diff --git a/components/eam/bld/build-namelist b/components/eam/bld/build-namelist index e77bc9fa4e04..4bd3e91c6380 100755 --- a/components/eam/bld/build-namelist +++ b/components/eam/bld/build-namelist @@ -3848,6 +3848,16 @@ if ($nl->get_value('use_gw_convect') =~ /$TRUE/io) { add_default($nl, 'use_gw_energy_fix' , 'val'=>'.true.'); +add_default($nl, 'use_waterloading' , 'val'=>'.true.'); +add_default($nl, 'use_cpstar' , 'val'=>'.true.'); +add_default($nl, 'use_enthalpy_cpdry' , 'val'=>'.true.'); +add_default($nl, 'use_enthalpy_cl' , 'val'=>'.false.'); +add_default($nl, 'use_enthalpy_cpwv' , 'val'=>'.false.'); +add_default($nl, 'use_enthalpy_theoretical' , 'val'=>'.false.'); +add_default($nl, 'use_global_cpterms_dme' , 'val'=>'.true.'); +add_default($nl, 'dycore_fixer_only' , 'val'=>'.false.'); + + # Some complexity to unpack. # 1. In WACCM, we want tau_0_ubc = .false., because it hasn't been tuned # for that option. diff --git a/components/eam/bld/namelist_files/namelist_defaults_eam.xml b/components/eam/bld/namelist_files/namelist_defaults_eam.xml index e2f9144e6748..aa0958b85ee5 100755 --- a/components/eam/bld/namelist_files/namelist_defaults_eam.xml +++ b/components/eam/bld/namelist_files/namelist_defaults_eam.xml @@ -1678,6 +1678,31 @@ with se_tstep, dt_remap_factor, dt_tracer_factor set to -1 10.0 0.375 .true. + + + .true. + + + + .true. + + + + .true. + + + .false. + + + .false. + + + .true. + + + .false. + + 2.5D0 268.15D0 1.50D0 diff --git a/components/eam/bld/namelist_files/namelist_definition.xml b/components/eam/bld/namelist_files/namelist_definition.xml index 362374c4f716..d489f2501950 100644 --- a/components/eam/bld/namelist_files/namelist_definition.xml +++ b/components/eam/bld/namelist_files/namelist_definition.xml @@ -1028,6 +1028,52 @@ IEFLX fixer options Default: 0 + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + +Default: set by build-namelist. + + + + + + shr_kind_r8 use ppgrid, only: pcols, pver, begchunk, endchunk use spmd_utils, only: masterproc use phys_gmean, only: gmean - use physconst, only: gravit, latvap, latice, cpair, cpairv + use physconst, only: gravit, latvap, latice, cpair, cpairv, cpwv, cpliq, cpice use physics_types, only: physics_state, physics_tend, physics_ptend, physics_ptend_init use constituents, only: cnst_get_ind, pcnst, cnst_name, cnst_get_type_byind, & icldliq, icldice, irain, isnow use time_manager, only: is_first_step use cam_logfile, only: iulog use cam_abortutils, only: endrun - use phys_control, only: ieflx_opt + use phys_control, only: ieflx_opt, use_cpstar, use_waterloading implicit none private @@ -69,6 +70,12 @@ module check_energy public :: check_ieflx_fix ! add ieflx to sensible heat flux public :: energy_helper_eam_def + public :: energy_helper_eam_def_column + public :: fixer_semi_pb + public :: fixer_pb + public :: fixer_pb_simple + public :: dme_adjust + public :: dme_adjust_wl ! Private module data @@ -269,14 +276,42 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type) integer i,k ! column, level indices real(r8) :: wr(state%ncol) ! vertical integral of rain real(r8) :: ws(state%ncol) ! vertical integral of snow + + real(r8) :: qdryini !----------------------------------------------------------------------- lchnk = state%lchnk ncol = state%ncol +! qini (:ncol,:pver) = state%q(:ncol,:pver, 1) +! cldliqini(:ncol,:pver) = state%q(:ncol,:pver,icldliq) +! cldiceini(:ncol,:pver) = state%q(:ncol,:pver,icldice) +! rainini(:ncol,:pver) = state%q(:ncol,:pver,irain) +! snowini(:ncol,:pver) = state%q(:ncol,:pver,isnow) + + do i=1,ncol + do k=1,pver + + qdryini = 1.0 - state%q(i,k, 1) - state%q(i,k,icldliq) & + - state%q(i,k,icldice) - state%q(i,k,irain) & + - state%q(i,k,isnow) + + ! a lot of infrastructure for the energy fixer is called before qini* values are init-ed in tphysbc + !therefore, init cpstar here (this will be called in dp_coupling layer). not the best solution overall. + + !do it before dycore energy fixer + state%cpstar(i,k) = cpair*qdryini + cpwv*state%q(i,k,1) + cpliq*( state%q(i,k,icldliq) + state%q(i,k,irain) ) + & + cpice*( state%q(i,k,icldice) + state%q(i,k,isnow) ) + + enddo + enddo + + + call energy_helper_eam_def(state%u,state%v,state%T,state%q,state%ps,state%pdel,state%phis, & ke,se,wv,wl,wi,wr,ws,te,tw, & - ncol) + ncol, & + state%cpstar) state%te_ini(:ncol) = te(:ncol) state%tw_ini(:ncol) = tw(:ncol) @@ -298,6 +333,7 @@ subroutine check_energy_timestep_init(state, tend, pbuf, col_type) end subroutine check_energy_timestep_init + !=============================================================================== subroutine check_energy_chng(state, tend, name, nstep, ztodt, & @@ -361,7 +397,8 @@ subroutine check_energy_chng(state, tend, name, nstep, ztodt, & call energy_helper_eam_def(state%u,state%v,state%T,state%q,state%ps,state%pdel,state%phis, & ke,se,wv,wl,wi,wr,ws,te,tw, & - ncol) + ncol, & + cpstar=state%cpstar) ! compute expected values and tendencies do i = 1, ncol @@ -459,10 +496,11 @@ subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) integer :: ncol ! number of active columns integer :: lchnk ! chunk index - real(r8) :: te(pcols,begchunk:endchunk,3) + real(r8) :: te(pcols,begchunk:endchunk,13) ! total energy of input/output states (copy) - real(r8) :: te_glob(3) ! global means of total energy + real(r8) :: te_glob(13) ! global means of total energy real(r8), pointer :: teout(:) + integer :: i !----------------------------------------------------------------------- ! Copy total energy out of input and output states @@ -476,14 +514,27 @@ subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) ! output energy call pbuf_get_field(pbuf_get_chunk(pbuf2d,lchnk),teout_idx, teout) - te(:ncol,lchnk,2) = teout(1:ncol) + te(:ncol,lchnk,2) = teout(:ncol) ! surface pressure for heating rate te(:ncol,lchnk,3) = state(lchnk)%pint(:ncol,pver+1) + + te(:ncol,lchnk,4) = state(lchnk)%cptermp(:ncol) + te(:ncol,lchnk,5) = state(lchnk)%cpterme(:ncol) + te(:ncol,lchnk,6) = state(lchnk)%pw(:ncol) + + te(:ncol,lchnk,7) = state(lchnk)%qflx(:ncol) + te(:ncol,lchnk,8) = state(lchnk)%liqflx(:ncol) + te(:ncol,lchnk,9) = state(lchnk)%iceflx(:ncol) + te(:ncol,lchnk,10) = state(lchnk)%dvapor(:ncol) + te(:ncol,lchnk,11) = state(lchnk)%dliquid(:ncol) + te(:ncol,lchnk,12) = state(lchnk)%dice(:ncol) + te(:ncol,lchnk,13) = state(lchnk)%pwvapor(:ncol) + end do ! Compute global means of input and output energies and of ! surface pressure for heating rate (assume uniform ptop) - call gmean(te, te_glob, 3) + call gmean(te, te_glob, 13) if (begchunk .le. endchunk) then teinp_glob = te_glob(1) @@ -497,6 +548,13 @@ subroutine check_energy_gmean(state, pbuf2d, dtime, nstep) if (masterproc) then write(iulog,'(1x,a9,1x,i8,4(1x,e25.17))') "nstep, te", nstep, teinp_glob, teout_glob, heat_glob, psurf_glob + write(iulog,'(1x,a13,1x,i8,2(1x,e25.17))') "nstep, pw, cp", nstep, te_glob(6), (te_glob(4)-te_glob(5))*dtime + write(iulog,'(1x,a15,1x,i8,2(1x,e25.17))') "nstep, cpp, cpe", nstep, te_glob(4)*dtime, te_glob(5)*dtime + + write(iulog,'(1x,a19,1x,i8,2(1x,e25.17))') "nstep, qflx, dqv/dt", nstep, te_glob(7), te_glob(10)/dtime + write(iulog,'(1x,a20,1x,i8,2(1x,e25.17))') "nstep, liqflx, dql/dt", nstep, te_glob(8), te_glob(11)/dtime + write(iulog,'(1x,a20,1x,i8,2(1x,e25.17))') "nstep, iceflx, dqi/dt", nstep, te_glob(9), te_glob(12)/dtime + write(iulog,'(1x,a19,1x,i8,2(1x,e25.17))') "nstep, pw, pwvapor", nstep, te_glob(6), te_glob(13) end if else heat_glob = 0._r8 @@ -1115,7 +1173,10 @@ end subroutine check_prect subroutine energy_helper_eam_def(u,v,T,q,ps,pdel,phis, & ke,se,wv,wl,wi,wr,ws,te,tw, & - ncol,teloc,psterm) + ncol, & + cpstar, & + qini,cldliqini,cldiceini,rainini,snowini,qdryini, & + teloc,psterm) !state vars are of size psetcols,pver, so, not exactly correct real(r8), intent(in) :: u(pcols,pver) @@ -1127,6 +1188,15 @@ subroutine energy_helper_eam_def(u,v,T,q,ps,pdel,phis, & real(r8), intent(in) :: phis(pcols) + real(r8), intent(in ), optional :: qini(pcols,pver) ! initial specific humidity + real(r8), intent(in ), optional :: cldliqini(pcols,pver) ! initial specific humidity + real(r8), intent(in ), optional :: cldiceini(pcols,pver) ! initial specific humidity + real(r8), intent(in ), optional :: rainini(pcols,pver) ! initial specific humidity + real(r8), intent(in ), optional :: snowini(pcols,pver) ! initial specific humidity + real(r8), intent(in ), optional :: qdryini(pcols,pver) ! initial specific humidity + + real(r8), intent(in ), optional :: cpstar(pcols,pver) + real(r8), intent(inout) :: ke(ncol) ! vertical integral of kinetic energy real(r8), intent(inout) :: se(ncol) ! vertical integral of static energy real(r8), intent(inout) :: wv(ncol) ! vertical integral of water (vapor) @@ -1138,6 +1208,7 @@ subroutine energy_helper_eam_def(u,v,T,q,ps,pdel,phis, & real(r8), intent(inout) :: ws(ncol) ! vertical integral of snow ! do not use in this version +!why are these pcol variables? real(r8), intent(inout), optional :: teloc(pcols,pver) real(r8), intent(inout), optional :: psterm(pcols) @@ -1145,19 +1216,58 @@ subroutine energy_helper_eam_def(u,v,T,q,ps,pdel,phis, & integer :: i,k if (icldliq > 1 .and. icldice > 1 .and. irain > 1 .and. isnow > 1) then + + if(present(cpstar))then + +!use cpstar +!print *, "OG we use cpstar in energy_helper_eam_def" + do i = 1, ncol + call energy_helper_eam_def_column(u(i,:),v(i,:),T(i,:),q(i,1:pver,1:pcnst),& + ps(i),pdel(i,:),phis(i), & + ke(i),se(i),wv(i),wl(i),wi(i),wr(i),ws(i),te(i),tw(i), & + cpstar=cpstar(i,:)) + enddo + + elseif(present(qini))then + +!use QINI + +!print *, "OG we use new def of energy with init fields" + do i = 1, ncol + call energy_helper_eam_def_column(u(i,:),v(i,:),T(i,:),q(i,1:pver,1:pcnst),& + ps(i),pdel(i,:),phis(i), & + ke(i),se(i),wv(i),wl(i),wi(i),wr(i),ws(i),te(i),tw(i), & + qini=qini(i,:),cldliqini=cldliqini(i,:),cldiceini=cldiceini(i,:),& + rainini=rainini(i,:),snowini=snowini(i,:),qdryini=qdryini(i,:) ) + enddo + else + +!use cpdry + do i = 1, ncol call energy_helper_eam_def_column(u(i,:),v(i,:),T(i,:),q(i,1:pver,1:pcnst),& ps(i),pdel(i,:),phis(i), & - ke(i),se(i),wv(i),wl(i),wi(i),wr(i),ws(i),te(i),tw(i) ) + ke(i),se(i),wv(i),wl(i),wi(i),wr(i),ws(i),te(i),tw(i) ) enddo + + endif + else call endrun('energy_helper...column is not implemented if water forms do not exist') endif end subroutine energy_helper_eam_def + + + + + + subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & ke,se,wv,wl,wi,wr,ws,te,tw, & + cpstar, & + qini,cldliqini,cldiceini,rainini,snowini,qdryini, & teloc,psterm) !state vars are of size psetcols,pver, so, not exactly correct @@ -1169,6 +1279,14 @@ subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & real(r8), intent(in) :: pdel(pver) real(r8), intent(in) :: phis + real(r8), intent(in ), optional :: qini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: cldliqini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: cldiceini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: rainini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: snowini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: qdryini(pver) ! initial specific humidity + real(r8), intent(in ), optional :: cpstar(pver) ! initial specific humidity + real(r8), intent(inout) :: ke ! vertical integral of kinetic energy real(r8), intent(inout) :: se ! vertical integral of static energy real(r8), intent(inout) :: wv ! vertical integral of water (vapor) @@ -1183,6 +1301,10 @@ subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & real(r8), intent(inout), optional :: psterm integer :: i,k + real(r8) :: cpstar_loc, qdry + + !default value, if not cpstar or computing cpstar from qini + cpstar_loc = cpair ! Compute vertical integrals of dry static energy and water (vapor, liquid, ice) ke = 0._r8 @@ -1197,6 +1319,8 @@ subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & if (present(teloc) .and. present(psterm))then teloc = 0.0; psterm = 0.0 do k = 1, pver +! ADD CPSTAR here too +!right now this code is not used, could be used in local fixers teloc(k) = 0.5_r8*(u(k)**2 + v(k)**2)*pdel(k)/gravit & + t(k)*cpair*pdel(k)/gravit & + (latvap+latice)*q(k,1 )*pdel(k)/gravit @@ -1208,9 +1332,33 @@ subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & do k = 1, pver ke = ke + 0.5_r8*(u(k)**2 + v(k)**2)*pdel(k)/gravit - se = se + t(k)*cpair*pdel(k)/gravit + + if (.not.use_cpstar) then + se = se + t(k)*cpair*pdel(k)/gravit + else + + if(present(cpstar))then + cpstar_loc = cpstar(k) + elseif (present(qini))then + !cpstar based on qini etc + cpstar_loc = cpair*qdryini(k) + cpwv*qini(k) + cpliq*( cldliqini(k) + rainini(k) ) + & + cpice*( cldiceini(k) + snowini(k) ) +#if 1 +!do not use such option + else + !cpstar based on current q + qdry = 1.0_r8 - q(k,1) - q(k,icldliq) - q(k,icldice) - q(k,irain) - q(k,isnow) + cpstar_loc = cpair*qdry + cpwv*q(k,1) + cpliq*( q(k,icldliq) + q(k,irain) ) + & + cpice*( q(k,icldice) + q(k,isnow) ) +#endif + endif + + se = se + t(k)*cpstar_loc*pdel(k)/gravit + endif + wv = wv + q(k,1 )*pdel(k)/gravit end do + se = se + phis*ps/gravit do k = 1, pver @@ -1224,6 +1372,7 @@ subroutine energy_helper_eam_def_column(u,v,T,q,ps,pdel,phis, & end do ! Compute vertical integrals of frozen static energy and total water. + te = se + ke + (latvap+latice)*wv + latice*( wl + wr ) tw = wv + wl + wi + wr + ws @@ -1231,5 +1380,152 @@ end subroutine energy_helper_eam_def_column +!use: state2 with TE2, +!state2 += ttend/dtime ==> state2 now has energy TE1 + subroutine fixer_semi_pb(teloc1, teloc2, psterm1, psterm2, pdel, ttend) + +!state vars are of size psetcols,pver, so, not exactly correct + real(r8), intent(inout) :: ttend(pver) + real(r8), intent(in) :: pdel(pver) + + real(r8), intent(in) :: teloc1(pver) + real(r8), intent(in) :: teloc2(pver) + real(r8), intent(in) :: psterm1 + real(r8), intent(in) :: psterm2 + + integer :: i,k + real(r8) :: fq + + !compute temperature tend from PW adjustment + !keep is as local as possible + ttend(1:pver)=0.0 + !first, tendency from terms at each vertical level + fq=0.0 + do k=1,pver + ttend(k)=(teloc1(k)-teloc2(k))*gravit/cpair/pdel(k) + !sum pdel into fq for the next, boundary (or ps) term + fq=fq+pdel(k) + enddo + + !second, tendency from ps term is peanutbuttered to each level + ttend(1:pver) = ttend(1:pver) + (psterm1-psterm2)*gravit/cpair/fq + + end subroutine fixer_semi_pb + + subroutine fixer_pb(teloc1, teloc2, psterm1, psterm2, pdel, ttend) + +!state vars are of size psetcols,pver, so, not exactly correct + real(r8), intent(inout) :: ttend(pver) + real(r8), intent(in) :: pdel(pver) + real(r8), intent(in) :: teloc1(pver) + real(r8), intent(in) :: teloc2(pver) + real(r8), intent(in) :: psterm1 + real(r8), intent(in) :: psterm2 + + real(r8) :: zeroarray(pver) + + zeroarray(:) = 0.0 + call fixer_semi_pb(zeroarray, zeroarray, & + psterm1+sum(teloc1), psterm2+sum(teloc2), pdel, ttend) + + end subroutine fixer_pb + +!state1 with te1 +!state2 with pdel +!returns ttend st state2 with ttend and pdel adds deltate + subroutine fixer_pb_simple(deltate, pdel, ttend) + +!state vars are of size psetcols,pver, so, not exactly correct + real(r8), intent(inout) :: ttend(pver) + real(r8), intent(in) :: pdel(pver) + real(r8), intent(in) :: deltate + + real(r8) :: zeroarray(pver) + + zeroarray(:) = 0.0 + call fixer_semi_pb(zeroarray, zeroarray, & + deltate, 0.0_r8, pdel, ttend) + + end subroutine fixer_pb_simple + + + subroutine dme_adjust(state, qini, dt) + + implicit none + type(physics_state), intent(inout) :: state + real(r8), intent(in ) :: qini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: dt ! model physics timestep + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: i,k,m ! Longitude, level indices + real(r8) :: fdq(pcols) ! mass adjustment factor + + lchnk = state%lchnk + ncol = state%ncol + + ! adjust dry mass in each layer back to input value, while conserving + ! constituents, momentum, and total energy + do k = 1, pver + +!!! ! adjusment factor is just change in water vapor + fdq(:ncol) = 1._r8 + state%q(:ncol,k,1) - qini(:ncol,k) + + ! adjust constituents to conserve mass in each layer + do m = 1, pcnst + state%q(:ncol,k,m) = state%q(:ncol,k,m) / fdq(:ncol) + end do + +! compute new total pressure variables + state%pdel (:ncol,k ) = state%pdel(:ncol,k ) * fdq(:ncol) + state%pint (:ncol,k+1) = state%pint(:ncol,k ) + state%pdel(:ncol,k) + end do + state%ps(:ncol) = state%pint (:ncol,pver+1) + + end subroutine dme_adjust + + + subroutine dme_adjust_wl(state, qini, cldliqini, cldiceini, rainini, snowini, dt) + + implicit none + type(physics_state), intent(inout) :: state + real(r8), intent(in ) :: qini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: cldliqini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: cldiceini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: rainini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: snowini(pcols,pver) ! initial specific humidity + real(r8), intent(in ) :: dt ! model physics timestep + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: i,k,m ! Longitude, level indices + real(r8) :: fdq(pcols) ! mass adjustment factor + + lchnk = state%lchnk + ncol = state%ncol + + ! adjust dry mass in each layer back to input value, while conserving + ! constituents, momentum, and total energy + do k = 1, pver + + fdq(:ncol) = 1._r8 + state%q(:ncol,k,1) - qini(:ncol,k) & + + state%q(:ncol,k,icldliq) - cldliqini(:ncol,k) & + + state%q(:ncol,k,icldice) - cldiceini(:ncol,k) & + + state%q(:ncol,k,irain) - rainini(:ncol,k) & + + state%q(:ncol,k,isnow) - snowini(:ncol,k) + ! adjust constituents to conserve mass in each layer + do m = 1, pcnst + state%q(:ncol,k,m) = state%q(:ncol,k,m) / fdq(:ncol) + end do + +! compute new total pressure variables + state%pdel (:ncol,k ) = state%pdel(:ncol,k ) * fdq(:ncol) + state%pint (:ncol,k+1) = state%pint(:ncol,k ) + state%pdel(:ncol,k) + end do + state%ps(:ncol) = state%pint (:ncol,pver+1) + + end subroutine dme_adjust_wl + + + + end module check_energy diff --git a/components/eam/src/physics/cam/phys_control.F90 b/components/eam/src/physics/cam/phys_control.F90 index 82a351319fbf..a1ec8f5fc92e 100644 --- a/components/eam/src/physics/cam/phys_control.F90 +++ b/components/eam/src/physics/cam/phys_control.F90 @@ -145,6 +145,17 @@ module phys_control !GW energy fix logical, public, protected :: use_gw_energy_fix = .false. +!WL, cpstar, enthalpies +!all set in eam guts and via namelists, their values are printed at init +logical, public, protected :: use_waterloading = .false. +logical, public, protected :: use_cpstar = .false. +logical, public, protected :: use_enthalpy_cpdry = .false. +logical, public, protected :: use_enthalpy_cl = .false. +logical, public, protected :: use_enthalpy_cpwv = .false. +logical, public, protected :: use_enthalpy_theoretical = .false. +logical, public, protected :: use_global_cpterms_dme = .false. +logical, public, protected :: dycore_fixer_only = .false. + ! Switches that turn on/off individual parameterizations. ! ! Comment by Hui Wan (PNNL, 2014-12): @@ -179,7 +190,7 @@ subroutine phys_ctl_readnl(nlfile) use namelist_utils, only: find_group_name use units, only: getunit, freeunit use mpishorthand - use cam_control_mod, only: cam_ctrl_set_physics_type + use cam_control_mod, only: cam_ctrl_set_physics_type, cam_ctrl_waterloading use physconst, only: pi character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input @@ -204,6 +215,8 @@ subroutine phys_ctl_readnl(nlfile) print_fixer_message, & use_hetfrz_classnuc, use_gw_oro, use_gw_front, use_gw_convect, & use_gw_energy_fix, & + use_waterloading, use_cpstar, use_enthalpy_cpdry, use_enthalpy_cl, use_enthalpy_cpwv, & + use_enthalpy_theoretical, use_global_cpterms_dme, dycore_fixer_only, & cld_macmic_num_steps, micro_do_icesupersat, & fix_g1_err_ndrop, ssalt_tuning, resus_fix, convproc_do_aer, & convproc_do_gas, convproc_method_activate, liqcf_fix, regen_fix, demott_ice_nuc, pergro_mods, pergro_test_active, & @@ -276,7 +289,15 @@ subroutine phys_ctl_readnl(nlfile) call mpibcast(use_gw_oro, 1 , mpilog, 0, mpicom) call mpibcast(use_gw_front, 1 , mpilog, 0, mpicom) call mpibcast(use_gw_convect, 1 , mpilog, 0, mpicom) - call mpibcast(use_gw_energy_fix, 1 , mpilog, 0, mpicom) + call mpibcast(use_gw_energy_fix, 1 , mpilog, 0, mpicom) + call mpibcast(use_waterloading, 1 , mpilog, 0, mpicom) + call mpibcast(use_cpstar, 1 , mpilog, 0, mpicom) + call mpibcast(use_enthalpy_cpdry, 1 , mpilog, 0, mpicom) + call mpibcast(use_enthalpy_cl, 1 , mpilog, 0, mpicom) + call mpibcast(use_enthalpy_cpwv, 1 , mpilog, 0, mpicom) + call mpibcast(use_enthalpy_theoretical, 1 , mpilog, 0, mpicom) + call mpibcast(use_global_cpterms_dme, 1 , mpilog, 0, mpicom) + call mpibcast(dycore_fixer_only, 1 , mpilog, 0, mpicom) call mpibcast(fix_g1_err_ndrop, 1 , mpilog, 0, mpicom) call mpibcast(ssalt_tuning, 1 , mpilog, 0, mpicom) call mpibcast(resus_fix, 1 , mpilog, 0, mpicom) @@ -312,6 +333,8 @@ subroutine phys_ctl_readnl(nlfile) call cam_ctrl_set_physics_type(cam_physpkg) + cam_ctrl_waterloading = use_waterloading + ! Error checking: ! Defaults for PBL and microphysics are set in build-namelist. Check here that diff --git a/components/eam/src/physics/cam/physics_types.F90 b/components/eam/src/physics/cam/physics_types.F90 index ce9282375f4f..5a408b9a3e5c 100644 --- a/components/eam/src/physics/cam/physics_types.F90 +++ b/components/eam/src/physics/cam/physics_types.F90 @@ -7,14 +7,14 @@ module physics_types use shr_kind_mod, only: r8 => shr_kind_r8 use ppgrid, only: pcols, pver, psubcols - use constituents, only: pcnst, qmin, cnst_name, icldliq, icldice + use constituents, only: pcnst, qmin, cnst_name, icldliq, icldice, irain, isnow use geopotential, only: geopotential_t use physconst, only: zvir, gravit, cpair, rair, cpairv, rairv use dycore, only: dycore_is use phys_grid, only: get_ncols_p, get_rlon_all_p, get_rlat_all_p, get_gcol_all_p use cam_logfile, only: iulog use cam_abortutils, only: endrun - use phys_control, only: waccmx_is, use_mass_borrower + use phys_control, only: waccmx_is, use_mass_borrower, use_waterloading, use_cpstar use shr_const_mod,only: shr_const_rwv use perf_mod, only: t_startf, t_stopf @@ -68,6 +68,7 @@ module physics_types lat, &! latitude (radians) lon, &! longitude (radians) ps, &! surface pressure + oldps, &! surface pressure psdry, &! dry surface pressure phis, &! surface geopotential ulat, &! unique latitudes (radians) @@ -81,16 +82,19 @@ module physics_types pmid, &! midpoint pressure (Pa) pmiddry, &! midpoint pressure dry (Pa) pdel, &! layer thickness (Pa) + oldpdel, &! layer thickness (Pa) pdeldry, &! layer thickness dry (Pa) rpdel, &! reciprocal of layer thickness (Pa) rpdeldry,&! recipricol layer thickness dry (Pa) lnpmid, &! ln(pmid) lnpmiddry,&! log midpoint pressure dry (Pa) exner, &! inverse exner function w.r.t. surface pressure (ps/p)^(R/cp) + cpstar, & zm ! geopotential height above surface at midpoints (m) real(r8), dimension(:,:,:),allocatable :: & - q ! constituent mixing ratio (kg/kg moist or dry air depending on type) + q, & ! constituent mixing ratio (kg/kg moist or dry air depending on type) + oldq ! constituent mixing ratio (kg/kg moist or dry air depending on type) real(r8), dimension(:,:),allocatable :: & pint, &! interface pressure (Pa) @@ -103,6 +107,8 @@ module physics_types te_ini, &! vertically integrated total (kinetic + static) energy of initial state te_cur, &! vertically integrated total (kinetic + static) energy of current state tw_ini, &! vertically integrated total water of initial state + cpterme, cptermp, pw, pwvapor, & + dvapor, dliquid, dice, qflx, liqflx, iceflx, & tw_cur ! vertically integrated total water of new state integer :: count ! count of values with significant energy or water imbalances integer, dimension(:),allocatable :: & @@ -203,7 +209,7 @@ subroutine physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, psetcol end subroutine physics_type_alloc !=============================================================================== - subroutine physics_update_main(state, ptend, dt, tend) + subroutine physics_update_main(state, ptend, dt, tend, isinterface) !----------------------------------------------------------------------- ! Update the state and or tendency structure with the parameterization tendencies !----------------------------------------------------------------------- @@ -223,6 +229,7 @@ subroutine physics_update_main(state, ptend, dt, tend) type(physics_tend ), intent(inout), optional :: tend ! Physics tendencies over timestep ! This is usually only needed by calls from physpkg. + logical, intent(in) :: isinterface ! !---------------------------Local storage------------------------------- integer :: i,k,m ! column,level,constituent indices @@ -236,6 +243,7 @@ subroutine physics_update_main(state, ptend, dt, tend) real(r8) :: zvirv(state%psetcols,pver) ! Local zvir array pointer real(r8) :: rairv_loc(state%psetcols,pver) + real(r8) :: cpstar_loc(state%psetcols) ! PERGRO limits cldliq/ice for macro/microphysics: character(len=24), parameter :: pergro_cldlim_names(4) = & @@ -247,7 +255,13 @@ subroutine physics_update_main(state, ptend, dt, tend) ! Whether to do validation of state on each call. logical :: state_debug_checks + logical :: use_cpdry_loc + ncol = state%ncol + + use_cpdry_loc = .true. + cpstar_loc(:ncol) = cpair + !----------------------------------------------------------------------- ! The column radiation model does not update the state @@ -281,8 +295,6 @@ subroutine physics_update_main(state, ptend, dt, tend) !----------------------------------------------------------------------- call phys_getopts(state_debug_checks_out=state_debug_checks) - ncol = state%ncol - ! Update u,v fields if(ptend%lu) then do k = ptend%top_level, ptend%bot_level @@ -370,13 +382,27 @@ subroutine physics_update_main(state, ptend, dt, tend) ! Update dry static energy(moved from above for WACCM-X so updating after cpairv_loc update) !------------------------------------------------------------------------------------------- if(ptend%ls) then + + !use_cpstar=true and only for interface calls + use_cpdry_loc = .not. ( isinterface .and. use_cpstar) + + !do k = ptend%top_level, ptend%bot_level + ! if (present(tend)) & + ! tend%dtdt(:ncol,k) = tend%dtdt(:ncol,k) + ptend%s(:ncol,k)/cpair + ! state%t(:ncol,k) = state%t(:ncol,k) + ptend%s(:ncol,k)/cpair * dt + !end do + do k = ptend%top_level, ptend%bot_level + + if(.not. use_cpdry_loc) then + cpstar_loc(:ncol) = state%cpstar(:ncol,k) + endif + if (present(tend)) & - tend%dtdt(:ncol,k) = tend%dtdt(:ncol,k) + ptend%s(:ncol,k)/cpair -! we first assume that dS is really dEn, En=enthalpy=c_p*T, then -! dT = dEn/c_p, so, state%t += ds/c_p. - state%t(:ncol,k) = state%t(:ncol,k) + ptend%s(:ncol,k)/cpair * dt + tend%dtdt(:ncol,k) = tend%dtdt(:ncol,k) + ptend%s(:ncol,k)/cpstar_loc(:ncol) + state%t(:ncol,k) = state%t(:ncol,k) + ptend%s(:ncol,k)/cpstar_loc(:ncol) * dt end do + end if if (ptend%ls .or. ptend%lq(1)) then @@ -388,6 +414,7 @@ subroutine physics_update_main(state, ptend, dt, tend) state%zi , state%zm ,& ncol) +!keep SE defined via cpdry. It is used in physics and should remain 'dry'. do k = ptend%top_level, ptend%bot_level state%s(:ncol,k) = state%t(:ncol,k )*cpair & + gravit*state%zm(:ncol,k) + state%phis(:ncol) @@ -502,6 +529,28 @@ subroutine physics_state_check(state, name) varname="state%te_ini", msg=msg) call shr_assert_in_domain(state%te_cur(:ncol), is_nan=.false., & varname="state%te_cur", msg=msg) + + call shr_assert_in_domain(state%cpterme(:ncol), is_nan=.false., & + varname="state%cpterme", msg=msg) + call shr_assert_in_domain(state%cptermp(:ncol), is_nan=.false., & + varname="state%cptermp", msg=msg) + call shr_assert_in_domain(state%pw(:ncol), is_nan=.false., & + varname="state%pw", msg=msg) + call shr_assert_in_domain(state%pwvapor(:ncol), is_nan=.false., & + varname="state%pwvapor", msg=msg) + + call shr_assert_in_domain(state%cpstar(:ncol,:), is_nan=.false., & + varname="state%cpstar", msg=msg) + + call shr_assert_in_domain(state%oldps(:ncol), is_nan=.false., & + varname="state%ps", msg=msg) + call shr_assert_in_domain(state%oldpdel(:ncol,:), is_nan=.false., & + varname="state%pdel", msg=msg) + call shr_assert_in_domain(state%oldq(:ncol,:,:), is_nan=.false., & + varname="state%q", msg=msg) + + + call shr_assert_in_domain(state%tw_ini(:ncol), is_nan=.false., & varname="state%tw_ini", msg=msg) call shr_assert_in_domain(state%tw_cur(:ncol), is_nan=.false., & @@ -576,6 +625,33 @@ subroutine physics_state_check(state, name) varname="state%te_ini", msg=msg) call shr_assert_in_domain(state%te_cur(:ncol), lt=posinf_r8, gt=neginf_r8, & varname="state%te_cur", msg=msg) + +!not sure which criteria to use here + call shr_assert_in_domain(state%cpterme(:ncol), lt=posinf_r8, gt=neginf_r8, & + varname="state%cpterme", msg=msg) + call shr_assert_in_domain(state%cptermp(:ncol), lt=posinf_r8, gt=neginf_r8, & + varname="state%cptermp", msg=msg) + call shr_assert_in_domain(state%pw(:ncol), lt=posinf_r8, gt=neginf_r8, & + varname="state%pw", msg=msg) + call shr_assert_in_domain(state%pwvapor(:ncol), lt=posinf_r8, gt=neginf_r8, & + varname="state%pwvapor", msg=msg) + + call shr_assert_in_domain(state%cpstar(:ncol,:), lt=posinf_r8, gt=neginf_r8, & + varname="state%cpstar", msg=msg) + +!set them where ps and pdel are set + call shr_assert_in_domain(state%oldps(:ncol), lt=posinf_r8, gt=0._r8, & + varname="state%ps", msg=msg) + call shr_assert_in_domain(state%oldpdel(:ncol,:), lt=posinf_r8, gt=neginf_r8, & + varname="state%pdel", msg=msg) + ! 3-D variables + do m = 1,pcnst + call shr_assert_in_domain(state%oldq(:ncol,:,m), lt=posinf_r8, gt=neginf_r8, & + varname="state%q ("//trim(cnst_name(m))//")", msg=msg) + end do + + + call shr_assert_in_domain(state%tw_ini(:ncol), lt=posinf_r8, gt=neginf_r8, & varname="state%tw_ini", msg=msg) call shr_assert_in_domain(state%tw_cur(:ncol), lt=posinf_r8, gt=neginf_r8, & @@ -1235,9 +1311,20 @@ subroutine physics_state_copy(state_in, state_out) state_out%lat(i) = state_in%lat(i) state_out%lon(i) = state_in%lon(i) state_out%ps(i) = state_in%ps(i) + state_out%oldps(i) = state_in%oldps(i) state_out%phis(i) = state_in%phis(i) state_out%te_ini(i) = state_in%te_ini(i) state_out%te_cur(i) = state_in%te_cur(i) + state_out%cptermp(i) = state_in%cptermp(i) + state_out%cpterme(i) = state_in%cpterme(i) + state_out%dvapor(i) = state_in%dvapor(i) + state_out%dliquid(i) = state_in%dliquid(i) + state_out%dice(i) = state_in%dice(i) + state_out%qflx(i) = state_in%qflx(i) + state_out%liqflx(i) = state_in%liqflx(i) + state_out%iceflx(i) = state_in%iceflx(i) + state_out%pw(i) = state_in%pw(i) + state_out%pwvapor(i) = state_in%pwvapor(i) state_out%tw_ini(i) = state_in%tw_ini(i) state_out%tw_cur(i) = state_in%tw_cur(i) end do @@ -1251,6 +1338,8 @@ subroutine physics_state_copy(state_in, state_out) state_out%omega(i,k) = state_in%omega(i,k) state_out%pmid(i,k) = state_in%pmid(i,k) state_out%pdel(i,k) = state_in%pdel(i,k) + state_out%oldpdel(i,k) = state_in%oldpdel(i,k) + state_out%cpstar(i,k) = state_in%cpstar(i,k) state_out%rpdel(i,k) = state_in%rpdel(i,k) state_out%lnpmid(i,k) = state_in%lnpmid(i,k) state_out%exner(i,k) = state_in%exner(i,k) @@ -1289,6 +1378,7 @@ subroutine physics_state_copy(state_in, state_out) do k = 1, pver do i = 1, ncol state_out%q(i,k,m) = state_in%q(i,k,m) + state_out%oldq(i,k,m) = state_in%oldq(i,k,m) end do end do end do @@ -1374,9 +1464,17 @@ subroutine set_state_pdry (state,pdeld_calc) state%pintdry(:ncol,1) = state%pint(:ncol,1) if (do_pdeld_calc) then + if(use_waterloading)then + do k = 1, pver + state%pdeldry(:ncol,k) = state%pdel(:ncol,k)*& + (1._r8-state%q(:ncol,k,1)-state%q(:ncol,k,icldliq)-state%q(:ncol,k,icldice)& + -state%q(:ncol,k,irain)-state%q(:ncol,k,isnow)) + end do + else do k = 1, pver state%pdeldry(:ncol,k) = state%pdel(:ncol,k)*(1._r8-state%q(:ncol,k,1)) end do + endif endif do k = 1, pver state%pintdry(:ncol,k+1) = state%pintdry(:ncol,k)+state%pdeldry(:ncol,k) @@ -1484,7 +1582,12 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%ps(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%ps') - + + + allocate(state%oldps(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%ps') + + allocate(state%psdry(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%psdry') @@ -1520,7 +1623,15 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%pdel(psetcols,pver), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%pdel') - + + + allocate(state%oldpdel(psetcols,pver), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%pdel') + + allocate(state%cpstar(psetcols,pver), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cpstar') + + allocate(state%pdeldry(psetcols,pver), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%pdeldry') @@ -1544,7 +1655,12 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%q(psetcols,pver,pcnst), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%q') - + + + allocate(state%oldq(psetcols,pver,pcnst), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%q') + + allocate(state%pint(psetcols,pver+1), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%pint') @@ -1566,6 +1682,34 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%te_cur(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%te_cur') + + allocate(state%cptermp(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cpp') + + allocate(state%cpterme(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cpe') + + allocate(state%dvapor(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + allocate(state%dliquid(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + allocate(state%dice(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + allocate(state%qflx(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + allocate(state%liqflx(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + allocate(state%iceflx(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation') + + + allocate(state%pw(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cpe') + + allocate(state%pwvapor(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%cpe') + + allocate(state%tw_ini(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%tw_ini') @@ -1586,6 +1730,7 @@ subroutine physics_state_alloc(state,lchnk,psetcols) state%ulat(:) = inf state%ulon(:) = inf state%ps(:) = inf + state%oldps(:) = 0.0 state%psdry(:) = inf state%phis(:) = inf state%t(:,:) = inf @@ -1596,6 +1741,8 @@ subroutine physics_state_alloc(state,lchnk,psetcols) state%pmid(:,:) = inf state%pmiddry(:,:) = inf state%pdel(:,:) = inf + state%oldpdel(:,:) = 0.0 + state%cpstar(:,:) = cpair state%pdeldry(:,:) = inf state%rpdel(:,:) = inf state%rpdeldry(:,:) = inf @@ -1604,6 +1751,7 @@ subroutine physics_state_alloc(state,lchnk,psetcols) state%exner(:,:) = inf state%zm(:,:) = inf state%q(:,:,:) = inf + state%oldq(:,:,:) = 0.0 state%pint(:,:) = inf state%pintdry(:,:) = inf @@ -1616,6 +1764,18 @@ subroutine physics_state_alloc(state,lchnk,psetcols) state%tw_ini(:) = inf state%tw_cur(:) = inf + state%cpterme(:) = 0.0 + state%cptermp(:) = 0.0 + + state%dvapor(:) = 0.0 + state%dliquid(:) = 0.0 + state%dice(:) = 0.0 + state%qflx(:) = 0.0 + state%liqflx(:) = 0.0 + state%iceflx(:) = 0.0 + state%pw(:) = 0.0 + state%pwvapor(:) = 0.0 + end subroutine physics_state_alloc !=============================================================================== @@ -1636,6 +1796,11 @@ subroutine physics_state_dealloc(state) deallocate(state%ps, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%ps') + + deallocate(state%oldps, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%ps') + + deallocate(state%psdry, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%psdry') @@ -1672,6 +1837,13 @@ subroutine physics_state_dealloc(state) deallocate(state%pdel, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%pdel') + + deallocate(state%oldpdel, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%pdel') + + deallocate(state%cpstar, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%cpstar') + deallocate(state%pdeldry, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%pdeldry') @@ -1696,6 +1868,11 @@ subroutine physics_state_dealloc(state) deallocate(state%q, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%q') + + deallocate(state%oldq, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%q') + + deallocate(state%pint, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%pint') @@ -1717,6 +1894,33 @@ subroutine physics_state_dealloc(state) deallocate(state%te_cur, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%te_cur') + + deallocate(state%cpterme, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%cpterme') + deallocate(state%cptermp, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%cptermp') + + deallocate(state%dvapor, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + deallocate(state%dliquid, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + deallocate(state%dice, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + deallocate(state%qflx, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + deallocate(state%liqflx, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + deallocate(state%iceflx, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation') + + + + deallocate(state%pw, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%cptermp') + deallocate(state%pwvapor, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%cptermp') + + deallocate(state%tw_ini, stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_dealloc error: deallocation error for state%tw_ini') diff --git a/components/eam/src/physics/cam/physics_update_mod.F90 b/components/eam/src/physics/cam/physics_update_mod.F90 index d023f0bd237a..3459451b2d2e 100644 --- a/components/eam/src/physics/cam/physics_update_mod.F90 +++ b/components/eam/src/physics/cam/physics_update_mod.F90 @@ -84,7 +84,7 @@ end subroutine physics_update_init !---------------------------------------------------------------------------- !---------------------------------------------------------------------------- - subroutine physics_update(state, ptend, dt, tend) + subroutine physics_update(state, ptend, dt, tend, isinterface) !purpose: This subroutine calls physics_update_main (old physics_update) !and also output variables for pergro test @@ -98,7 +98,8 @@ subroutine physics_update(state, ptend, dt, tend) !optional arguments type(physics_tend ), intent(inout), optional :: tend ! Physics tendencies over timestep - + logical, intent(in), optional :: isinterface + !Local vars character(len = fieldname_len) :: pname, varname, vsuffix @@ -120,8 +121,12 @@ subroutine physics_update(state, ptend, dt, tend) if (.not. (any(ptend%lq(:)) .or. ptend%ls .or. ptend%lu .or. ptend%lv)) outfld_active = .false. !call the old physics update call - call physics_update_main (state, ptend, dt, tend) - + if(present(isinterface).and.isinterface)then + call physics_update_main (state, ptend, dt, tend, .true.) + else + call physics_update_main (state, ptend, dt, tend, .false.) + endif + if (pergro_test_active .and. outfld_active) then !write text file to be used for the post processing diff --git a/components/eam/src/physics/cam/physpkg.F90 b/components/eam/src/physics/cam/physpkg.F90 index 11f238ffe32a..594c215dc75b 100644 --- a/components/eam/src/physics/cam/physpkg.F90 +++ b/components/eam/src/physics/cam/physpkg.F90 @@ -12,10 +12,10 @@ module physpkg ! July 2015 B. Singh Added code for unified convective transport !----------------------------------------------------------------------- - + use shr_infnan_mod, only: shr_infnan_isnan use shr_kind_mod, only: i8 => shr_kind_i8, r8 => shr_kind_r8 use spmd_utils, only: masterproc - use physconst, only: latvap, latice, rh2o + use physconst, only: latvap, latice, rh2o, cpair, cpwv, cpliq, cpice, gravit use physics_types, only: physics_state, physics_tend, physics_state_set_grid, & physics_ptend, physics_tend_init, & physics_type_alloc, physics_ptend_dealloc,& @@ -29,7 +29,10 @@ module physpkg use camsrfexch, only: cam_out_t, cam_in_t use cam_control_mod, only: ideal_phys, adiabatic - use phys_control, only: phys_do_flux_avg, phys_getopts, waccmx_is + use phys_control, only: phys_do_flux_avg, phys_getopts, waccmx_is, & + use_waterloading, use_cpstar, use_enthalpy_cpdry, use_enthalpy_cpwv, & + use_enthalpy_cl, use_enthalpy_theoretical, use_global_cpterms_dme, & + dycore_fixer_only use zm_conv, only: do_zmconv_dcape_ull => trigdcape_ull, & do_zmconv_dcape_only => trig_dcape_only use scamMod, only: single_column, scm_crm_mode @@ -47,6 +50,10 @@ module physpkg use modal_aero_wateruptake, only: modal_aero_wateruptake_init, & modal_aero_wateruptake_reg + use check_energy, only: dme_adjust, dme_adjust_wl, energy_helper_eam_def, energy_helper_eam_def_column,& + fixer_pb_simple + + implicit none private @@ -57,6 +64,9 @@ module physpkg integer :: qini_idx = 0 integer :: cldliqini_idx = 0 integer :: cldiceini_idx = 0 + integer :: rainini_idx = 0 + integer :: snowini_idx = 0 + integer :: qdryini_idx = 0 integer :: static_ener_ac_idx = 0 integer :: water_vap_ac_idx = 0 @@ -197,6 +207,9 @@ subroutine phys_register call pbuf_add_field('QINI', 'physpkg', dtype_r8, (/pcols,pver/), qini_idx) call pbuf_add_field('CLDLIQINI', 'physpkg', dtype_r8, (/pcols,pver/), cldliqini_idx) call pbuf_add_field('CLDICEINI', 'physpkg', dtype_r8, (/pcols,pver/), cldiceini_idx) + call pbuf_add_field('RAININI', 'physpkg', dtype_r8, (/pcols,pver/), rainini_idx) + call pbuf_add_field('SNOWINI', 'physpkg', dtype_r8, (/pcols,pver/), snowini_idx) + call pbuf_add_field('QDRYINI', 'physpkg', dtype_r8, (/pcols,pver/), qdryini_idx) call pbuf_add_field('static_ener_ac', 'global', dtype_r8, (/pcols/), static_ener_ac_idx) call pbuf_add_field('water_vap_ac', 'global', dtype_r8, (/pcols/), water_vap_ac_idx) @@ -932,6 +945,9 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out ) !BSINGH - addfld and adddefault calls for perturb growth testing if(pergro_test_active)call add_fld_default_calls() + !print WL- and cpstar-related variables + call print_wl_cpstar_info() + end subroutine phys_init ! @@ -1170,6 +1186,7 @@ subroutine phys_run1_adiabatic_or_ideal(ztodt, phys_state, phys_tend, pbuf2d) if (dycore_is('LR') .or. dycore_is('SE') ) then call check_energy_fix(phys_state(c), ptend(c), nstep, flx_heat) + !not ready for cpstar call physics_update(phys_state(c), ptend(c), ztodt, phys_tend(c)) call check_energy_chng(phys_state(c), phys_tend(c), "chkengyfix", nstep, ztodt, & zero, zero, zero, flx_heat) @@ -1503,6 +1520,9 @@ subroutine tphysac (ztodt, cam_in, & real(r8), pointer, dimension(:,:) :: qini real(r8), pointer, dimension(:,:) :: cldliqini real(r8), pointer, dimension(:,:) :: cldiceini + real(r8), pointer, dimension(:,:) :: rainini + real(r8), pointer, dimension(:,:) :: snowini + real(r8), pointer, dimension(:,:) :: qdryini real(r8), pointer, dimension(:,:) :: dtcore real(r8), pointer, dimension(:,:) :: ast ! relative humidity cloud fraction @@ -1523,6 +1543,18 @@ subroutine tphysac (ztodt, cam_in, & logical :: l_gw_drag logical :: l_ac_energy_chk +!ogdef + real(r8) :: ke(pcols), se(pcols), wv(pcols), wl(pcols), wi(pcols), te(pcols),& + tw(pcols), wr(pcols), ws(pcols), cpterm(pcols), & + te_before_pw(pcols), te_after_pw(pcols), deltat(pcols) + + real(r8) :: small_ttend(pcols,pver) + + real(r8) :: keloc, seloc, wvloc, wlloc, wiloc, teloc1, teloc2,& + twloc, wrloc, wsloc + integer :: ic + real(r8) :: qdry(pver), cpstar(pver), qvmass(pcols) + ! !----------------------------------------------------------------------- ! @@ -1562,6 +1594,9 @@ subroutine tphysac (ztodt, cam_in, & call pbuf_get_field(pbuf, qini_idx, qini) call pbuf_get_field(pbuf, cldliqini_idx, cldliqini) call pbuf_get_field(pbuf, cldiceini_idx, cldiceini) + call pbuf_get_field(pbuf, rainini_idx, rainini) + call pbuf_get_field(pbuf, snowini_idx, snowini) + call pbuf_get_field(pbuf, qdryini_idx, qdryini) ifld = pbuf_get_index('CLD') call pbuf_get_field(pbuf, ifld, cld, start=(/1,1,itim_old/),kount=(/pcols,pver,1/)) @@ -1620,25 +1655,25 @@ subroutine tphysac (ztodt, cam_in, & ! Test tracers call tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_tracers_chng(state, tracerint, "tracers_timestep_tend", nstep, ztodt, & cam_in%cflx) call aoa_tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_tracers_chng(state, tracerint, "aoa_tracers_timestep_tend", nstep, ztodt, & cam_in%cflx) ! add tendency from aircraft emissions call co2_cycle_set_ptend(state, pbuf, ptend) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! Chemistry calculation if (chem_is_active()) then call chem_timestep_tend(state, ptend, cam_in, cam_out, ztodt, & pbuf, fh2o, fsds) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, & cam_in%cflx) @@ -1660,7 +1695,7 @@ subroutine tphysac (ztodt, cam_in, & call clubb_surface ( state, ptend, ztodt, cam_in, surfric, obklen) ! Update surface flux constituents - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) else if (l_vdiff) then @@ -1678,7 +1713,7 @@ subroutine tphysac (ztodt, cam_in, & call mspd_intr (ztodt ,state ,ptend) endif - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call t_stopf ('vertical_diffusion_tend') end if ! l_vdiff @@ -1691,7 +1726,7 @@ subroutine tphysac (ztodt, cam_in, & !=================================================== call t_startf('rayleigh_friction') call rayleigh_friction_tend( ztodt, state, ptend) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call t_stopf('rayleigh_friction') if (do_clubb_sgs) then @@ -1710,7 +1745,7 @@ subroutine tphysac (ztodt, cam_in, & ! aerosol dry deposition processes call t_startf('aero_drydep') call aero_model_drydep( state, pbuf, obklen, surfric, cam_in, ztodt, cam_out, ptend ) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call t_stopf('aero_drydep') !--------------------------------------------------------------------------------- @@ -1728,14 +1763,14 @@ subroutine tphysac (ztodt, cam_in, & call gw_tend(state, sgh, pbuf, ztodt, ptend, cam_in) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! Check energy integrals call check_energy_chng(state, tend, "gwdrag", nstep, ztodt, zero, zero, zero, zero) call t_stopf('gw_tend') ! QBO relaxation call qbo_relax(state, pbuf, ptend) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! Check energy integrals call check_energy_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) @@ -1754,7 +1789,7 @@ subroutine tphysac (ztodt, cam_in, & call ionos_intr(state, ptend, pbuf, ztodt) endif - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! Check energy integrals call check_energy_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) call t_stopf ( 'iondrag' ) @@ -1766,13 +1801,13 @@ subroutine tphysac (ztodt, cam_in, & !=================================================== if((Nudge_Model).and.(Nudge_ON)) then call nudging_timestep_tend(state,ptend) - call physics_update(state,ptend,ztodt,tend) + call physics_update(state,ptend,ztodt,tend, isinterface=.true.) endif if (l_ac_energy_chk) then !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv - call pbuf_set_field(pbuf, teout_idx, state%te_cur, (/1,itim_old/),(/pcols,1/)) +! call pbuf_set_field(pbuf, teout_idx, state%te_cur, (/1,itim_old/),(/pcols,1/)) tmp_t(:ncol,:pver) = state%t(:ncol,:pver) @@ -1789,11 +1824,203 @@ subroutine tphysac (ztodt, cam_in, & if ( dycore_is('LR') .or. dycore_is('SE')) call set_dry_to_wet(state) ! Physics had dry, dynamics wants moist + +!!!!!!!!!!!!!!!!!!!!! CODE related to EFLX +! cam_out%precc (i) = prec_dp(i) + prec_sh(i) +! cam_out%precl (i) = prec_sed(i) + prec_pcw(i) +! cam_out%precsc(i) = snow_dp(i) + snow_sh(i) +! cam_out%precsl(i) = snow_sed(i) + snow_pcw(i) +! 263 PRECL m/s 1 A Large-scale (stable) precipitation rate (liq + ice) +! 264 PRECC m/s 1 A Convective precipitation rate (liq + ice) +! 265 PRECT m/s 1 A Total (convective and large-scale) precipitation rate (liq + ice) +! 269 PRECSL m/s 1 A Large-scale (stable) snow rate (water equivalent) +! 270 PRECSC m/s 1 A Convective snow rate (water equivalent) + + + if (use_enthalpy_cpdry) then + state%cpterme(:ncol) = cpair * state%t(:ncol,pver) * cam_in%cflx(:ncol,1) + state%cptermp(:ncol) = 1000.0 * cpair * state%t(:ncol,pver) * cam_out%precl(:ncol) + & + 1000.0 * cpair * state%t(:ncol,pver) * cam_out%precc(:ncol) + + elseif(use_enthalpy_theoretical) then + state%cpterme(:ncol) = cpwv * state%t(:ncol,pver) * cam_in%cflx(:ncol,1) + state%cptermp(:ncol) = & + 1000.0*cpliq*state%t(:ncol,pver)*(cam_out%precl(:ncol)+cam_out%precc(:ncol)-cam_out%precsc(:ncol)-cam_out%precsl(:ncol))+& + 1000.0*cpice*state%t(:ncol,pver)*( cam_out%precsc(:ncol)+cam_out%precsl(:ncol)) + + elseif(use_enthalpy_cl) then + state%cpterme(:ncol) = cpliq * state%t(:ncol,pver) * cam_in%cflx(:ncol,1) + state%cptermp(:ncol) = 1000.0 * cpliq * state%t(:ncol,pver) * cam_out%precl(:ncol) + & + 1000.0 * cpliq * state%t(:ncol,pver) * cam_out%precc(:ncol) + + elseif(use_enthalpy_cpwv) then + state%cpterme(:ncol) = cpwv * state%t(:ncol,pver) * cam_in%cflx(:ncol,1) + state%cptermp(:ncol) = 1000.0 * cpwv * state%t(:ncol,pver) * cam_out%precl(:ncol) + & + 1000.0 * cpwv * state%t(:ncol,pver) * cam_out%precc(:ncol) + endif + + !collect all deltas and fluxes for glob mean + state%qflx(:ncol) = cam_in%cflx(:ncol,1) + state%dvapor(:ncol) = 0.0 + state%dliquid(:ncol) = 0.0 + state%dice(:ncol) = 0.0 + do k=1,pver + state%dvapor(:ncol) = state%dvapor(:ncol) + state%pdel(:ncol,k)*(state%q(:ncol,k,1) - qini(:ncol,k))/gravit + state%dliquid(:ncol) = state%dliquid(:ncol) & + + state%pdel(:ncol,k)*(state%q(:ncol,k,icldliq) - cldliqini(:ncol,k))/gravit & + + state%pdel(:ncol,k)*(state%q(:ncol,k,irain) - rainini(:ncol,k))/gravit + state%dice(:ncol) = state%dice(:ncol) & + + state%pdel(:ncol,k)*(state%q(:ncol,k,icldice) - cldiceini(:ncol,k))/gravit & + + state%pdel(:ncol,k)*(state%q(:ncol,k,isnow) - snowini(:ncol,k))/gravit + enddo + + state%liqflx(:ncol) = 1000.0*(cam_out%precl(:ncol)+cam_out%precc(:ncol)-cam_out%precsc(:ncol)-cam_out%precsl(:ncol)) + state%iceflx(:ncol) = 1000.0*( cam_out%precsc(:ncol)+cam_out%precsl(:ncol)) + + call energy_helper_eam_def(state%u,state%v,state%T,state%q,state%ps,state%pdel,state%phis, & + ke(:ncol),se(:ncol),wv(:ncol),wl(:ncol),& + wi(:ncol),wr(:ncol),ws(:ncol),te_before_pw(:ncol),tw(:ncol), & + ncol, & + cpstar=state%cpstar) + + state%te_cur(:ncol) = te_before_pw(:ncol) + + !compute energy of PW (DME adjust) terms + + !save dp, ps, q first + !did mot save pint -- seems to not be used? + state%oldq = state%q + state%oldpdel = state%pdel + state%oldps = state%ps + + !first, compute PW of loading only vapor + call dme_adjust(state, qini, ztodt) + + !compute energy after PW, TE with PW is in te_after_pw + !if cpstar is on, this version should use cpstar, but it does not need to use qini and other init tracer values + call energy_helper_eam_def(state%u,state%v,state%T,state%q,state%ps,state%pdel,state%phis, & + ke(:ncol),se(:ncol),wv(:ncol),wl(:ncol),& + wi(:ncol),wr(:ncol),ws(:ncol),te_after_pw(:ncol),tw(:ncol), & + ncol) + + !compute DME adjust energy vapor only + state%pwvapor(:ncol) = state%te_cur(:ncol) - te_after_pw(:ncol) + !now we can reset all variables after DME adjust + state%q = state%oldq + state%pdel = state%oldpdel + state%ps = state%oldps + + + if(use_waterloading) then + !WL version + call dme_adjust_wl(state, qini, cldliqini, cldiceini, rainini, snowini, ztodt) + else + !no WL version + call dme_adjust(state, qini, ztodt) + endif + + !compute energy after PW, TE with PW is in te_after_pw + !if cpstar is on, this version should use cpstar, but it does not need to use qini and other init tracer values + call energy_helper_eam_def(state%u,state%v,state%T,state%q,state%ps,state%pdel,state%phis, & + ke(:ncol),se(:ncol),wv(:ncol),wl(:ncol),& + wi(:ncol),wr(:ncol),ws(:ncol),te_after_pw(:ncol),tw(:ncol), & + ncol) + + !finally, compute DME adjust energy + state%pw(:ncol) = state%te_cur(:ncol) - te_after_pw(:ncol) + + !units and dt: + !for cpdry, PW has to match cp terms, and it does + !pw ~= (cptermp - cpterme)*ztodt + + !define total cp term for simplicity + cpterm(:ncol) = state%cptermp(:ncol)*ztodt - state%cpterme(:ncol)*ztodt + + !delta is for using a local fixer after transfer of enthalpies + deltat(:ncol) = cpterm(:ncol) - state%pw(:ncol) + +!debug around local fixer, ignore for now +#if 0 + do ic=1,ncol +!what we want, small local fixer + call fixer_pb_simple(deltat(ic), state%pdel(ic,:), small_ttend(ic,:) ) +!pw as a local fixer +! call fixer_pb_simple(state%pw(ic), state%pdel(ic,:), small_ttend(ic,:) ) +!cpterms as local fixer +! call fixer_pb_simple((state%cptermp(ic)-state%cpterme(ic))*ztodt, & +! state%pdel(ic,:), small_ttend(ic,:) ) + +!check + call energy_helper_eam_def_column(state%u(ic,:),state%v(ic,:),& + state%T(ic,:),& + state%q(ic,:,:),state%ps(ic),state%pdel(ic,:),state%phis(ic), & + keloc,seloc,wvloc,wlloc,wiloc,wrloc,wsloc,teloc1,twloc) + call energy_helper_eam_def_column(state%u(ic,:),state%v(ic,:),& + state%T(ic,:)+small_ttend(ic,:),& + state%q(ic,:,:),state%ps(ic),state%pdel(ic,:),state%phis(ic), & + keloc,seloc,wvloc,wlloc,wiloc,wrloc,wsloc,teloc2,twloc) +!this one gets close results +!print *, 'ic, deltat,teloc2-teloc1',deltat(ic), teloc2-teloc1 +!this one gets error of 1e-9 or less +!if(abs(deltat(ic)) > 0.5)then +!print *, 'ic, deltat - (teloc2-teloc1) / ..', (deltat(ic) - (teloc2-teloc1))/deltat(ic) +!endif + + enddo +#endif + + + !now we can reset all variables after DME adjust + state%q = state%oldq + state%pdel = state%oldpdel + state%ps = state%oldps + + !check how fixer looks like without PW in it. + ! that is, when fixer is for dycore only + if(dycore_fixer_only)then + state%te_cur(:ncol) = state%te_cur(:ncol) - state%pw(:ncol) + endif + + !remove cp terms from fixer to mimic enthalpy transfers + if(use_global_cpterms_dme)then + !take CP term out of te_cur + state%te_cur(:ncol) = state%te_cur(:ncol) - state%cptermp(:ncol)*ztodt & + + state%cpterme(:ncol)*ztodt + endif + +!ignore +#if 0 +!LOCAL fixer + tend%dtdt(:ncol,1:pver) = tend%dtdt(:ncol,1:pver) - small_ttend(:ncol,1:pver)/ztodt +#endif + +!ignore +!rescale TTEND with cpstar +#if 0 + do ic=1,ncol + qdry(:) = 1.0 - state%q(ic,:,1) - state%q(ic,:,icldice) - state%q(ic,:,icldliq) & + - state%q(ic,:,irain) - state%q(ic,:,isnow) + cpstar(:) = cpair*qdry(:) + cpwv*state%q(ic,:,1) + cpliq*( state%q(ic,:,icldliq) + state%q(ic,:,irain) ) & + + cpice*( state%q(ic,:,icldice) + state%q(ic,:,isnow) ) + tend%dtdt(ic,1:pver) = tend%dtdt(ic,1:pver) / cpstar(1:pver) * cpair + enddo +#endif + + +!!!!!!!!!!!!!! end of EFLX work + + + + + + call pbuf_set_field(pbuf, teout_idx, state%te_cur, (/1,itim_old/),(/pcols,1/)) + ! Scale dry mass and energy (does nothing if dycore is EUL or SLD) tmp_q (:ncol,:pver) = state%q(:ncol,:pver,1) tmp_cldliq(:ncol,:pver) = state%q(:ncol,:pver,icldliq) tmp_cldice(:ncol,:pver) = state%q(:ncol,:pver,icldice) - call physics_dme_adjust(state, tend, qini, ztodt) + +! call physics_dme_adjust(state, tend, qini, ztodt) !!! REMOVE THIS CALL, SINCE ONLY Q IS BEING ADJUSTED. WON'T BALANCE ENERGY. TE IS SAVED BEFORE THIS !!! call check_energy_chng(state, tend, "drymass", nstep, ztodt, zero, zero, zero, zero) @@ -1992,6 +2219,9 @@ subroutine tphysbc (ztodt, & real(r8), pointer, dimension(:,:) :: qini real(r8), pointer, dimension(:,:) :: cldliqini real(r8), pointer, dimension(:,:) :: cldiceini + real(r8), pointer, dimension(:,:) :: rainini + real(r8), pointer, dimension(:,:) :: snowini + real(r8), pointer, dimension(:,:) :: qdryini real(r8), pointer, dimension(:,:) :: dtcore real(r8), pointer, dimension(:,:,:) :: fracis ! fraction of transported species that are insoluble @@ -2154,6 +2384,9 @@ subroutine tphysbc (ztodt, & call pbuf_get_field(pbuf, qini_idx, qini) call pbuf_get_field(pbuf, cldliqini_idx, cldliqini) call pbuf_get_field(pbuf, cldiceini_idx, cldiceini) + call pbuf_get_field(pbuf, rainini_idx, rainini) + call pbuf_get_field(pbuf, snowini_idx, snowini) + call pbuf_get_field(pbuf, qdryini_idx, qdryini) ifld = pbuf_get_index('DTCORE') call pbuf_get_field(pbuf, ifld, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) @@ -2168,6 +2401,33 @@ subroutine tphysbc (ztodt, & tend %dvdt(:ncol,:pver) = 0._r8 !!== KZ_WCON + + + + + qini (:ncol,:pver) = state%q(:ncol,:pver, 1) + cldliqini(:ncol,:pver) = state%q(:ncol,:pver,icldliq) + cldiceini(:ncol,:pver) = state%q(:ncol,:pver,icldice) + rainini(:ncol,:pver) = state%q(:ncol,:pver,irain) + snowini(:ncol,:pver) = state%q(:ncol,:pver,isnow) + + do i=1,ncol + do k=1,pver + + qdryini(i,k) = 1.0 - state%q(i,k, 1) - state%q(i,k,icldliq) & + - state%q(i,k,icldice) - state%q(i,k,irain) & + - state%q(i,k,isnow) + + !with updated *ini, update cpstar, too + !do it before dycore energy fixer + +!update cpstar in dp_coupling instead +! state%cpstar(i,k) = cpair*qdryini(i,k) + cpwv*qini(i,k) + cpliq*( cldliqini(i,k) + rainini(i,k) ) + & +! cpice*( cldiceini(i,k) + snowini(i,k) ) + + enddo + enddo + call check_qflx (state, tend, "PHYBC01", nstep, ztodt, cam_in%cflx(:,1)) call check_water(state, tend, "PHYBC01", nstep, ztodt) @@ -2256,24 +2516,22 @@ subroutine tphysbc (ztodt, & !=================================================== ! Global mean total energy fixer !=================================================== + if (l_bc_energy_fix) then call t_startf('energy_fixer') tini(:ncol,:pver) = state%t(:ncol,:pver) if (dycore_is('LR') .or. dycore_is('SE')) then + !if(.not. is_first_step())then call check_energy_fix(state, ptend, nstep, flx_heat) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) + !endif call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) end if ! Save state for convective tendency calculations. call diag_conv_tend_ini(state, pbuf) - qini (:ncol,:pver) = state%q(:ncol,:pver, 1) - cldliqini(:ncol,:pver) = state%q(:ncol,:pver,icldliq) - cldiceini(:ncol,:pver) = state%q(:ncol,:pver,icldice) - - call outfld('TEOUT', teout , pcols, lchnk ) call outfld('TEINP', state%te_ini, pcols, lchnk ) call outfld('TEFIX', state%te_cur, pcols, lchnk ) @@ -2311,7 +2569,9 @@ subroutine tphysbc (ztodt, & ptend%s, ptend%q(1,1,1)) ptend%s(:ncol,:) = (ptend%s(:ncol,:) - state%t(:ncol,:) )/ztodt * cpair ptend%q(:ncol,:,1) = (ptend%q(:ncol,:,1) - state%q(:ncol,:,1))/ztodt - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) + + call check_energy_chng(state, tend, "dadadj", nstep, ztodt, zero, zero, zero, zero) call t_stopf('dry_adjustment') @@ -2335,7 +2595,7 @@ subroutine tphysbc (ztodt, & dsubcld, jt, maxg, ideep, lengath) call t_stopf('convect_deep_tend') - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call pbuf_get_field(pbuf, prec_dp_idx, prec_dp ) call pbuf_get_field(pbuf, snow_dp_idx, snow_dp ) @@ -2367,7 +2627,7 @@ subroutine tphysbc (ztodt, & state , ptend , pbuf , sh_e_ed_ratio , sgh, sgh30, cam_in) call t_stopf ('convect_shallow_tend') - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) flx_cnd(:ncol) = prec_sh(:ncol) + rliq2(:ncol) call check_energy_chng(state, tend, "convect_shallow", nstep, ztodt, zero, flx_cnd, snow_sh, zero) @@ -2404,7 +2664,7 @@ subroutine tphysbc (ztodt, & cmfmc, cmfmc2, & cam_in%ts, cam_in%sst, zdu) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "cldwat_tend", nstep, ztodt, zero, prec_str(:ncol), snow_str(:ncol), zero) call t_stopf('stratiform_tend') @@ -2435,7 +2695,7 @@ subroutine tphysbc (ztodt, & call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "mp_aero_tend", nstep, ztodt, zero, zero, zero, zero) endif @@ -2468,7 +2728,7 @@ subroutine tphysbc (ztodt, & ! ptend down by the number of substeps, then applying it for ! the full time (ztodt). call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "macrop_tend", nstep, ztodt, & zero, flx_cnd/cld_macmic_num_steps, & det_ice/cld_macmic_num_steps, flx_heat/cld_macmic_num_steps) @@ -2511,7 +2771,7 @@ subroutine tphysbc (ztodt, & call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) ! Update physics tendencies and copy state to state_eq, because that is ! input for microphysics - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "clubb_tend", nstep, ztodt, & cam_in%cflx(:,1)/cld_macmic_num_steps, flx_cnd/cld_macmic_num_steps, & det_ice/cld_macmic_num_steps, flx_heat/cld_macmic_num_steps) @@ -2564,7 +2824,7 @@ subroutine tphysbc (ztodt, & ! (see above note for macrophysics). call physics_ptend_scale(ptend_sc, 1._r8/cld_macmic_num_steps, ncol) - call physics_update (state_sc, ptend_sc, ztodt, tend_sc) + call physics_update (state_sc, ptend_sc, ztodt, tend_sc, isinterface=.true.) call check_energy_chng(state_sc, tend_sc, "microp_tend_subcol", & nstep, ztodt, zero_sc, prec_str_sc(:ncol)/cld_macmic_num_steps, & snow_str_sc(:ncol)/cld_macmic_num_steps, zero_sc) @@ -2585,7 +2845,7 @@ subroutine tphysbc (ztodt, & ! (see above note for macrophysics). call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) - call physics_update (state, ptend, ztodt, tend) + call physics_update (state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "microp_tend", nstep, ztodt, & zero, prec_str(:ncol)/cld_macmic_num_steps, & snow_str(:ncol)/cld_macmic_num_steps, zero) @@ -2637,13 +2897,13 @@ subroutine tphysbc (ztodt, & sh_e_ed_ratio, mu, md, du, eu, ed, dp, dsubcld, & jt, maxg, ideep, lengath, species_class, & cam_out, pbuf, ptend ) ! outputs - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! deep convective aerosol transport call convect_deep_tend_2( state, ptend, ztodt, pbuf, & mu, eu, du, md, ed, dp, dsubcld, jt, maxg, & ideep, lengath, species_class ) - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) ! check tracer integrals call check_tracers_chng(state, tracerint, "cmfmca", nstep, ztodt, zero_tracers) @@ -2698,7 +2958,7 @@ subroutine tphysbc (ztodt, & do i=1,ncol tend%flx_net(i) = net_flx(i) end do - call physics_update(state, ptend, ztodt, tend) + call physics_update(state, ptend, ztodt, tend, isinterface=.true.) call check_energy_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) call t_stopf('radiation') @@ -2886,4 +3146,43 @@ subroutine add_fld_default_calls() end subroutine add_fld_default_calls + + +subroutine print_wl_cpstar_info() + use spmd_utils, only : masterproc + use cam_abortutils, only : endrun + + if (masterproc) write(iulog,*) 'use_waterloading = ', use_waterloading + if (masterproc) write(iulog,*) 'use_cpstar = ', use_cpstar + if (masterproc) write(iulog,*) 'use_enthalpy_cpdry = ', use_enthalpy_cpdry + if (masterproc) write(iulog,*) 'use_enthalpy_cl = ', use_enthalpy_cl + if (masterproc) write(iulog,*) 'use_enthalpy_cpwv = ', use_enthalpy_cpwv + if (masterproc) write(iulog,*) 'use_enthalpy_theoretical = ', use_enthalpy_theoretical + if (masterproc) write(iulog,*) 'use_global_cpterms_dme = ', use_global_cpterms_dme + if (masterproc) write(iulog,*) 'dycore_fixer_only = ', dycore_fixer_only + + if ((use_cpstar).and.(.not. use_waterloading)) then + call endrun ('PHYS init error: use_cpstar=true requires use_waterloading=true') + endif + + if ((use_enthalpy_cpdry).and.(use_enthalpy_cl.or.use_enthalpy_theoretical.or.use_enthalpy_cpwv)) then + call endrun ('PHYS init error: use_enthalpy* -- only one should be set to true') + endif + + if ((use_enthalpy_cl).and.(use_enthalpy_cpdry.or.use_enthalpy_theoretical.or.use_enthalpy_cpwv)) then + call endrun ('PHYS init error: use_enthalpy* -- only one should be set to true') + endif + + if ((use_enthalpy_theoretical).and.(use_enthalpy_cl.or.use_enthalpy_cpdry.or.use_enthalpy_cpwv)) then + call endrun ('PHYS init error: use_enthalpy* -- only one should be set to true') + endif + + if ((use_enthalpy_cpwv).and.(use_enthalpy_cl.or.use_enthalpy_cpdry.or.use_enthalpy_theoretical)) then + call endrun ('PHYS init error: use_enthalpy* -- only one should be set to true') + endif + +end subroutine print_wl_cpstar_info + + + end module physpkg diff --git a/components/homme/src/share/control_mod.F90 b/components/homme/src/share/control_mod.F90 index f27454f47ecc..978c9025cb44 100644 --- a/components/homme/src/share/control_mod.F90 +++ b/components/homme/src/share/control_mod.F90 @@ -137,6 +137,7 @@ module control_mod ! internally the code should use logical "use_moisture" character(len=MAX_STRING_LEN) , public :: moisture +! OG: PREQX cpstar, not used in eflx work integer, public :: use_cpstar=0 ! use cp or cp* in thermodynamics logical, public :: use_moisture=.false. ! use Q(:,:,:,1) to compute T_v diff --git a/components/homme/src/share/prim_driver_base.F90 b/components/homme/src/share/prim_driver_base.F90 index 36c827abd481..e7a340e16efd 100644 --- a/components/homme/src/share/prim_driver_base.F90 +++ b/components/homme/src/share/prim_driver_base.F90 @@ -1562,6 +1562,11 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) use bfb_mod, only : bfb_pow #endif #endif + +#ifdef CAM + use cam_control_mod, only : cam_ctrl_waterloading +#endif + implicit none type (element_t), intent(inout) :: elem real (kind=real_kind), intent(in) :: dt @@ -1586,6 +1591,31 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) real (kind=real_kind) :: dpnh_dp_i(np,np,nlevp) #endif + integer :: q1ind, q2ind, q3ind, q6ind, q7ind + + logical :: local_waterloading + +#if 0 +1 Q Specific humidity +2 CLDLIQ Grid box averaged cloud liquid amount +3 CLDICE Grid box averaged cloud ice amount +4 NUMLIQ Grid box averaged cloud liquid number +5 NUMICE Grid box averaged cloud ice number +6 RAINQM Grid box averaged rain amount +7 SNOWQM Grid box averaged snow amount +8 NUMRAI Grid box averaged rain number +9 NUMSNO Grid box averaged snow number +#endif + +#ifdef CAM +local_waterloading = cam_ctrl_waterloading +q1ind=1; q2ind=2; q3ind=3; q6ind=6; q7ind=7; +#else +local_waterloading = .false. +q1ind=-1; q2ind=-1; q3ind=-1; q6ind=-1; q7ind=-1; +#endif + + #ifdef HOMMEXX_BFB_TESTING ! BFB comparison with C++ requires to perform the reduction ! of FQ over the whole column *before* adding to ps @@ -1597,11 +1627,11 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) if (dt_remap_factor==0) then adjust_ps=.true. ! stay on reference levels for Eulerian case else -#ifdef SCREAM - adjust_ps=.false. ! Lagrangian case can support adjusting dp3d or ps -#else - adjust_ps=.true. ! Lagrangian case can support adjusting dp3d or ps -#endif +!#ifdef SCREAM +! adjust_ps=.false. ! Lagrangian case can support adjusting dp3d or ps +!#else + adjust_ps=.false. ! Lagrangian case can support adjusting dp3d or ps +!#endif endif #else adjust_ps=.true. ! preqx requires forcing to stay on reference levels @@ -1644,7 +1674,20 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) ! dyn_in%elem(ie)%state%Qdp(i,j,k,q,tl_fQdp) + fq elem%state%Qdp(i,j,k,q,np1_qdp) = & dp(i,j,k)*elem%derived%FQ(i,j,k,q) - + + +!this will only run with EAM, no need to check qsize +if(local_waterloading) then + + if (q==q1ind .or. q==q2ind .or. q==q3ind .or. q==q6ind .or. q==q7ind) then +!repeated code + fq = dp(i,j,k)*( elem%derived%FQ(i,j,k,q) -& + elem%state%Q(i,j,k,q)) + ! force ps to conserve mass: + ps(i,j)=ps(i,j) + fq + dp_adj(i,j,k)=dp_adj(i,j,k) + fq ! ps = ps0+sum(dp(k)) + endif +else if (q==1) then fq = dp(i,j,k)*( elem%derived%FQ(i,j,k,q) -& elem%state%Q(i,j,k,q)) @@ -1656,6 +1699,8 @@ subroutine applyCAMforcing_tracers(elem,hvcoord,np1,np1_qdp,dt,adjustment) #endif dp_adj(i,j,k)=dp_adj(i,j,k) + fq ! ps = ps0+sum(dp(k)) endif +endif + enddo end do end do diff --git a/share/util/shr_reprosum_mod.F90 b/share/util/shr_reprosum_mod.F90 index 9d10b9ebb0d9..05167a301a29 100644 --- a/share/util/shr_reprosum_mod.F90 +++ b/share/util/shr_reprosum_mod.F90 @@ -432,6 +432,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & integer omp_get_max_threads external omp_get_max_threads #endif + + integer :: idxx(nflds), ii + ! !----------------------------------------------------------------------- ! @@ -443,6 +446,8 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & repro_sum_both = 0 nonrepro_sum = 0 + idxx(:) = -1 + ! set MPI communicator if ( present(commid) ) then mpi_comm = commid @@ -471,6 +476,10 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & if (nan_check .or. inf_check) then +do ii=1,nflds +if (any(shr_infnan_isnan(arr(:,ii)))) idxx(ii)=ii +enddo + nan_count = count(shr_infnan_isnan(arr)) inf_count = count(shr_infnan_isinf(arr)) @@ -479,6 +488,9 @@ subroutine shr_reprosum_calc (arr, arr_gsum, nsummands, dsummands, & write(s_logunit,37) real(nan_count,r8), real(inf_count,r8), mypid 37 format("SHR_REPROSUM_CALC: Input contains ",e12.5, & " NaNs and ", e12.5, " INFs on process ", i7) + + print *, 'indices are', idxx + call shr_sys_abort("shr_reprosum_calc ERROR: NaNs or INFs in input") endif