From 11d061cdbbcc4dea9f4c7f5ea7425c2d42f461a8 Mon Sep 17 00:00:00 2001 From: cacraigucar Date: Wed, 28 Jun 2023 17:09:41 -0600 Subject: [PATCH] Merge pull request #537 from cacraigucar/frierson_v2 cam6_3_117: Bring in Frierson Simple Model ESCOMP commit: 6019d2de13e9eb5c82b233da18f5a2286f374313 --- bld/build-namelist | 20 +- .../use_cases/dctest_frierson.xml | 31 + doc/ChangeLog | 78 ++ src/control/cam_control_mod.F90 | 8 +- src/control/runtime_opts.F90 | 6 + src/physics/simple/frierson.F90 | 1174 +++++++++++++++++ src/physics/simple/frierson_cam.F90 | 1046 +++++++++++++++ src/physics/simple/physpkg.F90 | 89 +- src/physics/simple/restart_physics.F90 | 14 + 9 files changed, 2456 insertions(+), 10 deletions(-) create mode 100644 bld/namelist_files/use_cases/dctest_frierson.xml create mode 100644 src/physics/simple/frierson.F90 create mode 100644 src/physics/simple/frierson_cam.F90 diff --git a/bld/build-namelist b/bld/build-namelist index 0de0ead4cb..907886d28c 100755 --- a/bld/build-namelist +++ b/bld/build-namelist @@ -473,10 +473,26 @@ if ($phys eq 'adiabatic') { ++$phys_mode_flags; } my $ideal_mode = 0; -if ($phys eq 'kessler' or $phys eq 'held_suarez' or $phys eq 'tj2016') { +if ($phys eq 'kessler' or $phys eq 'held_suarez' or $phys eq 'tj2016' or $phys eq 'grayrad' ) { $ideal_mode = 1; ++$phys_mode_flags; } +if ($phys eq 'grayrad' ) { + add_default($nl, 'frierson_albedo'); + add_default($nl, 'frierson_c0'); + add_default($nl, 'frierson_deltas'); + add_default($nl, 'frierson_fb'); + add_default($nl, 'frierson_linfrac'); + add_default($nl, 'frierson_ri_c'); + add_default($nl, 'frierson_tau_eqtr'); + add_default($nl, 'frierson_tau_pole'); + add_default($nl, 'frierson_tdlt'); + add_default($nl, 'frierson_tmin'); + add_default($nl, 'frierson_twidth'); + add_default($nl, 'frierson_wetdrycoef'); + add_default($nl, 'frierson_wind_min'); + add_default($nl, 'frierson_z0'); +} if ($phys_mode_flags > 1) { die "$ProgName - ERROR: Only one of the variables atm_adiabatic, atm_ideal_phys, and aqua_planet can be set .true. \n"; } @@ -4858,6 +4874,8 @@ sub check_snapshot_settings { push (@validList_bc, ("'kessler_tend'")); } elsif ($phys eq 'tj2016') { push (@validList_bc, ("'thatcher_jablonowski_precip_tend'")); + } elsif ($phys eq 'grayrad') { + push (@validList_bc, ("'frierson_tend'")); } if ($chem ne 'none') { push (@validList_bc, ("'chem_timestep_tend'")); diff --git a/bld/namelist_files/use_cases/dctest_frierson.xml b/bld/namelist_files/use_cases/dctest_frierson.xml new file mode 100644 index 0000000000..340052bc84 --- /dev/null +++ b/bld/namelist_files/use_cases/dctest_frierson.xml @@ -0,0 +1,31 @@ + + + + + +0. +0. +0. +fixed_parameters + + +.false. + + +.true. +-720 +'A' + + + 'U:A','T:A','V:A','Q:A','Z3:A','PRECL:A','PRECC:A','PS:A','SST:A','TS:A', 'gray_DTCOND', 'gray_DQCOND', 'gray_EVAPDT', 'gray_EVAPDQ', 'gray_PRECL' , 'gray_PRECC' , 'gray_Tsurf' , 'gray_Qsurf' , 'gray_Cdrag' , 'gray_Zpbl' , 'gray_KVH' , 'gray_KVM' , 'gray_VSE' , 'gray_Zm' , 'gray_Rf' , 'gray_DTV' , 'gray_DUV' , 'gray_DVV' , 'gray_VD01' , 'gray_SHflux', 'gray_LHflux', 'gray_TauU' , 'gray_TauV' , 'gray_QRL' , 'gray_QRS' , 'gray_SWflux', 'gray_LUflux', 'gray_LDflux', 'gray_LWflux', 'gray_LUflux_TOA', 'gray_LDflux_TOA', 'gray_LWflux_TOA' + + +'moist_baroclinic_wave_dcmip2016' + +atm/cam/inic/fv/FGRAYRAD_f19.cam.i.0051-01-01-00000_c20230510.nc + + + 'TT_SLOT','TT_GBALL','TT_TANH','TT_EM8','TT_Y2_2','TT_Y32_16' + + + diff --git a/doc/ChangeLog b/doc/ChangeLog index 8462a0d7c1..4131d5c1ba 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,6 +1,84 @@ =============================================================== +Tag name: cam6_3_117 +Originator(s): cacraig, patc, islas +Date: June 28, 2023 +One-line Summary: Bring in Frierson Simple Model (gray radiation) +Github PR URL: https://github.com/ESCOMP/CAM/pull/537 + +Purpose of changes (include the issue number and title text for each relevant GitHub issue): + - Frierson gray radiation: https://github.com/ESCOMP/CAM/issues/8 + +Describe any changes made to build system: + - Add physics option "grayrad" + +Describe any changes made to the namelist: + - Introduced namelist settings for Frierson Gray radiation (all prepended with "frierson") + +List any changes to the defaults for the boundary datasets: N/A + +Describe any substantial timing or memory changes: N/A + +Code reviewed by: fvitt, peverwhee, nusbaume + +List all files eliminated: N/A + +List all files added and what they do: +A bld/namelist_files/use_cases/dctest_frierson.xml + - Add use_case for Gray radiation model + +A src/physics/simple/frierson.F90 +A src/physics/simple/frierson_cam.F90 + - Bring in Frierson code along with the CAM interface code to it + +List all existing files that have been modified, and describe the changes: +M bld/build-namelist + - introduce grayrad physics option and accompanying namelist values + +M bld/config_files/definition.xml +M bld/configure + - introduce grayrad physics option + +M bld/namelist_files/namelist_defaults_cam.xml +M bld/namelist_files/namelist_definition.xml + - introduce grayrad namelist values + +M cime_config/config_component.xml + - introduce %GRAYRAD compset qualifier + +M cime_config/config_compsets.xml + - introduce FGRAYRAD compset + +M cime_config/testdefs/testlist_cam.xml + - introduce prebeta test for FGRAYRAD compset + +M src/control/cam_control_mod.F90 +M src/control/runtime_opts.F90 +M src/physics/simple/physpkg.F90 +M src/physics/simple/restart_physics.F90 + - mods to support gray radiation + +If there were any failures reported from running test_driver.sh on any test +platform, and checkin with these failures has been OK'd by the gatekeeper, +then copy the lines from the td.*.status files for the failed tests to the +appropriate machine below. All failed tests must be justified. + +cheyenne/intel/aux_cam: all BFB except: + ERP_Ln9_Vnuopc.ne30pg3_ne30pg3_mg17.FW2000climo.cheyenne_intel.cam-outfrq9s_wcm_ne30 (Overall: FAIL) details: + - pre-existing failure + +izumi/nag/aux_cam: all BFB except: + DAE_Vnuopc.f45_f45_mg37.FHS94.izumi_nag.cam-dae (Overall: FAIL) details: + - pre-existing failure + +izumi/gnu/aux_cam: all BFB + +Isla Simpson analyzed the results from these code mods and approved the release of the Gray Radiation option. + +=============================================================== +=============================================================== + Tag name:cam6_3_116 Originator(s): cacraig, hannay, jedwards, eaton Date: June 23, 2023 diff --git a/src/control/cam_control_mod.F90 b/src/control/cam_control_mod.F90 index ce6b3deaad..3d954f68ce 100644 --- a/src/control/cam_control_mod.F90 +++ b/src/control/cam_control_mod.F90 @@ -33,8 +33,9 @@ module cam_control_mod logical, protected :: ideal_phys ! true => run Held-Suarez (1994) physics logical, protected :: kessler_phys ! true => run Kessler physics logical, protected :: tj2016_phys ! true => run tj2016 physics +logical, protected :: frierson_phys ! true => run frierson physics logical, protected :: simple_phys ! true => adiabatic or ideal_phys or kessler_phys - ! or tj2016 + ! or tj2016 or frierson logical, protected :: aqua_planet ! Flag to run model in "aqua planet" mode logical, protected :: moist_physics ! true => moist physics enabled, i.e., ! (.not. ideal_phys) .and. (.not. adiabatic) @@ -135,8 +136,9 @@ subroutine cam_ctrl_set_physics_type(phys_package) ideal_phys = trim(phys_package) == 'held_suarez' kessler_phys = trim(phys_package) == 'kessler' tj2016_phys = trim(phys_package) == 'tj2016' + frierson_phys = trim(phys_package) == 'grayrad' - simple_phys = adiabatic .or. ideal_phys .or. kessler_phys .or. tj2016_phys + simple_phys = adiabatic .or. ideal_phys .or. kessler_phys .or. tj2016_phys .or. frierson_phys moist_physics = .not. (adiabatic .or. ideal_phys) @@ -154,6 +156,8 @@ subroutine cam_ctrl_set_physics_type(phys_package) write(iulog,*) 'Run model with Kessler warm-rain physics forcing' else if (tj2016_phys) then write(iulog,*) 'Run model with Thatcher-Jablonowski (2016) physics forcing (moist Held-Suarez)' + else if (frierson_phys) then + write(iulog,*) 'Run model with Frierson (2006) physics' end if end if diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index f09554244d..2d707c4839 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -91,6 +91,9 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use rate_diags, only: rate_diags_readnl use tracers, only: tracers_readnl use nudging, only: nudging_readnl +#if ( defined SIMPLE ) + use frierson_cam, only: frierson_readnl +#endif use dyn_comp, only: dyn_readnl use ionosphere_interface,only: ionosphere_readnl @@ -195,6 +198,9 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call rate_diags_readnl(nlfilename) call scam_readnl(nlfilename, single_column, scmlat, scmlon) call nudging_readnl(nlfilename) +#if ( defined SIMPLE ) + call frierson_readnl(nlfilename) +#endif call dyn_readnl(nlfilename) call ionosphere_readnl(nlfilename) diff --git a/src/physics/simple/frierson.F90 b/src/physics/simple/frierson.F90 new file mode 100644 index 0000000000..08e524bbdc --- /dev/null +++ b/src/physics/simple/frierson.F90 @@ -0,0 +1,1174 @@ +module frierson +!------------------------------------------------------------------------------------ +! +! Purpose: Implement idealized forcings described in +! Frierson, et al. (2006), "A Gray-Radiation +! Aquaplanet Moist GCM. Part I. Static Stability and Eddy Scale" +! J. Atmos. Sci., Vol. 63, 2548-2566. +! +! DOI: https://doi.org/10.1175/JAS3753.1 +! +!==================================================================================== + ! + ! The only modules that are permitted + !-------------------------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use shr_const_mod, only: pi => shr_const_pi + + ! Set all Global values and routine to private by default + ! and then explicitly set their exposure + !--------------------------------------------------------- + implicit none + private + save + + public:: frierson_set_const + public:: frierson_condensate_NONE + public:: frierson_condensate + public:: frierson_condensate_TJ16 + public:: frierson_condensate_USER + public:: frierson_pbl + public:: frierson_pbl_USER + public:: frierson_radiation + public:: frierson_radiation_USER + + ! Global Tuning Parameters: + ! T0 and E0 are the temperature and saturation vapor pressure used + ! to calculate qsat values, the saturation value for Q (kg/kg) + !-------------------------------------------------------------------- + real(r8):: T0 + real(r8):: E0 + real(r8):: Erad + real(r8):: Wind_min + real(r8):: Z0 + real(r8):: Ri_c + real(r8):: Karman + real(r8):: Fb + real(r8):: Rs0 + real(r8):: DeltaS + real(r8):: Tau_eqtr + real(r8):: Tau_pole + real(r8):: LinFrac + real(r8):: Boltz + real(r8):: C_ocn + + ! Private data + !---------------------- + real(r8),private :: gravit ! g: gravitational acceleration (m/s2) + real(r8),private :: cappa ! Rd/cp + real(r8),private :: rair ! Rd: dry air gas constant (J/K/kg) + real(r8),private :: cpair ! cp: specific heat of dry air (J/K/kg) + real(r8),private :: latvap ! L: latent heat of vaporization (J/kg) + real(r8),private :: rh2o ! Rv: water vapor gas constant (J/K/kg) + real(r8),private :: epsilo ! Rd/Rv: ratio of h2o to dry air molecular weights + real(r8),private :: rhoh2o ! density of liquid water (kg/m3) + real(r8),private :: zvir ! (rh2o/rair) - 1, needed for virtual temperature + real(r8),private :: ps0 ! Base state surface pressure (Pa) + + real(r8),private :: latvap_div_cpair ! latvap/cpair + real(r8),private :: latvap_div_rh2o ! latvap/rh2o + + real(r8),private,allocatable:: etamid(:) ! hybrid coordinate - midpoints + + +contains + !======================================================================= + subroutine frierson_set_const(I_gravit,I_cappa ,I_rair ,I_cpair ,I_latvap , & + I_rh2o ,I_epsilo ,I_rhoh2o ,I_zvir ,I_ps0 , & + I_etamid,I_T0 ,I_E0 ,I_Erad ,I_Wind_min, & + I_Z0 ,I_Ri_c ,I_Karman ,I_Fb ,I_Rs0 , & + I_DeltaS,I_Tau_eqtr,I_Tau_pole,I_LinFrac,I_Boltz , & + I_Cocn ) + ! + ! frierson_set_const: Set parameters and constants for the Frierson + ! Model fomulation. Optional inputs can be provided + ! to over-ride the model defaults. + !===================================================================== + + use cam_abortutils, only: handle_allocate_error + + ! + ! Passed Variables + !------------------- + real(r8),intent(in):: I_gravit + real(r8),intent(in):: I_cappa + real(r8),intent(in):: I_rair + real(r8),intent(in):: I_cpair + real(r8),intent(in):: I_latvap + real(r8),intent(in):: I_rh2o + real(r8),intent(in):: I_epsilo + real(r8),intent(in):: I_rhoh2o + real(r8),intent(in):: I_zvir + real(r8),intent(in):: I_ps0 + real(r8),intent(in):: I_etamid(:) + + real(r8),intent(in) :: I_T0 + real(r8),intent(in) :: I_E0 + real(r8),intent(in) :: I_Erad + real(r8),intent(in) :: I_Wind_min + real(r8),intent(in) :: I_Z0 + real(r8),intent(in) :: I_Ri_c + real(r8),intent(in) :: I_Karman + real(r8),intent(in) :: I_Fb + real(r8),intent(in) :: I_Rs0 + real(r8),intent(in) :: I_DeltaS + real(r8),intent(in) :: I_Tau_eqtr + real(r8),intent(in) :: I_Tau_pole + real(r8),intent(in) :: I_LinFrac + real(r8),intent(in) :: I_Boltz + real(r8),intent(in) :: I_Cocn + + integer :: ierr + + ! Set global constants for later use + !------------------------------------ + gravit = I_gravit + cappa = I_cappa + rair = I_rair + cpair = I_cpair + latvap = I_latvap + rh2o = I_rh2o + epsilo = I_epsilo + rhoh2o = I_rhoh2o + zvir = I_zvir + ps0 = I_ps0 + T0 = I_T0 + E0 = I_E0 + Erad = I_Erad + Wind_min = I_Wind_min + Z0 = I_Z0 + Ri_c = I_Ri_c + Karman = I_Karman + Fb = I_Fb + Rs0 = I_Rs0 + DeltaS = I_DeltaS + Tau_eqtr = I_Tau_eqtr + Tau_pole = I_Tau_pole + LinFrac = I_LinFrac + Boltz = I_Boltz + C_ocn = I_Cocn + + latvap_div_cpair = latvap/cpair + latvap_div_rh2o = latvap/rh2o + + ! allocate space and set the level information + !---------------------------------------------- + allocate(etamid(size(I_etamid)),stat=ierr) + if (ierr /= 0) then + call handle_allocate_error(ierr, 'frierson_set_const', 'etamid') + end if + + etamid = I_etamid + + end subroutine frierson_set_const + !======================================================================= + + + !======================================================================= + subroutine frierson_condensate_NONE(ncol,pver,pmid,T,qv,relhum,precl) + ! + ! Precip_process: Implement NO large-scale condensation/precipitation + !======================================================================= + ! + ! Passed Variables + !--------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(inout):: qv (ncol,pver) ! specific humidity Q (kg/kg) + real(r8),intent(out) :: relhum(ncol,pver) ! relative humidity + real(r8),intent(out) :: precl (ncol) ! large-scale precipitation rate (m/s) + ! + ! Local Values + !------------- + real(r8):: qsat + integer :: i, k + + ! Set large scale precipitation rates to zero + !-------------------------------------------------------------------------- + precl(:) = 0.0_r8 + + ! Large-Scale Condensation and Precipitation without cloud stage + !--------------------------------------------------------------- + do k = 1, pver + do i = 1, ncol + ! calculate saturation value for Q + !---------------------------------- + qsat = epsilo*E0/pmid(i,k)*exp(-latvap_div_rh2o*((1._r8/T(i,k))-1._r8/T0)) + + ! Set percent relative humidity + !------------------------------- + relhum(i,k) = (qv(i,k)/qsat)*100._r8 + end do + end do + + end subroutine frierson_condensate_NONE + !======================================================================= + + + !======================================================================= + subroutine frierson_condensate(ncol,pver,dtime,pmid,pdel,T,qv,relhum,precl,evapdt,evapdq) + ! + ! Precip_process: Implement large-scale condensation/precipitation + ! from Frierson 2006. + ! + !======================================================================= + ! + ! Passed Variables + !--------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: dtime ! time step (s) + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(in) :: pdel (ncol,pver) ! layer thickness (Pa) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(inout):: qv (ncol,pver) ! specific humidity Q (kg/kg) + real(r8),intent(out) :: relhum(ncol,pver) ! relative humidity + real(r8),intent(out) :: precl (ncol) ! large-scale precipitation rate (m/s) + real(r8),intent(out) :: evapdt(ncol,pver) ! T tendency due to re-evaporation + real(r8),intent(out) :: evapdq(ncol,pver) ! Q tendency due to re-evaporation + ! + ! Local Values + !------------- + logical,parameter :: do_evap = .true. + + real(r8):: esat (ncol,pver) + real(r8):: qsat (ncol,pver) + real(r8):: dqsat(ncol,pver) + real(r8):: qdel (ncol,pver) + real(r8):: tdel (ncol,pver) + real(r8):: qnew (ncol,pver) + real(r8):: tnew (ncol,pver) + real(r8):: qext (ncol) + real(r8):: qdef (ncol) + + integer :: i, k + + ! Large-Scale Condensation and Precipitation + !-------------------------------------------- + do k = 1,pver + + ! calculate saturation vapor pressure + !------------------------------------- + esat(:,k) = E0*exp(-(latvap_div_rh2o)*((1._r8/T(:,k))-1._r8/T0)) + + ! calculate saturation value for Q + !---------------------------------- + do i = 1,ncol + if(pmid(i,k) > (1._r8-epsilo)*esat(i,k)) then + qsat (i,k) = epsilo*esat(i,k)/pmid(i,k) + dqsat(i,k) = (latvap_div_rh2o)*qsat(i,k)/(T(i,k)**2) + else + qsat (i,k) = 0._r8 + dqsat(i,k) = 0._r8 + endif + end do + + ! if > 100% relative humidity, rain falls out + !--------------------------------------------- + where(((qv(:,k)-qsat(:,k))*qsat(:,k)) > 0._r8) + qdel (:,k) = (qsat(:,k)-qv(:,k))/(1._r8+(latvap_div_cpair)*dqsat(:,k)) + tdel (:,k) = -(latvap_div_cpair)*qdel(:,k) + else where + qdel (:,k) = 0._r8 + tdel (:,k) = 0._r8 + end where + + ! Update temperature and water vapor + !----------------------------------- + qnew(:,k) = qv(:,k) + qdel(:,k) + tnew(:,k) = T(:,k) + tdel(:,k) + end do + + ! optionally allow for re-evaporation of falling precip + !------------------------------------------------------- + if(do_evap) then + ! Initialize work array for excess Q + !-------------------------------------- + qext(:) = 0._r8 + + ! Loop down thru the model levels + !--------------------------------- + do k = 1, pver + + ! Add qdel for the current level to the excess + !---------------------------------------------- + where(qdel(:,k) < 0._r8) qext(:) = qext(:) - qdel(:,k)*pdel(:,k)/gravit + + ! Evaporate excess Q where needed + !---------------------------------- + qdef(:) = 0._r8 + where((qdel(:,k) >= 0._r8).and.(qext(:) > 0._r8)) + qext(:) = qext(:)*gravit/pdel(:,k) + qdef(:) = (qsat(:,k)-qv(:,k))/(1._r8+(latvap_div_cpair)*dqsat(:,k)) + qdef(:) = min(qext(:),max(qdef(:),0._r8)) + qdel(:,k) = qdel(:,k) + qdef(:) + tdel(:,k) = tdel(:,k) -(latvap_div_cpair)*qdef(:) + qext(:) = (qext(:)-qdef(:))*pdel(:,k)/gravit + + ! Update temperature and water vapor + !----------------------------------- + qnew(:,k) = qv(:,k) + qdel(:,k) + tnew(:,k) = T(:,k) + tdel(:,k) + end where + + ! Save T/Q tendencies due to re-evaporation + !-------------------------------------------- + evapdq(:,k) = qdef(:)/dtime + evapdt(:,k) = -qdef(:)*latvap_div_cpair/dtime + end do ! k = 1, pver + else + ! Set T/Q re-evaporation tendencies to 0 + !-------------------------------------------- + evapdt(:,:) = 0._r8 + evapdq(:,:) = 0._r8 + endif + + ! Set large scale precipitation rates to zero + !-------------------------------------------------------------------------- + precl(:) = 0.0_r8 + + ! Calculate resulting precip value and relative humidity + !-------------------------------------------------------- + do k = 1, pver + precl (:) = precl(:) - (qdel(:,k)*pdel(:,k))/(gravit*rhoh2o) + qsat (:,k) = (epsilo/pmid(:,k))*E0*exp(-latvap_div_rh2o*((1._r8/tnew(:,k))-1._r8/T0)) + relhum(:,k) = (qnew(:,k)/qsat (:,k))*100._r8 + end do + precl(:) = max(precl(:),0._r8)/dtime + + ! Update T and qv values due to precipitation + !-------------------------------------------- + qv(:,:) = qnew(:,:) + T (:,:) = tnew(:,:) + + end subroutine frierson_condensate + !======================================================================= + + + !======================================================================= + subroutine frierson_condensate_TJ16(ncol,pver,dtime,pmid,pdel,T,qv,relhum,precl) + ! + ! Precip_process: Implement large-scale condensation/precipitation + ! from TJ16. + ! + !======================================================================= + ! + ! Passed Variables + !--------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: dtime ! time step (s) + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(in) :: pdel (ncol,pver) ! layer thickness (Pa) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(inout):: qv (ncol,pver) ! specific humidity Q (kg/kg) + real(r8),intent(out) :: relhum(ncol,pver) ! relative humidity + real(r8),intent(out) :: precl (ncol) ! large-scale precipitation rate (m/s) + ! + ! Local Values + !------------- + real(r8):: qsat + real(r8):: Crate + integer :: i, k + + ! Set large scale precipitation rates to zero + !-------------------------------------------------------------------------- + precl(:) = 0.0_r8 + + ! Large-Scale Condensation and Precipitation without cloud stage + !--------------------------------------------------------------- + do k = 1, pver + do i = 1, ncol + ! calculate saturation value for Q + !---------------------------------- + qsat = epsilo*E0/pmid(i,k)*exp(-latvap_div_rh2o*((1._r8/T(i,k))-1._r8/T0)) + + ! if > 100% relative humidity rain falls out + !------------------------------------------- + if(qv(i,k) > qsat) then + ! calc the condensation and large-scale precipitation(m/s) rates + !------------------------------------------------------------------- + Crate = ((qv(i,k)-qsat)/dtime) & + /(1._r8+(latvap_div_cpair)*(epsilo*latvap*qsat/(rair*T(i,k)**2))) + precl(i) = precl(i) + (Crate*pdel(i,k))/(gravit*rhoh2o) + + ! Update T and qv values due to precipitation + !-------------------------------------------- + T (i,k) = T (i,k) + Crate*(latvap_div_cpair)*dtime + qv(i,k) = qv(i,k) - Crate*dtime + + ! recompute qsat with updated T + !------------------------------- + qsat = epsilo*E0/pmid(i,k)*exp(-latvap_div_rh2o*((1._r8/T(i,k))-1._r8/T0)) + endif + + ! Set percent relative humidity + !------------------------------- + relhum(i,k) = (qv(i,k)/qsat)*100._r8 + end do + end do + + end subroutine frierson_condensate_TJ16 + !======================================================================= + + + !======================================================================= + subroutine frierson_condensate_USER(ncol,pver,dtime,pmid,pdel,T,qv,relhum,precl) + ! + ! frierson_condensate_USER: This routine is a stub which users can use + ! to develop and test their own large scale + ! condensation scheme + !======================================================================= + ! + ! Passed Variables + !--------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: dtime ! time step (s) + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(in) :: pdel (ncol,pver) ! layer thickness (Pa) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(inout):: qv (ncol,pver) ! specific humidity Q (kg/kg) + real(r8),intent(out) :: relhum(ncol,pver) ! relative humidity + real(r8),intent(out) :: precl (ncol) ! large-scale precipitation rate (m/s) + ! + ! Local Values + !------------- + real(r8):: qsat + integer :: i, k + + ! Set large scale precipitation rates to zero + !-------------------------------------------------------------------------- + precl(:) = 0.0_r8 + + ! Large-Scale Condensation and Precipitation without cloud stage + !--------------------------------------------------------------- + do k = 1, pver + do i = 1, ncol + ! calculate saturation value for Q + !---------------------------------- + qsat = epsilo*E0/pmid(i,k)*exp(-latvap_div_rh2o*((1._r8/T(i,k))-1._r8/T0)) + + ! Set percent relative humidity + !------------------------------- + relhum(i,k) = (qv(i,k)/qsat)*100._r8 + end do + end do + + end subroutine frierson_condensate_USER + !======================================================================= + + + !======================================================================= + subroutine frierson_pbl(ncol, pver, dtime, pmid, pint, Zm, Zi, & + Psfc, Tsfc, Qsfc, T, U, V, Q, Fsw, Fdn, & + Cdrag, Km, Ke, VSE, Z_pbl, Rf, dQa, dTa, dUa, dVa, & + LHflux, SHflux, TUflux, TVflux ) + ! + ! The implicit PBL parameterization based on Frierson, et al. 2006. + ! + ! frierson_pbl: This is an implementation of the implicit computation + ! derived from the code of the Frierson model. The + ! calculations are roughly divided up into sections of + ! the model where they should be carried out. + ! + !========================================================================== + ! + ! Passed Variables + !------------------ + integer ,intent(in) :: ncol ! Number of columns + integer ,intent(in) :: pver ! Number of levels + real(r8),intent(in) :: dtime ! Time Step + real(r8),intent(in) :: pmid (ncol,pver) ! Pressure at model levels + real(r8),intent(in) :: pint (ncol,pver+1) ! Pressure at interface levels. + real(r8),intent(in) :: Zm (ncol,pver) ! Height at model levels. + real(r8),intent(in) :: Zi (ncol,pver) ! Height at interface levels. + real(r8),intent(in) :: Psfc (ncol) ! Surface Pressure. + real(r8),intent(inout):: Tsfc (ncol) ! SST temperature K + real(r8),intent(inout):: Qsfc (ncol) ! sea surface water vapor (kg/kg) + real(r8),intent(inout):: T (ncol,pver) ! ATM Temperature values. + real(r8),intent(inout):: U (ncol,pver) ! ATM Zonal Wind values. + real(r8),intent(inout):: V (ncol,pver) ! ATM Meridional Wind values. + real(r8),intent(inout):: Q (ncol,pver) ! ATM Water vapor values. + real(r8),intent(in) :: Fdn (ncol) ! Downward LW flux at surface + real(r8),intent(in) :: Fsw (ncol) ! Net SW flux at surface from gray radiation + real(r8),intent(out) :: Cdrag (ncol) ! Surface drage coef. + real(r8),intent(out) :: Km (ncol,pver+1) ! Eddy diffusivity for PBL (momentum) + real(r8),intent(out) :: Ke (ncol,pver+1) ! Eddy diffusivity for PBL + real(r8),intent(out) :: VSE (ncol,pver) ! Virtual-Dry Static energy + real(r8),intent(out) :: Z_pbl (ncol) ! Height of PBL layer. + real(r8),intent(out) :: Rf (ncol,pver) + real(r8),intent(out) :: dTa (ncol,pver) + real(r8),intent(out) :: dQa (ncol,pver) + real(r8),intent(out) :: dUa (ncol,pver) + real(r8),intent(out) :: dVa (ncol,pver) + real(r8),intent(out) :: LHflux(ncol) ! Latent Heat Flux + real(r8),intent(out) :: SHflux(ncol) ! Sensible Heat Flux + real(r8),intent(out) :: TUflux(ncol) ! U Surface stress + real(r8),intent(out) :: TVflux(ncol) ! V Surface stress + ! + ! Local Values + !--------------- + real(r8):: Tv (ncol,pver) + real(r8):: Thv (ncol,pver) + real(r8):: Ws (ncol,pver) + real(r8):: rho (ncol) ! Air density near the ground (kg/m3) + real(r8):: Z_sfc (ncol) + real(r8):: Rf_sfc(ncol) + real(r8):: Ri_a (ncol) + real(r8):: Ri (ncol,pver) + integer :: K_sfc (ncol) + integer :: K_pbl (ncol) + real(r8):: Ke_pbl(ncol) + real(r8):: Km_pbl(ncol) + real(r8):: Z_a (ncol) ! Height at midpoint of the lowest model level (m) + real(r8):: Ws_a (ncol) ! wind speed at the lowest model level (m/s) + real(r8):: Thv_a (ncol) + real(r8):: Thv_s (ncol) + real(r8):: Ustar (ncol) + real(r8):: Bstar (ncol) + real(r8):: ZETA,PHI + + real(r8):: MU (ncol,pver) + real(r8):: NUe(ncol,pver) + real(r8):: NUm(ncol,pver) + real(r8):: Am (ncol,pver) + real(r8):: Bm (ncol,pver) + real(r8):: Cm (ncol,pver) + real(r8):: Ae (ncol,pver) + real(r8):: Be (ncol,pver) + real(r8):: Ce (ncol,pver) + real(r8):: FLu(ncol,pver) + real(r8):: FLv(ncol,pver) + real(r8):: FLq(ncol,pver) + real(r8):: FLt(ncol,pver) + real(r8):: Et (ncol,pver) + real(r8):: Eq (ncol,pver) + real(r8):: Eu (ncol,pver) + real(r8):: Ev (ncol,pver) + + real(r8):: Fval_t(ncol,pver) + real(r8):: Fval_q(ncol,pver) + real(r8):: Fval_u(ncol,pver) + real(r8):: Fval_v(ncol,pver) + real(r8):: Eval_m(ncol,pver) + real(r8):: Eval_e(ncol,pver) + integer :: i, k + + real(r8):: Su(ncol,pver) + real(r8):: Sv(ncol,pver) + real(r8):: St(ncol,pver) + real(r8):: Sq(ncol,pver) + + real(r8):: Estar_u(ncol) + real(r8):: Estar_v(ncol) + real(r8):: Estar_q(ncol) + real(r8):: Estar_t(ncol) + real(r8):: dFa_dTa(ncol) + real(r8):: dFa_dQa(ncol) + real(r8):: dFa_dUa(ncol) + real(r8):: dFa_dVa(ncol) + + real(r8):: Th_a (ncol) + real(r8):: Th_s (ncol) + real(r8):: Ft (ncol) + real(r8):: dFt_dTa (ncol) + real(r8):: dFt_dTs (ncol) + real(r8):: Fq (ncol) + real(r8):: dFq_dQa (ncol) + real(r8):: dFq_dTs (ncol) + real(r8):: Fu (ncol) + real(r8):: dFu_dUa (ncol) + real(r8):: Fv (ncol) + real(r8):: dFv_dVa (ncol) + real(r8):: Fup (ncol) + real(r8):: dFup_dTs(ncol) + + real(r8):: FN_u (ncol) + real(r8):: FN_v (ncol) + real(r8):: EN_t (ncol) + real(r8):: FN_t (ncol) + real(r8):: EN_q (ncol) + real(r8):: FN_q (ncol) + real(r8):: Flux (ncol) + real(r8):: dFlux(ncol) + real(r8):: dTs (ncol) + + real(r8):: Tsfc_bc(ncol) + + !============================================================================ + ! tphysbc(): + ! + ! Required Values: + ! T(:,:),Q(:,:),U(:,:),V(:,:) + ! Pmid(:,:),Pint(:,:),Zm(:,:),Zi(:,:) + ! Tsfc(:),Qsfc(:),Psfc(:) + !============================================================================ + + ! Sx() values allow for explicit source tendencies to be passed to + ! implicit PBL calculation. Set all values to 0. for now. + !------------------------------------------------------------------------- + Su(:,:) = 0._r8 + Sv(:,:) = 0._r8 + St(:,:) = 0._r8 + Sq(:,:) = 0._r8 + + ! Calc some values we will need later on + !------------------------------------------ + do k = 1, pver + Ws (:,k) = sqrt(U(:,k)**2 + V(:,k)**2 + Wind_min) + Tv (:,k) = T (:,k)*(1._r8+zvir*Q(:,k)) + Thv(:,k) = Tv(:,k)*((ps0/pmid(:,k))**cappa) + VSE(:,k) = Tv(:,k)+gravit*Zm(:,k)/cpair + end do + + ! Calculate Drag Coef and related values + !----------------------------------------- + do i = 1,ncol + Z_a (i) = Zm (i,pver) + Ws_a (i) = Ws (i,pver) + Thv_a(i) = Thv(i,pver) + Thv_s(i) = Tsfc(i)*(1._r8+zvir*Qsfc(i) )*((ps0/Psfc(i))**cappa) + Ri_a (i) = (gravit*Z_a(i)/(Ws_a(i)**2))*(Thv_a(i)-Thv_s(i))/Thv_s(i) + if(Ri_a(i) <= 0._r8) then + Cdrag(i) = (Karman/log((Z_a(i)/Z0)))**2 + elseif(Ri_a(i) >= Ri_c) then + Cdrag(i) = 0._r8 + else + Cdrag(i) = ((1._r8-(Ri_a(i)/Ri_c))*Karman/log((Z_a(i)/Z0)))**2 + endif + Ustar(i) = sqrt(Cdrag(i))*Ws_a(i) + Bstar(i) = sqrt(Cdrag(i))*(gravit*(Thv_a(i)-Thv_s(i))/Thv_s(i)) + end do + + ! Calculate a bulk Richardson number and determine + ! depths of boundary/surface layers. + !---------------------------------------------------- + do k = 1,pver + Ri(:,k) = (gravit*Zm(:,k)/(Ws(:,k)**2))*(VSE(:,k)-VSE(:,pver))/VSE(:,pver) + Rf(:,k) = Ri(:,k)/Ri_c + end do + + do i =1,ncol + Z_pbl(i) = Zm(i,pver) + K_pbl(i) = pver + do k = (pver-1),1,-1 + if(Rf(i,k) > 1._r8) then + K_pbl(i) = k + 1 + Z_pbl(i) = (Zm(i,k+1)*(Rf(i,k)- 1._r8 ) & + +Zm(i,k )*( 1._r8 - Rf(i,k+1)))/(Rf(i,k)-Rf(i,k+1)) + exit + endif + end do + + ! surface layer height is a fixed fraction of the PBL + ! determine the corresponding level index and Rf value + !----------------------------------------------------- + Z_sfc(i) = Fb*Z_pbl(i) + K_sfc(i) = pver + do k = (pver-1),1,-1 + if(Zm(i,k) > Z_sfc(i)) then + K_sfc (i) = k + 1 + Rf_sfc(i) = (Rf(i,k+1)*(Zm(i,k) - Z_sfc(i) ) & + + Rf(i,k )*(Z_sfc(i) - Zm(i,k+1)))/(Zm(i,k)-Zm(i,k+1)) + exit + endif + end do + end do ! i =1,ncol + + ! Compute diffusion coefs + !------------------------- + Ke(:,:) = 0._r8 + Ke_pbl(:) = 0._r8 + do i = 1,ncol + if (Cdrag(i) /= 0._r8) then + do k = pver,K_pbl(i),-1 + ZETA = Zi(i,k)*Karman*Bstar(i)/(Ustar(i)*Ustar(i)) + if(ZETA < 0._r8) then + if(k >= K_sfc(i)) then + Ke(i,k) = Karman*Ustar(i)*Zi(i,k) + else + Ke(i,k) = Karman*Ustar(i)*Zi(i,k) & + *(((Z_pbl(i)-Zi(i,k))/(Z_pbl(i)-Z_sfc(i)))**2) + endif + elseif (ZETA < Ri_c) then + PHI = 1._r8 + ZETA + if(k >= K_sfc(i)) then + Ke(i,k) = Karman*Ustar(i)*Zi(i,k)/PHI + else + Ke(i,k) = Karman*Ustar(i)*Zi(i,k) & + *(((Z_pbl(i)-Zi(i,k))/(Z_pbl(i)-Z_sfc(i)))**2)/PHI + endif + endif + end do + Ke_pbl(i) = Ke(i,K_sfc(i))*Z_sfc(i)/Zi(i,K_sfc(i)) + endif + end do + + ! The Same coefs used for momentum + !----------------------------------- + Km(:,:) = Ke(:,:) + Km_pbl(:) = Ke_pbl(:) + + ! Compute downward values for the implicit PBL scheme + !----------------------------------------------------- + do k = 1,pver + MU (:,k) = gravit*dtime/(Pint(:,k+1) - Pint(:,k)) + end do + + NUe(:,:) = 0._r8 + NUm(:,:) = 0._r8 + do k = 2,pver + rho(:) = 2._r8*Pint(:,k)/(rair*(Tv(:,k)+Tv(:,k-1))) + NUe(:,k) = rho(:)*Ke(:,k)/(Zm(:,k)-Zm(:,k-1)) + NUm(:,k) = rho(:)*Km(:,k)/(Zm(:,k)-Zm(:,k-1)) + end do + + Am(:,1 ) = MU(:,1)*NUm(:,2) + Cm(:,1 ) = 0._r8 + Am(:,pver) = 0._r8 + Cm(:,pver) = MU(:,pver)*NUm(:,pver) + Ae(:,1 ) = MU(:,1 )*NUe(:,2) + Ce(:,1 ) = 0._r8 + Ae(:,pver) = 0._r8 + Ce(:,pver) = MU(:,pver)*NUe(:,pver) + do k = 2,(pver-1) + Am(:,k) = MU(:,k)*NUm(:,k+1) + Cm(:,k) = MU(:,k)*NUm(:,k ) + Ae(:,k) = MU(:,k)*NUe(:,k+1) + Ce(:,k) = MU(:,k)*NUe(:,k ) + end do + Bm(:,:) = 1._r8 - Am(:,:) - Cm(:,:) + Be(:,:) = 1._r8 - Ae(:,:) - Ce(:,:) + + FLu(:,1) = 0._r8 + FLv(:,1) = 0._r8 + FLq(:,1) = 0._r8 + FLt(:,1) = 0._r8 + do k = 2,pver + FLu(:,k) = NUm(:,k)*(U (:,k)-U (:,k-1)) + FLv(:,k) = NUm(:,k)*(V (:,k)-V (:,k-1)) + FLq(:,k) = NUe(:,k)*(Q (:,k)-Q (:,k-1)) + FLt(:,k) = NUe(:,k)*(VSE(:,k)-VSE(:,k-1)) + end do + do k = 1,(pver-1) + Eu(:,k) = Su(:,k) + MU(:,k)*(FLu(:,k)-FLu(:,k+1)) + Ev(:,k) = Sv(:,k) + MU(:,k)*(FLv(:,k)-FLv(:,k+1)) + Eq(:,k) = Sq(:,k) + MU(:,k)*(FLq(:,k)-FLq(:,k+1)) + Et(:,k) = St(:,k) + MU(:,k)*(FLt(:,k)-FLt(:,k+1)) + end do + Eu(:,pver) = Su(:,pver) + MU(:,pver)*FLu(:,pver) + Ev(:,pver) = Sv(:,pver) + MU(:,pver)*FLv(:,pver) + Eq(:,pver) = Sq(:,pver) + MU(:,pver)*FLq(:,pver) + Et(:,pver) = St(:,pver) + MU(:,pver)*FLt(:,pver) + + Eval_m(:,1) = -Am(:,1)/Bm(:,1) + Eval_e(:,1) = -Ae(:,1)/Be(:,1) + Fval_u(:,1) = Eu(:,1)/Bm(:,1) + Fval_v(:,1) = Ev(:,1)/Bm(:,1) + Fval_q(:,1) = Eq(:,1)/Be(:,1) + Fval_t(:,1) = Et(:,1)/Be(:,1) + do k = 2,(pver-1) + Eval_m(:,k) = -Am(:,k)/(Bm(:,k)+Cm(:,k)*Eval_m(:,k-1)) + Eval_e(:,k) = -Ae(:,k)/(Be(:,k)+Ce(:,k)*Eval_e(:,k-1)) + Fval_u(:,k) = (Eu(:,k)-Cm(:,k)*Fval_u(:,k-1)) & + /(Bm(:,k)+Cm(:,k)*Eval_m(:,k-1)) + Fval_v(:,k) = (Ev(:,k)-Cm(:,k)*Fval_v(:,k-1)) & + /(Bm(:,k)+Cm(:,k)*Eval_m(:,k-1)) + Fval_q(:,k) = (Eq(:,k)-Ce(:,k)*Fval_q(:,k-1)) & + /(Be(:,k)+Ce(:,k)*Eval_e(:,k-1)) + Fval_t(:,k) = (Et(:,k)-Ce(:,k)*Fval_t(:,k-1)) & + /(Be(:,k)+Ce(:,k)*Eval_e(:,k-1)) + end do + Eval_m(:,pver) = 0._r8 + Eval_e(:,pver) = 0._r8 + Fval_u(:,pver) = 0._r8 + Fval_v(:,pver) = 0._r8 + Fval_q(:,pver) = 0._r8 + Fval_t(:,pver) = 0._r8 + + Estar_u(:) = (Eu(:,pver)-Cm(:,pver)*Fval_u(:,pver-1)) + Estar_v(:) = (Ev(:,pver)-Cm(:,pver)*Fval_v(:,pver-1)) + Estar_q(:) = (Eq(:,pver)-Ce(:,pver)*Fval_q(:,pver-1)) + Estar_t(:) = (Et(:,pver)-Ce(:,pver)*Fval_t(:,pver-1)) + + dFa_dTa(:) = NUe(:,pver)*(1._r8-Eval_e(:,pver-1)) + dFa_dQa(:) = NUe(:,pver)*(1._r8-Eval_e(:,pver-1)) + dFa_dUa(:) = NUm(:,pver)*(1._r8-Eval_m(:,pver-1)) + dFa_dVa(:) = NUm(:,pver)*(1._r8-Eval_m(:,pver-1)) + + !============================================================================ + ! flux calculation(): + ! + ! Required Values: + ! Passed from : + ! T(:,pver),Q(:,pver),U(:,pver),V(:,pver),Cdrag(:),Pmid(:,pver) + ! MU(:,pver),dFa_dTa(:),dFa_dQa(:),dFa_dUa(:),dFa_dVa(:) + ! Estar_t(:),Estar_q(:),Estar_u(:),Estar_v(:) + ! Passed from : + ! Tsfc(:),Qsfc(:),Psfc(:) + ! + !============================================================================ + + ! Calculate Surface flux values and their derivatives + !-------------------------------------------------------- + do i = 1, ncol + Th_a(i) = T (i,pver)*((ps0/pmid(i,pver))**cappa) + Th_s(i) = Tsfc(i) *((ps0/Psfc (i) )**cappa) + rho (i) = pmid (i,pver)/(rair*Tv(i,pver)) + + Ft (i) = rho(i)*Cdrag(i)*Ws_a(i)*(Th_s (i) - Th_a(i)) + Fq (i) = rho(i)*Cdrag(i)*Ws_a(i)*(Qsfc(i) - Q(i,pver)) + Fu (i) = -rho(i)*Cdrag(i)*Ws_a(i)*U(i,pver) + Fv (i) = -rho(i)*Cdrag(i)*Ws_a(i)*V(i,pver) + Fup (i) = Boltz*Tsfc(i)**4 + + dFt_dTa(i) = -rho(i)*Cdrag(i)*Ws_a(i)*((ps0/pmid(i,pver))**cappa) + dFq_dQa(i) = -rho(i)*Cdrag(i)*Ws_a(i) + dFu_dUa(i) = -rho(i)*Cdrag(i)*Ws_a(i) + dFv_dVa(i) = -rho(i)*Cdrag(i)*Ws_a(i) + + dFt_dTs(i) = rho(i)*Cdrag(i)*Ws_a(i)*((ps0/Psfc (i) )**cappa) + dFq_dTs(i) = rho(i)*Cdrag(i)*Ws_a(i)*Qsfc(i)*latvap/(rh2o*(Tsfc(i)**2)) + dFup_dTs(i) = 4._r8*Boltz*Tsfc(i)**3 + end do + + ! Incorporate surface fluxes into implicit scheme, then + ! update flux values and derivatives + !------------------------------------------ + FN_u (:) = (Estar_u(:) + MU(:,pver)*Fu(:)) & + /(1._r8-MU(:,pver)*(dFa_dUa(:)+dFu_dUa(:))) + FN_v (:) = (Estar_v(:) + MU(:,pver)*Fv(:)) & + /(1._r8-MU(:,pver)*(dFa_dVa(:)+dFv_dVa(:))) + FN_t (:) = (Estar_t(:) + MU(:,pver)*Ft(:)) & + /(1._r8-MU(:,pver)*(dFa_dTa(:)+dFt_dTa(:))) + FN_q (:) = (Estar_q(:) + MU(:,pver)*Fq(:)) & + /(1._r8-MU(:,pver)*(dFa_dQa(:)+dFq_dQa(:))) + + EN_t (:) = ( MU(:,pver)*dFt_dTs(:) ) & + /(1._r8-MU(:,pver)*(dFa_dTa(:)+dFt_dTa(:))) + EN_q (:) = ( MU(:,pver)*dFq_dTs(:) ) & + /(1._r8-MU(:,pver)*(dFa_dQa(:)+dFq_dQa(:))) + + Ft (:) = Ft(:) + dFt_dTa(:)*FN_t(:) + Fq (:) = Fq(:) + dFq_dQa(:)*FN_q(:) + + dFt_dTs(:) = dFt_dTs(:) + dFt_dTa(:)*EN_t(:) + dFq_dTs(:) = dFq_dTs(:) + dFq_dQa(:)*EN_q(:) + + !============================================================================ + ! surface calculation(): + ! + ! Required Values: + ! Passed from : + ! Fup(:),Ft(:),Fq(:) + ! dFup_dTs(:),dFt_dTs(:),dFq_dTs(:) + ! Fsw(:) + ! Passed from : + ! Fdn(:) + ! + !============================================================================ + + ! Update surface values + !----------------------- + Tsfc_bc(:) = Tsfc(:) + + Flux (:) = (dtime/C_ocn)*(Fdn(:) - Fup(:) + Fsw(:) & + -cpair*Ft(:) -latvap*Fq(:) ) + dFlux(:) = (dtime/C_ocn)*(-dFup_dTs(:) -cpair*dFt_dTs(:) -latvap*dFq_dTs(:)) + Tsfc (:) = Tsfc(:) + (Flux(:)/(1._r8-dFlux(:))) + Qsfc (:) = epsilo*E0/Psfc(:)*exp(-latvap_div_rh2o*((1._r8/Tsfc(:))-1._r8/T0)) + dTs (:) = Tsfc(:) - Tsfc_bc(:) + + LHflux(:) = latvap*Fq(:) + SHflux(:) = cpair *Ft(:) + TUflux(:) = Fu(:) + TVflux(:) = Fv(:) + + !============================================================================ + ! tphysac(): + ! + ! Required Values: + ! Passed from : + ! FN_t(:),FN_q(:),FN_u(:),FN_v(:) + ! EN_t(:),EN_q(:),dTs(:) + ! Passed from : + ! Fval_t(:),Fval_q(:),Fval_u(:),Fval_v(:) + ! Eval_e(:),Eval_m(:) + ! + !============================================================================ + + ! Compute upward values for the implicit PBL scheme + !----------------------------------------------------- + dTa(:,pver) = FN_t(:) + EN_t(:)*dTs(:) + dQa(:,pver) = FN_q(:) + EN_q(:)*dTs(:) + dUa(:,pver) = FN_u(:) + dVa(:,pver) = FN_v(:) + do k=(pver-1),1,-1 + dTa(:,k) = Fval_t(:,k) + Eval_e(:,k)*dTa(:,k+1) + dQa(:,k) = Fval_q(:,k) + Eval_e(:,k)*dQa(:,k+1) + dUa(:,k) = Fval_u(:,k) + Eval_m(:,k)*dUa(:,k+1) + dVa(:,k) = Fval_v(:,k) + Eval_m(:,k)*dVa(:,k+1) + end do + + ! Update atmosphere values + !-------------------------- + U(:,:) = U(:,:) + dUa(:,:) + V(:,:) = V(:,:) + dVa(:,:) + Q(:,:) = Q(:,:) + dQa(:,:) + T(:,:) = T(:,:) + dTa(:,:) + + ! Return resulting Tendency values + !---------------------------------- + dUa(:,:) = dUa(:,:)/dtime + dVa(:,:) = dVa(:,:)/dtime + dQa(:,:) = dQa(:,:)/dtime + dTa(:,:) = dTa(:,:)/dtime + + end subroutine frierson_pbl + !======================================================================= + + + !======================================================================= + subroutine frierson_pbl_USER(ncol, pver, dtime, pmid, pint, Zm, Zi, & + Psfc, Tsfc, Qsfc, T, U, V, Q, Fsw, Fdn, & + Cdrag, Km, Ke, VSE, Z_pbl, Rf, dQa, dTa, dUa, dVa, & + LHflux, SHflux, TUflux, TVflux ) + ! + ! frierson_pbl_USER: This routine is a stub which users can use + ! to develop and test their own PBL scheme + !========================================================================== + ! + ! Passed Variables + !------------------ + integer ,intent(in) :: ncol ! Number of columns + integer ,intent(in) :: pver ! Number of levels + real(r8),intent(in) :: dtime ! Time Step + real(r8),intent(in) :: pmid (ncol,pver) ! Pressure at model levels + real(r8),intent(in) :: pint (ncol,pver+1) ! Pressure at interface levels. + real(r8),intent(in) :: Zm (ncol,pver) ! Height at model levels. + real(r8),intent(in) :: Zi (ncol,pver) ! Height at interface levels. + real(r8),intent(in) :: Psfc (ncol) ! Surface Pressure. + real(r8),intent(inout):: Tsfc (ncol) ! SST temperature K + real(r8),intent(inout):: Qsfc (ncol) ! sea surface water vapor (kg/kg) + real(r8),intent(inout):: T (ncol,pver) ! ATM Temperature values. + real(r8),intent(inout):: U (ncol,pver) ! ATM Zonal Wind values. + real(r8),intent(inout):: V (ncol,pver) ! ATM Meridional Wind values. + real(r8),intent(inout):: Q (ncol,pver) ! ATM Water vapor values. + real(r8),intent(in) :: Fdn (ncol) ! Downward LW flux at surface + real(r8),intent(in) :: Fsw (ncol) ! Net SW flux at surface from gray radiation + real(r8),intent(out) :: Cdrag (ncol) ! Surface drage coef. + real(r8),intent(out) :: Km (ncol,pver+1) ! Eddy diffusivity for PBL + real(r8),intent(out) :: Ke (ncol,pver+1) ! Eddy diffusivity for PBL + real(r8),intent(out) :: VSE (ncol,pver) ! Virtual-Dry Static energy.(huh?) + real(r8),intent(out) :: Z_pbl (ncol) ! Height of PBL layer. + real(r8),intent(out) :: Rf (ncol,pver) + real(r8),intent(out) :: dTa (ncol,pver) + real(r8),intent(out) :: dQa (ncol,pver) + real(r8),intent(out) :: dUa (ncol,pver) + real(r8),intent(out) :: dVa (ncol,pver) + real(r8),intent(out) :: LHflux(ncol) ! Latent Heat Flux + real(r8),intent(out) :: SHflux(ncol) ! Sensible Heat Flux + real(r8),intent(out) :: TUflux(ncol) ! U Surface stress + real(r8),intent(out) :: TVflux(ncol) ! V Surface stress + ! + ! Local Values + !--------------- + + + Cdrag = 0._r8 + Km = 0._r8 + Ke = 0._r8 + VSE = 0._r8 + Z_pbl = 0._r8 + Rf = 0._r8 + dTa = 0._r8 + dQa = 0._r8 + dUa = 0._r8 + dVa = 0._r8 + LHflux = 0._r8 + SHflux = 0._r8 + TUflux = 0._r8 + TVflux = 0._r8 + + end subroutine frierson_pbl_USER + !======================================================================= + + + !======================================================================= + subroutine frierson_radiation(ncol,pver,dtime,clat,pint,pmid, & + Psfc,Tsfc,Qsfc,T,qv,dtdt_rad, & + Fsolar,Fup_s,Fdown_s,Fup_toa,Fdown_toa) + ! + ! The gray radiation parameterization based on Frierson, et al. 2006. + ! + ! frierson_radiation: This is an implementation of the gray radiation + ! scheme used in the Frierson model. + !========================================================================== + ! + ! Passed Variables + !------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: dtime ! time step (s) + real(r8),intent(in) :: clat (ncol) ! latitude + real(r8),intent(in) :: pint (ncol,pver+1) ! mid-point pressure (Pa) + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(in) :: Psfc (ncol) ! surface pressure + real(r8),intent(in) :: Tsfc (ncol) ! surface temperature (K) + real(r8),intent(in) :: Qsfc (ncol) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(in) :: qv (ncol,pver) ! Q (kg/kg) + real(r8),intent(out) :: dtdt_rad(ncol,pver) ! temperature tendency in K/s from relaxation + real(r8),intent(out) :: Fsolar (ncol) ! + real(r8),intent(out) :: Fup_s (ncol) ! + real(r8),intent(out) :: Fdown_toa(ncol) ! + real(r8),intent(out) :: Fup_toa (ncol) ! + real(r8),intent(out) :: Fdown_s (ncol) ! + ! + ! Local Values + !------------- + real(r8):: sinsq (ncol) ! sinlat**2 + real(r8):: Tv_srf (ncol) + real(r8):: Tv (ncol,pver ) + real(r8):: Zm (ncol,pver ) + real(r8):: Zscl (ncol) + real(r8):: Tbar (ncol) + real(r8):: Tdif (ncol) + real(r8):: Zfrac (ncol) + real(r8):: Pfrac (ncol) + real(r8):: Tau_lat(ncol) + real(r8):: Tau (ncol,pver+1) + real(r8):: Zi (ncol,pver+1) + real(r8):: Fup (ncol,pver+1) + real(r8):: Fdown (ncol,pver+1) + real(r8):: Bval (ncol,pver) + real(r8):: Etau (ncol,pver) + integer :: k + + ! Calc current Tv values + !--------------------------------- + Tv_srf(:) = Tsfc(:)*(1._r8+zvir*Qsfc(:)) + do k = 1, pver + Tv(:,k) = T(:,k)*(1._r8+zvir*qv(:,k)) + end do + + ! Calc Geopotential Heights at model interface + ! levels and at model layer levels + !---------------------------------------------- + Zi(:,pver+1) = 0._r8 + Zm(:,pver ) = rair*((Tv(:,pver)+Tv_srf(:))/2._r8)*log(Psfc(:)/pmid(:,pver))/gravit + do k = pver-1,1,-1 + Zscl (:) = rair*log(pmid(:,k+1)/pmid(:,k))/gravit + Tbar (:) = (Tv(:,k)+Tv(:,k+1))/2._r8 + Tdif (:) = (Tv(:,k)-Tv(:,k+1))/2._r8 + Zfrac (:) = log(pint(:,k+1)/pmid(:,k+1))/log(pmid(:,k)/pmid(:,k+1)) + Zm(:,k ) = Zm(:,k+1) + Zscl(:)*Tbar(:) + Zi(:,k+1) = Zm(:,k+1) + Zscl(:)*((Tv(:,k+1)-2._r8*Tdif(:))*Zfrac(:) & + + Tdif(:) *Zfrac(:)**2) + end do + Zfrac(:) = log(pint(:,1)/pmid(:,2))/log(pmid(:,1)/pmid(:,2)) + Zi(:,1) = Zm(:,2) + Zscl(:)*((Tv(:,2)-2._r8*Tdif(:))*Zfrac(:) & + + Tdif(:) *Zfrac(:)**2) + + ! Set Solar flux + !------------------------ + sinsq (:) = sin(clat(:))*sin(clat(:)) + Fsolar(:) = (Rs0/4._r8)*(1._r8 + DeltaS*(1._r8 - 3._r8*sinsq(:))/4._r8) + + ! Calc optical depths + !------------------------ + Tau_lat(:) = Tau_eqtr + (Tau_pole-Tau_eqtr)*sinsq(:) + do k = 1,pver+1 + Pfrac(:) = pint(:,k)/Psfc(:) + Tau(:,k) = Tau_lat(:)*(LinFrac*Pfrac(:) + (1._r8-LinFrac)*Pfrac(:)**4) + end do + + ! Lowest order solution for up/down flux assumes B is constant for the layer + !---------------------------------------------------------------------------- + do k=1,pver + Bval(:,k) = Boltz*T(:,k)**4 + Etau(:,k) = exp(Tau(:,k)-Tau(:,k+1)) + end do + + Fup(:,pver+1) = Boltz*Tsfc(:)**4 + do k=pver,1,-1 + Fup(:,k) = Fup(:,k+1)*Etau(:,k) + Bval(:,k)*(1._r8-Etau(:,k)) + end do + + Fdown(:,1) = 0._r8 + do k=1,pver + Fdown(:,k+1) = Fdown(:,k)*Etau(:,k) + Bval(:,k)*(1._r8-Etau(:,k)) + end do + + ! Calc Radiative heating + !------------------------- + do k=1,pver + dtdt_rad(:,k) = -(cappa*Tv(:,k)/pmid(:,k)) & + *((Fup(:,k+1)-Fdown(:,k+1)) - (Fup(:,k)-Fdown(:,k))) & + /( Zi(:,k+1) - Zi(:,k) ) + end do + + ! Return Upward/Downward long wave radiation at Surface and TOA + !---------------------------------------------------------- + Fup_s (:) = Fup (:,pver+1) + Fdown_s(:) = Fdown(:,pver+1) + Fup_toa (:) = Fup (:,1) + Fdown_toa(:) = Fdown(:,1) + + ! Update T values + !------------------- + do k=1,pver + T(:,k) = T(:,k) + dtdt_rad(:,k)*dtime + end do + + end subroutine frierson_radiation + !======================================================================= + + + !======================================================================= + subroutine frierson_radiation_USER(ncol,pver,dtime,clat,pint,pmid, & + Psfc,Tsfc,Qsfc,T,qv,dtdt_rad, & + Fsolar,Fup_s,Fdown_s,Fup_toa,Fdown_toa) + ! + ! frierson_radiation_USER: This routine is a stub which users can use + ! to develop and test their own radiation scheme + !========================================================================== + ! + ! Passed Variables + !------------------- + integer ,intent(in) :: ncol ! number of columns + integer ,intent(in) :: pver ! number of vertical levels + real(r8),intent(in) :: dtime ! time step (s) + real(r8),intent(in) :: clat (ncol) ! latitude + real(r8),intent(in) :: pint (ncol,pver+1) ! mid-point pressure (Pa) + real(r8),intent(in) :: pmid (ncol,pver) ! mid-point pressure (Pa) + real(r8),intent(in) :: Psfc (ncol) ! surface pressure + real(r8),intent(in) :: Tsfc (ncol) ! surface temperature (K) + real(r8),intent(in) :: Qsfc (ncol) + real(r8),intent(inout):: T (ncol,pver) ! temperature (K) + real(r8),intent(in) :: qv (ncol,pver) ! Q (kg/kg) + real(r8),intent(out) :: dtdt_rad(ncol,pver) ! temperature tendency in K/s from relaxation + real(r8),intent(out) :: Fsolar (ncol) ! + real(r8),intent(out) :: Fup_s (ncol) ! + real(r8),intent(out) :: Fdown_s (ncol) ! + real(r8),intent(out) :: Fup_toa (ncol) ! + real(r8),intent(out) :: Fdown_toa(ncol) ! + ! + ! Local Values + !------------- + + dtdt_rad = 0._r8 + Fsolar = 0._r8 + Fup_s = 0._r8 + Fdown_s = 0._r8 + Fup_toa = 0._r8 + Fdown_toa = 0._r8 + + + end subroutine frierson_radiation_USER + !======================================================================= + +end module frierson diff --git a/src/physics/simple/frierson_cam.F90 b/src/physics/simple/frierson_cam.F90 new file mode 100644 index 0000000000..0fe9d824bb --- /dev/null +++ b/src/physics/simple/frierson_cam.F90 @@ -0,0 +1,1046 @@ +module frierson_cam +!----------------------------------------------------------------------- +! +! Purpose: Implement idealized forcings described in +! Frierson, et al. (2006), " A Gray-Radiation Aquaplanet +! Moist GCM, Part I. Static Stability and Eddy Scale" +! J. Atmos. Sci, Vol 63, 2548-2566. +! doi: 10.1175/JAS3753.1 +! +!============================================================================ + ! Useful modules + !------------------- + use shr_kind_mod, only: r8 => shr_kind_r8 + use shr_const_mod, only: pi => shr_const_pi + use physconst, only: gravit, cappa, rair, cpair, latvap, rh2o, epsilo, rhoh2o, zvir + use ppgrid, only: pcols, pver, pverp, begchunk, endchunk + use constituents, only: pcnst + use physics_buffer, only: dtype_r8, pbuf_add_field, physics_buffer_desc, & + pbuf_set_field, pbuf_get_field + use camsrfexch, only: cam_in_t,cam_out_t + use cam_history, only: outfld + use time_manager, only: is_first_step + use cam_abortutils, only: endrun + use spmd_utils, only: masterproc + use cam_logfile, only: iulog + use hycoef, only: ps0, etamid + use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_real8 + + use pio ,only: file_desc_t, var_desc_t, io_desc_t, pio_double, pio_def_var + use pio ,only: pio_write_darray, pio_read_darray, pio_inq_varid + use cam_grid_support,only: cam_grid_id, cam_grid_dimensions, cam_grid_get_decomp + use shr_const_mod, only: SHR_CONST_STEBOL, SHR_CONST_REARTH, SHR_CONST_KARMAN, SHR_CONST_TKTRIP + + ! Set all Global values and routines to private by default + ! and then explicitly set their exposure. + !--------------------------------------------------------- + implicit none + private + save + + public :: frierson_register + public :: frierson_readnl + public :: frierson_init + public :: frierson_condensate_tend + public :: frierson_pbl_tend + public :: frierson_radiative_tend + public :: frierson_restart_init + public :: frierson_restart_write + public :: frierson_restart_read + + private :: frierson_surface_init + + ! PBL Configuatons + !------------------ + integer,parameter :: PBL_FRIERSON = 0 ! Implementation of Frierson PBL + integer,parameter :: PBL_USER = 1 ! Optional call for user defined PBL + + ! Tags to identify optional model formulations + !------------------------------------------------ + integer,parameter :: CONDENSATE_NONE = 0 ! No Condensation, PRECL=0 + integer,parameter :: CONDENSATE_FRIERSON = 1 ! Frierson condensation w/ re-evaporation + integer,parameter :: CONDENSATE_TJ16 = 2 ! Condensation from TJ2016 model. + integer,parameter :: CONDENSATE_USER = 3 ! Optional user defined Condensation scheme + + integer,parameter :: RADIATION_FRIERSON = 0 ! Frierson Gray radiation. + integer,parameter :: RADIATION_USER = 1 ! Optional user defined Radiation scheme + + ! Options selecting which PRECIP, PBL, RADIATION, etc.. formulations to use. + !--------------------------------------------------------------------------------- + integer,parameter :: PBL_OPT = PBL_FRIERSON + integer,parameter :: CONDENSATE_OPT = CONDENSATE_FRIERSON + integer,parameter :: RADIATION_OPT = RADIATION_FRIERSON + + ! Global Constants + !--------------------- + real(r8),parameter :: frierson_T0 = SHR_CONST_TKTRIP ! Reference Temperature for E0 + real(r8),parameter :: frierson_E0 = 610.78_r8 ! Saturation Vapor pressure @ T0 + real(r8),parameter :: frierson_Rs0 = 1360.0_r8 ! Solar Constant + real(r8),parameter :: frierson_Erad = SHR_CONST_REARTH ! Earth Radius + real(r8),parameter :: frierson_Karman = SHR_CONST_KARMAN ! Von Karman constant + real(r8),parameter :: frierson_Boltz = SHR_CONST_STEBOL ! Stefan-Boltzmann constant + + ! Some Physics buffer indices + !------------------------------- + integer :: prec_pcw_idx = 0 + integer :: prec_dp_idx = 0 + integer :: relhum_idx = 0 + + ! Global values for Surface Temp, surface fluxes, and radiative heating + !---------------------------------------------------------------------- + type(var_desc_t) :: Tsurf_desc ! Vardesc for restarts + type(var_desc_t) :: Qsurf_desc ! Vardesc for restarts + real(r8),allocatable :: Tsurf (:,:) ! Surface Temp + real(r8),allocatable :: Qsurf (:,:) ! Surface Q + real(r8),allocatable :: Fsolar(:,:) ! Net Solar Heating + real(r8),allocatable :: Fup (:,:) ! Upward Longwave flux + real(r8),allocatable :: Fdown (:,:) ! Downward Longwave flux + real(r8),allocatable :: Fup_toa (:,:) ! Upward Longwave flux at TOA + real(r8),allocatable :: Fdown_toa(:,:) ! Downward Longwave flux at TOA + real(r8),allocatable :: SHflux(:,:) ! Sensible Heat flux + real(r8),allocatable :: LHflux(:,:) ! Latent Heat Flux + real(r8),allocatable :: TUflux(:,:) ! U stress momentum flux + real(r8),allocatable :: TVflux(:,:) ! V stress momentum flux + real(r8),allocatable :: Cd (:,:) ! Surface Drag + real(r8),allocatable :: clat (:,:) ! latitudes(radians) for columns + real(r8),allocatable :: Fnet (:,:) ! Net Radiative Surface Heating + real(r8),allocatable :: Fnet_toa(:,:) ! Net Radiative Surface Heating at TOA + + real(r8), parameter :: unset_r8 = huge(1.0_r8) + + ! Global Tuning values + !------------------------ + real(r8) :: frierson_Wind_min = unset_r8 ! Minimum wind threshold + real(r8) :: frierson_Z0 = unset_r8 ! Roughness Length + real(r8) :: frierson_Ri_c = unset_r8 ! Crit. Richardson # for stable mixing + real(r8) :: frierson_Fb = unset_r8 ! Surface layer Fraction + real(r8) :: frierson_Albedo = unset_r8 ! Frierson Albedo + real(r8) :: frierson_DeltaS = unset_r8 ! Lat variation of shortwave radiation + real(r8) :: frierson_Tau_eqtr = unset_r8 ! Longwave optical depth at Equator + real(r8) :: frierson_Tau_pole = unset_r8 ! Longwave optical depth at poles. + real(r8) :: frierson_LinFrac = unset_r8 ! Stratosphere Linear optical depth param + real(r8) :: frierson_C0 = unset_r8 ! Ocean mixed layer heat capacity + real(r8) :: frierson_WetDryCoef = unset_r8 ! E0 Scale factor to control moisture + real(r8) :: frierson_Tmin = unset_r8 ! IC: Minimum sst (K) + real(r8) :: frierson_Tdlt = unset_r8 ! IC: eq-polar difference sst (K) + real(r8) :: frierson_Twidth = unset_r8 ! IC: Latitudinal width parameter for sst (degrees latitude) + +contains + !============================================================================== + subroutine frierson_register() + ! + ! frierson_register: Register physics buffer values + !===================================================================== + + call pbuf_add_field('PREC_PCW','physpkg',dtype_r8, (/pcols/), prec_pcw_idx) + call pbuf_add_field('PREC_DP' ,'physpkg',dtype_r8, (/pcols/), prec_dp_idx ) + call pbuf_add_field('RELHUM' ,'physpkg',dtype_r8, (/pcols,pver/),relhum_idx ) + + end subroutine frierson_register + !============================================================================== + + + !============================================================================== + subroutine frierson_readnl(nlfile) + ! + ! frierson_readnl: Read in parameters controlling Frierson parameterizations. + !===================================================================== + use namelist_utils,only: find_group_name + use units ,only: getunit, freeunit + ! + ! Passed Variables + !------------------ + character(len=*),intent(in):: nlfile + ! + ! Local Values + !-------------- + integer:: ierr,unitn + + character(len=*), parameter :: sub = 'frierson_readnl' + + namelist /frierson_nl/ frierson_Wind_min, frierson_Z0 , frierson_Ri_c , & + frierson_Fb , frierson_Albedo , frierson_DeltaS , & + frierson_Tau_eqtr, frierson_Tau_pole , frierson_LinFrac, & + frierson_C0 , frierson_WetDryCoef, frierson_Tmin , & + frierson_Tdlt , frierson_Twidth + + ! Read in namelist values + !------------------------- + if(masterproc) then + unitn = getunit() + open(unitn,file=trim(nlfile),status='old') + call find_group_name(unitn,'frierson_nl',status=ierr) + if(ierr == 0) then + read(unitn,frierson_nl,iostat=ierr) + if(ierr /= 0) then + call endrun(sub//': ERROR reading namelist') + endif + endif + close(unitn) + call freeunit(unitn) + endif + + ! Broadcast namelist values + !--------------------------- + call mpi_bcast(frierson_Wind_min , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Wind_min") + call mpi_bcast(frierson_Z0 , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Z0") + call mpi_bcast(frierson_Ri_c , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Ri_c") + call mpi_bcast(frierson_Fb , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Fb") + call mpi_bcast(frierson_Albedo , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Albedo") + call mpi_bcast(frierson_DeltaS , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_DeltaS") + call mpi_bcast(frierson_Tau_eqtr , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Tau_eqtr") + call mpi_bcast(frierson_Tau_pole , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Tau_pole") + call mpi_bcast(frierson_LinFrac , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_LinFrac") + call mpi_bcast(frierson_C0 , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_C0") + call mpi_bcast(frierson_Tmin , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Tmin") + call mpi_bcast(frierson_Tdlt , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Tdlt") + call mpi_bcast(frierson_Twidth , 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_Twidth") + call mpi_bcast(frierson_WetDryCoef, 1, mpi_real8 , mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: frierson_WetDryCoef") + + end subroutine frierson_readnl + !============================================================================== + + + !============================================================================== + subroutine frierson_init(phys_state,pbuf2d) + ! + ! frierson_init: allocate space for global arrays and initialize values. + ! Add variables to history outputs + !===================================================================== + use physics_types, only: physics_state + use error_messages,only: alloc_err + use cam_history, only: addfld, add_default,horiz_only + use phys_grid, only: get_ncols_p, get_rlat_p + use frierson, only: frierson_set_const + ! + ! Passed Variables + !------------------ + type(physics_state) ,pointer:: phys_state(:) + type(physics_buffer_desc),pointer:: pbuf2d (:,:) + ! + ! Local Values + !--------------- + integer :: istat,lchnk,icol,ncol + real(r8):: adjusted_E0 + real(r8):: frierson_Rs + + ! Initialize constants in Frierson module + !------------------------------------------ + adjusted_E0 = frierson_WetDryCoef*frierson_E0 + frierson_Rs = frierson_Rs0*(1._r8-frierson_Albedo) + call frierson_set_const(gravit,cappa,rair,cpair,latvap,rh2o,epsilo,rhoh2o,zvir,ps0,etamid, & + frierson_T0 ,adjusted_E0 ,frierson_Erad ,frierson_Wind_min, & + frierson_Z0 ,frierson_Ri_c ,frierson_Karman ,frierson_Fb , & + frierson_Rs ,frierson_DeltaS,frierson_Tau_eqtr,frierson_Tau_pole, & + frierson_LinFrac,frierson_Boltz ,frierson_C0 ) + + ! Add values for history output + !--------------------------------- + call addfld('gray_QRL' ,(/'lev' /),'A','K/s' ,'Longwave heating rate for gray atmosphere' ) + call addfld('gray_QRS' ,(/'lev' /),'A','K/s' ,'Solar heating rate for gray atmosphere' ) + call addfld('gray_DTCOND',(/'lev' /),'A','K/s' ,'T tendency - gray atmosphere moist process' ) + call addfld('gray_DQCOND',(/'lev' /),'A','kg/kg/s','Q tendency - gray atmosphere moist process' ) + call addfld('gray_EVAPDT',(/'lev' /),'A','K/s' ,'T tendency due to re-evaporation' ) + call addfld('gray_EVAPDQ',(/'lev' /),'A','kg/kg/s','Q tendency due to re-evaporation' ) + call addfld('gray_KVH' ,(/'ilev'/),'A','m2/s' ,'Vertical diffusion diffusivities (heat/moisture)') + call addfld('gray_KVM' ,(/'ilev'/),'A','m2/s' ,'Vertical diffusion diffusivities (momentum)' ) + call addfld('gray_VSE' ,(/'lev' /),'A','K' ,'VSE: (Tv + gZ/Cp)' ) + call addfld('gray_Zm' ,(/'lev' /),'A','m' ,'Geopotential height' ) + call addfld('gray_Rf' ,(/'lev' /),'A','1' ,'Bulk Richardson number (Frierson et al 2006, eq 16) / Ri_c' ) + call addfld('gray_DTV' ,(/'lev' /),'A','K/s' ,'T tendency due to vertical diffusion' ) + call addfld('gray_DUV' ,(/'lev' /),'A','m/s2' ,'U tendency due to vertical diffusion' ) + call addfld('gray_DVV' ,(/'lev' /),'A','m/s2' ,'V tendency due to vertical diffusion' ) + call addfld('gray_VD01' ,(/'lev' /),'A','kg/kg/s','Q tendency (vertical diffusion)' ) + call addfld('gray_PRECL' ,horiz_only,'A','m/s' ,'Large-scale precipitation rate' ) + call addfld('gray_PRECC' ,horiz_only,'A','m/s' ,'Convective precipitation rate' ) + call addfld('gray_Tsurf ',horiz_only,'I','K' ,'Surface Temperature' ) + call addfld('gray_Qsurf ',horiz_only,'I','kg/kg' ,'Surface Water Vapor' ) + call addfld('gray_Cdrag' ,horiz_only,'A','1' ,'Surface Drag Coefficient' ) + call addfld('gray_Zpbl' ,horiz_only,'I','m' ,'PBL Height' ) + call addfld('gray_SWflux',horiz_only,'I','W/m2' ,'SW Solar Flux' ) + call addfld('gray_LUflux',horiz_only,'I','W/m2' ,'LW Upward Radiative Flux at Surface' ) + call addfld('gray_LDflux',horiz_only,'I','W/m2' ,'LW Downward Radiative Flux at Surface' ) + call addfld('gray_LWflux',horiz_only,'I','W/m2' ,'LW Net Radiative Flux at Surface' ) + call addfld('gray_LUflux_TOA',horiz_only,'I','W/m2' ,'LW Upward Radiative Flux at TOA' ) + call addfld('gray_LDflux_TOA',horiz_only,'I','W/m2' ,'LW Downward Radiative Flux at TOA' ) + call addfld('gray_LWflux_TOA',horiz_only,'I','W/m2' ,'LW Net Radiative Flux at TOA' ) + call addfld('gray_SHflux',horiz_only,'I','W/m2' , 'Sensible Heat Flux' ) + call addfld('gray_LHflux',horiz_only,'I','W/m2' , 'Latent Heat Flux' ) + call addfld('gray_TauU' ,horiz_only,'I','N/m2' , 'U Surface Stress' ) + call addfld('gray_TauV' ,horiz_only,'I','N/m2' , 'V Surface Stress' ) + + call add_default('gray_QRL' ,1,' ') + call add_default('gray_QRS' ,1,' ') + call add_default('gray_DTCOND',1,' ') + call add_default('gray_DQCOND',1,' ') + call add_default('gray_EVAPDT',1,' ') + call add_default('gray_EVAPDQ',1,' ') + call add_default('gray_KVH' ,1,' ') + call add_default('gray_KVM' ,1,' ') + call add_default('gray_VSE' ,1,' ') + call add_default('gray_Zm' ,1,' ') + call add_default('gray_Rf' ,1,' ') + call add_default('gray_DTV' ,1,' ') + call add_default('gray_DUV' ,1,' ') + call add_default('gray_DVV' ,1,' ') + call add_default('gray_VD01' ,1,' ') + call add_default('gray_PRECC' ,1,' ') + call add_default('gray_PRECL' ,1,' ') + call add_default('gray_Tsurf' ,1,' ') + call add_default('gray_Qsurf' ,1,' ') + call add_default('gray_Cdrag' ,1,' ') + call add_default('gray_Zpbl' ,1,' ') + call add_default('gray_SWflux',1,' ') + call add_default('gray_LUflux',1,' ') + call add_default('gray_LDflux',1,' ') + call add_default('gray_LWflux',1,' ') + call add_default('gray_LUflux_TOA',1,' ') + call add_default('gray_LDflux_TOA',1,' ') + call add_default('gray_LWflux_TOA',1,' ') + call add_default('gray_SHflux',1,' ') + call add_default('gray_LHflux',1,' ') + call add_default('gray_TauU' ,1,' ') + call add_default('gray_TauV' ,1,' ') + + ! Allocate Global arrays + !------------------------- + allocate(Fsolar(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fsolar',pcols*(endchunk-begchunk+1)) + allocate(Fup (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fup' ,pcols*(endchunk-begchunk+1)) + allocate(Fdown (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fdown' ,pcols*(endchunk-begchunk+1)) + allocate(Fup_toa (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fup_toa' ,pcols*(endchunk-begchunk+1)) + allocate(Fdown_toa (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fdown_toa' ,pcols*(endchunk-begchunk+1)) + allocate(Fnet(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fnet',pcols*(endchunk-begchunk+1)) + allocate(Fnet_toa(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Fnet_toa' ,pcols*(endchunk-begchunk+1)) + + allocate(SHflux(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','SHflux',pcols*(endchunk-begchunk+1)) + allocate(LHflux(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','LHflux',pcols*(endchunk-begchunk+1)) + allocate(TUflux(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','TUflux',pcols*(endchunk-begchunk+1)) + allocate(TVflux(pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','TVflux',pcols*(endchunk-begchunk+1)) + allocate(Cd (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Cd' ,pcols*(endchunk-begchunk+1)) + allocate(clat (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','clat' ,pcols*(endchunk-begchunk+1)) + + ! Initialize time indices and latitudes + !---------------------------------------- + do lchnk = begchunk,endchunk + ncol = get_ncols_p(lchnk) + do icol = 1,ncol + clat(icol,lchnk) = get_rlat_p(lchnk,icol) + end do + end do + + ! At first model step, initialize some values + !----------------------------------------------- + if(is_first_step()) then + ! Initialize physics buffer values + !---------------------------------- + call pbuf_set_field(pbuf2d, prec_pcw_idx, 0._r8) + call pbuf_set_field(pbuf2d, prec_dp_idx , 0._r8) + + ! Allocate Surface fields + !------------------------- + allocate(Tsurf (pcols,begchunk:endchunk),stat=istat) + call alloc_err(istat,'Frierson INIT','Tsurf' ,pcols*(endchunk-begchunk+1)) + allocate(Qsurf (pcols,begchunk:endchunk) ,stat=istat) + call alloc_err(istat,'Frierson INIT','Qsurf' ,pcols*(endchunk-begchunk+1)) + + ! Initialize Surface temperatures and Q + !----------------------------------------------------------------------- + do lchnk = begchunk,endchunk + ncol = get_ncols_p(lchnk) + + ! Set to reference values for initialization + !------------------------------------------------------------ + phys_state(lchnk)%ps(:ncol) = ps0 + + call frierson_surface_init(ncol, clat(:ncol,lchnk), & + phys_state(lchnk)%ps(:ncol), & + Tsurf(:ncol,lchnk), & + Qsurf(:ncol,lchnk) ) + end do + endif + + ! Initialize radiation and flux values to 0.0 + !--------------------------------------------------------------------------- + do lchnk = begchunk,endchunk + Fsolar(:,lchnk) = 0._r8 + Fup (:,lchnk) = 0._r8 + Fdown (:,lchnk) = 0._r8 + Fup_toa (:,lchnk) = 0._r8 + Fdown_toa (:,lchnk) = 0._r8 + SHflux(:,lchnk) = 0._r8 + LHflux(:,lchnk) = 0._r8 + TUflux(:,lchnk) = 0._r8 + TVflux(:,lchnk) = 0._r8 + Cd (:,lchnk) = 0._r8 + Fnet (:,lchnk) = 0._r8 + end do + + ! Informational Output + !---------------------- + if(masterproc) then + write(iulog,*) ' ' + write(iulog,*) '-----------------------------------------------------------' + write(iulog,*) ' FRIERSON MODULE INITIALIZED WITH THE FOLLOWING SETTINGS: ' + write(iulog,*) '-----------------------------------------------------------' + write(iulog,*) 'FRIERSON: gravit=' , gravit + write(iulog,*) 'FRIERSON: cappa=' , cappa + write(iulog,*) 'FRIERSON: rair =' , rair + write(iulog,*) 'FRIERSON: cpair=' , cpair + write(iulog,*) 'FRIERSON: latvap=' , latvap + write(iulog,*) 'FRIERSON: rh2o=' , rh2o + write(iulog,*) 'FRIERSON: epsilo=' , epsilo + write(iulog,*) 'FRIERSON: rhoh2o=' , rhoh2o + write(iulog,*) 'FRIERSON: zvir=' , zvir + write(iulog,*) 'FRIERSON: ps0=' , ps0 + write(iulog,*) 'FRIERSON: etamid=' , etamid + write(iulog,*) 'FRIERSON: T0=' , frierson_T0 + write(iulog,*) 'FRIERSON: E0=' , frierson_E0 + write(iulog,*) 'FRIERSON: Erad=' , frierson_Erad + write(iulog,*) 'FRIERSON: Wind_min=' , frierson_Wind_min + write(iulog,*) 'FRIERSON: Z0=' , frierson_Z0 + write(iulog,*) 'FRIERSON: Ri_c=' , frierson_Ri_c + write(iulog,*) 'FRIERSON: Karman=' , frierson_Karman + write(iulog,*) 'FRIERSON: Fb=' , frierson_Fb + write(iulog,*) 'FRIERSON: Rs0=' , frierson_Rs0 + write(iulog,*) 'FRIERSON: Albedo=' , frierson_Albedo + write(iulog,*) 'FRIERSON: Rs=' , frierson_Rs + write(iulog,*) 'FRIERSON: DeltaS=' , frierson_DeltaS + write(iulog,*) 'FRIERSON: Tau_eqtr=' , frierson_Tau_eqtr + write(iulog,*) 'FRIERSON: Tau_pole=' , frierson_Tau_pole + write(iulog,*) 'FRIERSON: LinFrac=' , frierson_LinFrac + write(iulog,*) 'FRIERSON: Boltz=' , frierson_Boltz + write(iulog,*) 'FRIERSON: C0=' , frierson_C0 + write(iulog,*) 'FRIERSON: Tmin=' , frierson_Tmin + write(iulog,*) 'FRIERSON: Tdlt=' , frierson_Tdlt + write(iulog,*) 'FRIERSON: Twidth=' , frierson_Twidth + write(iulog,*) 'FRIERSON: WetDryCoef=', frierson_WetDryCoef + write(iulog,*) ' ' + endif + + end subroutine frierson_init + !============================================================================== + + + !============================================================================== + + + !============================================================================== + subroutine frierson_condensate_tend(state, ptend, ztodt, pbuf) + ! + ! frierson_condensate_tend: Run the selected process to compute precipitation + ! due to large scale condensation. + !===================================================================== + use physics_types,only: physics_state, physics_ptend + use physics_types,only: physics_ptend_init + use frierson, only: frierson_condensate_NONE,frierson_condensate + use frierson, only: frierson_condensate_USER,frierson_condensate_TJ16 + ! + ! Passed Variables + !------------------ + type(physics_state) ,intent(inout):: state + real(r8) ,intent(in) :: ztodt + type(physics_ptend) ,intent(out) :: ptend + type(physics_buffer_desc),pointer :: pbuf(:) + ! + ! Local Values + !----------------- + real(r8),pointer:: relhum (:,:) + real(r8),pointer:: prec_pcw(:) ! large scale precip + real(r8) :: prec_cnv(state%ncol) ! Convective Precip + real(r8) :: evapdt(state%ncol, pver) ! T tendency due to re-evaporation of condensation + real(r8) :: evapdq(state%ncol, pver) ! Q tendency due to re-evaporation of condensation + real(r8) :: dtcond(state%ncol, pver) ! Temperature tendency due to condensation + real(r8) :: dqcond(state%ncol, pver) ! Q tendency due to condensation + real(r8) :: T (state%ncol, pver) ! T temporary + real(r8) :: qv (state%ncol, pver) ! Q temporary + logical :: lq(pcnst) ! Calc tendencies? + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: k + + ! Set local copies of values + !--------------------------------- + lchnk = state%lchnk + ncol = state%ncol + T (:ncol,:) = state%T(:ncol,:) + qv(:ncol,:) = state%Q(:ncol,:,1) + + ! initialize individual parameterization tendencies + !--------------------------------------------------- + lq = .false. + lq(1) = .true. + call physics_ptend_init(ptend, state%psetcols, 'Frierson condensate', & + ls=.true., lu=.true., lv=.true., lq=lq) + + ! Get values from the physics buffer + !------------------------------------ + call pbuf_get_field(pbuf,prec_pcw_idx,prec_pcw) + call pbuf_get_field(pbuf, relhum_idx,relhum ) + + ! Initialize values for condensate tendencies + !--------------------------------------------- + do k = 1, pver + dtcond(:ncol,k) = state%T(:ncol,k) + dqcond(:ncol,k) = state%q(:ncol,k,1) + end do + + ! Call the Selected condensation routine ~~DEVO style~~ + !-------------------------------------------------------- + if(CONDENSATE_OPT == CONDENSATE_NONE) then + prec_cnv(:ncol) = 0._r8 + evapdt (:ncol,:) = 0._r8 + evapdq (:ncol,:) = 0._r8 + call frierson_condensate_NONE(ncol,pver,state%pmid(:ncol,:), & + T(:ncol,:), & + qv(:ncol,:), & + relhum(:ncol,:), & + prec_pcw(:ncol) ) + elseif(CONDENSATE_OPT == CONDENSATE_FRIERSON) then + prec_cnv(:ncol) = 0._r8 + call frierson_condensate(ncol,pver,ztodt,state%pmid(:ncol,:), & + state%pdel(:ncol,:), & + T(:ncol,:), & + qv(:ncol,:), & + relhum(:ncol,:), & + prec_pcw(:ncol) , & + evapdt(:ncol,:), & + evapdq(:ncol,:) ) + elseif(CONDENSATE_OPT == CONDENSATE_TJ16) then + prec_cnv(:ncol) = 0._r8 + evapdt (:ncol,:) = 0._r8 + evapdq (:ncol,:) = 0._r8 + call frierson_condensate_TJ16(ncol,pver,ztodt,state%pmid(:ncol,:), & + state%pdel(:ncol,:), & + T(:ncol,:), & + qv(:ncol,:), & + relhum(:ncol,:), & + prec_pcw(:ncol) ) + elseif(CONDENSATE_OPT == CONDENSATE_USER) then + prec_cnv(:ncol) = 0._r8 + evapdt (:ncol,:) = 0._r8 + evapdq (:ncol,:) = 0._r8 + call frierson_condensate_USER(ncol,pver,ztodt,state%pmid(:ncol,:), & + state%pdel(:ncol,:), & + T(:ncol,:), & + qv(:ncol,:), & + relhum(:ncol,:), & + prec_pcw(:ncol) ) + else + ! ERROR: Unknown CONDENSATE_OPT value + !------------------------------------- + write(iulog,*) 'ERROR: unknown CONDENSATE_OPT=',CONDENSATE_OPT + call endrun('frierson_condensate_tend() CONDENSATE_OPT ERROR') + endif + + ! Back out temperature and specific humidity + ! tendencies from updated fields + !-------------------------------------------- + do k = 1, pver + ptend%s(:ncol,k) = (T (:,k)-state%T(:ncol,k) )/ztodt*cpair + ptend%q(:ncol,k,1) = (qv(:,k)-state%q(:ncol,k,1))/ztodt + end do + + ! Output condensate tendencies + !------------------------------ + do k = 1, pver + dtcond(:ncol,k) = (T (:ncol,k) - dtcond(:ncol,k))/ztodt + dqcond(:ncol,k) = (qv(:ncol,k) - dqcond(:ncol,k))/ztodt + end do + call outfld('gray_EVAPDT',evapdt ,ncol,lchnk) + call outfld('gray_EVAPDQ',evapdq ,ncol,lchnk) + call outfld('gray_DTCOND',dtcond ,ncol,lchnk) + call outfld('gray_DQCOND',dqcond ,ncol,lchnk) + call outfld('gray_PRECL' ,prec_pcw,ncol,lchnk) + call outfld('gray_PRECC' ,prec_cnv,ncol,lchnk) + + end subroutine frierson_condensate_tend + !============================================================================== + + + !============================================================================ + subroutine frierson_pbl_tend(state, ptend, ztodt, cam_in) + ! + ! frierson_pbl_tend: Run the selected PBL process. + !========================================================================= + use physics_types,only: physics_state, physics_ptend + use physics_types,only: physics_ptend_init + use phys_grid, only: get_rlat_all_p + use frierson, only: frierson_pbl,frierson_pbl_USER + ! + ! Passed Variables + !------------------- + type(physics_state),intent(in) :: state + real(r8), intent(in) :: ztodt + type(physics_ptend),intent(out) :: ptend + type(cam_in_t), intent(inout):: cam_in + ! + ! Local Values + !---------------- + real(r8) :: T (state%ncol,pver) ! T temporary + real(r8) :: qv (state%ncol,pver) ! Q temporary (specific humidity) + real(r8) :: U (state%ncol,pver) ! U temporary + real(r8) :: V (state%ncol,pver) ! V temporary + real(r8) :: dqdt_vdiff(state%ncol,pver) ! PBL Q vertical diffusion tend kg/kg/s + real(r8) :: dtdt_vdiff(state%ncol,pver) ! PBL T vertical diffusion tend K/s + real(r8) :: dudt_vdiff(state%ncol,pver) ! PBL U vertical diffusion tend m/s/s + real(r8) :: dvdt_vdiff(state%ncol,pver) ! PBL V vertical diffusion tend m/s/s + real(r8) :: Km (state%ncol,pverp) ! Eddy diffusivity at layer interfaces (m2/s) + real(r8) :: Ke (state%ncol,pverp) ! Eddy diffusivity at layer interfaces (m2/s) + real(r8) :: VSE (state%ncol,pver) ! Dry Static Energy divided by Cp (K) + real(r8) :: Zm (state%ncol,pver) ! + real(r8) :: Zi (state%ncol,pver) ! + real(r8) :: Z_pbl (state%ncol) ! + real(r8) :: Rf (state%ncol,pver) ! + real(r8) :: Tsfc (state%ncol) ! Surface T + real(r8) :: Qsfc (state%ncol) ! Surface Q (saturated) + real(r8) :: Cdrag (state%ncol) ! Cdrag coef from surface calculation + + logical :: lq (pcnst) ! Calc tendencies? + real(r8) :: dTs (state%ncol) + real(r8) :: dUa (state%ncol,pver) + real(r8) :: dVa (state%ncol,pver) + real(r8) :: dTa (state%ncol,pver) + real(r8) :: dQa (state%ncol,pver) + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: kk ! loop index + + ! Set local copies of values + !--------------------------------- + lchnk = state%lchnk + ncol = state%ncol + Zm (:ncol,:) = state%zm (:ncol,:) + Zi (:ncol,1:pver) = state%zi (:ncol,1:pver) + T (:ncol,:) = state%T (:ncol,:) + U (:ncol,:) = state%U (:ncol,:) + V (:ncol,:) = state%V (:ncol,:) + qv (:ncol,:) = state%Q (:ncol,:,1) + + ! Initialize individual parameterization tendencies + !----------------------------------------------------- + lq = .false. + lq(1) = .true. + call physics_ptend_init(ptend,state%psetcols,'Frierson pbl_tend', & + ls=.true., lu=.true., lv=.true., lq=lq) + + ! Call the Selected PBL routine + !-------------------------------------------------------- + Tsfc(:ncol) = Tsurf(:ncol,lchnk) + Qsfc(:ncol) = Qsurf(:ncol,lchnk) + if(PBL_OPT == PBL_FRIERSON) then + ! Call Frierson PBL scheme + !-------------------------------------------------- + call frierson_pbl(ncol, pver, ztodt,state%pmid (:ncol,:), & + state%pint (:ncol,:), & + Zm(:ncol,:), & + Zi(:ncol,:), & + state%ps(:ncol) , & + Tsfc(:ncol) , & + Qsfc(:ncol) , & + T(:ncol,:), & + U(:ncol,:), & + V(:ncol,:), & + qv(:ncol,:), & + Fsolar(:ncol,lchnk), & + Fdown(:ncol,lchnk), & + Cdrag(:ncol) , & + Km(:ncol,:), & + Ke(:ncol,:), & + VSE(:ncol,:), & + Z_pbl(:ncol) , & + Rf(:ncol,:), & + dqdt_vdiff(:ncol,:), & + dtdt_vdiff(:ncol,:), & + dudt_vdiff(:ncol,:), & + dvdt_vdiff(:ncol,:), & + LHflux(:ncol,lchnk), & + SHflux(:ncol,lchnk), & + TUflux(:ncol,lchnk), & + TVflux(:ncol,lchnk) ) + elseif(PBL_OPT == PBL_USER) then + ! Call USER implemented routine in frierson module + !-------------------------------------------------- + call frierson_pbl_USER(ncol, pver, ztodt,state%pmid (:ncol,:), & + state%pint (:ncol,:), & + Zm(:ncol,:), & + Zi(:ncol,:), & + state%ps(:ncol) , & + Tsfc(:ncol) , & + Qsfc(:ncol) , & + T(:ncol,:), & + U(:ncol,:), & + V(:ncol,:), & + qv(:ncol,:), & + Fsolar(:ncol,lchnk), & + Fdown(:ncol,lchnk), & + Cdrag(:ncol) , & + Km(:ncol,:), & + Ke(:ncol,:), & + VSE(:ncol,:), & + Z_pbl(:ncol) , & + Rf(:ncol,:), & + dqdt_vdiff(:ncol,:), & + dtdt_vdiff(:ncol,:), & + dudt_vdiff(:ncol,:), & + dvdt_vdiff(:ncol,:), & + LHflux(:ncol,lchnk), & + SHflux(:ncol,lchnk), & + TUflux(:ncol,lchnk), & + TVflux(:ncol,lchnk) ) + else + ! ERROR: Unknown PBL_OPT value + !------------------------------------- + write(iulog,*) 'ERROR: unknown PBL_OPT=',PBL_OPT + call endrun('frierson_pbl_tend() PBL_OPT ERROR') + endif + Tsurf(:ncol,lchnk) = Tsfc (:ncol) + Qsurf(:ncol,lchnk) = Qsfc (:ncol) + Cd (:ncol,lchnk) = Cdrag(:ncol) + + ! Back out tendencies from updated fields + !----------------------------------------- + do kk = 1, pver + ptend%s(:ncol,kk ) = (T (:,kk)-state%T(:ncol,kk ))/ztodt*cpair + ptend%u(:ncol,kk ) = (U (:,kk)-state%U(:ncol,kk ))/ztodt + ptend%v(:ncol,kk ) = (V (:,kk)-state%V(:ncol,kk ))/ztodt + ptend%q(:ncol,kk,1) = (qv(:,kk)-state%q(:ncol,kk,1))/ztodt + end do + + ! Archive diagnostic fields + !---------------------------- + call outfld('gray_Tsurf' ,Tsurf(:ncol,lchnk) ,ncol,lchnk) + call outfld('gray_Qsurf' ,Qsurf(:ncol,lchnk) ,ncol,lchnk) + call outfld('gray_Cdrag' ,Cd (:ncol,lchnk) ,ncol,lchnk) + call outfld('gray_Zpbl' ,Z_pbl ,ncol,lchnk) ! + call outfld('gray_KVH' ,Ke ,ncol,lchnk) ! Eddy diffusivity (heat and moisture,m2/s) + call outfld('gray_KVM' ,Km ,ncol,lchnk) ! Eddy diffusivity (momentum, m2/s) + call outfld('gray_VSE' ,VSE ,ncol,lchnk) ! Virtual Dry Static Energy divided by Cp (K) + call outfld('gray_Zm' ,Zm ,ncol,lchnk) ! + call outfld('gray_Rf' ,Rf ,ncol,lchnk) ! + call outfld('gray_DTV' ,dtdt_vdiff ,ncol,lchnk) ! PBL + surface flux T tendency (K/s) + call outfld('gray_DUV' ,dudt_vdiff ,ncol,lchnk) ! PBL u tendency (m/s2) + call outfld('gray_DVV' ,dvdt_vdiff ,ncol,lchnk) ! PBL v tendency (m/s2) + call outfld('gray_VD01' ,dqdt_vdiff ,ncol,lchnk) ! PBL + surface flux Q tendency (kg/kg/s) + call outfld('gray_SHflux',SHflux(:ncol,lchnk),ncol,lchnk) ! Sensible Heat Flux + call outfld('gray_LHflux',LHflux(:ncol,lchnk),ncol,lchnk) ! Latent Heat Flux + call outfld('gray_TauU' ,TUflux(:ncol,lchnk),ncol,lchnk) ! U Surface Stress + call outfld('gray_TauV' ,TVflux(:ncol,lchnk),ncol,lchnk) ! V Surface Stress + + end subroutine frierson_pbl_tend + !============================================================================ + + + !============================================================================ + subroutine frierson_radiative_tend(state, ptend, ztodt,cam_in,cam_out) + ! + ! frierson_radiative_tend: Run the radiative process + !========================================================================= + use physics_types,only: physics_state, physics_ptend + use physics_types,only: physics_ptend_init + use phys_grid, only: get_rlat_all_p + use frierson, only: frierson_radiation,frierson_radiation_USER + ! + ! Passed Variables + !------------------ + type(physics_state),intent(in) :: state + real(r8) ,intent(in) :: ztodt + type(physics_ptend),intent(out) :: ptend + type(cam_in_t), intent(inout):: cam_in + type(cam_out_t), intent(inout):: cam_out + ! + ! Local Values + !--------------- + real(r8):: T (state%ncol,pver) ! T temporary + real(r8):: qv (state%ncol,pver) ! Q temporary + real(r8):: dtdt_heating(state%ncol,pver) ! Longwave heating tendency K/s + real(r8):: dtdt_solar (state%ncol,pver) ! Shortwave heating tendency K/s + real(r8):: Tsfc (state%ncol) ! Surface T + real(r8):: Qsfc (state%ncol) ! Surface Q (saturated) + logical :: lq(pcnst) ! Calc tendencies? + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer :: k ! loop index + + ! Copy to local values + !------------------------------------------------- + lchnk = state%lchnk + ncol = state%ncol + T (:ncol,:) = state%T(:ncol,:) + qv (:ncol,:) = state%Q(:ncol,:,1) + + !-------------------------------------- + Tsfc(:ncol) = Tsurf(:ncol,lchnk) + Qsfc(:ncol) = Qsurf(:ncol,lchnk) + + ! initialize individual parameterization tendencies + !--------------------------------------------------- + lq(:) = .false. + call physics_ptend_init(ptend, state%psetcols, 'Frierson radiative_tend', & + ls=.true., lu=.false., lv=.false., lq=lq) + + ! Call the Selected radiative routine + !-------------------------------------------------------- + if(RADIATION_OPT == RADIATION_FRIERSON) then + call frierson_radiation(ncol,pver,ztodt,clat(:ncol,lchnk), & + state%pint(:ncol,:), & + state%pmid(:ncol,:), & + state%ps(:ncol), & + Tsfc(:ncol), & + Qsfc(:ncol), & + T(:ncol,:), & + qv(:ncol,:), & + dtdt_heating(:ncol,:), & + Fsolar(:ncol,lchnk), & + Fup(:ncol,lchnk), & + Fdown(:ncol,lchnk), & + Fup_toa(:ncol,lchnk), & + Fdown_toa(:ncol,lchnk) ) + dtdt_solar(:ncol,:) = 0._r8 + elseif(RADIATION_OPT == RADIATION_USER) then + call frierson_radiation_USER(ncol,pver,ztodt,clat(:ncol,lchnk), & + state%pint(:ncol,:), & + state%pmid(:ncol,:), & + state%ps(:ncol), & + Tsfc(:ncol), & + Qsfc(:ncol), & + T(:ncol,:), & + qv(:ncol,:), & + dtdt_heating(:ncol,:), & + Fsolar(:ncol,lchnk), & + Fup(:ncol,lchnk), & + Fdown(:ncol,lchnk), & + Fup_toa(:ncol,lchnk), & + Fdown_toa(:ncol,lchnk) ) + dtdt_solar(:ncol,:) = 0._r8 + else + ! ERROR: Unknown RADIATION_OPT value + !------------------------------------- + write(iulog,*) 'ERROR: unknown RADIATION_OPT=',RADIATION_OPT + call endrun('frierson_pbl_tend() RADIATION_OPT ERROR') + endif + + Fnet (:ncol,lchnk) = Fup(:ncol,lchnk) - Fdown (:ncol,lchnk) + Fnet_toa (:ncol,lchnk) = Fup_toa(:ncol,lchnk) - Fdown_toa (:ncol,lchnk) + + ! Copy downward LW radiative heating values to cam_out% + !--------------------------------------------------------- + cam_out%flwds(:ncol) = Fdown (:ncol,lchnk) + cam_out%netsw(:ncol) = Fsolar(:ncol,lchnk) + cam_out%sols (:ncol) = Fsolar(:ncol,lchnk) + cam_out%solsd(:ncol) = Fsolar(:ncol,lchnk) + cam_out%soll (:ncol) = Fsolar(:ncol,lchnk) + cam_out%solld(:ncol) = Fsolar(:ncol,lchnk) + + ! Back out tendencies from updated T field + !-------------------------------------------- + do k = 1, pver + ptend%s(:ncol,k) = (T(:,k)-state%T(:ncol,k))/ztodt*cpair + end do + + ! Archive T tendency from temperature relaxation (mimics radiation, K/s) + !----------------------------------------------------------------------- + call outfld('gray_QRL' ,dtdt_heating, ncol,lchnk) + call outfld('gray_QRS' ,dtdt_solar , ncol,lchnk) + call outfld('gray_SWflux',Fsolar(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LUflux',Fup(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LDflux',Fdown(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LWflux',Fnet(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LUflux_TOA',Fup_toa(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LDflux_TOA',Fdown_toa(:ncol,lchnk) , ncol,lchnk) + call outfld('gray_LWflux_TOA',Fnet_toa(:ncol,lchnk) , ncol,lchnk) + + end subroutine frierson_radiative_tend + !============================================================================ + + + !======================================================================= + subroutine frierson_surface_init(ncol, clat, PS, Tsfc, Qsfc) + ! + ! + !========================================================================== + ! + ! Passed variables + !-------------------- + integer ,intent(in) :: ncol + real(r8),intent(in) :: clat (ncol) + real(r8),intent(in) :: PS (ncol) + real(r8),intent(out):: Tsfc(ncol) + real(r8),intent(out):: Qsfc(ncol) + ! + ! Local values + !-------------- + integer :: ii + real(r8):: T_width + + ! set SST profile + !------------------ + T_width = frierson_Twidth*pi/180.0_r8 + do ii = 1, ncol + Tsfc(ii) = frierson_Tmin + frierson_Tdlt*exp(-((clat(ii)/T_width)**2)/2.0_r8) + Qsfc(ii) = epsilo*frierson_E0/PS(ii) & + *exp(-latvap/rh2o*((1._r8/Tsfc(ii))-1._r8/frierson_T0)) + end do + + end subroutine frierson_surface_init + !======================================================================= + + + !======================================================================= + subroutine frierson_restart_init(File,hdimids,hdimcnt) + ! + ! frierson_restart_init: + !========================================================================== + ! + ! Passed variables + !-------------------- + type(file_desc_t),intent(inout):: File + integer ,intent(in) :: hdimcnt + integer ,intent(in) :: hdimids(1:hdimcnt) + ! + ! Local values + !-------------- + integer:: ierr + + ierr = pio_def_var(File,'Frierson_Tsfc',pio_double, hdimids, Tsurf_desc) + if (ierr /= 0) then + call endrun('frierson_restart_init: ERROR defining Frierson_Tsfc') + end if + + ierr = pio_def_var(File,'Frierson_Qsfc',pio_double, hdimids, Qsurf_desc) + if (ierr /= 0) then + call endrun('frierson_restart_init: ERROR defining Frierson_Qsfc') + end if + + end subroutine frierson_restart_init + !======================================================================= + + + !======================================================================= + subroutine frierson_restart_write(File) + ! + ! frierson_restart_write: + !========================================================================== + ! + ! Passed variables + !-------------------- + type(file_desc_t),intent(inout):: File + ! + ! Local values + !-------------- + type(io_desc_t),pointer:: iodesc + integer:: dims(3),gdims(3),nhdims + integer:: physgrid + integer:: ierr + + ! Get the iodesc for write calls + !--------------------------------- + dims(1) = pcols + dims(2) = endchunk - begchunk + 1 + physgrid = cam_grid_id('physgrid') + call cam_grid_dimensions(physgrid, gdims(1:2), nhdims) + call cam_grid_get_decomp(physgrid, dims(1:2), gdims(1:nhdims), pio_double, iodesc) + + ! Write Surface values + !--------------------- + call pio_write_darray(File, Tsurf_desc, iodesc, Tsurf, ierr) + if (ierr /= 0) then + call endrun('frierson_restart_write: ERROR writing Tsurf') + end if + + call pio_write_darray(File, Qsurf_desc, iodesc, Qsurf, ierr) + if (ierr /= 0) then + call endrun('frierson_restart_write: ERROR writing Qsurf') + end if + + end subroutine frierson_restart_write + !======================================================================= + + + !======================================================================= + subroutine frierson_restart_read(File) + ! + ! frierson_restart_read: + !========================================================================== + use error_messages,only: alloc_err + ! + ! Passed variables + !-------------------- + type(file_desc_t),intent(inout):: File + ! + ! Local values + !-------------- + type( io_desc_t),pointer:: iodesc + type(var_desc_t) :: vardesc + integer:: dims(3),gdims(3),nhdims + integer:: physgrid + integer:: ierr + + ! Allocate space for the restart fields + !----------------------------------------- + allocate(Tsurf (pcols,begchunk:endchunk),stat=ierr) + call alloc_err(ierr,'Frierson RESTART','Tsurf' ,pcols*(endchunk-begchunk+1)) + allocate(Qsurf (pcols,begchunk:endchunk) ,stat=ierr) + call alloc_err(ierr,'Frierson RESTART','Qsurf' ,pcols*(endchunk-begchunk+1)) + + ! Get the iodesc for read calls + !--------------------------------- + dims(1) = pcols + dims(2) = endchunk - begchunk + 1 + physgrid = cam_grid_id('physgrid') + call cam_grid_dimensions(physgrid, gdims(1:2), nhdims) + call cam_grid_get_decomp(physgrid, dims(1:2), gdims(1:nhdims), pio_double, iodesc) + + ! Read Surface values + !--------------------- + ierr = pio_inq_varid(File,'Frierson_Tsfc',vardesc) + if (ierr /= 0) then + call endrun('frierson_restart_read: ERROR PIO unable to find variable Frierson_Tsfc') + end if + + call pio_read_darray(File, vardesc, iodesc, Tsurf, ierr) + if (ierr /= 0) then + call endrun('frierson_restart_read: ERROR PIO unable to read variable Tsurf') + end if + + ierr = pio_inq_varid(File,'Frierson_Qsfc',vardesc) + if (ierr /= 0) then + call endrun('frierson_restart_read: ERROR PIO unable to find variable Frierson_Qsfc') + end if + + call pio_read_darray(File, vardesc, iodesc, Qsurf, ierr) + if (ierr /= 0) then + call endrun('frierson_restart_read: ERROR PIO unable to read variable Qsurf') + end if + + end subroutine frierson_restart_read + !======================================================================= + +end module frierson_cam + diff --git a/src/physics/simple/physpkg.F90 b/src/physics/simple/physpkg.F90 index 9c1e4c61bf..31c41def58 100644 --- a/src/physics/simple/physpkg.F90 +++ b/src/physics/simple/physpkg.F90 @@ -19,7 +19,7 @@ module physpkg use camsrfexch, only: cam_out_t, cam_in_t, cam_export ! Note: ideal_phys is true for Held-Suarez (1994) physics - use cam_control_mod, only: moist_physics, adiabatic, ideal_phys, kessler_phys, tj2016_phys + use cam_control_mod, only: moist_physics, adiabatic, ideal_phys, kessler_phys, tj2016_phys, frierson_phys use phys_control, only: phys_getopts use perf_mod, only: t_barrierf, t_startf, t_stopf, t_adj_detailf use cam_logfile, only: iulog @@ -80,6 +80,7 @@ subroutine phys_register use check_energy, only: check_energy_register use kessler_cam, only: kessler_register use tj2016_cam, only: thatcher_jablonowski_register + use frierson_cam, only: frierson_register !---------------------------Local variables----------------------------- ! @@ -111,6 +112,8 @@ subroutine phys_register call kessler_register() else if (tj2016_phys) then call thatcher_jablonowski_register() + else if (frierson_phys) then + call frierson_register() end if ! Fields for physics package diagnostics @@ -197,10 +200,12 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) use held_suarez_cam, only: held_suarez_init use kessler_cam, only: kessler_cam_init use tj2016_cam, only: thatcher_jablonowski_init + use frierson_cam, only: frierson_init use tracers, only: tracers_init use wv_saturation, only: wv_sat_init use phys_debug_util, only: phys_debug_init use qneg_module, only: qneg_init + use nudging, only: Nudge_Model, nudging_init use cam_snapshot, only: cam_snapshot_init use cam_budget, only: cam_budget_init @@ -243,7 +248,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) ! wv_saturation is relatively independent of everything else and ! low level, so init it early. Must at least do this before radiation. - if (kessler_phys .or. tj2016_phys) then + if (kessler_phys .or. tj2016_phys .or. frierson_phys) then call wv_sat_init() end if @@ -259,8 +264,14 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) call kessler_cam_init(pbuf2d) else if (tj2016_phys) then call thatcher_jablonowski_init(pbuf2d) + else if (frierson_phys) then + call frierson_init(phys_state,pbuf2d) end if + ! Initialize Nudging Parameters + !-------------------------------- + if(Nudge_Model) call nudging_init + if (chem_is_active()) then ! Prognostic chemistry. call chem_init(phys_state,pbuf2d) @@ -476,6 +487,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) use constituents, only: cnst_get_ind, pcnst use cam_diagnostics, only: diag_phys_tend_writeout, diag_surf use tj2016_cam, only: thatcher_jablonowski_sfc_pbl_hs_tend + use frierson_cam, only: frierson_pbl_tend use dycore, only: dycore_is use check_energy, only: tot_energy_phys use cam_history, only: hist_fld_active @@ -483,6 +495,10 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) use cam_budget, only: thermo_budget_history use dyn_tests_utils, only: vc_dycore, vc_height, vc_dry_pressure use air_composition, only: cpairv, cp_or_cv_dycore + use time_manager, only: get_nstep + use nudging, only: Nudge_Model, Nudge_ON, nudging_timestep_tend + use check_energy, only: check_energy_chng + ! Arguments ! real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) @@ -495,6 +511,9 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) !---------------------------Local workspace----------------------------- + integer :: nstep ! current timestep number + real(r8):: zero(pcols) ! array of zeros + type(physics_ptend) :: ptend ! indivdual parameterization tendencies real(r8) :: tmp_q(pcols, pver) real(r8) :: tmp_cldliq(pcols, pver) @@ -518,6 +537,10 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) real(r8) :: scaling(pcols,pver) !-------------------------------------------------------------------------- + ! get nstep and zero array for energy checker + zero = 0._r8 + nstep = get_nstep() + ! number of active atmospheric columns ncol = state%ncol lchnk = state%lchnk @@ -555,14 +578,28 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) call physics_update(state, ptend, ztodt, tend) end if + if (frierson_phys) then + ! Update surface, PBL + call frierson_pbl_tend(state, ptend, ztodt, cam_in) + call physics_update(state, ptend, ztodt, tend) + end if + + ! Update Nudging values, if needed + !---------------------------------- + if (Nudge_Model .and. Nudge_ON) then + call nudging_timestep_tend(state,ptend) + call physics_update(state, ptend, ztodt, tend) + call check_energy_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) + endif + call tot_energy_phys(state, 'phAP') call tot_energy_phys(state, 'dyAP',vc=vc_dycore) - + ! FV: convert dry-type mixing ratios to moist here because ! physics_dme_adjust assumes moist. This is done in p_d_coupling for ! other dynamics. Bundy, Feb 2004. ! - moist_mixing_ratio_dycore = dycore_is('LR').or. dycore_is('FV3') + moist_mixing_ratio_dycore = dycore_is('LR').or. dycore_is('FV3') ! ! update cp/cv for energy computation based in updated water variables ! @@ -608,7 +645,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) else ! ! for moist-mixing ratio based dycores - ! + ! ! Note: this operation will NOT be reverted with set_wet_to_dry after set_dry_to_wet call ! call set_dry_to_wet(state) @@ -619,7 +656,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) if (vc_dycore == vc_height.or.vc_dycore == vc_dry_pressure) then ! ! MPAS and SE specific scaling of temperature for enforcing energy consistency - ! (and to make sure that temperature dependent diagnostic tendencies + ! (and to make sure that temperature dependent diagnostic tendencies ! are computed correctly; e.g. dtcore) ! scaling(1:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk) @@ -630,7 +667,7 @@ subroutine tphysac (ztodt, cam_in, cam_out, state, tend, pbuf) ! else: do nothing for dycores with energy consistent with CAM physics ! end if - + else tmp_q (:ncol,:pver) = 0.0_r8 tmp_cldliq(:ncol,:pver) = 0.0_r8 @@ -692,6 +729,8 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) use held_suarez_cam, only: held_suarez_tend use kessler_cam, only: kessler_tend use tj2016_cam, only: thatcher_jablonowski_precip_tend + use frierson_cam, only: frierson_condensate_tend + use frierson_cam, only: frierson_radiative_tend use dycore, only: dycore_is use cam_snapshot_common,only: cam_snapshot_all_outfld use cam_snapshot_common,only: cam_snapshot_ptend_outfld @@ -886,6 +925,37 @@ subroutine tphysbc (ztodt, state, tend, pbuf, cam_out, cam_in ) if (trim(cam_take_snapshot_after) == "thatcher_jablonowski_precip_tend") then call cam_snapshot_all_outfld(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf) end if + else if (frierson_phys) then + ! Compute the large-scale precipitation + !---------------------------------------- + if (trim(cam_take_snapshot_before) == "frierson_condensate_tend") then + call cam_snapshot_all_outfld(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf) + end if + call frierson_condensate_tend(state, ptend, ztodt, pbuf) + if ( (trim(cam_take_snapshot_after) == "frierson_condensate_tend") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + if (trim(cam_take_snapshot_after) == "frierson_condensate_tend") then + call cam_snapshot_all_outfld(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf) + end if + + ! Compute the radiative tendencies + !----------------------------------- + if (trim(cam_take_snapshot_before) == "frierson_radiative_tend") then + call cam_snapshot_all_outfld(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf) + end if + call frierson_radiative_tend(state, ptend, ztodt, cam_in, cam_out) + if ( (trim(cam_take_snapshot_after) == "frierson_radiative_tend") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + if (trim(cam_take_snapshot_after) == "frierson_radiative_tend") then + call cam_snapshot_all_outfld(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf) + end if + end if ! Can't turn on conservation error messages unless the appropriate heat @@ -945,6 +1015,7 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) !-------------------------------------------------------------------------- use physics_types, only: physics_state use physics_buffer, only: physics_buffer_desc + use nudging, only: Nudge_Model, nudging_timestep_init implicit none @@ -956,6 +1027,10 @@ subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) !-------------------------------------------------------------------------- + ! Update Nudging values, if needed + !---------------------------------- + if(Nudge_Model) call nudging_timestep_init(phys_state) + end subroutine phys_timestep_init !====================================================================================== diff --git a/src/physics/simple/restart_physics.F90 b/src/physics/simple/restart_physics.F90 index ef8f8795ef..fb40c5921b 100644 --- a/src/physics/simple/restart_physics.F90 +++ b/src/physics/simple/restart_physics.F90 @@ -13,6 +13,9 @@ module restart_physics pio_inq_varid, pio_def_var, pio_def_dim, & pio_put_var, pio_get_var + use cam_control_mod, only: frierson_phys + use frierson_cam,only: frierson_restart_init, frierson_restart_write, frierson_restart_read + implicit none private save @@ -59,6 +62,10 @@ subroutine init_restart_physics ( File, pbuf2d) call pbuf_init_restart(File, pbuf2d) + if (frierson_phys) then + call frierson_restart_init(File,hdimids,hdimcnt) + end if + end subroutine init_restart_physics subroutine write_restart_physics (File, cam_in, cam_out, pbuf2d) @@ -85,6 +92,10 @@ subroutine write_restart_physics (File, cam_in, cam_out, pbuf2d) ! Physics buffer call pbuf_write_restart(File, pbuf2d) + if (frierson_phys) then + call frierson_restart_write(File) + end if + end subroutine write_restart_physics !####################################################################### @@ -110,6 +121,9 @@ subroutine read_restart_physics(File, cam_in, cam_out, pbuf2d) call pbuf_read_restart(File, pbuf2d) + if (frierson_phys) then + call frierson_restart_read(File) + end if end subroutine read_restart_physics end module restart_physics