diff --git a/physics/CONV/SAMF/samfdeepcnv.f b/physics/CONV/SAMF/samfdeepcnv.f index 18be662f0..2acf3b55b 100644 --- a/physics/CONV/SAMF/samfdeepcnv.f +++ b/physics/CONV/SAMF/samfdeepcnv.f @@ -220,10 +220,9 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! ! parameters for prognostic sigma closure real(kind=kind_phys) omega_u(im,km),zdqca(im,km),tmfq(im,km), - & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km), - & tentr(im,km) - real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins - parameter(sigmind=0.01,sigmins=0.03,sigminm=0.01) + & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im),qadv(im,km) + real(kind=kind_phys) gravinv,invdelt,sigmind,sigminm,sigmins, + & wc_min, wc_eff logical flag_shallow, flag_mid c physical parameters ! parameter(grav=grav,asolfac=0.958) @@ -347,6 +346,21 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !************************************************************************ ! ! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif + km1 = km - 1 !> - Initialize column-integrated and other single-value-per-column variable arrays. c @@ -1137,7 +1151,6 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do i=1,im if(cnvflg(i) .and. & (k > kbcon(i) .and. k < kmax(i))) then - tentr(i,k)=xlamue(i,k)*fent1(i,k) tem = cxlamet(i) * frh(i,k) * fent2(i,k) xlamue(i,k) = xlamue(i,k)*fent1(i,k) + tem tem1 = cxlamdt(i) * frh(i,k) @@ -1788,11 +1801,11 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,tentr,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1878,7 +1891,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -2907,27 +2920,25 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, the convective adjustment time (dtconv) is set to be proportional to the convective turnover time, which is computed using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - if(hwrf_samfdeep) then - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - else - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = min(dtconv(i),dtmax) - endif - enddo - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif + !grid spacing scaling (disabled for HWRF SAMF deep) + if (.not. hwrf_samfdeep) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif + !bounds + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = min(dtconv(i), dtmax) + endif + enddo ! !> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im diff --git a/physics/CONV/SAMF/samfshalcnv.f b/physics/CONV/SAMF/samfshalcnv.f index 1ec8747f0..43d7eeec2 100644 --- a/physics/CONV/SAMF/samfshalcnv.f +++ b/physics/CONV/SAMF/samfshalcnv.f @@ -171,7 +171,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im),qadv(im,km) real(kind=kind_phys) gravinv,dxcrtas,invdelt,sigmind,sigmins, - & sigminm + & sigminm,wc_min,wc_eff logical flag_shallow,flag_mid c physical parameters ! parameter(g=grav,asolfac=0.89) @@ -205,8 +205,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! parameter(bet1=1.875,cd1=.506,f1=2.0,gam1=.5) parameter(betaw=.03,dxcrtc0=9.e3) parameter(h1=0.33333333) -! progsigma - parameter(dxcrtas=500.e3,sigmind=0.01,sigmins=0.03,sigminm=0.01) c local variables and arrays real(kind=kind_phys) pfld(im,km), to(im,km), qo(im,km), & uo(im,km), vo(im,km), qeso(im,km), @@ -254,7 +252,6 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & real(kind=kind_phys) tf, tcr, tcrf parameter (tf=233.16, tcr=263.16, tcrf=1.0/(tcr-tf)) - c----------------------------------------------------------------------- ! ! Initialize CCPP error handling variables @@ -276,8 +273,10 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progsigma) then dxcrt=10.e3 + dxcrtas=500.e3 else dxcrt=15.e3 + dxcrtas=500.e3 endif c----------------------------------------------------------------------- @@ -295,7 +294,21 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & prsl = prslp * 0.001 del = delp * 0.001 !************************************************************************ -! +! - Initialize parameters related to prognostic closure + if (progsigma) then + if (progomega) then + sigmind = 0.03 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + else + sigmind = 0.01 + sigmins = 0.03 + sigminm = 0.03 + wc_min = 0.2 + endif + endif +! km1 = km - 1 c c initialize arrays @@ -1520,11 +1533,11 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & if (progomega) then call progomega_calc(first_time_step,restart,im,km, & kbcon1,ktcon,omegain,delt,del,zi,cnvflg,omegaout, - & grav,buo,drag,wush,xlamue,bb1,bb2) + & grav,buo,drag,wush,bb1,bb2) do k = 1, km do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then omega_u(i,k)=omegaout(i,k) omega_u(i,k)=MAX(omega_u(i,k),-80.) ! Convert to m/s for use in convective time-scale: @@ -1611,7 +1624,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then dp = 1000. * del(i,k) tem = 0.5 * (omega_u(i,k) + omega_u(i,k-1)) omegac(i) = omegac(i) + tem * dp @@ -1636,7 +1649,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & do k = 2, km1 do i = 1, im if (cnvflg(i)) then - if(k > kbcon1(i) .and. k < ktcon(i)) then + if(k >= kbcon1(i) .and. k < ktcon(i)) then if(omega_u(i,k) .ne. 0.)then zeta(i,k)=eta(i,k)*(omegac(i)/omega_u(i,k)) else @@ -1955,21 +1968,28 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & ! compute convective turn-over time ! !> - Following Bechtold et al. (2008) \cite bechtold_et_al_2008, calculate the convective turnover time using the mean updraft velocity (wc) and the cloud depth. It is also proportional to the grid size (gdx). - do i= 1, im - if(cnvflg(i)) then - tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) - dtconv(i) = tem / wc(i) - if (.not.hwrf_samfshal) then - tfac = 1. + gdx(i) / 75000. - dtconv(i) = tfac * dtconv(i) - endif - dtconv(i) = max(dtconv(i),dtmin) - dtconv(i) = max(dtconv(i),dt2) - dtconv(i) = min(dtconv(i),dtmax) - endif + do i = 1, im + if (cnvflg(i)) then + tem = zi(i,ktcon1(i)) - zi(i,kbcon1(i)) + if (progomega) then + wc_eff = max(wc(i), wc_min) + dtconv(i) = tem / wc_eff + else + dtconv(i) = tem / wc(i) + endif +! - grid spacing scaling (disabled for HWRF shallow option) + if (.not. hwrf_samfshal) then + tfac = 1. + gdx(i) / 75000. + dtconv(i) = tfac * dtconv(i) + endif +! - limits + dtconv(i) = max(dtconv(i), dtmin) + dtconv(i) = max(dtconv(i),dt2) + dtconv(i) = min(dtconv(i), dtmax) + endif enddo -! -!> - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. +! +! > - Calculate advective time scale (tauadv) using a mean cloud layer wind speed. do i= 1, im if(cnvflg(i)) then sumx(i) = 0. diff --git a/physics/CONV/progomega_calc.f90 b/physics/CONV/progomega_calc.f90 index 52d99d801..d22659156 100644 --- a/physics/CONV/progomega_calc.f90 +++ b/physics/CONV/progomega_calc.f90 @@ -20,7 +20,7 @@ module progomega !!\section gen_progomega progomega_calc General Algorithm subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegain,delt,del, & - zi,cnvflg,omegaout,grav,buo,drag,wush,tentr,bb1,bb2) + zi,cnvflg,omegaout,grav,buo,drag,wush,bb1,bb2) use machine, only : kind_phys use funcphys, only : fpvs @@ -30,20 +30,17 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai integer, intent(in) :: kbcon1(im),ktcon(im) real(kind=kind_phys), intent(in) :: delt,grav,bb1,bb2 real(kind=kind_phys), intent(in) :: omegain(im,km), del(im,km),zi(im,km) - real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km),tentr(im,km) + real(kind=kind_phys), intent(in) :: drag(im,km),buo(im,km),wush(im,km) real(kind=kind_phys), intent(inout) :: omegaout(im,km) logical, intent(in) :: cnvflg(im),first_time_step,flag_restart real(kind=kind_phys) :: termA(im,km),termB(im,km),termC(im,km),omega(im,km) real(kind=kind_phys) :: RHS(im,km),Kd(im,km) - real(kind=kind_phys) :: dp,dz,entrn,Kdn,discr,wush_pa,lbb1,lbb2,lbb3 + real(kind=kind_phys) :: dp,dz,discr,wush_pa,lbb1,lbb2,lbb3 integer :: i,k - entrn = 0.8E-4 !0.5E-4 !m^-1 - Kdn = 0.5E-4 !2.9E-4 !m^-1 - lbb1 = 0.5 !1.0 - lbb2 = 3.2 !3.0 - lbb3 = 0.5 !0.5 - + lbb1 = 1.5 + lbb2 = 0.6 + lbb3 = 1.2 !Initialization 2D do k = 1,km @@ -56,6 +53,16 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai enddo enddo + do k = 1,km + do i = 1,im + if(cnvflg(i))then + if(omega(i,k) < 1.0E-5) then + omega(i,k) = 0. + endif + endif + enddo + enddo + if(first_time_step .and. .not. flag_restart)then do k = 1,km do i = 1,im @@ -76,10 +83,7 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai do k = 2, km do i = 1, im if (cnvflg(i)) then - if (k > kbcon1(i) .and. k < ktcon(i)) then - - ! Aerodynamic drag parameter - Kd(i,k) = (tentr(i,k)/entrn)*Kdn + if (k >= kbcon1(i) .and. k < ktcon(i)) then ! Scale by dp/dz to have equation in Pa/s !(dp/dz > 0) @@ -91,8 +95,8 @@ subroutine progomega_calc(first_time_step,flag_restart,im,km,kbcon1,ktcon,omegai !termC - Adds buoyancy forcing !Coefficients for the quadratic equation - termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz)) + (Kd(i,k) * (dp/dz))) - termB(i,k) = -1.0 - delt * lbb3 * wush(i,k) * dp/dz + termA(i,k) = delt * ((lbb1 * drag(i,k) * (dp/dz))) + termB(i,k) = 1.0 - delt * lbb3 * wush(i,k) * dp/dz termC(i,k) = omega(i,k) - delt * lbb2 * buo(i,k) * (dp/dz) & - delt * omega(i,k) * (omega(i,k-1) - omega(i,k)) / dp !Compute the discriminant diff --git a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 index 648f6bf81..24bc6f86c 100644 --- a/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 +++ b/physics/Interstitials/UFS_SCM_NEPTUNE/GFS_surface_generic_post.F90 @@ -117,7 +117,7 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, cpl v1(i) = vgrs_1(i) enddo - if (cplflx .or. cplchm .or. cplwav) then + if (cplflx .or. cplchm .or. cplwav .or. cpl_fire) then do i=1,im u10mi_cpl(i) = u10m(i) v10mi_cpl(i) = v10m(i) diff --git a/physics/Radiation/radiation_clouds.f b/physics/Radiation/radiation_clouds.f index 03ce83b7c..727b896b7 100644 --- a/physics/Radiation/radiation_clouds.f +++ b/physics/Radiation/radiation_clouds.f @@ -518,8 +518,8 @@ subroutine radiation_clouds_prop & & iovr_exp, ! Flag for exponential cloud overlap method & iovr_exprand, ! Flag for exponential-random cloud overlap method & idcor, - & idcor_con, - & idcor_hogan, + & idcor_con, + & idcor_hogan, & idcor_oreopoulos @@ -600,7 +600,7 @@ subroutine radiation_clouds_prop & end do end do - if (imp_physics == imp_physics_mg) then ! + if (imp_physics == imp_physics_mg) then ! ! unified cloud and/or with MG microphysics if (uni_cld .and. ncndl >= 2) then @@ -795,7 +795,7 @@ subroutine radiation_clouds_prop & ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) enddo enddo - + ! Compute cloud decorrelation length if (idcor == idcor_hogan) then call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) @@ -1271,7 +1271,7 @@ subroutine progcld_fer_hires & do k = 1, NLAY-1 do i = 1, IX cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & - & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) end do end do else @@ -1280,7 +1280,7 @@ subroutine progcld_fer_hires & do k = 1, NLAY-1 do i = 1, IX cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & - & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) end do end do endif @@ -1548,7 +1548,7 @@ subroutine progcld_thompson_wsm6 & else rei(i,k)=reice_def endif - endif + endif crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) csp(i,k) = max(0.0, (1.-snow2ice)*clw(i,k,ntsw) * & gfac * delp(i,k)) @@ -1586,7 +1586,7 @@ subroutine progcld_thompson_wsm6 & do k = 1, NLAY-1 do i = 1, IX cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & - & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) + & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0.) end do end do else @@ -1597,7 +1597,7 @@ subroutine progcld_thompson_wsm6 & do i = 1, IX cldtot(i,k) = cld_frac_XuRandall(plyr(i,k), qstl(i,k), & & rhly(i,k), clwf(i,k), xrc3, xr_exp, 0., & - & cond_cfrac_onRH) + & cond_cfrac_onRH) end do end do endif @@ -1898,8 +1898,8 @@ subroutine progcld_thompson & rew(i,k) = 9.5 else !--- Land rew(i,k) = 5.5 - endif - endif + endif + endif if (qi1d(k).gt.clwmin .and. cldfra1d(k).lt.ovcst) then cip(i,k) = qi1d(k) * dz1d(k)*1000. idx_rei = int(t1d(k)-179.) @@ -1907,7 +1907,7 @@ subroutine progcld_thompson & corr = t1d(k) - int(t1d(k)) rei(i,K) = max(5.0, retab(idx_rei)*(1.-corr) + & & retab(idx_rei+1)*corr) - endif + endif enddo enddo @@ -3131,7 +3131,7 @@ END SUBROUTINE adjust_cloudFinal !> This function computes the cloud-fraction following !! Xu-Randall(1996) \cite xu_and_randall_1996 !! - function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha, & + function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha, & ! --- inputs & lambda, factor, cond_cfrac_onRH) implicit none ! Inputs @@ -3155,7 +3155,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha, & ! ! Parameters ! real(kind_phys) :: & -! lambda = 0.50 ! , & +! lambda = 0.50 ! , & ! P = 0.25 clwt = 1.0e-6 * (p_lay*0.001) diff --git a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 index a5f855f11..b0c539f60 100644 --- a/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 +++ b/physics/SFC_Models/Land/Noahmp/noahmpdrv.F90 @@ -187,9 +187,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & integer :: lsoil_incr integer, allocatable :: mask_tile(:) - integer,allocatable :: stc_updated(:), slc_updated(:) + integer,allocatable :: stc_updated(:) logical :: soil_freeze, soil_ice - integer :: soiltype, n_stc, n_slc + integer :: soiltype real(kind=kind_phys) :: slc_new integer :: i, j, ij, l, k, ib @@ -197,10 +197,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & real(kind=kind_phys) :: smp !< for computing supercooled water real(kind=kind_phys) :: hc_incr - - integer :: nother, nsnowupd - integer :: nstcupd, nslcupd, nfrozen, nfrozen_upd - logical :: print_update_stats = .False. + real(kind=kind_phys), parameter :: tfreez_noahmp=273.16 ! tfreez used in NoahMP to determine frozen ground ! --- Initialize CCPP error handling variables errmsg = '' @@ -223,7 +220,7 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & endif if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) then - print*, "adding land iau increments " + print*, "adding land iau increments" endif if (Land_IAU_Control%lsoil .ne. km) then @@ -236,11 +233,9 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & allocate(stc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols allocate(slc_inc_flat(Land_IAU_Control%nx * Land_IAU_Control%ny, km)) !GFS_Control%ncols allocate(stc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny)) - allocate(slc_updated(Land_IAU_Control%nx * Land_IAU_Control%ny)) !copy background stc stc_updated = 0 - slc_updated = 0 ib = 1 do j = 1, Land_IAU_Control%ny do k = 1, km @@ -260,53 +255,37 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & lensfc = Land_IAU_Control%nx * Land_IAU_Control%ny if(Land_IAU_Control%me == Land_IAU_Control%mpi_root) print*,' adjusting first ', lsoil_incr, ' surface layers only, delt ', delt - ! initialize variables for counts statitics to be zeros - nother = 0 ! grid cells not land - nsnowupd = 0 ! grid cells with snow (temperature not yet updated) - nstcupd = 0 ! grid cells that are updated stc - nslcupd = 0 ! grid cells that are updated slc - nfrozen = 0 ! not update as frozen soil - nfrozen_upd = 0 ! not update as frozen soil - -!TODO---if only fv3 increment files are used, this can be read from file + allocate(mask_tile(lensfc)) call calculate_landinc_mask(weasd, vegtype, soiltyp, lensfc, isice_table, mask_tile) - !IAU increments are in units of 1/sec !Land_IAU_Control%dtp - !* only updating soil temp for now + dz(1) = -zsoil(1) + do k = 2, km + dz(k) = -zsoil(k) + zsoil(k-1) + enddo + !IAU increments are in units of 1/sec ij_loop : do ij = 1, lensfc ! mask: 1 - soil, 2 - snow, 0 - land-ice, -1 - not land if (mask_tile(ij) == 1) then soil_freeze=.false. soil_ice=.false. - do k = 1, lsoil_incr ! k = 1, km - if ( stc(ij,k) < con_t0c) soil_freeze=.true. + do k = 1, lsoil_incr + if ( stc(ij,k) < tfreez_noahmp ) soil_freeze=.true. if ( smc(ij,k) - slc(ij,k) > 0.001 ) soil_ice=.true. if (Land_IAU_Control%upd_stc) then - stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt !Land_IAU_Control%dtp - if (k==1) then - stc_updated(ij) = 1 - nstcupd = nstcupd + 1 - endif + if (k==1) stc_updated(ij) = 1 + stc(ij,k) = stc(ij,k) + stc_inc_flat(ij,k)*delt endif - if ( (stc(ij,k) < con_t0c) .and. (.not. soil_freeze) .and. (k==1) ) nfrozen_upd = nfrozen_upd + 1 - - ! do not do updates if this layer or any above is frozen + ! do not do SLC updates if this layer or any above is frozen if ( (.not. soil_freeze ) .and. (.not. soil_ice ) ) then if (Land_IAU_Control%upd_slc) then - if (k==1) then - nslcupd = nslcupd + 1 - slc_updated(ij) = 1 - endif - ! apply zero limit here (higher, model-specific limits are later) - slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0) - smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, 0.0) + ! if soil moisture is <0.1 mm in layer, prevent DA from further reducing it + slc(ij,k) = max(slc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), slc(ij,k))) + smc(ij,k) = max(smc(ij,k) + slc_inc_flat(ij,k)*delt, min(0.0001/dz(k), smc(ij,k))) endif - else - if (k==1) nfrozen = nfrozen+1 endif enddo endif ! if soil/snow point @@ -318,27 +297,25 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & errmsg = 'FATAL ERROR in noahmpdrv_timestep_init: problem in set_soilveg_noahmp' return endif - n_stc = 0 - n_slc = 0 if (Land_IAU_Control%do_stcsmc_adjustment) then if (Land_IAU_Control%upd_stc) then do i=1,lensfc if (stc_updated(i) == 1 ) then ! soil-only location - n_stc = n_stc+1 soiltype = soiltyp(i) do l = 1, lsoil_incr if (abs(stc_inc_flat(i,l)) > Land_IAU_Control%min_T_increment) then !the following if case applies when updated stc > melting point, it handles both !case 1: frz ==> frz, recalculate slc, smc remains !case 2: unfrz ==> frz, recalculate slc, smc remains - if (stc(i,l) .LT. con_t0c )then + if (stc(i,l) .LE. tfreez_noahmp )then !recompute supercool liquid water,smc_anl remain unchanged - smp = con_hfus*(con_t0c-stc(i,l))/(con_g*stc(i,l)) !(m) + smp = con_hfus*(tfreez_noahmp-stc(i,l))/(con_g*stc(i,l)) !(m) slc_new=maxsmc(soiltype)*(smp/satpsi(soiltype))**(-1./bb(soiltype)) slc(i,l) = max( min( slc_new, smc(i,l)), 0.0 ) endif - !case 3: frz ==> unfrz (or unfrz ==> unfrz), melt all soil ice (if any) - if (stc(i,l) .GT. con_t0c )then !do not rely on stc_bck + !case 3: frz ==> unfrz melt all soil ice (if any) + if ( stc(i,l) .GT. tfreez_noahmp .and. & + stc(i,l) - stc_inc_flat(i,l)*delt .LE. tfreez_noahmp ) then slc(i,l)=smc(i,l) endif endif @@ -347,34 +324,10 @@ subroutine noahmpdrv_timestep_init (itime, fhour, delt, km, ncols, & enddo endif - if (Land_IAU_Control%upd_slc) then - dz(1) = -zsoil(1) - do l = 2, km - dz(l) = -zsoil(l) + zsoil(l-1) - enddo - do i=1,lensfc - if (slc_updated(i) == 1 ) then - n_slc = n_slc+1 - ! apply SM bounds (later: add upper SMC limit) - do l = 1, lsoil_incr - if (abs(slc_inc_flat(i, l)) > Land_IAU_Control%min_SLC_increment) then - ! noah-mp minimum is 1 mm per layer (in SMC) - ! no need to maintain frozen amount, would be v. small. - slc(i,l) = max( 0.001/dz(l), slc(i,l) ) - smc(i,l) = max( 0.001/dz(l), smc(i,l) ) - endif - enddo - endif - enddo - endif endif - - deallocate(stc_inc_flat, slc_inc_flat, stc_updated, slc_updated) - deallocate(mask_tile) + deallocate(stc_inc_flat, slc_inc_flat, stc_updated) + deallocate(mask_tile) - !Remove non-warning, non-error log write - !write(*,'(a,i4,a,i8)') 'noahmpdrv_timestep_init rank ', Land_IAU_Control%me, ' # of cells with stc update ', nstcupd - end subroutine noahmpdrv_timestep_init