Skip to content

Commit

Permalink
Merge branch 'noahmp' into 'main' (#15)
Browse files Browse the repository at this point in the history
NOAH MP update - work completed by Kai-Yuan Cheng at GFDL

See merge request fv3team/fv3_gfsphysics!48

Co-authored-by: Lucas Harris <[email protected]>
  • Loading branch information
laurenchilutti and lharris4 authored Sep 8, 2022
1 parent 8c46d4f commit 6ad4e72
Show file tree
Hide file tree
Showing 9 changed files with 3,007 additions and 2,787 deletions.
2 changes: 1 addition & 1 deletion FV3GFS/FV3GFS_io.F90
Original file line number Diff line number Diff line change
Expand Up @@ -799,7 +799,7 @@ subroutine register_sfc_prop_restart_vars(Model, nx, ny, nvar_s2m, action)

!--- register the NOAH-MP 3D fields
if (Model%lsm == Model%lsm_noahmp) then
opt = .false.
opt = .true.
do num = nvar_s3+1,nvar_s3+3
var3_p => sfc_var3sn(:,:,:,num)
dim_names_3d(3) = "zaxis_2"
Expand Down
5 changes: 3 additions & 2 deletions GFS_layer/GFS_physics_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1278,7 +1278,7 @@ subroutine GFS_physics_driver &
! --- inputs:
(im, Model%lsoil, kdt, Statein%pgr, Statein%ugrs, Statein%vgrs, &
Statein%tgrs, Statein%qgrs, soiltyp, vegtype, sigmaf, &
Radtend%semis, gabsbdlw, adjsfcdsw, adjsfcnsw, dtf, &
Radtend%semis, adjsfcdlw, adjsfcdsw, adjsfcnsw, dtf, &
Sfcprop%tg3, cd, cdq, Statein%prsl(:,1), work3, &
Diag%zlvl, dry, wind, slopetyp, &
Sfcprop%shdmin, Sfcprop%shdmax, Sfcprop%snoalb, &
Expand Down Expand Up @@ -1306,9 +1306,10 @@ subroutine GFS_physics_driver &
Sfcprop%fastcpxy, Sfcprop%xlaixy, Sfcprop%xsaixy, &
Sfcprop%taussxy, Sfcprop%smoiseq, Sfcprop%smcwtdxy, &
Sfcprop%deeprechxy, Sfcprop%rechxy, &
Sfcprop%albdvis, Sfcprop%albdnir, Sfcprop%albivis, Sfcprop%albinir,&
! --- outputs:
Sfcprop%sncovr, qss, gflx, drain, evap, hflx, ep1d, runof, &
Diag%cmm, Diag%chh, evbs, evcw, sbsno, snowc, Diag%soilm, &
Diag%cmm, Diag%chh, evbs, evcw, sbsno, snowc, Diag%soilm, &
snohf, Diag%smcwlt2, Diag%smcref2, Diag%wet1, t2mmp, q2mp)

endif
Expand Down
7 changes: 5 additions & 2 deletions GFS_layer/GFS_radiation_driver.F90
Original file line number Diff line number Diff line change
Expand Up @@ -1643,7 +1643,9 @@ subroutine GFS_radiation_driver &
tsfg, tsfa, Sfcprop%hprim, Sfcprop%alvsf, &
Sfcprop%alnsf, Sfcprop%alvwf, Sfcprop%alnwf, &
Sfcprop%facsf, Sfcprop%facwf, Sfcprop%fice, &
Sfcprop%tisfc, IM, &
Sfcprop%tisfc,Sfcprop%albdvis, Sfcprop%albdnir, &
Sfcprop%albivis, Sfcprop%albinir, &
IM, &
sfcalb) ! --- outputs

!> -# Approximate mean surface albedo from vis- and nir- diffuse values.
Expand Down Expand Up @@ -1747,7 +1749,8 @@ subroutine GFS_radiation_driver &

call setemis (Grid%xlon, Grid%xlat, Sfcprop%slmsk, & ! --- inputs
Sfcprop%snowd, Sfcprop%sncovr, Sfcprop%zorl, &
tsfg, tsfa, Sfcprop%hprim, IM, &
tsfg, tsfa, Sfcprop%hprim, &
Sfcprop%emiss, IM, &
Radtend%semis) ! --- outputs

!> - Call module_radlw_main::lwrad(), to compute LW heating rates and
Expand Down
1,342 changes: 674 additions & 668 deletions gsmphys/module_sf_noahmp_glacier.f90

Large diffs are not rendered by default.

3,975 changes: 2,000 additions & 1,975 deletions gsmphys/module_sf_noahmplsm.f90

Large diffs are not rendered by default.

248 changes: 124 additions & 124 deletions gsmphys/noahmp_tables.f90

Large diffs are not rendered by default.

177 changes: 170 additions & 7 deletions gsmphys/radiation_surface.f
Original file line number Diff line number Diff line change
Expand Up @@ -183,6 +183,7 @@ subroutine sfc_init &
!! \n physparam::ialbflg
!! - = 0: using climatology surface albedo scheme for SW
!! - = 1: using MODIS based land surface albedo for SW
!! - = 2: using land surface model albedo for SW

