Skip to content

Commit

Permalink
Merge pull request #71 from joeolson42/mynn_snowfix
Browse files Browse the repository at this point in the history 
MYNN fix for numerical stability issues with mixing snow
  • Loading branch information
@grantfirl
grantfirl authored Jun 12, 2023
2 parents 8cb1643 + db322a1 commit 90c7089
Show file tree
Hide file tree
Showing 5 changed files with 60 additions and 76 deletions.
4 changes: 2 additions & 2 deletions CODEOWNERS
View file
Validating CODEOWNERS rules …
Original file line number Diff line number Diff line change
Expand Up @@ -116,7 +116,7 @@ physics/module_gfdl_cloud_microphys.* @RuiyuSun
physics/module_MP_FER_HIRES.* @ericaligo-NOAA @grantfirl @Qingfu-Liu @dustinswales
physics/module_mp_nssl_2mom.F90 @grantfirl @Qingfu-Liu @dustinswales
physics/module_mp_radar.* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales
physics/module_mp_thompson* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales
physics/module_mp_thompson* @gthompsnWRF @RuiyuSun @AndersJensen-NOAA @grantfirl @Qingfu-Liu @dustinswales
physics/module_nst* @XuLi-NOAA @grantfirl @Qingfu-Liu @dustinswales
physics/module_sf_exchcoef.f90 @grantfirl @Qingfu-Liu @dustinswales
physics/module_SF_JSFC.F90 @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales
Expand All @@ -126,7 +126,7 @@ physics/module_soil_pre.* @tanyasmirnova
physics/moninshoc.* @SMoorthi-emc @grantfirl @Qingfu-Liu @dustinswales
physics/mp_fer_hires.* @ericaligo-NOAA @grantfirl @Qingfu-Liu @dustinswales
physics/mp_nssl.* @grantfirl @Qingfu-Liu @dustinswales
physics/mp_thompson* @gthompsnWRF @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales
physics/mp_thompson* @gthompsnWRF @RuiyuSun @AndersJensen-NOAA @grantfirl @Qingfu-Liu @dustinswales
physics/multi_gases.F90 @RuiyuSun @grantfirl @Qingfu-Liu @dustinswales
physics/myjpbl_wrapper.* @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales
physics/myjsfc_wrapper.* @Qingfu-Liu @grantfirl @Qingfu-Liu @dustinswales
Expand Down
10 changes: 4 additions & 6 deletions physics/module_bl_mynn.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1003,19 +1003,17 @@ SUBROUTINE mynn_bl_driver( &
if ( mix_chem ) then
do ic = 1,ndvel
vd1(ic) = vdep(i,ic) ! dry deposition velocity
chem1(kts,ic) = chem3d(i,kts,ic)
enddo
do k = kts+1,kte
do k = kts,kte
do ic = 1,nchem
chem1(k,ic) = chem3d(i,k,ic)
enddo
enddo
else
do ic = 1,ndvel
vd1(ic) = 0. ! dry deposition velocity
chem1(kts,ic) = 0.
enddo
do k = kts+1,kte
do k = kts,kte
do ic = 1,nchem
chem1(k,ic) = 0.
enddo
Expand Down Expand Up @@ -7521,7 +7519,7 @@ FUNCTION phim(zet)

dummy_0=(1.-am_unst*zet) ! parentesis arg
dummy_1=dummy_0**0.333333 ! y
dummy_11=-0.33333*am_unst*dummy_0**-0.6666667 ! dy/dzet
dummy_11=-0.33333*am_unst*dummy_0**(-0.6666667) ! dy/dzet
dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.) ! f
dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.) ! df/dzet
dummy_3 = 0.57735*(2.*dummy_1+1.) ! g
Expand Down Expand Up @@ -7573,7 +7571,7 @@ FUNCTION phih(zet)

