diff --git a/.gitmodules b/.gitmodules index 11bdef17c3..f49ccae8d1 100644 --- a/.gitmodules +++ b/.gitmodules @@ -32,8 +32,8 @@ [submodule "atmos_phys"] path = src/atmos_phys - url = https://github.com/NorESMhub/atmospheric_physics - fxtag = atmos_phys0_14_001_noresm_v0 + url = https://github.com/mvertens/atmospheric_physics + fxtag = 339524f fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/NorESMhub/atmospheric_physics @@ -76,8 +76,8 @@ [submodule "oslo_aero"] path = src/chemistry/oslo_aero - url = https://github.com/NorESMhub/OSLO_AERO - fxtag = oslo_aero_3_0a007 + url = https://github.com/mvertens/OSLO_AERO + fxtag = 55cd786 fxrequired = AlwaysRequired fxDONOTUSEurl = https://github.com/NorESMhub/OSLO_AERO.git diff --git a/bld/build-namelist b/bld/build-namelist index e128e5079c..c9bad05851 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -455,8 +455,7 @@ if ($print>=2) { # Composition of air add_default($nl, 'dry_air_species'); add_default($nl, 'water_species_in_air'); -# Enthalpy flux -add_default($nl, 'compute_enthalpy_flux'); + # Spectral Element dycore my $dyn = $cfg->get('dyn'); @@ -4553,6 +4552,7 @@ add_default($nl, 'history_gas'); add_default($nl, 'history_aerosol_forcing'); add_default($nl, 'history_aerosol_radiation'); add_default($nl, 'history_aerosol_debug_output'); +add_default($nl, 'history_enthalpy_flux'); # The history output for the AMWG variability diagnostics assumes that auxilliary history # files h1, h2, and h3 contain daily, 6-hrly, and 3-hrly output respectively. If this output diff --git a/bld/namelist_files/namelist_defaults_cam.xml b/bld/namelist_files/namelist_defaults_cam.xml index 1d8c9d3cee..39c189f11e 100644 --- a/bld/namelist_files/namelist_defaults_cam.xml +++ b/bld/namelist_files/namelist_defaults_cam.xml @@ -2787,6 +2787,14 @@ See https://github.com/NorESMhub/noresm3_dev_simulations/discussions/78 .false. 0.5 + + 1.0 + .false. + .false. + 2e-4 + 0.1 + 6.e2 + 1.0 @@ -2993,6 +3001,7 @@ See https://github.com/NorESMhub/noresm3_dev_simulations/discussions/78 .false. .false. .false. + .false. -1 -1 @@ -3012,9 +3021,6 @@ See https://github.com/NorESMhub/noresm3_dev_simulations/discussions/78 'Q','CLDLIQ','CLDICE','RAINQM','SNOWQM' 'Q','CLDLIQ','CLDICE','RAINQM','SNOWQM','GRAUQM' - - .false. - diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index a3ce04d21c..63ad1021c1 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -2275,7 +2275,7 @@ Default: 0,-24,-24,-24,-24,-24,-24,-24,-24,-24 group="cam_history_nl" valid_values=""> If interpolate_output(k) = .true., then the k'th history file will be interpolated to a lat/lon grid before output. -Default: .false.,.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false.,.false. +Default: .false. + +Output enthalpy flux diagnostics if CAM computes enthalpy fluxes +Default: .false. + + True when model is configured to use an offline driver. @@ -8282,15 +8288,6 @@ Switch to apply lunar tidal tendencies to neutral winds. Default: FALSE - - - -Enthalpy flux terms explicitly computed and added in the atmosphere and -passed to an active ocean component. -Default: FALSE - - + + + + + + + + + diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/shell_commands b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/shell_commands new file mode 100644 index 0000000000..eb40ad83e0 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/shell_commands @@ -0,0 +1,2 @@ +./xmlchange ROF_NCPL=\$ATM_NCPL +./xmlchange GLC_NCPL=\$ATM_NCPL diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cam b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cam new file mode 100644 index 0000000000..77424c653b --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cam @@ -0,0 +1,5 @@ +mfilt=1,1,1,1,1,1 +ndens=1,1,1,1,1,1 +nhtfrq=9,9,9,9,9,9 +write_nstep0=.true. +inithist='ENDOFRUN' diff --git a/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cpl b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cpl new file mode 100644 index 0000000000..58bfdcc086 --- /dev/null +++ b/cime_config/testdefs/testmods_dirs/cam/outfrq9s_enthalpy_from_cam/user_nl_cpl @@ -0,0 +1 @@ +component_computes_enthalpy_flux = 'atm' \ No newline at end of file diff --git a/src/atmos_phys b/src/atmos_phys index 3b9898008c..339524f0d7 160000 --- a/src/atmos_phys +++ b/src/atmos_phys @@ -1 +1 @@ -Subproject commit 3b9898008cc3bf580403c674a92ed48509de7d0e +Subproject commit 339524f0d7bae50669de9cc651fd5bbc2252ae87 diff --git a/src/chemistry/oslo_aero b/src/chemistry/oslo_aero index 0c99adb327..55cd786ca0 160000 --- a/src/chemistry/oslo_aero +++ b/src/chemistry/oslo_aero @@ -1 +1 @@ -Subproject commit 0c99adb32794b3cec92766b4efbce9f0241dde44 +Subproject commit 55cd786ca0c6b0f3e2603b8d9b2464626ff5f743 diff --git a/src/control/cam_comp.F90 b/src/control/cam_comp.F90 index c07d342923..c55a0d6a70 100644 --- a/src/control/cam_comp.F90 +++ b/src/control/cam_comp.F90 @@ -28,6 +28,7 @@ module cam_comp use perf_mod use cam_logfile, only: iulog use cam_abortutils, only: endrun +use air_composition, only: air_composition_register implicit none private @@ -56,7 +57,8 @@ module cam_comp subroutine cam_init( & caseid, ctitle, model_doi_url, & initial_run_in, restart_run_in, branch_run_in, post_assim_in, & - calendar, brnch_retain_casename, aqua_planet, dms_from_ocn, & + calendar, brnch_retain_casename, aqua_planet, dms_from_ocn, & + compute_enthalpy_flux, & single_column, scmlat, scmlon, & eccen, obliqr, lambm0, mvelpp, & perpetual_run, perpetual_ymd, & @@ -103,6 +105,7 @@ subroutine cam_init( & logical, intent(in) :: single_column logical, intent(in) :: dms_from_ocn + logical, intent(in) :: compute_enthalpy_flux real(r8), intent(in) :: scmlat real(r8), intent(in) :: scmlon @@ -168,14 +171,19 @@ subroutine cam_init( & ! Register zonal average grid for phys TEM diagnostics call phys_grid_ctem_reg() + ! Need to call this before phys_register - sets module variable + ! compute_enthalpy_flux in air_composition_register + call air_composition_register(compute_enthalpy_flux) + ! Register advected tracers and physics buffer fields - call phys_register () + call phys_register() ! Initialize ghg surface values before default initial distributions ! are set in dyn_init call chem_surfvals_init() call air_composition_init() + ! initialize ionosphere call ionosphere_init() diff --git a/src/control/camsrfexch.F90 b/src/control/camsrfexch.F90 index fbbb7a20c2..6380d0a79e 100644 --- a/src/control/camsrfexch.F90 +++ b/src/control/camsrfexch.F90 @@ -16,6 +16,7 @@ module camsrfexch active_Fall_flxdst1, active_Fall_flxvoc, active_Fall_flxfire use cam_control_mod, only: aqua_planet, simple_phys + implicit none private @@ -25,7 +26,7 @@ module camsrfexch public atm2hub_deallocate public hub2atm_deallocate public cam_export - + public get_prec_vars ! Public data types public cam_out_t ! Data from atmosphere public cam_in_t ! Merged surface data @@ -52,6 +53,11 @@ module camsrfexch real(r8) :: precsl(pcols) ! real(r8) :: precc(pcols) ! real(r8) :: precl(pcols) ! + real(r8) :: hrain(pcols) ! material enth. flx for liquid precip + real(r8) :: hsnow(pcols) ! material enth. flx for frozen precip + real(r8) :: hevap(pcols) ! material enth. flx for evaporation + real(r8) :: hmat (pcols) ! material enth. flx at surface, total + real(r8) :: hlat (pcols) ! variable latent heat component of hmat real(r8) :: soll(pcols) ! real(r8) :: sols(pcols) ! real(r8) :: solld(pcols) ! @@ -115,6 +121,8 @@ module camsrfexch real(r8) :: icefrac(pcols) ! sea-ice areal fraction real(r8) :: ocnfrac(pcols) ! ocean areal fraction real(r8) :: cflx(pcols,pcnst) ! constituent flux (emissions) + real(r8) :: evap_ocn(pcols) ! evaporation over ocean + real(r8) :: hrof(pcols) ! enthalpy from river runoff real(r8) :: ustar(pcols) ! atm/ocn saved version of ustar real(r8) :: re(pcols) ! atm/ocn saved version of re real(r8) :: ssq(pcols) ! atm/ocn saved version of ssq @@ -250,6 +258,7 @@ subroutine hub2atm_alloc( cam_in ) cam_in(c)%meganflx(:,:) = 0.0_r8 cam_in(c)%cflx (:,:) = 0._r8 + cam_in(c)%evap_ocn (:) = 0._r8 cam_in(c)%ustar (:) = 0._r8 cam_in(c)%re (:) = 0._r8 cam_in(c)%ssq (:) = 0._r8 @@ -329,20 +338,17 @@ subroutine atm2hub_alloc( cam_out ) cam_out(c)%dstwet3(:) = 0._r8 cam_out(c)%dstdry4(:) = 0._r8 cam_out(c)%dstwet4(:) = 0._r8 + cam_out(c)%hevap(:) = 0._r8 nullify(cam_out(c)%nhx_nitrogen_flx) nullify(cam_out(c)%noy_nitrogen_flx) - - if (.not. (simple_phys .or. aqua_planet)) then - - allocate (cam_out(c)%nhx_nitrogen_flx(pcols), stat=ierror) - if ( ierror /= 0 ) call endrun(sub//': allocation error nhx_nitrogen_flx') - cam_out(c)%nhx_nitrogen_flx(:) = 0._r8 - - allocate (cam_out(c)%noy_nitrogen_flx(pcols), stat=ierror) - if ( ierror /= 0 ) call endrun(sub//': allocation error noy_nitrogen_flx') - cam_out(c)%noy_nitrogen_flx(:) = 0._r8 - + if (.not.(simple_phys .or. aqua_planet)) then + allocate (cam_out(c)%nhx_nitrogen_flx(pcols), stat=ierror) + if ( ierror /= 0 ) call endrun(sub//': allocation error nhx_nitrogen_flx') + cam_out(c)%nhx_nitrogen_flx(:) = 0._r8 + allocate (cam_out(c)%noy_nitrogen_flx(pcols), stat=ierror) + if ( ierror /= 0 ) call endrun(sub//': allocation error noy_nitrogen_flx') + cam_out(c)%noy_nitrogen_flx(:) = 0._r8 endif end do @@ -410,7 +416,7 @@ end subroutine hub2atm_deallocate !====================================================================== -subroutine cam_export(state,cam_out,pbuf) +subroutine cam_export(state, cam_out, pbuf, cam_in) ! Transfer atmospheric fields into necessary surface data structures @@ -419,18 +425,22 @@ subroutine cam_export(state,cam_out,pbuf) use cam_history, only: outfld use chem_surfvals, only: chem_surfvals_get use co2_cycle, only: co2_transport, c_i - use physconst, only: rair, mwdry, mwco2, gravit, mwo3 + use physconst, only: rair, mwdry, mwco2, gravit, mwo3, cpliq, cpice, cpwv, tmelt use constituents, only: pcnst - use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc + use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_set_field use rad_constituents, only: rad_cnst_get_gas use cam_control_mod, only: simple_phys - + use air_composition, only: t00a, t00o, h00a, h00o + use air_composition, only: hliq_idx, hice_idx, fliq_idx, fice_idx + use air_composition, only: compute_enthalpy_flux, num_enthalpy_vars + use cam_history, only: outfld !xxx debug implicit none ! Input arguments - type(physics_state), intent(in) :: state - type (cam_out_t), intent(inout) :: cam_out - type(physics_buffer_desc), pointer :: pbuf(:) + type(physics_state), intent(in) :: state + type (cam_out_t), intent(inout) :: cam_out + type(physics_buffer_desc), pointer :: pbuf(:) + type (cam_in_t ), optional, intent(in) :: cam_in ! Local variables @@ -439,23 +449,23 @@ subroutine cam_export(state,cam_out,pbuf) integer :: lchnk ! Chunk index integer :: ncol integer :: psl_idx - integer :: prec_dp_idx, snow_dp_idx, prec_sh_idx, snow_sh_idx - integer :: prec_sed_idx,snow_sed_idx,prec_pcw_idx,snow_pcw_idx integer :: srf_ozone_idx, lightning_idx + integer :: enthalpy_prec_bc_idx, enthalpy_prec_ac_idx, enthalpy_evop_idx + real(r8):: ubot, vbot real(r8), pointer :: psl(:) - real(r8), pointer :: prec_dp(:) ! total precipitation from ZM convection - real(r8), pointer :: snow_dp(:) ! snow from ZM convection - real(r8), pointer :: prec_sh(:) ! total precipitation from Hack convection - real(r8), pointer :: snow_sh(:) ! snow from Hack convection - real(r8), pointer :: prec_sed(:) ! total precipitation from ZM convection - real(r8), pointer :: snow_sed(:) ! snow from ZM convection - real(r8), pointer :: prec_pcw(:) ! total precipitation from Hack convection - real(r8), pointer :: snow_pcw(:) ! snow from Hack convection real(r8), pointer :: o3_ptr(:,:), srf_o3_ptr(:) real(r8), pointer :: lightning_ptr(:) + + ! enthalpy variables (if applicable) + real(r8), dimension(:,:), pointer :: enthalpy_prec_ac + real(r8), dimension(:) , pointer :: hevap_ocn + real(r8), dimension(pcols) :: fliq_tot, fice_tot + real(r8), dimension(pcols,num_enthalpy_vars) :: enthalpy_prec_bc + + character(len=*), parameter :: sub = 'cam_export' !----------------------------------------------------------------------- lchnk = state%lchnk @@ -464,42 +474,73 @@ subroutine cam_export(state,cam_out,pbuf) psl_idx = pbuf_get_index('PSL') call pbuf_get_field(pbuf, psl_idx, psl) - prec_dp_idx = pbuf_get_index('PREC_DP', errcode=i) - snow_dp_idx = pbuf_get_index('SNOW_DP', errcode=i) - prec_sh_idx = pbuf_get_index('PREC_SH', errcode=i) - snow_sh_idx = pbuf_get_index('SNOW_SH', errcode=i) - prec_sed_idx = pbuf_get_index('PREC_SED', errcode=i) - snow_sed_idx = pbuf_get_index('SNOW_SED', errcode=i) - prec_pcw_idx = pbuf_get_index('PREC_PCW', errcode=i) - snow_pcw_idx = pbuf_get_index('SNOW_PCW', errcode=i) - srf_ozone_idx = pbuf_get_index('SRFOZONE', errcode=i) - lightning_idx = pbuf_get_index('LGHT_FLASH_FREQ', errcode=i) + if (compute_enthalpy_flux) then + enthalpy_prec_bc_idx = pbuf_get_index('ENTHALPY_PREC_BC', errcode=i) + enthalpy_prec_ac_idx = pbuf_get_index('ENTHALPY_PREC_AC', errcode=i) + if (enthalpy_prec_bc_idx==0.or.enthalpy_prec_ac_idx==0) then + call endrun(sub//": pbufs for enthalpy flux not allocated") + end if + call pbuf_get_field(pbuf, enthalpy_prec_ac_idx, enthalpy_prec_ac) + + !------------------------------------------------------------------ + ! + ! compute precipitation fluxes and set associated physics buffers + ! + !------------------------------------------------------------------ + call get_prec_vars(ncol,pbuf,fliq=fliq_tot,fice=fice_tot,& + precc_out=cam_out%precc,precl_out=cam_out%precl,& + precsc_out=cam_out%precsc,precsl_out=cam_out%precsl) + + ! fliq_tot holds liquid precipitation from tphysbc and + ! tphysac from previous physics time-step: back out fliq_bc + ! Idem for ice + enthalpy_prec_bc(:ncol,fice_idx) = fice_tot(:ncol) -enthalpy_prec_ac(:ncol,fice_idx) ! out of atm + enthalpy_prec_bc(:ncol,fliq_idx) = fliq_tot(:ncol) -enthalpy_prec_ac(:ncol,fliq_idx) ! out of atm + + ! compute precipitation enthalpy fluxes from tphysbc + ! correct for reference T of latent heats (liquid reference state), and use tbot (=T(pver), updated later below) + enthalpy_prec_bc(:ncol,hice_idx) = -enthalpy_prec_bc(:ncol,fice_idx)*(cpice*(state%T(:ncol,pver)-t00a)+(cpliq*t00a+h00a)) + enthalpy_prec_bc(:ncol,hliq_idx) = -enthalpy_prec_bc(:ncol,fliq_idx)*(cpliq*(state%T(:ncol,pver)-t00a)+(cpliq*t00a+h00a)) + + ! export all prec_bc to pbuf + call pbuf_set_field(pbuf, enthalpy_prec_bc_idx, enthalpy_prec_bc) + + ! Compute enthalpy fluxes for the coupler: + cam_out%hsnow(:ncol) = enthalpy_prec_bc(:ncol,hice_idx)+enthalpy_prec_ac(:ncol,hice_idx) ! into atm + cam_out%hrain(:ncol) = enthalpy_prec_bc(:ncol,hliq_idx)+enthalpy_prec_ac(:ncol,hliq_idx) ! into atm + + ! change enthalpy flux to sign convention of ocean model and change zero points + cam_out%hsnow(:ncol) = -cam_out%hsnow(:ncol) + fice_tot(:ncol)*((h00o-h00a)+(cpliq-cpice)*(t00o-t00a)) ! into ocn; fice_tot is out of atm + cam_out%hrain(:ncol) = -cam_out%hrain(:ncol) + fliq_tot(:ncol)* (h00o-h00a)! +0. ! into ocn; fliq_tot is out of atm + + if (present(cam_in)) then + ! hevap is one time-step old, consistently with rest of enthalpy_prec_ac + enthalpy_evop_idx = pbuf_get_index('ENTHALPY_EVOP', errcode=i) + if (enthalpy_evop_idx==0) then + call endrun(sub//": pbuf for enthalpy evop not allocated") + end if + call pbuf_get_field(pbuf, enthalpy_evop_idx, hevap_ocn) + cam_out%hevap(:ncol) = -hevap_ocn(:ncol) - cam_in%evap_ocn(:ncol)*((h00o-h00a)+(cpliq-cpwv )*(t00o-t00a)) ! into ocn; cflux is into atm + end if + + cam_out%hmat(:ncol) = cam_out%hsnow(:ncol) + cam_out%hrain(:ncol) + cam_out%hevap(:ncol) ! this is into ocean + ! variable latent heat component + ! N.B.: approximate due to difference between ts and tbot, also note lagged SST + cam_out%hlat(:ncol) = cam_in%evap_ocn(:ncol)*(cpliq-cpwv )*(cam_in%sst(:ncol)-t00a) & + -fice_tot (:ncol)*(cpliq-cpice)*(cam_in%sst(:ncol)-t00a) + else + + call get_prec_vars(ncol,pbuf,& + precc_out=cam_out%precc,precl_out=cam_out%precl,& + precsc_out=cam_out%precsc,precsl_out=cam_out%precsl) + cam_out%hmat(:ncol) = 0._r8 + cam_out%hlat(:ncol) = 0._r8 - if (prec_dp_idx > 0) then - call pbuf_get_field(pbuf, prec_dp_idx, prec_dp) - end if - if (snow_dp_idx > 0) then - call pbuf_get_field(pbuf, snow_dp_idx, snow_dp) - end if - if (prec_sh_idx > 0) then - call pbuf_get_field(pbuf, prec_sh_idx, prec_sh) - end if - if (snow_sh_idx > 0) then - call pbuf_get_field(pbuf, snow_sh_idx, snow_sh) - end if - if (prec_sed_idx > 0) then - call pbuf_get_field(pbuf, prec_sed_idx, prec_sed) - end if - if (snow_sed_idx > 0) then - call pbuf_get_field(pbuf, snow_sed_idx, snow_sed) - end if - if (prec_pcw_idx > 0) then - call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw) - end if - if (snow_pcw_idx > 0) then - call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw) end if + srf_ozone_idx = pbuf_get_index('SRFOZONE', errcode=i) + lightning_idx = pbuf_get_index('LGHT_FLASH_FREQ', errcode=i) + do i=1,ncol cam_out%tbot(i) = state%t(i,pver) cam_out%thbot(i) = state%t(i,pver) * state%exner(i,pver) @@ -510,7 +551,6 @@ subroutine cam_export(state,cam_out,pbuf) cam_out%pbot(i) = state%pmid(i,pver) cam_out%psl(i) = psl(i) cam_out%rho(i) = cam_out%pbot(i)/(rair*cam_out%tbot(i)) - ! Direction of bottom level wind ubot = state%u(i,pver) vbot = state%v(i,pver) @@ -547,51 +587,120 @@ subroutine cam_export(state,cam_out,pbuf) call pbuf_get_field(pbuf, lightning_idx, lightning_ptr) cam_out%lightning_flash_freq(:ncol) = lightning_ptr(:ncol) end if - - ! - ! Precipation and snow rates from shallow convection, deep convection and stratiform processes. - ! Compute total convective and stratiform precipitation and snow rates - ! - do i=1,ncol - cam_out%precc (i) = 0._r8 - cam_out%precl (i) = 0._r8 - cam_out%precsc(i) = 0._r8 - cam_out%precsl(i) = 0._r8 - if (prec_dp_idx > 0) then - cam_out%precc (i) = cam_out%precc (i) + prec_dp(i) - end if - if (prec_sh_idx > 0) then - cam_out%precc (i) = cam_out%precc (i) + prec_sh(i) - end if - if (prec_sed_idx > 0) then - cam_out%precl (i) = cam_out%precl (i) + prec_sed(i) - end if - if (prec_pcw_idx > 0) then - cam_out%precl (i) = cam_out%precl (i) + prec_pcw(i) - end if - if (snow_dp_idx > 0) then - cam_out%precsc(i) = cam_out%precsc(i) + snow_dp(i) - end if - if (snow_sh_idx > 0) then - cam_out%precsc(i) = cam_out%precsc(i) + snow_sh(i) - end if - if (snow_sed_idx > 0) then - cam_out%precsl(i) = cam_out%precsl(i) + snow_sed(i) - end if - if (snow_pcw_idx > 0) then - cam_out%precsl(i) = cam_out%precsl(i) + snow_pcw(i) - end if - - ! jrm These checks should not be necessary if they exist in the parameterizations - if (cam_out%precc(i) .lt.0._r8) cam_out%precc(i)=0._r8 - if (cam_out%precl(i) .lt.0._r8) cam_out%precl(i)=0._r8 - if (cam_out%precsc(i).lt.0._r8) cam_out%precsc(i)=0._r8 - if (cam_out%precsl(i).lt.0._r8) cam_out%precsl(i)=0._r8 - if (cam_out%precsc(i).gt.cam_out%precc(i)) cam_out%precsc(i)=cam_out%precc(i) - if (cam_out%precsl(i).gt.cam_out%precl(i)) cam_out%precsl(i)=cam_out%precl(i) - - end do - end subroutine cam_export +! +! Precipation and snow rates from shallow convection, deep convection and stratiform processes. +! Compute total convective and stratiform precipitation and snow rates +! +subroutine get_prec_vars(ncol,pbuf,fliq,fice, precc_out,precl_out,precsc_out,precsl_out) + use ppgrid, only: pcols + use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc + + integer, intent(in) :: ncol + type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), dimension(pcols) , optional, intent(out):: fliq!rain flux (out of atm) in SI units + real(r8), dimension(pcols) , optional, intent(out):: fice!snow flux (out of atm) in SI units + + real(r8), dimension(pcols), optional, intent(out):: precc_out !total precipitation from convection + real(r8), dimension(pcols), optional, intent(out):: precl_out !total large scale precipitation + real(r8), dimension(pcols), optional, intent(out):: precsc_out!frozen precipitation from convection + real(r8), dimension(pcols), optional, intent(out):: precsl_out!frozen large scale precipitation + + integer :: i + + real(r8), pointer :: prec_dp(:) !total precipitation from from deep convection + real(r8), pointer :: snow_dp(:) !frozen precipitation from deep convection + real(r8), pointer :: prec_sh(:) !total precipitation from shallow convection + real(r8), pointer :: snow_sh(:) !frozen precipitation from from shallow convection + real(r8), pointer :: prec_sed(:) !total precipitation from cloud sedimentation + real(r8), pointer :: snow_sed(:) !frozen precipitation from sedimentation + real(r8), pointer :: prec_pcw(:) !total precipitation from from microphysics + real(r8), pointer :: snow_pcw(:) !frozen precipitation from from microphysics + + real(r8), dimension(pcols):: precc, precl, precsc, precsl + integer :: prec_dp_idx, snow_dp_idx, prec_sh_idx, snow_sh_idx + integer :: prec_sed_idx,snow_sed_idx,prec_pcw_idx,snow_pcw_idx + ! + ! get fields from pbuf + ! + prec_dp_idx = pbuf_get_index('PREC_DP', errcode=i) + snow_dp_idx = pbuf_get_index('SNOW_DP', errcode=i) + prec_sh_idx = pbuf_get_index('PREC_SH', errcode=i) + snow_sh_idx = pbuf_get_index('SNOW_SH', errcode=i) + prec_sed_idx = pbuf_get_index('PREC_SED', errcode=i) + snow_sed_idx = pbuf_get_index('SNOW_SED', errcode=i) + prec_pcw_idx = pbuf_get_index('PREC_PCW', errcode=i) + snow_pcw_idx = pbuf_get_index('SNOW_PCW', errcode=i) + + if (prec_dp_idx > 0) then + call pbuf_get_field(pbuf, prec_dp_idx, prec_dp) + end if + if (snow_dp_idx > 0) then + call pbuf_get_field(pbuf, snow_dp_idx, snow_dp) + end if + if (prec_sh_idx > 0) then + call pbuf_get_field(pbuf, prec_sh_idx, prec_sh) + end if + if (snow_sh_idx > 0) then + call pbuf_get_field(pbuf, snow_sh_idx, snow_sh) + end if + if (prec_sed_idx > 0) then + call pbuf_get_field(pbuf, prec_sed_idx, prec_sed) + end if + if (snow_sed_idx > 0) then + call pbuf_get_field(pbuf, snow_sed_idx, snow_sed) + end if + if (prec_pcw_idx > 0) then + call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw) + end if + if (snow_pcw_idx > 0) then + call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw) + end if + + precc = 0._r8 + precl = 0._r8 + precsc = 0._r8 + precsl = 0._r8 + if (prec_dp_idx > 0) then + precc(:ncol) = precc(:ncol) + prec_dp(:ncol) + end if + if (prec_sh_idx > 0) then + precc(:ncol) = precc(:ncol) + prec_sh(:ncol) + end if + if (prec_sed_idx > 0) then + precl(:ncol) = precl(1:ncol) + prec_sed(:ncol) + end if + if (prec_pcw_idx > 0) then + precl(:ncol) = precl(1:ncol) + prec_pcw(:ncol) + end if + if (snow_dp_idx > 0) then + precsc(:ncol) = precsc(:ncol) + snow_dp(:ncol) + end if + if (snow_sh_idx > 0) then + precsc(:ncol) = precsc(:ncol) + snow_sh(:ncol) + end if + if (snow_sed_idx > 0) then + precsl(:ncol) = precsl(:ncol) + snow_sed(:ncol) + end if + if (snow_pcw_idx > 0) then + precsl(:ncol)= precsl(:ncol) + snow_pcw(:ncol) + end if + + do i=1,ncol + precc(i) = MAX(precc(i), 0.0_r8) + precl(i) = MAX(precl(i), 0.0_r8) + precsc(i) = MAX(precsc(i),0.0_r8) + precsl(i) = MAX(precsl(i),0.0_r8) + if (precsc(i).gt.precc(i)) precsc(i)=precc(i) + if (precsl(i).gt.precl(i)) precsl(i)=precl(i) + end do + if (present(precc_out )) precc_out (:ncol) = precc (:ncol) + if (present(precl_out )) precl_out (:ncol) = precl (:ncol) + if (present(precsc_out)) precsc_out(:ncol) = precsc(:ncol) + if (present(precsl_out)) precsl_out(:ncol) = precsl(:ncol) + + if (present(fice)) fice(:ncol) = 1000.0_r8*(precsc(:ncol)+precsl(:ncol)) !snow flux + if (present(fliq)) fliq(:ncol) = 1000.0_r8*(precc (:ncol)-precsc(:ncol)+precl(:ncol)-precsl(:ncol))!rain flux + end subroutine get_prec_vars end module camsrfexch diff --git a/src/cpl/nuopc/atm_comp_nuopc.F90 b/src/cpl/nuopc/atm_comp_nuopc.F90 index 10150e4dc2..42a459a548 100644 --- a/src/cpl/nuopc/atm_comp_nuopc.F90 +++ b/src/cpl/nuopc/atm_comp_nuopc.F90 @@ -127,6 +127,8 @@ module atm_comp_nuopc character(len=*) , parameter :: orb_variable_year = 'variable_year' character(len=*) , parameter :: orb_fixed_parameters = 'fixed_parameters' + logical :: compute_enthalpy_flux ! If true, CAM computes enthalpy flux + real(R8) , parameter :: grid_tol = 1.e-2_r8 ! tolerance for calculated lat/lon vs read in type(ESMF_Mesh) :: model_mesh ! model_mesh @@ -301,6 +303,29 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) call shr_sys_abort(subname//'Need to set attribute ScalarFieldIdxNextSwCday') endif + call NUOPC_CompAttributeGet(gcomp, name='atm_computes_enthalpy_flux', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + compute_enthalpy_flux = .false. + if (isPresent .and. isSet) then + if (trim(cvalue) == 'atm') then + compute_enthalpy_flux = .true. + end if + end if + + call NUOPC_CompAttributeGet(gcomp, name='component_computes_enthalpy_flux', value=cvalue, & + isPresent=isPresent, isSet=isSet, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (isPresent .and. isSet) then + if (trim(cvalue) == 'atm') then + compute_enthalpy_flux = .true. + else + compute_enthalpy_flux = .false. + end if + else + compute_enthalpy_flux = .false. + end if + ! read mediator fields namelists call read_surface_fields_namelists() @@ -309,7 +334,7 @@ subroutine InitializeAdvertise(gcomp, importState, exportState, clock, rc) if (isPresent .and. isSet) then read (cvalue,*) mediator_present if (mediator_present) then - call advertise_fields(gcomp, flds_scalar_name, rc) + call advertise_fields(gcomp, flds_scalar_name, compute_enthalpy_flux, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return end if else @@ -650,6 +675,7 @@ subroutine InitializeRealize(gcomp, importState, exportState, clock, rc) branch_run_in=branch_run, post_assim_in=dart_mode, & calendar=calendar, brnch_retain_casename=brnch_retain_casename, & aqua_planet=aqua_planet, dms_from_ocn=dms_from_ocn, & + compute_enthalpy_flux=compute_enthalpy_flux, & single_column=single_column, scmlat=scol_lat, scmlon=scol_lon, & eccen=eccen, obliqr=obliqr, lambm0=lambm0, mvelpp=mvelpp, & perpetual_run=perpetual_run, perpetual_ymd=perpetual_ymd, & diff --git a/src/cpl/nuopc/atm_import_export.F90 b/src/cpl/nuopc/atm_import_export.F90 index 907c4d2f1e..011be4e9ef 100644 --- a/src/cpl/nuopc/atm_import_export.F90 +++ b/src/cpl/nuopc/atm_import_export.F90 @@ -13,7 +13,7 @@ module atm_import_export use shr_mpi_mod , only : shr_mpi_min, shr_mpi_max use nuopc_shr_methods , only : chkerr use cam_logfile , only : iulog - use cam_history , only: outfld + use cam_history , only : outfld use spmd_utils , only : masterproc, mpicom use srf_field_check , only : set_active_Sl_ram1 use srf_field_check , only : set_active_Sl_fv @@ -62,9 +62,12 @@ module atm_import_export integer :: megan_nflds = -huge(1) ! number of MEGAN voc fields from lnd-> atm integer :: emis_nflds = -huge(1) ! number of fire emission fields from lnd-> atm logical :: atm_provides_lightning = .false. ! cld to grnd lightning flash freq (min-1) + logical, public :: dms_from_ocn = .false. ! dms is obtained from ocean as atm import data logical, public :: brf_from_ocn = .false. ! brf is obtained from ocean as atm import data logical, public :: n2o_from_ocn = .false. ! n2o is obtained from ocean as atm import data logical, public :: nh3_from_ocn = .false. ! nh3 is obtained from ocean as atm import data + logical, protected :: compute_enthalpy_flux = .false. + character(*),parameter :: F01 = "('(cam_import_export) ',a,i8,2x,i8,2x,d21.14)" character(*),parameter :: F02 = "('(cam_import_export) ',a,i8,2x,i8,2x,i8,2x,d21.14)" character(*),parameter :: u_FILE_u = __FILE__ @@ -98,11 +101,12 @@ end subroutine read_surface_fields_namelists !----------------------------------------------------------- ! advertise fields !----------------------------------------------------------- - subroutine advertise_fields(gcomp, flds_scalar_name, rc) + subroutine advertise_fields(gcomp, flds_scalar_name, compute_enthalpy_flux_in, rc) ! input/output variables type(ESMF_GridComp) :: gcomp character(len=*) , intent(in) :: flds_scalar_name + logical , intent(in) :: compute_enthalpy_flux_in integer , intent(out) :: rc ! local variables @@ -115,7 +119,6 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) logical :: flds_co2b ! use case logical :: flds_co2c ! use case character(len=128) :: fldname - logical :: dms_from_ocn ! dms is obtained from ocean as atm import data logical :: ispresent logical :: isset character(len=*), parameter :: subname='(atm_import_export:advertise_fields): ' @@ -130,6 +133,12 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) ! determine necessary toggles for below !-------------------------------- + ! Set module variable + compute_enthalpy_flux = compute_enthalpy_flux_in + if (masterproc) then + write(iulog,'(2a,l)') trim(subname), 'compute_enthalpy_flux = ',compute_enthalpy_flux + end if + call NUOPC_CompAttributeGet(gcomp, name='flds_co2a', value=cvalue, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return read(cvalue,*) flds_co2a @@ -219,6 +228,10 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_rainl' ) call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_snowc' ) call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_snowl' ) + if (compute_enthalpy_flux) then + call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_hmat' ) ! enthalpy flux computed by cam + call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_hlat' ) ! var.lat.ht.part + end if call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_lwdn' ) call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_swndr' ) call fldlist_add(fldsFrAtm_num, fldsFrAtm, 'Faxa_swvdr' ) @@ -301,6 +314,10 @@ subroutine advertise_fields(gcomp, flds_scalar_name, rc) call fldlist_add(fldsToAtm_num, fldsToAtm, 'Faxx_sen' ) call fldlist_add(fldsToAtm_num, fldsToAtm, 'Faxx_lwup' ) call fldlist_add(fldsToAtm_num, fldsToAtm, 'Faxx_evap' ) + if (compute_enthalpy_flux) then + call fldlist_add(fldsToAtm_num, fldsToAtm, 'Faox_evap' ) + call fldlist_add(fldsToAtm_num, fldsToAtm, 'Faxx_hrof' ) + end if ! dust fluxes from land (4 sizes) call fldlist_add(fldsToAtm_num, fldsToAtm, 'Fall_flxdst', ungridded_lbound=1, ungridded_ubound=4) @@ -583,6 +600,8 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) real(r8), pointer :: fldptr_tauy(:) real(r8), pointer :: fldptr_sen(:) real(r8), pointer :: fldptr_evap(:) + real(r8), pointer :: fldptr_evop(:) + real(r8), pointer :: fldptr_hrof(:) logical, save :: first_time = .true. character(len=*), parameter :: subname='(atm_import_export:import_fields)' !--------------------------------------------------------------------------- @@ -611,6 +630,31 @@ subroutine import_fields( gcomp, cam_in, restart_init, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getfldptr(importState, 'Faxx_evap', fldptr=fldptr_evap, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + + ! ***NOTE:*** if atm_compute_enthalpy_flux is .false. and if in + ! CMEPS med_computes_enthalpy_flux is .true., then the mediator + ! will compute it if the ocean requests it and add a correction + ! to the sensible heat sent to cam. + ! This is the case if CAM is coupled to MOM6. + ! However, it is not the case if CAM is coupled to BLOM. + + if (compute_enthalpy_flux) then + ! ocean-point hevap + call state_getfldptr(importState, 'Faox_evap', fldptr=fldptr_evop, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + ! enthalpy of runoff + call state_getfldptr(importState, 'Faxx_hrof', fldptr=fldptr_hrof, rc=rc) + if (ChkErr(rc,__LINE__,u_FILE_u)) return + g = 1 + do c = begchunk,endchunk + do i = 1,get_ncols_p(c) + cam_in(c)%evap_ocn(i) = -fldptr_evop(g) * med2mod_areacor(g) + cam_in(c)%hrof(i) = -fldptr_hrof(g) * med2mod_areacor(g) + g = g + 1 + end do + end do + end if ! end of compute_enthalpy_flux + g = 1 do c = begchunk,endchunk do i = 1,get_ncols_p(c) @@ -1052,6 +1096,7 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) real(r8), pointer :: fldptr_soll(:) , fldptr_sols(:) real(r8), pointer :: fldptr_solld(:) , fldptr_solsd(:) real(r8), pointer :: fldptr_snowc(:) , fldptr_snowl(:) + real(r8), pointer :: fldptr_hmat (:) , fldptr_hlat (:) ! enthalpy flux computed by cam real(r8), pointer :: fldptr_rainc(:) , fldptr_rainl(:) real(r8), pointer :: fldptr_lwdn(:) , fldptr_swnet(:) real(r8), pointer :: fldptr_topo(:) , fldptr_zbot(:) @@ -1150,6 +1195,20 @@ subroutine export_fields( gcomp, model_mesh, model_clock, cam_out, rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return call state_getfldptr(exportState, 'Faxa_swvdf', fldptr=fldptr_solsd, rc=rc) if (ChkErr(rc,__LINE__,u_FILE_u)) return + if (compute_enthalpy_flux) then + call state_getfldptr(exportState, 'Faxa_hmat' , fldptr=fldptr_hmat , rc=rc) ! enthalpy + if (ChkErr(rc,__LINE__,u_FILE_u)) return + call state_getfldptr(exportState, 'Faxa_hlat' , fldptr=fldptr_hlat , rc=rc) ! variable latent heat part + if (ChkErr(rc,__LINE__,u_FILE_u)) return + g = 1 + do c = begchunk,endchunk + do i = 1,get_ncols_p(c) + fldptr_hmat (g) = cam_out(c)%hmat(i) * mod2med_areacor(g) ! enthalpy + fldptr_hlat (g) = cam_out(c)%hlat(i) * mod2med_areacor(g) ! variable latent heat part + g = g + 1 + end do + end do + end if g = 1 do c = begchunk,endchunk do i = 1,get_ncols_p(c) diff --git a/src/physics/cam/cam_diagnostics.F90 b/src/physics/cam/cam_diagnostics.F90 index 5f7e7d9a60..9376e729fe 100644 --- a/src/physics/cam/cam_diagnostics.F90 +++ b/src/physics/cam/cam_diagnostics.F90 @@ -17,7 +17,7 @@ module cam_diagnostics use constituents, only: pcnst, cnst_name, cnst_longname, cnst_cam_outfld use constituents, only: ptendnam, apcnst, bpcnst, cnst_get_ind use dycore, only: dycore_is -use phys_control, only: phys_getopts +use phys_control, only: phys_getopts, history_enthalpy_flux use wv_saturation, only: qsat, qsat_water, svp_ice_vect use time_manager, only: is_first_step @@ -191,10 +191,12 @@ subroutine diag_init_dry(pbuf2d) use cam_history, only: register_vector_field use tidal_diag, only: tidal_diag_init use cam_budget, only: cam_budget_em_snapshot, cam_budget_em_register, thermo_budget_history + use air_composition, only: compute_enthalpy_flux type(physics_buffer_desc), pointer, intent(in) :: pbuf2d(:,:) integer :: istage + ! outfld calls in diag_phys_writeout call addfld (cnst_name(1), (/ 'lev' /), 'A', 'kg/kg', cnst_longname(1)) call addfld ('NSTEP', horiz_only, 'A', 'timestep', 'Model timestep') @@ -224,6 +226,19 @@ subroutine diag_init_dry(pbuf2d) call addfld ('TFIX', horiz_only, 'A', 'K/s', 'T fixer (T equivalent of Energy correction)') call addfld ('TTEND_TOT', (/ 'lev' /), 'A', 'K/s', 'Total temperature tendency') + if (compute_enthalpy_flux) then + call addfld('EBREAK' , horiz_only, 'A','W/m2', & + 'Global-mean energy-nonconservation (W/m2)' ) + call addfld('PTTEND_DME', (/ 'lev' /), 'A', 'K/s ', & + 'T-tendency due to water fluxes (end of tphysac)' ) + call addfld('IETEND_DME', horiz_only, 'A','W/m2 ', & + 'Column enthalpy tendency due to water fluxes (end of tphysac)' ) + call addfld('EFLX ' , horiz_only, 'A','W/m2 ', & + 'Surface water material enthalpy flux (end of tphysac)' ) + call addfld('MFLX ' , horiz_only, 'A','W/m2 ', & + 'Mass flux due to dry mass adjustment / water changes (end of tphysac)') + end if + ! outfld calls in diag_phys_tend_writeout call addfld ('UTEND_TOT', (/ 'lev' /), 'A', 'm/s2', 'Total zonal wind tendency') call addfld ('VTEND_TOT', (/ 'lev' /), 'A', 'm/s2', 'Total meridional wind tendency') @@ -392,6 +407,44 @@ subroutine diag_init_dry(pbuf2d) call addfld( 'CPAIRV', (/ 'lev' /), 'I', 'J/K/kg', 'Variable specific heat cap air' ) call addfld( 'RAIRV', (/ 'lev' /), 'I', 'J/K/kg', 'Variable dry air gas constant' ) + if (compute_enthalpy_flux) then + if (history_enthalpy_flux) then + call addfld('enth_prec_ac_hice',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_ac_hliq',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_bc_hice',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_bc_hliq',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_ac_fice',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_ac_fliq',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_bc_fice',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_prec_bc_fliq',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fevap' ,horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_frain_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fsnow_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fwatr_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_frain_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fsnow_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fwatr_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_frain_tt_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fsnow_tt_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_fwatr_tt_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hevap_atm' ,horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hevap_ocn' ,horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hrain_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hsnow_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hwatr_bc_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hrain_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hsnow_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hwatr_ac_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hrain_tt_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hsnow_tt_err',horiz_only, 'A', 'W/m2', '' ) + call addfld('enth_hwatr_tt_err',horiz_only, 'A', 'W/m2', '' ) + endif + call addfld('te_tnd' , horiz_only, 'A', 'W/m2', 'Total column integrated energy tendency from CAM physics' ) + call addfld('dEdt_dme' , horiz_only, 'A', 'W/m2', 'Column integrated dEdt from water update') + call addfld('dEdt_physics' , horiz_only, 'A', 'W/m2', '' )!xxx diags will remove + call addfld('dEdt_efix_physics', horiz_only, 'A', 'W/m2', 'Column integrated physics energy fixer dEdt from enthalpy fixer' ) + endif + if (thermo_budget_history) then ! ! energy diagnostics addflds for vars_stage combinations plus e_m_snapshots @@ -1579,14 +1632,14 @@ subroutine diag_conv(state, ztodt, pbuf) type(physics_buffer_desc), pointer :: pbuf(:) ! convective precipitation variables - real(r8), pointer :: prec_dp(:) ! total precipitation from ZM convection - real(r8), pointer :: snow_dp(:) ! snow from ZM convection - real(r8), pointer :: prec_sh(:) ! total precipitation from Hack convection - real(r8), pointer :: snow_sh(:) ! snow from Hack convection - real(r8), pointer :: prec_sed(:) ! total precipitation from ZM convection - real(r8), pointer :: snow_sed(:) ! snow from ZM convection - real(r8), pointer :: prec_pcw(:) ! total precipitation from Hack convection - real(r8), pointer :: snow_pcw(:) ! snow from Hack convection + real(r8), pointer :: prec_dp(:) ! total precipitation from ZM convection + real(r8), pointer :: snow_dp(:) ! snow from ZM convection + real(r8), pointer :: prec_sh(:) ! total precipitation from Hack convection + real(r8), pointer :: snow_sh(:) ! snow from Hack convection + real(r8), pointer :: prec_sed(:) ! total precipitation from MG sedimentation + real(r8), pointer :: snow_sed(:) ! snow from MG sedimentation + real(r8), pointer :: prec_pcw(:) ! total precipitation from MG prog. cloud + real(r8), pointer :: snow_pcw(:) ! snow from MG prog. cloud ! Local variables: @@ -2028,6 +2081,7 @@ subroutine diag_phys_tend_writeout_dry(state, pbuf, tend, ztodt) use check_energy, only: check_energy_get_integrals use physconst, only: cpair + use air_composition, only: compute_enthalpy_flux ! Arguments @@ -2044,6 +2098,7 @@ subroutine diag_phys_tend_writeout_dry(state, pbuf, tend, ztodt) real(r8) :: ftem2(pcols) ! Temporary workspace for outfld variables real(r8) :: ftem3(pcols,pver) ! Temporary workspace for outfld variables real(r8) :: heat_glob ! global energy integral (FV only) + real(r8) :: tedif_glob ! energy flux from fixer ! CAM pointers to get variables from the physics buffer real(r8), pointer, dimension(:,:) :: t_ttend real(r8), pointer, dimension(:,:) :: t_utend @@ -2064,9 +2119,18 @@ subroutine diag_phys_tend_writeout_dry(state, pbuf, tend, ztodt) ! Total physics tendency for Temperature ! (remove global fixer tendency from total for FV and SE dycores) - call check_energy_get_integrals( heat_glob_out=heat_glob ) - ftem2(:ncol) = heat_glob/cpair - call outfld('TFIX', ftem2, pcols, lchnk ) + if (compute_enthalpy_flux) then + call check_energy_get_integrals(heat_glob_out=heat_glob, tedif_glob_out=tedif_glob) + ftem2(:ncol) = tedif_glob/ztodt + call outfld('EBREAK', ftem2, pcols, lchnk) + ftem2(:ncol) = heat_glob/cpair + call outfld('TFIX', ftem2, pcols, lchnk) + else + call check_energy_get_integrals( heat_glob_out=heat_glob ) + ftem2(:ncol) = heat_glob/cpair + call outfld('TFIX', ftem2, pcols, lchnk) + end if + ftem3(:ncol,:pver) = tend%dtdt(:ncol,:pver) - heat_glob/cpair call outfld('PTTEND',ftem3, pcols, lchnk ) ftem3(:ncol,:pver) = tend%dudt(:ncol,:pver) diff --git a/src/physics/cam/check_energy.F90 b/src/physics/cam/check_energy.F90 index d1d59e173f..358925cb3e 100644 --- a/src/physics/cam/check_energy.F90 +++ b/src/physics/cam/check_energy.F90 @@ -21,12 +21,12 @@ module check_energy !--------------------------------------------------------------------------------- use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: pcols, pver + use ppgrid, only: pcols, pver, begchunk, endchunk use spmd_utils, only: masterproc - use physconst, only: rga + use physconst, only: gravit, rga, latvap, latice, cpair, rair use air_composition, only: cpairv, cp_or_cv_dycore - use physics_types, only: physics_state + 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 use cam_logfile, only: iulog @@ -55,6 +55,9 @@ module check_energy public :: check_energy_cam_fix ! add heating rate required for global mean total energy conservation + ! This routine adjusts enthalpy if compute_enthalpy_flux = .true. + public :: enthalpy_adjustment + ! Private module data logical :: print_energy_errors = .false. @@ -67,6 +70,7 @@ module check_energy real(r8) :: heat_glob ! global mean heating rate ! Physics buffer indices + integer, public :: teout_idx = 0 ! teout index in physics buffer integer, public :: dtcore_idx = 0 ! dtcore index in physics buffer integer, public :: dqcore_idx = 0 ! dqcore index in physics buffer @@ -793,9 +797,7 @@ subroutine check_energy_cam_chng(state, tend, name, nstep, ztodt, & if(.not. all(cpairv(:,:,:) == cpair)) then call endrun('check_energy_chng: cpairv is not allowed to vary when subcolumns are turned on') endif - local_cp_phys(:,:) = cpair - ! Note: cp_or_cv set above for pressure coordinate if (vc_dycore == ENERGY_FORMULA_DYCORE_MPAS) then ! compute cv if vertical coordinate is height: cv = cp - R @@ -920,4 +922,255 @@ subroutine check_energy_cam_fix(state, ptend, nstep, eshflx) ) end subroutine check_energy_cam_fix + +!=============================================================================== + + subroutine enthalpy_adjustment(ncol, lchnk, state, cam_in, cam_out, pbuf, ztodt, itim_old,& + qini,totliqini,toticeini,tend) + + ! This routine is called by routine tphysac and is only called if compute_enthalpy_flux is .true. + + use camsrfexch, only: cam_in_t, cam_out_t, get_prec_vars + use physics_buffer, only: pbuf_get_index, physics_buffer_desc, pbuf_set_field, pbuf_get_field + use cam_abortutils, only: endrun + use air_composition, only: hliq_idx, hice_idx, fliq_idx, fice_idx, num_enthalpy_vars + use air_composition, only: cpairv, cp_or_cv_dycore, te_init + use air_composition, only: thermodynamic_active_species_liq_num,thermodynamic_active_species_liq_idx + use air_composition, only: thermodynamic_active_species_ice_num,thermodynamic_active_species_ice_idx + use physconst, only: cpliq, cpice, cpwv, tmelt + use air_composition, only: t00a, h00a + use physconst, only: rga, latvap, latice + use dyn_tests_utils, only: vc_dycore + use cam_thermo, only: get_hydrostatic_energy + use physics_types, only: physics_dme_adjust_camnor, dyn_te_idx + use cam_thermo, only: cam_thermo_water_update + use cam_history, only: outfld + use cam_budget, only: thermo_budget_history + use time_manager, only: get_nstep + + ! Arguments + integer, intent(in) :: ncol, lchnk + type(physics_state), intent(inout) :: state + type(cam_in_t), intent(in ) :: cam_in + type(cam_out_t), intent(inout) :: cam_out + type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), intent(in) :: ztodt + integer, intent(in) :: itim_old + real(r8), dimension(pcols,pver), intent(in) :: qini, totliqini, toticeini + type(physics_tend ) , intent(inout) :: tend + + ! Local variables + integer:: enthalpy_prec_bc_idx, enthalpy_prec_ac_idx, enthalpy_evop_idx + real(r8), dimension(:,:), pointer :: enthalpy_prec_bc + real(r8), dimension(pcols,num_enthalpy_vars) :: enthalpy_prec_ac + real(r8), dimension(pcols) :: fliq_tot, fice_tot + + integer:: dp_ntprp_idx, dp_ntsnp_idx + real(r8), dimension(:,:), pointer :: dp_ntprp, dp_ntsnp + integer:: qrain_mg_idx,qsnow_mg_idx + real(r8), dimension(:,:), pointer :: qrain_mg, qsnow_mg + + real(r8), dimension(pcols) :: te , se , po , ke + real(r8), dimension(pcols) :: te_endphys, se_endphys, po_endphys, ke_endphys + real(r8), dimension(pcols) :: te_dme , se_dme , po_dme , ke_dme + real(r8), dimension(pcols) :: te_enth_fix , se_enth_fix , po_enth_fix , ke_enth_fix + real(r8), dimension(pcols) :: fct_bc_tot, fct_ac_tot + real(r8), dimension(pcols) :: enthalpy_heating_fix_bc, enthalpy_heating_fix_ac + + real(r8), dimension(pcols) :: dEdt_physics + real(r8), dimension(pcols) :: dEdt_dme + real(r8), dimension(pcols) :: dEdt_cpdycore + real(r8), dimension(pcols) :: dEdt_enth_fix, dEdt_efix + real(r8), dimension(pcols) :: constant_latent_heat_surface !xxx diagnostics + real(r8), dimension(pcols) :: variable_latent_heat_surface_cpice_term !xxx diagnostics + real(r8), dimension(pcols) :: variable_latent_heat_surface_ls_term !xxx diagnostics + real(r8), dimension(pcols) :: variable_latent_heat_surface_lf_term !xxx diagnostics + real(r8), dimension(pcols) :: enthalpy_flux_atm, enthalpy_flux_ocn + real(r8), dimension(pcols,pver) :: tmp_t, pdel_rf, qinp, totliqinp, toticeinp + real(r8), dimension(pcols) :: zero, dsema, dcp_heat, iedme + real(r8), dimension(pcols) :: water_flux_bc, water_flux_ac, enthalpy_flux_bc, enthalpy_flux_ac + real(r8), dimension(pcols) :: eflx_out + real(r8), dimension(pcols) :: mflx_out + real(r8), dimension(pcols) :: hevap_atm, hevap_ocn + real(r8), dimension(pcols) :: tevp, tprc, nocnfrc + + real(r8), dimension(pcols,pver) :: rnsrc_pbc, snsrc_pbc + real(r8), dimension(pcols,pver) :: rnsrc_pac, snsrc_pac + real(r8), dimension(pcols,pver) :: rnsrc_tot, snsrc_tot + real(r8), dimension(pcols) :: watrerr,rainerr,snowerr + + integer nstep, ixq, m, m_cnst + real(r8), dimension(pcols,pver) :: fct_bc, fct_ac + real(r8), dimension(pcols,pver) :: scale_cpdry_cpdycore, ttend_hfix + real(r8), parameter :: eps=1.E-10_r8 + logical , parameter :: debug_enthalpy=.false. + integer :: i, k + real(r8):: tot, wgt_bc, wgt_ac + !----------------------------------------------------------------------------- + + nstep = get_nstep() + zero(:)=0._r8 + + ! scale temperature for consistency with dycore (partial adj. after cp update done implicitly in dme) + do k = 1, pver + do i = 1, ncol + scale_cpdry_cpdycore(i,k) = cpairv(i,k,lchnk)/cp_or_cv_dycore(i,k,lchnk) + state%T (i,k) = state%temp_ini(i,k)+scale_cpdry_cpdycore(i,k)*(state%T(i,k)- state%temp_ini(i,k)) + tend%dtdt(i,k) = scale_cpdry_cpdycore(i,k)*tend%dtdt(i,k) + end do + end do + + !------------------------------------------------------------------------------------------- + ! from this point onwards computation consistent with variable latent heat total energy formula + ! Equation 78 in https://agupubs.onlinelibrary.wiley.com/doi/full/10.1029/2022MS003117 + !------------------------------------------------------------------------------------------- + + !=== start computation of material enthalpy fluxes === + ! evaporation enthalpy flux + enthalpy_evop_idx = pbuf_get_index('ENTHALPY_EVOP' , errcode=i) + if (enthalpy_evop_idx==0) then + call endrun("pbufs for enthalpy evap flux not allocated") + end if + + ! using merged quantities, for atmospheric mat.enthalpy flux (used in check_energy) + if (minval(cam_in%ts(:ncol)).gt.0._r8) then + hevap_atm(:ncol) = cam_in%cflx(:ncol,1)*(cpwv*(cam_in%ts (:ncol)-t00a)+(cpliq*t00a+h00a)) ! into atm + tevp(:ncol)= cam_in%ts(:ncol) + ! for ocean-only mat.enthalpy flux (passed to ocean) + hevap_ocn (:ncol)= cam_in%evap_ocn(:ncol)*(cpwv*(cam_in%sst(:ncol)-t00a)+(cpliq*t00a+h00a)) + else ! not great but better than zeros + hevap_atm (:ncol)= cam_in%cflx(:ncol,1)*(cpwv*(state%t(:ncol,pver)-t00a)+(cpliq*t00a+h00a)) ! into atm + tevp (:ncol)= state%t(:ncol,pver) + hevap_ocn (:ncol)= hevap_atm(:ncol) ! out of ocn + endif + call pbuf_set_field(pbuf, enthalpy_evop_idx, hevap_ocn) + + !------------------------------------------------------------------ + ! compute precipitation fluxes and set associated physics buffers + !------------------------------------------------------------------ + enthalpy_prec_bc_idx = pbuf_get_index('ENTHALPY_PREC_BC', errcode=i) + enthalpy_prec_ac_idx = pbuf_get_index('ENTHALPY_PREC_AC', errcode=i) + if (enthalpy_prec_bc_idx==0.or.enthalpy_prec_ac_idx==0) then + call endrun("pbufs for enthalpy precip flux not allocated") + end if + call pbuf_get_field(pbuf, enthalpy_prec_bc_idx, enthalpy_prec_bc) + call get_prec_vars(ncol,pbuf,fliq=fliq_tot,fice=fice_tot) + + ! fliq_tot holds liquid precipitation from tphysbc and tphysac; idem for ice + enthalpy_prec_ac(:ncol,fice_idx) = fice_tot(:ncol)-enthalpy_prec_bc(:ncol,fice_idx) + enthalpy_prec_ac(:ncol,fliq_idx) = fliq_tot(:ncol)-enthalpy_prec_bc(:ncol,fliq_idx) + + ! compute precipitation enthalpy fluxes from tphysbc + tprc(:ncol) = cam_out%tbot(:ncol) + + ! correct for reference T of latent heats (liquid reference state) + enthalpy_prec_ac(:ncol,hice_idx) = -enthalpy_prec_ac(:ncol,fice_idx)*(cpice*(tprc(:ncol)-t00a)+(cpliq*t00a+h00a)) + enthalpy_prec_ac(:ncol,hliq_idx) = -enthalpy_prec_ac(:ncol,fliq_idx)*(cpliq*(tprc(:ncol)-t00a)+(cpliq*t00a+h00a)) + call pbuf_set_field(pbuf, enthalpy_prec_ac_idx, enthalpy_prec_ac) + + ! compute total enthalpy flux + enthalpy_flux_bc (:ncol) = enthalpy_prec_bc(:ncol,hliq_idx)+enthalpy_prec_bc(:ncol,hice_idx) + enthalpy_flux_ac (:ncol) = enthalpy_prec_ac(:ncol,hliq_idx)+enthalpy_prec_ac(:ncol,hice_idx) & + + hevap_atm(:ncol) + water_flux_bc(:ncol) = enthalpy_prec_bc(:ncol,fliq_idx)+enthalpy_prec_bc(:ncol,fice_idx) + water_flux_ac(:ncol) = enthalpy_prec_ac(:ncol,fliq_idx)+enthalpy_prec_ac(:ncol,fice_idx) & + - cam_in%cflx(:ncol,1) + enthalpy_flux_atm(:ncol) = enthalpy_prec_bc(:ncol,hliq_idx)+enthalpy_prec_bc(:ncol,hice_idx) & + + enthalpy_prec_ac(:ncol,hliq_idx)+enthalpy_prec_ac(:ncol,hice_idx) & + + hevap_atm(:ncol) + enthalpy_flux_ocn(:ncol) = enthalpy_prec_bc(:ncol,hliq_idx)+enthalpy_prec_bc(:ncol,hice_idx) & + + enthalpy_prec_ac(:ncol,hliq_idx)+enthalpy_prec_ac(:ncol,hice_idx) & + + hevap_ocn(:ncol) + enthalpy_flux_ocn(:ncol) = cam_in%ocnfrac(:ncol)*enthalpy_flux_ocn(:ncol) + + if (debug_enthalpy) then + call outfld("enth_prec_ac_hice" , enthalpy_prec_ac(:,hice_idx) , pcols ,lchnk ) + call outfld("enth_prec_ac_hliq" , enthalpy_prec_ac(:,hliq_idx) , pcols ,lchnk ) + call outfld("enth_prec_bc_hice" , enthalpy_prec_bc(:,hice_idx) , pcols ,lchnk ) + call outfld("enth_prec_bc_hliq" , enthalpy_prec_bc(:,hliq_idx) , pcols ,lchnk ) + call outfld("enth_prec_ac_fice" , enthalpy_prec_ac(:,fice_idx) , pcols ,lchnk ) + call outfld("enth_prec_ac_fliq" , enthalpy_prec_ac(:,fliq_idx) , pcols ,lchnk ) + call outfld("enth_prec_bc_fice" , enthalpy_prec_bc(:,fice_idx) , pcols ,lchnk ) + call outfld("enth_prec_bc_fliq" , enthalpy_prec_bc(:,fliq_idx) , pcols ,lchnk ) + call outfld("enth_hevap_atm" , hevap_atm (:) , pcols ,lchnk ) + call outfld("enth_hevap_ocn" , hevap_ocn (:) , pcols ,lchnk ) + endif + !=== end computation of material enthalpy fluxes === + + !+++ diags + ! compute total energy after physics using equation 78 + call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & + state%pdel(1:ncol,1:pver), cp_or_cv_dycore(:ncol,:,lchnk), & + state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver),& + vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & + te = te_endphys(:ncol), se=se_endphys(:ncol), po=po_endphys(:ncol), ke=ke_endphys(:ncol)) + ! the column integrated total energy change should match accumlated te_tnd: + ! dEdt_physics=te_tnd + call outfld ('te_tnd',tend%te_tnd , pcols, lchnk) + dEdt_physics(:ncol) = (te_endphys(:ncol)-te_init(:ncol,1,lchnk))/ztodt + call outfld ('dEdt_physics', dEdt_physics, pcols, lchnk) + !--- sgaid + + !+ get pbuf fields for precip + dp_ntprp_idx = pbuf_get_index('dp_ntprp',errcode=i) !prec production from ZM + dp_ntsnp_idx = pbuf_get_index('dp_ntsnp',errcode=i) !snow production from ZM + call pbuf_get_field(pbuf, dp_ntprp_idx , dp_ntprp) + call pbuf_get_field(pbuf, dp_ntsnp_idx , dp_ntsnp) + qrain_mg_idx = pbuf_get_index('qrain_mg',errcode=i) !rain production from MG + qsnow_mg_idx = pbuf_get_index('qsnow_mg',errcode=i) !snow production from MG + call pbuf_get_field(pbuf, qrain_mg_idx, qrain_mg) + call pbuf_get_field(pbuf, qsnow_mg_idx, qsnow_mg) + rnsrc_pbc(:ncol,:) = dp_ntprp(:ncol,:)-dp_ntsnp(:ncol,:) + snsrc_pbc(:ncol,:) = dp_ntsnp(:ncol,:) + rnsrc_pac(:ncol,:) = qrain_mg(:ncol,:) + snsrc_pac(:ncol,:) = qsnow_mg(:ncol,:) + rnsrc_tot(:ncol,:) = rnsrc_pbc(:ncol,:)+rnsrc_pac(:ncol,:) + snsrc_tot(:ncol,:) = snsrc_pbc(:ncol,:)+snsrc_pac(:ncol,:) + !- picerp rof sdleif fubp teg + + ! Adjust the dry mass in each layer back to the value of physics input state + ! Adjust air specific enthalpy accordingly. Diagnose boundary enthalpy flux. + ! Author: Thomas Toniazzo (17.07.21) + call physics_dme_adjust_camnor(state, tend, qini, totliqini, toticeini, ztodt, & + ntrnprd=rnsrc_tot*ztodt, & + ntsnprd=snsrc_tot*ztodt, & + tevap=tevp, & + tprec=tprc, & + mflx=water_flux_bc+water_flux_ac, & + eflx=enthalpy_flux_atm, & + mflx_out=mflx_out, & + eflx_out=eflx_out, & + ent_tnd=dsema, & + pdel_rf=pdel_rf ) + + call outfld('IETEND_DME', dsema , pcols, lchnk) + call outfld('EFLX' , enthalpy_flux_atm , pcols, lchnk) + call outfld('MFLX' , water_flux_bc+water_flux_ac , pcols, lchnk) + + ! compute and store new column-integrated enthalpy and associated tendency + call get_hydrostatic_energy(state%q(1:ncol,1:pver,1:pcnst),.true., & + state%pdel(1:ncol,1:pver), cp_or_cv_dycore(:ncol,:,lchnk), & + state%u(1:ncol,1:pver), state%v(1:ncol,1:pver), state%T(1:ncol,1:pver),& + vc_dycore, ptop=state%pintdry(1:ncol,1), phis = state%phis(1:ncol), & + te = te(:ncol), se=se(:ncol), po=po(:ncol), ke=ke(:ncol)) + + ! Save final energy for use with global fixer in next timestep -- note sign conventions, and coupling-dependent options + ! subtract from te the h flux (sign: into atm) that is *not* passed to surface components + ! and also remove enthalpy of run-off (if added to BLOM) + state%te_cur(:ncol,dyn_te_idx) = te(:ncol) & + - ztodt*(enthalpy_flux_atm(:ncol) - enthalpy_flux_ocn(:ncol) - cam_in%hrof(:ncol)) + tend%te_tnd(:ncol) = tend%te_tnd(:ncol) + (enthalpy_flux_ocn(:ncol) + cam_in%hrof(:ncol)) ! B. with run-off + + if (thermo_budget_history) then + call tot_energy_phys(state, 'phAM') + call tot_energy_phys(state, 'dyAM', vc=vc_dycore) + endif + + call pbuf_set_field(pbuf, teout_idx, state%te_cur(:,dyn_te_idx), (/1,itim_old/),(/pcols,1/)) + ! the amount of total energy we need energy fixer to fix (associated with enthalpy flux) + dEdt_efix(:ncol) = (state%te_cur(:ncol,dyn_te_idx)-te (:ncol))/ztodt + call outfld("dEdt_efix_physics" , dEdt_efix , pcols ,lchnk ) + + end subroutine enthalpy_adjustment + end module check_energy diff --git a/src/physics/cam/phys_control.F90 b/src/physics/cam/phys_control.F90 index 210c2a1d66..7528bd11d4 100644 --- a/src/physics/cam/phys_control.F90 +++ b/src/physics/cam/phys_control.F90 @@ -71,12 +71,15 @@ module phys_control logical :: history_dust = .false. logical :: history_scwaccm_forcing = .false. logical :: history_chemspecies_srf = .false. + logical, public, protected :: history_aerosol_base = .true. logical, public, protected :: history_aerosol_decomposed = .false. logical, public, protected :: history_gas = .false. logical, public, protected :: history_aerosol_forcing = .false. logical, public, protected :: history_aerosol_radiation = .false. logical, public, protected :: history_aerosol_debug_output = .false. +logical, public, protected :: history_enthalpy_flux = .false. + logical :: do_clubb_sgs logical :: do_hb_above_clubb = .false. ! enable HB vertical mixing above clubb top @@ -143,6 +146,7 @@ subroutine phys_ctl_readnl(nlfile) history_clubb, history_dust, & history_cesm_forcing, history_scwaccm_forcing, history_chemspecies_srf, history_aerosol_base, history_aerosol_debug_output, & history_aerosol_decomposed, history_gas, history_aerosol_forcing, history_aerosol_radiation, & + history_enthalpy_flux, & do_clubb_sgs, state_debug_checks, use_hetfrz_classnuc, use_gw_oro, use_gw_front, & use_gw_front_igw, use_gw_convect_dp, use_gw_convect_sh, use_gw_movmtn_pbl, cld_macmic_num_steps, & offline_driver, convproc_do_aer, cam_snapshot_before_num, cam_snapshot_after_num, & @@ -200,6 +204,7 @@ subroutine phys_ctl_readnl(nlfile) call mpi_bcast(history_aerosol_debug_output,1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(history_dust, 1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(history_scwaccm_forcing, 1, mpi_logical, masterprocid, mpicom, ierr) + call mpi_bcast(history_enthalpy_flux, 1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(do_clubb_sgs, 1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(state_debug_checks, 1, mpi_logical, masterprocid, mpicom, ierr) call mpi_bcast(use_hetfrz_classnuc, 1, mpi_logical, masterprocid, mpicom, ierr) diff --git a/src/physics/cam/physics_types.F90 b/src/physics/cam/physics_types.F90 index 9fddfdd811..4fb5efe0ca 100644 --- a/src/physics/cam/physics_types.F90 +++ b/src/physics/cam/physics_types.F90 @@ -7,6 +7,7 @@ module physics_types use ppgrid, only: pcols, pver use constituents, only: pcnst, qmin, cnst_name, cnst_get_ind use geopotential, only: geopotential_t + use physconst, only: cpliq, cpwv use physconst, only: zvir, gravit, cpair, rair use air_composition, only: cpairv, rairv use phys_grid, only: get_ncols_p, get_rlon_all_p, get_rlat_all_p, get_gcol_all_p @@ -14,6 +15,7 @@ module physics_types use cam_abortutils, only: endrun use phys_control, only: waccmx_is use shr_const_mod, only: shr_const_rwv + use spmd_utils, only: masterproc implicit none private ! Make default type private to the module @@ -32,6 +34,7 @@ module physics_types public physics_ptend_init public physics_state_set_grid public physics_dme_adjust ! adjust dry mass and energy for change in water + public physics_dme_adjust_camnor ! adjust dry mass and energy for change in water public physics_state_copy ! copy a physics_state object public physics_ptend_copy ! copy a physics_ptend object public physics_ptend_sum ! accumulate physics_ptend objects @@ -53,7 +56,14 @@ module physics_types public physics_cnst_limit ! apply limiters to constituents (waccmx) !------------------------------------------------------------------------------- integer, parameter, public :: phys_te_idx = 1 - integer ,parameter, public :: dyn_te_idx = 2 + integer, parameter, public :: dyn_te_idx = 2 + + integer, parameter, public :: num_hflx = 4 + + integer, parameter, public :: ihrain = 1 ! index for enthalpy flux associated with liquid precipitation + integer, parameter, public :: ihsnow = 2 ! index for enthalpy flux associated with frozen precipiation + integer, parameter, public :: ifrain = 3 ! index for flux of liquid precipitation + integer, parameter, public :: ifsnow = 4 ! index for flux of frozen precipitation type physics_state integer :: & @@ -101,7 +111,7 @@ module physics_types ! (dyn_te_idx) dycore total energy computed in physics te_ini, &! vertically integrated total (kinetic + static) energy of initial state te_cur ! vertically integrated total (kinetic + static) energy of current state - real(r8), dimension(:), allocatable :: & + real(r8), dimension(: ),allocatable :: & tw_ini, &! vertically integrated total water of initial state tw_cur ! vertically integrated total water of new state real(r8), dimension(:,:),allocatable :: & @@ -123,9 +133,11 @@ module physics_types integer :: psetcols=0 ! max number of columns set- if subcols = pcols*psubcols, else = pcols real(r8), dimension(:,:),allocatable :: dtdt, dudt, dvdt + real(r8), dimension(:,:),allocatable :: s_dme, qt_dme real(r8), dimension(:), allocatable :: flx_net real(r8), dimension(:), allocatable :: & te_tnd, &! cumulative boundary flux of total energy + te_sen, &! cumulative sensible heat flux tw_tnd ! cumulative boundary flux of total water end type physics_tend @@ -169,6 +181,7 @@ module physics_types end type physics_ptend + logical :: levels_are_moist=.true. ! TODO: put in namelist? !=============================================================================== contains @@ -204,14 +217,17 @@ subroutine physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, psetcol end subroutine physics_type_alloc !=============================================================================== - subroutine physics_update(state, ptend, dt, tend) + subroutine physics_update(state, ptend, dt, tend ) !----------------------------------------------------------------------- ! Update the state and or tendency structure with the parameterization tendencies !----------------------------------------------------------------------- use scamMod, only: scm_crm_mode, single_column use phys_control, only: phys_getopts - use cam_thermo, only: cam_thermo_dry_air_update ! Routine which updates physconst variables (WACCM-X) - use air_composition, only: dry_air_species_num, thermodynamic_active_species_num, thermodynamic_active_species_idx + use cam_thermo, only: cam_thermo_dry_air_update ! Routine which updates physconst variables (WACCM-X) + use cam_thermo, only: get_conserved_energy, inv_conserved_energy + use air_composition, only: dry_air_species_num + use air_composition, only: thermodynamic_active_species_num, thermodynamic_active_species_idx + use air_composition, only: compute_enthalpy_flux use qneg_module , only: qneg3 !------------------------------Arguments-------------------------------- @@ -233,6 +249,8 @@ subroutine physics_update(state, ptend, dt, tend) integer :: ixh, ixh2 ! constituent indices for H, H2 logical :: derive_new_geopotential ! derive new geopotential fields? + real(r8) :: te(state%psetcols,pver),t_tmp(state%psetcols,pver),pdel(state%psetcols,pver) + real(r8) :: zvirv(state%psetcols,pver) ! Local zvir array pointer real(r8),allocatable :: cpairv_loc(:,:) @@ -411,16 +429,36 @@ subroutine physics_update(state, ptend, dt, tend) !------------------------------------------------------------------------------------------------------------- ! Update temperature from dry static energy (moved from above for WACCM-X so updating after cpairv_loc update) !------------------------------------------------------------------------------------------------------------- - if(ptend%ls) then - do k = ptend%top_level, ptend%bot_level - state%t(:ncol,k) = state%t(:ncol,k) + ptend%s(:ncol,k)*dt/cpairv_loc(:ncol,k) - if (present(tend)) & - tend%dtdt(:ncol,k) = tend%dtdt(:ncol,k) + ptend%s(:ncol,k)/cpairv_loc(:ncol,k) - end do + + if(compute_enthalpy_flux) then + !use conserved energy (pe and te are output variables in get_conserved_energy call) + call get_conserved_energy(levels_are_moist, ptend%top_level, ptend%bot_level, & + cpairv_loc(:ncol,:), state%T(:ncol,:), state%q(:ncol,:,:), state%pdel(:ncol,:), & + pdel(:ncol,:), te(:ncol,:)) + te(:ncol,ptend%top_level:ptend%bot_level) = te(:ncol,ptend%top_level:ptend%bot_level) + & + ptend%s(:ncol,ptend%top_level:ptend%bot_level)*dt + call inv_conserved_energy(levels_are_moist, ptend%top_level, ptend%bot_level, & + te(:ncol,:), cpairv_loc(:ncol,:), state%q(:ncol,:,:), state%pdel(:ncol,:), & + pdel(:ncol,:), t_tmp(:ncol,:)) + if (present(tend)) then + tend%dtdt(:ncol,ptend%top_level:ptend%bot_level) = tend%dtdt(:ncol,ptend%top_level:ptend%bot_level) + & + (T_tmp(:ncol,ptend%top_level:ptend%bot_level) - & + state%t(:ncol,ptend%top_level:ptend%bot_level))/dt + end if + state%T(:ncol,ptend%top_level:ptend%bot_level) = T_tmp(:ncol,ptend%top_level:ptend%bot_level) + else + do k = ptend%top_level, ptend%bot_level + state%t(:ncol,k) = state%t(:ncol,k) + ptend%s(:ncol,k)*dt/cpairv_loc(:ncol,k) + if (present(tend)) then + tend%dtdt(:ncol,k) = tend%dtdt(:ncol,k) + ptend%s(:ncol,k)/cpairv_loc(:ncol,k) + end if + end do + endif + end if - ! Derive new geopotential fields if heating or water species tendency not 0. + ! Derive new geopotential fields if heating or water tendency not 0. derive_new_geopotential = .false. if(ptend%ls) then ! Heating tendency not 0 @@ -552,9 +590,9 @@ 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%tw_ini(:ncol), is_nan=.false., & + 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., & + call shr_assert_in_domain(state%tw_cur(:ncol ), is_nan=.false., & varname="state%tw_cur", msg=msg) call shr_assert_in_domain(state%temp_ini(:ncol,:), is_nan=.false., & varname="state%temp_ini", msg=msg) @@ -630,9 +668,9 @@ 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) - call shr_assert_in_domain(state%tw_ini(:ncol), lt=posinf_r8, gt=neginf_r8, & + 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, & + call shr_assert_in_domain(state%tw_cur(:ncol ), lt=posinf_r8, gt=neginf_r8, & varname="state%tw_cur", msg=msg) call shr_assert_in_domain(state%temp_ini(:ncol,:), lt=posinf_r8, gt=neginf_r8, & varname="state%temp_ini", msg=msg) @@ -1319,9 +1357,56 @@ subroutine physics_dme_adjust(state, tend, qini, liqini, iceini, dt) end subroutine physics_dme_adjust -!----------------------------------------------------------------------- +!=============================================================================== + + subroutine physics_dme_adjust_camnor(state, tend, qini, liqini, iceini, dt, & + ntrnprd, ntsnprd, tevap, tprec, mflx, eflx, eflx_out, mflx_out, & + ent_tnd, pdel_rf) + + ! Purpose: Diagnose boundary enthalpy flux and local heating rates associated to + ! atmospheric moisture change: Author: Thomas Toniazzo (17.07.21) + + use dme_adjust_camnor, only: dme_adjust_camnor_run + ! + ! Arguments + ! + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + real(r8), intent(in) :: qini(pcols,pver) ! initial specific humidity + real(r8), intent(in) :: liqini(pcols,pver) ! initial total liquid + real(r8), intent(in) :: iceini(pcols,pver) ! initial total ice + real(r8), intent(in) :: dt + real(r8), intent(in) :: ntrnprd(pcols,pver) ! net precip (liq+ice) production in layer + real(r8), intent(in) :: ntsnprd(pcols,pver) ! net snow production in layer + real(r8), intent(in) :: tevap(pcols) ! temperature of surface evaporation + real(r8), intent(in) :: tprec(pcols) ! temperature of surface precipitation + real(r8), intent(in) :: mflx(pcols) ! mass flux for use in check_energy + real(r8), intent(in) :: eflx(pcols) ! energy flux for use in check_energy + real(r8), intent(out) :: mflx_out(pcols) ! column (surfce) enthalpy flux from bflx (sanity check) + real(r8), intent(out) :: eflx_out(pcols) ! column (surfce) enthalpy flux from bflx (sanity check) + real(r8), intent(out) :: ent_tnd(pcols) ! column-integrated enthalpy tendency + real(r8), intent(out) :: pdel_rf(pcols,pver) ! ratio old pdel / new pdel + !----------------------------------------------------------------------- + + if (state%psetcols /= pcols) then + call endrun('physics_dme_adjust_camnor: cannot pass in a state which has sub-columns') + end if + + call dme_adjust_camnor_run(state%lchnk, state%ncol, & + state%psetcols, state%pint, state%pmid, & + state%pdel, state%rpdel, state%pdeldry, state%lnpint, state%lnpmid, & + state%ps, state%phis, state%zm, state%zi, & + state%t, state%u, state%v, state%q, state%s, & + tend%dudt, tend%dvdt, tend%dtdt, & + qini, liqini, iceini, dt, & + ntrnprd, ntsnprd, tevap, tprec, mflx, eflx, eflx_out, mflx_out, & + ent_tnd, pdel_rf) + + end subroutine physics_dme_adjust_camnor !=============================================================================== + + subroutine physics_state_copy(state_in, state_out) use ppgrid, only: pver, pverp @@ -1357,10 +1442,10 @@ subroutine physics_state_copy(state_in, state_out) state_out%ps(i) = state_in%ps(i) state_out%phis(i) = state_in%phis(i) end do - state_out%te_ini(:ncol,:) = state_in%te_ini(:ncol,:) - state_out%te_cur(:ncol,:) = state_in%te_cur(:ncol,:) - state_out%tw_ini(:ncol) = state_in%tw_ini(:ncol) - state_out%tw_cur(:ncol) = state_in%tw_cur(:ncol) + state_out%te_ini (:ncol,:) = state_in%te_ini (:ncol,:) + state_out%te_cur (:ncol,:) = state_in%te_cur (:ncol,:) + state_out%tw_ini (:ncol ) = state_in%tw_ini (:ncol ) + state_out%tw_cur (:ncol ) = state_in%tw_cur (:ncol ) do k = 1, pver do i = 1, ncol @@ -1435,11 +1520,14 @@ subroutine physics_tend_init(tend) call endrun('physics_tend_init: tend must be allocated before it can be initialized') end if + tend%s_dme = 0._r8 + tend%qt_dme = 0._r8 tend%dtdt = 0._r8 tend%dudt = 0._r8 tend%dvdt = 0._r8 tend%flx_net = 0._r8 tend%te_tnd = 0._r8 + tend%te_sen = 0._r8 tend%tw_tnd = 0._r8 end subroutine physics_tend_init @@ -1675,10 +1763,10 @@ subroutine physics_state_alloc(state,lchnk,psetcols) allocate(state%te_cur(psetcols,2), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%te_cur') - allocate(state%tw_ini(psetcols), stat=ierr) + allocate(state%tw_ini(psetcols ), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%tw_ini') - allocate(state%tw_cur(psetcols), stat=ierr) + allocate(state%tw_cur(psetcols ), stat=ierr) if ( ierr /= 0 ) call endrun('physics_state_alloc error: allocation error for state%tw_cur') allocate(state%temp_ini(psetcols,pver), stat=ierr) @@ -1726,12 +1814,12 @@ subroutine physics_state_alloc(state,lchnk,psetcols) state%lnpintdry(:,:) = inf state%zi(:,:) = inf - state%te_ini(:,:) = inf - state%te_cur(:,:) = inf - state%tw_ini(:) = inf - state%tw_cur(:) = inf + state%te_ini (:,:) = inf + state%te_cur (:,:) = inf + state%tw_ini (: ) = inf + state%tw_cur (: ) = inf state%temp_ini(:,:) = inf - state%z_ini(:,:) = inf + state%z_ini (:,:) = inf end subroutine physics_state_alloc @@ -1872,6 +1960,11 @@ subroutine physics_tend_alloc(tend,psetcols) tend%psetcols = psetcols + allocate(tend%s_dme(psetcols,pver), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%s_dme') + allocate(tend%qt_dme(psetcols,pver), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%qt_dme') + allocate(tend%dtdt(psetcols,pver), stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%dtdt') @@ -1887,14 +1980,24 @@ subroutine physics_tend_alloc(tend,psetcols) allocate(tend%te_tnd(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%te_tnd') + allocate(tend%te_sen(psetcols), stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%te_sen') + + !allocate(tend%te_lat(psetcols), stat=ierr) + !if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%te_lat') + allocate(tend%tw_tnd(psetcols), stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_alloc error: allocation error for tend%tw_tnd') + tend%s_dme (:,:)= inf + tend%qt_dme(:,:)= inf tend%dtdt(:,:) = inf tend%dudt(:,:) = inf tend%dvdt(:,:) = inf tend%flx_net(:) = inf tend%te_tnd(:) = inf + tend%te_sen(:) = inf + !tend%te_lat(:) = inf tend%tw_tnd(:) = inf end subroutine physics_tend_alloc @@ -1908,6 +2011,12 @@ subroutine physics_tend_dealloc(tend) type(physics_tend), intent(inout) :: tend integer :: ierr = 0 + deallocate(tend%s_dme, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%s_dme') + + deallocate(tend%qt_dme, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%qt_dme') + deallocate(tend%dtdt, stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%dtdt') @@ -1923,6 +2032,12 @@ subroutine physics_tend_dealloc(tend) deallocate(tend%te_tnd, stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%te_tnd') + deallocate(tend%te_sen, stat=ierr) + if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%te_sen') + + !deallocate(tend%te_lat, stat=ierr) + !if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%te_lat') + deallocate(tend%tw_tnd, stat=ierr) if ( ierr /= 0 ) call endrun('physics_tend_dealloc error: deallocation error for tend%tw_tnd') end subroutine physics_tend_dealloc diff --git a/src/physics/cam/qneg_module.F90 b/src/physics/cam/qneg_module.F90 index 773bf220a5..432324e1bd 100644 --- a/src/physics/cam/qneg_module.F90 +++ b/src/physics/cam/qneg_module.F90 @@ -309,7 +309,7 @@ subroutine qneg3 (subnam, idx, ncol, ncold, lver, lconst_beg, & end subroutine qneg3 subroutine qneg4 (subnam, lchnk, ncol, ztodt, & - qbot, srfrpdel, shflx, lhflx, qflx) + qbot, srfrpdel, shflx, lhflx, qflx, seflx) !----------------------------------------------------------------------- ! ! Purpose: @@ -325,7 +325,7 @@ subroutine qneg4 (subnam, lchnk, ncol, ztodt, & ! Author: J. Olson ! !----------------------------------------------------------------------- - use physconst, only: gravit, latvap + use physconst, only: gravit, latvap, latice use constituents, only: qmin use cam_history, only: outfld @@ -343,9 +343,10 @@ subroutine qneg4 (subnam, lchnk, ncol, ztodt, & ! ! Input/Output arguments ! - real(r8), intent(inout) :: shflx(ncol) ! Surface sensible heat flux (J/m2/s) - real(r8), intent(inout) :: lhflx(ncol) ! Surface latent heat flux (J/m2/s) - real(r8), intent(inout) :: qflx (ncol,pcnst) ! surface water flux (kg/m^2/s) + real(r8), intent(inout) :: shflx(ncol) ! Surface sensible heat flux (J/m2/s) + real(r8), intent(inout) :: lhflx(ncol) ! Surface latent heat flux (J/m2/s) + real(r8), intent(inout) :: qflx (ncol,pcnst) ! surface water flux (kg/m^2/s) + real(r8), intent(inout), optional :: seflx(ncol) ! heat flux for energy checker (ice ref.state) ! !---------------------------Local workspace----------------------------- ! @@ -390,11 +391,15 @@ subroutine qneg4 (subnam, lchnk, ncol, ztodt, & qflx (i,1) = qflx (i,1) - excess(i) lhflx(i) = lhflx(i) - excess(i)*latvap shflx(i) = shflx(i) + excess(i)*latvap + if (present(seflx)) then + seflx(i) = seflx(i) + excess(i)*(latvap+latice) + end if if (index > 0) then qneg4_warn_num(index) = qneg4_warn_num(index) + 1 end if end if end do + ! Maybe output bad values if ((cnst_outfld((2*pcnst)+1)) .and. (worst < worst_reset)) then do i = 1, ncol diff --git a/src/physics/cam/zm_conv_intr.F90 b/src/physics/cam/zm_conv_intr.F90 index 5079daa7f8..65c2a3243f 100644 --- a/src/physics/cam/zm_conv_intr.F90 +++ b/src/physics/cam/zm_conv_intr.F90 @@ -60,7 +60,9 @@ module zm_conv_intr dlfzm_idx, & ! detrained convective cloud water mixing ratio. prec_dp_idx, & snow_dp_idx, & - mconzm_idx ! convective mass flux + mconzm_idx, & ! convective mass flux + dp_ntprp_idx, & ! needed in check_energy for new enthalpy + dp_ntsnp_idx ! needed in check_energy for new enthalpy real(r8), parameter :: unset_r8 = huge(1.0_r8) real(r8) :: zmconv_c0_lnd = unset_r8 @@ -134,10 +136,14 @@ subroutine zm_conv_register ! map gathered points to chunk index call pbuf_add_field('ZM_IDEEP', 'physpkg', dtype_i4, (/pcols/), zm_ideep_idx) -! Flux of precipitation from deep convection (kg/m2/s) + ! Flux of precipitation from deep convection (kg/m2/s) call pbuf_add_field('DP_FLXPRC','global',dtype_r8,(/pcols,pverp/),dp_flxprc_idx) -! Flux of snow from deep convection (kg/m2/s) + ! Needed for check_energy for new enthalpy computations + call pbuf_add_field('dp_ntprp','physpkg',dtype_r8,(/pcols,pver /),dp_ntprp_idx) + call pbuf_add_field('dp_ntsnp','physpkg',dtype_r8,(/pcols,pver /),dp_ntsnp_idx) + + ! Flux of snow from deep convection (kg/m2/s) call pbuf_add_field('DP_FLXSNW','global',dtype_r8,(/pcols,pverp/),dp_flxsnw_idx) call pbuf_add_field('ICWMRDP', 'physpkg',dtype_r8,(/pcols,pver/),icwmrdp_idx) @@ -754,6 +760,9 @@ subroutine zm_conv_tend(pblh ,mcon ,cme , & evapcdp(:ncol,:pver) = ptend_loc%q(:ncol,:pver,1) + ! Needed in check_energy for new enthalpy computations + call pbuf_set_field(pbuf, dp_ntprp_idx, ntprprd) + call pbuf_set_field(pbuf, dp_ntsnp_idx, ntsnprd) ! ! Write out variables from zm_conv_evap_run ! diff --git a/src/physics/cam7/micro_pumas_cam.F90 b/src/physics/cam7/micro_pumas_cam.F90 index a33949034e..b11ff71c8a 100644 --- a/src/physics/cam7/micro_pumas_cam.F90 +++ b/src/physics/cam7/micro_pumas_cam.F90 @@ -219,7 +219,7 @@ module micro_pumas_cam qrain_idx=-1, qsnow_idx=-1, & nrain_idx=-1, nsnow_idx=-1, & qcsedten_idx=-1, qrsedten_idx=-1, & - qisedten_idx=-1, qssedten_idx=-1, & + qisedten_idx=-1, qssedten_idx=-1, qgsedten_idx=-1, & vtrmc_idx=-1, umr_idx=-1, & vtrmi_idx=-1, ums_idx=-1, & qcsevap_idx=-1, qisevap_idx=-1 @@ -816,6 +816,12 @@ subroutine micro_pumas_cam_register call pbuf_add_field('UMS', 'physpkg', dtype_r8, (/pcols,pver/), ums_idx) call pbuf_add_field('QCSEVAP', 'physpkg', dtype_r8, (/pcols,pver/), qcsevap_idx) call pbuf_add_field('QISEVAP', 'physpkg', dtype_r8, (/pcols,pver/), qisevap_idx) + else + call pbuf_add_field('QCSEDTEN', 'physpkg', dtype_r8, (/pcols,pver/), qcsedten_idx) + call pbuf_add_field('QRSEDTEN', 'physpkg', dtype_r8, (/pcols,pver/), qrsedten_idx) + call pbuf_add_field('QISEDTEN', 'physpkg', dtype_r8, (/pcols,pver/), qisedten_idx) + call pbuf_add_field('QSSEDTEN', 'physpkg', dtype_r8, (/pcols,pver/), qssedten_idx) + call pbuf_add_field('QGSEDTEN', 'physpkg', dtype_r8, (/pcols,pver/), qgsedten_idx) end if end subroutine micro_pumas_cam_register @@ -1511,6 +1517,7 @@ subroutine micro_pumas_cam_init(pbuf2d) if (qrsedten_idx > 0) call pbuf_set_field(pbuf2d, qrsedten_idx, 0._r8) if (qisedten_idx > 0) call pbuf_set_field(pbuf2d, qisedten_idx, 0._r8) if (qssedten_idx > 0) call pbuf_set_field(pbuf2d, qssedten_idx, 0._r8) + if (qgsedten_idx > 0) call pbuf_set_field(pbuf2d, qgsedten_idx, 0._r8) if (vtrmc_idx > 0) call pbuf_set_field(pbuf2d, vtrmc_idx, 0._r8) if (umr_idx > 0) call pbuf_set_field(pbuf2d, umr_idx, 0._r8) if (vtrmi_idx > 0) call pbuf_set_field(pbuf2d, vtrmi_idx, 0._r8) @@ -1919,6 +1926,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) real(r8) :: qrsedtenout_grid(pcols,pver) real(r8) :: qisedtenout_grid(pcols,pver) real(r8) :: qssedtenout_grid(pcols,pver) + real(r8) :: qgsedtenout_grid(pcols,pver) real(r8) :: vtrmcout_grid(pcols,pver) real(r8) :: umrout_grid(pcols,pver) real(r8) :: vtrmiout_grid(pcols,pver) @@ -1993,6 +2001,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) real(r8), pointer :: qrsedtenout_grid_ptr(:,:) real(r8), pointer :: qisedtenout_grid_ptr(:,:) real(r8), pointer :: qssedtenout_grid_ptr(:,:) + real(r8), pointer :: qgsedtenout_grid_ptr(:,:) real(r8), pointer :: vtrmcout_grid_ptr(:,:) real(r8), pointer :: umrout_grid_ptr(:,:) real(r8), pointer :: vtrmiout_grid_ptr(:,:) @@ -2258,6 +2267,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) if (qrsedten_idx > 0) call pbuf_get_field(pbuf, qrsedten_idx, qrsedtenout_grid_ptr) if (qisedten_idx > 0) call pbuf_get_field(pbuf, qisedten_idx, qisedtenout_grid_ptr) if (qssedten_idx > 0) call pbuf_get_field(pbuf, qssedten_idx, qssedtenout_grid_ptr) + if (qgsedten_idx > 0) call pbuf_get_field(pbuf, qgsedten_idx, qgsedtenout_grid_ptr) if (vtrmc_idx > 0) call pbuf_get_field(pbuf, vtrmc_idx, vtrmcout_grid_ptr) if (umr_idx > 0) call pbuf_get_field(pbuf, umr_idx, umrout_grid_ptr) if (vtrmi_idx > 0) call pbuf_get_field(pbuf, vtrmi_idx, vtrmiout_grid_ptr) @@ -2987,6 +2997,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) qisevapout_grid(:ncol,:top_lev-1) = 0._r8 qrsedtenout_grid(:ncol,:top_lev-1) = 0._r8 qssedtenout_grid(:ncol,:top_lev-1) = 0._r8 + qgsedtenout_grid(:ncol,:top_lev-1) = 0._r8 umrout_grid(:ncol,:top_lev-1) = 0._r8 umsout_grid(:ncol,:top_lev-1) = 0._r8 psacro_grid(:ncol,:top_lev-1) = 0._r8 @@ -3077,6 +3088,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) ns_grid = state_loc%q(:,:,ixnumsnow) qrsedtenout_grid(:ncol,top_lev:) = proc_rates%qrsedten qssedtenout_grid(:ncol,top_lev:) = proc_rates%qssedten + qgsedtenout_grid(:ncol,top_lev:) = proc_rates%qgsedten umrout_grid(:ncol,top_lev:) = proc_rates%umr umsout_grid(:ncol,top_lev:) = proc_rates%ums @@ -3569,6 +3581,7 @@ subroutine micro_pumas_cam_tend(state, ptend, dtime, pbuf) if (qrsedten_idx > 0) qrsedtenout_grid_ptr = qrsedtenout_grid if (qisedten_idx > 0) qisedtenout_grid_ptr = qisedtenout_grid if (qssedten_idx > 0) qssedtenout_grid_ptr = qssedtenout_grid + if (qgsedten_idx > 0) qgsedtenout_grid_ptr = qgsedtenout_grid if (vtrmc_idx > 0) vtrmcout_grid_ptr = vtrmcout_grid if (umr_idx > 0) umrout_grid_ptr = umrout_grid if (vtrmi_idx > 0) vtrmiout_grid_ptr = vtrmiout_grid diff --git a/src/utils/air_composition.F90 b/src/utils/air_composition.F90 index ed2e1a9afc..440e034f4a 100644 --- a/src/utils/air_composition.F90 +++ b/src/utils/air_composition.F90 @@ -1,7 +1,8 @@ -! air_composition module defines major species of the atmosphere and manages -! the physical properties that are dependent on the composition of air module air_composition + ! air_composition module defines major species of the atmosphere and manages + ! the physical properties that are dependent on the composition of air + use shr_kind_mod, only: r8 => shr_kind_r8 use cam_abortutils, only: endrun @@ -9,6 +10,7 @@ module air_composition private save + public :: air_composition_register ! sets module variable compute_enthalpy_flux public :: air_composition_readnl public :: air_composition_init public :: dry_air_composition_update @@ -24,12 +26,21 @@ module air_composition public :: get_R ! get_mbarv: molecular weight of dry air public :: get_mbarv + ! + ! enthalpy variables in physics buffer + ! + integer, parameter, public :: num_enthalpy_vars = 4 ! index for enthalpy flux associated with liquid precipitation + integer, parameter, public :: hliq_idx = 1 ! index for enthalpy flux associated with liquid precipitation + integer, parameter, public :: hice_idx = 2 ! index for enthalpy flux associated with frozen precipiation + integer, parameter, public :: fliq_idx = 3 ! index for flux of liquid precipitation + integer, parameter, public :: fice_idx = 4 ! index for flux of frozen precipitation private :: air_species_info integer, parameter :: unseti = -HUGE(1) real(r8), parameter :: unsetr = HUGE(1.0_r8) + ! composition of air ! integer, parameter :: num_names_max = 20 ! Should match namelist definition @@ -38,7 +49,6 @@ module air_composition integer, protected, public :: dry_air_species_num integer, protected, public :: water_species_in_air_num - logical, protected, public :: compute_enthalpy_flux ! Thermodynamic variables integer, protected, public :: thermodynamic_active_species_num = unseti @@ -68,8 +78,9 @@ module air_composition integer, allocatable, protected, public :: thermodynamic_active_species_ice_idx(:) ! thermodynamic_active_species_ice_idx_dycore: index of ice water species integer, allocatable, public :: thermodynamic_active_species_ice_idx_dycore(:) - ! enthalpy_reference_state: choices: 'ice', 'liq', 'wv' - character(len=3), public, protected :: enthalpy_reference_state = 'xxx' + ! enthalpy_reference_state: choices: 'ice', 'liq', 'vap' + ! 'wv'->'vap' (stick to three characters, 'water' is presumably implicit in all of these...) + character(len=3), public, protected :: enthalpy_reference_state = 'ice' integer, protected, public :: wv_idx = -1 ! Water vapor index @@ -83,6 +94,16 @@ module air_composition real(r8), public, protected :: n2_mwi = unsetr ! Inverse mol. weight of N2 real(r8), public, protected :: mbar = unsetr ! Mean mass at mid level + ! explicitly declare reference enthalpies and temperatures for atmosphere and ocean + ! only used if compute_enthalpy_flux is true + logical , public, protected :: compute_enthalpy_flux = .false. ! obtained from nuopc mediator + real(r8), public, protected :: t00o = unsetr ! Water enthalpy reference temperature, ocean (K) + real(r8), public, protected :: t00a = unsetr ! Water enthalpy reference temperature, atmosphere (K) + real(r8), public, protected :: h00o = unsetr ! Material enthalpy zero, liquid reference state, ocean water (J/kg) + real(r8), public, protected :: h00a = unsetr ! Material enthalpy zero, liquid reference state, atmos water (J/kg) + real(r8), public, protected :: h00a_vap = unsetr ! Material enthalpy zero, vapor reference state, atmos (J/kg) + real(r8), public, protected :: h00a_ice = unsetr ! Material enthalpy zero, vapor reference state, atmos (J/kg) + ! coefficients in expressions for molecular diffusion coefficients ! kv1,..,kv3 are coefficients for kmvis calculation ! kc1,..,kc3 are coefficients for kmcnd calculation @@ -109,6 +130,7 @@ module air_composition ! cp_or_cv_dycore: enthalpy or internal energy scaling factor for ! energy consistency real(r8), public, protected, allocatable :: cp_or_cv_dycore(:,:,:) + real(r8), public , allocatable :: te_init(:,:,:)!xxx to be removed ! ! Interfaces for public routines interface get_cp_dry @@ -137,8 +159,29 @@ module air_composition CONTAINS - ! Read namelist variables. + subroutine air_composition_register(compute_enthalpy_flux_in) + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + + logical, intent(in) :: compute_enthalpy_flux_in + + ! Set module variable compute_enthalpy_flux + compute_enthalpy_flux = compute_enthalpy_flux_in + if (masterproc) then + if (compute_enthalpy_flux) then + write(iulog, *) ' ' + write(iulog, *) 'CAM computes enthalpy flux and sends it to surface.' + else + write(iulog, *) 'CAM does not compute enthalpy flux. ' + end if + end if + + end subroutine air_composition_register + + !=========================================================================== + subroutine air_composition_readnl(nlfile) + ! Read namelist variables. use namelist_utils, only: find_group_name use spmd_utils, only: masterproc, mpicom, masterprocid use spmd_utils, only: mpi_character, mpi_logical @@ -156,7 +199,6 @@ subroutine air_composition_readnl(nlfile) ! Variable components of dry air and water species in air namelist /air_composition_nl/ dry_air_species, water_species_in_air - namelist /air_composition_nl/ compute_enthalpy_flux !----------------------------------------------------------------------- banner = repeat('*', lsize) @@ -165,7 +207,6 @@ subroutine air_composition_readnl(nlfile) ! Read variable components of dry air and water species in air dry_air_species = (/ (' ', indx = 1, num_names_max) /) water_species_in_air = (/ (' ', indx = 1, num_names_max) /) - compute_enthalpy_flux = .false. if (masterproc) then open(newunit=unitn, file=trim(nlfile), status='old') @@ -186,9 +227,6 @@ subroutine air_composition_readnl(nlfile) len(water_species_in_air)*num_names_max, mpi_character, & masterprocid, mpicom, ierr) if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: water_species_in_air") - call mpi_bcast(compute_enthalpy_flux, 1, mpi_logical, & - masterprocid, mpicom, ierr) - if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: compute_enthalpy_flux") dry_air_species_num = 0 water_species_in_air_num = 0 @@ -226,10 +264,6 @@ subroutine air_composition_readnl(nlfile) do indx = 1, water_species_in_air_num write(iulog, *) ' ', trim(water_species_in_air(indx)) end do - if (compute_enthalpy_flux) then - write(iulog, *) ' ' - write(iulog, *) 'CAM computes enthalpy flux and sends to surface.' - end if write(iulog, *) bline write(iulog, *) banner end if @@ -239,22 +273,24 @@ end subroutine air_composition_readnl !=========================================================================== subroutine air_composition_init() + use string_utils, only: int2str use spmd_utils, only: masterproc use cam_logfile, only: iulog - use physconst, only: r_universal, cpair, rair, cpwv, rh2o, cpliq, cpice, mwdry + use physconst, only: r_universal, cpair, rair, cpwv, rh2o, cpliq, cpice, mwdry, cpwv, latice, latvap, tmelt use constituents, only: cnst_get_ind, cnst_mw use ppgrid, only: pcols, pver, begchunk, endchunk + + ! Local variables integer :: icnst, ix, isize, ierr, idx integer :: liq_num, ice_num integer :: liq_idx(water_species_in_air_num) integer :: ice_idx(water_species_in_air_num) logical :: has_liq, has_ice real(r8) :: mw - + ! character(len=*), parameter :: subname = 'composition_init' character(len=*), parameter :: errstr = subname//": failed to allocate " - ! ! define cp and R for species in species_name ! @@ -276,6 +312,7 @@ subroutine air_composition_init() real(r8), parameter :: dof3 = 6._r8 real(r8), parameter :: cv3 = 0.5_r8 * r_universal * dof3 real(r8), parameter :: cp3 = 0.5_r8 * r_universal * (2._r8 + dof3) + !----------------------------------------------------------------------- liq_num = 0 ice_num = 0 @@ -348,7 +385,7 @@ subroutine air_composition_init() if (ierr /= 0) then call endrun(errstr//"cp_or_cv_dycore") end if - + allocate(te_init(pcols,4,begchunk:endchunk), stat=ierr)!xxx to be removed thermodynamic_active_species_idx = -HUGE(1) thermodynamic_active_species_idx_dycore = -HUGE(1) thermodynamic_active_species_cp = 0.0_r8 @@ -629,11 +666,50 @@ subroutine air_composition_init() (1 + liq_num + ice_num), " (1 + liq_num + ice_num)" call endrun(subname//': water_species_in_air_num /= 1+liq_num+ice_num') end if - enthalpy_reference_state = 'ice' - if (masterproc) then - write(iulog, *) 'Enthalpy reference state : ', & - TRIM(enthalpy_reference_state) + + if (compute_enthalpy_flux) then + + ! Initialising t00's and h00's + ! N.B. latent heats should be adjusted to t00a, but unless t00a=tmelt, this will break all physics + ! physics and SE dycore make different, mutually inconsistent, + ! hard-wired assumptions on t00 and h00: + ! physics : t00=tmelt, h00(ice)=L(ice; liq, T=tmelt) + ! dynamics (SE): t00=0, h00=0 + ! As a result, any water non-conservation in the dycore results in fixer + ! increments, proportional to h00a as set below. + + ! ocean choice for enthalpy at T=0 (liquid reference phase) + t00o = tmelt + h00o = -cpliq*t00o + + ! atmo choices for enthalpy at T=0 (liquid reference phase): + t00a = tmelt + h00a = -cpliq*t00a + + ! hard-wiring here + enthalpy_reference_state = 'ice' ! TODO (mvertens): should this be a namelist variable? + if (enthalpy_reference_state == 'ice') then + h00a = -cpliq*t00a ! conserve single formula for global energy + else if (enthalpy_reference_state == 'vap') then + h00a =-((cpliq-cpwv )*t00a + latvap) + endif + + ! the following ensure that the value of atmospheric enthalpy is independent of reference state + h00a_vap = h00a + ((cpliq-cpwv )*t00a + latvap) + h00a_ice = h00a + ((cpliq-cpice)*t00a - latice) + + if (masterproc) then + write(iulog, *) ' ocean t00o: ', t00o + write(iulog, *) ' ocean h00o: ', h00o + write(iulog, *) 'atmos. enthalpy_reference_state: ', trim(enthalpy_reference_state) + write(iulog, *) ' t00a: ', t00a + write(iulog, *) ' h00a: ', h00a + write(iulog, *) ' h00a_ice: ', h00a_ice + write(iulog, *) ' h00a_vap: ', h00a_vap + endif + end if + end subroutine air_composition_init !=========================================================================== diff --git a/src/utils/cam_thermo.F90 b/src/utils/cam_thermo.F90 index f65649c4ef..565eda9553 100644 --- a/src/utils/cam_thermo.F90 +++ b/src/utils/cam_thermo.F90 @@ -31,6 +31,8 @@ module cam_thermo ! DOI: 10.1029/2017MS001257 ! https://opensky.ucar.edu/islandora/object/articles:21929 + public :: get_conserved_energy, inv_conserved_energy + ! cam_thermo_init: Initialize constituent dependent properties public :: cam_thermo_init ! cam_thermo_dry_air_update: Update dry air composition dependent properties @@ -79,6 +81,7 @@ module cam_thermo ! mixing_ratio options integer, public, parameter :: DRY_MIXING_RATIO = 1 integer, public, parameter :: MASS_MIXING_RATIO = 2 + !--------------- Variables below here are for WACCM-X --------------------- ! kmvis: molecular viscosity kg/m/s real(r8), public, protected, allocatable :: kmvis(:,:,:) @@ -285,29 +288,29 @@ subroutine cam_thermo_water_update(mmr, lchnk, ncol, vcoord, to_dry_factor) !------------------------------Arguments---------------------------------- real(r8), intent(in) :: mmr(:,:,:) ! constituents array - integer, intent(in) :: lchnk ! Chunk number - integer, intent(in) :: ncol ! number of columns - integer, intent(in) :: vcoord - real(r8), optional, intent(in) :: to_dry_factor(:,:) + integer, intent(in) :: lchnk ! Chunk number + integer, intent(in) :: ncol ! number of columns + integer, intent(in) :: vcoord + real(r8), optional, intent(in) :: to_dry_factor(:,:) ! logical :: lcp call water_composition_update(mmr, lchnk, ncol, vcoord, to_dry_factor=to_dry_factor) end subroutine cam_thermo_water_update - !=========================================================================== + !=========================================================================== - ! - !*********************************************************************** - ! - ! Compute enthalpy = cp*T*dp, where dp is pressure level thickness, - ! cp is generalized cp and T temperature - ! - ! Note: tracer is in units of m*dp_dry ("mass") - ! - !*********************************************************************** - ! - subroutine get_enthalpy_1hd(tracer_mass, temp, dp_dry, & + ! + !*********************************************************************** + ! + ! Compute enthalpy = cp*T*dp, where dp is pressure level thickness, + ! cp is generalized cp and T temperature + ! + ! Note: tracer is in units of m*dp_dry ("mass") + ! + !*********************************************************************** + ! + subroutine get_enthalpy_1hd(tracer_mass, temp, dp_dry, & enthalpy, active_species_idx_dycore) use air_composition, only: dry_air_species_num, get_cp_dry ! Dummy arguments @@ -567,7 +570,7 @@ subroutine get_sum_species_1hd(tracer, active_species_idx, & real(r8), optional, intent(in) :: dp_dry(:, :) ! sum_species: sum species real(r8), intent(out) :: sum_species(:, :) - ! factor: to moist factor + ! factor: to moist factor real(r8), optional, intent(out) :: factor(:, :) ! Local variables real(r8) :: factor_loc(SIZE(tracer, 1), SIZE(tracer, 2)) @@ -722,7 +725,7 @@ end subroutine get_dp_2hd ! compute mid-level (full level) pressure from dry pressure and water tracers ! !************************************************************************************************************************* - ! + ! subroutine get_pmid_from_dpdry_1hd(tracer, mixing_ratio, active_species_idx, dp_dry, ptop, pmid, pint, dp) real(r8), intent(in) :: tracer(:,:,:) ! tracers; quantity specified by mixing_ratio arg @@ -883,7 +886,7 @@ subroutine get_gz_from_dp_dry_ptop_temp_1hd(tracer, mixing_ratio, active_species real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2)) :: pmid_local, t_v_local, dp_local, R_dry real(r8), dimension(SIZE(tracer, 1), SIZE(tracer, 2) + 1) :: pint character(len=*), parameter :: subname = 'get_gz_from_dp_dry_ptop_temp_1hd: ' - + call get_pmid_from_dp(tracer, mixing_ratio, active_species_idx, & dp_dry, ptop, pmid_local, pint=pint, dp=dp_local) @@ -1024,7 +1027,7 @@ end subroutine get_Richardson_number_1hd ! subroutine get_ps_1hd(tracer_mass, active_species_idx, dp_dry, ps, ptop) use air_composition, only: dry_air_species_num - + real(r8), intent(in) :: tracer_mass(:,:,:) ! Tracer array (q*dp) real(r8), intent(in) :: dp_dry(:,:) ! dry pressure level thickness real(r8), intent(out) :: ps(:) ! surface pressure @@ -1571,7 +1574,7 @@ end subroutine cam_thermo_calc_kappav_2hd ! if subroutine is asked to compute "te" then the latent heat terms are ! added to the kinetic (ke), internal + geopotential (se) energy terms ! - ! subroutine assumes that enthalpy term (rho*cp*T) uses dry air heat capacity + ! subroutine assumes that enthalpy term (rho*cp*T) uses dry air heat capacity !tht: why? not true ! !*************************************************************************** ! @@ -1583,6 +1586,9 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure use air_composition, only: wv_idx use physconst, only: rga, latvap, latice + use physconst, only: cpliq, cpice, cpwv, tmelt + use air_composition, only: t00a, h00a, h00a_vap, h00a_ice + use air_composition, only: compute_enthalpy_flux ! Dummy arguments ! tracer: tracer mixing ratio @@ -1612,7 +1618,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & real(r8), intent(out), optional :: te (:) ! KE: vertically integrated kinetic energy real(r8), intent(out), optional :: ke (:) - ! SE: vertically integrated enthalpy (pressure coordinate) + ! SE: vertically integrated enthalpy (pressure coordinate) ! or internal energy (z coordinate) real(r8), intent(out), optional :: se (:) ! PO: vertically integrated PHIS term (pressure coordinate) @@ -1632,6 +1638,7 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & real(r8) :: wv_vint(SIZE(tracer, 1)) ! Vertical integral of wv real(r8) :: liq_vint(SIZE(tracer, 1)) ! Vertical integral of liq real(r8) :: ice_vint(SIZE(tracer, 1)) ! Vertical integral of ice + real(r8) :: wtot_vint(SIZE(tracer, 1))! Vertical integral of water real(r8) :: pdel(SIZE(tracer, 1),SIZE(tracer, 2)) !moist pressure level thickness real(r8) :: latsub ! latent heat of sublimation @@ -1787,22 +1794,48 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & end do end do if (present(ice)) ice = ice_vint + ! Compute vertical integrals of total water. if (present(H2O)) then H2O = wv_vint + liq_vint + ice_vint end if - ! + ! latent heat terms depend on enthalpy reference state - ! + ! note choices in physconst however, ensuring they actually + wtot_vint = wv_vint + liq_vint + ice_vint latsub = latvap + latice if (present(te)) then select case (TRIM(enthalpy_reference_state)) case('ice') te = te + (latsub * wv_vint) + (latice * liq_vint) + if (compute_enthalpy_flux) then + if (vcoord .ne. vc_moist_pressure) then + ! add t00 and h00 terms + te = te + wv_vint*(cpice-cpwv )*t00a + te = te + liq_vint*(cpice-cpliq)*t00a + te = te + wtot_vint*h00a_ice + endif + end if case('liq') te = te + (latvap * wv_vint) - (latice * ice_vint) - case('wv') + if (compute_enthalpy_flux) then + if (vcoord .ne. vc_moist_pressure) then + ! add t00 and h00 terms + te = te + wv_vint*(cpliq-cpwv )*t00a + te = te + ice_vint*(cpliq-cpice)*t00a + te = te + wtot_vint*h00a + endif + end if + case('vap') te = te - (latvap * liq_vint) - (latsub * ice_vint) + if (compute_enthalpy_flux) then + if (vcoord .ne. vc_moist_pressure) then + ! add t00 and h00 terms + te = te + liq_vint*(cpwv -cpliq)*t00a + te = te + ice_vint*(cpwv -cpice)*t00a + te = te + wtot_vint*h00a_vap + endif + end if case default write(iulog, *) subname, ' enthalpy reference state not ', & 'supported: ', TRIM(enthalpy_reference_state) @@ -1812,4 +1845,566 @@ subroutine get_hydrostatic_energy_1hd(tracer, moist_mixing_ratio, pdel_in, & deallocate(species_idx, species_liq_idx, species_ice_idx) end subroutine get_hydrostatic_energy_1hd + !=========================================================================== + + subroutine get_conserved_energy(moist_mixing_ratio, ktop, kbot & + , cp_or_cv, T, tracer, pdel_in & + , pdel, te & + , qini, liqini, iceini & + , phis & + , gph & + , U, V, W, rairv & + , flatent,latent,potential,kinetic,temce & + , refstate, vcoord, dycore_idx) + + use dycore, only: dycore_is + use cam_logfile, only: iulog + use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + use air_composition, only: wv_idx + use physconst, only: rga, latvap, latice + use physconst, only: cpliq, cpice, cpwv, tmelt + use air_composition, only: t00a, h00a, h00a_vap, h00a_ice + + ! arguments in: + ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry + logical , intent(in) :: moist_mixing_ratio + integer , intent(in) :: ktop, kbot + ! cp_or_cv: dry air heat capacity under constant pressure or + ! constant volume (depends on vcoord) + real(r8), intent(in) :: cp_or_cv(:,:) + real(r8), intent(in) :: T(:,:) + real(r8), intent(in) :: tracer(:,:,:) + ! pdel: pressure level thickness + real(r8), intent(in) :: pdel_in(:,:) !N.B. this should be g*\rho*dz for MPAS + + ! arguments out: + ! conserved total energy/enthalpy per unit mass + real(r8), intent(out) :: te (:,:) + ! pdel: layer mass + real(r8), intent(out) :: pdel(:,:) !N.B. this should be g*\rho*dz for MPAS + + ! arguments optional: + real(r8), intent(in), optional :: qini(:,:), liqini(:,:), iceini(:,:) + ! surface geopotential -- should be made mandatory arg + real(r8), intent(in), optional :: phis(:) + ! geopotential height, required for MPAS: te=u_m:=c_v*T+latent+gz+KE + ! dycore_is('MPAS') and gph not present -> stop + real(r8), intent(in), optional :: gph(:,:) + ! N.B. either PHIS or GPH must be present + ! horizontal winds --> add KE (should be made mandatory arguments) + real(r8), intent(in), optional :: U(:,:) + real(r8), intent(in), optional :: V(:,:) + ! vertical wind --> add to KE (non-hydrostatic) + real(r8), intent(in), optional :: W(:,:) + real(r8), intent(in), optional :: Rairv(:,:) + character(len=3),intent(in),optional :: refstate + integer, intent(in), optional :: vcoord ! vertical coordinate + ! dycore_idx: use dycore index for thermodynamic active species + logical, intent(in) , optional :: dycore_idx + real(r8), intent(out), optional :: flatent(:,:) + real(r8), intent(out), optional :: latent(:,:) + real(r8), intent(out), optional :: potential(:,:) + real(r8), intent(out), optional :: kinetic(:,:) + real(r8), intent(out), optional :: temce(:,:) ! Total Enthalpy Minus Conserved Energy + + ! Local variables + real(r8) :: qwv (SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qliq(SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qice(SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qtot(SIZE(tracer, 1),SIZE(tracer, 2)), latsub + real(r8) :: work(SIZE(tracer, 1),SIZE(tracer, 2)) + + integer :: ierr + integer :: kdx, idx, nkd, nid ! coord indices + integer :: qdx ! tracer index + integer :: wvidx ! water vapor index + integer, allocatable :: species_idx(:) + integer, allocatable :: species_liq_idx(:) + integer, allocatable :: species_ice_idx(:) + character(len=3) :: loc_refstate + character(len=*), parameter :: subname = 'get_conserved_energy' + + allocate(species_idx(thermodynamic_active_species_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_idx array') + end if + allocate(species_liq_idx(thermodynamic_active_species_liq_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_liq_idx array') + end if + allocate(species_ice_idx(thermodynamic_active_species_ice_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_ice_idx array') + end if + + nkd=SIZE(tracer, 2) + nid=SIZE(tracer, 1) + + if(present(refstate))then + loc_refstate=trim(refstate) + else + loc_refstate=trim(enthalpy_reference_state) + endif + + if (present(dycore_idx))then + if (dycore_idx) then + species_idx(:) = thermodynamic_active_species_idx_dycore(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:) + else + species_idx(:) = thermodynamic_active_species_idx(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) + end if + else + species_idx(:) = thermodynamic_active_species_idx(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) + end if + + if (moist_mixing_ratio) then + pdel = pdel_in*rga + else + pdel = pdel_in*rga + if (present(qini).and.present(liqini).and.present(iceini))then + pdel(:,:) = pdel(:,:) + pdel_in(:, :)*(qini(:,:)+liqini(:,:)+iceini(:,:))*rga + else + do qdx = dry_air_species_num+1, thermodynamic_active_species_num + pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx))*rga + end do + endif + end if + + do kdx = ktop, kbot + do idx = 1, nid + te(idx,kdx) = T(idx,kdx)*cp_or_cv(idx, kdx) + end do + end do + + work(:,:)=0._r8 + if(present(phis))then + do kdx = ktop, kbot + do idx = 1, nid + work(idx,kdx) = phis(idx) + end do + end do + endif + if(dycore_is('MPAS')) then + if(.not.present(gph)) call endrun(subname//': conserved_energy function'// & + ' requires GPH in input for non-hydrostatic case') + do kdx = ktop, kbot + do idx = 1, nid + work(idx,kdx) = work(idx,kdx) + gph(idx,kdx)/rga + end do + end do + endif + if (present(potential)) then + do kdx = ktop, kbot + do idx = 1, nid + potential(idx,kdx) = work(idx,kdx) + end do + end do + else + do kdx = ktop, kbot + do idx = 1, nid + te(idx,kdx) = te(idx,kdx) + work(idx,kdx) + end do + end do + endif + + if(present(qini).and.present(liqini).and.present(iceini))then + qwv (:,:)=qini (:,:) + qliq(:,:)=liqini(:,:) + qice(:,:)=iceini(:,:) + else + qwv (:,:) = tracer(:,:,wv_idx) + qliq(:,:) = 0._r8 + do qdx = 1, thermodynamic_active_species_liq_num + qliq(:,:) = qliq(:,:) + tracer(:,:,species_liq_idx(qdx)) + enddo + qice(:,:) = 0._r8 + do qdx = 1, thermodynamic_active_species_ice_num + qice(:,:) = qice(:,:) + tracer(:,:,species_ice_idx(qdx)) + enddo + endif + + latsub = latvap + latice + select case (TRIM(loc_refstate)) + case('ice') + work(:,:) = (latsub * qwv ) + (latice * qliq) + case('liq') + work(:,:) = (latvap * qwv ) - (latice * qice) + case('vap') + work(:,:) =-(latvap * qliq) - (latsub * qice) + case default + write(iulog, *) subname, ' enthalpy reference state not ', & + 'supported: ', TRIM(loc_refstate) + call endrun(subname//': enthalpy reference state not supported') + end select + if (present(latent).or.present(flatent)) then + if (present(flatent)) then + do kdx = ktop, kbot + do idx = 1, nid + flatent(idx,kdx) = work(idx,kdx) + end do + end do + endif + if (present(latent)) then + do kdx = ktop, kbot + do idx = 1, nid + latent(idx,kdx) = work(idx,kdx) + end do + end do + endif + else + do kdx = ktop, kbot + do idx = 1, nid + te(idx,kdx) = te(idx,kdx) + work(idx,kdx) + end do + end do + endif + + ! add t00 and h00 terms + if(present(vcoord))then + if(vcoord.ne.vc_moist_pressure) then + qtot(:,:) = qice(:,:) + qliq(:,:) + qwv (:,:) + select case (TRIM(loc_refstate)) + case('ice') + work(:,:) = qwv (:,:)*(cpice-cpwv )*t00a & + + qliq(:,:)*(cpice-cpliq)*t00a & + + qtot(:,:)*h00a_ice + case('liq') + work(:,:) = qwv (:,:)*(cpliq-cpwv )*t00a & + + qice(:,:)*(cpliq-cpice)*t00a & + + qtot(:,:)*h00a + case('vap') + work(:,:) = qliq(:,:)*(cpwv -cpliq)*t00a & + + qice(:,:)*(cpwv -cpice)*t00a & + + qtot(:,:)*h00a_vap + end select + if (present(latent)) then + do kdx = ktop, kbot + do idx = 1, nid + latent(idx,kdx) = latent(idx,kdx)+work(idx,kdx) + end do + end do + else + do kdx = ktop, kbot + do idx = 1, nid + te(idx,kdx) = te(idx,kdx) + work(idx,kdx) + end do + end do + endif + endif + endif + + if(present(U).and.present(V)) then + do kdx = ktop, kbot + do idx = 1, nid + work(idx,kdx) = .5_r8*(u(idx,kdx)**2+v(idx,kdx)**2) + enddo + enddo + if (present(kinetic)) then + do kdx = ktop, kbot + do idx = 1, nid + kinetic(idx,kdx)= work(idx,kdx) + end do + end do + else + do kdx = ktop, kbot + do idx = 1, nid + te(idx,kdx) = te(idx,kdx) + work(idx,kdx) + end do + end do + endif + endif + + if(present(temce)) then + if(dycore_is('MPAS'))then + if(.not.(present(rairv))) call endrun(subname//': TEMCE required but'// & + ' Rairv not provided in non-hydrostatic case') + do kdx = ktop, kbot + do idx = 1, nid + temce(idx,kdx) = T(idx,kdx)*rairv(idx, kdx) + end do + end do + else + if(.not.(present(gph))) call endrun(subname//': TEMCE required but'// & + ' GPH not provided in hydrostatic case') + do kdx = ktop, kbot + do idx = 1, nid + temce(idx,kdx) = gph(idx,kdx)/rga + end do + end do + endif + endif + + deallocate(species_idx, species_liq_idx, species_ice_idx) + + end subroutine get_conserved_energy + + !=========================================================================== + + subroutine inv_conserved_energy(moist_mixing_ratio & + , ktop, kbot & + , te, cp_or_cv, tracer, pdel_in & + , pdel, T & + , phis & + , gph & + , U, V, W & + , flatent,latent,potential,kinetic & + , refstate, vcoord, dycore_idx) + + use cam_logfile, only: iulog + use dycore, only: dycore_is + use dyn_tests_utils, only: vc_height, vc_moist_pressure, vc_dry_pressure + use air_composition, only: wv_idx + use physconst, only: rga, latvap, latice + use physconst, only: cpliq, cpice, cpwv, tmelt + use air_composition, only: t00a, h00a, h00a_vap, h00a_ice + + ! arguments in: + ! note - if pdeldry passed to subroutine then tracer mixing ratio must be dry + logical , intent(in) :: moist_mixing_ratio + integer , intent(in) :: ktop, kbot + ! conserved energy/enthalpy + real(r8), intent(in) :: te(:,:) + ! cp_or_cv: dry air heat capacity under constant pressure or + ! constant volume (depends on vcoord) + real(r8), intent(in) :: cp_or_cv(:,:) + real(r8), intent(in) :: tracer(:,:,:) + ! pdel: pressure level thickness + real(r8), intent(in) :: pdel_in(:,:) !N.B. this should be g*\rho*dz for MPAS + + ! arguments out: + ! temperature + real(r8), intent(out) :: T(:,:) + ! pdel: layer mass + real(r8), intent(out) :: pdel(:,:) !N.B. this should be g*\rho*dz for MPAS + + ! arguments optional: + ! surface geopotential --> compute te=e_m:=c_p*T+latent+phis+KE (hydrostatic) + real(r8), intent(in), optional :: phis(:) + ! geopotential height --> compute te=u_m:=c_v*T+latent+gz+KE (MPAS) + ! should be =z_mid in output os subroutine geopotential_t + real(r8), intent(in), optional :: gph(:,:) + character(len=3),intent(in),optional :: refstate + integer, intent(in), optional :: vcoord ! vertical coordinate + !N.B. either PHIS or GPH must be present + ! dycore_idx: use dycore index for thermodynamic active species + logical, intent(in), optional :: dycore_idx + ! horizontal winds --> add KE (will be made mandatory arguments later) + real(r8), intent(in), optional :: U(:,:) + real(r8), intent(in), optional :: V(:,:) + ! vertical wind --> add to KE (MPAS) + real(r8), intent(in), optional :: W(:,:) + real(r8), intent(in), optional :: flatent(:,:) + real(r8), intent(in), optional :: latent(:,:) + real(r8), intent(in), optional :: potential(:,:) + real(r8), intent(in), optional :: kinetic(:,:) + + ! Local variables + real(r8) ::tetmp(SIZE(tracer, 1),SIZE(tracer, 2)) + real(r8) :: qwv (SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qliq(SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qice(SIZE(tracer, 1),SIZE(tracer, 2)) & + ,qtot(SIZE(tracer, 1),SIZE(tracer, 2)), latsub + + integer :: ierr + integer :: kdx, idx, nkd, nid ! coord indices + integer :: qdx ! tracer index + integer :: wvidx ! water vapor index + integer, allocatable :: species_idx(:) + integer, allocatable :: species_liq_idx(:) + integer, allocatable :: species_ice_idx(:) + character(len=3) :: loc_refstate + character(len=*), parameter :: subname = 'get_conserved_energy' + + allocate(species_idx(thermodynamic_active_species_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_idx array') + end if + allocate(species_liq_idx(thermodynamic_active_species_liq_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_liq_idx array') + end if + allocate(species_ice_idx(thermodynamic_active_species_ice_num), stat=ierr) + if ( ierr /= 0 ) then + call endrun(subname//': allocation error for species_ice_idx array') + end if + + nkd=SIZE(tracer, 2) + nid=SIZE(tracer, 1) + + if(present(refstate))then + loc_refstate=trim(refstate) + else + loc_refstate=trim(enthalpy_reference_state) + endif + + if (present(dycore_idx))then + if (dycore_idx) then + species_idx(:) = thermodynamic_active_species_idx_dycore(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx_dycore(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx_dycore(:) + else + species_idx(:) = thermodynamic_active_species_idx(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) + end if + else + species_idx(:) = thermodynamic_active_species_idx(:) + species_liq_idx(:) = thermodynamic_active_species_liq_idx(:) + species_ice_idx(:) = thermodynamic_active_species_ice_idx(:) + end if + + if (moist_mixing_ratio) then + pdel = pdel_in*rga + else + pdel = pdel_in*rga + do qdx = dry_air_species_num+1, thermodynamic_active_species_num + pdel(:,:) = pdel(:,:) + pdel_in(:, :)*tracer(:,:,species_idx(qdx))*rga + end do + end if + + if(present(kinetic)) then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = te(idx,kdx) - kinetic(idx,kdx) + enddo + enddo + else if(present(U).and.present(V)) then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = te(idx,kdx) - .5_r8*(u(idx,kdx)**2+v(idx,kdx)**2) + enddo + enddo + else + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = te(idx,kdx) + end do + end do + endif + + if(present(potential)) then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = tetmp(idx,kdx) - potential(idx,kdx) + end do + end do + else + if(present(phis))then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = tetmp(idx,kdx) - phis(idx) + end do + end do + endif + if(dycore_is('MPAS')) then + if(.not.present(gph)) call endrun(subname//': conserved_energy function'// & + ' requires GPH in input for non-hydrostatic case') + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = tetmp(idx,kdx) - gph(idx,kdx)/rga + end do + end do + endif + endif + + if (present(latent)) then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = tetmp(idx,kdx) - latent(idx,kdx) + end do + end do + else + qwv (:,:) = tracer(:,:,wv_idx) + qliq(:,:) = 0._r8 + do qdx = 1, thermodynamic_active_species_liq_num + qliq(:,:) = qliq(:,:) + tracer(:,:,species_liq_idx(qdx)) + enddo + qice(:,:) = 0._r8 + do qdx = 1, thermodynamic_active_species_ice_num + qice(:,:) = qice(:,:) + tracer(:,:,species_ice_idx(qdx)) + enddo + qtot(:,:) = qice(:,:) + qliq(:,:) + qwv (:,:) + if (present(flatent)) then + do kdx = ktop, kbot + do idx = 1, nid + tetmp(idx,kdx) = tetmp(idx,kdx) - flatent(idx,kdx) + end do + end do + if(present(vcoord))then + if(vcoord.ne.vc_moist_pressure) then + ! add t00 and h00 terms + select case (TRIM(loc_refstate)) + case('ice') + tetmp(:,:) = tetmp(:,:) -(qwv (:,:)*(cpice-cpwv )*t00a & + +qliq(:,:)*(cpice-cpliq)*t00a & + +qtot(:,:)*h00a_ice ) + case('liq') + tetmp(:,:) = tetmp(:,:) -(qwv (:,:)*(cpliq-cpwv )*t00a & + +qice(:,:)*(cpliq-cpice)*t00a & + +qtot(:,:)*h00a ) + case('vap') + tetmp(:,:) = tetmp(:,:) -(qliq(:,:)*(cpwv -cpliq)*t00a & + +qice(:,:)*(cpwv -cpice)*t00a & + +qtot(:,:)*h00a_vap ) + case default + write(iulog, *) subname, ' enthalpy reference state not ', & + 'supported: ', TRIM(loc_refstate) + call endrun(subname//': enthalpy reference state not supported') + end select + endif + endif + else + latsub = latvap + latice + select case (TRIM(loc_refstate)) + case('ice') + tetmp(:,:) = tetmp(:,:) - (latsub * qwv ) - (latice * qliq) + if(present(vcoord))then + if(vcoord.ne.vc_moist_pressure) then + tetmp(:,:) = tetmp(:,:) -(qwv (:,:)*(cpice-cpwv )*t00a & + +qliq(:,:)*(cpice-cpliq)*t00a & + +qtot(:,:)*h00a_ice ) + endif + endif + case('liq') + tetmp(:,:) = tetmp(:,:) - (latvap * qwv ) + (latice * qice) + if(present(vcoord))then + if(vcoord.ne.vc_moist_pressure) then + tetmp(:,:) = tetmp(:,:) -(qwv (:,:)*(cpliq-cpwv )*t00a & + +qice(:,:)*(cpliq-cpice)*t00a & + +qtot(:,:)*h00a ) + endif + endif + case('vap') + tetmp(:,:) = tetmp(:,:) + (latvap * qliq) + (latsub * qice) + if(present(vcoord))then + if(vcoord.ne.vc_moist_pressure) then + tetmp(:,:) = tetmp(:,:) -(qliq(:,:)*(cpwv -cpliq)*t00a & + +qice(:,:)*(cpwv -cpice)*t00a & + +qtot(:,:)*h00a_vap ) + endif + endif + case default + write(iulog, *) subname, ' enthalpy reference state not ', & + 'supported: ', TRIM(loc_refstate) + call endrun(subname//': enthalpy reference state not supported') + end select + endif + endif + + do kdx = ktop, kbot + do idx = 1, nid + T(idx,kdx) = tetmp(idx,kdx)/cp_or_cv(idx, kdx) + end do + end do + + deallocate(species_idx, species_liq_idx, species_ice_idx) + + end subroutine inv_conserved_energy + +!------------------------------------------------------------------------------- end module cam_thermo