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

Jump to bottom

Add initialization of water vapor mixing ratio at the surface #8

Closed
wants to merge 3 commits into from
Closed

Add initialization of water vapor mixing ratio at the surface #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
c4e0873
Add initialization of water vapor mixing ratio at the surface to
tanyasmirnova Apr 16, 2021
cbbafc5
Bug fix in the composite for ialb=2 option.
tanyasmirnova Apr 16, 2021
d83c1a1
Add fractional code to ialb=1 option used with the Noah LSM.
tanyasmirnova Apr 19, 2021
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
219 changes: 107 additions & 112 deletions physics/radiation_surface.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -449,118 +449,113 @@ subroutine setalb &

do i = 1, IMAX

!> - Calculate snow cover input directly for land model, no
!! conversion needed.
!-- water albedo
asevd_wat = 0.06
asend_wat = 0.06
asevb_wat = asevd_wat
asenb_wat = asevd_wat

! direct albedo CZA dependence over water
if (fraco(i) > f_zero .and. coszf(i) > 0.0001) then
if (tsknf(i) >= con_t0c) then
asevb_wat = max (asevd_wat, 0.026/(coszf(i)**1.7 + 0.065) &
& + 0.15 * (coszf(i)-0.1) * (coszf(i)-0.5) &
& * (coszf(i)-f_one))
asenb_wat = asevb_wat
endif

fsno0 = sncovr(i) ! snow fraction on land

if (nint(slmsk(i))==0 .and. tsknf(i)>con_tice) fsno0 = f_zero

if (nint(slmsk(i)) == 2) then
if(lsm == lsm_ruc) then
!-- use RUC LSM's snow-cover fraction for ice
fsno0 = sncovr_ice(i) ! snow fraction on ice
else
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
endif

fsno1 = f_one - fsno0 ! snow-free fraction (land or ice), 1-sea
flnd0 = min(f_one, facsf(i)+facwf(i)) ! 1-land, 0-sea/ice
fsea0 = max(f_zero, f_one-flnd0)! ! 1-sea/ice, 0-land
fsno = fsno0 ! snow cover, >0 - land/ice
fsea = fsea0 * fsno1 ! 1-sea/ice, 0-land
flnd = flnd0 * fsno1 ! <=1-land,0-sea/ice

!> - 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
if (icy(i)) then
!-- Computation of ice albedo
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
! diffused
if (tsknf(i) < 271.1) then
asevd_ice = 0.70
asend_ice = 0.65
else
b1 = f_zero
a1 = (tsknf(i) - 271.1)**2
asevd_ice = 0.7 - 4.0*a1
asend_ice = 0.65 - 3.6875*a1
endif
! direct
asevb_ice = asevd_ice
asenb_ice = asend_ice

if (fsno0 > f_zero) then
! Snow on ice
dtgd = max(f_zero, min(5.0, (con_ttp-tisfc(i)) ))
b1 = 0.03 * dtgd
asnvd = (asevd_ice + b1) ! diffused snow albedo
asnnd = (asend_ice + b1)
if (coszf(i) > 0.0001 .and. coszf(i) < 0.5) then ! direct snow albedo
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

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
else
asnvb = snoalb(i)
asnnb = snoalb(i)
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))
! composite ice and snow albedos
asevd_ice = asevd_ice * (1. - fsno0) + asnvd * fsno0
asend_ice = asend_ice * (1. - fsno0) + asnnd * fsno0
asevb_ice = asevb_ice * (1. - fsno0) + asnvb * fsno0
asenb_ice = asenb_ice * (1. - fsno0) + asnnb * fsno0
endif ! snow
else
! icy = false, fill in values
asevd_ice = 0.70
asend_ice = 0.65
asevb_ice = 0.70
asenb_ice = 0.65
endif ! end icy

if (tsknf(i) >= con_t0c) then
!- sea
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
if (fracl(i) > f_zero) then
!> - Calculate snow cover input directly for land model, no
!! conversion needed.

fsno0 = sncovr(i) ! snow fraction on land

fsno1 = f_one - fsno0
flnd0 = min(f_one, facsf(i)+facwf(i))
flnd = flnd0 * fsno1 ! snow-free fraction
fsno = f_one - flnd ! snow-covered fraction