dummy_0=(1.-ah_unst*zet) ! parentesis arg
dummy_1=dummy_0**0.333333 ! y
dummy_11=-0.33333*ah_unst*dummy_0**-0.6666667 ! dy/dzet
dummy_11=-0.33333*ah_unst*dummy_0**(-0.6666667) ! dy/dzet
dummy_2 = 0.33333*(dummy_1**2.+dummy_1+1.) ! f
dummy_22 = 0.3333*dummy_11*(2.*dummy_1+1.) ! df/dzet
dummy_3 = 0.57735*(2.*dummy_1+1.) ! g
Expand Down
96 changes: 41 additions & 55 deletions physics/module_sf_ruclsm.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -415,7 +415,7 @@ SUBROUTINE LSMRUC(xlat,xlon, &
curat, &
INFILTRP
real (kind_phys) :: cq,r61,r273,arp,brp,x,evs,eis
real (kind_phys) :: cropsm
real (kind_phys) :: cropfr, cropsm, newsm, factor

real (kind_phys) :: meltfactor, ac,as, wb,rovcp
INTEGER :: NROOT
Expand Down Expand Up @@ -445,8 +445,8 @@ SUBROUTINE LSMRUC(xlat,xlon, &
NDDZS=2*(nzs-2)

!--
testptlat = 48.7074_kind_phys !39.958 !42.05 !39.0 !74.12 !29.5
testptlon = 289.03_kind_phys !271.622 !286.75 !280.6 !164.0 !283.0
testptlat = 35.55 !48.7074_kind_phys !39.958 !42.05 !39.0 !74.12 !29.5
testptlon = 278.66 !289.03_kind_phys !271.622 !286.75 !280.6 !164.0 !283.0
!--


Expand Down Expand Up @@ -983,63 +983,49 @@ SUBROUTINE LSMRUC(xlat,xlon, &

! Fraction of cropland category in the grid box should not have soil moisture below
! wilting point during the growing season.
! Let's keep soil moisture 20% above wilting point for the fraction of grid box under
! croplands.
! Let's keep soil moisture 5% above wilting point for the crop fraction of grid box.
! This change violates LSM moisture budget, but
! can be considered as a compensation for irrigation not included into LSM.
!tgs - "irrigation" uses fractional landuse, therefore mosaic_lu=1.
! "Irigation" could be applied when landuse fractional information
! is available and mosaic_lu=1.
if(mosaic_lu == 1) then
IF (lufrac(crop) > zero .and. lai(i,j) > 1.1_kind_phys) THEN
! cropland
do k=1,nroot
cropsm=1.1_kind_phys*wilt - qmin
if(soilm1d(k) < cropsm*lufrac(crop)) then
IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.2 .and. &
abs(xlon-testptlon).lt.0.2)then
print * ,'Soil moisture is below wilting in cropland category at time step',ktau
print*,' lat,lon=',xlat,xlon &
,'lufrac(crop),k,soilm1d(k),wilt,cropsm', &
lufrac(crop),k,soilm1d(k),wilt,cropsm
endif
ENDIF
soilm1d(k) = cropsm*lufrac(crop)
IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.2 .and. &
abs(xlon-testptlon).lt.0.2)then
print*,' lat,lon=',xlat,xlon
print * ,'Added soil water to cropland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k)
endif
ENDIF
endif
enddo
! greenness factor: between 0 for min greenness and 1 for max greenness.
factor = max(zero,min(one,(vegfra(i,j)-shdmin(i,j))/max(one,(shdmax(i,j)-shdmin(i,j)))))
if (abs(xlat-testptlat).lt.0.1 .and. &
abs(xlon-testptlon).lt.0.1)then
print *,' lat,lon=',xlat,xlon,' factor=',factor
endif

ELSEIF (ivgtyp(i,j) == natural .and. lai(i,j) > 0.7) THEN
! grassland: assume that 40% of grassland is irrigated cropland
do k=1,nroot
cropsm=1.2_kind_phys*wilt - qmin
if(soilm1d(k) < cropsm*lufrac(natural)*0.4) then
IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.2 .and. &
abs(xlon-testptlon).lt.0.2)then
print * ,'Soil moisture is below wilting in mixed grassland/cropland category at time step',ktau
print*,' lat,lon=',xlat,xlon, &
'lufrac(natural),k,soilm1d(k),wilt', &
lufrac(natural),k,soilm1d(k),wilt
endif
ENDIF
soilm1d(k) = cropsm * lufrac(natural)*0.4_kind_phys
if((ivgtyp(i,j) == natural .or. ivgtyp(i,j) == crop) .and. factor > 0.75) then
! cropland or grassland, apply irrigation during the growing seaspon when fraction
! of greenness is > 0.75.

IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.2 .and. &
abs(xlon-testptlon).lt.0.2)then
print*,' lat,lon=',xlat,xlon
print * ,'Added soil water to grassland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k)
endif
ENDIF
endif
do k=1,nroot
cropsm=1.05_kind_phys*wilt - qmin
cropfr = min(one,lufrac(crop) + 0.4*lufrac(natural)) ! assume that 40% of natural is cropland
newsm = cropsm*cropfr + (1.-cropfr)*soilm1d(k)
if(soilm1d(k) < newsm) then
IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.1 .and. &
abs(xlon-testptlon).lt.0.1)then
print * ,'Soil moisture is below wilting in cropland areas at time step',ktau
print * ,' lat,lon=',xlat,xlon
print * ,' lufrac=',lufrac,'factor=',factor &
,'lai,ivgtyp,lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm,', &
lai(i,j),ivgtyp(i,j),lufrac(crop),k,soilm1d(k),cropfr,wilt,cropsm,newsm
endif
ENDIF
soilm1d(k) = newsm
IF (debug_print ) THEN
if (abs(xlat-testptlat).lt.0.1 .and. &
abs(xlon-testptlon).lt.0.1)then
print*,' lat,lon=',xlat,xlon
print * ,'Added soil water to cropland areas, k,soilm1d(k)',k,soilm1d(k)
endif
ENDIF
endif ! < cropsm
enddo
ENDIF
endif ! crop
endif ! mosaic_lu

