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

Add RUC LSM option to land perturbations, add snow albedo and albedo to perturbed parameters #1

Merged
merged 9 commits into from
Jan 9, 2021
Merged

Add RUC LSM option to land perturbations, add snow albedo and albedo to perturbed parameters #1

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
d26fec3
Add RUC LSM option to land perturbations. Also, added snow albedo and
tanyasmirnova Dec 17, 2020
a44daed
Use GFS_Control max and min soil moisture to bound the perturbed soil
tanyasmirnova Dec 18, 2020
57afa4d
Removed the old comment line.
tanyasmirnova Dec 28, 2020
545fd51
Add semis to lndp_apply_perts
climbfuji Dec 30, 2020
d6ed599
Merge pull request #1 from climbfuji/lndp_apply_perts_dom_for_tanya
tanyasmirnova Dec 30, 2020
ac026f6
Cleanup work in lndp_apply_perts.F90: define zs_ruc locally, remove t…
climbfuji Dec 30, 2020
81d091f
Merge pull request #2 from climbfuji/more_cleanup_for_tanya
tanyasmirnova Dec 30, 2020
ac5e241
Switch from level depths to layer thicknesses in zs_ruc
tanyasmirnova Dec 30, 2020
4391df4
1. Added time step kdt and lndp_each_step to parameter list.
tanyasmirnova Jan 7, 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
2 changes: 1 addition & 1 deletion compns_stochy.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -294,7 +294,7 @@ subroutine compns_stochy (me,sz_nml,input_nml_file,fn_nml,nlunit,deltim,iret)
'land perturbations will be applied to selected paramaters, using newer scheme designed for DA ens spread'
do k =1,n_var_lndp
select case (lndp_var_list(k))
case('vgf','smc','stc')
case('vgf','smc','stc','alb', 'sal','emi')
if (me==0) print*, 'land perturbation will be applied to ', lndp_var_list(k)
case default
print*, 'ERROR: land perturbation requested for new parameter - will need to be coded in lndp_apply_pert', lndp_var_list(k)
Expand Down
212 changes: 148 additions & 64 deletions lndp_apply_perts.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,60 +14,104 @@ module lndp_apply_perts_mod
! lndp_apply_perts
!====================================================================
! Driver for applying perturbations to sprecified land states or parameters
! Draper, July 2020.
! Note on location: requires access to namelist_soilveg
! Draper, July 2020.

subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmc,param_update_flag, &
smc, slc, stc, vfrac, ierr)
subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, &
dtf, kdt, lndp_each_step, &
n_var_lndp, lndp_var_list, lndp_prt_list, &
sfc_wts, xlon, xlat, stype, smcmax, smcmin, param_update_flag, &
smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, &
snoalb, semis, ierr)

implicit none

! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, lsm
character(len=3), intent(in) :: lndp_var_list(:)
! intent(in)
integer, intent(in) :: blksz(:)
integer, intent(in) :: n_var_lndp, lsoil, kdt
logical, intent(in) :: lndp_each_step
integer, intent(in) :: lsm, lsm_noah, lsm_ruc
character(len=3), intent(in) :: lndp_var_list(:)
real(kind=kind_dbl_prec), intent(in) :: lndp_prt_list(:)
real(kind=kind_dbl_prec), intent(in) :: dtf
real(kind=kind_dbl_prec), intent(in) :: sfc_wts(:,:,:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
logical, intent(in) :: param_update_flag
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
logical, intent(in) :: param_update_flag
! true = parameters have been updated, apply perts
real(kind=kind_dbl_prec), intent(in) :: stype(:,:)
real(kind=kind_dbl_prec), intent(in) :: maxsmc(:)
real(kind=kind_dbl_prec), intent(in) :: smcmax(:)
real(kind=kind_dbl_prec), intent(in) :: smcmin(:)

! intent(inout)
! intent(inout)
real(kind=kind_dbl_prec), intent(inout) :: smc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: slc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: stc(:,:,:)
real(kind=kind_dbl_prec), intent(inout) :: vfrac(:,:)

! intent(out)
real(kind=kind_dbl_prec), intent(inout) :: snoalb(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alvsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alnsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alvwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: alnwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: facsf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:)
real(kind=kind_dbl_prec), intent(inout) :: semis(:,:)

! intent(out)
integer, intent(out) :: ierr

! local
integer :: nblks, print_i, print_nb, i, nb
integer :: this_im, v, soiltyp, k
logical :: print_flag
integer :: this_im, v, soiltyp, k
logical :: print_flag

real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert
real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert, factor
real(kind=kind_dbl_prec), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale

! decrease in applied pert with depth
real(kind=kind_dbl_prec), dimension(4), parameter :: smc_vertscale = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: stc_vertscale = (/1.0,0.5,0.25,0.125/)

! model-dependent values, hard-wired in noah code.
!-- Noah lsm
real(kind=kind_dbl_prec), dimension(4), parameter :: smc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: stc_vertscale_noah = (/1.0,0.5,0.25,0.125/)
real(kind=kind_dbl_prec), dimension(4), parameter :: zs_noah = (/0.1, 0.3, 0.6, 1.0/)
real(kind=kind_dbl_prec), parameter :: minsmc = 0.02
!-- RUC lsm
real(kind=kind_dbl_prec), dimension(9), parameter :: smc_vertscale_ruc = (/1.0,0.9,0.8,0.6,0.4,0.2,0.1,0.05,0./)
real(kind=kind_dbl_prec), dimension(9), parameter :: stc_vertscale_ruc = (/1.0,0.9,0.8,0.6,0.4,0.2,0.1,0.05,0./)
real(kind=kind_dbl_prec), dimension(9), parameter :: zs_ruc = (/0.05, 0.15, 0.20, 0.20, 0.40, 0.60, 0.60, 0.80, 1.00/)

ierr = 0

if (lsm/=lsm_noah .and. lsm/=lsm_ruc) then
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is Noah or RUC,', &
' may need to adapt variable names for a different LSM'
ierr=10
return
endif

ierr = 0
!write (0,*) 'Input to lndp_apply_pert'
!write (0,*) 'lsm, lsoil, lsm_ruc, lsoil_lsm =', lsm, lsoil, lsm_ruc, lsoil_lsm
!write (0,*) 'zs_lsm =', zs_lsm
!write (0,*) 'n_var_lndp, lndp_var_list =', n_var_lndp, lndp_var_list
!write (0,*) 'smcmin =',smcmin

! lndp_prt_list input is per hour, factor converts to per timestep
! Do conversion only when variables are perturbed at every time step
if(lndp_each_step) then
factor = dtf/3600.
else
factor = 1.
endif

if (lsm .NE. 1 ) then
write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah, ', &
' may need to adapt variable names for a different LSM'
ierr=10
return
if (lsm == lsm_noah) then
do k = 1, lsoil
zslayer(k) = zs_noah(k)
smc_vertscale(k) = smc_vertscale_noah(k)
stc_vertscale(k) = stc_vertscale_noah(k)
enddo
elseif (lsm == lsm_ruc) then
do k = 1, lsoil
zslayer(k) = zs_ruc(k)
smc_vertscale(k) = smc_vertscale_ruc(k)
stc_vertscale(k) = stc_vertscale_ruc(k)
enddo
endif

nblks = size(blksz)
Expand All @@ -84,39 +128,43 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &

if ( smc(nb,i,1) .EQ. 1.) cycle ! skip non-soil (land-ice and non-land)
! set printing
if ( (i==print_i) .and. (nb==print_nb) ) then
if ( (i==print_i) .and. (nb==print_nb) ) then
print_flag = .true.
else
print_flag=.false.
else
print_flag=.false.
endif

do v = 1,n_var_lndp
do v = 1,n_var_lndp
select case (trim(lndp_var_list(v)))
!=================================================================
! State updates - performed every cycle
!=================================================================
case('smc')
p=5.
case('smc')
p=5.
soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc
min_bound = minsmc
max_bound = maxsmc(soiltyp)
min_bound = smcmin(soiltyp)
DomHeinzeller marked this conversation as resolved.
Show resolved Hide resolved
max_bound = smcmax(soiltyp)

if ((lsm /= lsm_ruc) .or. (lsm == lsm_ruc .and. kdt == 2)) then
! with RUC LSM perturb smc only at time step = 2, as in HRRR
do k=1,lsoil
!store frozen soil moisture
tmp_sic= smc(nb,i,k) - slc(nb,i,k)

! perturb total soil moisture
! perturb total soil moisture
! factor of sldepth*1000 converts from mm to m3/m3
pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zs_noah(k)*1000.)
pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep
! lndp_prt_list(v) = 0.3 in input.nml
pert = sfc_wts(nb,i,v)*smc_vertscale(k)*lndp_prt_list(v)/(zslayer(k)*1000.)
pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep
! (necessary for state vars only)
call apply_pert('smc',pert,print_flag, smc(nb,i,k),ierr,p,min_bound, max_bound)

