diff --git a/sorc/global_cycle.fd/CMakeLists.txt b/sorc/global_cycle.fd/CMakeLists.txt index 0f73d4ab6..0af117441 100644 --- a/sorc/global_cycle.fd/CMakeLists.txt +++ b/sorc/global_cycle.fd/CMakeLists.txt @@ -14,7 +14,9 @@ set(lib_src set(exe_src cycle.f90) if(CMAKE_Fortran_COMPILER_ID MATCHES "^(Intel)$") - set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") +# set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -r8") +# Yuan Xue: set flag to suppress long name detector + set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -diag-disable:5462 -r8") elseif(CMAKE_Fortran_COMPILER_ID MATCHES "^(GNU)$") set(CMAKE_Fortran_FLAGS "${CMAKE_Fortran_FLAGS} -fdefault-real-8") if(CMAKE_Fortran_COMPILER_VERSION VERSION_GREATER_EQUAL 10) diff --git a/sorc/global_cycle.fd/cycle.f90 b/sorc/global_cycle.fd/cycle.f90 index 4511416c9..4607ac987 100644 --- a/sorc/global_cycle.fd/cycle.f90 +++ b/sorc/global_cycle.fd/cycle.f90 @@ -70,8 +70,9 @@ !! -DO_SFCCYCLE Call sfccycle routine to update surface fields !! -DO_LNDINC Read in land increment files, and add increments to !! relevant states. -!! -DO_SOI_INC Do land increments to soil states. -!! -DO_SNO_INC Do land increments to snow states. +!! -DO_SOI_INC_GSI Do land increments (on the Guassian grid) to soil states. +!! -DO_SOI_INC_JEDI Do land increments (on the cubed sphere tiles) to soil. +!! -DO_SNO_INC Do land increments (on the cubed sphere tiles) to snow states. !! - ISOT Use statsgo soil type when '1'. Use zobler when '0'. !! - IVEGSRC Use igbp veg type when '1'. Use sib when '2'. !! - ZSEA1/2_MM When running with NSST model, this is the lower/ @@ -302,10 +303,14 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & USE READ_WRITE_DATA use machine USE MPI - USE LAND_INCREMENTS, ONLY: ADD_INCREMENT_SOIL, & + USE NETCDF + USE LAND_INCREMENTS, ONLY: ADD_GSI_INCREMENT_SOIL, & + ADD_JEDI_INCREMENT_SOIL_TEMP, & + ADD_JEDI_INCREMENT_SOIL_MOISTURE, & ADD_INCREMENT_SNOW, & CALCULATE_LANDINC_MASK, & - APPLY_LAND_DA_ADJUSTMENTS_SOIL, & + APPLY_LAND_DA_ADJUSTMENTS_SOIL_TEMP, & + APPLY_LAND_DA_ADJUSTMENTS_SOIL_MOISTURE, & APPLY_LAND_DA_ADJUSTMENTS_SND, & LSM_NOAH, LSM_NOAHMP @@ -328,7 +333,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & CHARACTER(LEN=500) :: LND_SOI_FILE CHARACTER(LEN=4) :: INPUT_NML_FILE(SZ_NML) - INTEGER :: I, IERR + INTEGER :: I, IERR, NCID, ID_VAR INTEGER :: I_INDEX(LENSFC), J_INDEX(LENSFC) INTEGER :: IDUM(IDIM,JDIM) integer :: num_parthds, num_threads @@ -366,6 +371,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & REAL, ALLOCATABLE :: SLIFCS_FG(:) INTEGER, ALLOCATABLE :: LANDINC_MASK_FG(:), LANDINC_MASK(:) REAL, ALLOCATABLE :: SND_BCK(:), SND_INC(:), SWE_BCK(:) + REAL, ALLOCATABLE :: SLCINC(:,:), STCINC(:,:) REAL(KIND=KIND_IO8), ALLOCATABLE :: SLMASKL(:), SLMASKW(:) TYPE(NSST_DATA) :: NSST @@ -374,7 +380,9 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & INTEGER :: veg_type_landice INTEGER, DIMENSION(LENSFC) :: STC_UPDATED, SLC_UPDATED - LOGICAL :: FILE_EXISTS, DO_SOI_INC, DO_SNO_INC + REAL, DIMENSION(LENSFC,LSOIL) :: STC_INC_TMP, SLC_INC_TMP + + LOGICAL :: FILE_EXISTS, DO_SOI_INC_GSI, DO_SOI_INC_JEDI, DO_SNO_INC !-------------------------------------------------------------------------------- ! NST_FILE is the path/name of the gaussian GSI file which contains NSST @@ -387,7 +395,8 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & DATA LND_SOI_FILE/'NULL'/ DO_SNO_INC = .FALSE. - DO_SOI_INC = .FALSE. + DO_SOI_INC_GSI = .FALSE. + DO_SOI_INC_JEDI = .FALSE. SIG1T = 0.0 ! Not a dead start! @@ -446,13 +455,28 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & IF (DO_LNDINC) THEN ! identify variables to be updates, and allocate arrays. IF (TRIM(LND_SOI_FILE) .NE. "NULL") THEN - DO_SOI_INC = .TRUE. + DO_SOI_INC_GSI = .TRUE. PRINT* PRINT*," APPLYING SOIL INCREMENTS FROM THE GSI" - ALLOCATE(STC_BCK(LENSFC, LSOIL), SMC_BCK(LENSFC, LSOIL), SLC_BCK(LENSFC,LSOIL)) - ALLOCATE(LANDINC_MASK_FG(LENSFC)) + ELSE + INQUIRE(FILE='xainc.001', EXIST=file_exists) + IF (file_exists) then + IERR=NF90_OPEN('xainc.001',NF90_NOWRITE,NCID) + IERR=NF90_INQ_VARID(NCID,"stc_inc",ID_VAR) + IF (IERR == 0) THEN + DO_SOI_INC_JEDI = .TRUE. + PRINT* + PRINT*," APPLYING SOIL INCREMENTS FROM THE JEDI" + ALLOCATE(SLCINC(LENSFC, LSOIL), STCINC(LENSFC, LSOIL)) + ENDIF + ENDIF ENDIF + ALLOCATE(STC_BCK(LENSFC, LSOIL), SMC_BCK(LENSFC, LSOIL), SLC_BCK(LENSFC,LSOIL)) + ALLOCATE(LANDINC_MASK_FG(LENSFC)) ! FOR NOW, CODE SO CAN DO BOTH, BUT MIGHT NEED TO THINK ABOUT THIS SOME MORE. + ! NOTE: BOTH snow and soil related increment file are named as "xainc.00N" + ! This may create issues if doing soil and snow DA together + ! For now, just put an additional inquiry above to distinguish - Yuan Xue IF (DO_SNO_INC) THEN ! ideally, would check here that sfcsub snow DA update is not also requested ! but latter is controlled by fnsol, which is read in within that routine. @@ -474,6 +498,8 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & ! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE. !-------------------------------------------------------------------------------- +IF (.NOT. DO_SOI_INC_JEDI) THEN + CALL READ_DATA(LSOIL,LENSFC,DO_NSST,.false.,IS_NOAHMP=IS_NOAHMP, & TSFFCS=TSFFCS,SMCFCS=SMCFCS, & SWEFCS=SWEFCS,STCFCS=STCFCS,TG3FCS=TG3FCS,ZORFCS=ZORFCS, & @@ -485,6 +511,13 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & VMNFCS=VMNFCS,VMXFCS=VMXFCS,SLCFCS=SLCFCS,SLPFCS=SLPFCS, & ABSFCS=ABSFCS,T2M=T2M ,Q2M=Q2M ,SLMASK=SLMASK, & ZSOIL=ZSOIL, NSST=NSST) +ELSE + CALL READ_DATA(LSOIL,LENSFC,DO_NSST,.false.,IS_NOAHMP=IS_NOAHMP, & + VEGFCS=VEGFCS, SOTFCS=SOTFCS, & + SWEFCS=SWEFCS,SLIFCS=SLIFCS,VETFCS=VETFCS,SNDFCS=SNDFCS,& + SMCFCS=SMCFCS,SLCFCS=SLCFCS,STCFCS=STCFCS,SLMASK=SLMASK,ZSOIL=ZSOIL) + +ENDIF IF (IS_NOAHMP) THEN LSM=LSM_NOAHMP @@ -519,7 +552,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & ! CALCULATE MASK FOR LAND INCREMENTS IF (DO_LNDINC) & - CALL CALCULATE_LANDINC_MASK(SMCFCS(:,1),SWEFCS, VETFCS, & + CALL CALCULATE_LANDINC_MASK(SWEFCS, VETFCS, SOTFCS, & LENSFC,VEG_TYPE_LANDICE, LANDINC_MASK) !-------------------------------------------------------------------------------- @@ -648,7 +681,7 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & ENDIF ! SOIL INCREMENTS - IF (DO_SOI_INC) THEN + IF (DO_SOI_INC_GSI .OR. DO_SOI_INC_JEDI) THEN !-------------------------------------------------------------------------------- ! re-calculate soilsnow mask if snow has been updated. @@ -657,23 +690,26 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & LANDINC_MASK_FG = LANDINC_MASK IF (DO_SFCCYCLE .OR. DO_SNO_INC) THEN - CALL CALCULATE_LANDINC_MASK(SMCFCS(:,1),SWEFCS, VETFCS, LENSFC, & + CALL CALCULATE_LANDINC_MASK(SWEFCS, VETFCS, SOTFCS, LENSFC, & VEG_TYPE_LANDICE, LANDINC_MASK ) ENDIF - + !-------------------------------------------------------------------------------- ! read increments in !-------------------------------------------------------------------------------- - INQUIRE(FILE=trim(LND_SOI_FILE), EXIST=file_exists) - IF (.not. file_exists) then + IF (DO_SOI_INC_GSI) THEN + INQUIRE(FILE=trim(LND_SOI_FILE), EXIST=file_exists) + IF (.not. file_exists) then print *, 'FATAL ERROR: land increment update requested, but file does not exist: ', & trim(lnd_soi_file) call MPI_ABORT(MPI_COMM_WORLD, 10, IERR) + ENDIF + CALL READ_GSI_DATA(LND_SOI_FILE, 'LND', LSOIL=LSOIL) + ELSE IF (DO_SOI_INC_JEDI) THEN + CALL READ_DATA(LSOIL,LENSFC,.false.,.true.,SLCINC=SLCINC,STCINC=STCINC) ENDIF - CALL READ_GSI_DATA(LND_SOI_FILE, 'LND', LSOIL=LSOIL) - !-------------------------------------------------------------------------------- ! add increments to state vars !-------------------------------------------------------------------------------- @@ -688,17 +724,58 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & SLC_BCK = SLCFCS ! below updates [STC/SMC/STC]FCS to hold the analysis - CALL ADD_INCREMENT_SOIL(RLA,RLO,STCFCS,SMCFCS,SLCFCS,STC_UPDATED, SLC_UPDATED, & + ! ONLY isolate soil temperature and soil moisture ingest/adjust routines in + ! the jedi-related portion. In gsi-related part, it is *not* updated yet. + ! The gsi-related part is still using the old assimilation sequence + ! (i.e., soil moisture incr. ingest based on background state, not the + ! analysis state as edited in the new assimilation sequence) + ! The new assimilation/adjustment sequence in jedi is: + ! a) ingest soil temp incr. + ! b) adjust soil ice/liquid moisture as needed + ! c) ingest soil moisure incr. based on newly-adjusted soil state due to b) + ! d) ingest soil moisture incr. + ! e) adjust/bound soil moisture + IF (DO_SOI_INC_GSI) THEN + CALL ADD_GSI_INCREMENT_SOIL(RLA,RLO,STCFCS,SMCFCS,SLCFCS,STC_UPDATED, SLC_UPDATED, & + STC_INC_TMP, SLC_INC_TMP, & LANDINC_MASK,LANDINC_MASK_FG,LENSFC,LSOIL,IDIM,JDIM,LSM,MYRANK) !-------------------------------------------------------------------------------- ! make any necessary adjustments to dependent variables !-------------------------------------------------------------------------------- + CALL APPLY_LAND_DA_ADJUSTMENTS_SOIL_TEMP(LSM, ISOT, IVEGSRC, LENSFC,& + LSOIL, SOTFCS, LANDINC_MASK_FG,& + STC_BCK, STCFCS, SMCFCS, SLCFCS,& + STC_UPDATED) + CALL APPLY_LAND_DA_ADJUSTMENTS_SOIL_MOISTURE(LSM, LENSFC, & + LSOIL, SMCFCS, SLCFCS,& + SLC_UPDATED, ZSOIL) + + ! if read in GSI-based increments, save interpolated fields onto cubed sphere tiles + CALL WRITE_DATA(LENSFC,IDIM,JDIM,LSOIL,DO_NSST,.true.,NSST, & + SLC_INC_TMP=SLC_INC_TMP,STC_INC_TMP=STC_INC_TMP) + + ELSE IF (DO_SOI_INC_JEDI) THEN + CALL ADD_JEDI_INCREMENT_SOIL_TEMP(STCINC, STCFCS, & + STC_UPDATED, & + LANDINC_MASK, LANDINC_MASK_FG, LENSFC,& + LSOIL, LSM, MYRANK) + CALL APPLY_LAND_DA_ADJUSTMENTS_SOIL_TEMP(LSM, ISOT, IVEGSRC, LENSFC,& + LSOIL, SOTFCS, LANDINC_MASK_FG,& + STC_BCK, STCFCS, SMCFCS, SLCFCS,& + STC_UPDATED) + + CALL ADD_JEDI_INCREMENT_SOIL_MOISTURE(SLCINC, STCFCS, SMCFCS, SLCFCS,& + SLC_UPDATED, & + LANDINC_MASK, LANDINC_MASK_FG, LENSFC,& + LSOIL, LSM, MYRANK) + CALL APPLY_LAND_DA_ADJUSTMENTS_SOIL_MOISTURE(LSM, LENSFC, & + LSOIL, SMCFCS, SLCFCS,& + SLC_UPDATED, ZSOIL) + + ENDIF ! switch between reading gsi and jedi based increments - CALL APPLY_LAND_DA_ADJUSTMENTS_SOIL(LSM, ISOT, IVEGSRC,LENSFC, LSOIL, & - SOTFCS, LANDINC_MASK_FG, STC_BCK, STCFCS, SMCFCS, SLCFCS, STC_UPDATED, & - SLC_UPDATED,ZSOIL) ENDIF ! soil increments @@ -706,7 +783,6 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & ! clean up !-------------------------------------------------------------------------------- - ! to do - save and write out STC_INC? (soil temperature increments) IF(ALLOCATED(LANDINC_MASK_FG)) DEALLOCATE(LANDINC_MASK_FG) IF(ALLOCATED(LANDINC_MASK)) DEALLOCATE(LANDINC_MASK) IF(ALLOCATED(STC_BCK)) DEALLOCATE(STC_BCK) @@ -715,6 +791,8 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & IF(ALLOCATED(SND_BCK)) DEALLOCATE(SND_BCK) IF(ALLOCATED(SWE_BCK)) DEALLOCATE(SWE_BCK) IF(ALLOCATED(SND_INC)) DEALLOCATE(SND_INC) + IF(ALLOCATED(SLCINC)) DEALLOCATE(SLCINC) + IF(ALLOCATED(STCINC)) DEALLOCATE(STCINC) ENDIF !-------------------------------------------------------------------------------- @@ -723,13 +801,13 @@ SUBROUTINE SFCDRV(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, & IF (LSM==LSM_NOAHMP) THEN - CALL WRITE_DATA(LENSFC,IDIM,JDIM,LSOIL,DO_NSST,NSST,VEGFCS=VEGFCS, & + CALL WRITE_DATA(LENSFC,IDIM,JDIM,LSOIL,DO_NSST,.false.,NSST,VEGFCS=VEGFCS, & SLCFCS=SLCFCS,SMCFCS=SMCFCS,STCFCS=STCFCS) ELSEIF (LSM==LSM_NOAH) THEN CALL WRITE_DATA(LENSFC,IDIM,JDIM,LSOIL, & - DO_NSST,NSST,SLIFCS=SLIFCS,TSFFCS=TSFFCS,VEGFCS=VEGFCS, & + DO_NSST,.false.,NSST,SLIFCS=SLIFCS,TSFFCS=TSFFCS,VEGFCS=VEGFCS, & SWEFCS=SWEFCS,TG3FCS=TG3FCS,ZORFCS=ZORFCS, & ALBFCS=ALBFCS,ALFFCS=ALFFCS,CNPFCS=CNPFCS, & F10M=F10M,T2M=T2M,Q2M=Q2M,VETFCS=VETFCS, & @@ -831,7 +909,7 @@ SUBROUTINE ADJUST_NSST(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,& INTEGER :: ISTART, IEND, JSTART, JEND INTEGER :: MASK_TILE, MASK_FG_TILE INTEGER :: ITILE, JTILE - INTEGER :: MAX_SEARCH, J, IERR + INTEGER :: MAX_SEARCH, J, IERR INTEGER :: IGAUSP1, JGAUSP1 integer :: nintp,nsearched,nice,nland integer :: nfill,nfill_tice,nfill_clm diff --git a/sorc/global_cycle.fd/land_increments.f90 b/sorc/global_cycle.fd/land_increments.f90 index b9738250c..9cab1cc9b 100644 --- a/sorc/global_cycle.fd/land_increments.f90 +++ b/sorc/global_cycle.fd/land_increments.f90 @@ -6,10 +6,13 @@ module land_increments private - public add_increment_soil + public add_gsi_increment_soil + public add_jedi_increment_soil_temp + public add_jedi_increment_soil_moisture public add_increment_snow public calculate_landinc_mask - public apply_land_da_adjustments_soil + public apply_land_da_adjustments_soil_temp + public apply_land_da_adjustments_soil_moisture public apply_land_da_adjustments_snd public lsm_noah, lsm_noahmp @@ -18,7 +21,15 @@ module land_increments !! copied from GFS_typedefs.F90 ! control state for soil analysis: - integer, parameter :: lsoil_incr=3 !< number of layers to add incrments to + ! Yuan Xue: ONLY updating top 2 layers because analysis shown that + ! Pan-Arctic permafrost is too sensitive to slc adjustment due to soil temperature ingest in 3rd layer + ! whenever soil temperature approaches freezing point during a synthetic + ! twin experiment. + ! This is indeed a generic issue due to opt_frz=1 (Niu and Yang, 2006 jhm) + ! A +/-0.05K difference in soil temp (slightly below tfreez) can lead to + ! 0.05-0.06m3/m3 difference in liquid soil moisture, which is also seen in the + ! surface layers. + integer, parameter :: lsoil_incr=2 !< number of layers to add incrments to real, parameter :: tfreez=273.16 !< con_t0c in physcons contains @@ -49,10 +60,14 @@ module land_increments !! @param[in] JDIM 'J' dimension of a tile !! @param[in] lsm Integer flag indicating which land model is used (1-Noah, 2-Noah-MP) !! @param[in] MYRANK MPI rank number - !! + !! @param[out] stc_inc_tmp Soil temperature increments on the cubed-sphere tile + !! @param[out] slc_inc_tmp Liquid soil moisture increments on the cubed-sphere tile !! @author Clara Draper. @date March 2021 + !! @author Yuan Xue. @date 11/08/2023: add capability to save interpolated + !! GSI increments onto cubed-sphere tiles -subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, slc_updated, & +subroutine add_gsi_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, slc_updated, & + stc_inc_tmp, slc_inc_tmp, & soilsnow_tile,soilsnow_fg_tile,lensfc,lsoil,idim,jdim,lsm, myrank) use utils @@ -72,6 +87,9 @@ subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, real, intent(inout) :: smc_state(lensfc, lsoil) integer, intent(out) :: stc_updated(lensfc), slc_updated(lensfc) + real :: slc_inc_tmp(lensfc, lsoil) + real :: stc_inc_tmp(lensfc, lsoil) + integer :: iopt, nret, kgds_gaus(200) integer :: igaus, jgaus, ij integer :: mask_tile, mask_fg_tile @@ -203,6 +221,9 @@ subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, nfrozen = 0 ! not update as frozen soil nfrozen_upd = 0 ! not update as frozen soil + ! initialize matrix + stc_inc_tmp = 0.0 + slc_inc_tmp = 0.0 ij_loop : do ij = 1, lensfc @@ -301,10 +322,12 @@ subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, wsum = wsum + s2c(itile,jtile,4) endif - ! normalize increments + ! normalize increments and save them do k = 1, lsoil_incr stc_inc(k) = stc_inc(k) / wsum + stc_inc_tmp(ij,k) = stc_inc(k) slc_inc(k) = slc_inc(k) / wsum + slc_inc_tmp(ij,k) = slc_inc(k) enddo !---------------------------------------------------------------------- ! add the interpolated increment to the background @@ -361,7 +384,245 @@ subroutine add_increment_soil(rla,rlo,stc_state,smc_state,slc_state,stc_updated, deallocate(id1, id2, jdc, s2c) -end subroutine add_increment_soil +end subroutine add_gsi_increment_soil + + !> Similar to gsi based increment ingest, this routine is to + !! add jedi based stc increments to model state + !! Without any follow-up slc adjustment, add_gsi_increment_soil is identical to + !! add_jedi_increment_soil_temp plus add_jedi_increment_soil_moisture -- sanity checks are all done + !! @param[in] STCINC Soil temperature increments on the cubed-sphere tile + !! @param[inout] STC_STATE Soil temperature state vector + !! @param[out] STC_UPDATED Integer to record if STC in each grid is updated + !! @param[in] SOILSNOW_TILE Land mask for increments on the cubed-sphere tile + !! @param[in] SOILSNOW_FG_TILE First guess land mask for increments on the + !! cubed-sphere tile + !! @param[in] LENSFC Number of land points on a tile + !! @param[in] LSOIL Number of soil layers + !! @param[in] lsm Integer flag indicating which land model is used (1-Noah, + !! 2-Noah-MP) + !! @param[in] MYRANK MPI rank number + !! @author Yuan Xue @date: 11/08/2023 +subroutine add_jedi_increment_soil_temp(stcinc,stc_state,& + stc_updated,& + soilsnow_tile,soilsnow_fg_tile,lensfc,lsoil,& + lsm,myrank) + use mpi + + implicit none + + integer, intent(in) :: lensfc, lsoil, myrank, lsm + integer, intent(in) :: soilsnow_tile(lensfc), soilsnow_fg_tile(lensfc) + real, intent(in) :: stcinc(lensfc,lsoil) + real, intent(inout) :: stc_state(lensfc, lsoil) + integer, intent(out) :: stc_updated(lensfc) + + integer :: ij + integer :: mask_tile, mask_fg_tile + logical :: upd_stc + + integer :: k, nother, nsnowupd + integer :: nstcupd + + stc_updated=0 + + if (lsm==lsm_noah) then + upd_stc=.true. + elseif (lsm==lsm_noahmp) then + upd_stc=.true. + endif + + print* + print*,'adjust soil temp using jedi increments on cubed sphere tiles' + print*,'updating soil temps', upd_stc + print*,'adjusting first ', lsoil_incr, ' surface layers only' + + ! initialize variables for counts statitics to be zeros + nother = 0 ! grid cells not land + nsnowupd = 0 ! grid cells with snow (temperature not yet updated) + nstcupd = 0 ! grid cells that are updated + + ij_loop : do ij = 1, lensfc + + mask_tile = soilsnow_tile(ij) + mask_fg_tile = soilsnow_fg_tile(ij) + + !---------------------------------------------------------------------- + ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land + !---------------------------------------------------------------------- + + if (mask_tile <= 0) then ! skip if neither soil nor snow + nother = nother + 1 + cycle ij_loop + endif + + !---------------------------------------------------------------------- + ! if snow is present before or after snow update, skip soil analysis + !---------------------------------------------------------------------- + + if (mask_fg_tile == 2 .or. mask_tile == 2) then + nsnowupd = nsnowupd + 1 + cycle ij_loop + endif + + if (mask_tile == 1) then + !---------------------------------------------------------------------- + ! add the interpolated increment to the background + !---------------------------------------------------------------------- + do k = 1, lsoil_incr + + if (upd_stc) then + stc_state(ij,k) = stc_state(ij,k) + stcinc(ij,k) + if (k==1) then + stc_updated(ij) = 1 + nstcupd = nstcupd + 1 + endif + endif + + enddo + + endif ! if soil/snow point + + enddo ij_loop + + write(*,'(a,i2)') ' statistics of grids number processed for rank : ',myrank + write(*,'(a,i8)') ' soil grid total', lensfc + write(*,'(a,i8)') ' soil grid cells stc updated = ',nstcupd + write(*,'(a,i8)') ' (not updated yet) snow grid cells = ', nsnowupd + write(*,'(a,i8)') ' grid cells, without soil or snow = ', nother + +end subroutine add_jedi_increment_soil_temp + +!> Add jedi based slc (liquid soil moisture) increments to model state +!! @param[in] SLCINC Liquid soil moisture increments on the cubed-sphere tile +!! @param[in] STC_STATE Soil temperature state vector +!! @param[inout] SMC_STATE Total soil moisture state vector +!! @param[inout] SLC_STATE Liquid soil moisture state vector +!! @param[out] SLC_UPDATED Integer to record if SLC in each grid is updated +!! @param[in] SOILSNOW_TILE Land mask for increments on the cubed-sphere tile +!! @param[in] SOILSNOW_FG_TILE First guess land mask for increments on the +!! cubed-sphere tile +!! @param[in] LENSFC Number of land points on a tile +!! @param[in] LSOIL Number of soil layers +!! @param[in] lsm Integer flag indicating which land model is used (1-Noah, +!! 2-Noah-MP) +!! @param[in] MYRANK MPI rank number +!! @author Yuan Xue @date: 11/08/2023 +subroutine add_jedi_increment_soil_moisture(slcinc,stc_state,smc_state,slc_state,& + slc_updated, & + soilsnow_tile,soilsnow_fg_tile,lensfc,lsoil,& + lsm,myrank) + use mpi + + implicit none + + integer, intent(in) :: lensfc, lsoil, myrank, lsm + integer, intent(in) :: soilsnow_tile(lensfc), soilsnow_fg_tile(lensfc) + real, intent(in) :: slcinc(lensfc,lsoil) + real, intent(in) :: stc_state(lensfc, lsoil) + real, intent(inout) :: slc_state(lensfc, lsoil) + real, intent(inout) :: smc_state(lensfc, lsoil) + integer, intent(out) :: slc_updated(lensfc) + + integer :: ij + integer :: mask_tile, mask_fg_tile + logical :: upd_slc + + integer :: k, nother, nsnowupd + integer :: nslcupd, nfrozen, nfrozen_upd + logical :: soil_freeze, soil_ice + + + slc_updated=0 + + if (lsm==lsm_noah) then + upd_slc=.false. ! not coded + elseif (lsm==lsm_noahmp) then + upd_slc=.true. + endif + + print* + print*,'adjust soil moisture using jedi increments on cubed sphere tiles' + print*,'updating soil moisture', upd_slc + print*,'adjusting first ', lsoil_incr, ' surface layers only' + + ! initialize variables for counts statitics to be zeros + nother = 0 ! grid cells not land + nsnowupd = 0 ! grid cells with snow (not updated) + nslcupd = 0 ! grid cells that are updated + nfrozen = 0 ! not update as frozen soil + nfrozen_upd = 0 ! not update as frozen soil + + ij_loop : do ij = 1, lensfc + + mask_tile = soilsnow_tile(ij) + mask_fg_tile = soilsnow_fg_tile(ij) + + !---------------------------------------------------------------------- + ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land + !---------------------------------------------------------------------- + + if (mask_tile <= 0) then ! skip if neither soil nor snow + nother = nother + 1 + cycle ij_loop + endif + + !---------------------------------------------------------------------- + ! if snow is present before or after snow update, skip soil analysis + !---------------------------------------------------------------------- + + if (mask_fg_tile == 2 .or. mask_tile == 2) then + nsnowupd = nsnowupd + 1 + cycle ij_loop + endif + + if (mask_tile == 1) then + !---------------------------------------------------------------------- + ! add the interpolated increment to the background + !---------------------------------------------------------------------- + + soil_freeze=.false. + soil_ice=.false. + do k = 1, lsoil_incr + + if ( stc_state(ij,k) < tfreez) soil_freeze=.true. + if ( smc_state(ij,k) - slc_state(ij,k) > 0.001 ) soil_ice=.true. + + if ( (stc_state(ij,k) < tfreez) .and. (.not. soil_freeze) .and. & + (k==1) ) & + nfrozen_upd = nfrozen_upd + 1 + + ! do not do updates if this layer or any above is frozen + if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then + if (upd_slc) then + if (k==1) then + nslcupd = nslcupd + 1 + slc_updated(ij) = 1 + endif + ! apply zero limit here (higher, model-specific limits are + ! later) + slc_state(ij,k) = max(slc_state(ij,k) + slcinc(ij,k), 0.0) + smc_state(ij,k) = max(smc_state(ij,k) + slcinc(ij,k), 0.0) + endif + else + if (k==1) nfrozen = nfrozen+1 + endif + + enddo + + endif ! if soil/snow point + + enddo ij_loop + + write(*,'(a,i2)') ' statistics of grids number processed for rank : ',myrank + write(*,'(a,i8)') ' soil grid total', lensfc + write(*,'(a,i8)') ' soil grid cells slc updated = ',nslcupd + write(*,'(a,i8)') ' soil grid cells not updated, frozen = ',nfrozen + write(*,'(a,i8)') ' soil grid cells update, became frozen = ',nfrozen_upd + write(*,'(a,i8)') ' (not updated yet) snow grid cells = ', nsnowupd + write(*,'(a,i8)') ' grid cells, without soil or snow = ', nother + +end subroutine add_jedi_increment_soil_moisture + !> Add snow depth increment to model snow depth state, !! and limit output to be non-negative. JEDI increments are @@ -400,18 +661,19 @@ end subroutine add_increment_snow !! !! @param[in] lensfc Number of land points for this tile !! @param[in] veg_type_landice Value of vegetion class that indicates land-ice -!! @param[in] smc Model soil moisture. +!! @param[in] stype Soil type !! @param[in] swe Model snow water equivalent !! @param[in] vtype Model vegetation type !! @param[out] mask Land mask for increments !! @author Clara Draper @date March 2021 -subroutine calculate_landinc_mask(smc,swe,vtype,lensfc,veg_type_landice,mask) +!! @author Yuan Xue: introduce stype to make the mask calculation more generic +subroutine calculate_landinc_mask(swe,vtype,stype,lensfc,veg_type_landice,mask) implicit none integer, intent(in) :: lensfc, veg_type_landice - real, intent(in) :: smc(lensfc), swe(lensfc) - real, intent(in) :: vtype(lensfc) + real, intent(in) :: swe(lensfc) + real, intent(in) :: vtype(lensfc),stype(lensfc) integer, intent(out) :: mask(lensfc) integer :: i @@ -420,7 +682,7 @@ subroutine calculate_landinc_mask(smc,swe,vtype,lensfc,veg_type_landice,mask) ! land (but not land-ice) do i=1,lensfc - if (smc(i) .LT. 0.99) then + if (nint(stype(i)) .GT. 0) then if (swe(i) .GT. 0.001) then ! snow covered land mask(i) = 2 else ! non-snow covered land @@ -434,18 +696,17 @@ subroutine calculate_landinc_mask(smc,swe,vtype,lensfc,veg_type_landice,mask) end subroutine calculate_landinc_mask -!> Make adjustments to dependent variables after applying land increments. -!! These adjustments are model-dependent, and are currently only coded -!! if full for Noah LSM. +!> Make adjustments to dependent variables after applying soil temp increments +!! These adjustments are model-dependent !! For Noah LSM, copy relevent code blocks from model code (same as has !! been done in sfc_sub). -!! For Noah-MP, have inserted place-holders to simply reset the model -!! variables back to the analysis if adjustments are needed. Later, will replace -!! this with appropriate adjustmenets (in summary, for now we do not -!! make STC updates if soils are frozen, and are also not applying the -!! appropriate max. values for SMC). -!! Here: adjust (frozen) soil moisture to be consistent with changes in -!! soil temperature from DA +!! For Noah-MP, the adjustment scheme shown below as of 11/09/2023: +!! Case 1: frozen ==> frozen, recalculate slc following opt_frz=1, smc remains +!! Case 2: unfrozen ==> frozen, recalculate slc following opt_frz=1, smc remains +!! Case 3: frozen ==> unfrozen, melt all soil ice (if any) +!! Case 4: unfrozen ==> unfrozen along with other cases (e.g., temp=tfreez), do nothing +!! Note: For Case 3, Yuan Xue thoroughly evaluated a total of four options and the +!! current option is found to be the best as of 11/09/2023 !! @param[in] lsm Integer code for the LSM !! @param[in] isot Integer code for the soil type data set !! @param[in] ivegsrc Integer code for the vegetation type data set @@ -454,72 +715,73 @@ end subroutine calculate_landinc_mask !! @param[in] rsoiltype Array of input soil types !! @param[in] mask Mask indicating surface type !! @param[in] stc_bck Background soil temperature states -!! @param[in] stc_adj Analysis soil temperature states +!! @param[inout] stc_adj Analysis soil temperature states !! @param[inout] smc_adj Analysis soil moisture states !! @param[inout] slc_adj Analysis liquid soil moisture states !! @param[in] stc_updated Integer to record whether STC in each grid cell was updated -!! @param[in] slc_updated Integer to record whether SLC in each grid cell was updated -!! @param[in] zsoil Depth of bottom of each soil layer !! @author Clara Draper @date April 2021 +!! @author Yuan Xue: isolate soil temperature and soil moisture adjustment +!! routines @ date Sept 2023 -subroutine apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, & +subroutine apply_land_da_adjustments_soil_temp(lsm, isot, ivegsrc,lensfc, & lsoil, rsoiltype, mask, stc_bck, stc_adj, smc_adj, slc_adj, & - stc_updated, slc_updated, zsoil) + stc_updated) use mpi - use set_soilveg_snippet_mod, only: set_soilveg + use set_soilveg_snippet_mod, only: set_soilveg_noah, set_soilveg_noahmp use sflx_snippet, only: frh2o implicit none - + integer, intent(in) :: lsm, lensfc, lsoil, isot, ivegsrc real, intent(in) :: rsoiltype(lensfc) ! soil types, as real integer, intent(in) :: mask(lensfc) real, intent(in) :: stc_bck(lensfc, lsoil) - integer, intent(in) :: stc_updated(lensfc), slc_updated(lensfc) - real, intent(inout) :: smc_adj(lensfc,lsoil), slc_adj(lensfc,lsoil) + integer, intent(in) :: stc_updated(lensfc) + real, intent(inout) :: smc_adj(lensfc,lsoil),slc_adj(lensfc,lsoil) real, intent(inout) :: stc_adj(lensfc, lsoil) - real(kind=4), intent(in) :: zsoil(lsoil) - + logical :: frzn_bck, frzn_anl logical :: soil_freeze, soil_ice - integer :: i, l, n_freeze, n_thaw, ierr, n_revert - integer :: myrank, soiltype, iret, n_stc, n_slc - logical :: upd_slc, upd_stc + integer :: i, l, n_freeze, n_thaw, ierr + integer :: myrank, soiltype, iret, n_stc + logical :: upd_stc real :: slc_new real, parameter :: tfreez=273.16 !< con_t0c in physcons + real, dimension(30) :: maxsmc, bb, satpsi - real, dimension(4) :: dz ! layer thickness + + real, parameter :: hfus=0.3336e06 !< latent heat of fusion(j/kg) + real, parameter :: grav=9.80616 !< gravity accel.(m/s2) + real :: smp !< for computing supercooled water within + !!frozen soil regime call mpi_comm_rank(mpi_comm_world, myrank, ierr) - if (lsm==lsm_noah) then + if (lsm==lsm_noah) then upd_stc=.true. - upd_slc=.false. - elseif (lsm==lsm_noahmp) then + elseif (lsm==lsm_noahmp) then upd_stc=.true. - upd_slc=.true. endif - select case (lsm ) - case(lsm_noah) - ! initialise soil properties - call set_soilveg(isot, ivegsrc, maxsmc, bb, satpsi, iret) - if (iret < 0) then - print *, 'FATAL ERROR: problem in set_soilveg' + select case (lsm ) + case(lsm_noah) + call set_soilveg_noah(isot, ivegsrc, maxsmc, bb, satpsi, iret) + if (iret < 0) then + print *, 'FATAL ERROR: problem in set_soilveg_noah' call mpi_abort(mpi_comm_world, 10, ierr) - endif + endif - print *, 'Adjusting noah model smc after stc DA update' + print *, 'Adjusting noah model smc after stc DA update' - n_freeze = 0 - n_thaw = 0 - - do i=1,lensfc + n_freeze = 0 + n_thaw = 0 + + do i=1,lensfc if (mask(i) > 0) then ! if soil location do l = 1, lsoil frzn_bck = (stc_bck(i,l) .LT. tfreez ) @@ -536,8 +798,8 @@ subroutine apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, & ! make adjustment (same routine for both) soiltype = nint(rsoiltype(i)) ! bb and maxsmc are in the namelist_soilveg, need soiltype index - call frh2o(stc_adj(i,l), smc_adj(i,l),slc_adj(i,l), maxsmc(soiltype), & - bb(soiltype), satpsi(soiltype),slc_new) + call frh2o(stc_adj(i,l), smc_adj(i,l),slc_adj(i,l),& + maxsmc(soiltype), bb(soiltype), satpsi(soiltype),slc_new) slc_adj(i,l) = max( min( slc_new, smc_adj(i,l)), 0.0 ) enddo @@ -545,44 +807,109 @@ subroutine apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, & enddo print *, 'adjusted: ', n_thaw,' thawed,', n_freeze, ' frozen' - case (lsm_noahmp) + + case (lsm_noahmp) - if (upd_stc) then + if (upd_stc) then + n_stc=0 + call set_soilveg_noahmp(isot, ivegsrc, maxsmc, bb, satpsi, iret) + if (iret < 0) then + print *, 'FATAL ERROR: problem in set_soilveg_noahmp' + call mpi_abort(mpi_comm_world, 10, ierr) + endif - print *, 'Reverting frozen noah-mp model stc back to background' - n_revert=0 - n_stc = 0 - n_slc = 0 + print *, 'Adjusting noah-mp model slc/stc after stc DA update' - do i=1,lensfc + do i=1,lensfc if (stc_updated(i) == 1 ) then n_stc = n_stc+1 - ! remove soil temperature increments if frozen - soil_freeze=.false. - soil_ice=.false. - do l = 1, lsoil_incr - if ( min(stc_bck(i,l),stc_adj(i,l)) < tfreez) soil_freeze=.true. - if ( smc_adj(i,l) - slc_adj(i,l) > 0.001 ) soil_ice=.true. - if ( soil_freeze .or. soil_ice ) then - ! for now, revert update. Later, adjust SMC/SLC for update. - if (l==1) n_revert = n_revert+1 - stc_adj(i,l)=stc_bck(i,l) - endif - enddo - endif - enddo + soiltype = nint(rsoiltype(i)) + if (mask(i) == 1) then + do l = 1, lsoil_incr !Only adjusts top two layers + !case 1: frz ==> frz, recalculate slc, smc remains + !case 2: unfrz ==> frz, recalculate slc, smc remains + !both cases are considered in the following if case + if (stc_adj(i,l) .LT. tfreez )then + !recompute supercool liquid water,smc_anl remain unchanged + smp = hfus*(tfreez-stc_adj(i,l))/(grav*stc_adj(i,l)) !(m) + slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype)) !(m3/m3) + slc_adj(i,l) = max( min( slc_new, smc_adj(i,l)), 0.0 ) + endif + !case 3: frz ==> unfrz, melt all soil ice (if any) + if ((stc_bck(i,l) .LT. tfreez) .and. & + (stc_adj(i,l) .GT. tfreez))then + slc_adj(i,l)=smc_adj(i,l) + endif + enddo + endif + endif + enddo + + endif + + case default + print *, 'FATAL ERROR: unrecognised LSM,', lsm + call mpi_abort(mpi_comm_world, 10, ierr) + end select - endif - if (upd_slc) then + write(*,'(a,i2)') 'statistics of grids number processed for rank : ', myrank + write(*,'(a,i8)') ' soil grid total', lensfc + write(*,'(a,i8)') ' soil grid cells with stc update', n_stc - dz(1) = -zsoil(1) - do l = 2,lsoil - dz(l) = -zsoil(l) + zsoil(l-1) - enddo - print *, 'Applying soil moisture mins ' +end subroutine apply_land_da_adjustments_soil_temp + +!> Make adjustments to dependent variables after applying soil moisture increments +!! These adjustments are model-dependent +!! For Noah LSM, no soil moisture adjustment is applied as of 11/09/2023 +!! For Noah-MP, the adjustment is applying soil moisture mins as of 11/09/2023 +!! @param[in] lsm Integer code for the LSM +!! @param[in] lensfc Number of land points for this tile +!! @param[in] lsoil Number of soil layers +!! @param[inout] smc_adj Analysis soil moisture states +!! @param[inout] slc_adj Analysis liquid soil moisture states +!! @param[in] slc_updated Integer to record if SLC in each grid cell was updated +!! @param[in] zsoil Depth of bottom of each soil layer + +subroutine apply_land_da_adjustments_soil_moisture(lsm,lensfc, & + lsoil, smc_adj, slc_adj, & + slc_updated, zsoil) + + use mpi + + implicit none + + integer, intent(in) :: lsm, lensfc, lsoil + integer, intent(in) :: slc_updated(lensfc) + real, intent(inout) :: smc_adj(lensfc,lsoil) + real, intent(inout) :: slc_adj(lensfc,lsoil) + real(kind=4), intent(in) :: zsoil(lsoil) + + integer :: i, l, ierr + integer :: myrank, n_slc + logical :: upd_slc + real, dimension(4) :: dz ! layer thickness + + call mpi_comm_rank(mpi_comm_world, myrank, ierr) + + if (lsm==lsm_noah) then + upd_slc=.false. + elseif (lsm==lsm_noahmp) then + upd_slc=.true. + endif + + select case (lsm) + case(lsm_noahmp) + + if (upd_slc) then + n_slc = 0 + dz(1) = -zsoil(1) + do l = 2,lsoil + dz(l) = -zsoil(l) + zsoil(l-1) + enddo + print *, 'Applying soil moisture mins ' do i=1,lensfc - if (slc_updated(i) == 1 ) then + if (slc_updated(i) == 1 ) then n_slc = n_slc+1 ! apply SM bounds (later: add upper SMC limit) do l = 1, lsoil_incr @@ -593,20 +920,19 @@ subroutine apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, & enddo endif enddo - endif + endif - case default + case default print *, 'FATAL ERROR: unrecognised LSM,', lsm call mpi_abort(mpi_comm_world, 10, ierr) end select - write(*,'(a,i2)') 'statistics of grids number processed for rank : ', myrank + write(*,'(a,i2)') 'statistics of grids number processed for rank : ', myrank write(*,'(a,i8)') ' soil grid total', lensfc write(*,'(a,i8)') ' soil grid cells with slc update', n_slc - write(*,'(a,i8)') ' soil grid cells with stc update', n_stc - write(*,'(a,i8)') ' soil grid cells reverted', n_revert -end subroutine apply_land_da_adjustments_soil +end subroutine apply_land_da_adjustments_soil_moisture + !> Make adjustments to dependent variables after applying land increments. !! These adjustments are model-dependent, and are currently only coded diff --git a/sorc/global_cycle.fd/read_write_data.f90 b/sorc/global_cycle.fd/read_write_data.f90 index 1cc9ad90f..24c75ac10 100644 --- a/sorc/global_cycle.fd/read_write_data.f90 +++ b/sorc/global_cycle.fd/read_write_data.f90 @@ -79,6 +79,7 @@ MODULE READ_WRITE_DATA !! @param[in] lensfc Total number of points on a tile. !! @param[in] lsoil Number of soil layers. !! @param[in] do_nsst When true, nsst fields were processed. + !! @param[in] inc_file When true, write out increments to files !! @param[in] nsst Data structure containing nsst fields. !! @param[in] slifcs Land-sea mask. !! @param[in] tsffcs Skin temperature. @@ -114,17 +115,21 @@ MODULE READ_WRITE_DATA !! @param[in] slcfcs Liquid portion of volumetric soil moisture. !! @param[in] smcfcs Total volumetric soil moisture. !! @param[in] stcfcs Soil temperature. - !! + !! @param[in] stc_inc_tmp Soil temperature increments on the cubed-sphere tiles + !! @param[in] slc_inc_tmp Liquid soil moisture increments on the cubed-sphere tiles !! @author George Gayno NOAA/EMC - + !! @author Yuan Xue: gsi-based soil increments are interpolated onto the + !! cubed sphered tiles for assimilation, hence, add capability to write out + !! interpolated increments into files subroutine write_data(lensfc,idim,jdim,lsoil, & - do_nsst,nsst,slifcs,tsffcs,vegfcs,swefcs, & + do_nsst,inc_file,nsst,slifcs,tsffcs,vegfcs,swefcs, & tg3fcs,zorfcs,albfcs,alffcs, & cnpfcs,f10m,t2m,q2m,vetfcs, & sotfcs,ustar,fmm,fhh,sicfcs, & sihfcs,sitfcs,tprcp,srflag, & swdfcs,vmnfcs,vmxfcs,slpfcs, & - absfcs,slcfcs,smcfcs,stcfcs) + absfcs,slcfcs,smcfcs,stcfcs, & + stc_inc_tmp,slc_inc_tmp) use mpi @@ -134,6 +139,7 @@ subroutine write_data(lensfc,idim,jdim,lsoil, & integer, intent(in) :: idim, jdim logical, intent(in) :: do_nsst + logical, intent(in) :: inc_file real, intent(in), optional :: slifcs(lensfc),tsffcs(lensfc) real, intent(in), optional :: swefcs(lensfc),tg3fcs(lensfc) @@ -150,6 +156,8 @@ subroutine write_data(lensfc,idim,jdim,lsoil, & real, intent(in), optional :: vmxfcs(lensfc), slpfcs(lensfc) real, intent(in), optional :: absfcs(lensfc), slcfcs(lensfc,lsoil) real, intent(in), optional :: smcfcs(lensfc,lsoil), stcfcs(lensfc,lsoil) + real, intent(in), optional :: stc_inc_tmp(lensfc,lsoil) + real, intent(in), optional :: slc_inc_tmp(lensfc,lsoil) type(nsst_data), intent(in) :: nsst @@ -161,18 +169,21 @@ subroutine write_data(lensfc,idim,jdim,lsoil, & character(len=3) :: rankch integer :: myrank, error, ncid, id_var + integer :: varid_slc, varid_stc + integer :: dim_soil call mpi_comm_rank(mpi_comm_world, myrank, error) write(rankch, '(i3.3)') (myrank+1) - fnbgso = "./fnbgso." // rankch + if (.NOT.(inc_file)) then - print* - print*,"update OUTPUT SFC DATA TO: ",trim(fnbgso) + fnbgso = "./fnbgso." // rankch + print* + print*,"update OUTPUT SFC DATA TO: ",trim(fnbgso) - ERROR=NF90_OPEN(TRIM(fnbgso),NF90_WRITE,NCID) - CALL NETCDF_ERR(ERROR, 'OPENING FILE: '//TRIM(fnbgso) ) + ERROR=NF90_OPEN(TRIM(fnbgso),NF90_WRITE,NCID) + CALL NETCDF_ERR(ERROR, 'OPENING FILE: '//TRIM(fnbgso) ) if(present(slifcs)) then error=nf90_inq_varid(ncid, "slmsk", id_var) @@ -472,6 +483,47 @@ subroutine write_data(lensfc,idim,jdim,lsoil, & call remove_checksum(ncid, id_var) endif + else + + fnbgso = "./xainc." // rankch + print* + print*,"Write increments onto cubed sphere tiles to: ", trim(fnbgso) + + error=nf90_create(trim(fnbgso),NF90_64BIT_OFFSET,ncid) + CALL netcdf_err(error, 'OPENING FILE: '//trim(fnbgso) ) + + ! Define dimensions in the file. + error = nf90_def_dim(ncid, "xaxis_1", idim, dim_x) + call netcdf_err(error, 'defining xaxis_1') + + error = nf90_def_dim(ncid, "yaxis_1", jdim, dim_y) + call netcdf_err(error, 'defining yaxis_1') + + error = nf90_def_dim(ncid, "soil_levels",lsoil, dim_soil) + call netcdf_err(error, 'defining soil_levels') + + ! Define variables in the file. + error=nf90_def_var(ncid, "slc_inc", NF90_DOUBLE, & + (/dim_x,dim_y,dim_soil/),varid_slc) + call netcdf_err(error, 'defining slc_inc'); + + error=nf90_def_var(ncid, "stc_inc", NF90_DOUBLE, & + (/dim_x,dim_y,dim_soil/),varid_stc) + call netcdf_err(error, 'defining stc_inc'); + + error = nf90_enddef(ncid) + + ! Put variables in the file. + dum3d = reshape(stc_inc_tmp, (/idim,jdim,lsoil/)) + error = nf90_put_var( ncid, varid_stc, dum3d) + call netcdf_err(error, 'writing stc_inc record' ) + + dum3d = reshape(slc_inc_tmp, (/idim,jdim,lsoil/)) + error = nf90_put_var( ncid, varid_slc, dum3d) + call netcdf_err(error, 'writing slc_inc record' ) + + endif + if(do_nsst) then error=nf90_inq_varid(ncid, "tref", id_var) @@ -1014,7 +1066,11 @@ END SUBROUTINE READ_GSI_DATA !! @param[out] SLMASK Land-sea mask without ice flag. !! @param[out] ZSOIL Soil layer thickness. !! @param[out] NSST Data structure containing nsst fields. + !! @param[in] SLCINC Liquid soil moisture increments on the cubed-sphere tiles + !! @param[in] STCINC Soil temperature increments on the cubed-sphere tiles !! @author George Gayno NOAA/EMC + !! @author Yuan Xue: add capability to read soil related increments on the cubed + !! sphere tiles from files SUBROUTINE READ_DATA(LSOIL,LENSFC,DO_NSST,INC_FILE,IS_NOAHMP, & TSFFCS,SMCFCS,SWEFCS,STCFCS, & TG3FCS,ZORFCS, & @@ -1025,6 +1081,7 @@ SUBROUTINE READ_DATA(LSOIL,LENSFC,DO_NSST,INC_FILE,IS_NOAHMP, & SIHFCS,SICFCS,SITFCS, & TPRCP,SRFLAG,SNDFCS, & VMNFCS,VMXFCS,SLCFCS, & + SLCINC,STCINC,& SLPFCS,ABSFCS,T2M,Q2M,SLMASK, & ZSOIL,NSST) USE MPI @@ -1054,6 +1111,8 @@ SUBROUTINE READ_DATA(LSOIL,LENSFC,DO_NSST,INC_FILE,IS_NOAHMP, & REAL, OPTIONAL, INTENT(OUT) :: SLCFCS(LENSFC,LSOIL) REAL, OPTIONAL, INTENT(OUT) :: SMCFCS(LENSFC,LSOIL) REAL, OPTIONAL, INTENT(OUT) :: STCFCS(LENSFC,LSOIL) + REAL, OPTIONAL, INTENT(OUT) :: SLCINC(LENSFC,LSOIL) + REAL, OPTIONAL, INTENT(OUT) :: STCINC(LENSFC,LSOIL) REAL(KIND=4), OPTIONAL, INTENT(OUT) :: ZSOIL(LSOIL) TYPE(NSST_DATA), OPTIONAL :: NSST ! intent(out) will crash @@ -1106,7 +1165,7 @@ SUBROUTINE READ_DATA(LSOIL,LENSFC,DO_NSST,INC_FILE,IS_NOAHMP, & ERROR=NF90_INQ_VARID(NCID, "canliqxy", ID_VAR) ERROR2=NF90_INQ_VARID(NCID, "tsnoxy", ID_VAR) IS_NOAHMP=.FALSE. - IF(ERROR == 0 .AND. ERROR2 == 0) THEN + IF(ERROR == 0 .OR. ERROR2 == 0) THEN IS_NOAHMP=.TRUE. print*,"- WILL PROCESS FOR NOAH-MP LSM." ENDIF @@ -1501,6 +1560,23 @@ SUBROUTINE READ_DATA(LSOIL,LENSFC,DO_NSST,INC_FILE,IS_NOAHMP, & STCFCS = RESHAPE(DUMMY3D, (/LENSFC,LSOIL/)) ENDIF + IF (PRESENT(SLCINC)) THEN + ERROR=NF90_INQ_VARID(NCID, "slc_inc", ID_VAR) + CALL NETCDF_ERR(ERROR, 'READING slc_inc ID' ) + ERROR=NF90_GET_VAR(NCID, ID_VAR, dummy3d) + CALL NETCDF_ERR(ERROR, 'READING slc increments' ) + SLCINC = RESHAPE(DUMMY3D, (/LENSFC,LSOIL/)) + ENDIF + + IF (PRESENT(STCINC)) THEN + ERROR=NF90_INQ_VARID(NCID, "stc_inc", ID_VAR) + CALL NETCDF_ERR(ERROR, 'READING stc_inc ID' ) + ERROR=NF90_GET_VAR(NCID, ID_VAR, dummy3d) + CALL NETCDF_ERR(ERROR, 'READING stc increments' ) + STCINC = RESHAPE(DUMMY3D, (/LENSFC,LSOIL/)) + ENDIF + + DEALLOCATE(DUMMY3D) ! cloud fields not in warm restart files. set to zero? diff --git a/sorc/lsm_routines.fd/noah.fd/set_soilveg_snippet.f90 b/sorc/lsm_routines.fd/noah.fd/set_soilveg_snippet.f90 index 0ba432ba1..a0743d847 100644 --- a/sorc/lsm_routines.fd/noah.fd/set_soilveg_snippet.f90 +++ b/sorc/lsm_routines.fd/noah.fd/set_soilveg_snippet.f90 @@ -1,9 +1,13 @@ !> @file !> @brief Routine to set Noah LSM soil and veg params needed for sflx_snippet !> @author Clara Draper - !> Below was extracted from namelist_soilveg.f and set_soilveg.f !! (couldn't get above to compile for doxygen) +!> @author Yuan Xue: add Noah-MP LSM soil and veg params needed for global_cycle +!> Noah-MP related parameters were extracted from noahmp_table.f +!> isot (soil type) = 1: STATSGO must be selected if NoahMP is used +!> ivet (vegetation type) = 1: IBGP is used by UFS offline Land DA for Noah-MP +!> as of 07/13/2023 module set_soilveg_snippet_mod @@ -11,19 +15,20 @@ module set_soilveg_snippet_mod private - public set_soilveg + public set_soilveg_noah + public set_soilveg_noahmp contains !> This subroutine initializes soil and vegetation -!! parameters needed in global_cycle/land_increment.f90 +!! parameters needed in global_cycle/land_increment.f90 for noah !! @param[in] isot Soil type !! @param[in] ivet Vegetation type !! @param[out] maxsmc Maximum soil moisture for each soil type !! @param[out] bb B exponent for each soil type !! @param[out] satpsi Saturated matric potential for each soil type !! @param[out] iret Return integer -subroutine set_soilveg(isot,ivet, maxsmc, bb, satpsi, iret) +subroutine set_soilveg_noah(isot,ivet, maxsmc, bb, satpsi, iret) implicit none integer, intent(in) :: isot,ivet @@ -85,6 +90,51 @@ subroutine set_soilveg(isot,ivet, maxsmc, bb, satpsi, iret) iret = 0 -end subroutine set_soilveg +end subroutine set_soilveg_noah + +!> This subroutine initializes soil and vegetation +!! parameters needed in global_cycle/land_increment.f90 for noah-mp +!! @param[in] isot Soil type +!! @param[in] ivet Vegetation type +!! @param[out] maxsmc Maximum soil moisture for each soil type +!! @param[out] bb B exponent for each soil type +!! @param[out] satpsi Saturated matric potential for each soil type +!! @param[out] iret Return integer +subroutine set_soilveg_noahmp(isot,ivet, maxsmc, bb, satpsi,iret) + + implicit none + + integer, intent(in) :: isot,ivet !ivet is *not* used for now + real, dimension(30), intent(out) :: maxsmc, bb, satpsi + integer, intent(out) :: iret + + if (isot .eq. 1) then + +! set soil-dependent params (STATSGO is the only option for UFS, 07/13/2023) + maxsmc= (/0.339, 0.421, 0.434, 0.476, 0.484,& + & 0.439, 0.404, 0.464, 0.465, 0.406, 0.468, 0.468, & + & 0.439, 1.000, 0.200, 0.421, 0.468, 0.200, & + & 0.339, 0.339, 0.000, 0.000, 0.000, 0.000, & + & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/) + bb= (/2.79, 4.26, 4.74, 5.33, 3.86, 5.25,& + & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, & + & 5.25, 0.0, 2.79, 4.26, 11.55, 2.79, & + & 2.79, 0.00, 0.00, 0.00, 0.00, 0.00, & + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) + satpsi= (/0.069, 0.036, 0.141, 0.759, 0.955, & + & 0.355, 0.135, 0.617, 0.263, 0.098, 0.324, 0.468, & + & 0.355, 0.00, 0.069, 0.036, 0.468, 0.069, & + & 0.069, 0.00, 0.00, 0.00, 0.00, 0.00, & + & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/) + + else + print*, 'For Noah-MP, set_soilveg is not supported for soil type ', isot + iret = -1 + return + + endif + + iret = 0 +end subroutine set_soilveg_noahmp end module set_soilveg_snippet_mod diff --git a/tests/global_cycle/CMakeLists.txt b/tests/global_cycle/CMakeLists.txt index 2427551f2..ba297034e 100644 --- a/tests/global_cycle/CMakeLists.txt +++ b/tests/global_cycle/CMakeLists.txt @@ -16,10 +16,28 @@ include_directories(${PROJECT_SOURCE_DIR}) execute_process( COMMAND ${CMAKE_COMMAND} -E copy ${CMAKE_CURRENT_SOURCE_DIR}/LSanSuppress.supp ${CMAKE_CURRENT_BINARY_DIR}/LSanSuppress.supp) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.001 ${CMAKE_CURRENT_BINARY_DIR}/xainc.001) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.002 ${CMAKE_CURRENT_BINARY_DIR}/xainc.002) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.003 ${CMAKE_CURRENT_BINARY_DIR}/xainc.003) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.004 ${CMAKE_CURRENT_BINARY_DIR}/xainc.004) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.005 ${CMAKE_CURRENT_BINARY_DIR}/xainc.005) +execute_process( COMMAND ${CMAKE_COMMAND} -E copy + ${CMAKE_CURRENT_SOURCE_DIR}/data/xainc.006 ${CMAKE_CURRENT_BINARY_DIR}/xainc.006) add_executable(ftst_land_increments ftst_land_increments.F90) target_link_libraries(ftst_land_increments global_cycle_lib) +add_executable(ftst_read_increments ftst_read_increments.F90) +target_link_libraries(ftst_read_increments global_cycle_lib) add_mpi_test(global_cycle-ftst_land_increments EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_land_increments NUMPROCS 1 TIMEOUT 60) +add_mpi_test(global_cycle-ftst_read_increments + EXECUTABLE ${CMAKE_CURRENT_BINARY_DIR}/ftst_read_increments + NUMPROCS 6 TIMEOUT 60) + diff --git a/tests/global_cycle/ftst_land_increments.F90 b/tests/global_cycle/ftst_land_increments.F90 index 8519e178b..d4ebded2e 100644 --- a/tests/global_cycle/ftst_land_increments.F90 +++ b/tests/global_cycle/ftst_land_increments.F90 @@ -1,8 +1,12 @@ program ftst_land_increments -! Test "apply_land_da_adjustments" using sample points. +! Test "apply_land_da_adjustments_soil_temp" using sample points. ! ! author: George Gayno (george.gayno@noaa.gov) +! Yuan Xue: add a total of four testing points to test out +! the newly added frozen soil ice fraction calculations +! under different scenarios after ingesting +! soil temp increments use mpi use land_increments @@ -14,14 +18,13 @@ program ftst_land_increments real, parameter :: EPSILON=0.001 - real(kind=4), allocatable :: zsoil(:) real, allocatable :: rsoiltype(:) integer, allocatable :: mask(:) real, allocatable :: stc_bck(:,:) real, allocatable :: smc_anl(:,:) real, allocatable :: slc_anl(:,:) real, allocatable :: stc_anl(:,:) - integer, allocatable :: stc_updated(:), slc_updated(:) + integer, allocatable :: stc_updated(:) call mpi_init(ierr) @@ -32,9 +35,8 @@ program ftst_land_increments ivegsrc = 1 ! IGBP vegetation type. lsoil = 4 ! Number of soil layers. - lensfc= 3 ! Number of test points. + lensfc= 4 ! Number of test points. - allocate(zsoil(lsoil)) allocate(rsoiltype(lensfc)) ! Soil type. allocate(mask(lensfc)) ! Land mask allocate(stc_bck(lensfc,lsoil)) ! Background soil temperature (K). @@ -42,18 +44,13 @@ program ftst_land_increments allocate(slc_anl(lensfc,lsoil)) ! Analyzed liquid soil moisture. allocate(stc_anl(lensfc,lsoil)) ! Analyzed soil temperature (K). allocate(stc_updated(lensfc)) - allocate(slc_updated(lensfc)) - zsoil(1) = -0.1 - zsoil(2) = -0.4 - zsoil(3) = -1.0 - zsoil(4) = -2.0 - slc_updated=1 stc_updated=1 ! Point 1 is above freezing before the adjustment ! and above freezing after the adjustment. Therefore, -! the increments to STC and SLC will be retained. +! the increments to STC will be retained. +! temp: unfrozen ==> unfrozen rsoiltype(1) = 5. mask(1) = 1 @@ -61,55 +58,77 @@ program ftst_land_increments smc_anl(1,:) = .25 slc_anl(1,:) = .25 - stc_anl(1,1:3) = 281.0 ! DA only updates 3 layers + stc_anl(1,1:2) = 281.0 ! DA only updates 2 layers ! Point 2 is below freezing before the adjustment ! and above freezing after the adjustment. Therefore, -! the increment to STC will be removed, and SMC / SLC -! are unchanged. +! all soil ice will be melted +! temp: frozen ==> unfrozen rsoiltype(2) = 5. mask(2) = 1 - stc_bck(2,:) = 270.0 + stc_bck(2,:) = 271.0 smc_anl(2,:) = .25 slc_anl(2,:) = .20 - stc_anl(2,1:3) = 274.0 + stc_anl(2,1:2) = 275.0 -! Point 3 freezes. Therefore, -! the increment to STC will be removed, and SMC / SLC -! are unchanged. +! Point 3 freezes before and after the adjustment. Therefore, +! SLC will be recomputed, SMC remained unchanged +! temp: frozen ==> frozen rsoiltype(3) = 5. mask(3) = 1 - stc_bck(3,:) = 274.0 + stc_bck(3,:) = 272.0 smc_anl(3,:) = .25 - slc_anl(3,:) = .25 - stc_anl(3,1:3) = 271.0 + slc_anl(3,:) = .20 + stc_anl(3,1:2) = 271.0 - call apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, & +! Point 4 is above freezing before and below freezing after the adjustment +! Therfore, SLC will be recomputed, SMC remained unchanged +! temp: unfrozen ==> frozen + + rsoiltype(4) = 5. + mask(4) = 1 + stc_bck(4,:) = 280.0 + + smc_anl(4,:) = .25 + slc_anl(4,:) = .25 + stc_anl(4,1:2) = 271.0 + + + call apply_land_da_adjustments_soil_temp(lsm, isot, ivegsrc,lensfc, & lsoil, rsoiltype, mask, stc_bck, stc_anl, smc_anl, slc_anl, & - stc_updated, slc_updated,zsoil) + stc_updated) - do l = 1, 3 - if (abs(smc_anl(1,l) - 0.25) > EPSILON) stop 2 - if (abs(slc_anl(1,l) - 0.25) > EPSILON) stop 4 - if (abs(stc_anl(1,l) - 281.) > EPSILON) stop 3 + do l = 1,2 + !print*, 'Testing Point#1' + if (abs(smc_anl(1,l) - 0.25) > EPSILON) stop 2 + if (abs(slc_anl(1,l) - 0.25) > EPSILON) stop 3 + if (abs(stc_anl(1,l) - 281.) > EPSILON) stop 4 enddo - do l = 1, 3 - if (abs(smc_anl(2,l) - 0.25) > EPSILON) stop 6 - if (abs(slc_anl(2,l) - 0.20) > EPSILON) stop 8 - if (abs(stc_anl(2,l) - 270.) > EPSILON) stop 5 + do l = 1,2 + !print*, 'Testing Point#2' + if (abs(smc_anl(2,l) - 0.25) > EPSILON) stop 5 + if (abs(slc_anl(2,l) - 0.25) > EPSILON) stop 6 + if (abs(stc_anl(2,l) - 275.) > EPSILON) stop 7 enddo - do l = 1, 3 - if (abs(smc_anl(3,l) - 0.25) > EPSILON) stop 10 - if (abs(slc_anl(3,l) - 0.25) > EPSILON) stop 12 - if (abs(stc_anl(3,l) - 274.) > EPSILON) stop 11 + do l = 1,2 + !print*, 'Testing Point#3' + if (abs(smc_anl(3,l) - 0.25) > EPSILON) stop 8 + if (abs(slc_anl(3,l) - 0.112) > EPSILON) stop 9 + if (abs(stc_anl(3,l) - 271.) > EPSILON) stop 10 enddo + do l = 1,2 + !print*, 'Testing Point#4' + if (abs(smc_anl(4,l) - 0.25) > EPSILON) stop 11 + if (abs(slc_anl(4,l) - 0.112) > EPSILON) stop 12 + if (abs(stc_anl(4,l) - 271.) > EPSILON) stop 13 + enddo call mpi_finalize(ierr) deallocate(rsoiltype,stc_bck,smc_anl,slc_anl,stc_anl,mask) diff --git a/tests/global_cycle/ftst_read_increments.F90 b/tests/global_cycle/ftst_read_increments.F90 new file mode 100644 index 000000000..555489976 --- /dev/null +++ b/tests/global_cycle/ftst_read_increments.F90 @@ -0,0 +1,139 @@ +! Unit test for global_cycle routine "read_data". +! +! Reads a sample increment file from file xainc.$NNN +! NOTE: $NNN corresponds to (mpi rank + 1) +! Total number of processes= 6 (corresponds to 6 tiles) +! Each process is assigned one increment file on one tile to read +! If any portion of the variable does not match expected values, +! the test fails. +! +! Author: Yuan Xue, 07/17/2023 + + program read_increments + + use read_write_data, only : read_data + use mpi + + implicit none + + integer:: ierr,my_rank + integer:: lsoil, l + integer:: lensfc, ij_input + + integer, parameter:: NUM_VALUES=4 + real, parameter :: EPSILON=0.0001 + + real, allocatable:: SLCINC(:,:),STCINC(:,:) + real:: stc_inc_expected_values_tile1(NUM_VALUES) + real:: stc_inc_expected_values_tile2(NUM_VALUES) + real:: stc_inc_expected_values_tile3(NUM_VALUES) + real:: stc_inc_expected_values_tile4(NUM_VALUES) + real:: stc_inc_expected_values_tile5(NUM_VALUES) + real:: stc_inc_expected_values_tile6(NUM_VALUES) + real:: slc_inc_expected_values_tile1(NUM_VALUES) + real:: slc_inc_expected_values_tile2(NUM_VALUES) + real:: slc_inc_expected_values_tile3(NUM_VALUES) + real:: slc_inc_expected_values_tile4(NUM_VALUES) + real:: slc_inc_expected_values_tile5(NUM_VALUES) + real:: slc_inc_expected_values_tile6(NUM_VALUES) + + !expected values were extracted from MATLAB, which directly reads in xainc file + !each tile is examined separately here + data stc_inc_expected_values_tile1 / -0.6302, -0.1116, 0.0341, 0.0 / + data stc_inc_expected_values_tile2 / 0.0825, 0.0071, -0.0255, 0.0 / + data stc_inc_expected_values_tile3 / 0.2070, 0.0608, 0.0001, 0.0 / + data stc_inc_expected_values_tile4 / 0.0, 0.0, 0.0, 0.0 / + data stc_inc_expected_values_tile5 / -0.1031, -0.0386, -0.0356, 0.0/ + data stc_inc_expected_values_tile6 / 0.0, 0.0, 0.0, 0.0 / + data slc_inc_expected_values_tile1 / -0.0007285, 0.0000055, 0.0000003, 0.0 / + data slc_inc_expected_values_tile2 / -0.0006, -0.0059, -0.0087, 0.0 / + data slc_inc_expected_values_tile3 / -0.0015, 0.0030, -0.0007, 0.0 / + data slc_inc_expected_values_tile4 / 0.0, 0.0, 0.0, 0.0 / + data slc_inc_expected_values_tile5 / 0.0014, 0.0012, 0.0006, 0.0 / + data slc_inc_expected_values_tile6 / 0.0, 0.0, 0.0, 0.0 / + + call mpi_init(ierr) + call MPI_Comm_rank(MPI_COMM_WORLD, my_rank, ierr) + + lsoil=4 + lensfc=36864 + ij_input = 28541 + + allocate(SLCINC(lensfc,lsoil)) + allocate(STCINC(lensfc,lsoil)) + + if (my_rank .eq. 0) print*,"Starting test of global_cycle routine read_data." + if (my_rank .eq. 0) print*,"Call routine read_data" + + call read_data(lsoil,lensfc,.false.,.true.,SLCINC=SLCINC,STCINC=STCINC) + + if (my_rank .eq. 0) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile1(l))& + > EPSILON) stop 20 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile1(l))& + > EPSILON) stop 40 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + if (my_rank .eq. 1) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile2(l))& + > EPSILON) stop 21 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile2(l))& + > EPSILON) stop 41 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + if (my_rank .eq. 2) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile3(l))& + > EPSILON) stop 22 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile3(l))& + > EPSILON) stop 42 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + if (my_rank .eq. 3) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile4(l))& + > EPSILON) stop 23 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile4(l))& + > EPSILON) stop 43 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + if (my_rank .eq. 4) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile5(l))& + > EPSILON) stop 24 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile5(l))& + > EPSILON) stop 44 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + if (my_rank .eq. 5) then + do l = 1,4 + if (abs(STCINC(ij_input,l) - stc_inc_expected_values_tile6(l))& + > EPSILON) stop 25 + if (abs(SLCINC(ij_input,l) - slc_inc_expected_values_tile6(l))& + > EPSILON) stop 45 + enddo + print*, "tile#", my_rank+1, "reads OK" + endif + + call MPI_Barrier(MPI_COMM_WORLD, ierr) + + if (my_rank .eq. 0) print*, "ALL is OK" + if (my_rank .eq. 0) print*, "SUCCESS!" + + deallocate(SLCINC,STCINC) + call mpi_finalize(ierr) + + end program read_increments +