!*** DIAGNOSTICS
Expand Down Expand Up @@ -6599,7 +6585,7 @@ SUBROUTINE TRANSF( debug_print, &
TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
ENDIF
!-- uncomment next line for non-linear root distribution
! TRANF(1)=part(1)
TRANF(1)=part(1)

DO K=2,NROOT
totliq=soiliqw(k)+qmin
Expand Down
20 changes: 10 additions & 10 deletions physics/mynnedmf_wrapper.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -418,7 +418,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
FLAG_QI = .true.
FLAG_QNI= .true.
FLAG_QC = .true.
FLAG_QS = .true.
FLAG_QS = .false.
FLAG_QNC= .true.
FLAG_QNWFA= .true.
FLAG_QNIFA= .true.
Expand All @@ -428,7 +428,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
sqv(i,k) = qgrs_water_vapor(i,k)
sqc(i,k) = qgrs_liquid_cloud(i,k)
sqi(i,k) = qgrs_ice(i,k)
sqs(i,k) = qgrs_snow(i,k)
sqs(i,k) = 0. !qgrs_snow(i,k)
qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k)
qni(i,k) = qgrs_cloud_ice_num_conc(i,k)
ozone(i,k) = qgrs_ozone(i,k)
Expand All @@ -441,7 +441,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
FLAG_QI = .true.
FLAG_QNI= .true.
FLAG_QC = .true.
FLAG_QS = .true.
FLAG_QS = .false.
FLAG_QNC= .true.
FLAG_QNWFA= .false.
FLAG_QNIFA= .false.
Expand All @@ -451,7 +451,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
sqv(i,k) = qgrs_water_vapor(i,k)
sqc(i,k) = qgrs_liquid_cloud(i,k)
sqi(i,k) = qgrs_ice(i,k)
sqs(i,k) = qgrs_snow(i,k)
sqs(i,k) = 0. !qgrs_snow(i,k)
qnc(i,k) = qgrs_cloud_droplet_num_conc(i,k)
qni(i,k) = qgrs_cloud_ice_num_conc(i,k)
ozone(i,k) = qgrs_ozone(i,k)
Expand All @@ -464,7 +464,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
FLAG_QI = .true.
FLAG_QNI= .true.
FLAG_QC = .true.
FLAG_QS = .true.
FLAG_QS = .false.
FLAG_QNC= .false.
FLAG_QNWFA= .false.
FLAG_QNIFA= .false.
Expand All @@ -474,7 +474,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
sqv(i,k) = qgrs_water_vapor(i,k)
sqc(i,k) = qgrs_liquid_cloud(i,k)
sqi(i,k) = qgrs_ice(i,k)
sqs(i,k) = qgrs_snow(i,k)
sqs(i,k) = 0. !qgrs_snow(i,k)
qnc(i,k) = 0.
qni(i,k) = qgrs_cloud_ice_num_conc(i,k)
ozone(i,k) = qgrs_ozone(i,k)
Expand Down Expand Up @@ -834,7 +834,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k)
dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k)
dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k))
!dqdt_ozone(i,k) = 0.0
dqdt_water_aer_num_conc(i,k) = RQNWFABLTEN(i,k)
dqdt_ice_aer_num_conc(i,k) = RQNIFABLTEN(i,k)
Expand Down Expand Up @@ -869,7 +869,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
dqdt_cloud_droplet_num_conc(i,k) = RQNCBLTEN(i,k)
dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k)
dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k))
enddo
enddo
if(ldiag3d .and. .not. flag_for_pbl_generic_tend) then
Expand All @@ -887,7 +887,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
dqdt_liquid_cloud(i,k) = RQCBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_ice(i,k) = RQIBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_ice_num_conc(i,k) = RQNIBLTEN(i,k)
dqdt_snow(i,k) = RQSBLTEN(i,k) !/(1.0 + qv(i,k))
dqdt_snow(i,k) = 0.0 !RQSBLTEN(i,k) !/(1.0 + qv(i,k))
!dqdt_ozone(i,k) = 0.0
enddo
enddo
Expand All @@ -896,7 +896,7 @@ SUBROUTINE mynnedmf_wrapper_run( &
call dtend_helper(100+ntcw,RQCBLTEN)
call dtend_helper(100+ntiw,RQIBLTEN)
call dtend_helper(100+ntinc,RQNIBLTEN)
call dtend_helper(100+ntsw,RQSBLTEN)
!call dtend_helper(100+ntsw,RQSBLTEN)
endif
!do k=1,levs
! do i=1,im
Expand Down
6 changes: 3 additions & 3 deletions physics/set_soilveg_ruc.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -190,9 +190,9 @@ subroutine set_soilveg_ruc(me,isot,ivet,nlunit,errmsg,errflg)
ifortbl =(/1, 2, 4, 3, 2, 4, 4, 5, 5, 5, 4, 7, 9, 7, &
& 9, 9, 9, 5, 5, 5, 9, 0, 0, 0, 0, 0, 0, 0, 0, 0/)

laitbl =(/6.40, 6.48, 5.16, 3.31, 5.50, 3.66, 2.60, &
& 3.66, 3.66, 2.90, 5.72, 5.68, 1.00, 4.29, &
& 0.01, 0.75, 0.01, 3.35, 3.35, 3.35, 0.01, &
laitbl =(/2.80, 5.18, 4.16, 4.81, 4.20, 1.16, 0.90, &
& 3.00, 3.00, 1.10, 1.72, 2.58, 1.00, 2.29, &
& 0.01, 0.75, 0.01, 1.00, 1.00, 1.00, 0.01, &
& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, &
& 0.00, 0.00/)

Expand Down

0 comments on commit 90c7089

Please sign in to comment.