diff --git a/lis/core/LIS_histDataMod.F90 b/lis/core/LIS_histDataMod.F90 index 30b7f15ad..9386e3262 100644 --- a/lis/core/LIS_histDataMod.F90 +++ b/lis/core/LIS_histDataMod.F90 @@ -487,6 +487,26 @@ module LIS_histDataMod public :: LIS_MOC_SNOWSURFACEQ !public :: LIS_MOC_SNOWWIND_DIR + ! MMF SW + !public :: LIS_MOC_SMCWTD + public :: LIS_MOC_QRF + public :: LIS_MOC_QRFS + public :: LIS_MOC_QSPRING + public :: LIS_MOC_QSPRINGS + public :: LIS_MOC_QSLAT + public :: LIS_MOC_RECH + public :: LIS_MOC_WTD + public :: LIS_MOC_DEEPRECH + !integer :: LIS_MOC_SMCWTD = -9999 + integer :: LIS_MOC_QRF = -9999 + integer :: LIS_MOC_QRFS = -9999 + integer :: LIS_MOC_DEEPRECH = -9999 + integer :: LIS_MOC_QSPRING = -9999 + integer :: LIS_MOC_QSPRINGS = -9999 + integer :: LIS_MOC_QSLAT = -9999 + integer :: LIS_MOC_RECH = -9999 + integer :: LIS_MOC_WTD = -9999 + ! SnowModel outputs: public :: LIS_MOC_SWE_SM public :: LIS_MOC_SNOWDEPTH_SM @@ -1197,6 +1217,115 @@ subroutine LIS_histDataInit(n, ntiles) !------------------------------------------------------------------------- ! read the meta data attributes for each variable !------------------------------------------------------------------------- + !!! MMF scheme + !call ESMF_ConfigFindLabel(modelSpecConfig,"SmcWtd:",rc=rc) + !call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + ! "SmcWtd",& + ! "soil_moisture_between_bottom_of_the_soil_and_the_water_table",& + ! "soil moisture between bottom of the soil and the water table",rc) + !if ( rc == 1 ) then + ! call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_SMCWTD,& + ! LIS_histData(n)%head_lsm_list,& + ! n,1,ntiles,(/"m3/m3"/),1,(/"-"/),1,112,0,& + ! model_patch=.true.) + !endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Qrf:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Qrf",& + "groundwater_river_water_flux",& + "groundwater river water flux",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_QRF,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Qrfs:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Qrfs",& + "cumulated_groundwater_river_water_flux",& + "cumulated groundwater river water flux",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_QRFS,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Qspring:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Qspring",& + "water_springing_at_the_surface_from_groundwater_convergence_in_the_column",& + "water springing at the surface from groundwater convergence in the column",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_QSPRING,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Qsprings:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Qsprings",& + "accumulated_water_springing_at_the_surface_from_groundwater_convergence_in_the_column",& + "accumulated water springing at the surface from groundwater convergence in the column",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_QSPRINGS,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Qslat:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Qslat",& + "accumulated_lateral_flow", & + "accumulated lateral flow",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_QSLAT,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Rech:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Rech",& + "accumulated_recharge", & + "accumulated recharge",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_RECH,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"DeepRech:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "DeepRech",& + "accumulated_recharge_when_deep", & + "accumulated recharge when deep",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_DEEPRECH,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"mm"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + + call ESMF_ConfigFindLabel(modelSpecConfig,"Wtd:",rc=rc) + call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & + "Wtd",& + "accumulated_groundwater_river_water_flux",& + "accumulated groundwater river water flux",rc) + if ( rc == 1 ) then + call register_dataEntry(LIS_MOC_LSM_COUNT,LIS_MOC_WTD,& + LIS_histData(n)%head_lsm_list,& + n,1,ntiles,(/"m"/),1,(/"-"/),1,112,0,& + model_patch=.true.) + endif + !!! JULES call ESMF_ConfigFindLabel(modelSpecConfig,"SnowSoot:",rc=rc) call get_moc_attributes(modelSpecConfig, LIS_histData(n)%head_lsm_list, & diff --git a/lis/core/LIS_historyMod.F90 b/lis/core/LIS_historyMod.F90 index eafc2f03a..d9291a060 100644 --- a/lis/core/LIS_historyMod.F90 +++ b/lis/core/LIS_historyMod.F90 @@ -118,6 +118,7 @@ module LIS_historyMod public :: LIS_gather_1dgrid_to_2dgrid public :: LIS_scatter_global_to_local_grid public :: LIS_gather_2d_local_to_global + public :: LIS_gather_masterproc_2d_local_to_global !EOP !BOP @@ -396,6 +397,24 @@ module LIS_historyMod !EOP end interface +!BOP +! +! !ROUTINE: LIS_gather_masterproc_2d_local_to_global +! \label{LIS_gather_masterproc_2d_local_to_global} +! +! !INTERFACE: + interface LIS_gather_masterproc_2d_local_to_global +! !PRIVATE MEMBER FUNCTIONS: + module procedure gather_masterproc_2d_local_to_global_int + module procedure gather_masterproc_2d_local_to_global_real +! +! !DESCRIPTION: +! This interface provides routines for gathering 2d local arrays +! to 2d global array on masterproc. +! +!EOP + end interface + contains !BOP @@ -9803,7 +9822,7 @@ end subroutine LIS_gather_1dgrid_to_2dgrid ! 30 Jan 2009: Sujay Kumar, Initial Code ! ! !INTERFACE: -subroutine LIS_scatter_global_to_local_grid(n, gtmp,ltmp) +subroutine LIS_scatter_global_to_local_grid(n, gtmp, ltmp) ! !USES: ! !ARGUMENTS: @@ -9815,7 +9834,7 @@ subroutine LIS_scatter_global_to_local_grid(n, gtmp,ltmp) ! ! !DESCRIPTION: -! This routine gathers a gridded 1d array into a 2d global array. +! This routine scatters a gridded 2d global array into a 2d local array. ! ! This process accounts for the halo. ! @@ -9824,15 +9843,16 @@ subroutine LIS_scatter_global_to_local_grid(n, gtmp,ltmp) ! \item [n] ! index of the current nest ! \item [gtmp] -! return array for the gridded output data -! \item [var] +! array for the gridded global input data +! \item [ltmp] +! array for the gridded local output data ! output data to process ! \end{description} !EOP real, allocatable :: gtmp1d(:) real, allocatable :: gtmp2d(:,:) - integer :: i,c,r,m,t,l + integer :: c,r integer :: gid integer :: ierr @@ -10002,4 +10022,194 @@ subroutine LIS_gather_2d_local_to_global(n, lvar, gvar) deallocate(var1) end subroutine LIS_gather_2d_local_to_global +subroutine gather_masterproc_2d_local_to_global_int(n, lvar, gvar) +! !USES: + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + integer, intent(in) :: lvar(LIS_rc%lnc(n), LIS_rc%lnr(n)) + integer, intent(out) :: gvar(LIS_rc%gnc(n), LIS_rc%gnr(n)) + +! !DESCRIPTION: +! This routine gathers local (sub-domain) gridded 2d arrays into a +! 2d global array on the masterproc. +! +! The arguments are: +! \begin{description} +! \item [n] +! index of the current nest +! \item [lvar] +! local (per-process) 2d array +! \item [gvar] +! global 2d array +! \end{description} +!EOP + + integer, allocatable :: var1(:) + integer, allocatable :: gtmp1(:) + integer :: c, r, l + integer :: c1, r1 + integer :: count1 + integer :: ierr + integer :: deltas + + allocate(var1(LIS_rc%lnc_red(n)*LIS_rc%lnr_red(n))) + allocate(gtmp1(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + + count1 = 1 + do r = LIS_nss_halo_ind(n,LIS_localPet+1), LIS_nse_halo_ind(n,LIS_localPet+1) + do c = LIS_ews_halo_ind(n,LIS_localPet+1), LIS_ewe_halo_ind(n,LIS_localPet+1) + if ( r .ge. LIS_nss_ind(n,LIS_localPet+1) .and. & + r .le. LIS_nse_ind(n,LIS_localPet+1) .and. & + c .ge. LIS_ews_ind(n,LIS_localPet+1) .and. & + c .le. LIS_ewe_ind(n,LIS_localPet+1) ) then + r1 = r - LIS_nss_halo_ind(n,LIS_localPet+1) + 1 + c1 = c - LIS_ews_halo_ind(n,LIS_localPet+1) + 1 + var1(count1) = lvar(c1,r1) + count1 = count1 + 1 + endif + enddo + enddo + + +#if (defined SPMD) + deltas = LIS_deltas(n,LIS_localPet) + call MPI_GATHERV(var1, deltas, MPI_INTEGER, & + gtmp1, LIS_deltas(n,:), LIS_offsets(n,:), MPI_INTEGER, & + 0, LIS_mpi_comm, ierr) +#else + gtmp1 = var1 +#endif + +! Looping over halo is unnecessary here because gtmp does not contain +! any halo points. I am keeping this here for reference. +! if ( LIS_masterproc ) then +! count1 = 1 +! do l = 1, LIS_npes +! do r = LIS_nss_halo_ind(n,l), LIS_nse_halo_ind(n,l) +! do c = LIS_ews_halo_ind(n,l), LIS_ewe_halo_ind(n,l) +! if ( r .ge. LIS_nss_ind(n,l) .and. & +! r .le. LIS_nse_ind(n,l) .and. & +! c .ge. LIS_ews_ind(n,l) .and. & +! c .le. LIS_ewe_ind(n,l) ) then +! gvar(c,r) = gtmp1(count1) +! count1 = count1 + 1 +! endif +! enddo +! enddo +! enddo +! endif + if ( LIS_masterproc ) then + count1 = 1 + do l = 1, LIS_npes + do r = LIS_nss_ind(n,l), LIS_nse_ind(n,l) + do c = LIS_ews_ind(n,l), LIS_ewe_ind(n,l) + gvar(c,r) = gtmp1(count1) + count1 = count1 + 1 + enddo + enddo + enddo + endif + + + deallocate(gtmp1) + deallocate(var1) +end subroutine gather_masterproc_2d_local_to_global_int + +subroutine gather_masterproc_2d_local_to_global_real(n, lvar, gvar) +! !USES: + + implicit none +! !ARGUMENTS: + integer, intent(in) :: n + real, intent(in) :: lvar(LIS_rc%lnc(n), LIS_rc%lnr(n)) + real, intent(out) :: gvar(LIS_rc%gnc(n), LIS_rc%gnr(n)) + +! !DESCRIPTION: +! This routine gathers local (sub-domain) gridded 2d arrays into a +! 2d global array on the masterproc. +! +! The arguments are: +! \begin{description} +! \item [n] +! index of the current nest +! \item [lvar] +! local (per-process) 2d array +! \item [gvar] +! global 2d array +! \end{description} +!EOP + + real, allocatable :: var1(:) + real, allocatable :: gtmp1(:) + integer :: c, r, l + integer :: c1, r1 + integer :: count1 + integer :: ierr + integer :: deltas + + allocate(var1(LIS_rc%lnc_red(n)*LIS_rc%lnr_red(n))) + allocate(gtmp1(LIS_rc%gnc(n)*LIS_rc%gnr(n))) + + count1 = 1 + do r = LIS_nss_halo_ind(n,LIS_localPet+1), LIS_nse_halo_ind(n,LIS_localPet+1) + do c = LIS_ews_halo_ind(n,LIS_localPet+1), LIS_ewe_halo_ind(n,LIS_localPet+1) + if ( r .ge. LIS_nss_ind(n,LIS_localPet+1) .and. & + r .le. LIS_nse_ind(n,LIS_localPet+1) .and. & + c .ge. LIS_ews_ind(n,LIS_localPet+1) .and. & + c .le. LIS_ewe_ind(n,LIS_localPet+1) ) then + r1 = r - LIS_nss_halo_ind(n,LIS_localPet+1) + 1 + c1 = c - LIS_ews_halo_ind(n,LIS_localPet+1) + 1 + var1(count1) = lvar(c1,r1) + count1 = count1 + 1 + endif + enddo + enddo + + +#if (defined SPMD) + deltas = LIS_deltas(n,LIS_localPet) + call MPI_GATHERV(var1, deltas, MPI_REAL, & + gtmp1, LIS_deltas(n,:), LIS_offsets(n,:), MPI_REAL, & + 0, LIS_mpi_comm, ierr) +#else + gtmp1 = var1 +#endif + +! Looping over halo is unnecessary here because gtmp does not contain +! any halo points. I am keeping this here for reference. +! if ( LIS_masterproc ) then +! count1 = 1 +! do l = 1, LIS_npes +! do r = LIS_nss_halo_ind(n,l), LIS_nse_halo_ind(n,l) +! do c = LIS_ews_halo_ind(n,l), LIS_ewe_halo_ind(n,l) +! if ( r .ge. LIS_nss_ind(n,l) .and. & +! r .le. LIS_nse_ind(n,l) .and. & +! c .ge. LIS_ews_ind(n,l) .and. & +! c .le. LIS_ewe_ind(n,l) ) then +! gvar(c,r) = gtmp1(count1) +! count1 = count1 + 1 +! endif +! enddo +! enddo +! enddo +! endif + if ( LIS_masterproc ) then + count1 = 1 + do l = 1, LIS_npes + do r = LIS_nss_ind(n,l), LIS_nse_ind(n,l) + do c = LIS_ews_ind(n,l), LIS_ewe_ind(n,l) + gvar(c,r) = gtmp1(count1) + count1 = count1 + 1 + enddo + enddo + enddo + endif + + + deallocate(gtmp1) + deallocate(var1) +end subroutine gather_masterproc_2d_local_to_global_real + end module LIS_historyMod diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_coldstart.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_coldstart.F90 index e6c64a28f..61eb38918 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_coldstart.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_coldstart.F90 @@ -68,7 +68,7 @@ subroutine NoahMP401_coldstart(mtype) real, dimension(4) :: DZS ! Thickness of the soil layers [m] real :: dx, dy real, dimension( 1, 1 ) :: msftx, msfty - real :: wtddt, dtbl + real :: wtddt integer :: stepwtd real, dimension( 1, 1 ) :: & @@ -178,6 +178,8 @@ subroutine NoahMP401_coldstart(mtype) real, dimension(1,60,1) :: gecros_state ! Optional gecros crop + + ! PLANTING(1,1) = 126 ! default planting date ! HARVEST(1,1) = 290 ! default harvest date ! SEASON_GDD(1,1) = 1605 ! default total seasonal growing degree days @@ -343,13 +345,12 @@ subroutine NoahMP401_coldstart(mtype) rechxy(1,1) = 0.0 deeprechxy(1,1) = 0.0 areaxy(1,1) = 100. - dx = 10.0 - dy = 10.0 - msftx(1,1) = 0.0 - msfty(1,1) = 0.0 - wtddt = 0.0 + dx = 10.0 ! SW, not used + dy = 10.0 ! SW, not used + msftx(1,1) = 1.0 + msfty(1,1) = 1.0 + wtddt = NOAHMP401_struc(n)%ts/60.0 ! wtddt in minute? stepwtd = 0 - dtbl = 0.0 qrfsxy(1,1) = 0.0 qslatxy(1,1) = 0.0 fdepthxy(1,1) = 0.0 @@ -384,7 +385,7 @@ subroutine NoahMP401_coldstart(mtype) ims,ime, jms,jme, kms,kme, & ! memory its,ite, jts,jte, kts,kte & ! tile ,smoiseqxy,smcwtdxy ,rechxy ,deeprechxy, areaxy ,dx, dy, msftx, msfty,& - wtddt ,stepwtd ,dtbl ,qrfsxy ,qspringsxy ,qslatxy, & + wtddt ,stepwtd ,NOAHMP401_struc(n)%ts ,qrfsxy ,qspringsxy ,qslatxy, & fdepthxy ,HT ,riverbedxy ,eqzwt ,rivercondxy ,pexpxy, & rechclim ,gecros_state & ) @@ -452,6 +453,11 @@ subroutine NoahMP401_coldstart(mtype) enddo ! t=1,1 endif ! coldstart + + !!! MMF initialization + if(NOAHMP401_struc(n)%run_opt==5) then + call mmf_start(n) + endif deallocate(zsnso) deallocate(tsnow) @@ -477,4 +483,6 @@ subroutine NoahMP401_coldstart(mtype) write(LIS_logunit,*) "[INFO] NoahMP401_coldstart -- ", & "Using the specified start time ", LIS_rc%time enddo ! nnest + end subroutine NoahMP401_coldstart + diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_finalize.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_finalize.F90 index 9448f5068..a9bc0af19 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_finalize.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_finalize.F90 @@ -61,7 +61,32 @@ subroutine NoahMP401_finalize() deallocate(NOAHMP401_struc(n)%init_smc) deallocate(NOAHMP401_struc(n)%init_tslb) deallocate(NOAHMP401_struc(n)%init_gecros_state) - + ! SW, MMF variables, 10/20/2021 + if(NOAHMP401_struc(n)%run_opt==5) then + deallocate(NOAHMP401_struc(n)%smoiseq) + deallocate(NOAHMP401_struc(n)%smcwtd) + deallocate(NOAHMP401_struc(n)%deeprech) + deallocate(NOAHMP401_struc(n)%rech) + deallocate(NOAHMP401_struc(n)%wtd) + deallocate(NOAHMP401_struc(n)%qslat) + deallocate(NOAHMP401_struc(n)%qrfs) + deallocate(NOAHMP401_struc(n)%qsprings) + deallocate(NOAHMP401_struc(n)%fdepth) + deallocate(NOAHMP401_struc(n)%area) + deallocate(NOAHMP401_struc(n)%topo) + deallocate(NOAHMP401_struc(n)%eqwtd) + deallocate(NOAHMP401_struc(n)%riverbed) + deallocate(NOAHMP401_struc(n)%rivercond) + deallocate(NOAHMP401_struc(n)%qrf) + deallocate(NOAHMP401_struc(n)%qspring) + deallocate(NOAHMP401_struc(n)%rechclim) + deallocate(NOAHMP401_struc(n)%rct_idx) + deallocate(NOAHMP401_struc(n)%soil3d) + deallocate(NOAHMP401_struc(n)%soil2d) + deallocate(NOAHMP401_struc(n)%vege3d) + deallocate(NOAHMP401_struc(n)%vege2d) + endif + end do ! nest loop deallocate(NOAHMP401_struc) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 index d38c5c959..a2b96ea02 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 @@ -210,6 +210,29 @@ module NoahMP401_lsmMod integer :: pedo_opt integer :: crop_opt integer :: urban_opt + integer :: row_min, row_max, col_min, col_max ! SW, for MMF + real, pointer :: smoiseq(:,:,:) + real, pointer :: smcwtd(:,:) !MMF state + real, pointer :: deeprech(:,:) !MMF state + real, pointer :: rech(:,:) !MMF state + real, pointer :: wtd(:,:) !MMF state, SW, 09142021 + real, pointer :: qslat(:,:) !MMF state, SW + real, pointer :: qrfs(:,:) !MMF state, SW + real, pointer :: qsprings(:,:) !MMF state, SW + real, pointer :: fdepth(:,:) !SW for MMF 09142021 + real, pointer :: area(:,:) !SW for MMF + real, pointer :: topo(:,:) !SW for MMF + real, pointer :: eqwtd(:,:) !SW for MMF + real, pointer :: riverbed(:,:) !SW for MMF + real, pointer :: rivercond(:,:) !SW for MMF + real, pointer :: qrf(:,:) !MMF output SW for MMF + real, pointer :: qspring(:,:) !MMF output, SW + real, pointer :: rechclim(:,:) !SW + real, pointer :: soil3d(:,:,:) !SW, for MMF + integer, pointer :: soil2d(:,:) !SW, for MMF + real, pointer :: vege3d(:,:,:) !SW, for MMF + integer, pointer :: vege2d(:,:) !SW, for MMF + integer,allocatable:: rct_idx(:,:) ! col, row, tile type(NoahMP401dec), pointer :: noahmp401(:) end type NoahMP401_type_dec @@ -225,7 +248,8 @@ module NoahMP401_lsmMod ! !INTERFACE: subroutine NoahMP401_ini() ! !USES: - use LIS_coreMod, only : LIS_rc + !use LIS_coreMod, only : LIS_rc + use LIS_coreMod use LIS_logMod, only : LIS_verify, LIS_logunit use LIS_timeMgrMod, only : LIS_clock, LIS_calendar, & LIS_update_timestep, LIS_registerAlarm @@ -243,10 +267,9 @@ subroutine NoahMP401_ini() ! \end{description} !EOP implicit none - integer :: n, t + integer :: n, t, num_t, col, row integer :: status character*3 :: fnest ! EMK for RHMin - ! allocate memory for nest allocate(NOAHMP401_struc(LIS_rc%nnest)) @@ -277,6 +300,7 @@ subroutine NoahMP401_ini() allocate(NOAHMP401_struc(n)%noahmp401(t)%snowliq(NOAHMP401_struc(n)%nsnow)) allocate(NOAHMP401_struc(n)%noahmp401(t)%smoiseq(NOAHMP401_struc(n)%nsoil)) allocate(NOAHMP401_struc(n)%noahmp401(t)%gecros_state(60)) + allocate(NOAHMP401_struc(n)%noahmp401(t)%relsmc(NOAHMP401_struc(n)%nsoil)) !SW enddo ! initialize forcing variables to zeros @@ -344,6 +368,101 @@ subroutine NoahMP401_ini() LIS_sfmodel_struc(n)%lyrthk(:) = & 100*NOAHMP401_struc(n)%sldpth(:) LIS_sfmodel_struc(n)%ts = NOAHMP401_struc(n)%ts + + !SW, for MMF, 09202021 + if(NOAHMP401_struc(n)%run_opt .eq. 5) then + !NOAHMP401_struc(n)%row_min = huge(1) + !NOAHMP401_struc(n)%row_max = 0 + !NOAHMP401_struc(n)%col_min = huge(1) + !NOAHMP401_struc(n)%col_max = 0 + !do t = 1, LIS_rc%npatch(n, LIS_rc%lsm_index) + ! NOAHMP401_struc(n)%row_min = min(NOAHMP401_struc(n)%row_min, LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row) + ! NOAHMP401_struc(n)%row_max = max(NOAHMP401_struc(n)%row_max, LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row) + ! NOAHMP401_struc(n)%col_min = min(NOAHMP401_struc(n)%col_min, LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col) + ! NOAHMP401_struc(n)%col_max = max(NOAHMP401_struc(n)%col_max, LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col) + !enddo + + + NOAHMP401_struc(n)%row_min = 1 + NOAHMP401_struc(n)%row_max = LIS_rc%lnr(n) + NOAHMP401_struc(n)%col_min = 1 + NOAHMP401_struc(n)%col_max = LIS_rc%lnc(n) + + allocate(NOAHMP401_struc(n)%rct_idx(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + + NOAHMP401_struc(n)%rct_idx(:,:) = LIS_rc%udef + do t = 1, LIS_rc%npatch(n, LIS_rc%lsm_index) + NOAHMP401_struc(n)%rct_idx(LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col, & + LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row) = t + enddo + + !num_t = (NOAHMP401_struc(n)%row_max - NOAHMP401_struc(n)%row_min + 1) * & + ! (NOAHMP401_struc(n)%col_max - NOAHMP401_struc(n)%col_min + 1) + !col = NOAHMP401_struc(n)%col_min + !row = NOAHMP401_struc(n)%row_min + !do t = 1, num_t + ! NOAHMP401_struc(n)%rct_idx(col,row) = 1; + ! col = col + 1; + ! if (col .gt. NOAHMP401_struc(n)%col_max) then + ! col = NOAHMP401_struc(n)%col_min + ! row = row + 1; + ! if (col .gt. NOAHMP401_struc(n)%col_max + 1) then + ! print *, "WARNING: rct_idx array exceeds declared size" + ! end if + ! end if + !enddo + !print*, t + + ! allocate 2D data structure for MMF 10 DEC 2021; Now all col x row TML + allocate(NOAHMP401_struc(n)%smoiseq(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max, & + NOAHMP401_struc(n)%nsoil)) + + allocate(NOAHMP401_struc(n)%smcwtd(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%deeprech(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%rech(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%wtd(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%qslat(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%qrfs(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%qsprings(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%fdepth(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%area(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%topo(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%eqwtd(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%riverbed(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%rivercond(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%qrf(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%qspring(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%rechclim(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + + allocate(NOAHMP401_struc(n)%soil3d(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max, & + LIS_rc%nsoiltypes)) ! 16 is the number of soil categories, hardcoded temporarily + allocate(NOAHMP401_struc(n)%soil2d(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(NOAHMP401_struc(n)%vege3d(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max, & + LIS_rc%nsurfacetypes)) ! 20 is the number of surfacetype categories, hardcoded temporarily + allocate(NOAHMP401_struc(n)%vege2d(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, & + NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + endif enddo end subroutine NoahMP401_ini end module NoahMP401_lsmMod diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_main.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_main.F90 index b87ee3c0b..5c68406d2 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_main.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_main.F90 @@ -36,6 +36,7 @@ subroutine NoahMP401_main(n) use LIS_logMod, only : LIS_logunit, LIS_endrun use LIS_FORC_AttributesMod use NoahMP401_lsmMod + use module_sf_noahmp_groundwater_401, only : WTABLE_mmf_noahmp implicit none ! !ARGUMENTS: @@ -45,7 +46,7 @@ subroutine NoahMP401_main(n) real :: dt real :: lat, lon real :: tempval - integer :: row, col, tid + integer :: row, col, tid, ridx, cidx integer :: year, month, day, hour, minute, second logical :: alarmCheck @@ -246,6 +247,9 @@ subroutine NoahMP401_main(n) real, dimension(1,1) :: tmp_infxs1rt ! extra output for WRF-HYDRO [m] real, dimension(1,1) :: tmp_soldrain1rt ! extra output for WRF-HYDRO [m] + ! TML: Debugging term to print model variables if set to 1. + integer :: tmp_printdebug ! print model output if true + !ag (05Jan2021) real :: tmp_rivsto real :: tmp_fldsto @@ -256,6 +260,55 @@ subroutine NoahMP401_main(n) character*3 :: fnest REAL, PARAMETER:: LVH2O = 2.501000E+6 ! Latent heat for evapo for water + !!!! MMF scheme + integer :: min_row, max_row, min_col, max_col, nsoil + integer :: isurban ! urban type + real :: wtddt ! + real :: ddz ! + real :: totwater, deltat + integer, allocatable :: isltyp(:, :) ! soil type + integer, allocatable :: ivgtyp(:, :) ! vegetation type + real, allocatable :: smoiseq(:, :, :) ! equilibrium soil moisture content [m3/m3] + real, allocatable :: dzs(:) ! thickness of soil layers [m] + real, allocatable :: fdepth(:, :) ! efolding depth for transmissivity (m) + real, allocatable :: area(:, :) ! grid cell area in m^2 ? + real, allocatable :: topo(:, :) + real, allocatable :: eqwtd(:, :) + real, allocatable :: riverbed(:, :) + real, allocatable :: rivercond(:, :) + real, allocatable :: rechclim(:,:) + + !inout + real, allocatable :: smois(:, :, :) ! SMC + real, allocatable :: sh2oxy(:, :, :) ! SH2O + real, allocatable :: wtd(:, :) ! the depth to water table [m] + real, allocatable :: smcwtd(:, :) ! soil moisture between bottom of the soil and the water table [m3/m3] + real, allocatable :: deeprech(:, :) ! + real, allocatable :: qslat(:, :) ! lateral flow + real, allocatable :: qrfs(:, :) + real, allocatable :: qsprings(:, :) + real, allocatable :: rech(:, :) ! groundwater recharge (net vertical flux across the water table), positive up + real, allocatable :: xland(:, :) ! XLAND ! =2 ocean; =1 land/seaice + real, allocatable :: xice(:, :) ! sea ice fraction + real, allocatable :: pexp(:, :) + real, allocatable :: qrf(:, :) + real, allocatable :: qspring(:, :) + + ! local variables, need to be very careful + integer:: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + + integer :: tmp_row, tmp_col + ! local + real :: xice_threshold + integer :: isice + + + ! used to be parameters but set to constant number + + + allocate( tmp_sldpth( NOAHMP401_struc(n)%nsoil ) ) allocate( tmp_shdfac_monthly( 12 ) ) allocate( tmp_soilcomp( 8 ) ) @@ -283,6 +336,15 @@ subroutine NoahMP401_main(n) lat = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lat lon = LIS_domain(n)%grid(LIS_domain(n)%gindex(col, row))%lon + ! TML print statement for model debugging. + tmp_printdebug = 0 + if (row .eq. 12) then + if (col .eq. 25) then + tmp_printdebug = 1 + !print *, "PRINTED DEBUG VARIABLES ON" + endif + endif + ! retrieve forcing data from NOAHMP401_struc(n)%noahmp401(t) and assign to local variables ! tair: air temperature tmp_tair = NOAHMP401_struc(n)%noahmp401(t)%tair / NOAHMP401_struc(n)%forc_count @@ -308,9 +370,9 @@ subroutine NoahMP401_main(n) tmp_lwdown = NOAHMP401_struc(n)%noahmp401(t)%lwdown / NOAHMP401_struc(n)%forc_count ! prcp: total precipitation (rainfall+snowfall) - ! Both Noah-MP-3.6 and Noah-MP-4.0.1 require total precipitation as forcing input. - ! In Noah-MP-3.6, the forcing is required to be precipitation rate [kg m-2 sec-1]. - ! In Noah-MP-4.0.1, the forcing is required to be precipitation amount [kg m-2]. + ! Both NoahMP-3.6.1 and NoahMP-4.0.1 requires total precipitation as forcing input. + ! In LIS/NoahMP-3.6.1, the input forcing is total precipitation [mm], but in + ! LIS/NoahMP-4.0.1, the forcing data provides precipitation rate [mm/s] !!! ! T. Lahmers: Correct total precip for cases when model time step > forcing timestep. ! Edit suggested by D. Mocko and K. Arsenault @@ -755,7 +817,8 @@ subroutine NoahMP401_main(n) NOAHMP401_struc(n)%noahmp401(t)%param, & ! out - relative soil moisture [-] tmp_sfcheadrt , & tmp_infxs1rt , & - tmp_soldrain1rt ) ! out - extra output for WRF-HYDRO [m] + tmp_soldrain1rt , & !) ! out - extra output for WRF-HYDRO [m] + tmp_printdebug ) ! TML: Added debug print statement ! save state variables from local variables to global variables NOAHMP401_struc(n)%noahmp401(t)%sfcrunoff = tmp_sfcrunoff @@ -783,6 +846,7 @@ subroutine NoahMP401_main(n) NOAHMP401_struc(n)%noahmp401(t)%qsnow = tmp_qsnow NOAHMP401_struc(n)%noahmp401(t)%wslake = tmp_wslake NOAHMP401_struc(n)%noahmp401(t)%zwt = tmp_zwt + NOAHMP401_struc(n)%noahmp401(t)%wtd = tmp_zwt !!! wtd should be the same as zwt NOAHMP401_struc(n)%noahmp401(t)%wa = tmp_wa NOAHMP401_struc(n)%noahmp401(t)%wt = tmp_wt NOAHMP401_struc(n)%noahmp401(t)%tsno(:) = tmp_tsno(:) @@ -867,6 +931,11 @@ subroutine NoahMP401_main(n) NOAHMP401_struc(n)%noahmp401(t)%chuc = tmp_chuc NOAHMP401_struc(n)%noahmp401(t)%chv2 = tmp_chv2 NOAHMP401_struc(n)%noahmp401(t)%chb2 = tmp_chb2 + !SW + noahmp401_struc(n)%noahmp401(t)%subsnow = tmp_subsnow + noahmp401_struc(n)%noahmp401(t)%qsnbot = tmp_qsnbot + noahmp401_struc(n)%noahmp401(t)%pah = tmp_pah + noahmp401_struc(n)%noahmp401(t)%relsmc(:) = tmp_relsmc(:) NOAHMP401_struc(n)%noahmp401(t)%infxs1rt = tmp_infxs1rt(1,1) NOAHMP401_struc(n)%noahmp401(t)%soldrain1rt = tmp_soldrain1rt(1,1) @@ -880,6 +949,265 @@ subroutine NoahMP401_main(n) tmp_q2sat = 0.622*tmp_es/(tmp_psurf-(1.-0.622)*tmp_es) noahmp401_struc(n)%noahmp401(t)%rhmin = tmp_qair / tmp_q2sat endif + enddo ! end of tile (t) loop, for Noah-MP physics + + !!! MMF scheme, SW + if(NOAHMP401_struc(n)%run_opt .eq. 5) then + ! determine nrow, ncol + min_row = noahmp401_struc(n)%row_min + max_row = noahmp401_struc(n)%row_max + min_col = noahmp401_struc(n)%col_min + max_col = noahmp401_struc(n)%col_max + !do t = 1, LIS_rc%npatch(n, LIS_rc%lsm_index) + ! dt = LIS_rc%ts + ! row = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row + ! col = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col + ! min_row = min(min_row, row) + ! max_row = max(max_row, row) + ! min_col = min(min_col, col) + ! max_col = max(max_col, col) + !enddo + + !print *, 'max row' + !print *, max_row + !print *, 'max col' + !print *, max_col + ims = min_col + ime = max_col + jms = min_row + jme = max_row + ids = min_col + ide = max_col + jds = min_row + jde = max_row + its = min_col + ite = max_col + jts = min_row + jte = max_row + + kds = 1 + kde = 1 + kms = 1 + kme = 2 + kts = 1 + kte = 2 + nsoil = NOAHMP401_struc(n)%nsoil + + ! assume i corresponds to col and j corresponds to row + ! this assumption should not matter to results + allocate( fdepth(min_col:max_col, min_row:max_row)) ! P + allocate( area(min_col:max_col, min_row:max_row)) ! P + allocate( topo(min_col:max_col, min_row:max_row)) ! P + allocate( eqwtd(min_col:max_col, min_row:max_row)) ! P + allocate( pexp(min_col:max_col, min_row:max_row)) ! P + allocate( riverbed(min_col:max_col, min_row:max_row)) ! P + allocate( rivercond(min_col:max_col, min_row:max_row)) ! P + allocate( rechclim(min_col:max_col, min_row:max_row)) ! P + + allocate( isltyp(min_col:max_col, min_row:max_row)) ! P + allocate( ivgtyp(min_col:max_col, min_row:max_row)) ! P + + allocate( smoiseq(min_col:max_col, 1:NOAHMP401_struc(n)%nsoil, min_row:max_row)) ! S + allocate( smois(min_col:max_col, 1:NOAHMP401_struc(n)%nsoil, min_row:max_row)) ! S + allocate( sh2oxy(min_col:max_col, 1:NOAHMP401_struc(n)%nsoil, min_row:max_row)) ! S + allocate( dzs(NOAHMP401_struc(n)%nsoil)) ! P + allocate( wtd(min_col:max_col, min_row:max_row)) ! S + allocate( smcwtd(min_col:max_col, min_row:max_row)) ! S + allocate( deeprech(min_col:max_col, min_row:max_row)) ! S recharge to or from the water table when deep [m] + allocate( qslat(min_col:max_col, min_row:max_row)) ! S cumulative term of qlat + allocate( qrfs(min_col:max_col, min_row:max_row)) ! S cumulative term of groundwater-river water flux + allocate( qsprings(min_col:max_col, min_row:max_row)) ! S cumulative term of water springing at the surface from groundwater convergence in the column + allocate( rech(min_col:max_col, min_row:max_row)) ! S recharge to or from the water table when shallow [m] (diagnostic) + allocate( xland(min_col:max_col, min_row:max_row)) ! P + allocate( xice(min_col:max_col, min_row:max_row)) ! P + allocate( qrf(min_col:max_col, min_row:max_row)) ! O groundwater - river water flux + allocate( qspring(min_col:max_col, min_row:max_row)) ! O water springing at the surface from groundwater convergence in the column + + ! we may not run sub-grid tiling with MMF + + ! land, ice, sea ice setup + xland(:,:) = 1.0 ! the grid is fully filled with land + xice(:,:) = 0.0 ! ice fraction is 0 + xice_threshold = 0.5 ! be consistent with hrldas + isice = 100 ! this number should not be one of vegetation type id + isurban = LIS_rc%urbanclass ! be consistent with hrldas + ! soil layer thickness + dzs(:) = tmp_sldpth(:) + + do row = min_row, max_row + do col = min_col, max_col + t = NOAHMP401_struc(n)%rct_idx(col,row) ! rct_idx is col x row TML + if(t .eq. LIS_rc%udef) then ! undefined for land, but maybe still valid for MMF + xland(col,row) = 2 ! 2 is used for the undefined grid cells in the HRLDAS test case + isltyp(col,row) = NOAHMP401_struc(n)%soil2d(col,row) ! required by MMF, soil2d is col x row TML + ivgtyp(col,row) = NOAHMP401_struc(n)%vege2d(col,row) ! vege2d is col x row TML + else + ! soil type + isltyp(col,row) = NOAHMP401_struc(n)%noahmp401(t)%soiltype + ! vegetation type + ivgtyp(col,row) = NOAHMP401_struc(n)%noahmp401(t)%vegetype + ! soil moisture + smois(col, :, row) = NOAHMP401_struc(n)%noahmp401(t)%smc(:) + sh2oxy(col, :, row) = NOAHMP401_struc(n)%noahmp401(t)%sh2o(:) + smoiseq(col, :, row) = NOAHMP401_struc(n)%noahmp401(t)%smoiseq(:) + smcwtd(col, row) = NOAHMP401_struc(n)%noahmp401(t)%smcwtd + deeprech(col, row) = NOAHMP401_struc(n)%noahmp401(t)%deeprech + rechclim(col, row) = NOAHMP401_struc(n)%noahmp401(t)%rechclim + wtd(col,row) = noahmp401_struc(n)%noahmp401(t)%wtd + + rech(col,row) = noahmp401_struc(n)%noahmp401(t)%rech + qslat(col,row) = noahmp401_struc(n)%noahmp401(t)%qslat + qrfs(col,row) = noahmp401_struc(n)%noahmp401(t)%qrfs + qsprings(col,row) = noahmp401_struc(n)%noahmp401(t)%qsprings + qrf(col,row) = noahmp401_struc(n)%noahmp401(t)%qrf + qspring(col,row) = noahmp401_struc(n)%noahmp401(t)%qspring + + endif + + if(isltyp(col,row) .eq. 14) then + smcwtd(col, row) = 1.0 + wtd(col,row) = 0.0 + endif + + ! MMF parameters may be needed for undefined grid cells + ridx = row - NOAHMP401_struc(n)%row_min + 1 + cidx = col - NOAHMP401_struc(n)%col_min + 1 + fdepth(col,row) = NOAHMP401_struc(n)%fdepth(cidx, ridx) !! all variables col x row TML + topo(col,row) = NOAHMP401_struc(n)%topo(cidx, ridx) + area(col,row) = NOAHMP401_struc(n)%area(cidx, ridx) + riverbed(col,row) = NOAHMP401_struc(n)%riverbed(cidx, ridx) + eqwtd(col,row) = NOAHMP401_struc(n)%eqwtd(cidx, ridx) + rivercond(col,row) = NOAHMP401_struc(n)%rivercond(cidx, ridx) + rechclim(col,row) = NOAHMP401_struc(n)%rechclim(cidx, ridx) + pexp(col,row) = 1.0 + enddo ! col loop + enddo ! row loop + + wtddt = LIS_rc%ts/60.0 ! time step in minutes? + + !print*, 'WTABLE_mmf_noahmp Initial Values' + !print*, 'FDEPTH = ',fdepth(25,12) + !print*, 'AREA = ',area(25,12) + !print*, 'TOPO = ',topo(25,12) + !print*, 'ISURBAN = ',isurban(25,12) + !print*, 'IVGTYP = ',ivgtyp(25,12) + !print*, 'RIVERCOND = ',rivercond(25,12) + !print*, 'RIVERBED = ',riverbed(25,12) + !print*, 'EQWTD = ',eqwtd(25,12) + !print*, 'PEXP = ',pexp(25,12) + !print*, 'SMC-1 = ',smois(25,1,12) + !print*, 'SMC-2 = ',smois(25,2,12) + !print*, 'SMC-3 = ',smois(25,3,12) + !print*, 'SMC-4 = ',smois(25,4,12) + !print*, 'SH2O-1 = ',sh2oxy(25,1,12) + !print*, 'SH2O-2 = ',sh2oxy(25,2,12) + !print*, 'SH2O-3 = ',sh2oxy(25,3,12) + !print*, 'SH2O-4 = ',sh2oxy(25,4,12) + !print*, 'SMCWTD = ',smcwtd(25,12) + !print*, 'WTD = ',wtd(25,12) + !print*, 'QRF = ',qrf(25,12) + !print*, 'DEEPRECH = ',deeprech(25,12) + !print*, 'QSPRING = ',qspring(25,12) + !print*, 'QSLAT = ',qslat(25,12) + !print*, 'QRFS = ',qrfs(25,12) + !print*, 'QSPRINGS = ',qsprings(25,12) + !print*, 'RECH = ',rech(25,12) + + !!! call MMF physics + call WTABLE_mmf_noahmp (nsoil ,xland ,xice ,xice_threshold ,isice ,& !in + isltyp ,smoiseq ,dzs ,wtddt ,& !in + fdepth ,area ,topo ,isurban ,ivgtyp ,& !in + rivercond ,riverbed ,eqwtd ,pexp ,& !in + smois ,sh2oxy ,smcwtd ,wtd ,qrf ,& !inout + deeprech ,qspring ,qslat ,qrfs ,qsprings ,rech ,& !inout + ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte ) + + !print*, 'WTABLE_mmf_noahmp Final Values' + !print*, 'FDEPTH = ',fdepth(25,12) + !print*, 'AREA = ',area(25,12) + !print*, 'TOPO = ',topo(25,12) + !print*, 'ISURBAN = ',isurban(25,12) + !print*, 'IVGTYP = ',ivgtyp(25,12) + !print*, 'RIVERCOND = ',rivercond(25,12) + !print*, 'RIVERBED = ',riverbed(25,12) + !print*, 'EQWTD = ',eqwtd(25,12) + !print*, 'PEXP = ',pexp(25,12) + !print*, 'SMC-1 = ',smois(25,1,12) + !print*, 'SMC-2 = ',smois(25,2,12) + !print*, 'SMC-3 = ',smois(25,3,12) + !print*, 'SMC-4 = ',smois(25,4,12) + !print*, 'SH2O-1 = ',sh2oxy(25,1,12) + !print*, 'SH2O-2 = ',sh2oxy(25,2,12) + !print*, 'SH2O-3 = ',sh2oxy(25,3,12) + !print*, 'SH2O-4 = ',sh2oxy(25,4,12) + !print*, 'SMCWTD = ',smcwtd(25,12) + !print*, 'WTD = ',wtd(25,12) + !print*, 'QRF = ',qrf(25,12) + !print*, 'DEEPRECH = ',deeprech(25,12) + !print*, 'QSPRING = ',qspring(25,12) + !print*, 'QSLAT = ',qslat(25,12) + !print*, 'QRFS = ',qrfs(25,12) + !print*, 'QSPRINGS = ',qsprings(25,12) + !print*, 'RECH = ',rech(25,12) + + ! get state and output + do row = min_row, max_row + do col = min_col, max_col + NOAHMP401_struc(n)%eqwtd(col, row) = eqwtd(col,row) + + t = NOAHMP401_struc(n)%rct_idx(col,row) ! rct_idx is row major, Shugong Wang + if (t .eq. LIS_rc%udef) cycle + noahmp401_struc(n)%noahmp401(t)%smc(:) = smois(col, :, row) + noahmp401_struc(n)%noahmp401(t)%sh2o(:) = sh2oxy(col, :, row) + noahmp401_struc(n)%noahmp401(t)%smoiseq(:) = smoiseq(col, :, row) + noahmp401_struc(n)%noahmp401(t)%smcwtd = smcwtd(col,row) + noahmp401_struc(n)%noahmp401(t)%wtd = wtd(col,row) !TML added wtd and zwt. + noahmp401_struc(n)%noahmp401(t)%zwt = wtd(col,row) !Same variable... + noahmp401_struc(n)%noahmp401(t)%deeprech = deeprech(col,row) + noahmp401_struc(n)%noahmp401(t)%rech = rech(col,row) + noahmp401_struc(n)%noahmp401(t)%qslat = qslat(col,row) + noahmp401_struc(n)%noahmp401(t)%qrfs = qrfs(col,row) + noahmp401_struc(n)%noahmp401(t)%qsprings = qsprings(col,row) + noahmp401_struc(n)%noahmp401(t)%qrf = qrf(col,row) + noahmp401_struc(n)%noahmp401(t)%qspring = qspring(col,row) + enddo ! col loop + enddo ! row loop + + ! free memory + deallocate( fdepth) + deallocate( area) + deallocate( topo) + deallocate( eqwtd) + deallocate( pexp) + deallocate( riverbed) + deallocate(rivercond) + deallocate( rechclim) + + deallocate( isltyp) + deallocate( ivgtyp) + + deallocate( smoiseq) + deallocate( smois) + deallocate( sh2oxy) + deallocate( dzs) + deallocate( wtd) + deallocate( smcwtd) + deallocate( deeprech) + deallocate( qslat) + deallocate( qrfs) + deallocate( qsprings) + deallocate( rech) + deallocate( xland) + deallocate( xice) + deallocate( qrf) + deallocate( qspring) + endif ! if run MMF + + + ! write output + do t = 1, LIS_rc%npatch(n, LIS_rc%lsm_index) call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_RHMIN, & value=noahmp401_struc(n)%noahmp401(t)%rhmin, & @@ -914,7 +1242,8 @@ subroutine NoahMP401_main(n) call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_ALBEDO, value = NOAHMP401_struc(n)%noahmp401(t)%albedo, & vlevel=1, unit="-", direction="-", surface_type = LIS_rc%lsm_index) - if (tmp_albedo.ne.LIS_rc%udef) tmp_albedo = tmp_albedo * 100.0 + !if (tmp_albedo.ne.LIS_rc%udef) tmp_albedo = tmp_albedo * 100.0 + if (NOAHMP401_struc(n)%noahmp401(t)%albedo .ne. LIS_rc%udef) tmp_albedo = NOAHMP401_struc(n)%noahmp401(t)%albedo * 100.0 call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_ALBEDO, value = tmp_albedo, & vlevel=1, unit="%", direction="-", surface_type = LIS_rc%lsm_index) @@ -1310,11 +1639,13 @@ subroutine NoahMP401_main(n) ! Code added by Zhuo Wang on 02/28/2019 ![ 92] output variable: qsnbot (unit=kg m-2 s-1). *** melting water out of snow bottom - call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QSM, value = tmp_qsnbot, & + !call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QSM, value = tmp_qsnbot, & + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QSM, value = NOAHMP401_struc(n)%noahmp401(t)%qsnbot, & vlevel=1, unit="kg m-2 s-1", direction="S2L", surface_type = LIS_rc%lsm_index) ![ 93] output variable: subsnow (unit=kg m-2 s-1). *** snow sublimation - call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_SUBSNOW, value = tmp_subsnow, & + !call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_SUBSNOW, value = tmp_subsnow, & + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_SUBSNOW, value = NOAHMP401_struc(n)%noahmp401(t)%subsnow, & vlevel=1, unit="kg m-2 s-1", direction="-", surface_type = LIS_rc%lsm_index) ![ 94] output variable: AvgSurfT (unit=K). *** average surface temperature @@ -1331,10 +1662,12 @@ subroutine NoahMP401_main(n) (NOAHMP401_struc(n)%noahmp401(t)%canliq + & NOAHMP401_struc(n)%noahmp401(t)%canice) endif + do i = 1,NOAHMP401_struc(n)%nsoil TWS_out = TWS_out + & (NOAHMP401_struc(n)%noahmp401(t)%smc(i) * & - tmp_sldpth(i)*LIS_CONST_RHOFW) + !tmp_sldpth(i)*lis_const_rhofw) + NOAHMP401_struc(n)%sldpth(i)*lis_const_rhofw) enddo if (NOAHMP401_struc(n)%noahmp401(t)%wa.ge.0.0) then TWS_out = TWS_out + NOAHMP401_struc(n)%noahmp401(t)%wa @@ -1344,9 +1677,30 @@ subroutine NoahMP401_main(n) ![ 96] Qa - Advective energy - Heat transferred to a snow cover by rain ! - (unit=W m-2) - added by David Mocko - call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QA, value = tmp_pah, & + !call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QA, value = tmp_pah, & + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QA, value = NOAHMP401_struc(n)%noahmp401(t)%pah, & vlevel=1, unit="W m-2",direction="DN",surface_type=LIS_rc%lsm_index) + ![ 97] Qslat - Accumulated Lateral Flow (for MMF Groundwater) - (unit=mm) + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QSLAT, value = NOAHMP401_struc(n)%noahmp401(t)%qslat, & + vlevel=1, unit="mm", direction="-", surface_type = LIS_rc%lsm_index) + + ![ 98] Qsprings - Accumulated Seeping Water (for MMF Groundwater) - (unit=mm) + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QSPRINGS, value = NOAHMP401_struc(n)%noahmp401(t)%qsprings, & + vlevel=1, unit="mm", direction="-", surface_type = LIS_rc%lsm_index) + + ![ 99] Qrfs - Accumulated Flux From Groundwater to Rivers (for MMF Groundwater) - (unit=mm) + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_QRFS, value = NOAHMP401_struc(n)%noahmp401(t)%qrfs, & + vlevel=1, unit="mm", direction="-", surface_type = LIS_rc%lsm_index) + + ![ 100] Wtd - Water Table Depth (for MMF Groundwater) - (unit=m) + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_WTD, value = NOAHMP401_struc(n)%noahmp401(t)%wtd, & + vlevel=1, unit="m", direction="-", surface_type = LIS_rc%lsm_index) + + ![ 101] Rech - Accumulated Recharge (for MMF Groundwater) - (unit=mm) + call LIS_diagnoseSurfaceOutputVar(n, t, LIS_MOC_RECH, value = NOAHMP401_struc(n)%noahmp401(t)%rech, & + vlevel=1, unit="mm", direction="-", surface_type = LIS_rc%lsm_index) + ! Added water balance change terms - David Mocko endsm = 0.0 do i = 1,tmp_nsoil @@ -1379,25 +1733,25 @@ subroutine NoahMP401_main(n) surface_type=LIS_rc%lsm_index) ! David Mocko (10/29/2019) - Copy RELSMC calculation from Noah-3.X - do i = 1,tmp_nsoil - if (tmp_relsmc(i).gt.1.0) then - tmp_relsmc(i) = 1.0 + do i = 1,NOAHMP401_struc(n)%nsoil + if (NOAHMP401_struc(n)%noahmp401(t)%relsmc(i).gt.1.0) then + NOAHMP401_struc(n)%noahmp401(t)%relsmc(i) = 1.0 endif - if (tmp_relsmc(i).lt.0.01) then - tmp_relsmc(i) = 0.01 + if (NOAHMP401_struc(n)%noahmp401(t)%relsmc(i).lt.0.01) then + NOAHMP401_struc(n)%noahmp401(t)%relsmc(i) = 0.01 endif ! J.Case (9/11/2014) -- Set relative soil moisture to missing (LIS_rc%udef) ! if the vegetation type is urban class. - if (tmp_vegetype.eq.tmp_urban_vegetype) then - tmp_relsmc(i) = LIS_rc%udef + if (NOAHMP401_struc(n)%noahmp401(t)%vegetype.eq.tmp_urban_vegetype) then + NOAHMP401_struc(n)%noahmp401(t)%relsmc(i) = LIS_rc%udef endif call LIS_diagnoseSurfaceOutputVar(n,t,LIS_MOC_RELSMC,vlevel=i, & - value=tmp_relsmc(i),unit='-',direction="-",surface_type=LIS_rc%lsm_index) - if (tmp_relsmc(i).eq.LIS_rc%udef) then - tempval = tmp_relsmc(i) + value=NOAHMP401_struc(n)%noahmp401(t)%relsmc(i),unit='-',direction="-",surface_type=LIS_rc%lsm_index) + if (NOAHMP401_struc(n)%noahmp401(t)%relsmc(i).eq.LIS_rc%udef) then + tempval = NOAHMP401_struc(n)%noahmp401(t)%relsmc(i) else - tempval = tmp_relsmc(i)*100.0 + tempval = NOAHMP401_struc(n)%noahmp401(t)%relsmc(i)*100.0 endif call LIS_diagnoseSurfaceOutputVar(n,t,LIS_MOC_RELSMC,vlevel=i, & value=tempval,unit='%',direction="-",surface_type=LIS_rc%lsm_index) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_module.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_module.F90 index 1b47bd72d..be4703c47 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_module.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_module.F90 @@ -434,9 +434,22 @@ module NoahMP401_module real :: sai real :: tauss real, pointer :: smoiseq(:) - real :: smcwtd - real :: deeprech - real :: rech + real :: smcwtd !MMF state + real :: deeprech !MMF state + real :: rech !MMF state + real :: wtd !MMF state, SW, 09142021 + real :: qslat !MMF state, SW + real :: qrfs !MMF state, SW + real :: qsprings !MMF state, SW + real :: fdepth !SW for MMF 09142021 + real :: area !SW for MMF + real :: topo !SW for MMF + real :: eqwtd !SW for MMF + real :: riverbed !SW for MMF + real :: rivercond !SW for MMF + real :: qrf !MMF output SW for MMF + real :: qspring !MMF output, SW + real :: rechclim !SW real :: grain real :: gdd integer :: pgs @@ -506,7 +519,12 @@ module NoahMP401_module real :: chuc real :: chv2 real :: chb2 - + + !SW + real :: qsnbot + real :: subsnow + real :: pah + real, pointer :: relsmc(:) !EMK for 557WW real :: tair_agl_min real :: rhmin diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readrst.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readrst.F90 index ca55e253e..021e3dbe5 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readrst.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readrst.F90 @@ -24,7 +24,7 @@ ! !INTERFACE: subroutine NoahMP401_readrst() ! !USES: - use LIS_coreMod, only : LIS_rc, LIS_masterproc + use LIS_coreMod, only : LIS_rc, LIS_masterproc, LIS_surface use LIS_historyMod, only : LIS_readvar_restart use LIS_logMod, only : LIS_logunit, LIS_endrun, & LIS_getNextUnitNumber, & @@ -110,6 +110,7 @@ subroutine NoahMP401_readrst() implicit none integer :: t, l + integer :: row, col, ridx, cidx integer :: nc, nr, npatch integer :: n integer :: ftn @@ -428,13 +429,42 @@ subroutine NoahMP401_readrst() varname="PGS", wformat=wformat) ! read: optional gecros crop -! do l=1, 60 ! TODO: check loop -! call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, tmptilen, varname="GECROS_STATE", & -! dim=l, vlevels = 60, wformat=wformat) -! do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) -! NOAHMP401_struc(n)%noahmp401(t)%gecros_state(l) = tmptilen(t) -! enddo -! enddo + !do l=1, 60 ! TODO: check loop + ! call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, tmptilen, varname="GECROS_STATE", & + ! dim=l, vlevels = 60, wformat=wformat) + ! do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) + ! NOAHMP401_struc(n)%noahmp401(t)%gecros_state(l) = tmptilen(t) + ! enddo + !enddo + + !SW for MMF + if(NOAHMP401_struc(n)%run_opt==5) then + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qslat, & + varname="QSLAT", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qrfs, & + varname="QRFS", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qsprings, & + varname="QSPRINGS", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%wtd, & + varname="WTD", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%eqwtd, & + varname="EQWTD", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%rivercond, & + varname="RIVERCOND", wformat=wformat) + call LIS_readvar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%riverbed, & + varname="RIVERBED", wformat=wformat) + do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) + col = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col + row = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row + cidx = col - NOAHMP401_struc(n)%col_min + 1 + ridx = row - NOAHMP401_struc(n)%row_min + 1 + + NOAHMP401_struc(n)%eqwtd(cidx, ridx) = NOAHMP401_struc(n)%noahmp401(t)%eqwtd + NOAHMP401_struc(n)%riverbed(cidx, ridx) = NOAHMP401_struc(n)%noahmp401(t)%riverbed + NOAHMP401_struc(n)%rivercond(cidx, ridx) = NOAHMP401_struc(n)%noahmp401(t)%rivercond + enddo + + endif ! close restart file if(wformat .eq. "binary") then diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_setup.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_setup.F90 index 7799d9952..5d7a31ead 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_setup.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_setup.F90 @@ -61,8 +61,8 @@ subroutine NoahMP401_setup() implicit none integer :: mtype integer :: t, k, n - integer :: col, row - real, allocatable :: placeholder(:,:) + integer :: col, row, ridx, cidx + real, allocatable :: placeholder(:,:), tmp integer :: soilcolor, vegtyp, soiltyp(4), slopetyp, croptype mtype = LIS_rc%lsm_index @@ -228,6 +228,174 @@ subroutine NoahMP401_setup() enddo enddo endif + + ! MMF scheme support, SW, 09142021 + if(NOAHMP401_struc(n)%run_opt .eq. 5) then + ! read: fdepth + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter FDEPTH from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "MMF_FDEPTH", placeholder) + !! Note: the tile based data structure introduce missing values to MMF parameters where landmask=0 + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%fdepth = placeholder(col, row) + enddo + !!! the 2-D array should solve the missing values of MMF parameters where landmask=0 + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%fdepth(cidx, ridx) = placeholder(col, row) + enddo + enddo + + ! read: area + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter AREA from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "AREAXY", placeholder) + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%area = placeholder(col, row)*1E6 ! m2 + enddo + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%area(cidx, ridx) = placeholder(col, row)*1E6 ! m2 + enddo + enddo + ! read: topo + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter TOPO/ELEVATION from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "MMF_HGTM", placeholder) + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%topo = placeholder(col, row) + enddo + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%topo(cidx, ridx) = placeholder(col, row) + enddo + enddo + ! read: eqwtd + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter EQWTD from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "MMF_EQWTD", placeholder) + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%eqwtd = placeholder(col, row) + enddo + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%eqwtd(cidx, ridx) = placeholder(col, row) + enddo + enddo + ! read: RECHCLIM + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter RECHCLIM from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "MMF_RECHCLIM", placeholder) + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%rechclim = placeholder(col, row) + enddo + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%rechclim(cidx, ridx) = placeholder(col, row) + enddo + enddo + ! read: riverbed + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter RIVERBED from ", & + trim(LIS_rc%paramfile(n)) + call LIS_read_param(n, "MMF_RIVERBED", placeholder) + do t = 1, LIS_rc%npatch(n, mtype) + col = LIS_surface(n, mtype)%tile(t)%col + row = LIS_surface(n, mtype)%tile(t)%row + NOAHMP401_struc(n)%noahmp401(t)%riverbed = placeholder(col, row) + enddo + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%riverbed(cidx, ridx) = placeholder(col, row) + enddo + enddo + !print*, 'LOADED RIVERBED = ',NOAHMP401_struc(n)%riverbed(25,12) + + + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter TEXTURE from ", & + trim(LIS_rc%paramfile(n)) + do k = 1, LIS_rc%nsoiltypes ! 16 is hardcoded temporarily + call NOAHMP401_read_MULTILEVEL_param(n, "TEXTURE", k, placeholder) + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%soil3d(cidx, ridx, k) = placeholder(col, row) + enddo + enddo + enddo + + ! determine soil type + NOAHMP401_struc(n)%soil2d(:,:) = 1 + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + + do k=2, LIS_rc%nsoiltypes + if(NOAHMP401_struc(n)%soil3d(cidx, ridx, k) > NOAHMP401_struc(n)%soil3d(cidx, ridx, NOAHMP401_struc(n)%soil2d(cidx,ridx))) then + NOAHMP401_struc(n)%soil2d(cidx,ridx) = k + endif + enddo + enddo + enddo + + write(LIS_logunit,*) "[INFO] Noah-MP.4.0.1 reading parameter LANDCOVER from ", & + trim(LIS_rc%paramfile(n)) + do k = 1, LIS_rc%nsurfacetypes ! 20 is hardcoded temporarily + call NOAHMP401_read_MULTILEVEL_param(n, "LANDCOVER", k, placeholder) + ! 2-D array + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + row = ridx - NOAHMP401_struc(n)%row_min + 1 + col = cidx - NOAHMP401_struc(n)%col_min + 1 + NOAHMP401_struc(n)%vege3d(cidx, ridx, k) = placeholder(col, row) + enddo + enddo + enddo + ! determine vegetype + NOAHMP401_struc(n)%vege2d(:,:) = 1 ! Crrently, we have 20 layer of land cover for LIS. However, there could be a 21 layer in HRLDAS. vegetation is set to 17 for such case + do ridx = NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do cidx = NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + + do k=2, LIS_rc%nsurfacetypes + if(NOAHMP401_struc(n)%vege3d(cidx, ridx, k) > NOAHMP401_struc(n)%vege3d(cidx, ridx, NOAHMP401_struc(n)%vege2d(cidx,ridx))) then + NOAHMP401_struc(n)%vege2d(cidx,ridx) = k + endif + enddo + + ! if the total fraction from layer 1 to layer 20 is no more than 0.5, this means there is the fraction of the 21st layer is >= 0.5, set the landcover type to be 17 as HRLDAS + if(.not. sum(NOAHMP401_struc(n)%vege3d(cidx, ridx, :))>0.5) then + NOAHMP401_struc(n)%vege2d(cidx,ridx) = LIS_rc%waterclass + endif + enddo + enddo + endif deallocate(placeholder) !!!! read Noah-MP parameter tables Shugong Wang 11/06/2018 diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_writerst.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_writerst.F90 index eb44a25fc..80c27bcc2 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_writerst.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_writerst.F90 @@ -23,7 +23,7 @@ ! !INTERFACE: subroutine NoahMP401_writerst(n) ! !USES: - use LIS_coreMod, only : LIS_rc, LIS_masterproc + use LIS_coreMod, only : LIS_rc, LIS_masterproc, LIS_surface use LIS_timeMgrMod, only : LIS_isAlarmRinging use LIS_logMod, only : LIS_logunit, LIS_getNextUnitNumber, & LIS_releaseUnitNumber , LIS_verify @@ -123,7 +123,7 @@ end subroutine NoahMP401_writerst subroutine NoahMP401_dump_restart(n, ftn, wformat) ! !USES: - use LIS_coreMod, only : LIS_rc, LIS_masterproc + use LIS_coreMod, only : LIS_rc, LIS_masterproc, LIS_surface use LIS_logMod, only : LIS_logunit use LIS_historyMod use NoahMP401_lsmMod @@ -217,7 +217,8 @@ subroutine NoahMP401_dump_restart(n, ftn, wformat) ! !EOP - integer :: l, t + integer :: l, t + integer :: row, col, ridx, cidx real :: tmptilen(LIS_rc%npatch(n, LIS_rc%lsm_index)) integer :: dimID(11) integer :: sfcrunoff_ID @@ -262,6 +263,15 @@ subroutine NoahMP401_dump_restart(n, ftn, wformat) integer :: tauss_ID integer :: smoiseq_ID integer :: smcwtd_ID + integer :: wtd_ID + integer :: eqwtd_ID + integer :: rivercond_ID + integer :: riverbed_ID + integer :: qslat_ID + integer :: qrfs_ID + integer :: qsprings_ID + integer :: qrf_ID + integer :: qspring_ID integer :: deeprech_ID integer :: rech_ID integer :: grain_ID @@ -779,12 +789,70 @@ subroutine NoahMP401_dump_restart(n, ftn, wformat) varid=pgs_ID, dim=1, wformat=wformat) ! optional gecros crop -! do l=1, 60 ! TODO: check loop -! tmptilen = 0 -! do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) -! tmptilen(t) = NOAHMP401_struc(n)%noahmp401(t)%gecros_state(l) -! enddo -! call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, tmptilen, & -! varid=gecros_state_ID, dim=l, wformat=wformat) -! enddo + !do l=1, 60 ! TODO: check loop + ! tmptilen = 0 + ! do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) + ! tmptilen(t) = NOAHMP401_struc(n)%noahmp401(t)%gecros_state(l) + ! enddo + ! call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, tmptilen, & + ! varid=gecros_state_ID, dim=l, wformat=wformat) + !enddo + + !SW for MMF + if(NOAHMP401_struc(n)%run_opt==5) then + call LIS_writeHeader_restart(ftn, n, dimID, qslat_ID, "QSLAT", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qslat, & + varid=qslat_ID, dim=1, wformat=wformat) + + call LIS_writeHeader_restart(ftn, n, dimID, qrfs_ID, "QRFS", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qrfs, & + varid=qrfs_ID, dim=1, wformat=wformat) + + call LIS_writeHeader_restart(ftn, n, dimID, qsprings_ID, "QSPRINGS", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%qsprings, & + varid=qsprings_ID, dim=1, wformat=wformat) + + call LIS_writeHeader_restart(ftn, n, dimID, wtd_ID, "WTD", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%wtd, & + varid=wtd_ID, dim=1, wformat=wformat) + + !UPDATE 1D VARIABLES TO BE CONSITENT WITH 2D VARIABLES + do t=1, LIS_rc%npatch(n, LIS_rc%lsm_index) + col = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%col + row = LIS_surface(n, LIS_rc%lsm_index)%tile(t)%row + cidx = col - NOAHMP401_struc(n)%col_min + 1 + ridx = row - NOAHMP401_struc(n)%row_min + 1 + + NOAHMP401_struc(n)%noahmp401(t)%eqwtd = NOAHMP401_struc(n)%eqwtd(cidx, ridx) + NOAHMP401_struc(n)%noahmp401(t)%riverbed = NOAHMP401_struc(n)%riverbed(cidx, ridx) + NOAHMP401_struc(n)%noahmp401(t)%rivercond = NOAHMP401_struc(n)%rivercond(cidx, ridx) + enddo + + + call LIS_writeHeader_restart(ftn, n, dimID, eqwtd_ID, "EQWTD", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%eqwtd, & + varid=eqwtd_ID, dim=1, wformat=wformat) + + call LIS_writeHeader_restart(ftn, n, dimID, rivercond_ID, "RIVERCOND", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%rivercond, & + varid=rivercond_ID, dim=1, wformat=wformat) + + call LIS_writeHeader_restart(ftn, n, dimID, riverbed_ID, "RIVERBED", & + "soil moisture content in the layer to the water table when deep", & + "-", vlevels=1, valid_min=-99999.0, valid_max=99999.0) + call LIS_writevar_restart(ftn, n, LIS_rc%lsm_index, NOAHMP401_struc(n)%noahmp401%riverbed, & + varid=riverbed_ID, dim=1, wformat=wformat) + endif end subroutine NoahMP401_dump_restart diff --git a/lis/surfacemodels/land/noahmp.4.0.1/mmf_start.F90 b/lis/surfacemodels/land/noahmp.4.0.1/mmf_start.F90 new file mode 100644 index 000000000..f4f09f60f --- /dev/null +++ b/lis/surfacemodels/land/noahmp.4.0.1/mmf_start.F90 @@ -0,0 +1,383 @@ +subroutine mmf_start(n) + use LIS_coreMod + use NoahMP401_lsmMod + use module_sf_noahmpdrv_401 + use LIS_historyMod, only: LIS_gather_masterproc_2d_local_to_global, & + LIS_scatter_global_to_local_grid + use LIS_mpiMod + use LIS_logMod, only : LIS_logunit + implicit none + integer :: i,j + integer :: n, row, col, t, ridx, cidx, ierr + real :: wtddt + ! SW, MMF + integer, allocatable,dimension(:,:) :: isltyp, ivgtyp + !real, allocatable :: fdepth(:,:) + real, allocatable,dimension(:,:) :: fdepth, topo , area, rechclim, rivercond, & + wtd, riverbed, eqwtd, pexp, smcwtdxy, & + deeprechxy, rechxy, qslatxy, qrfsxy, qspringsxy + real, allocatable,dimension(:,:,:) :: smois, sh2o, smoiseq +#if (defined SPMD) + integer, allocatable, dimension(:,:) :: gisltyp, givgtyp + real, allocatable, dimension(:,:) :: gfdepth, gtopo , garea, grechclim, grivercond, & + gwtd, griverbed, geqwtd, gpexp, gsmcwtdxy, & + gdeeprechxy, grechxy, gqslatxy, gqrfsxy, gqspringsxy + real, allocatable,dimension(:,:,:) :: gsmois, gsh2o, gsmoiseq +#endif + wtddt = int(LIS_rc%ts/60) ! in minutes? + + allocate(isltyp(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(ivgtyp(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(fdepth(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(topo(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(area(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(rechclim(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(rivercond(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(wtd(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(riverbed(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(eqwtd(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(pexp(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(smcwtdxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(deeprechxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(rechxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(qslatxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(qrfsxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(qspringsxy(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(smois(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, 1:NOAHMP401_struc(n)%nsoil, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(sh2o(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, 1:NOAHMP401_struc(n)%nsoil, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + allocate(smoiseq(NOAHMP401_struc(n)%col_min:NOAHMP401_struc(n)%col_max, 1:NOAHMP401_struc(n)%nsoil, NOAHMP401_struc(n)%row_min:NOAHMP401_struc(n)%row_max)) + + do row=NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do col=NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + t = NOAHMP401_struc(n)%rct_idx(col,row) ! rct_idx is col x row TML + if(t .ne. LIS_rc%udef) then + isltyp(col,row) = NOAHMP401_struc(n)%noahmp401(t)%soiltype + ivgtyp(col,row) = NOAHMP401_struc(n)%noahmp401(t)%vegetype + else + isltyp(col,row) = NOAHMP401_struc(n)%soil2d(col,row) ! soil2d is col x row TML + ivgtyp(col,row) = NOAHMP401_struc(n)%vege2d(col,row) ! vege2d is col x row TML + endif + smois(col,:,row) = NOAHMP401_struc(n)%init_smc(:) + sh2o(col,:,row) = NOAHMP401_struc(n)%init_smc(:) + smoiseq(col,:,row) = 0.0 + enddo + enddo + !!! 2-D, MMF, SW + do row=NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do col=NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + ridx = row - NOAHMP401_struc(n)%row_min + 1 + cidx = col - NOAHMP401_struc(n)%col_min + 1 + fdepth(col,row) = NOAHMP401_struc(n)%fdepth(cidx, ridx) + topo(col,row) = NOAHMP401_struc(n)%topo(cidx, ridx) + area(col,row) = NOAHMP401_struc(n)%area(cidx, ridx) + riverbed(col,row) = NOAHMP401_struc(n)%riverbed(cidx, ridx) + eqwtd(col,row) = NOAHMP401_struc(n)%eqwtd(cidx, ridx) + rivercond(col,row) = NOAHMP401_struc(n)%rivercond(cidx, ridx) + rechclim(col,row) = NOAHMP401_struc(n)%rechclim(cidx, ridx) + pexp(col,row) = 1.0 + enddo + enddo + + !print*, 'mmf_start INIT VARIABLES: ' + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) + !print*, 'SMCWTD = ',smcwtdxy(18,12) + !print*, 'EQWTD = ',eqwtd(18,12) + !print*, 'WTD = ',wtd(18,12) + +#if (defined SPMD) + if ( LIS_masterproc ) then + ! Allocate global in and inout variables + write(LIS_logunit,*) "[INFO] Allocating Arrays for MMF Init." + allocate(gisltyp(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(givgtyp(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gfdepth(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gtopo(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(griverbed(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(geqwtd(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gpexp(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(garea(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gwtd(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gsmcwtdxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gdeeprechxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(grechxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gqslatxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gqrfsxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gqspringsxy(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(grechclim(LIS_rc%gnc(n), LIS_rc%gnr(n))) + allocate(gsmois(LIS_rc%gnc(n), NOAHMP401_struc(n)%nsoil, LIS_rc%gnr(n))) + allocate(gsh2o(LIS_rc%gnc(n), NOAHMP401_struc(n)%nsoil, LIS_rc%gnr(n))) + allocate(gsmoiseq(LIS_rc%gnc(n), NOAHMP401_struc(n)%nsoil, LIS_rc%gnr(n))) + else + ! Allocate dummy "global" in and inout variables + write(LIS_logunit,*) "[INFO] Allocating Dummy Arrays for MMF Init." + allocate(gisltyp(1,1)) + allocate(givgtyp(1,1)) + allocate(gfdepth(1,1)) + allocate(gtopo(1,1)) + allocate(griverbed(1,1)) + allocate(geqwtd(1,1)) + allocate(gpexp(1,1)) + allocate(garea(1,1)) + allocate(gwtd(1,1)) + allocate(gsmcwtdxy(1,1)) + allocate(gdeeprechxy(1,1)) + allocate(grechxy(1,1)) + allocate(gqslatxy(1,1)) + allocate(gqrfsxy(1,1)) + allocate(gqspringsxy(1,1)) + allocate(grechclim(1,1)) + allocate(gsmois(1,1,1)) + allocate(gsh2o(1,1,1)) + allocate(gsmoiseq(1,1,1)) + endif + + !TML ADD GATHER CALL HERE... + ! Gather in and inout variables + write(LIS_logunit,*) "[INFO] Gathering Arrays for MMF Init." + write(LIS_logunit,*) "[INFO] Gathering Array: isltyp" + call LIS_gather_masterproc_2d_local_to_global(n, isltyp, gisltyp) + write(LIS_logunit,*) "[INFO] Gathering Array: ivgtyp" + call LIS_gather_masterproc_2d_local_to_global(n, ivgtyp, givgtyp) + write(LIS_logunit,*) "[INFO] Gathering Array: fdepth" + call LIS_gather_masterproc_2d_local_to_global(n, fdepth, gfdepth) + write(LIS_logunit,*) "[INFO] Gathering Array: topo" + call LIS_gather_masterproc_2d_local_to_global(n, topo, gtopo) + call LIS_gather_masterproc_2d_local_to_global(n, riverbed, griverbed) + call LIS_gather_masterproc_2d_local_to_global(n, eqwtd, geqwtd) + call LIS_gather_masterproc_2d_local_to_global(n, pexp, gpexp) + call LIS_gather_masterproc_2d_local_to_global(n, area, garea) + call LIS_gather_masterproc_2d_local_to_global(n, wtd, gwtd) + call LIS_gather_masterproc_2d_local_to_global(n, smcwtdxy, gsmcwtdxy) + call LIS_gather_masterproc_2d_local_to_global(n, deeprechxy, gdeeprechxy) + call LIS_gather_masterproc_2d_local_to_global(n, rechxy, grechxy) + call LIS_gather_masterproc_2d_local_to_global(n, qslatxy, gqslatxy) + call LIS_gather_masterproc_2d_local_to_global(n, qrfsxy, gqrfsxy) + call LIS_gather_masterproc_2d_local_to_global(n, qspringsxy, gqspringsxy) + call LIS_gather_masterproc_2d_local_to_global(n, rechclim, grechclim) + call LIS_gather_masterproc_2d_local_to_global(n, smois(:,1,:), gsmois(:,1,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smois(:,2,:), gsmois(:,2,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smois(:,3,:), gsmois(:,3,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smois(:,4,:), gsmois(:,4,:)) + call LIS_gather_masterproc_2d_local_to_global(n, sh2o(:,1,:), gsh2o(:,1,:)) + call LIS_gather_masterproc_2d_local_to_global(n, sh2o(:,2,:), gsh2o(:,2,:)) + call LIS_gather_masterproc_2d_local_to_global(n, sh2o(:,3,:), gsh2o(:,3,:)) + call LIS_gather_masterproc_2d_local_to_global(n, sh2o(:,4,:), gsh2o(:,4,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smoiseq(:,1,:), gsmoiseq(:,1,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smoiseq(:,2,:), gsmoiseq(:,2,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smoiseq(:,3,:), gsmoiseq(:,3,:)) + call LIS_gather_masterproc_2d_local_to_global(n, smoiseq(:,4,:), gsmoiseq(:,4,:)) + + if ( LIS_masterproc ) then + ! Allocate out variables + write(LIS_logunit,*) "[INFO] Allocating Output Arrays for MMF Init." + allocate(grivercond(LIS_rc%gnc(n), LIS_rc%gnr(n))) + else + ! Allocate dummy "global" output variables + write(LIS_logunit,*) "[INFO] Allocating Dummy Output Arrays for MMF Init." + allocate(grivercond(1,1)) + endif + + if ( LIS_masterproc ) then + write(LIS_logunit,*) "[INFO] Initialization Groundwater for MMF." + + !write(*,*) "mmf_start GROUNDWATER VARIABLES: " + !do i = 50,60 + ! write(*,*) "EQWTD(",i,",35) = ",geqwtd(i,35) + !enddo + !do j = 30,40 + ! write(*,*) "EQWTD(55,",j,") = ",geqwtd(55,j) + !enddo + !do i = 50,60 + ! write(*,*) "ISLTYP(",i,",35) = ",gisltyp(i,35) + !enddo + !do j = 30,40 + ! write(*,*) "ISLTYP(55,",j,") = ",gisltyp(55,j) + !enddo + + call groundwater_init (noahmp401_struc(n)%nsoil, & !nsoil , + noahmp401_struc(n)%sldpth, & !dzs, + gisltyp, givgtyp, wtddt , & + gfdepth, gtopo, griverbed, geqwtd, grivercond, gpexp , garea ,gwtd , & + gsmois,gsh2o, gsmoiseq, gsmcwtdxy, gdeeprechxy, grechxy , & + gqslatxy, gqrfsxy, gqspringsxy, & + grechclim , & + 1, & !ids, + LIS_rc%gnc(n), & !ide, +1 for test + 1, & !jds, + LIS_rc%gnr(n), & !jde, + 1, 1, & !kds,kde, + 1, & !ims, + LIS_rc%gnc(n), & !ime, + 1, & !jms, + LIS_rc%gnr(n), & !jme, + 1, 1, & !kms,kme, + 1, & !ips, + LIS_rc%gnc(n), & !ipe, + 1, & !jps, + LIS_rc%gnr(n), & !jpe, + 1,1, & !kps,kpe, + 1, & !its, + LIS_rc%gnc(n), & !ite, + 1, & !jts, + LIS_rc%gnr(n), & !jte, + 1,1) !kts,kte + endif + write(LIS_logunit,*) "[INFO] Waiting for MMF Processors." + call MPI_Barrier(LIS_MPI_COMM, ierr) + + ! Scatter out and inout variables + write(LIS_logunit,*) "[INFO] Scattering MMF Init. Arrays." + call LIS_scatter_global_to_local_grid(n, gpexp, pexp) + call LIS_scatter_global_to_local_grid(n, geqwtd, eqwtd) + call LIS_scatter_global_to_local_grid(n, griverbed, riverbed) + call LIS_scatter_global_to_local_grid(n, gwtd, wtd) + call LIS_scatter_global_to_local_grid(n, gsmcwtdxy, smcwtdxy) + call LIS_scatter_global_to_local_grid(n, gdeeprechxy, deeprechxy) + call LIS_scatter_global_to_local_grid(n, grechxy, rechxy) + call LIS_scatter_global_to_local_grid(n, gqslatxy, qslatxy) + call LIS_scatter_global_to_local_grid(n, gqrfsxy, qrfsxy) + call LIS_scatter_global_to_local_grid(n, gqspringsxy, qspringsxy) + call LIS_scatter_global_to_local_grid(n, gsmois(:,1,:), smois(:,1,:)) + call LIS_scatter_global_to_local_grid(n, gsmois(:,2,:), smois(:,2,:)) + call LIS_scatter_global_to_local_grid(n, gsmois(:,3,:), smois(:,3,:)) + call LIS_scatter_global_to_local_grid(n, gsmois(:,4,:), smois(:,4,:)) + call LIS_scatter_global_to_local_grid(n, gsh2o(:,1,:), sh2o(:,1,:)) + call LIS_scatter_global_to_local_grid(n, gsh2o(:,2,:), sh2o(:,2,:)) + call LIS_scatter_global_to_local_grid(n, gsh2o(:,3,:), sh2o(:,3,:)) + call LIS_scatter_global_to_local_grid(n, gsh2o(:,4,:), sh2o(:,4,:)) + call LIS_scatter_global_to_local_grid(n, gsmoiseq(:,1,:), smoiseq(:,1,:)) + call LIS_scatter_global_to_local_grid(n, gsmoiseq(:,2,:), smoiseq(:,2,:)) + call LIS_scatter_global_to_local_grid(n, gsmoiseq(:,3,:), smoiseq(:,3,:)) + call LIS_scatter_global_to_local_grid(n, gsmoiseq(:,4,:), smoiseq(:,4,:)) + call LIS_scatter_global_to_local_grid(n, grivercond, rivercond) + + ! Deallocate temporary global variables + deallocate(gisltyp) + deallocate(givgtyp) + deallocate(gfdepth) + deallocate(gtopo) + deallocate(griverbed) + deallocate(geqwtd) + deallocate(grivercond) + deallocate(gpexp) + deallocate(garea) + deallocate(gwtd) + deallocate(gsmcwtdxy) + deallocate(gdeeprechxy) + deallocate(grechxy) + deallocate(gqslatxy) + deallocate(gqrfsxy) + deallocate(gqspringsxy) + deallocate(grechclim) + deallocate(gsmois) + deallocate(gsh2o) + deallocate(gsmoiseq) + +#else + + !write(*,*) "mmf_start GROUNDWATER VARIABLES: " + !do i = 50,60 + ! write(*,*) "EQWTD(",i,",35) = ",eqwtd(i,35) + !enddo + !do j = 30,40 + ! write(*,*) "EQWTD(55,",j,") = ",eqwtd(55,j) + !enddo + !do i = 50,60 + ! write(*,*) "ISLTYP(",i,",35) = ",isltyp(i,35) + !enddo + !do j = 30,40 + ! write(*,*) "ISLTYP(55,",j,") = ",isltyp(55,j) + !enddo + + call groundwater_init (noahmp401_struc(n)%nsoil, & !nsoil , + noahmp401_struc(n)%sldpth, & !dzs, + isltyp, ivgtyp, wtddt , & + fdepth, topo, riverbed, eqwtd, rivercond, pexp , area ,wtd , & + smois,sh2o, smoiseq, smcwtdxy, deeprechxy, rechxy , & + qslatxy, qrfsxy, qspringsxy, & + rechclim , & + NOAHMP401_struc(n)%col_min, & !ids, + NOAHMP401_struc(n)%col_max, & !ide, +1 for test + NOAHMP401_struc(n)%row_min, & !jds, + NOAHMP401_struc(n)%row_max, & !jde, + 1, 1, & !kds,kde, + NOAHMP401_struc(n)%col_min, & !ims, + NOAHMP401_struc(n)%col_max, & !ime, + NOAHMP401_struc(n)%row_min, & !jms, + NOAHMP401_struc(n)%row_max, & !jme, + 1, 1, & !kms,kme, + NOAHMP401_struc(n)%col_min, & !ips, + NOAHMP401_struc(n)%col_max, & !ipe, + NOAHMP401_struc(n)%row_min, & !jps, + NOAHMP401_struc(n)%row_max, & !jpe, + 1,1, & !kps,kpe, + NOAHMP401_struc(n)%col_min, & !its, + NOAHMP401_struc(n)%col_max, & !ite, + NOAHMP401_struc(n)%row_min, & !jts, + NOAHMP401_struc(n)%row_max, & !jte, + 1,1) !kts,kte +#endif + !print*, 'mmf_start GROUNDWATER VARIABLES: ' + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) + !print*, 'SMCWTD = ',smcwtdxy(18,12) + !print*, 'EQWTD = ',eqwtd(18,12) + !print*, 'WTD = ',wtd(18,12) + + do row=NOAHMP401_struc(n)%row_min, NOAHMP401_struc(n)%row_max + do col=NOAHMP401_struc(n)%col_min, NOAHMP401_struc(n)%col_max + t = NOAHMP401_struc(n)%rct_idx(col,row) + !print*, col + !print*, row + !print*, wtd(col,row) + !print*, t + ! Added if statement to deal with no-data values. TML + if(t .ne. LIS_rc%udef) then + NOAHMP401_struc(n)%noahmp401(t)%wtd = wtd(col,row) + NOAHMP401_struc(n)%noahmp401(t)%zwt = wtd(col,row) !!!! zwt should be the same as wtd + NOAHMP401_struc(n)%noahmp401(t)%rivercond = rivercond(col,row) + NOAHMP401_struc(n)%noahmp401(t)%smc(:) = smois(col,:,row) ! smois + NOAHMP401_struc(n)%noahmp401(t)%sh2o(:) = sh2o(col,:,row) + NOAHMP401_struc(n)%noahmp401(t)%smoiseq(:)= smoiseq(col,:,row) + if(isltyp(col,row) .eq. 14) then + NOAHMP401_struc(n)%noahmp401(t)%smcwtd = 1.0 + else + NOAHMP401_struc(n)%noahmp401(t)%smcwtd = smcwtdxy(col,row) + endif + NOAHMP401_struc(n)%noahmp401(t)%deeprech = deeprechxy(col,row) + NOAHMP401_struc(n)%noahmp401(t)%rech = rechxy(col,row) + NOAHMP401_struc(n)%noahmp401(t)%qslat = qslatxy(col,row) + NOAHMP401_struc(n)%noahmp401(t)%qrfs = qrfsxy(col,row) + NOAHMP401_struc(n)%noahmp401(t)%qsprings = qspringsxy(col,row) + endif + NOAHMP401_struc(n)%rivercond(col, row) = rivercond(col,row) !!! make a copy to the 2D paramter data structure + NOAHMP401_struc(n)%riverbed(col, row) = riverbed(col,row) !!! make a copy to the 2D paramter data structure + NOAHMP401_struc(n)%eqwtd(col, row) = eqwtd(col,row) !!! make a copy + enddo + enddo + deallocate(isltyp) + deallocate(ivgtyp) + deallocate(fdepth) + deallocate(topo) + deallocate(area) + deallocate(rechclim) + deallocate(rivercond) + deallocate(wtd) + deallocate(riverbed) + deallocate(eqwtd) + deallocate(pexp) + deallocate(smcwtdxy) + deallocate(deeprechxy) + deallocate(rechxy) + deallocate(qslatxy) + deallocate(qrfsxy) + deallocate(qspringsxy) + deallocate(smois) + deallocate(sh2o) + deallocate(smoiseq) +end subroutine mmf_start diff --git a/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 index 715264d19..cfa7f92b3 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 @@ -53,7 +53,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & !ag (12Sep2019) rivsto, fldsto, fldfrc,& parameters , & ! out Noah MP only - sfcheadrt , INFXSRT, soldrain) ! For WRF-Hydro + sfcheadrt , INFXSRT, soldrain, printdebug) ! For WRF-Hydro use module_sf_noahmpdrv_401, only: noahmplsm_401 use module_sf_noahmplsm_401 @@ -74,6 +74,8 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & integer, intent(in) :: minute ! minute of the current time step [-] real, intent(in) :: dt ! timestep [s] + integer, intent(in) :: printdebug !TML: Debuggging flag + ! Revised by Zhuo Wang and Shugong Wang real, intent(in) :: dz8w ! thickness of atmo layers [m] @@ -788,7 +790,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & #endif ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte) + its,ite, jts,jte, kts,kte, printdebug, loctime) ! TML: Added debugging term. ! Added by Zhuo Wang and Shugong on 10/30/2018 tsk = tskinout(1,1) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_groundwater_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_groundwater_401.F90 index b49f06979..4cf437bdc 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_groundwater_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_groundwater_401.F90 @@ -21,6 +21,7 @@ SUBROUTINE WTABLE_mmf_noahmp (NSOIL ,XLAND ,XICE ,XICE_THRESHOLD ,ISI ! ---------------------------------------------------------------------- USE NOAHMP_TABLES_401, ONLY: BEXP_TABLE, DKSAT_TABLE, SMCMAX_TABLE,PSISAT_TABLE, SMCWLT_TABLE + use LIS_logMod, only : LIS_logunit ! ---------------------------------------------------------------------- IMPLICIT NONE ! ---------------------------------------------------------------------- @@ -87,9 +88,12 @@ SUBROUTINE WTABLE_mmf_noahmp (NSOIL ,XLAND ,XICE ,XICE_THRESHOLD ,ISI INTEGER, DIMENSION( ims:ime, jms:jme ) :: LANDMASK !-1 for water (ice or no ice) and glacial areas, 1 for land where the LSM does its soil moisture calculations. REAL :: BEXP,DKSAT,PSISAT,SMCMAX,SMCWLT + !REAL :: QLAT_temp,QRF_temp,QSPRING_temp,DEEPRECH_temp + !print *,"WTABLE Called" + !print*, 'ZWT = ',WTD(25,12) DELTAT = WTDDT * 60. !timestep in seconds for this calculation - + ZSOIL(0) = 0. ZSOIL(1) = -DZS(1) DO K = 2, NSOIL @@ -109,7 +113,8 @@ SUBROUTINE WTABLE_mmf_noahmp (NSOIL ,XLAND ,XICE ,XICE_THRESHOLD ,ISI ,ids,ide,jds,jde,kds,kde & ,ims,ime,jms,jme,kms,kme & ,its,ite,jts,jte,kts,kte ) - + !print*, 'LATERALFLOW COMPLETE' + !print*, 'ZWT = ',WTD(25,12) !compute flux from grounwater to rivers in the cell @@ -171,29 +176,118 @@ SUBROUTINE WTABLE_mmf_noahmp (NSOIL ,XLAND ,XICE ,XICE_THRESHOLD ,ISI SH2O(1:NSOIL) = SH2OXY(I,1:NSOIL,J) SMCEQ(1:NSOIL) = SMOISEQ(I,1:NSOIL,J) + !IF (I == 25) THEN + ! IF (J == 12) THEN + ! print *, "CALLING UPDATEWTD" + ! print*, 'NSOIL = ',NSOIL + ! print*, 'DZS = ',DZS + ! print*, 'ZSOIL = ',ZSOIL + ! print*, 'SMCMAX = ',SMCMAX + ! print*, 'SMCWLT = ',SMCWLT + ! print*, 'PSISAT = ',PSISAT + !! print*, 'BEXP = ',BEXP + ! print*, 'TOTWATER = ',TOTWATER + ! print*, 'WTD(I,J) = ',WTD(I,J) + ! print*, 'SMC = ',SMC + ! print*, 'SH2O = ',SH2O + ! print*, 'SMCWTD = ',SMCWTD(I,J) + ! ENDIF + !ENDIF + + !Update the water table depth and soil moisture CALL UPDATEWTD ( NSOIL, DZS , ZSOIL, SMCEQ, SMCMAX, SMCWLT, PSISAT, BEXP ,I , J , &!in TOTWATER, WTD(I,J), SMC, SH2O, SMCWTD(I,J) , &!inout QSPRING(I,J) ) !out + !IF (I == 25) THEN + ! IF (J == 12) THEN + ! print *, "FINISHED UPDATEWTD" + ! print*, 'TOTWATER = ',TOTWATER + ! print*, 'WTD(I,J) = ',WTD(I,J) + ! print*, 'SMC = ',SMC + ! print*, 'SH2O = ',SH2O + ! print*, 'SMCWTD = ',SMCWTD(I,J) + ! print*, 'QSPRING = ',QSPRING(I,J) + ! ENDIF + !ENDIF + + !now update soil moisture SMOIS(I,1:NSOIL,J) = SMC(1:NSOIL) SH2OXY(I,1:NSOIL,J) = SH2O(1:NSOIL) + ELSE + ! Set Cumulative Variables to Zero over Water (i.e. LANDMASK < 0) + QLAT(I,J) = 0. + QRF(I,J) = 0. + QSPRING(I,J) = 0. + DEEPRECH(I,J) = 0. + ENDIF ENDDO ENDDO + !print*, 'UPDATEWTD COMPLETE' + !print*, 'ZWT = ',WTD(25,12) + !accumulate fluxes for output DO J=jts,jte DO I=its,ite - QSLAT(I,J) = QSLAT(I,J) + QLAT(I,J)*1.E3 - QRFS(I,J) = QRFS(I,J) + QRF(I,J)*1.E3 - QSPRINGS(I,J) = QSPRINGS(I,J) + QSPRING(I,J)*1.E3 - RECH(I,J) = RECH(I,J) + DEEPRECH(I,J)*1.E3 + !***** TML: Previous Debugging Lines to Find Floating Pt. Errors ***** + ! Remove This Later ... + !TML Round Down Small Values to Prevent Floating Point Errors + !IF (ABS(QLAT(I,J)).LT.1.E-08)THEN + ! QLAT(I,J) = 0. + !ENDIF + !IF (ABS(QRF(I,J)).LT.1.E-05)THEN + ! QRF(I,J) = 0. + !ENDIF + !IF (ABS(QSPRING(I,J)).LT.1.E-08)THEN + ! QSPRING(I,J) = 0. + !ENDIF + !IF (ABS(DEEPRECH(I,J)).LT.1.E-08)THEN + ! DEEPRECH(I,J) = 0. + !ENDIF + !write(LIS_logunit,*) 'QRF = ',QRF(I,J) + !QLAT_temp = QLAT(I,J)*1.E3 + !QRF_temp = QRF(I,J)*1.E3 + !QSPRING_temp = QSPRING(I,J)*1.E3 + !DEEPRECH_temp = DEEPRECH(I,J)*1.E3 + !Set Excessively large values to 0 (issue for parallel runs) + !IF (ABS(QLAT_temp).GT.1.E+01)THEN + ! print*, 'Warning: QLAT = ',QLAT_temp + ! QLAT_temp = 0. + !ENDIF + !IF (ABS(QRF_temp).GT.1.E+01)THEN + ! print*, 'Warning: QRF = ',QRF_temp + ! QRF_temp = 0. + !ENDIF + !IF (ABS(QSPRING_temp).GT.1.E+01)THEN + ! print*, 'Warning: QSPRING = ',QSPRING_temp + ! QSPRING_temp = 0. + !ENDIF + !IF (ABS(DEEPRECH_temp).GT.1.E+01)THEN + ! print*, 'Warning: DEEPRECH = ',DEEPRECH_temp + ! DEEPRECH_temp = 0. + !ENDIF + !write(LIS_logunit,*) 'DEEPRECH_temp = ',DEEPRECH_temp + !write(LIS_logunit,*) 'RECH = ',RECH(I,J) + !QSLAT(I,J) = QSLAT(I,J) + QLAT_temp + !QRFS(I,J) = QRFS(I,J) + QRF_temp + !QSPRINGS(I,J) = QSPRINGS(I,J) + QSPRING_temp + !RECH(I,J) = RECH(I,J) + DEEPRECH_temp + !***** TML: END DEBUGGING CODE; Returns to Original Code ***** + IF(LANDMASK(I,J).GT.0)THEN + QSLAT(I,J) = QSLAT(I,J) + QLAT(I,J)*1.E3 + QRFS(I,J) = QRFS(I,J) + QRF(I,J)*1.E3 + QSPRINGS(I,J) = QSPRINGS(I,J) + QSPRING(I,J)*1.E3 + RECH(I,J) = RECH(I,J) + DEEPRECH(I,J)*1.E3 + ENDIF !zero out DEEPRECH DEEPRECH(I,J) =0. + ENDDO ENDDO @@ -232,17 +326,22 @@ SUBROUTINE LATERALFLOW (ISLTYP,WTD,QLAT,FDEPTH,TOPO,LANDMASK,DELTAT,AREA & REAL, PARAMETER :: PI = 3.14159265 REAL, PARAMETER :: FANGLE = 0.22754493 ! = 0.5*sqrt(0.5*tan(pi/8)) -itsh=max(its-1,ids) -iteh=min(ite+1,ide-1) -jtsh=max(jts-1,jds) -jteh=min(jte+1,jde-1) +!itsh=max(its-1,ids) +!iteh=min(ite+1,ide-1) +!jtsh=max(jts-1,jds) +!jteh=min(jte+1,jde-1) +! TML: Exchanges not possible between boundaries; need to rely on halos ... +itsh=its +iteh=ite +jtsh=jts +jteh=jte DO J=jtsh,jteh DO I=itsh,iteh IF(FDEPTH(I,J).GT.0.)THEN KLAT = DKSAT_TABLE(ISLTYP(I,J)) * KLATFACTOR(ISLTYP(I,J)) - IF(WTD(I,J) < -1.5)THEN + IF(WTD(I,J) .LT. -1.5)THEN KCELL(I,J) = FDEPTH(I,J) * KLAT * EXP( (WTD(I,J) + 1.5) / FDEPTH(I,J) ) ELSE KCELL(I,J) = KLAT * ( WTD(I,J) + 1.5 + FDEPTH(I,J) ) @@ -255,10 +354,16 @@ SUBROUTINE LATERALFLOW (ISLTYP,WTD,QLAT,FDEPTH,TOPO,LANDMASK,DELTAT,AREA & ENDDO ENDDO -itsh=max(its,ids+1) -iteh=min(ite,ide-2) -jtsh=max(jts,jds+1) -jteh=min(jte,jde-2) +!itsh=max(its,ids+1) +!iteh=min(ite,ide-2) +!jtsh=max(jts,jds+1) +!jteh=min(jte,jde-2) + +! TML: Exchanges not possible between boundaries; need to rely on halos ... +itsh=its+1 +iteh=ite-1 +jtsh=jts+1 +jteh=jte-1 DO J=jtsh,jteh DO I=itsh,iteh @@ -439,6 +544,11 @@ SUBROUTINE UPDATEWTD (NSOIL, DZS, ZSOIL ,SMCEQ ,& !in if(totwater.le.maxwatup)then wtd = wtd + totwater/(smcmax-smcwtd) totwater=0. + !IF (ILOC == 25) THEN + ! IF (JLOC == 12) THEN + ! print *, "wtd updated" + ! ENDIF + !ENDIF else totwater=totwater-maxwatup wtd=zsoil(nsoil)-dzs(nsoil) @@ -605,7 +715,10 @@ SUBROUTINE UPDATEWTD (NSOIL, DZS, ZSOIL ,SMCEQ ,& !in ENDIF SH2O = SMC - SICE - + + !if(qspring.lt.-1.E30)then + ! print *, 'WARNING EXTREME GW EXCHANGE; Floating point error likely' + !endif END SUBROUTINE UPDATEWTD diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 index 5a2b25489..9e85cd7f4 100644 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 @@ -56,7 +56,7 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, #endif ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & - its,ite, jts,jte, kts,kte, & + its,ite, jts,jte, kts,kte, PRINTDEBUG, loctime, & MP_RAINC, MP_RAINNC, MP_SHCV, MP_SNOW, MP_GRAUP, MP_HAIL ) !---------------------------------------------------------------- USE MODULE_SF_NOAHMPLSM_401 @@ -74,6 +74,8 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, ! Added LIS undefined value as an input - David Mocko REAL, INTENT(IN ) :: LIS_undef_value + INTEGER, INTENT(IN ) :: PRINTDEBUG ! TML: Debugging Flag + real*8, INTENT(IN ) :: loctime INTEGER, INTENT(IN ) :: ITIMESTEP ! timestep number INTEGER, INTENT(IN ) :: YR ! 4-digit year REAL, INTENT(IN ) :: JULIAN ! Julian day @@ -783,6 +785,158 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, IF(VEGTYP == 27) FVEG = 0.0 IF(VEGTYP == 27) PLAI = 0.0 + !print*, 'I= ',I + !print*, 'J= ',J + + if (PRINTDEBUG .eq. 10) then + !print*, 'START VARIABLES: ' + !print*, 'TIME= ',loctime + !print*, 'XLAT = ', XLAT(I,J) + !print*, 'XLONG = ',XLONG(I,J) + !print*, 'T_ML = ',T_ML + !print*, 'Q_ML = ',Q_ML + !print*, 'U_ML = ',U_ML + !print*, 'V_ML = ',V_ML + !print*, 'SWDN = ',SWDN + !print*, 'LWDN = ',LWDN + !print*, 'P_ML = ',P_ML + !print*, 'PSFC = ',PSFC + !print*, 'PRCP = ',PRCP + !print*, 'TBOT =',TBOT + !print*, 'FOLN =',FOLN + !print*, 'FICEOLD =',FICEOLD + !print*, 'Z_ML =',Z_ML + !print*, 'ALBOLD =',ALBOLD + !print*, 'SNEQVO =',SNEQVO + !print*, 'SMCWTD = ',SMCWTD + !print*, 'ZWT = ',ZWT + !print*, 'SMC-1 = ',SMC(1) + !print*, 'SMC-2 = ',SMC(2) + !print*, 'SMC-3 = ',SMC(3) + !print*, 'SMC-4 = ',SMC(4) + !print*, 'STC-1 = ',STC(1) + !print*, 'STC-2 = ',STC(2) + !print*, 'STC-3 = ',STC(3) + !print*, 'STC-4 = ',STC(4) + !print*, 'SMH2O-1 =',SMH2O(1) + !print*, 'SMH2O-2 =',SMH2O(2) + !print*, 'SMH2O-3 =',SMH2O(3) + !print*, 'SMH2O-4 =',SMH2O(4) + !print*, 'TAH =',TAH + !print*, 'EAH =',EAH + !print*, 'FWET =',FWET + !print*, 'CANLIQ =',CANLIQ + !print*, 'CANICE =',CANICE + !print*, 'TV =',TV + !print*, 'TG =',TG + !print*, 'QSFC1D =',QSFC1D + !print*, 'QSNOW =',QSNOW + !print*, 'ISNOW =',ISNOW + !print*, 'ZSNSO =',ZSNSO + !print*, 'SNDPTH =',SNDPTH + !print*, 'SWE =',SWE + !print*, 'SNICE =',SNICE + !print*, 'SNLIQ =',SNLIQ + !print*, 'ZWT =',ZWT + !print*, 'WA =',WA + !print*, 'WT =',WT + !print*, 'WSLAKE =',WSLAKE + !print*, 'LFMASS =',LFMASS + !print*, 'RTMASS =',RTMASS + !print*, 'STMASS =',STMASS + !print*, 'WOOD =',WOOD + !print*, 'STBLCP =',STBLCP + !print*, 'FASTCP =',FASTCP + !print*, 'PLAI =',PLAI + !print*, 'PSAI =',PSAI + !print*, 'CM =',CM + !print*, 'CH =',CH + !print*, 'TAUSS =',TAUSS + !print*, 'GRAIN =',GRAIN + !print*, 'GDD =',GDD + !print*, 'PGS =',PGS + !print*, 'SMCWTD =',SMCWTD + !print*, 'DEEPRECH = ',DEEPRECH + !print*, 'RECH = ',RECH + !print*, 'GECROS1D = ',GECROS1D + !print*, 'Z0WRF =',Z0WRF + !print*, 'FSA =',FSA + !print*, 'FSR =',FSR + !print*, 'FIRA =',FIRA + !print*, 'FSH =',FSH + !print*, 'SSOIL =',SSOIL + !print*, 'FCEV =',FCEV + !print*, 'FGEV =',FGEV + !print*, 'FCTR =',FCTR + !print*, 'ECAN =',ECAN + !print*, 'ETRAN =',ETRAN + !print*, 'ESOIL =',ESOIL + !print*, 'TRAD =',TRAD + !print*, 'SUBSNOW =',SUBSNOW + !print*, 'RELSMC =',RELSMC + !print*, 'TGB =',TGB + !print*, 'TGV =',TGV + !print*, 'T2MV =',T2MV + !print*, 'T2MB =',T2MB + !print*, 'Q2MV =',Q2MV + !print*, 'Q2MB =',Q2MB + !print*, 'RUNSF =',RUNSF + !print*, 'RUNSB =',RUNSB + !print*, 'APAR =',APAR + !print*, 'PSN =',PSN + !print*, 'SAV =',SAV + !print*, 'SAG =',SAG + !print*, 'FSNO =',FSNO + !print*, 'NEE =',NEE + !print*, 'GPP =',GPP + !print*, 'NPP =',NPP + !print*, 'FVEGMP =',FVEGMP + !print*, 'SALB =',SALB + !print*, 'QSNBOT =',QSNBOT + !print*, 'PONDING =',PONDING + !print*, 'PONDING1 =',PONDING1 + !print*, 'PONDING2 =',PONDING2 + !print*, 'RSSUN =',RSSUN + !print*, 'RSSHA =',RSSHA + !print*, 'BGAP =',BGAP + !print*, 'WGAP =',WGAP + !print*, 'CHV =',CHV + !print*, 'CHB =',CHB + !print*, 'EMISSI =',EMISSI + !print*, 'SHG =',SHG + !print*, 'SHC =',SHC + !print*, 'SHB =',SHB + !print*, 'EVG =',EVG + !print*, 'EVB =',EVB + !print*, 'GHV =',GHV + !print*, 'GHB =',GHB + !print*, 'IRG =',IRG + !print*, 'IRC =',IRC + !print*, 'IRB =',IRB + !print*, 'TR =',TR + !print*, 'EVC =',EVC + !print*, 'FGEV_PET =',FGEV_PET + !print*, 'FCEV_PET =',FCEV_PET + !print*, 'FCTR_PET =',FCTR_PET + !print*, 'CHLEAF =',CHLEAF + !print*, 'CHUC =',CHUC + !print*, 'CHV2 =',CHV2 + !print*, 'CHB2 =',CHB2 + !print*, 'FPICE =',FPICE + !print*, 'PAHV =',PAHV + !print*, 'PAHG =',PAHG + !print*, 'PAHB =',PAHB + !print*, 'PAH =',PAH + !print*, 'LAISUN =',LAISUN + !print*, 'LAISHA =',LAISHA + !print*, 'RB =',RB + !print*, 'RIVSTO =',rivsto + !print*, 'FLDSTO =',fldsto + !print*, 'FLDFRC =',fldfrc + !print*, 'SFCHEADRT(i,j) =',sfcheadrt(i,j) + endif + + IF ( VEGTYP == ISICE_TABLE ) THEN ICE = -1 ! Land-ice point CALL NOAHMP_OPTIONS_GLACIER(IOPT_ALB ,IOPT_SNF ,IOPT_TBOT, IOPT_STC, IOPT_GLA , IOPT_SNDPTH) @@ -929,6 +1083,154 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, INFXSRT(i,j) = RUNSF*dt !mm , surface runoff #endif + ! Check of latitude is 0, effectively comments out this block... TML + if (PRINTDEBUG .eq. 10) then + !print*, 'END VARIABLES: ' + !print*, 'TIME= ',loctime + !print*, 'XLAT = ', XLAT(I,J) + !print*, 'XLONG = ',XLONG(I,J) + !print*, 'T_ML = ',T_ML + !print*, 'Q_ML = ',Q_ML + !print*, 'U_ML = ',U_ML + !print*, 'V_ML = ',V_ML + !print*, 'SWDN = ',SWDN + !print*, 'LWDN = ',LWDN + !print*, 'P_ML = ',P_ML + !print*, 'PSFC = ',PSFC + !print*, 'PRCP = ',PRCP + !print*, 'TBOT =',TBOT + !print*, 'FOLN =',FOLN + !print*, 'FICEOLD =',FICEOLD + !print*, 'Z_ML =',Z_ML + !print*, 'ALBOLD =',ALBOLD + !print*, 'SNEQVO =',SNEQVO + !print*, 'SMCWTD = ',SMCWTD + !print*, 'ZWT = ',ZWT + !print*, 'SMC-1 = ',SMC(1) + !print*, 'SMC-2 = ',SMC(2) + !print*, 'SMC-3 = ',SMC(3) + !print*, 'SMC-4 = ',SMC(4) + !print*, 'STC-1 = ',STC(1) + !print*, 'STC-2 = ',STC(2) + !print*, 'STC-3 = ',STC(3) + !print*, 'STC-4 = ',STC(4) + !print*, 'SMH2O-1 =',SMH2O(1) + !print*, 'SMH2O-2 =',SMH2O(2) + !print*, 'SMH2O-3 =',SMH2O(3) + !print*, 'SMH2O-4 =',SMH2O(4) + !print*, 'TAH =',TAH + !print*, 'EAH =',EAH + !print*, 'FWET =',FWET + !print*, 'CANLIQ =',CANLIQ + !print*, 'CANICE =',CANICE + !print*, 'TV =',TV + !print*, 'TG =',TG + !print*, 'QSFC1D =',QSFC1D + !print*, 'QSNOW =',QSNOW + !print*, 'ISNOW =',ISNOW + !print*, 'ZSNSO =',ZSNSO + !print*, 'SNDPTH =',SNDPTH + !print*, 'SWE =',SWE + !print*, 'SNICE =',SNICE + !print*, 'SNLIQ =',SNLIQ + !print*, 'ZWT =',ZWT + !print*, 'WA =',WA + !print*, 'WT =',WT + !print*, 'WSLAKE =',WSLAKE + !print*, 'LFMASS =',LFMASS + !print*, 'RTMASS =',RTMASS + !print*, 'STMASS =',STMASS + !print*, 'WOOD =',WOOD + !print*, 'STBLCP =',STBLCP + !print*, 'FASTCP =',FASTCP + !print*, 'PLAI =',PLAI + !print*, 'PSAI =',PSAI + !print*, 'CM =',CM + !print*, 'CH =',CH + !print*, 'TAUSS =',TAUSS + !print*, 'GRAIN =',GRAIN + !print*, 'GDD =',GDD + !print*, 'PGS =',PGS + !print*, 'SMCWTD =',SMCWTD + !print*, 'DEEPRECH = ',DEEPRECH + !print*, 'RECH = ',RECH + !print*, 'GECROS1D = ',GECROS1D + !print*, 'Z0WRF =',Z0WRF + !print*, 'FSA =',FSA + !print*, 'FSR =',FSR + !print*, 'FIRA =',FIRA + !print*, 'FSH =',FSH + !print*, 'SSOIL =',SSOIL + !print*, 'FCEV =',FCEV + !print*, 'FGEV =',FGEV + !print*, 'FCTR =',FCTR + !print*, 'ECAN =',ECAN + !print*, 'ETRAN =',ETRAN + !print*, 'ESOIL =',ESOIL + !print*, 'TRAD =',TRAD + !print*, 'SUBSNOW =',SUBSNOW + !print*, 'RELSMC =',RELSMC + !print*, 'TGB =',TGB + !print*, 'TGV =',TGV + !print*, 'T2MV =',T2MV + !print*, 'T2MB =',T2MB + !print*, 'Q2MV =',Q2MV + !print*, 'Q2MB =',Q2MB + !print*, 'RUNSF =',RUNSF + !print*, 'RUNSB =',RUNSB + !print*, 'APAR =',APAR + !print*, 'PSN =',PSN + !print*, 'SAV =',SAV + !print*, 'SAG =',SAG + !print*, 'FSNO =',FSNO + !print*, 'NEE =',NEE + !print*, 'GPP =',GPP + !print*, 'NPP =',NPP + !print*, 'FVEGMP =',FVEGMP + !print*, 'SALB =',SALB + !print*, 'QSNBOT =',QSNBOT + !print*, 'PONDING =',PONDING + !print*, 'PONDING1 =',PONDING1 + !print*, 'PONDING2 =',PONDING2 + !print*, 'RSSUN =',RSSUN + !print*, 'RSSHA =',RSSHA + !print*, 'BGAP =',BGAP + !print*, 'WGAP =',WGAP + !print*, 'CHV =',CHV + !print*, 'CHB =',CHB + !print*, 'EMISSI =',EMISSI + !print*, 'SHG =',SHG + !print*, 'SHC =',SHC + !print*, 'SHB =',SHB + !print*, 'EVG =',EVG + !print*, 'EVB =',EVB + !print*, 'GHV =',GHV + !print*, 'GHB =',GHB + !print*, 'IRG =',IRG + !print*, 'IRC =',IRC + !print*, 'IRB =',IRB + !print*, 'TR =',TR + !print*, 'EVC =',EVC + !print*, 'FGEV_PET =',FGEV_PET + !print*, 'FCEV_PET =',FCEV_PET + !print*, 'FCTR_PET =',FCTR_PET + !print*, 'CHLEAF =',CHLEAF + !print*, 'CHUC =',CHUC + !print*, 'CHV2 =',CHV2 + !print*, 'CHB2 =',CHB2 + !print*, 'FPICE =',FPICE + !print*, 'PAHV =',PAHV + !print*, 'PAHG =',PAHG + !print*, 'PAHB =',PAHB + !print*, 'PAH =',PAH + !print*, 'LAISUN =',LAISUN + !print*, 'LAISHA =',LAISHA + !print*, 'RB =',RB + !print*, 'RIVSTO =',rivsto + !print*, 'FLDSTO =',fldsto + !print*, 'FLDFRC =',fldfrc + !print*, 'SFCHEADRT(i,j) =',sfcheadrt(i,j) + endif ! INPUT/OUTPUT @@ -1715,7 +2017,7 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT else waxy (I,J) = 0. wtxy (I,J) = 0. - areaxy (I,J) = (DX * DY) / ( MSFTX(I,J) * MSFTY(I,J) ) + !areaxy (I,J) = (DX * DY) / ( MSFTX(I,J) * MSFTY(I,J) )! SW, we calculate areaxy outside endif IF(IVGTYP(I,J) == ISBARREN_TABLE .OR. IVGTYP(I,J) == ISICE_TABLE .OR. & @@ -2035,10 +2337,12 @@ SUBROUTINE GROUNDWATER_INIT ( & ZSOIL(NS) = ZSOIL(NS-1) - DZS(NS) END DO + ! SW,kludge solution to solve boundary divide by zero + !itf=min0(ite,ide-1) + !jtf=min0(jte,jde-1) - itf=min0(ite,ide-1) - jtf=min0(jte,jde-1) - + itf=min0(ite,ide) ! SW + jtf=min0(jte,jde) ! SW WHERE(IVGTYP.NE.ISWATER_TABLE.AND.IVGTYP.NE.ISICE_TABLE) LANDMASK=1 @@ -2057,6 +2361,25 @@ SUBROUTINE GROUNDWATER_INIT ( & NCOUNT=0 + !DO J=jts,jtf + ! DO I=its,itf + ! IF(WTD(I,J).EQ.-37.11)THEN + ! print*, 'WTD = ',WTD(I,J) + ! print*, 'I = ',I + ! print*, 'J = ',J + ! ENDIF + ! ENDDO + !ENDDO + + !print*, 'GROUNDWATER_INIT Initial Values' + !print*, 'WTD = ',WTD(18,12) + !print*, 'ISLTYP = ',ISLTYP(18,12) + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) + + DO NITER=1,500 !Calculate lateral flow @@ -2085,6 +2408,14 @@ SUBROUTINE GROUNDWATER_INIT ( & EQWTD=WTD + !print*, 'GROUNDWATER_INIT Lateral Flow (1)' + !print*, 'WTD = ',WTD(18,12) + !print*, 'ISLTYP = ',ISLTYP(18,12) + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) + !after adjusting, where qlat > 1cm/year now wtd is at the surface. !it may still happen that qlat + rech > 0 and eqwtd-rbed <0. There the wtd can !rise to the surface (poor drainage) but the et will then increase. @@ -2153,6 +2484,14 @@ SUBROUTINE GROUNDWATER_INIT ( & ENDDO ENDDO + !print*, 'GROUNDWATER_INIT Lateral Flow (2)' + !print*, 'WTD = ',WTD(18,12) + !print*, 'ISLTYP = ',ISLTYP(18,12) + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) + !now compute eq. soil moisture, change soil moisture to be compatible with the water table and compute deep soil moisture DO J = jts,jtf @@ -2198,6 +2537,8 @@ SUBROUTINE GROUNDWATER_INIT ( & DX = FUNC/DFUNC SMC = SMC - DX + SMC = MAX(SMC,1.E-4) !TML FIX TO PREVENT EXTREMELY LOW SMC + SMC = MIN(SMC,SMCMAX) !TML FIX TO PREVENT EXTREMELY HIGH SMC IF ( ABS (DX) < 1.E-6)EXIT ENDDO @@ -2245,7 +2586,13 @@ SUBROUTINE GROUNDWATER_INIT ( & ENDDO ENDDO - + !print*, 'GROUNDWATER_INIT Soil Update' + !print*, 'WTD = ',WTD(18,12) + !print*, 'ISLTYP = ',ISLTYP(18,12) + !print*, 'SMC-1 = ',smois(18,1,12) + !print*, 'SMC-2 = ',smois(18,2,12) + !print*, 'SMC-3 = ',smois(18,3,12) + !print*, 'SMC-4 = ',smois(18,4,12) END SUBROUTINE GROUNDWATER_INIT