if ( ialbflg == 0 ) then

Expand All @@ -196,6 +197,11 @@ subroutine sfc_init &
print *,' - Using MODIS based land surface albedo for sw'
endif

elseif ( ialbflg == 2 ) then ! use albedo from land model

if ( me == 0 ) then
print *,' - Using Albedo From Land Model'
endif
else
print *,' !! ERROR in Albedo Scheme Setting, IALB=',ialbflg
stop
Expand All @@ -205,6 +211,7 @@ subroutine sfc_init &
!! \n physparam::iemsflg
!! - = 0: fixed SFC emissivity at 1.0
!! - = 1: input SFC emissivity type map from "semis_file"
!! - = 2: using SFC emissivity from land model

iemslw = mod(iemsflg, 10) ! emissivity control
if ( iemslw == 0 ) then ! fixed sfc emis at 1.0
Expand Down Expand Up @@ -256,6 +263,11 @@ subroutine sfc_init &
close(NIRADSF)
endif ! end if_file_exist_block

elseif ( iemslw == 2 ) then ! use emiss from land model

if ( me == 0 ) then
print *,' - Using Surface Emissivity From Land Model'
endif
else
print *,' !! ERROR in Emissivity Scheme Setting, IEMS=',iemsflg
stop
Expand Down Expand Up @@ -308,6 +320,7 @@ end subroutine sfc_init
subroutine setalb &
& ( slmsk,snowf,sncovr,snoalb,zorlf,coszf,tsknf,tairf,hprif, & ! --- inputs:
& alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc, &
& lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir, &
& IMAX, &
& sfcalb & ! --- outputs:
& )
Expand Down Expand Up @@ -377,6 +390,7 @@ subroutine setalb &
real (kind=kind_phys), dimension(:), intent(in) :: &
& slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, &
& alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, &
& lsmalbdvis, lsmalbdnir, lsmalbivis, lsmalbinir, &
& sncovr, snoalb

! --- outputs
Expand All @@ -388,11 +402,11 @@ subroutine setalb &
real (kind=kind_phys) :: asnvb, asnnb, asnvd, asnnd, asevb &
&, asenb, asevd, asend, fsno, fsea, rfcs, rfcw, flnd &
&, asnow, argh, hrgh, fsno0, fsno1, flnd0, fsea0, csnow &
&, a1, a2, b1, b2, b3, ab1bm, ab2bm
&, a1, a2, b1, b2, b3, ab1bm, ab2bm, m, s, alpha, beta, albtmp

real (kind=kind_phys) ffw, dtgd

integer :: i, k
integer :: i, k, kk, iflag

!
!===> ... begin here
Expand Down Expand Up @@ -502,7 +516,7 @@ subroutine setalb &
enddo ! end_do_i_loop

!> -# If use modis based albedo for land area:
else
elseif ( ialbflg == 1 ) then

do i = 1, IMAX

Expand Down Expand Up @@ -612,8 +626,136 @@ subroutine setalb &
enddo ! end_do_i_loop
!> -# use land model output for land area:
elseif ( ialbflg == 2 ) then
do i = 1, IMAX
!> - Calculate snow cover input directly for land model, no
!! conversion needed.
fsno0 = f_zero
if (nint(slmsk(i)) == 2) then
asnow = 0.02*snowf(i)
argh = min(0.50, max(.025, 0.01*zorlf(i)))
hrgh = min(f_one, max(0.20, 1.0577-1.1538e-3*hprif(i) ) )
fsno0 = asnow / (argh + asnow) * hrgh
endif
fsno1 = f_one - fsno0
flnd0 = f_one
fsea0 = max(f_zero, f_one-flnd0)
fsno = fsno0
fsea = fsea0 * fsno1
flnd = flnd0 * fsno1
!> - Calculate diffused sea surface albedo.
if (tsknf(i) >= 271.5) then
asevd = 0.06
asend = 0.06
elseif (tsknf(i) < 271.1) then
asevd = 0.70
asend = 0.65
else
a1 = (tsknf(i) - 271.1)**2
asevd = 0.7 - 4.0*a1
asend = 0.65 - 3.6875*a1
endif
!> - Calculate diffused snow albedo, land area use input max snow
!! albedo.
if (nint(slmsk(i)) == 2) then
ffw = f_one - fice(i)
if (ffw < f_one) then
dtgd = max(f_zero, min(5.0, (con_ttp-tisfc(i)) ))
b1 = 0.03 * dtgd
else
b1 = f_zero
endif
b3 = 0.06 * ffw
asnvd = (0.70 + b1) * fice(i) + b3
asnnd = (0.60 + b1) * fice(i) + b3
asevd = 0.70 * fice(i) + b3
asend = 0.60 * fice(i) + b3
else
asnvd = snoalb(i)
asnnd = snoalb(i)
endif
!> - Calculate direct snow albedo.
if (nint(slmsk(i)) == 2) then
if (coszf(i) < 0.5) then
csnow = 0.5 * (3.0 / (f_one+4.0*coszf(i)) - f_one)
asnvb = min( 0.98, asnvd+(f_one-asnvd)*csnow )
asnnb = min( 0.98, asnnd+(f_one-asnnd)*csnow )
else
asnvb = asnvd
asnnb = asnnd
endif
endif
!> - Calculate direct sea surface albedo, use fanglin's zenith angle
!! treatment.

