diff --git a/lndp_apply_perts.F90 b/lndp_apply_perts.F90 index 055ffdf6..85d0e910 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