!> - use Fanglin's zenith angle treatment.
if (coszf(i) > 0.0001) then
rfcs = 1.775/(1.0+1.55*coszf(i))
else
!- ice
asevb = asevd
asenb = asend
!- no sun
rfcs = f_one
endif
else
!- no sun
rfcs = f_one
asevb = asevd
asenb = asend
endif

!- zenith dependence is applied only to direct beam albedo
ab1bm = min(0.99, alnsf(i)*rfcs)
ab2bm = min(0.99, alvsf(i)*rfcs)
sfcalb(i,1) = ab1bm *flnd + asenb*fsea + asnnb*fsno
sfcalb(i,2) = alnwf(i)*flnd + asend*fsea + asnnd*fsno
sfcalb(i,3) = ab2bm *flnd + asevb*fsea + asnvb*fsno
sfcalb(i,4) = alvwf(i)*flnd + asevd*fsea + asnvd*fsno
!- zenith dependence is applied only to direct beam albedo
ab1bm = min(0.99, alnsf(i)*rfcs)
ab2bm = min(0.99, alvsf(i)*rfcs)

alndnb = ab1bm *flnd + snoalb(i) * fsno
alndnd = alnwf(i)*flnd + snoalb(i) * fsno
alndvb = ab2bm *flnd + snoalb(i) * fsno
alndvd = alvwf(i)*flnd + snoalb(i) * fsno
else
!-- fill in values of land albedo
alndnb = 0.
alndnd = 0.
alndvb = 0.
alndvd = 0.
endif ! end land

!-- Composite mean surface albedo from land, open water and
!-- ice fractions
sfcalb(i,1) = min(0.99,max(0.01,alndnb))*fracl(i) & ! direct beam NIR
& + asenb_wat*fraco(i) + asenb_ice*fraci(i)
sfcalb(i,2) = min(0.99,max(0.01,alndnd))*fracl(i) & ! diffuse NIR
& + asend_wat*fraco(i) + asend_ice*fraci(i)
sfcalb(i,3) = min(0.99,max(0.01,alndvb))*fracl & ! direct beam visible
& + asevb_wat*fraco(i) + asevb_ice*fraci(i)
sfcalb(i,4) = min(0.99,max(0.01,alndvd))*fracl & ! diffuse visible
& + asevd_wat*fraco(i) + asevd_ice*fraci(i)

enddo ! end_do_i_loop

Expand Down Expand Up @@ -647,14 +642,14 @@ subroutine setalb &

!-- Composite mean surface albedo from land, open water and
!-- ice fractions
sfcalb(i,1) = min(0.99,max(0.01,lsmalbdnir(i)))*fracl(i) &
& + asenb_wat*fraco(i) + asenb_ice*fraci(i)
sfcalb(i,2) = min(0.99,max(0.01,lsmalbinir(i)))*fracl(i) &
sfcalb(i,1) = min(0.99,max(0.01,lsmalbdnir(i)))*fracl(i) & ! direct beam NIR
& + asenb_wat*fraco(i) + asenb_ice*fraci(i)
sfcalb(i,2) = min(0.99,max(0.01,lsmalbinir(i)))*fracl(i) & ! diffuse NIR
& + asend_wat*fraco(i) + asend_ice*fraci(i)
sfcalb(i,3) = min(0.99,max(0.01,lsmalbdvis(i)))*fracl(i) &
& + asevb_wat*fraco(i) + asenb_ice*fraci(i)
sfcalb(i,4) = min(0.99,max(0.01,lsmalbivis(i)))*fracl(i) &
& + asevd_wat*fraco(i) + asend_ice*fraci(i)
sfcalb(i,3) = min(0.99,max(0.01,lsmalbdvis(i)))*fracl(i) & ! direct beam visible
& + asevb_wat*fraco(i) + asevb_ice*fraci(i)
sfcalb(i,4) = min(0.99,max(0.01,lsmalbivis(i)))*fracl(i) & ! diffuse visible
& + asevd_wat*fraco(i) + asevd_ice*fraci(i)

enddo ! end_do_i_loop

