From cc5cf91d884f256f83eb4867be30e9a0cf5fefb9 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Mon, 5 Jun 2023 16:41:29 -0400 Subject: [PATCH 1/6] snow update additions for SnowModel --- .../da_snow/noahmp401_snow_update.F90 | 81 ++++++++++++++++--- 1 file changed, 68 insertions(+), 13 deletions(-) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 index 93e105177..c1a660d54 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 @@ -16,6 +16,7 @@ ! 14 Dec 2018: Yeosang Yoon; Modified code for NoahMP 4.0.1 and SNODEP ! 15 May 2019: Yeosang Yoon; Modified for NoahMP 4.0.1 and LDTSI ! 13 Dec 2019: Eric Kemp; Replaced LDTSI with SNOW +! 05 Jun 2023: Justin Pflug; fixes for SnowModel-defined snow updates ! ! !INTERFACE @@ -59,6 +60,7 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) integer :: snl_idx,i,j,iz integer :: iloc, jloc ! needed, but not use real :: smp,sneqv,snowh + real :: snoden real :: sneqv1,snowh1 real :: ponding1,ponding2 integer :: newnode @@ -137,6 +139,22 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) stc(1:nsoil) = & noahmp401_struc(n)%noahmp401(t)%tslb(1:nsoil) + ! NMP snow density calculation + if(snowh.gt.0) then + snoden = sneqv/(snowh*1000) + else + snoden = 0.55 + endif + + ! allow snow update even in cases where changes opp. directions + ! alter snow depth change to be in direction of SWE change + if((dsneqv.gt.0.and.dsnowh.le.0).or.& + (dsneqv.lt.0.and.dsnowh.ge.0)) then + dsnowh = (dsneqv/snoden)/1000 + ! set snow depth change to zero in instance where no SWE change + elseif(dsneqv.eq.0.and.dsnowh.ne.0) then + dsnowh = 0. + endif ! from snowfall routine ! creating a new layer @@ -147,12 +165,12 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) newnode = 0 - if(isnow == 0 .and. snowh >= 0.025.and.& - (dsneqv.gt.0.and.dsnowh.gt.0)) then !mb: change limit - isnow = -1 - newnode = 1 + if(isnow.eq.0.and.dsneqv.gt.0.and.dsnowh.gt.0) then + if(snowh.ge.0.025) then + isnow = -1 + newnode = 1 + endif dzsnso(0)= snowh - snowh = 0. stc(0) = min(273.16, noahmp401_struc(n)%noahmp401(t)%sfctmp) snice(0) = sneqv snliq(0) = 0. @@ -168,23 +186,60 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) if(dsneqv.lt.0.and.dsnowh.lt.0) then snowh1 = snowh + dsnowh sneqv1 = sneqv + dsneqv +! if dsnowh adjusted since dsneqv and dsnowh in opp. directions +! can cause one or other snowh1 or sneqv1 to be negative + if(sneqv1.gt.0.and.snowh1.le.0) then + snowh = ((sneqv1/snoden)/1000)-dsnowh + snowh1 = snowh + dsnowh +! if SWE disappears, also make sure snow depth disappears + elseif(sneqv.le.0) then + sneqv = -dsneqv + sneqv1 = sneqv + dsneqv + snowh = -dsnowh + snowh1 = snowh + dsnowh + endif +! make sure snow layers currently exist in decrease case + if(dzsnso(0).eq.0) then + if(snowh.ge.0.025) then + isnow = -1 + else + isnow = 0 + endif + dzsnso(0)= snowh + stc(0) = min(273.16, noahmp401_struc(n)%noahmp401(t)%sfctmp) + snice(0) = sneqv + snliq(0) = 0. + endif if(snowh1.ge.0.and.sneqv1.ge.0) then snowh = snowh + dsnowh sneqv = sneqv + dsneqv -! update dzsnso -! how do you determine the thickness of a layer? +! snow can no longer fill layer 1 if(snowh.le.dzsnso(0)) then isnow = 0 dzsnso(-nsnow+1:(isnow-1)) = 0 dzsnso(isnow) = snowh + snice(-nsnow+1:(isnow-1)) = 0 + snice(isnow) = sneqv + snliq(-nsnow+1:isnow) = 0 +! snow can no longer fill layer 1 and 2 elseif(snowh.le.(dzsnso(0)+dzsnso(-1))) then - isnow = -1 - dzsnso(-nsnow+1:(isnow-1)) = 0 - dzsnso(isnow) = snowh -dzsnso(isnow+1) - elseif(snowh.le.(dzsnso(0)+dzsnso(-1)+dzsnso(-2))) then isnow = -2 - dzsnso(-nsnow+1:(isnow-2)) = 0 - dzsnso(isnow) = snowh -dzsnso(isnow+2) + dzsnso(-nsnow+1:isnow) = 0 + dzsnso(isnow+1) = snowh -dzsnso(0) + ! scale swe in layers by ratio of depth to pack + do snl_idx=-snow+1,0 + snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) + enddo + snliq(-nsnow+1:isnow) = 0 +! all other cases + elseif(snowh.le.(dzsnso(0)+dzsnso(-1)+dzsnso(-2))) then + isnow = -3 + dzsnso(isnow+1) = snowh -dzsnso(-1) -dzsnso(0) + ! scale swe in layers by ratio of depth to pack + do snl_idx=-snow+1,0 + snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) + enddo + snliq(-nsnow+1:isnow) = 0 endif endif endif From 0bb0f2bcc1abacca0e66c76635379b077cc7c887 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Tue, 6 Jun 2023 11:25:16 -0400 Subject: [PATCH 2/6] fixed nsnow typo --- .../land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 index c1a660d54..8f9173232 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/da_snow/noahmp401_snow_update.F90 @@ -227,7 +227,7 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) dzsnso(-nsnow+1:isnow) = 0 dzsnso(isnow+1) = snowh -dzsnso(0) ! scale swe in layers by ratio of depth to pack - do snl_idx=-snow+1,0 + do snl_idx=-nsnow+1,0 snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) enddo snliq(-nsnow+1:isnow) = 0 @@ -236,7 +236,7 @@ subroutine noahmp401_snow_update(n, t, dsneqv, dsnowh) isnow = -3 dzsnso(isnow+1) = snowh -dzsnso(-1) -dzsnso(0) ! scale swe in layers by ratio of depth to pack - do snl_idx=-snow+1,0 + do snl_idx=-nsnow+1,0 snice(snl_idx) = sneqv*(dzsnso(snl_idx)/snowh) enddo snliq(-nsnow+1:isnow) = 0 From 224fe1867e913a102b16a7711e67086d9c2a1633 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Tue, 6 Jun 2023 16:29:39 -0400 Subject: [PATCH 3/6] added SM precip partitioning scheme --- .../land/noahmp.4.0.1/NoahMP401_lsmMod.F90 | 2 +- .../land/noahmp.4.0.1/NoahMP401_readcrd.F90 | 8 +-- .../land/noahmp.4.0.1/noahmp_driver_401.F90 | 5 +- .../phys/module_sf_noahmp_glacier_401.F90 | 30 ++++++-- .../phys/module_sf_noahmpdrv_401.F90 | 15 ++-- .../phys/module_sf_noahmplsm_401.F90 | 68 ++++++++++++++----- 6 files changed, 90 insertions(+), 38 deletions(-) mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 mode change 100644 => 100755 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 old mode 100644 new mode 100755 index e1b94f570..643292dfe --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 @@ -98,7 +98,7 @@ module NoahMP401_lsmMod ! \item[alb\_opt] ! snow surface albedo (1->BATS; 2->CLASS) ! \item[snf\_opt] -! rainfall \& snowfall (1->Jordan91; 2->BATS; 3->Noah) +! rainfall \& snowfall (1->Jordan91; 2->BATS; 3->Noah; 5->SnowModel+Dai(2008)) ! \item[tbot\_opt] ! lower boundary of soil temperature ! \item[stc\_opt] diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 old mode 100644 new mode 100755 index 9f6ef947d..7565dc933 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 @@ -283,7 +283,7 @@ subroutine NoahMP401_readcrd() NOAHMP401_struc(n)%alb_opt enddo - ! rainfall & snowfall (1->Jordan91; 2->BATS; 3->Noah) + ! rainfall & snowfall (1->Jordan91; 2->BATS; 3->Noah; 5->SnowModel+Dai(2008)) call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 rainfall & snowfall option:", rc = rc) do n=1, LIS_rc%nnest @@ -424,7 +424,7 @@ subroutine NoahMP401_readcrd() NOAHMP401_struc(n)%rformat = "netcdf" enddo ! restart run, read restart file - if (trim(LIS_rc%startcode) == "restart") then +! if (trim(LIS_rc%startcode) == "restart") then Call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 restart file:", rc=rc) do n=1,LIS_rc%nnest @@ -442,7 +442,7 @@ subroutine NoahMP401_readcrd() enddo ! coldstart run, read initial state variables - else +! else ! skin temperature call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 initial surface skin temperature:", rc = rc) @@ -563,7 +563,7 @@ subroutine NoahMP401_readcrd() " not defined") enddo endif - endif +! endif deallocate(nids) 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 old mode 100644 new mode 100755 index 88c747a5f..c4fa3f52e --- a/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 @@ -111,7 +111,8 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & integer, intent(in) :: snf_opt ! precipitation partitioning between snow and rain ! (1-Jordan91; 2->BATS: Snow when SFCTMP < TFRZ+2.2; ! 3->Noah: Snow when SFCTMP < TFRZ; - ! 4->Use WRF precipitation partitioning ) + ! 4->Use WRF precipitation partitioning; + ! 5->Use linear relationship from SnowModel, based on Dai (2008)) integer, intent(in) :: tbot_opt ! lower boundary of soil temperature (1->zero-flux; 2->Noah) integer, intent(in) :: stc_opt ! snow/soil temperature time scheme (1->semi-implicit; 2->fully implicit) @@ -742,7 +743,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & fldfrcin(1,1)=fldfrc - call noahmplsm_401 (LIS_rc%udef, & ! in : LIS undefined value (David Mocko) + call noahmplsm_401 (ttile,LIS_rc%udef, & ! in : LIS undefined value (David Mocko) itimestep,yearlen , julian , coszin , latin , lonin , & ! in : time/space-related dz8w3d(1), dt , zsoil , nsoil , dx , & ! in : model configuration vegetypein, soiltypein, vegfrain, vegmaxin, tbotin , & ! in : Vegetation/Soil characteristics diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 old mode 100644 new mode 100755 index 2a36e31e7..cbd845eea --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 @@ -33,7 +33,7 @@ MODULE NOAHMP_GLACIER_GLOBALS_401 INTEGER :: OPT_ALB != 2 !(suggested 2) ! options for partitioning precipitation into rainfall & snowfall -! 1 -> Jordan (1991); 2 -> BATS: when SFCTMP SFCTMP Jordan (1991); 2 -> BATS: when SFCTMP SFCTMPSnowModel+Dai(2008) INTEGER :: OPT_SNF != 1 !(suggested 1) @@ -280,10 +280,11 @@ SUBROUTINE NOAHMP_GLACIER (& FGEV = EDIR * LATHEA END IF - IF(MAXVAL(SICE) < 0.0001) THEN - WRITE(message,*) "GLACIER HAS MELTED AT:",ILOC,JLOC," ARE YOU SURE THIS SHOULD BE A GLACIER POINT?" - CALL wrf_debug(10,TRIM(message)) - END IF +! MLW remove message +! IF(MAXVAL(SICE) < 0.0001) THEN +! WRITE(message,*) "GLACIER HAS MELTED AT:",ILOC,JLOC," ARE YOU SURE THIS SHOULD BE A GLACIER POINT?" +! CALL wrf_debug(10,TRIM(message)) +! END IF ! water and energy balance check @@ -1208,7 +1209,7 @@ SUBROUTINE SFCDIF1_GLACIER(ITER ,ZLVL ,ZPD ,Z0H ,Z0M , & !in REAL, INTENT(INOUT) :: FH2 !sen heat stability correction, weighted by prior iters ! outputs - REAL, INTENT(INOUT) :: FV !friction velocity (m/s) + REAL, INTENT(OUT) :: FV !friction velocity (m/s) REAL, INTENT(OUT) :: CM !drag coefficient for momentum REAL, INTENT(OUT) :: CH !drag coefficient for heat REAL, INTENT(OUT) :: CH2 !drag coefficient for heat @@ -2065,6 +2066,10 @@ SUBROUTINE WATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,PRCP ,SFCTMP , & !in REAL, DIMENSION( 1:NSOIL) :: SH2O_SAVE !soil liquid water content [m3/m3] INTEGER :: ILEV +! JP added new precip partitioning + REAL, PARAMETER :: TAIR_C_CENTER = 274.26 ! center temperature [k] where FPICE = 0.5 + REAL, PARAMETER :: SLP = -0.30 ! change in FPICE per degree-change + REAL :: BINT !y-intercept, relationship btween air temperature and FPICE #ifdef WRF_HYDRO REAL , INTENT(INOUT) :: sfcheadrt #endif @@ -2112,6 +2117,17 @@ SUBROUTINE WATER_GLACIER (NSNOW ,NSOIL ,IMELT ,DT ,PRCP ,SFCTMP , & !in FPICE = 1.0 ENDIF ENDIF + + ! JP -- Adding precip partitioning option + ! Linear fit to Dai (2008), used in Liston's SnowModel + IF(OPT_SNF == 5) THEN + ! intercept, where 0.5 is the fraction for TAIR_C_CENTER + BINT = 0.5 - SLP * TAIR_C_CENTER + ! solve the equation in form y=mx+b + FPICE = SLP * SFCTMP + BINT + FPICE = MAX(0.0,FPICE) + FPICE = MIN(1.0,FPICE) + ENDIF ! print*, 'fpice: ',fpice ! Hedstrom NR and JW Pomeroy (1998), Hydrol. Processes, 12, 1611-1625 @@ -3028,7 +3044,7 @@ SUBROUTINE ERROR_GLACIER (ILOC ,JLOC ,SWDOWN ,FSA ,FSR ,FIRA , & call wrf_message(trim(message)) WRITE(message,'(i6,1x,i6,1x,5F10.4)')ILOC,JLOC,SAG,FIRA,FSH,FGEV,SSOIL call wrf_message(trim(message)) - call wrf_error_fatal("Energy budget problem in NOAHMP GLACIER") +! call wrf_error_fatal("Energy budget problem in NOAHMP GLACIER") END IF END_WB = SNEQV 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 old mode 100644 new mode 100755 index 5a2b25489..cea3dbad8 --- 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 @@ -11,7 +11,7 @@ MODULE module_sf_noahmpdrv_401 ! CONTAINS ! The subroutine name has been modifed for LIS implemenation Oct 22 2018 - SUBROUTINE noahmplsm_401(LIS_undef_value, & ! IN : LIS undefined value + SUBROUTINE noahmplsm_401(tid,LIS_undef_value, & ! IN : LIS undefined value ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN : Time/Space-related DZ8W, DT, DZS, NSOIL, DX, & ! IN : Model configuration IVGTYP, ISLTYP, VEGFRA, VEGMAX, TMN, & ! IN : Vegetation/Soil characteristics @@ -72,7 +72,8 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, ! IN only -! Added LIS undefined value as an input - David Mocko + ! Added LIS undefined value as an input - David Mocko + integer :: tid REAL, INTENT(IN ) :: LIS_undef_value INTEGER, INTENT(IN ) :: ITIMESTEP ! timestep number INTEGER, INTENT(IN ) :: YR ! 4-digit year @@ -102,7 +103,7 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, INTEGER, INTENT(IN ) :: IOPT_INF ! frozen soil permeability (1-> NY06; 2->Koren99) INTEGER, INTENT(IN ) :: IOPT_RAD ! radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg) INTEGER, INTENT(IN ) :: IOPT_ALB ! snow surface albedo (1->BATS; 2->CLASS) - INTEGER, INTENT(IN ) :: IOPT_SNF ! rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah) + INTEGER, INTENT(IN ) :: IOPT_SNF ! rainfall & snowfall (1-Jordan91; 2->BATS;3->Noah; 5->SnowModel+Dai(2008)) INTEGER, INTENT(IN ) :: IOPT_TBOT ! lower boundary of soil temperature (1->zero-flux; 2->Noah) INTEGER, INTENT(IN ) :: IOPT_STC ! snow/soil temperature time scheme INTEGER, INTENT(IN ) :: IOPT_GLA ! glacier option (1->phase change; 2->simple) @@ -878,7 +879,7 @@ SUBROUTINE noahmplsm_401(LIS_undef_value, ICE=0 ! Neither sea ice or land ice. CALL NOAHMP_SFLX (parameters, & - I , J , LAT , YEARLEN , JULIAN , COSZ , & ! IN : Time/Space-related + tid , J , LAT , YEARLEN , JULIAN , COSZ , & ! IN : Time/Space-related DT , DX , DZ8W1D , NSOIL , ZSOIL , NSNOW , & ! IN : Model configuration FVEG , FVGMAX , VEGTYP , ICE , IST , CROPTYPE, & ! IN : Vegetation/Soil characteristics SMCEQ , & ! IN : Vegetation/Soil characteristics @@ -1602,9 +1603,9 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT ,i,j,snow(i,j),snowh(i,j) CALL wrf_message(err_message) ENDIF - IF ( SNOW( i,j ) > 2000. ) THEN - SNOWH(I,J) = SNOWH(I,J) * 2000. / SNOW(I,J) ! SNOW in mm and SNOWH in m - SNOW (I,J) = 2000. ! cap SNOW at 2000, maintain density + IF ( SNOW( i,j ) > 10000. ) THEN + SNOWH(I,J) = SNOWH(I,J) * 10000. / SNOW(I,J) ! SNOW in mm and SNOWH in m + SNOW (I,J) = 10000. ! cap SNOW at 2000, maintain density ENDIF ENDDO ENDDO diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 old mode 100644 new mode 100755 index 001da5cc8..0686d867e --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 @@ -1,6 +1,7 @@ ! PET from Sujay has been added by Shugong 08/31/2021 MODULE MODULE_SF_NOAHMPLSM_401 + use LIS_coreMod use module_sf_gecros_401, only : gecros IMPLICIT NONE @@ -688,7 +689,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & DZSNSO(IZ) = - ZSNSO(IZ) ELSE DZSNSO(IZ) = ZSNSO(IZ-1) - ZSNSO(IZ) - END IF + END IF END DO ! root-zone temperature @@ -741,6 +742,14 @@ SUBROUTINE NOAHMP_SFLX (parameters, & PAHV ,PAHG ,PAHB ,QRAIN ,QSNOW ,SNOWHIN, & !out FWET ,CMC ) !out + ! KRA ADDED HERE THIS CHECK FOR VALUES BEFORE CALL TO ENERGY ... + IF(SNOWH <= 1.E-6 .OR. SNEQV <= 1.E-3) THEN + SNOWH = 0.0 + SNEQV = 0.0 + END IF + ! KRA + + ! compute energy budget (momentum & energy fluxes and phase changes) CALL ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in @@ -915,6 +924,10 @@ SUBROUTINE ATM (parameters,SFCPRS ,SFCTMP ,Q2 , REAL :: PRCP_FROZEN !total frozen precipitation [mm/s] ! MB/AN : v3.7 REAL, PARAMETER :: RHO_GRPL = 500.0 ! graupel bulk density [kg/m3] ! MB/AN : v3.7 REAL, PARAMETER :: RHO_HAIL = 917.0 ! hail bulk density [kg/m3] ! MB/AN : v3.7 +! JP added for new precip partitioning + REAL, PARAMETER :: TAIR_C_CENTER = 274.26 ! center temperature where FPICE = 0.5 + REAL, PARAMETER :: SLP = -0.30 ! change in FPICE per degree-change + REAL :: BINT !y-intercept, relationship btween air temperature and FPICE ! -------------------------------------------------------------------------------------------------- !jref: seems like PAIR should be P1000mb?? @@ -987,6 +1000,17 @@ SUBROUTINE ATM (parameters,SFCPRS ,SFCTMP ,Q2 , ENDIF ENDIF + ! JP -- Adding precip partitioning option + ! Linear fit to Dai (2008), used in Liston's SnowModel + IF(OPT_SNF == 5) THEN + ! intercept, where 0.5 is the fraction for TAIR_C_CENTER + BINT = 0.5 - SLP * TAIR_C_CENTER + ! solve the equation in form y=mx+b + FPICE = SLP * SFCTMP + BINT + FPICE = MAX(0.0,FPICE) + FPICE = MIN(1.0,FPICE) + ENDIF + ! Hedstrom NR and JW Pomeroy (1998), Hydrol. Processes, 12, 1611-1625 ! fresh snow density @@ -1100,13 +1124,12 @@ SUBROUTINE PHENOLOGY (parameters,VEGTYP ,croptype, SNOWH , TV , LAT , YEA FB = DB / MAX(1.E-06,parameters%HVT-parameters%HVB) IF(parameters%HVT> 0. .AND. parameters%HVT <= 1.0) THEN !MB: change to 1.0 and 0.2 to reflect - SNOWHC = parameters%HVT*EXP(-SNOWH/0.2) ! changes to HVT in MPTABLE - IF(SNOWHC>1.E-06) THEN !Wanshu: avoid very small SNOWHC induced floating invalid - FB = MIN(SNOWH,SNOWHC)/SNOWHC - ELSE - !print *,"small snowh in phenology=",snowh - FB = 1 - END IF + SNOWHC = parameters%HVT*EXP(-SNOWH/0.2) ! changes to HVT in MPTABLE + if(snowh.lt.5) then + FB = MIN(SNOWH,SNOWHC)/SNOWHC + else + FB = 1.0 + endif ENDIF ELAI = LAI*(1.-FB) @@ -1826,15 +1849,11 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in ! ground snow cover fraction [Niu and Yang, 2007, JGR] FSNO = 0. - IF(SNOWH.GT.0.) THEN +! IF(SNOWH.GT.0.) THEN + IF(SNOWH.GT.(0.001)) THEN ! Update by KRA when have very small snowdepth values BDSNO = SNEQV / SNOWH FMELT = (BDSNO/100.)**parameters%MFSNO - if (FMELT<0.000001) then !Bailing Li, added this for GRACE DA to catch smaller values - FSNO = 1 - !print *,"small FMELT due to snowh,fmelt,bdsno,para_mfsno",snowh,fmelt,bdsno,parameters%MFSNO - else FSNO = TANH( SNOWH /(2.5* Z0 * FMELT)) - end if ENDIF ! ground roughness length @@ -2134,7 +2153,7 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in IF(FIRE <=0.) THEN WRITE(6,*) 'emitted longwave <0; skin T may be wrong due to inconsistent' - WRITE(6,*) 'input of SHDFAC with LAI' + WRITE(6,*) 'input of SHDFAC with LAI', LIS_localPet WRITE(6,*) ILOC, JLOC, 'SHDFAC=',FVEG,'VAI=',VAI,'TV=',TV,'TG=',TG WRITE(6,*) 'LWDN=',LWDN,'FIRA=',FIRA,'SNOWH=',SNOWH call wrf_error_fatal("STOP in Noah-MP") @@ -6819,9 +6838,9 @@ SUBROUTINE SNOWWATER (parameters,NSNOW ,NSOIL ,IMELT ,DT ,ZSOIL , & !in !to obtain equilibrium state of snow in glacier region - IF(SNEQV > 2000.) THEN ! 2000 mm -> maximum water depth + IF(SNEQV > 10000.) THEN ! 2000 mm -> maximum water depth BDSNOW = SNICE(0) / DZSNSO(0) - SNOFLOW = (SNEQV - 2000.) + SNOFLOW = (SNEQV - 10000.) SNICE(0) = SNICE(0) - SNOFLOW DZSNSO(0) = DZSNSO(0) - SNOFLOW/BDSNOW SNOFLOW = SNOFLOW / DT @@ -6976,22 +6995,31 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ISNOW_OLD = ISNOW DO J = ISNOW_OLD+1,0 +! print *,'H1' IF (SNICE(J) <= .1) THEN +! print *,'H2' IF(J /= 0) THEN +! print *,'H3' SNLIQ(J+1) = SNLIQ(J+1) + SNLIQ(J) SNICE(J+1) = SNICE(J+1) + SNICE(J) ELSE +! print *,'H4' IF (ISNOW_OLD < -1) THEN ! MB/KM: change to ISNOW +! print *,'H5' SNLIQ(J-1) = SNLIQ(J-1) + SNLIQ(J) SNICE(J-1) = SNICE(J-1) + SNICE(J) ELSE +! print *,'H6' IF(SNICE(J) >= 0.) THEN +! print *,'H7' PONDING1 = SNLIQ(J) ! ISNOW WILL GET SET TO ZERO BELOW; PONDING1 WILL GET SNEQV = SNICE(J) ! ADDED TO PONDING FROM PHASECHANGE PONDING SHOULD BE SNOWH = DZSNSO(J) ! ZERO HERE BECAUSE IT WAS CALCULATED FOR THIN SNOW ELSE ! SNICE OVER-SUBLIMATED EARLIER +! print *,'H8' PONDING1 = SNLIQ(J) + SNICE(J) IF(PONDING1 < 0.) THEN ! IF SNICE AND SNLIQ SUBLIMATES REMOVE FROM SOIL +! print *,'H9' SICE(1) = MAX(0.0,SICE(1)+PONDING1/(DZSNSO(1)*1000.)) PONDING1 = 0.0 END IF @@ -7008,6 +7036,7 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! shift all elements above this down by one. IF (J > ISNOW+1 .AND. ISNOW < -1) THEN +! print *,'H10' DO I = J, ISNOW+2, -1 STC(I) = STC(I-1) SNLIQ(I) = SNLIQ(I-1) @@ -7022,6 +7051,7 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! to conserve water in case of too large surface sublimation IF(SICE(1) < 0.) THEN +! print *,'H11' SH2O(1) = SH2O(1) + SICE(1) SICE(1) = 0. END IF @@ -7034,6 +7064,7 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ZWLIQ = 0. DO J = ISNOW+1,0 +! print *,'H12' SNEQV = SNEQV + SNICE(J) + SNLIQ(J) SNOWH = SNOWH + DZSNSO(J) ZWICE = ZWICE + SNICE(J) @@ -7045,6 +7076,7 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in IF (SNOWH < 0.025 .AND. ISNOW < 0 ) THEN ! MB: change limit ! IF (SNOWH < 0.05 .AND. ISNOW < 0 ) THEN +! print *,'H13' ISNOW = 0 SNEQV = ZWICE PONDING2 = ZWLIQ ! LIMIT OF ISNOW < 0 MEANS INPUT PONDING @@ -7061,12 +7093,14 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! check the snow depth - snow layers combined IF (ISNOW < -1) THEN +! print *,'H14' ISNOW_OLD = ISNOW MSSI = 1 DO I = ISNOW_OLD+1,0 IF (DZSNSO(I) < DZMIN(MSSI)) THEN +! print *,'H15' IF (I == ISNOW+1) THEN NEIBOR = I + 1 From 7d35aad115ffdedd4b9e7b15355165ac190721c5 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Tue, 6 Jun 2023 16:51:37 -0400 Subject: [PATCH 4/6] fixed div. changes from local LIS config --- .../land/noahmp.4.0.1/noahmp_driver_401.F90 | 4 ++-- .../phys/module_sf_noahmp_glacier_401.F90 | 11 +++++------ .../noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 | 15 +++++++-------- .../noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 | 15 --------------- 4 files changed, 14 insertions(+), 31 deletions(-) mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmpdrv_401.F90 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 old mode 100755 new mode 100644 index c4fa3f52e..98c710ba1 --- a/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/noahmp_driver_401.F90 @@ -111,7 +111,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & integer, intent(in) :: snf_opt ! precipitation partitioning between snow and rain ! (1-Jordan91; 2->BATS: Snow when SFCTMP < TFRZ+2.2; ! 3->Noah: Snow when SFCTMP < TFRZ; - ! 4->Use WRF precipitation partitioning; + ! 4->Use WRF precipitation partitioning ! 5->Use linear relationship from SnowModel, based on Dai (2008)) integer, intent(in) :: tbot_opt ! lower boundary of soil temperature (1->zero-flux; 2->Noah) integer, intent(in) :: stc_opt ! snow/soil temperature time scheme (1->semi-implicit; 2->fully implicit) @@ -743,7 +743,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & fldfrcin(1,1)=fldfrc - call noahmplsm_401 (ttile,LIS_rc%udef, & ! in : LIS undefined value (David Mocko) + call noahmplsm_401 (LIS_rc%udef, & ! in : LIS undefined value (David Mocko) itimestep,yearlen , julian , coszin , latin , lonin , & ! in : time/space-related dz8w3d(1), dt , zsoil , nsoil , dx , & ! in : model configuration vegetypein, soiltypein, vegfrain, vegmaxin, tbotin , & ! in : Vegetation/Soil characteristics diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 index cbd845eea..6bd09569d 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 @@ -280,11 +280,10 @@ SUBROUTINE NOAHMP_GLACIER (& FGEV = EDIR * LATHEA END IF -! MLW remove message -! IF(MAXVAL(SICE) < 0.0001) THEN -! WRITE(message,*) "GLACIER HAS MELTED AT:",ILOC,JLOC," ARE YOU SURE THIS SHOULD BE A GLACIER POINT?" -! CALL wrf_debug(10,TRIM(message)) -! END IF + IF(MAXVAL(SICE) < 0.0001) THEN + WRITE(message,*) "GLACIER HAS MELTED AT:",ILOC,JLOC," ARE YOU SURE THIS SHOULD BE A GLACIER POINT?" + CALL wrf_debug(10,TRIM(message)) + END IF ! water and energy balance check @@ -3044,7 +3043,7 @@ SUBROUTINE ERROR_GLACIER (ILOC ,JLOC ,SWDOWN ,FSA ,FSR ,FIRA , & call wrf_message(trim(message)) WRITE(message,'(i6,1x,i6,1x,5F10.4)')ILOC,JLOC,SAG,FIRA,FSH,FGEV,SSOIL call wrf_message(trim(message)) -! call wrf_error_fatal("Energy budget problem in NOAHMP GLACIER") + call wrf_error_fatal("Energy budget problem in NOAHMP GLACIER") END IF END_WB = SNEQV 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 old mode 100755 new mode 100644 index cea3dbad8..5a2b25489 --- 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 @@ -11,7 +11,7 @@ MODULE module_sf_noahmpdrv_401 ! CONTAINS ! The subroutine name has been modifed for LIS implemenation Oct 22 2018 - SUBROUTINE noahmplsm_401(tid,LIS_undef_value, & ! IN : LIS undefined value + SUBROUTINE noahmplsm_401(LIS_undef_value, & ! IN : LIS undefined value ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN : Time/Space-related DZ8W, DT, DZS, NSOIL, DX, & ! IN : Model configuration IVGTYP, ISLTYP, VEGFRA, VEGMAX, TMN, & ! IN : Vegetation/Soil characteristics @@ -72,8 +72,7 @@ SUBROUTINE noahmplsm_401(tid,LIS_undef_value, ! IN only - ! Added LIS undefined value as an input - David Mocko - integer :: tid +! Added LIS undefined value as an input - David Mocko REAL, INTENT(IN ) :: LIS_undef_value INTEGER, INTENT(IN ) :: ITIMESTEP ! timestep number INTEGER, INTENT(IN ) :: YR ! 4-digit year @@ -103,7 +102,7 @@ SUBROUTINE noahmplsm_401(tid,LIS_undef_value, INTEGER, INTENT(IN ) :: IOPT_INF ! frozen soil permeability (1-> NY06; 2->Koren99) INTEGER, INTENT(IN ) :: IOPT_RAD ! radiation transfer (1->gap=F(3D,cosz); 2->gap=0; 3->gap=1-Fveg) INTEGER, INTENT(IN ) :: IOPT_ALB ! snow surface albedo (1->BATS; 2->CLASS) - INTEGER, INTENT(IN ) :: IOPT_SNF ! rainfall & snowfall (1-Jordan91; 2->BATS;3->Noah; 5->SnowModel+Dai(2008)) + INTEGER, INTENT(IN ) :: IOPT_SNF ! rainfall & snowfall (1-Jordan91; 2->BATS; 3->Noah) INTEGER, INTENT(IN ) :: IOPT_TBOT ! lower boundary of soil temperature (1->zero-flux; 2->Noah) INTEGER, INTENT(IN ) :: IOPT_STC ! snow/soil temperature time scheme INTEGER, INTENT(IN ) :: IOPT_GLA ! glacier option (1->phase change; 2->simple) @@ -879,7 +878,7 @@ SUBROUTINE noahmplsm_401(tid,LIS_undef_value, ICE=0 ! Neither sea ice or land ice. CALL NOAHMP_SFLX (parameters, & - tid , J , LAT , YEARLEN , JULIAN , COSZ , & ! IN : Time/Space-related + I , J , LAT , YEARLEN , JULIAN , COSZ , & ! IN : Time/Space-related DT , DX , DZ8W1D , NSOIL , ZSOIL , NSNOW , & ! IN : Model configuration FVEG , FVGMAX , VEGTYP , ICE , IST , CROPTYPE, & ! IN : Vegetation/Soil characteristics SMCEQ , & ! IN : Vegetation/Soil characteristics @@ -1603,9 +1602,9 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT ,i,j,snow(i,j),snowh(i,j) CALL wrf_message(err_message) ENDIF - IF ( SNOW( i,j ) > 10000. ) THEN - SNOWH(I,J) = SNOWH(I,J) * 10000. / SNOW(I,J) ! SNOW in mm and SNOWH in m - SNOW (I,J) = 10000. ! cap SNOW at 2000, maintain density + IF ( SNOW( i,j ) > 2000. ) THEN + SNOWH(I,J) = SNOWH(I,J) * 2000. / SNOW(I,J) ! SNOW in mm and SNOWH in m + SNOW (I,J) = 2000. ! cap SNOW at 2000, maintain density ENDIF ENDDO ENDDO diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 index 0686d867e..bb16e35ef 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 @@ -6995,31 +6995,22 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ISNOW_OLD = ISNOW DO J = ISNOW_OLD+1,0 -! print *,'H1' IF (SNICE(J) <= .1) THEN -! print *,'H2' IF(J /= 0) THEN -! print *,'H3' SNLIQ(J+1) = SNLIQ(J+1) + SNLIQ(J) SNICE(J+1) = SNICE(J+1) + SNICE(J) ELSE -! print *,'H4' IF (ISNOW_OLD < -1) THEN ! MB/KM: change to ISNOW -! print *,'H5' SNLIQ(J-1) = SNLIQ(J-1) + SNLIQ(J) SNICE(J-1) = SNICE(J-1) + SNICE(J) ELSE -! print *,'H6' IF(SNICE(J) >= 0.) THEN -! print *,'H7' PONDING1 = SNLIQ(J) ! ISNOW WILL GET SET TO ZERO BELOW; PONDING1 WILL GET SNEQV = SNICE(J) ! ADDED TO PONDING FROM PHASECHANGE PONDING SHOULD BE SNOWH = DZSNSO(J) ! ZERO HERE BECAUSE IT WAS CALCULATED FOR THIN SNOW ELSE ! SNICE OVER-SUBLIMATED EARLIER -! print *,'H8' PONDING1 = SNLIQ(J) + SNICE(J) IF(PONDING1 < 0.) THEN ! IF SNICE AND SNLIQ SUBLIMATES REMOVE FROM SOIL -! print *,'H9' SICE(1) = MAX(0.0,SICE(1)+PONDING1/(DZSNSO(1)*1000.)) PONDING1 = 0.0 END IF @@ -7036,7 +7027,6 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! shift all elements above this down by one. IF (J > ISNOW+1 .AND. ISNOW < -1) THEN -! print *,'H10' DO I = J, ISNOW+2, -1 STC(I) = STC(I-1) SNLIQ(I) = SNLIQ(I-1) @@ -7051,7 +7041,6 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! to conserve water in case of too large surface sublimation IF(SICE(1) < 0.) THEN -! print *,'H11' SH2O(1) = SH2O(1) + SICE(1) SICE(1) = 0. END IF @@ -7064,7 +7053,6 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ZWLIQ = 0. DO J = ISNOW+1,0 -! print *,'H12' SNEQV = SNEQV + SNICE(J) + SNLIQ(J) SNOWH = SNOWH + DZSNSO(J) ZWICE = ZWICE + SNICE(J) @@ -7076,7 +7064,6 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in IF (SNOWH < 0.025 .AND. ISNOW < 0 ) THEN ! MB: change limit ! IF (SNOWH < 0.05 .AND. ISNOW < 0 ) THEN -! print *,'H13' ISNOW = 0 SNEQV = ZWICE PONDING2 = ZWLIQ ! LIMIT OF ISNOW < 0 MEANS INPUT PONDING @@ -7093,14 +7080,12 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in ! check the snow depth - snow layers combined IF (ISNOW < -1) THEN -! print *,'H14' ISNOW_OLD = ISNOW MSSI = 1 DO I = ISNOW_OLD+1,0 IF (DZSNSO(I) < DZMIN(MSSI)) THEN -! print *,'H15' IF (I == ISNOW+1) THEN NEIBOR = I + 1 From 2e9230f7c2a37961d2772b8ac95754a1b66fbbd9 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Thu, 8 Jun 2023 16:54:36 -0400 Subject: [PATCH 5/6] first-round div NMP+SM development --- .../land/noahmp.4.0.1/NoahMP401_readcrd.F90 | 6 +++--- .../phys/module_sf_noahmp_glacier_401.F90 | 2 +- .../phys/module_sf_noahmplsm_401.F90 | 19 ++++++++++--------- 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 index 7565dc933..59e81cc5f 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 @@ -424,7 +424,7 @@ subroutine NoahMP401_readcrd() NOAHMP401_struc(n)%rformat = "netcdf" enddo ! restart run, read restart file -! if (trim(LIS_rc%startcode) == "restart") then + if (trim(LIS_rc%startcode) == "restart") then Call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 restart file:", rc=rc) do n=1,LIS_rc%nnest @@ -442,7 +442,7 @@ subroutine NoahMP401_readcrd() enddo ! coldstart run, read initial state variables -! else + else ! skin temperature call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 initial surface skin temperature:", rc = rc) @@ -563,7 +563,7 @@ subroutine NoahMP401_readcrd() " not defined") enddo endif -! endif + endif deallocate(nids) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 index 6bd09569d..a4a057ebd 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 @@ -1208,7 +1208,7 @@ SUBROUTINE SFCDIF1_GLACIER(ITER ,ZLVL ,ZPD ,Z0H ,Z0M , & !in REAL, INTENT(INOUT) :: FH2 !sen heat stability correction, weighted by prior iters ! outputs - REAL, INTENT(OUT) :: FV !friction velocity (m/s) + REAL, INTENT(INOUT) :: FV !friction velocity (m/s) REAL, INTENT(OUT) :: CM !drag coefficient for momentum REAL, INTENT(OUT) :: CH !drag coefficient for heat REAL, INTENT(OUT) :: CH2 !drag coefficient for heat diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 index bb16e35ef..8ce5f65f5 100755 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 @@ -1,7 +1,6 @@ ! PET from Sujay has been added by Shugong 08/31/2021 MODULE MODULE_SF_NOAHMPLSM_401 - use LIS_coreMod use module_sf_gecros_401, only : gecros IMPLICIT NONE @@ -689,7 +688,7 @@ SUBROUTINE NOAHMP_SFLX (parameters, & DZSNSO(IZ) = - ZSNSO(IZ) ELSE DZSNSO(IZ) = ZSNSO(IZ-1) - ZSNSO(IZ) - END IF + END IF END DO ! root-zone temperature @@ -742,11 +741,13 @@ SUBROUTINE NOAHMP_SFLX (parameters, & PAHV ,PAHG ,PAHB ,QRAIN ,QSNOW ,SNOWHIN, & !out FWET ,CMC ) !out - ! KRA ADDED HERE THIS CHECK FOR VALUES BEFORE CALL TO ENERGY ... - IF(SNOWH <= 1.E-6 .OR. SNEQV <= 1.E-3) THEN - SNOWH = 0.0 - SNEQV = 0.0 - END IF + IF(OPT_SNF == 5) THEN + ! KRA ADDED HERE THIS CHECK FOR VALUES BEFORE CALL TO ENERGY ... + IF(SNOWH <= 1.E-6 .OR. SNEQV <= 1.E-3) THEN + SNOWH = 0.0 + SNEQV = 0.0 + END IF + ENDIF ! KRA @@ -6838,9 +6839,9 @@ SUBROUTINE SNOWWATER (parameters,NSNOW ,NSOIL ,IMELT ,DT ,ZSOIL , & !in !to obtain equilibrium state of snow in glacier region - IF(SNEQV > 10000.) THEN ! 2000 mm -> maximum water depth + IF(SNEQV > 2000.) THEN ! 2000 mm -> maximum water depth BDSNOW = SNICE(0) / DZSNSO(0) - SNOFLOW = (SNEQV - 10000.) + SNOFLOW = (SNEQV - 2000.) SNICE(0) = SNICE(0) - SNOFLOW DZSNSO(0) = DZSNSO(0) - SNOFLOW/BDSNOW SNOFLOW = SNOFLOW / DT From cbb7f4bc71c94c6c0b6935e00759ee68a61d4774 Mon Sep 17 00:00:00 2001 From: "Justin A. Pflug" Date: Fri, 9 Jun 2023 15:50:10 -0400 Subject: [PATCH 6/6] reverted precip part. and intermediate files --- .../land/noahmp.4.0.1/NoahMP401_lsmMod.F90 | 2 +- .../land/noahmp.4.0.1/NoahMP401_readcrd.F90 | 2 +- .../land/noahmp.4.0.1/noahmp_driver_401.F90 | 3 +- .../phys/module_sf_noahmp_glacier_401.F90 | 17 +------ .../phys/module_sf_noahmplsm_401.F90 | 48 ++++++------------- 5 files changed, 18 insertions(+), 54 deletions(-) mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 mode change 100755 => 100644 lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmplsm_401.F90 diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 old mode 100755 new mode 100644 index 643292dfe..e1b94f570 --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_lsmMod.F90 @@ -98,7 +98,7 @@ module NoahMP401_lsmMod ! \item[alb\_opt] ! snow surface albedo (1->BATS; 2->CLASS) ! \item[snf\_opt] -! rainfall \& snowfall (1->Jordan91; 2->BATS; 3->Noah; 5->SnowModel+Dai(2008)) +! rainfall \& snowfall (1->Jordan91; 2->BATS; 3->Noah) ! \item[tbot\_opt] ! lower boundary of soil temperature ! \item[stc\_opt] diff --git a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 old mode 100755 new mode 100644 index 59e81cc5f..9f6ef947d --- a/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/NoahMP401_readcrd.F90 @@ -283,7 +283,7 @@ subroutine NoahMP401_readcrd() NOAHMP401_struc(n)%alb_opt enddo - ! rainfall & snowfall (1->Jordan91; 2->BATS; 3->Noah; 5->SnowModel+Dai(2008)) + ! rainfall & snowfall (1->Jordan91; 2->BATS; 3->Noah) call ESMF_ConfigFindLabel(LIS_config, & "Noah-MP.4.0.1 rainfall & snowfall option:", rc = rc) do n=1, LIS_rc%nnest 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 98c710ba1..88c747a5f 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 @@ -111,8 +111,7 @@ subroutine noahmp_driver_401(n, ttile, itimestep, & integer, intent(in) :: snf_opt ! precipitation partitioning between snow and rain ! (1-Jordan91; 2->BATS: Snow when SFCTMP < TFRZ+2.2; ! 3->Noah: Snow when SFCTMP < TFRZ; - ! 4->Use WRF precipitation partitioning - ! 5->Use linear relationship from SnowModel, based on Dai (2008)) + ! 4->Use WRF precipitation partitioning ) integer, intent(in) :: tbot_opt ! lower boundary of soil temperature (1->zero-flux; 2->Noah) integer, intent(in) :: stc_opt ! snow/soil temperature time scheme (1->semi-implicit; 2->fully implicit) diff --git a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 old mode 100755 new mode 100644 index a4a057ebd..2a36e31e7 --- a/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 +++ b/lis/surfacemodels/land/noahmp.4.0.1/phys/module_sf_noahmp_glacier_401.F90 @@ -33,7 +33,7 @@ MODULE NOAHMP_GLACIER_GLOBALS_401 INTEGER :: OPT_ALB != 2 !(suggested 2) ! options for partitioning precipitation into rainfall & snowfall -! 1 -> Jordan (1991); 2 -> BATS: when SFCTMP SFCTMPSnowModel+Dai(2008) +! 1 -> Jordan (1991); 2 -> BATS: when SFCTMP SFCTMP 0. .AND. parameters%HVT <= 1.0) THEN !MB: change to 1.0 and 0.2 to reflect - SNOWHC = parameters%HVT*EXP(-SNOWH/0.2) ! changes to HVT in MPTABLE - if(snowh.lt.5) then - FB = MIN(SNOWH,SNOWHC)/SNOWHC - else - FB = 1.0 - endif + SNOWHC = parameters%HVT*EXP(-SNOWH/0.2) ! changes to HVT in MPTABLE + IF(SNOWHC>1.E-06) THEN !Wanshu: avoid very small SNOWHC induced floating invalid + FB = MIN(SNOWH,SNOWHC)/SNOWHC + ELSE + !print *,"small snowh in phenology=",snowh + FB = 1 + END IF ENDIF ELAI = LAI*(1.-FB) @@ -1850,11 +1826,15 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in ! ground snow cover fraction [Niu and Yang, 2007, JGR] FSNO = 0. -! IF(SNOWH.GT.0.) THEN - IF(SNOWH.GT.(0.001)) THEN ! Update by KRA when have very small snowdepth values + IF(SNOWH.GT.0.) THEN BDSNO = SNEQV / SNOWH FMELT = (BDSNO/100.)**parameters%MFSNO + if (FMELT<0.000001) then !Bailing Li, added this for GRACE DA to catch smaller values + FSNO = 1 + !print *,"small FMELT due to snowh,fmelt,bdsno,para_mfsno",snowh,fmelt,bdsno,parameters%MFSNO + else FSNO = TANH( SNOWH /(2.5* Z0 * FMELT)) + end if ENDIF ! ground roughness length @@ -2154,7 +2134,7 @@ SUBROUTINE ENERGY (parameters,ICE ,VEGTYP ,IST ,NSNOW ,NSOIL , & !in IF(FIRE <=0.) THEN WRITE(6,*) 'emitted longwave <0; skin T may be wrong due to inconsistent' - WRITE(6,*) 'input of SHDFAC with LAI', LIS_localPet + WRITE(6,*) 'input of SHDFAC with LAI' WRITE(6,*) ILOC, JLOC, 'SHDFAC=',FVEG,'VAI=',VAI,'TV=',TV,'TG=',TG WRITE(6,*) 'LWDN=',LWDN,'FIRA=',FIRA,'SNOWH=',SNOWH call wrf_error_fatal("STOP in Noah-MP")