Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

MYNN fix for numerical stability issues with mixing snow #71

Merged
merged 8 commits into from
Jun 12, 2023
Merged

MYNN fix for numerical stability issues with mixing snow #71

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
90f4b65
temporary patch for snow mixing which causes numerical instabilities
joeolson42 May 13, 2023
17fc149
small change to table values of Leaf Area Index to match them better …
joeolson42 May 13, 2023
b766e81
fix Cray compiler error having to do with binary/unary operator prece…
grantfirl May 19, 2023
1591a27
add Anders Jensen to Thompson MP CODEOWNERS list
grantfirl May 19, 2023
da02a14
Merge pull request #3 from grantfirl/mynn_snowfix_gjf
joeolson42 May 19, 2023
5ffa961
remove benign bug from smoke/dust mixing when nchem /= ndvel
joeolson42 May 26, 2023
3a1c48d
changes to the irrigation algorithm in RUC to use real-time vegetatio…
joeolson42 May 26, 2023
db322a1
Merge branch 'ufs/dev' into mynn_snowfix
joeolson42 Jun 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions CODEOWNERS
View file
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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
Open in desktop
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