Expand Down
28 changes: 24 additions & 4 deletions physics/sfc_drv_ruc.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -31,17 +31,18 @@ module lsm_ruc
!! \htmlinclude lsm_ruc_init.html
!!
subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, &
flag_restart, flag_init, &
flag_restart, flag_init, con_fvirt, con_rd, &
im, lsoil_ruc, lsoil, kice, nlev, & ! in
lsm_ruc, lsm, slmsk, stype, vtype, & ! in
tsfc_lnd, tsfc_wat, & ! in
t1, q1, prsl1, tsfc_lnd, tsfc_ice, tsfc_wat, & ! in
tg3, smc, slc, stc, fice, min_seaice, & ! in
sncovr_lnd, sncovr_ice, snoalb, & ! in
facsf, facwf, alvsf, alvwf, alnsf, alnwf, & ! in
sfcqv_lnd, sfcqv_ice, & ! out
sfalb_lnd_bck, & ! out
semisbase, semis_lnd, semis_ice, & ! out
albdvis_lnd,albdnir_lnd,albivis_lnd,albinir_lnd, & ! out
albdvis_ice,albdnir_ice,albivis_ice,albinir_ice, & ! out
semisbase, semis_lnd, semis_ice, & ! out
zs, sh2o, smfrkeep, tslb, smois, wetness, & ! out
tsice, pores, resid, errmsg, errflg)

Expand All @@ -56,12 +57,18 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, &
integer, intent(in) :: kice
integer, intent(in) :: nlev
integer, intent(in) :: lsm_ruc, lsm
real (kind=kind_phys),intent(in) :: con_fvirt
real (kind=kind_phys),intent(in) :: con_rd


real (kind=kind_phys), dimension(im), intent(in) :: slmsk
real (kind=kind_phys), dimension(im), intent(in) :: stype
real (kind=kind_phys), dimension(im), intent(in) :: vtype
real (kind=kind_phys), dimension(im), intent(in) :: t1
real (kind=kind_phys), dimension(im), intent(in) :: q1
real (kind=kind_phys), dimension(im), intent(in) :: prsl1
real (kind=kind_phys), dimension(im), intent(in) :: tsfc_lnd
real (kind=kind_phys), dimension(im), intent(in) :: tsfc_ice
real (kind=kind_phys), dimension(im), intent(in) :: tsfc_wat
real (kind=kind_phys), dimension(im), intent(in) :: tg3
real (kind=kind_phys), dimension(im), intent(in) :: sncovr_lnd
Expand All @@ -87,7 +94,8 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, &
real (kind=kind_phys), dimension(im), intent(inout) :: semis_ice
real (kind=kind_phys), dimension(im), intent(inout) :: &
albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd, &
albdvis_ice, albdnir_ice, albivis_ice, albinir_ice
albdvis_ice, albdnir_ice, albivis_ice, albinir_ice, &
sfcqv_lnd, sfcqv_ice

! --- out
real (kind=kind_phys), dimension(:), intent(out) :: zs
Expand All @@ -102,6 +110,7 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, &
! --- local
real (kind=kind_phys), dimension(lsoil_ruc) :: dzs
real (kind=kind_phys) :: alb_lnd, alb_ice
real (kind=kind_phys) :: q0, qs1, rho
integer :: ipr, i, k
logical :: debug_print
integer, dimension(im) :: soiltyp, vegtype
Expand Down Expand Up @@ -193,6 +202,17 @@ subroutine lsm_ruc_init (me, master, isot, ivegsrc, nlunit, &
albivis_ice(i) = alb_ice
albinir_ice(i) = alb_ice

if (.not.flag_restart) then
!-- initialize QV mixing ratio at the surface from atm. 1st level
q0 = max(q1(i)/(1.-q1(i)), 1.e-8) ! q1=specific humidity at level 1 (kg/kg)
rho = prsl1(i) / (con_rd*t1(i)*(1.0+con_fvirt*q0))
qs1 = rslf(prsl1(i),tsfc_lnd(i)) !* qs1=sat. mixing ratio at level 1 (kg/kg)
q0 = min(qs1, q0)
sfcqv_lnd(i) = q0
qs1 = rslf(prsl1(i),tsfc_ice(i))
sfcqv_ice(i) = qs1
endif

enddo ! i

call init_soil_depth_3 ( zs , dzs , lsoil_ruc )
Expand Down
Loading