! assign all of applied pert to the liquid soil moisture
! assign all of applied pert to the liquid soil moisture
slc(nb,i,k) = smc(nb,i,k) - tmp_sic
enddo
endif

case('stc')
case('stc')

do k=1,lsoil
pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v)
Expand All @@ -128,22 +176,58 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, &
! Parameter updates - only if param_update_flag = TRUE
DomHeinzeller marked this conversation as resolved.
Show resolved Hide resolved
!=================================================================
case('vgf') ! vegetation fraction
if (param_update_flag) then
p =5.
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.
max_bound=1.

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound)
endif
case default
case('alb') ! albedo
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.0
max_bound=0.4

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
!call apply_pert ('alvsf',pert,print_flag, alvsf(nb,i), ierr,p,min_bound, max_bound)
call apply_pert ('alnsf',pert,print_flag, alnsf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('alvwf',pert,print_flag, alvwf(nb,i), ierr,p,min_bound, max_bound)
call apply_pert ('alnwf',pert,print_flag, alnwf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('facsf',pert,print_flag, facsf(nb,i), ierr,p,min_bound, max_bound)
!call apply_pert ('facwf',pert,print_flag, facwf(nb,i), ierr,p,min_bound, max_bound)
endif
case('sal') ! snow albedo
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.3
max_bound=0.85

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound)
endif
case('emi') ! emissivity
if (param_update_flag .or. lndp_each_step) then
p =5.
min_bound=0.8
max_bound=1.

pert = sfc_wts(nb,i,v)*lndp_prt_list(v)
pert = pert*factor
call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound)
endif
case default
print*, &
'ERROR: unrecognised lndp_prt_list option in lndp_apply_pert, exiting', trim(lndp_var_list(v))
ierr = 10
return
end select
enddo
enddo
'ERROR: unrecognised lndp_prt_list option in lndp_apply_pert, exiting', trim(lndp_var_list(v))
ierr = 10
return
end select
enddo
enddo
enddo
end subroutine lndp_apply_perts

Expand All @@ -160,14 +244,14 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
character(len=*), intent(in) :: vname ! name of variable being perturbed

real(kind=kind_dbl_prec), optional, intent(in) :: p ! flat-top paramater, 0 = no flat-top
! flat-top function is used for bounded variables
! flat-top function is used for bounded variables
! to reduce the magnitude of perturbations near boundaries.
real(kind=kind_dbl_prec), optional, intent(in) :: vmin, vmax ! min,max bounds of variable being perturbed

! intent (inout)
real(kind=kind_dbl_prec), intent(inout) :: state
! intent out

! intent out
integer :: ierr

!local
Expand All @@ -178,9 +262,9 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
endif

! apply perturbation
if (present(p) ) then
if (present(p) ) then
if ( .not. (present(vmin) .and. present(vmax) )) then
ierr=20
ierr=20
print*, 'error, flat-top function requires min & max to be specified'
endif

Expand All @@ -199,21 +283,21 @@ subroutine apply_pert(vname,pert,print_flag, state,ierr,p,vmin, vmax)
endif

end subroutine apply_pert


!====================================================================
! set_printing_nb_i
! set_printing_nb_i
!====================================================================
! routine to turn on print flag for selected location
!
!
subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb)

implicit none
implicit none

! intent (in)
integer, intent(in) :: blksz(:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlon(:,:)
real(kind=kind_dbl_prec), intent(in) :: xlat(:,:)


! intent (out)
Expand All @@ -236,7 +320,7 @@ subroutine set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb)
print_i=i
print_nb=nb
write(*,*) 'LNDP -print flag is on', xlon(nb,i)*57.29578, xlat(nb,i)*57.29578, nb, i
return
return
endif
enddo
enddo
Expand Down