if (coszf(i) > 0.0001) then

! rfcs = 1.89 - 3.34*coszf(i) + 4.13*coszf(i)*coszf(i) &
! & - 2.02*coszf(i)*coszf(i)*coszf(i)
rfcs = 1.775/(1.0+1.55*coszf(i))

if (tsknf(i) >= con_t0c) then
asevb = max(asevd, 0.026/(coszf(i)**1.7+0.065) &
& + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
& * (coszf(i)-f_one))
asenb = asevb
else
asevb = asevd
asenb = asend
endif
else
rfcs = f_one
asevb = asevd
asenb = asend
endif

sfcalb(i,1) = min(0.99,max(0.01,lsmalbdnir(i)))*flnd &
& + asenb*fsea + asnnb*fsno
sfcalb(i,2) = min(0.99,max(0.01,lsmalbinir(i)))*flnd &
& + asend*fsea + asnnd*fsno
sfcalb(i,3) = min(0.99,max(0.01,lsmalbdvis(i)))*flnd &
& + asevb*fsea + asnvb*fsno
sfcalb(i,4) = min(0.99,max(0.01,lsmalbivis(i)))*flnd &
& + asevd*fsea + asnvd*fsno
enddo ! end_do_i_loop

endif ! end if_ialbflg
!

!! sfc-perts, mgehne ***
!! perturb all 4 kinds of surface albedo, sfcalb(:,1:4)
! if (pertalb>0.0) then
! do i = 1, imax
! do kk=1, 4
! ! compute beta distribution parameters for all 4 albedos
! m = sfcalb(i,kk)
! s = pertalb*m*(1.-m)
! alpha = m*m*(1.-m)/(s*s)-m
! beta = alpha*(1.-m)/m
! ! compute beta distribution value corresponding
! ! to the given percentile albPpert to use as new albedo
! call ppfbet(albPpert(i),alpha,beta,iflag,albtmp)
! sfcalb(i,kk) = albtmp
! enddo
! enddo ! end_do_i_loop
! endif
!
!! *** sfc-perts, mgehne

return
!...................................
end subroutine setalb
Expand All @@ -639,7 +781,7 @@ end subroutine setalb
!-----------------------------------
subroutine setemis &
& ( xlon,xlat,slmsk,snowf,sncovr,zorlf,tsknf,tairf,hprif, & ! --- inputs:
& IMAX, &
& lsmemiss, IMAX, &
& sfcemis & ! --- outputs:
& )

Expand All @@ -665,6 +807,8 @@ subroutine setemis &
! tsknf (IMAX) - ground surface temperature in k !
! tairf (IMAX) - lowest model layer air temperature in k !
! hprif (IMAX) - topographic sdv in m !
! lsmemiss(IMAX)- emissivity from lsm !
! !
! IMAX - array horizontal dimension !
! !
! outputs: !
Expand All @@ -688,8 +832,8 @@ subroutine setemis &
integer, intent(in) :: IMAX

real (kind=kind_phys), dimension(:), intent(in) :: &
& xlon,xlat, slmsk, snowf,sncovr, zorlf, tsknf, tairf, hprif
& xlon,xlat, slmsk, snowf,sncovr, zorlf, tsknf, tairf, hprif,&
& lsmemiss
! --- outputs
real (kind=kind_phys), dimension(:), intent(out) :: sfcemis

Expand All @@ -715,7 +859,7 @@ subroutine setemis &
sfcemis(:) = f_one
return

else ! emiss set by sfc type and condition
elseif ( iemslw == 1 ) then ! emiss set by sfc type and condition

dltg = 360.0 / float(IMXEMS)
hdlt = 0.5 * dltg
Expand Down Expand Up @@ -796,6 +940,25 @@ subroutine setemis &

enddo lab_do_IMAX

elseif ( iemslw == 2 ) then ! sfc emiss updated in land model

do i = 1, IMAX

if ( nint(slmsk(i)) == 0 ) then ! sea point

sfcemis(i) = emsref(1)

else if ( nint(slmsk(i)) == 2 ) then ! sea-ice

sfcemis(i) = emsref(7)

else ! land

sfcemis(i) = lsmemiss(i)

endif ! end if_slmsk_block
enddo

endif ! end if_iemslw_block

!chk print *,' In setemis, iemsflg, sfcemis =',iemsflg,sfcemis
Expand Down
Loading

0 comments on commit 6ad4e72

Please sign in to comment.