From d26fec3d58f7b9ab99933ed8c1c0dca2d4939b30 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Thu, 17 Dec 2020 22:30:39 +0000 Subject: [PATCH 1/9] Add RUC LSM option to land perturbations. Also, added snow albedo and background albedo components to the perturbations. --- compns_stochy.F90 | 2 +- lndp_apply_perts.F90 | 124 ++++++++++++++++++++++++++++++++++++------- 2 files changed, 106 insertions(+), 20 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index b61e894..fe9fb21 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -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') 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) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index dabd48e..e74694c 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -15,17 +15,22 @@ module lndp_apply_perts_mod !==================================================================== ! Driver for applying perturbations to sprecified land states or parameters ! Draper, July 2020. -! Note on location: requires access to namelist_soilveg +! Note on location: requires access to namelist_soilveg or namelist_soilveg_ruc - 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, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, & + dtf, n_var_lndp, lndp_var_list, & + lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmcnoah, maxsmc, & + param_update_flag, & + smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, & + snoalb, ierr) + !snoalb, semis, ierr) implicit none ! intent(in) integer, intent(in) :: blksz(:) integer, intent(in) :: n_var_lndp, lsoil, lsm + integer, intent(in) :: lsoil_lsm, 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 @@ -35,41 +40,82 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & 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) :: maxsmcnoah(:) real(kind=kind_dbl_prec), intent(in) :: maxsmc(:) + real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:) ! 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(:,:) + 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 + integer :: this_im, v, soiltyp, k, nsoil logical :: print_flag real(kind=kind_dbl_prec) :: p, min_bound, max_bound, tmp_sic, pert + real(kind=kind_dbl_prec), dimension(max(lsoil,lsoil_lsm)) :: 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), parameter :: minsmc = 0.02 + real(kind=kind_dbl_prec), parameter :: zsbot = 6. ierr = 0 - if (lsm .NE. 1 ) then - write(6,*) 'ERROR: lndp_apply_pert assumes LSM is noah, ', & + if (lsm .NE. 1 .and. lsm .ne. 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 + !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,*) 'maxsmcnoah, maxsmc =', maxsmcnoah, maxsmc + + zslayer(:) = 0. + smc_vertscale(:) = 0. + stc_vertscale(:) = 0. + if (lsm == 1) then + nsoil = lsoil + do k = 1, nsoil + 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 + nsoil = lsoil_lsm + zslayer(nsoil) = zsbot - zs_lsm(nsoil) + do k = 1, nsoil-1 + zslayer(k) = zs_lsm(k+1) - zs_lsm(k) + smc_vertscale(k) = smc_vertscale_ruc(k) + stc_vertscale(k) = stc_vertscale_ruc(k) + enddo + endif + nblks = size(blksz) call set_printing_nb_i(blksz,xlon,xlat,print_i,print_nb) @@ -99,15 +145,17 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & p=5. soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc min_bound = minsmc - max_bound = maxsmc(soiltyp) - - do k=1,lsoil + max_bound = maxsmcnoah(soiltyp) + !write (*,*) 'max_bound=', max_bound + + do k=1,nsoil !store frozen soil moisture tmp_sic= smc(nb,i,k) - slc(nb,i,k) ! 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.) + ! 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) @@ -115,10 +163,11 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & ! assign all of applied pert to the liquid soil moisture slc(nb,i,k) = smc(nb,i,k) - tmp_sic enddo + !write (0,*) 'nb, i, smc(nb,i,:)', nb, i, smc(nb,i,:) case('stc') - do k=1,lsoil + do k=1,nsoil pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep ! (necessary for state vars only) @@ -128,14 +177,51 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil,dtf, n_var_lndp, lndp_var_list, & ! Parameter updates - only if param_update_flag = TRUE !================================================================= case('vgf') ! vegetation fraction - if (param_update_flag) then + !if (param_update_flag) then p =5. min_bound=0. max_bound=1. + ! lndp_prt_list(v) = 0.1 in input.nml pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) - endif + !endif + !write (0,*) 'nb, i, vfrac(nb,i)', nb, i, vfrac(nb,i) + case('alb') ! albedo + !if (param_update_flag) then + p =5. + min_bound=0.0 + max_bound=0.4 + + ! lndp_prt_list(v) = 0.4 (wrf) + pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + !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) then + p =5. + min_bound=0.3 + max_bound=0.85 + + ! lndp_prt_list(v) = 0.4 (wrf) + pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound) + !endif + !case('emi') ! emissivity + ! if (param_update_flag) then + ! p =5. + ! min_bound=0. + ! max_bound=1. + + ! lndp_prt_list(v) = 0.1 (wrf) + ! pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + ! 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)) From a44daed3dd16031af01067200109d8eb73207034 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Fri, 18 Dec 2020 19:56:23 +0000 Subject: [PATCH 2/9] Use GFS_Control max and min soil moisture to bound the perturbed soil moisture. --- lndp_apply_perts.F90 | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index e74694c..b939f41 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -19,7 +19,7 @@ module lndp_apply_perts_mod subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, & dtf, n_var_lndp, lndp_var_list, & - lndp_prt_list, sfc_wts, xlon, xlat, stype, maxsmcnoah, maxsmc, & + lndp_prt_list, sfc_wts, xlon, xlat, stype, smcmax, smcmin, & param_update_flag, & smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, & snoalb, ierr) @@ -40,8 +40,8 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, 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) :: maxsmcnoah(:) - real(kind=kind_dbl_prec), intent(in) :: maxsmc(:) + real(kind=kind_dbl_prec), intent(in) :: smcmax(:) + real(kind=kind_dbl_prec), intent(in) :: smcmin(:) real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:) ! intent(inout) @@ -94,7 +94,7 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, !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,*) 'maxsmcnoah, maxsmc =', maxsmcnoah, maxsmc + !write (0,*) 'smcmin =',smcmin zslayer(:) = 0. smc_vertscale(:) = 0. @@ -144,9 +144,8 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, case('smc') p=5. soiltyp = int( stype(nb,i)+0.5 ) ! also need for maxsmc - min_bound = minsmc - max_bound = maxsmcnoah(soiltyp) - !write (*,*) 'max_bound=', max_bound + min_bound = smcmin(soiltyp) + max_bound = smcmax(soiltyp) do k=1,nsoil !store frozen soil moisture @@ -163,7 +162,6 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, ! assign all of applied pert to the liquid soil moisture slc(nb,i,k) = smc(nb,i,k) - tmp_sic enddo - !write (0,*) 'nb, i, smc(nb,i,:)', nb, i, smc(nb,i,:) case('stc') @@ -186,7 +184,6 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) !endif - !write (0,*) 'nb, i, vfrac(nb,i)', nb, i, vfrac(nb,i) case('alb') ! albedo !if (param_update_flag) then p =5. From 57afa4d42b5cb9c2225b83def22f54bbcb9eba4d Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Mon, 28 Dec 2020 22:23:48 +0000 Subject: [PATCH 3/9] Removed the old comment line. --- compns_stochy.F90 | 2 +- lndp_apply_perts.F90 | 13 ++++++------- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index fe9fb21..dd9f3bf 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -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','alb', 'sal') + 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) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index b939f41..e7499db 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -15,7 +15,6 @@ module lndp_apply_perts_mod !==================================================================== ! Driver for applying perturbations to sprecified land states or parameters ! Draper, July 2020. -! Note on location: requires access to namelist_soilveg or namelist_soilveg_ruc subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, & dtf, n_var_lndp, lndp_var_list, & @@ -211,13 +210,13 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, !endif !case('emi') ! emissivity ! if (param_update_flag) then - ! p =5. - ! min_bound=0. - ! max_bound=1. + !p =5. + !min_bound=0. + !max_bound=1. - ! lndp_prt_list(v) = 0.1 (wrf) - ! pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - ! call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) + ! lndp_prt_list(v) = 0.1 (wrf) + !pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + !call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) ! endif case default print*, & From 545fd513951e2611e120446897b53a599f1676a6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 30 Dec 2020 07:56:02 -0700 Subject: [PATCH 4/9] Add semis to lndp_apply_perts --- lndp_apply_perts.F90 | 43 +++++++++++++++++++++---------------------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index e7499db..055ffdf 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -21,26 +21,25 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, lndp_prt_list, sfc_wts, xlon, xlat, stype, smcmax, smcmin, & param_update_flag, & smc, slc, stc, vfrac, alvsf, alnsf, alvwf, alnwf, facsf, facwf, & - snoalb, ierr) - !snoalb, semis, ierr) + snoalb, semis, ierr) implicit none ! intent(in) - integer, intent(in) :: blksz(:) - integer, intent(in) :: n_var_lndp, lsoil, lsm - integer, intent(in) :: lsoil_lsm, lsm_ruc - character(len=3), intent(in) :: lndp_var_list(:) + integer, intent(in) :: blksz(:) + integer, intent(in) :: n_var_lndp, lsoil, lsm + integer, intent(in) :: lsoil_lsm, 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) :: smcmax(:) - real(kind=kind_dbl_prec), intent(in) :: smcmin(:) + real(kind=kind_dbl_prec), intent(in) :: smcmax(:) + real(kind=kind_dbl_prec), intent(in) :: smcmin(:) real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:) ! intent(inout) @@ -55,7 +54,7 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, 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(:,:) + real(kind=kind_dbl_prec), intent(inout) :: semis(:,:) ! intent(out) integer, intent(out) :: ierr @@ -208,16 +207,16 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound) !endif - !case('emi') ! emissivity - ! if (param_update_flag) then - !p =5. - !min_bound=0. - !max_bound=1. - - ! lndp_prt_list(v) = 0.1 (wrf) - !pert = sfc_wts(nb,i,v)*lndp_prt_list(v) - !call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) - ! endif + case('emi') ! emissivity + !if (param_update_flag) then + p =5. + min_bound=0. + max_bound=1. + + !lndp_prt_list(v) = 0.1 (wrf) + pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + 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)) From ac026f6c04a104387d006995d32858b9ef28893a Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 30 Dec 2020 11:24:53 -0700 Subject: [PATCH 5/9] Cleanup work in lndp_apply_perts.F90: define zs_ruc locally, remove trailing whitespaces, pass only one value lsoil from wrapper --- lndp_apply_perts.F90 | 150 ++++++++++++++++++++----------------------- 1 file changed, 68 insertions(+), 82 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 055ffdf..85d0e91 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -14,35 +14,33 @@ module lndp_apply_perts_mod ! lndp_apply_perts !==================================================================== ! Driver for applying perturbations to sprecified land states or parameters -! Draper, July 2020. +! Draper, July 2020. - subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, & - dtf, 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) + subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & + dtf, 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) + ! intent(in) integer, intent(in) :: blksz(:) - integer, intent(in) :: n_var_lndp, lsoil, lsm - integer, intent(in) :: lsoil_lsm, lsm_ruc + integer, intent(in) :: n_var_lndp, lsoil + 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 + 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) :: smcmax(:) real(kind=kind_dbl_prec), intent(in) :: smcmin(:) - real(kind=kind_dbl_prec), intent(in) :: zs_lsm(:) - ! 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(:,:,:) @@ -56,16 +54,16 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:) real(kind=kind_dbl_prec), intent(inout) :: semis(:,:) - ! intent(out) + ! intent(out) integer, intent(out) :: ierr ! local integer :: nblks, print_i, print_nb, i, nb - integer :: this_im, v, soiltyp, k, nsoil - 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), dimension(max(lsoil,lsoil_lsm)) :: zslayer, smc_vertscale, stc_vertscale + real(kind=kind_dbl_prec), dimension(lsoil) :: zslayer, smc_vertscale, stc_vertscale ! decrease in applied pert with depth !-- Noah lsm @@ -75,17 +73,15 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, !-- 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), parameter :: minsmc = 0.02 - real(kind=kind_dbl_prec), parameter :: zsbot = 6. + real(kind=kind_dbl_prec), dimension(9), parameter :: zs_ruc = (/0.00, 0.05, 0.20, 0.40, 0.60, 1.00, 1.60, 2.20, 3.00/) - ierr = 0 + ierr = 0 - if (lsm .NE. 1 .and. lsm .ne. 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 + 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 !write (0,*) 'Input to lndp_apply_pert' @@ -94,21 +90,15 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, !write (0,*) 'n_var_lndp, lndp_var_list =', n_var_lndp, lndp_var_list !write (0,*) 'smcmin =',smcmin - zslayer(:) = 0. - smc_vertscale(:) = 0. - stc_vertscale(:) = 0. - if (lsm == 1) then - nsoil = lsoil - do k = 1, nsoil + 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 - nsoil = lsoil_lsm - zslayer(nsoil) = zsbot - zs_lsm(nsoil) - do k = 1, nsoil-1 - zslayer(k) = zs_lsm(k+1) - zs_lsm(k) + do k = 1, lsoil + zslayer(k) = zs_ruc(k) smc_vertscale(k) = smc_vertscale_ruc(k) stc_vertscale(k) = stc_vertscale_ruc(k) enddo @@ -128,42 +118,42 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, 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 = smcmin(soiltyp) max_bound = smcmax(soiltyp) - - do k=1,nsoil + + 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 ! 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 + 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 - case('stc') + case('stc') - do k=1,nsoil + do k=1,lsoil pert = sfc_wts(nb,i,v)*stc_vertscale(k)*lndp_prt_list(v) pert = pert*dtf/3600. ! lndp_prt_list input is per hour, convert to per timestep ! (necessary for state vars only) @@ -173,22 +163,20 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, ! Parameter updates - only if param_update_flag = TRUE !================================================================= case('vgf') ! vegetation fraction - !if (param_update_flag) then - p =5. + if (param_update_flag) then + p =5. min_bound=0. max_bound=1. - ! lndp_prt_list(v) = 0.1 in input.nml pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('vfrac',pert,print_flag, vfrac(nb,i), ierr,p,min_bound, max_bound) - !endif + endif case('alb') ! albedo - !if (param_update_flag) then + if (param_update_flag) then p =5. min_bound=0.0 max_bound=0.4 - ! lndp_prt_list(v) = 0.4 (wrf) pert = sfc_wts(nb,i,v)*lndp_prt_list(v) !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) @@ -196,35 +184,33 @@ subroutine lndp_apply_perts(blksz,lsm, lsoil, lsm_ruc, lsoil_lsm, zs_lsm, 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 + endif case('sal') ! snow albedo - !if (param_update_flag) then + if (param_update_flag) then p =5. min_bound=0.3 max_bound=0.85 - ! lndp_prt_list(v) = 0.4 (wrf) pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('snoalb',pert,print_flag, snoalb(nb,i), ierr,p,min_bound, max_bound) - !endif + endif case('emi') ! emissivity - !if (param_update_flag) then + if (param_update_flag) then p =5. min_bound=0. max_bound=1. - !lndp_prt_list(v) = 0.1 (wrf) pert = sfc_wts(nb,i,v)*lndp_prt_list(v) call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) - !endif - case default + 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 @@ -241,14 +227,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 @@ -259,9 +245,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 @@ -280,21 +266,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) @@ -317,7 +303,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 From ac5e2415ae04fcc8196a19d239cbcfd840a397e7 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Wed, 30 Dec 2020 22:07:51 +0000 Subject: [PATCH 6/9] Switch from level depths to layer thicknesses in zs_ruc --- lndp_apply_perts.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 85d0e91..51c2842 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -73,7 +73,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & !-- 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.00, 0.05, 0.20, 0.40, 0.60, 1.00, 1.60, 2.20, 3.00/) + 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 From 4391df4555470d3497c184bd4c64b5be999d3d9d Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Thu, 7 Jan 2021 21:07:37 +0000 Subject: [PATCH 7/9] 1. Added time step kdt and lndp_each_step to parameter list. 2. When lndp_each_step=.true. perform conversion from change/per hour to change per time step 3. Mimicking what is done in HRRR, apply soil moisture perturbations to RUC LSM at the beginning of the run when time step = 2. --- lndp_apply_perts.F90 | 33 +++++++++++++++++++++++++-------- 1 file changed, 25 insertions(+), 8 deletions(-) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 51c2842..9bbd1cd 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -17,7 +17,8 @@ module lndp_apply_perts_mod ! Draper, July 2020. subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & - dtf, n_var_lndp, lndp_var_list, lndp_prt_list, & + 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) @@ -26,7 +27,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! intent(in) integer, intent(in) :: blksz(:) - integer, intent(in) :: n_var_lndp, lsoil + 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(:) @@ -62,7 +64,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & 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 @@ -90,6 +92,14 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & !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 == lsm_noah) then do k = 1, lsoil zslayer(k) = zs_noah(k) @@ -135,6 +145,8 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & min_bound = smcmin(soiltyp) 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) @@ -150,6 +162,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! assign all of applied pert to the liquid soil moisture slc(nb,i,k) = smc(nb,i,k) - tmp_sic enddo + endif case('stc') @@ -163,21 +176,23 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & ! Parameter updates - only if param_update_flag = TRUE !================================================================= case('vgf') ! vegetation fraction - if (param_update_flag) then + 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('alb') ! albedo - if (param_update_flag) then + 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) @@ -186,21 +201,23 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & !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) then + 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) then + if (param_update_flag .or. lndp_each_step) then p =5. - min_bound=0. + 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 From ffc88dad1d659ffdd787e2fa9ae305340a5bcacf Mon Sep 17 00:00:00 2001 From: Tanya Smirnova Date: Tue, 12 Jan 2021 19:17:12 +0000 Subject: [PATCH 8/9] Added roughness length over land to the perturbed variables. --- compns_stochy.F90 | 2 +- lndp_apply_perts.F90 | 13 ++++++++++++- 2 files changed, 13 insertions(+), 2 deletions(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index dd9f3bf..c2c2f15 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -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','alb', 'sal','emi') + case('vgf','smc','stc','alb', 'sal','emi','zol','zoi','zow') 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) diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 9bbd1cd..bd3ef9f 100644 --- a/lndp_apply_perts.F90 +++ b/lndp_apply_perts.F90 @@ -21,7 +21,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & 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) + snoalb, semis, zorll, ierr) implicit none @@ -55,6 +55,7 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & real(kind=kind_dbl_prec), intent(inout) :: facsf(:,:) real(kind=kind_dbl_prec), intent(inout) :: facwf(:,:) real(kind=kind_dbl_prec), intent(inout) :: semis(:,:) + real(kind=kind_dbl_prec), intent(inout) :: zorll(:,:) ! intent(out) integer, intent(out) :: ierr @@ -220,6 +221,16 @@ subroutine lndp_apply_perts(blksz, lsm, lsm_noah, lsm_ruc, lsoil, & pert = pert*factor call apply_pert ('semis',pert,print_flag, semis(nb,i), ierr,p,min_bound, max_bound) endif + case('zol') ! land roughness length + if (param_update_flag .or. lndp_each_step) then + p =5. + min_bound=0. + max_bound=300. + + pert = sfc_wts(nb,i,v)*lndp_prt_list(v) + pert = pert*factor + call apply_pert ('zol',pert,print_flag, zorll(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)) From a4ac1335d00290f93e1460e01dd9b95346bc9490 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Thu, 14 Jan 2021 18:02:28 +0000 Subject: [PATCH 9/9] Removed roughnesses for ice (zoi) and water (zow) from the list of perturbed land variables. --- compns_stochy.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/compns_stochy.F90 b/compns_stochy.F90 index c2c2f15..992c18c 100644 --- a/compns_stochy.F90 +++ b/compns_stochy.F90 @@ -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','alb', 'sal','emi','zol','zoi','zow') + case('vgf','smc','stc','alb', 'sal','emi','zol') 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)