diff --git a/src/aero_model.F90 b/src/aero_model.F90 index fa5e51f..a0d8d11 100644 --- a/src/aero_model.F90 +++ b/src/aero_model.F90 @@ -43,7 +43,7 @@ module aero_model use oslo_aero_share, only: lifeCycleNumberMedianRadius, rhopart, lifeCycleSigma use oslo_aero_share, only: l_so4_a2, l_bc_n, l_bc_ax use oslo_aero_share, only: MODE_IDX_BC_NUC, MODE_IDX_BC_EXT_AC - use oslo_aero_control, only: oslo_aero_ctl_readnl + use oslo_aero_control, only: oslo_aero_ctl_readnl, use_aerocom use oslo_aero_depos, only: oslo_aero_depos_init use oslo_aero_depos, only: oslo_aero_depos_dry, oslo_aero_depos_wet, oslo_aero_wetdep_init use oslo_aero_coag, only: initializeCoagulation, coagtend, clcoag @@ -153,10 +153,10 @@ subroutine aero_model_init( pbuf2d ) call init_interp_constants() ! table initialization constants call initopt() ! table initialization call initlogn() ! table initialization -#ifdef AEROCOM - call initdry() ! table initialization - call initaeropt() ! table initialization -#endif + if (use_aerocom) then + call initdry() ! table initialization + call initaeropt() ! table initialization + end if call initializeCondensation() call oslo_aero_ocean_init() call oslo_aero_depos_init(pbuf2d) diff --git a/src/oslo_aero_control.F90 b/src/oslo_aero_control.F90 index 659a245..b5cef75 100644 --- a/src/oslo_aero_control.F90 +++ b/src/oslo_aero_control.F90 @@ -23,8 +23,11 @@ module oslo_aero_control integer, parameter :: unset_int = huge(1) integer, parameter, public :: dir_string_length=256 - ! Namelist variables: - real(r8) :: volc_fraction_coarse = 0.0_r8 !Fraction of volcanic aerosols in coarse mode + ! Public Namelist variables: + logical, public, protected :: use_aerocom = .false. ! If true, turn on aerocom output + + ! Private Namelist variables: + real(r8) :: volc_fraction_coarse = 0.0_r8 !Fraction of volcanic aerosols in coarse mode character(len=dir_string_length) :: aerotab_table_dir = unset_str ! DMS/Ocean namelist variables @@ -53,7 +56,8 @@ subroutine oslo_aero_ctl_readnl(nlfile) namelist /oslo_ctl_nl/ volc_fraction_coarse, aerotab_table_dir, dms_source, & dms_source_type, opom_source, opom_source_type, & - ocean_filename, ocean_filepath, dms_cycle_year, opom_cycle_year + ocean_filename, ocean_filepath, dms_cycle_year, opom_cycle_year, & + use_aerocom !----------------------------------------------------------------------------- if (masterproc) then @@ -69,6 +73,9 @@ subroutine oslo_aero_ctl_readnl(nlfile) end if ! Broadcast namelist variables + call mpi_bcast(use_aerocom, 1 , mpi_logical, mstrid, mpicom, ierr) + if (ierr /= mpi_success) call endrun(subname//" mpi_bcast: volc_fraction_coarse") + call mpi_bcast(volc_fraction_coarse, 1 , mpi_real8, mstrid, mpicom, ierr) if (ierr /= mpi_success) call endrun(subname//" mpi_bcast: volc_fraction_coarse") call mpi_bcast(aerotab_table_dir, len(aerotab_table_dir) , mpi_character, mstrid, mpicom, ierr) diff --git a/src/oslo_aero_diagnostics.F90 b/src/oslo_aero_diagnostics.F90 new file mode 100644 index 0000000..46d2e7f --- /dev/null +++ b/src/oslo_aero_diagnostics.F90 @@ -0,0 +1,362 @@ +module oslo_aero_diagnostics + +use shr_kind_mod , only: r8 => shr_kind_r8 +use cam_history , only: addfld, add_default, horiz_only +use oslo_aero_control , only: use_aerocom +use oslo_aero_share , only: nbmodes + +implicit none +private + +public :: oslo_aero_diagnostics_init + +contains + + subroutine oslo_aero_diagnostics_init() + + integer :: indx + character(len=20) :: varname + + call addfld ('AODVIS ',horiz_only, 'A','unitless' ,'Aerosol optical depth at 0.442-0.625um') ! CAM4-Oslo: 0.35-0.64um + call addfld ('ABSVIS ',horiz_only, 'A','unitless' ,'Aerosol absorptive optical depth at 0.442-0.625um') ! CAM4-Oslo: 0.35-0.64um + call addfld ('AODVVOLC ',horiz_only, 'A','unitless' ,'CMIP6 volcanic aerosol optical depth at 0.442-0.625um') ! CAM4-Oslo: 0.35-0.64um + call addfld ('ABSVVOLC ',horiz_only, 'A','unitless' ,'CMIP6 volcanic aerosol absorptive optical depth at 0.442-0.625um') ! CAM4-Oslo: 0.35-0.64um + call addfld ('CAODVIS ',horiz_only, 'A','unitless' ,'Clear air aerosol optical depth') + call addfld ('CABSVIS ',horiz_only, 'A','unitless' ,'Clear air aerosol absorptive optical depth') + call addfld ('CLDFREE ',horiz_only, 'A','unitless' ,'Cloud free fraction wrt CAODVIS and CABSVIS') + call addfld ('DAYFOC ',horiz_only, 'A','unitless' ,'Daylight fraction') + call addfld ('N_AER ',(/'lev'/), 'A', 'unitless','Aerosol number concentration') + call addfld ('SSAVIS ',(/'lev'/), 'A','unitless' ,'Aerosol single scattering albedo in visible wavelength band') + call addfld ('ASYMMVIS',(/'lev'/), 'A','unitless' ,'Aerosol assymetry factor in visible wavelength band') + call addfld ('EXTVIS ',(/'lev'/), 'A','1/km ' ,'Aerosol extinction') + call addfld ('BVISVOLC ',(/'lev'/), 'A','1/km ' ,'CMIP6 volcanic aerosol extinction at 0.442-0.625um') + + call addfld ('FSNT_DRF',horiz_only, 'A','W/m^2 ','Total column absorbed solar flux (DIRind)') + call addfld ('FSNTCDRF',horiz_only, 'A','W/m^2 ','Clear sky total column absorbed solar flux (DIRind)' ) + call addfld ('FSNS_DRF',horiz_only, 'A','W/m^2 ','Surface absorbed solar flux (DIRind)' ) + call addfld ('FSNSCDRF',horiz_only, 'A','W/m^2 ','Clear sky surface absorbed solar flux (DIRind)' ) + call addfld ('QRS_DRF ',(/'lev'/), 'A','K/s ','Solar heating rate (DIRind)') + call addfld ('QRSC_DRF',(/'lev'/), 'A','K/s ','Clearsky solar heating rate (DIRind)' ) + call addfld ('FLNT_DRF',horiz_only, 'A','W/m^2 ','Total column longwave flux (DIRind)' ) + call addfld ('FLNTCDRF',horiz_only, 'A','W/m^2 ','Clear sky total column longwave flux (DIRind)' ) + call addfld ('FSUTADRF',horiz_only, 'A','W/m^2 ','SW upwelling flux at TOA') + call addfld ('FSDS_DRF',horiz_only, 'A','W/m^2 ','SW downelling flux at surface') + call addfld ('FSUS_DRF',horiz_only, 'A','W/m^2 ','SW upwelling flux at surface') + call addfld ('FSDSCDRF',horiz_only, 'A','W/m^2 ','SW downwelling clear sky flux at surface') + call addfld ('FLUS ',horiz_only, 'A','W/m^2 ','LW surface upwelling flux') + + call add_default ('AODVIS ', 1, ' ') + call add_default ('ABSVIS ', 1, ' ') + call add_default ('AODVVOLC', 1, ' ') + call add_default ('ABSVVOLC', 1, ' ') + call add_default ('DAYFOC ', 1, ' ') + call add_default ('CAODVIS ', 1, ' ') + call add_default ('CABSVIS ', 1, ' ') + call add_default ('CLDFREE ', 1, ' ') + call add_default ('N_AER ', 1, ' ') + !call add_default ('N_AERORG', 1, ' ') + call add_default ('SSAVIS ', 1, ' ') + call add_default ('ASYMMVIS', 1, ' ') + call add_default ('EXTVIS ', 1, ' ') + call add_default ('BVISVOLC', 1, ' ') + call add_default ('FSNT_DRF', 1, ' ') + call add_default ('FSNTCDRF', 1, ' ') + call add_default ('FSNS_DRF', 1, ' ') + call add_default ('FSNSCDRF', 1, ' ') + call add_default ('QRS_DRF ', 1, ' ') + call add_default ('QRSC_DRF', 1, ' ') + call add_default ('FLNT_DRF', 1, ' ') + call add_default ('FLNTCDRF', 1, ' ') + call add_default ('FSUTADRF', 1, ' ') + call add_default ('FSDS_DRF', 1, ' ') + call add_default ('FSUS_DRF', 1, ' ') + call add_default ('FSDSCDRF', 1, ' ') + call add_default ('FLUS ', 1, ' ') + + if (use_aerocom) then + call addfld ('AKCXS ',horiz_only, 'A','mg/m2 ','Scheme excess aerosol mass burden') + call addfld ('PMTOT ',horiz_only, 'A','ug/m3 ','Aerosol PM, all sizes') + call addfld ('PM25 ',horiz_only, 'A','ug/m3 ','Aerosol PM2.5') + + call addfld ('PM2P5 ',(/'lev'/), 'A','ug/m3 ','3D aerosol PM2.5') + call addfld ('MMRPM2P5',(/'lev'/), 'A','kg/kg ','3D aerosol PM2.5 mass mixing ratio') + call addfld ('MMRPM1 ',(/'lev'/), 'A','kg/kg ','3D aerosol PM1.0 mass mixing ratio') + + call addfld ('MMRPM2P5_SRF' ,horiz_only, 'A','kg/kg ','Aerosol PM2.5 mass mixing ratio in bottom layer') + call addfld ('GRIDAREA' ,horiz_only, 'A','m2 ','Grid area for 1.9x2.5 horizontal resolution') + call addfld ('DAERH2O ' ,horiz_only, 'A','mg/m2 ','Aerosol water load') + + call addfld ('MMR_AH2O',(/'lev'/), 'A','kg/kg ','Aerosol water mmr') + call addfld ('ECDRYAER',(/'lev'/), 'A','kg/kg ','Dry aerosol extinction at 550nm') + call addfld ('ABSDRYAE',(/'lev'/), 'A','m-1 ','Dry aerosol absorption at 550nm') + call addfld ('ECDRY440',(/'lev'/), 'A','m-1 ','Dry aerosol extinction at 440nm') + call addfld ('ABSDR440',(/'lev'/), 'A','m-1 ','Dry aerosol absorption at 440nm') + call addfld ('ECDRY870',(/'lev'/), 'A','m-1 ','Dry aerosol extinction at 870nm') + call addfld ('ABSDR870',(/'lev'/), 'A','m-1 ','Dry aerosol absorption at 870nm') + call addfld ('ASYMMDRY',(/'lev'/), 'A','unitless','Dry asymmetry factor in visible wavelength band') + call addfld ('ECDRYLT1',(/'lev'/), 'A','m-1 ','Dry aerosol extinction at 550nm lt05') + call addfld ('ABSDRYBC',(/'lev'/), 'A','m-1 ','Dry BC absorption at 550nm') + call addfld ('ABSDRYOC',(/'lev'/), 'A','m-1 ','Dry OC absorption at 550nm') + call addfld ('ABSDRYSU',(/'lev'/), 'A','m-1 ','Dry sulfate absorption at 550nm') + call addfld ('ABSDRYSS',(/'lev'/), 'A','m-1 ','Dry sea-salt absorption at 550nm') + call addfld ('ABSDRYDU',(/'lev'/), 'A','m-1 ','Dry dust absorption at 550nm') + + call addfld ('OD550DRY',horiz_only, 'A','unitless','Dry aerosol optical depth at 550nm') + call addfld ('AB550DRY',horiz_only, 'A','unitless','Dry aerosol absorptive optical depth at 550nm') + call addfld ('DERLT05 ',horiz_only, 'A','um ','Effective aerosol dry radius<0.5um') + call addfld ('DERGT05 ',horiz_only, 'A','um ','Effective aerosol dry radius>0.5um') + call addfld ('DER ',horiz_only, 'A','um ','Effective aerosol dry radius') + call addfld ('DOD440 ',horiz_only, 'A','unitless','Aerosol optical depth at 440nm') + call addfld ('ABS440 ',horiz_only, 'A','unitless','Aerosol absorptive optical depth at 440nm') + call addfld ('DOD500 ',horiz_only, 'A','unitless','Aerosol optical depth at 500nm') + call addfld ('ABS500 ',horiz_only, 'A','unitless','Aerosol absorptive optical depth at 500nm') + call addfld ('DOD550 ',horiz_only, 'A','unitless','Aerosol optical depth at 550nm') + call addfld ('ABS550 ',horiz_only, 'A','unitless','Aerosol absorptive optical depth at 550nm') + call addfld ('ABS550AL',horiz_only, 'A','unitless','Alt. aerosol absorptive optical depth at 550nm') + call addfld ('DOD670 ',horiz_only, 'A','unitless','Aerosol optical depth at 670nm') + call addfld ('ABS670 ',horiz_only, 'A','unitless','Aerosol absorptive optical depth at 670nm') + call addfld ('DOD870 ',horiz_only, 'A','unitless','Aerosol optical depth at 870nm') + call addfld ('ABS870 ',horiz_only, 'A','unitless','Aerosol absorptive optical depth at 870nm') + + call addfld ('DLOAD_MI',horiz_only, 'A','mg/m2 ','mineral aerosol load') + call addfld ('DLOAD_SS',horiz_only, 'A','mg/m2 ','sea-salt aerosol load') + call addfld ('DLOAD_S4',horiz_only, 'A','mg/m2 ','sulfate aerosol load') + call addfld ('DLOAD_OC',horiz_only, 'A','mg/m2 ','OC aerosol load') + call addfld ('DLOAD_BC',horiz_only, 'A','mg/m2 ','BC aerosol load') + + call addfld ('LOADBCAC',horiz_only, 'A','mg/m2 ','BC aerosol coag load') + call addfld ('LOADBC0 ',horiz_only, 'A','mg/m2 ','BC aerosol mode 0 load') + call addfld ('LOADBC2 ',horiz_only, 'A','mg/m2 ','BC aerosol mode 2 load') + call addfld ('LOADBC4 ',horiz_only, 'A','mg/m2 ','BC aerosol mode 4 load') + call addfld ('LOADBC12',horiz_only, 'A','mg/m2 ','BC aerosol mode 12 load') + call addfld ('LOADBC14',horiz_only, 'A','mg/m2 ','BC aerosol mode 14 load') + call addfld ('LOADOCAC',horiz_only, 'A','mg/m2 ','OC aerosol coag load') + call addfld ('LOADOC3 ',horiz_only, 'A','mg/m2 ','OC aerosol mode 3 load') + call addfld ('LOADOC4 ',horiz_only, 'A','mg/m2 ','OC aerosol mode 4 load') + call addfld ('LOADOC13',horiz_only, 'A','mg/m2 ','OC aerosol mode 13 load') + call addfld ('LOADOC14',horiz_only, 'A','mg/m2 ','OC aerosol mode 14 load') + ! + call addfld ('EC550AER',(/'lev'/),'A','m-1 ','aerosol extinction coefficient') + call addfld ('ABS550_A',(/'lev'/),'A','m-1 ','aerosol absorption coefficient') + call addfld ('BS550AER',(/'lev'/),'A','m-1 sr-1','aerosol backscatter coefficient') + ! + call addfld ('EC550SO4',(/'lev'/), 'A','m-1 ','SO4 aerosol extinction coefficient') + call addfld ('EC550BC ',(/'lev'/), 'A','m-1 ','BC aerosol extinction coefficient') + call addfld ('EC550POM',(/'lev'/), 'A','m-1 ','POM aerosol extinction coefficient') + call addfld ('EC550SS ',(/'lev'/), 'A','m-1 ','SS aerosol extinction coefficient') + call addfld ('EC550DU ',(/'lev'/), 'A','m-1 ','DU aerosol extinction coefficient') + ! + call addfld ('CDOD440 ' ,horiz_only, 'A','unitless','Clear air Aerosol optical depth at 440nm') + call addfld ('CDOD550 ' ,horiz_only, 'A','unitless','Clear air Aerosol optical depth at 550nm') + call addfld ('CABS550 ' ,horiz_only, 'A','unitless','Clear air Aerosol abs optical depth at 550nm') + call addfld ('CABS550A ' ,horiz_only, 'A','unitless','Clear air Aerosol abs optical depth at 550nm') + call addfld ('CDOD870 ' ,horiz_only, 'A','unitless','Clear air Aerosol optical depth at 870nm') + call addfld ('A550_DU ' ,horiz_only, 'A','unitless', 'mineral abs. aerosol optical depth 550nm') + call addfld ('A550_SS ' ,horiz_only, 'A','unitless','sea-salt abs aerosol optical depth 550nm') + call addfld ('A550_SO4' ,horiz_only, 'A','unitless','SO4 aerosol abs. optical depth 550nm') + call addfld ('A550_POM' ,horiz_only, 'A','unitless', 'OC abs. aerosol optical depth 550nm') + call addfld ('A550_BC ' ,horiz_only, 'A','unitless', 'BC abs. aerosol optical depth 550nm') + call addfld ('D440_DU ' ,horiz_only, 'A','unitless','mineral aerosol optical depth 440nm') + call addfld ('D440_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 440nm') + call addfld ('D440_SO4' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 440nm') + call addfld ('D440_POM' ,horiz_only, 'A','unitless','OC aerosol optical depth 440nm') + call addfld ('D440_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 440nm') + call addfld ('D500_DU ' ,horiz_only, 'A','unitless','mineral aerosol optical depth 500nm') + call addfld ('D500_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 500nm') + call addfld ('D500_SO4' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 500nm') + call addfld ('D500_POM' ,horiz_only, 'A','unitless','OC aerosol optical depth 500nm') + call addfld ('D500_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 500nm') + call addfld ('D550_DU ' ,horiz_only, 'A','unitless','mineral aerosol optical depth 550nm') + call addfld ('D550_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 550nm') + call addfld ('D550_SO4' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 550nm') + call addfld ('D550_POM' ,horiz_only, 'A','unitless','OC aerosol optical depth 550nm') + call addfld ('D550_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 550nm') + call addfld ('D670_DU ' ,horiz_only, 'A','unitless','mineral aerosol optical depth 670nm') + call addfld ('D670_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 670nm') + call addfld ('D670_SO4' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 670nm') + call addfld ('D670_POM' ,horiz_only, 'A','unitless','OC aerosol optical depth 670nm') + call addfld ('D670_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 670nm') + call addfld ('D870_DU ' ,horiz_only, 'A','unitless','mineral aerosol optical depth 870nm') + call addfld ('D870_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 870nm') + call addfld ('D870_SO4' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 870nm') + call addfld ('D870_POM' ,horiz_only, 'A','unitless','OC aerosol optical depth 870nm') + call addfld ('D870_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 870nm') + call addfld ('DLT_DUST' ,horiz_only, 'A','unitless','mineral aerosol optical depth 550nm lt05') + call addfld ('DLT_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 550nm lt05') + call addfld ('DLT_SO4 ' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 550nm lt05') + call addfld ('DLT_POM ' ,horiz_only, 'A','unitless','OC aerosol optical depth 550nm lt05') + call addfld ('DLT_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 550nm lt05') + call addfld ('DGT_DUST' ,horiz_only, 'A','unitless','mineral aerosol optical depth 550nm gt05') + call addfld ('DGT_SS ' ,horiz_only, 'A','unitless','sea-salt aerosol optical depth 550nm gt05') + call addfld ('DGT_SO4 ' ,horiz_only, 'A','unitless','SO4 aerosol optical depth 550nm gt05') + call addfld ('DGT_POM ' ,horiz_only, 'A','unitless','OC aerosol optical depth 550nm gt05') + call addfld ('DGT_BC ' ,horiz_only, 'A','unitless','BC aerosol optical depth 550nm gt05') + call addfld ('AIRMASS ' ,horiz_only, 'A','kg/m2 ','Vertically integrated airmass') + + call addfld ('NNAT_0 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 0 number concentration') + call addfld ('NNAT_1 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 1 number concentration') + call addfld ('NNAT_2 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 2 number concentration') + call addfld ('NNAT_4 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 4 number concentration') + call addfld ('NNAT_5 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 5 number concentration') + call addfld ('NNAT_6 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 6 number concentration') + call addfld ('NNAT_7 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 7 number concentration') + call addfld ('NNAT_8 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 8 number concentration') + call addfld ('NNAT_9 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 9 number concentration') + call addfld ('NNAT_10 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 10 number concentration') + call addfld ('NNAT_12 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 12 number concentration') + call addfld ('NNAT_14 ',(/'lev'/),'A','1/cm3 ','Aerosol mode 14 number concentration') + call addfld ('AIRMASSL',(/'lev'/),'A','kg/m2 ','Layer airmass') + call addfld ('BETOTVIS',(/'lev'/),'A','1/km' ,'Aerosol 3d extinction at 0.442-0.625') ! CAM4-Oslo: 0.35-0.64um + call addfld ('BATOTVIS',(/'lev'/),'A','1/km' ,'Aerosol 3d absorption at 0.442-0.625') ! CAM4-Oslo: 0.35-0.64um + call addfld ('BATSW13 ',(/'lev'/),'A','1/km' ,'Aerosol 3d SW absorption at 3.077-3.846um') + call addfld ('BATLW01 ',(/'lev'/),'A','1/km' ,'Aerosol 3d LW absorption depth at 3.077-3.846um') + + do indx=1,nbmodes + if (indx /= 3) then + write(varName,'(A,I2.2)') "Camrel", indx + call addfld(varName, (/'lev'/),'A','unitless', 'relative added mass for mode'//trim(varName)) + write(varName,'(A,I2.2)') "Cxsrel", indx + call addfld(varName, horiz_only, 'A', 'unitless', 'relative exessive added mass column for mode'//trim(varname)) + end if + enddo + + call add_default ('AKCXS ', 1, ' ') + call add_default ('PMTOT ', 1, ' ') + call add_default ('PM25 ', 1, ' ') + call add_default ('PM2P5 ', 1, ' ') + call add_default ('MMRPM2P5', 1, ' ') + call add_default ('MMRPM1 ', 1, ' ') + call add_default ('GRIDAREA', 1, ' ') + call add_default ('DAERH2O ', 1, ' ') + call add_default ('MMR_AH2O', 1, ' ') + call add_default ('ECDRYAER', 1, ' ') + call add_default ('ABSDRYAE', 1, ' ') + call add_default ('ECDRY440', 1, ' ') + call add_default ('ABSDR440', 1, ' ') + call add_default ('ECDRY870', 1, ' ') + call add_default ('ABSDR870', 1, ' ') + call add_default ('ASYMMDRY', 1, ' ') + call add_default ('ECDRYLT1', 1, ' ') + call add_default ('ABSDRYBC', 1, ' ') + call add_default ('ABSDRYOC', 1, ' ') + call add_default ('ABSDRYSU', 1, ' ') + call add_default ('ABSDRYSS', 1, ' ') + call add_default ('ABSDRYDU', 1, ' ') + call add_default ('OD550DRY', 1, ' ') + call add_default ('AB550DRY', 1, ' ') + call add_default ('DERLT05 ', 1, ' ') + call add_default ('DERGT05 ', 1, ' ') + call add_default ('DER ', 1, ' ') + call add_default ('DOD440 ', 1, ' ') + call add_default ('ABS440 ', 1, ' ') + call add_default ('DOD500 ', 1, ' ') + call add_default ('ABS500 ', 1, ' ') + call add_default ('DOD550 ', 1, ' ') + call add_default ('ABS550 ', 1, ' ') + call add_default ('ABS550AL', 1, ' ') + call add_default ('DOD670 ', 1, ' ') + call add_default ('ABS670 ', 1, ' ') + call add_default ('DOD870 ', 1, ' ') + call add_default ('ABS870 ', 1, ' ') + call add_default ('DLOAD_MI', 1, ' ') + call add_default ('DLOAD_SS', 1, ' ') + call add_default ('DLOAD_S4', 1, ' ') + call add_default ('DLOAD_OC', 1, ' ') + call add_default ('DLOAD_BC', 1, ' ') + call add_default ('LOADBCAC', 1, ' ') + call add_default ('LOADBC0 ', 1, ' ') + call add_default ('LOADBC2 ', 1, ' ') + call add_default ('LOADBC4 ', 1, ' ') + call add_default ('LOADBC12', 1, ' ') + call add_default ('LOADBC14', 1, ' ') + call add_default ('LOADOCAC', 1, ' ') + call add_default ('LOADOC4 ', 1, ' ') + call add_default ('LOADOC14', 1, ' ') + ! + call add_default ('EC550AER', 1, ' ') + call add_default ('ABS550_A', 1, ' ') + call add_default ('BS550AER', 1, ' ') + ! + call add_default ('EC550SO4', 1, ' ') + call add_default ('EC550BC ', 1, ' ') + call add_default ('EC550POM', 1, ' ') + call add_default ('EC550SS ', 1, ' ') + call add_default ('EC550DU ', 1, ' ') + ! + call add_default ('CDOD440 ', 1, ' ') + call add_default ('CDOD550 ', 1, ' ') + call add_default ('CABS550 ', 1, ' ') + call add_default ('CABS550A', 1, ' ') + call add_default ('CDOD870 ', 1, ' ') + call add_default ('A550_DU ', 1, ' ') + call add_default ('A550_SS ', 1, ' ') + call add_default ('A550_SO4', 1, ' ') + call add_default ('A550_POM', 1, ' ') + call add_default ('A550_BC ', 1, ' ') + call add_default ('D440_DU ', 1, ' ') + call add_default ('D440_SS ', 1, ' ') + call add_default ('D440_SO4', 1, ' ') + call add_default ('D440_POM', 1, ' ') + call add_default ('D440_BC ', 1, ' ') + call add_default ('D500_DU ', 1, ' ') + call add_default ('D500_SS ', 1, ' ') + call add_default ('D500_SO4', 1, ' ') + call add_default ('D500_POM', 1, ' ') + call add_default ('D500_BC ', 1, ' ') + call add_default ('D550_DU ', 1, ' ') + call add_default ('D550_SS ', 1, ' ') + call add_default ('D550_SO4', 1, ' ') + call add_default ('D550_POM', 1, ' ') + call add_default ('D550_BC ', 1, ' ') + call add_default ('D670_DU ', 1, ' ') + call add_default ('D670_SS ', 1, ' ') + call add_default ('D670_SO4', 1, ' ') + call add_default ('D670_POM', 1, ' ') + call add_default ('D670_BC ', 1, ' ') + call add_default ('D870_DU ', 1, ' ') + call add_default ('D870_SS ', 1, ' ') + call add_default ('D870_SO4', 1, ' ') + call add_default ('D870_POM', 1, ' ') + call add_default ('D870_BC ', 1, ' ') + call add_default ('DLT_DUST', 1, ' ') + call add_default ('DLT_SS ', 1, ' ') + call add_default ('DLT_SO4 ', 1, ' ') + call add_default ('DLT_POM ', 1, ' ') + call add_default ('DLT_BC ', 1, ' ') + call add_default ('DGT_DUST', 1, ' ') + call add_default ('DGT_SS ', 1, ' ') + call add_default ('DGT_SO4 ', 1, ' ') + call add_default ('DGT_POM ', 1, ' ') + call add_default ('DGT_BC ', 1, ' ') + call add_default ('NNAT_0 ', 1, ' ') + call add_default ('NNAT_1 ', 1, ' ') + call add_default ('NNAT_2 ', 1, ' ') + call add_default ('NNAT_4 ', 1, ' ') + call add_default ('NNAT_5 ', 1, ' ') + call add_default ('NNAT_6 ', 1, ' ') + call add_default ('NNAT_7 ', 1, ' ') + call add_default ('NNAT_8 ', 1, ' ') + call add_default ('NNAT_9 ', 1, ' ') + call add_default ('NNAT_10 ', 1, ' ') + call add_default ('NNAT_12 ', 1, ' ') + call add_default ('NNAT_14 ', 1, ' ') + call add_default ('AIRMASSL', 1, ' ') + call add_default ('AIRMASS ', 1, ' ') + call add_default ('BETOTVIS', 1, ' ') + call add_default ('BATOTVIS', 1, ' ') + call add_default ('BATSW13 ', 1, ' ') + call add_default ('BATLW01 ', 1, ' ') + + do indx=1,nbmodes + if (indx /= 3) then + write(varName,'(A,I2.2)') "Camrel", indx + call add_default(varName, 1, ' ') + write(varName,'(A,I2.2)') "Cxsrel", indx + call add_default(varName, 1, ' ') + end if + enddo + end if ! end of if use_aerocom + + end subroutine oslo_aero_diagnostics_init + +end module oslo_aero_diagnostics diff --git a/src/oslo_aero_optical_params.F90 b/src/oslo_aero_optical_params.F90 index e548f19..bfe9ca1 100644 --- a/src/oslo_aero_optical_params.F90 +++ b/src/oslo_aero_optical_params.F90 @@ -18,6 +18,7 @@ module oslo_aero_optical_params use oslo_aero_conc, only: calculateBulkProperties, partitionMass use oslo_aero_sw_tables, only: interpol0, interpol1, interpol2to3, interpol4, interpol5to10 use oslo_aero_aerocom, only: aerocom1, aerocom2 + use oslo_aero_control, only: use_aerocom use perf_mod, only: t_startf, t_stopf implicit none @@ -153,19 +154,21 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, end do ! interpol-calculations only when daylight or not: -#ifdef AEROCOM ! always calculate optics (also at (polar) night) - do icol=1,ncol - daylight(icol) = .true. - end do -#else ! calculate optics only in daytime - do icol=1,ncol - if (coszrs(icol) > 0.0_r8) then + if (use_aerocom) then + ! always calculate optics (also at (polar) night) + do icol=1,ncol daylight(icol) = .true. - else - daylight(icol) = .false. - endif - end do -#endif + end do + else + ! calculate optics only in daytime + do icol=1,ncol + if (coszrs(icol) > 0.0_r8) then + daylight(icol) = .true. + else + daylight(icol) = .false. + endif + end do + end if ! Set SO4, BC and OC concentrations: @@ -246,12 +249,12 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, focm, fcm, xfac, ifac1, fbcm, xfbc, ifbc1, faqm, xfaq, ifaq1) call t_stopf('oslo_aero_input_interpol') -#ifdef AEROCOM - call aerocom1(lchnk, ncol, Cam, Nnatk, deltah_km, & - xct, ict1, xfac, ifac1, xfbc, ifbc1, xfaq, ifaq1, & - xfbcbg, ifbcbg1, xfbcbgn, ifbcbgn1, & - xfombg, ifombg1, Ctotdry) -#endif + if (use_aerocom) then + call aerocom1(lchnk, ncol, Cam, Nnatk, deltah_km, & + xct, ict1, xfac, ifac1, xfbc, ifbc1, xfaq, ifaq1, & + xfbcbg, ifbcbg1, xfbcbgn, ifbcbgn1, & + xfombg, ifombg1, Ctotdry) + end if call t_startf('oslo_aero_interpol') ! (Wet) Optical properties for each of the aerosol modes: @@ -475,18 +478,18 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, end do end do -#ifdef AEROCOM - do icol=1,ncol - do ilev=1,pver - batotsw13(icol,ilev)=betot(icol,ilev,13)*(1.0_r8-ssatot(icol,ilev,13)) - batotlw01(icol,ilev)=batotlw(icol,ilev,1) + if (use_aerocom) then + do icol=1,ncol + do ilev=1,pver + batotsw13(icol,ilev)=betot(icol,ilev,13)*(1.0_r8-ssatot(icol,ilev,13)) + batotlw01(icol,ilev)=batotlw(icol,ilev,1) + end do end do - end do - ! These two fields should be close to equal, both representing absorption - ! in the 3.077-3.846 um wavelenght band (i.e., a check of LUT for LW vs. SW). - call outfld('BATSW13 ',batotsw13,pcols,lchnk) - call outfld('BATLW01 ',batotlw01,pcols,lchnk) -#endif + ! These two fields should be close to equal, both representing absorption + ! in the 3.077-3.846 um wavelenght band (i.e., a check of LUT for LW vs. SW). + call outfld('BATSW13 ',batotsw13,pcols,lchnk) + call outfld('BATLW01 ',batotlw01,pcols,lchnk) + end if ! APPROXIMATE aerosol extinction and absorption at 550nm (0.442-0.625 um) ! (in the visible wavelength band) @@ -561,37 +564,37 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, call outfld('ABSVVOLC',absvisvolc ,pcols,lchnk) call outfld('BVISVOLC',bevisvolc ,pcols,lchnk) -#ifdef AEROCOM - ! Extinction and absorption for 0.55 um for the total aerosol, and AODs - call outfld('BETOTVIS',betotvis,pcols,lchnk) - call outfld('BATOTVIS',batotvis,pcols,lchnk) - call outfld('AIRMASSL',airmassl,pcols,lchnk) - call outfld('AIRMASS ',airmass ,pcols,lchnk) + if (use_aerocom) then + ! Extinction and absorption for 0.55 um for the total aerosol, and AODs + call outfld('BETOTVIS',betotvis,pcols,lchnk) + call outfld('BATOTVIS',batotvis,pcols,lchnk) + call outfld('AIRMASSL',airmassl,pcols,lchnk) + call outfld('AIRMASS ',airmass ,pcols,lchnk) - ! Mass concentration (ug/m3) and mmr (kg/kg) of aerosol condensed water - ! Condensed water mmr (kg/kg) - do ilev=1,pver - do icol=1,ncol - Cwater(icol,ilev) = Ctot(icol,ilev) - Ctotdry(icol,ilev) - mmr_aerh2o(icol,ilev)=1.e-9_r8*Cwater(icol,ilev)/rhoda(icol,ilev) - end do - enddo - call outfld('MMR_AH2O',mmr_aerh2o, pcols, lchnk) + ! Mass concentration (ug/m3) and mmr (kg/kg) of aerosol condensed water + ! Condensed water mmr (kg/kg) + do ilev=1,pver + do icol=1,ncol + Cwater(icol,ilev) = Ctot(icol,ilev) - Ctotdry(icol,ilev) + mmr_aerh2o(icol,ilev)=1.e-9_r8*Cwater(icol,ilev)/rhoda(icol,ilev) + end do + enddo + call outfld('MMR_AH2O',mmr_aerh2o, pcols, lchnk) - ! Condensed water loading (mg_m2) - daerh2o(:) = 0.0_r8 - do ilev=1,pver - do icol=1,ncol - daerh2o(icol) = daerh2o(icol) + Cwater(icol,ilev)*deltah_km(icol,ilev) + ! Condensed water loading (mg_m2) + daerh2o(:) = 0.0_r8 + do ilev=1,pver + do icol=1,ncol + daerh2o(icol) = daerh2o(icol) + Cwater(icol,ilev)*deltah_km(icol,ilev) + end do end do - end do - call outfld('DAERH2O ',daerh2o ,pcols,lchnk) + call outfld('DAERH2O ',daerh2o ,pcols,lchnk) - ! Aerocom second phase - call aerocom2(lchnk, ncol, Nnatk, pint, deltah_km, faitbc, f_soana, fnbc, rhoda, v_soana, & - xct, ict1, xfac, ifac1, xfbc, ifbc1, xfaq, ifaq1, xfbcbg, ifbcbg1, xfbcbgn, ifbcbgn1, & - xfombg, ifombg1, xrh, irh1) -#endif + ! Aerocom second phase + call aerocom2(lchnk, ncol, Nnatk, pint, deltah_km, faitbc, f_soana, fnbc, rhoda, v_soana, & + xct, ict1, xfac, ifac1, xfbc, ifbc1, xfaq, ifaq1, xfbcbg, ifbcbg1, xfbcbgn, ifbcbgn1, & + xfombg, ifombg1, xrh, irh1) + end if call t_stopf('oslo_aero_optical_params') diff --git a/src_cam/physpkg.F90 b/src_cam/physpkg.F90 new file mode 100644 index 0000000..9baa67b --- /dev/null +++ b/src_cam/physpkg.F90 @@ -0,0 +1,3025 @@ +module physpkg + !----------------------------------------------------------------------- + ! Purpose: + ! + ! Provides the interface to CAM physics package + ! + ! Module contains reordered physics to accomodate CLUBB + ! Modified after original physpkg module, Dec 2021, A. Herrington + !----------------------------------------------------------------------- + + use shr_kind_mod, only: r8 => shr_kind_r8 + use spmd_utils, only: masterproc + use physconst, only: latvap, latice + use physics_types, only: physics_state, physics_tend, physics_state_set_grid, & + physics_ptend, physics_tend_init, physics_update, & + physics_type_alloc, physics_ptend_dealloc,& + physics_state_alloc, physics_state_dealloc, physics_tend_alloc, physics_tend_dealloc + use phys_grid, only: get_ncols_p + use phys_gmean, only: gmean_mass + use ppgrid, only: begchunk, endchunk, pcols, pver, pverp, psubcols + use constituents, only: pcnst, cnst_name, cnst_get_ind + use camsrfexch, only: cam_out_t, cam_in_t + + use phys_control, only: use_hemco ! Use Harmonized Emissions Component (HEMCO) + + use cam_control_mod, only: ideal_phys, adiabatic + use phys_control, only: phys_do_flux_avg, phys_getopts, waccmx_is + use scamMod, only: single_column, scm_crm_mode + use flux_avg, only: flux_avg_init + use perf_mod + use cam_logfile, only: iulog + use camsrfexch, only: cam_export + + use modal_aero_calcsize, only: modal_aero_calcsize_init, modal_aero_calcsize_diag, modal_aero_calcsize_reg + use modal_aero_wateruptake, only: modal_aero_wateruptake_init, modal_aero_wateruptake_dr, modal_aero_wateruptake_reg + + implicit none + private + save + + ! Public methods + public phys_register ! was initindx - register physics methods + public phys_init ! Public initialization method + public phys_run1 ! First phase of the public run method + public phys_run2 ! Second phase of the public run method + public phys_final ! Public finalization method + + ! Private module data + + ! Physics package options + character(len=16) :: shallow_scheme + character(len=16) :: macrop_scheme + character(len=16) :: microp_scheme + character(len=16) :: subcol_scheme + character(len=32) :: cam_take_snapshot_before ! Physics routine to take a snapshot "before" + character(len=32) :: cam_take_snapshot_after ! Physics routine to take a snapshot "after" + integer :: cld_macmic_num_steps ! Number of macro/micro substeps + integer :: cam_snapshot_before_num ! tape number for before snapshots + integer :: cam_snapshot_after_num ! tape number for after snapshots + logical :: do_clubb_sgs + logical :: use_subcol_microp ! if true, use subcolumns in microphysics + logical :: state_debug_checks ! Debug physics_state. + logical :: clim_modal_aero ! climate controled by prognostic or prescribed modal aerosols + logical :: prog_modal_aero ! Prognostic modal aerosols present + + ! Physics buffer index + integer :: teout_idx = 0 + + integer :: landm_idx = 0 + integer :: sgh_idx = 0 + integer :: sgh30_idx = 0 + + integer :: qini_idx = 0 + integer :: cldliqini_idx = 0 + integer :: cldiceini_idx = 0 + integer :: totliqini_idx = 0 + integer :: toticeini_idx = 0 + + integer :: prec_str_idx = 0 + integer :: snow_str_idx = 0 + integer :: prec_sed_idx = 0 + integer :: snow_sed_idx = 0 + integer :: prec_pcw_idx = 0 + integer :: snow_pcw_idx = 0 + integer :: prec_dp_idx = 0 + integer :: snow_dp_idx = 0 + integer :: prec_sh_idx = 0 + integer :: snow_sh_idx = 0 + integer :: dlfzm_idx = 0 ! detrained convective cloud water mixing ratio. + integer :: ducore_idx = 0 ! ducore index in physics buffer + integer :: dvcore_idx = 0 ! dvcore index in physics buffer + integer :: dtcore_idx = 0 ! dtcore index in physics buffer + integer :: dqcore_idx = 0 ! dqcore index in physics buffer + integer :: cmfmczm_idx = 0 ! Zhang-McFarlane convective mass fluxes + integer :: rliqbc_idx = 0 ! tphysbc reserve liquid +!======================================================================= +contains +!======================================================================= + + subroutine phys_register + !----------------------------------------------------------------------- + ! + ! Purpose: Register constituents and physics buffer fields. + ! + ! Author: CSM Contact: M. Vertenstein, Aug. 1997 + ! B.A. Boville, Oct 2001 + ! A. Gettelman, Nov 2010 - put micro/macro physics into separate routines + ! + !----------------------------------------------------------------------- + use cam_abortutils, only: endrun + use physics_buffer, only: pbuf_init_time, pbuf_cam_snapshot_register + use physics_buffer, only: pbuf_add_field, dtype_r8, pbuf_register_subcol + use constituents, only: cnst_add, cnst_chk_dim + + use cam_control_mod, only: moist_physics + use chemistry, only: chem_register + use mo_lightning, only: lightning_register + use cloud_fraction, only: cldfrc_register + use microp_driver, only: microp_driver_register + use microp_aero, only: microp_aero_register + ! OSLO_AERO begin + use oslo_aero_microp, only: oslo_aero_microp_register + ! OSLO_AERO end + use macrop_driver, only: macrop_driver_register + use clubb_intr, only: clubb_register_cam + use conv_water, only: conv_water_register + use physconst, only: mwh2o, cpwv + use tracers, only: tracers_register + use check_energy, only: check_energy_register + use carma_intr, only: carma_register + use ghg_data, only: ghg_data_register + use vertical_diffusion, only: vd_register + use convect_deep, only: convect_deep_register + use convect_diagnostics,only: convect_diagnostics_register + use radiation, only: radiation_register + use co2_cycle, only: co2_register + use flux_avg, only: flux_avg_register + use iondrag, only: iondrag_register + use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_reg + use prescribed_ozone, only: prescribed_ozone_register + use prescribed_volcaero,only: prescribed_volcaero_register + use prescribed_strataero,only: prescribed_strataero_register + use prescribed_aero, only: prescribed_aero_register + use prescribed_ghg, only: prescribed_ghg_register + use aoa_tracers, only: aoa_tracers_register + use aircraft_emit, only: aircraft_emit_register + use cam_diagnostics, only: diag_register + use cloud_diagnostics, only: cloud_diagnostics_register + use cospsimulator_intr, only: cospsimulator_intr_register + use rad_constituents, only: rad_cnst_get_info ! Added to query if it is a modal aero sim or not + use radheat, only: radheat_register + use subcol, only: subcol_register + use subcol_utils, only: is_subcol_on, subcol_get_scheme + use dyn_comp, only: dyn_register + use offline_driver, only: offline_driver_reg + use hemco_interface, only: HCOI_Chunk_Init + + !---------------------------Local variables----------------------------- + ! + integer :: m ! loop index + integer :: mm ! constituent index + integer :: nmodes + !----------------------------------------------------------------------- + + ! Get physics options + call phys_getopts(shallow_scheme_out = shallow_scheme, & + macrop_scheme_out = macrop_scheme, & + microp_scheme_out = microp_scheme, & + cld_macmic_num_steps_out = cld_macmic_num_steps, & + do_clubb_sgs_out = do_clubb_sgs, & + use_subcol_microp_out = use_subcol_microp, & + state_debug_checks_out = state_debug_checks, & + cam_take_snapshot_before_out= cam_take_snapshot_before, & + cam_take_snapshot_after_out = cam_take_snapshot_after, & + cam_snapshot_before_num_out = cam_snapshot_before_num, & + cam_snapshot_after_num_out = cam_snapshot_after_num) + + subcol_scheme = subcol_get_scheme() + + ! Initialize dyn_time_lvls + call pbuf_init_time() + + ! Register the subcol scheme + call subcol_register() + + ! Register water vapor. + ! ***** N.B. ***** This must be the first call to cnst_add so that + ! water vapor is constituent 1. + if (moist_physics) then + call cnst_add('Q', mwh2o, cpwv, 1.E-12_r8, mm, & + longname='Specific humidity', readiv=.true., is_convtran1=.true.) + else + call cnst_add('Q', mwh2o, cpwv, 0.0_r8, mm, & + longname='Specific humidity', readiv=.false., is_convtran1=.true.) + end if + + ! Topography file fields. + call pbuf_add_field('LANDM', 'global', dtype_r8, (/pcols/), landm_idx) + call pbuf_add_field('SGH', 'global', dtype_r8, (/pcols/), sgh_idx) + call pbuf_add_field('SGH30', 'global', dtype_r8, (/pcols/), sgh30_idx) + + ! Fields for physics package diagnostics + call pbuf_add_field('QINI', 'physpkg', dtype_r8, (/pcols,pver/), qini_idx) + call pbuf_add_field('CLDLIQINI', 'physpkg', dtype_r8, (/pcols,pver/), cldliqini_idx) + call pbuf_add_field('CLDICEINI', 'physpkg', dtype_r8, (/pcols,pver/), cldiceini_idx) + call pbuf_add_field('TOTLIQINI', 'physpkg', dtype_r8, (/pcols,pver/), totliqini_idx) + call pbuf_add_field('TOTICEINI', 'physpkg', dtype_r8, (/pcols,pver/), toticeini_idx) + + ! check energy package + call check_energy_register + + ! If using a simple physics option (e.g., held_suarez, adiabatic), + ! the normal CAM physics parameterizations are not called. + if (moist_physics) then + + ! register fluxes for saving across time + if (phys_do_flux_avg()) call flux_avg_register() + + call cldfrc_register() + + ! cloud water + if (.not. do_clubb_sgs) call macrop_driver_register() + ! OSLO_AERO begin + call oslo_aero_microp_register() + ! OSLO_AERO end + call microp_driver_register() + + ! Register CLUBB_SGS here + if (do_clubb_sgs) call clubb_register_cam() + + call pbuf_add_field('PREC_STR', 'global',dtype_r8,(/pcols/),prec_str_idx) + call pbuf_add_field('SNOW_STR', 'global',dtype_r8,(/pcols/),snow_str_idx) + call pbuf_add_field('PREC_PCW', 'global',dtype_r8,(/pcols/),prec_pcw_idx) + call pbuf_add_field('SNOW_PCW', 'global',dtype_r8,(/pcols/),snow_pcw_idx) + call pbuf_add_field('PREC_SED', 'global',dtype_r8,(/pcols/),prec_sed_idx) + call pbuf_add_field('SNOW_SED', 'global',dtype_r8,(/pcols/),snow_sed_idx) + + if (is_subcol_on()) then + call pbuf_register_subcol('PREC_STR', 'phys_register', prec_str_idx) + call pbuf_register_subcol('SNOW_STR', 'phys_register', snow_str_idx) + call pbuf_register_subcol('PREC_PCW', 'phys_register', prec_pcw_idx) + call pbuf_register_subcol('SNOW_PCW', 'phys_register', snow_pcw_idx) + call pbuf_register_subcol('PREC_SED', 'phys_register', prec_sed_idx) + call pbuf_register_subcol('SNOW_SED', 'phys_register', snow_sed_idx) + end if + + ! Reserve liquid at end of tphysbc + call pbuf_add_field('RLIQBC','physpkg',dtype_r8,(/pcols/),rliqbc_idx) + + ! Who should add FRACIS? + ! -- It does not seem that aero_intr should add it since FRACIS is used in convection + ! even if there are no prognostic aerosols ... so do it here for now + call pbuf_add_field('FRACIS','physpkg',dtype_r8,(/pcols,pver,pcnst/),m) + + call conv_water_register() + + ! Determine whether its a 'modal' aerosol simulation or not + ! OSLO_AERO begin + clim_modal_aero = .false. + ! OSLO_AERO end + + ! register chemical constituents including aerosols ... + call chem_register() + + ! add prognostic lightning flash freq pbuf fld + call lightning_register() + + ! co2 constituents + call co2_register() + + ! register other constituents + call prescribed_volcaero_register() + call prescribed_strataero_register() + call prescribed_ozone_register() + call prescribed_aero_register() + call prescribed_ghg_register() + + ! register various data model gasses with pbuf + call ghg_data_register() + + ! carma microphysics + ! + call carma_register() + + ! Register iondrag variables with pbuf + call iondrag_register() + + ! Register ionosphere variables with pbuf if mode set to ionosphere + if( waccmx_is('ionosphere') ) then + call waccmx_phys_ion_elec_temp_reg() + endif + + call aircraft_emit_register() + + ! deep convection + call convect_deep_register + + ! convection diagnostics + call convect_diagnostics_register + + ! radiation + call radiation_register + call cloud_diagnostics_register + call radheat_register + + ! COSP + call cospsimulator_intr_register + + ! vertical diffusion + call vd_register() + else + ! held_suarez/adiabatic physics option should be in simple_physics + call endrun('phys_register: moist_physics configuration error') + end if + + ! Register diagnostics PBUF + call diag_register() + + ! Register age of air tracers + call aoa_tracers_register() + + ! Register test tracers + call tracers_register() + + call dyn_register() + + ! All tracers registered, check that the dimensions are correct + call cnst_chk_dim() + + ! ***NOTE*** No registering constituents after the call to cnst_chk_dim. + + call offline_driver_reg() + + if (use_hemco) then + ! initialize harmonized emissions component (HEMCO) + call HCOI_Chunk_Init() + endif + + ! This needs to be last as it requires all pbuf fields to be added + if (cam_snapshot_before_num > 0 .or. cam_snapshot_after_num > 0) then + call pbuf_cam_snapshot_register() + end if + + end subroutine phys_register + + + + !======================================================================= + + subroutine phys_inidat( cam_out, pbuf2d ) + use cam_abortutils, only: endrun + + use physics_buffer, only: pbuf_get_index, physics_buffer_desc, pbuf_set_field, dyn_time_lvls + + + use cam_initfiles, only: initial_file_get_id, topo_file_get_id + use cam_grid_support, only: cam_grid_check, cam_grid_id + use cam_grid_support, only: cam_grid_get_dim_names + use pio, only: file_desc_t + use ncdio_atm, only: infld + use dycore, only: dycore_is + use polar_avg, only: polar_average + use short_lived_species, only: initialize_short_lived_species + use cam_control_mod, only: aqua_planet + use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_inidat + + type(cam_out_t), intent(inout) :: cam_out(begchunk:endchunk) + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + integer :: lchnk, m, n, i, k, ncol + type(file_desc_t), pointer :: fh_ini, fh_topo + character(len=8) :: fieldname + real(r8), pointer :: tptr(:,:), tptr_2(:,:), tptr3d(:,:,:), tptr3d_2(:,:,:) + real(r8), pointer :: qpert(:,:) + + character(len=11) :: subname='phys_inidat' ! subroutine name + integer :: tpert_idx, qpert_idx, pblh_idx + + logical :: found=.false., found2=.false. + integer :: ierr + character(len=8) :: dim1name, dim2name + integer :: ixcldice, ixcldliq + integer :: grid_id ! grid ID for data mapping + + nullify(tptr,tptr_2,tptr3d,tptr3d_2) + + fh_ini => initial_file_get_id() + fh_topo => topo_file_get_id() + + ! dynamics variables are handled in dyn_init - here we read variables needed for physics + ! but not dynamics + + grid_id = cam_grid_id('physgrid') + if (.not. cam_grid_check(grid_id)) then + call endrun(trim(subname)//': Internal error, no "physgrid" grid') + end if + call cam_grid_get_dim_names(grid_id, dim1name, dim2name) + + allocate(tptr(1:pcols,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr(1:pcols,begchunk:endchunk)') + end if + + if (associated(fh_topo) .and. .not. aqua_planet) then + call infld('SGH', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr, found, gridname='physgrid') + if(.not. found) call endrun('ERROR: SGH not found on topo file') + + call pbuf_set_field(pbuf2d, sgh_idx, tptr) + + allocate(tptr_2(1:pcols,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr_2(1:pcols,begchunk:endchunk)') + end if + call infld('SGH30', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr_2, found, gridname='physgrid') + if(found) then + call pbuf_set_field(pbuf2d, sgh30_idx, tptr_2) + else + if (masterproc) write(iulog,*) 'Warning: Error reading SGH30 from topo file.' + if (masterproc) write(iulog,*) 'The field SGH30 will be filled using data from SGH.' + call pbuf_set_field(pbuf2d, sgh30_idx, tptr) + end if + + deallocate(tptr_2) + + call infld('LANDM_COSLAT', fh_topo, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr, found, gridname='physgrid') + + if(.not.found) call endrun(' ERROR: LANDM_COSLAT not found on topo dataset.') + + call pbuf_set_field(pbuf2d, landm_idx, tptr) + + else + call pbuf_set_field(pbuf2d, sgh_idx, 0._r8) + call pbuf_set_field(pbuf2d, sgh30_idx, 0._r8) + call pbuf_set_field(pbuf2d, landm_idx, 0._r8) + end if + + call infld('PBLH', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr(:,:), found, gridname='physgrid') + if(.not. found) then + tptr(:,:) = 0._r8 + if (masterproc) write(iulog,*) 'PBLH initialized to 0.' + end if + pblh_idx = pbuf_get_index('pblh') + + call pbuf_set_field(pbuf2d, pblh_idx, tptr) + + call infld('TPERT', fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr(:,:), found, gridname='physgrid') + if(.not. found) then + tptr(:,:) = 0._r8 + if (masterproc) write(iulog,*) 'TPERT initialized to 0.' + end if + tpert_idx = pbuf_get_index( 'tpert') + call pbuf_set_field(pbuf2d, tpert_idx, tptr) + + fieldname='QPERT' + qpert_idx = pbuf_get_index( 'qpert',ierr) + if (qpert_idx > 0) then + call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr, found, gridname='physgrid') + if(.not. found) then + tptr=0_r8 + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.' + end if + + allocate(tptr3d_2(pcols,pcnst,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr3d_2(pcols,pcnst,begchunk:endchunk)') + end if + tptr3d_2 = 0_r8 + tptr3d_2(:,1,:) = tptr(:,:) + + call pbuf_set_field(pbuf2d, qpert_idx, tptr3d_2) + deallocate(tptr3d_2) + end if + + fieldname='CUSH' + m = pbuf_get_index('cush', ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, dim2name, 1, pcols, begchunk, endchunk, & + tptr, found, gridname='physgrid') + if(.not.found) then + if(masterproc) write(iulog,*) trim(fieldname), ' initialized to 1000.' + tptr=1000._r8 + end if + do n=1,dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr, start=(/1,n/), kount=(/pcols,1/)) + end do + deallocate(tptr) + end if + + ! + ! 3-D fields + ! + + allocate(tptr3d(pcols,pver,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr3d(pcols,pver,begchunk:endchunk)') + end if + + fieldname='CLOUD' + m = pbuf_get_index('CLD') + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(found) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + else + call pbuf_set_field(pbuf2d, m, 0._r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.' + end if + + fieldname='QCWAT' + m = pbuf_get_index(fieldname,ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(.not. found) then + call infld('Q',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if (found) then + if (masterproc) write(iulog,*) trim(fieldname), ' initialized with Q' + if(dycore_is('LR')) call polar_average(pver, tptr3d) + else + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to huge()' + tptr3d = huge(1.0_r8) + end if + end if + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + end if + + fieldname = 'ICCWAT' + m = pbuf_get_index(fieldname, ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(found) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + else + call cnst_get_ind('CLDICE', ixcldice) + call infld('CLDICE',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(found) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + else + call pbuf_set_field(pbuf2d, m, 0._r8) + end if + if (masterproc) then + if (found) then + write(iulog,*) trim(fieldname), ' initialized with CLDICE' + else + write(iulog,*) trim(fieldname), ' initialized to 0.0' + end if + end if + end if + end if + + fieldname = 'LCWAT' + m = pbuf_get_index(fieldname,ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(found) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + else + allocate(tptr3d_2(pcols,pver,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr3d_2(pcols,pver,begchunk:endchunk)') + end if + call cnst_get_ind('CLDICE', ixcldice) + call cnst_get_ind('CLDLIQ', ixcldliq) + call infld('CLDICE',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + call infld('CLDLIQ',fh_ini,dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d_2, found2, gridname='physgrid') + if(found .and. found2) then + do lchnk = begchunk, endchunk + ncol = get_ncols_p(lchnk) + tptr3d(:ncol,:,lchnk)=tptr3d(:ncol,:,lchnk)+tptr3d_2(:ncol,:,lchnk) + end do + if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDICE + CLDLIQ' + else if (found) then ! Data already loaded in tptr3d + if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDICE only' + else if (found2) then + tptr3d(:,:,:)=tptr3d_2(:,:,:) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized with CLDLIQ only' + end if + + if (found .or. found2) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + if(dycore_is('LR')) call polar_average(pver, tptr3d) + else + call pbuf_set_field(pbuf2d, m, 0._r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.0' + end if + deallocate(tptr3d_2) + end if + end if + + fieldname = 'TCWAT' + m = pbuf_get_index(fieldname,ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(.not.found) then + call infld('T', fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if (found) then + if(dycore_is('LR')) call polar_average(pver, tptr3d) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized with T' + else + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to huge()' + tptr3d = huge(1._r8) + end if + end if + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + end if + + fieldname = 'CONCLD' + m = pbuf_get_index('CONCLD',ierr) + if (m > 0) then + call infld(fieldname, fh_ini, dim1name, 'lev', dim2name, 1, pcols, 1, pver, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if(found) then + do n = 1, dyn_time_lvls + call pbuf_set_field(pbuf2d, m, tptr3d, (/1,1,n/),(/pcols,pver,1/)) + end do + else + call pbuf_set_field(pbuf2d, m, 0._r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.' + end if + end if + + deallocate(tptr3d) + allocate(tptr3d(pcols,pverp,begchunk:endchunk), stat=ierr) + if (ierr /= 0) then + call endrun(subname//': Failed to allocate tptr3d(pcols,pver,begchunk:endchunk)') + end if + + fieldname = 'TKE' + m = pbuf_get_index( 'tke') + call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if (found) then + call pbuf_set_field(pbuf2d, m, tptr3d) + else + call pbuf_set_field(pbuf2d, m, 0.01_r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.01' + end if + + + fieldname = 'KVM' + m = pbuf_get_index('kvm') + call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if (found) then + call pbuf_set_field(pbuf2d, m, tptr3d) + else + call pbuf_set_field(pbuf2d, m, 0._r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.' + end if + + + fieldname = 'KVH' + m = pbuf_get_index('kvh') + call infld(fieldname, fh_ini, dim1name, 'ilev', dim2name, 1, pcols, 1, pverp, begchunk, endchunk, & + tptr3d, found, gridname='physgrid') + if (found) then + call pbuf_set_field(pbuf2d, m, tptr3d) + else + call pbuf_set_field(pbuf2d, m, 0._r8) + if (masterproc) write(iulog,*) trim(fieldname), ' initialized to 0.' + end if + + call initialize_short_lived_species(fh_ini, pbuf2d) + + !--------------------------------------------------------------------------------- + ! If needed, get ion and electron temperature fields from initial condition file + !--------------------------------------------------------------------------------- + + call waccmx_phys_ion_elec_temp_inidat(fh_ini,pbuf2d) + + end subroutine phys_inidat + + + subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_in, cam_out ) + + !----------------------------------------------------------------------- + ! + ! Initialization of physics package. + ! + !----------------------------------------------------------------------- + + use physics_buffer, only: physics_buffer_desc, pbuf_initialize, pbuf_get_index + use physconst, only: rair, cpair, gravit, zvir, & + karman + use cam_thermo, only: cam_thermo_init + use ref_pres, only: pref_edge, pref_mid + + use carma_intr, only: carma_init + use cam_control_mod, only: initial_run + use check_energy, only: check_energy_init + use chemistry, only: chem_init + use mo_lightning, only: lightning_init + use prescribed_ozone, only: prescribed_ozone_init + use prescribed_ghg, only: prescribed_ghg_init + use prescribed_aero, only: prescribed_aero_init + use aerodep_flx, only: aerodep_flx_init + use aircraft_emit, only: aircraft_emit_init + use prescribed_volcaero,only: prescribed_volcaero_init + use prescribed_strataero,only: prescribed_strataero_init + use cloud_fraction, only: cldfrc_init + use cldfrc2m, only: cldfrc2m_init + use co2_cycle, only: co2_init, co2_transport + use convect_deep, only: convect_deep_init + use convect_diagnostics,only: convect_diagnostics_init + use cam_diagnostics, only: diag_init + ! OSLO_AERO begin + use oslo_aero_diagnostics, only: oslo_aero_diagnostics_init + ! OSLO_AERO end + use gw_drag, only: gw_init + use radheat, only: radheat_init + use radiation, only: radiation_init + use cloud_diagnostics, only: cloud_diagnostics_init + use wv_saturation, only: wv_sat_init + use microp_driver, only: microp_driver_init + use microp_aero, only: microp_aero_init + ! OSLO_AERO begin + use oslo_aero_microp, only: oslo_aero_microp_init + ! OSLO_AERO end + use macrop_driver, only: macrop_driver_init + use conv_water, only: conv_water_init + use tracers, only: tracers_init + use aoa_tracers, only: aoa_tracers_init + use rayleigh_friction, only: rayleigh_friction_init + use pbl_utils, only: pbl_utils_init + use vertical_diffusion, only: vertical_diffusion_init + use phys_debug_util, only: phys_debug_init + use phys_debug, only: phys_debug_state_init + use rad_constituents, only: rad_cnst_init + use aer_rad_props, only: aer_rad_props_init + use subcol, only: subcol_init + use qbo, only: qbo_init + use qneg_module, only: qneg_init + use lunar_tides, only: lunar_tides_init + use iondrag, only: iondrag_init +#if ( defined OFFLINE_DYN ) + use metdata, only: metdata_phys_init +#endif + use epp_ionization, only: epp_ionization_init, epp_ionization_active + use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_init ! Initialization of ionosphere module (WACCM-X) + use waccmx_phys_intr, only: waccmx_phys_mspd_init ! Initialization of major species diffusion module (WACCM-X) + use clubb_intr, only: clubb_ini_cam + use tropopause, only: tropopause_init + use solar_data, only: solar_data_init + use dadadj_cam, only: dadadj_init + use cam_abortutils, only: endrun + use nudging, only: Nudge_Model, nudging_init + use cam_snapshot, only: cam_snapshot_init + use cam_history, only: addfld, register_vector_field, add_default + use cam_budget, only: cam_budget_init + use phys_grid_ctem, only: phys_grid_ctem_init + + use ccpp_constituent_prop_mod, only: ccpp_const_props_init + + ! Input/output arguments + type(physics_state), pointer :: phys_state(:) + type(physics_tend ), pointer :: phys_tend(:) + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + type(cam_in_t), intent(in) :: cam_in(begchunk:endchunk) + type(cam_out_t),intent(inout) :: cam_out(begchunk:endchunk) + + ! local variables + integer :: lchnk + integer :: ierr + + logical :: history_budget ! output tendencies and state variables for + ! temperature, water vapor, cloud + ! ice, cloud liquid, U, V + integer :: history_budget_histfile_num ! output history file number for budget fields + + !----------------------------------------------------------------------- + + call physics_type_alloc(phys_state, phys_tend, begchunk, endchunk, pcols) + + do lchnk = begchunk, endchunk + call physics_state_set_grid(lchnk, phys_state(lchnk)) + end do + + !------------------------------------------------------------------------------------------- + ! Initialize any variables in cam_thermo which are not temporally and/or spatially constant + !------------------------------------------------------------------------------------------- + call cam_thermo_init() + + ! Initialize debugging a physics column + call phys_debug_init() + + call pbuf_initialize(pbuf2d) + + ! Initialize subcol scheme + call subcol_init(pbuf2d) + + ! diag_init makes addfld calls for dynamics fields that are output from + ! the physics decomposition + call diag_init(pbuf2d) + ! OSLO_AERO begin + call oslo_aero_diagnostics_init() + ! OSLO_AERO end + + call check_energy_init() + + call tracers_init() + + ! age of air tracers + call aoa_tracers_init() + + teout_idx = pbuf_get_index( 'TEOUT') + + ! adiabatic or ideal physics should be only used if in simple_physics + if (adiabatic .or. ideal_phys) then + if (adiabatic) then + call endrun('phys_init: adiabatic configuration error') + else + call endrun('phys_init: ideal_phys configuration error') + end if + end if + + if (initial_run) then + call phys_inidat(cam_out, pbuf2d) + end if + + ! wv_saturation is relatively independent of everything else and + ! low level, so init it early. Must at least do this before radiation. + call wv_sat_init + + ! solar irradiance data modules + call solar_data_init() + + ! Initialize rad constituents and their properties + call rad_cnst_init() + + call radiation_init(pbuf2d) + + call aer_rad_props_init() + + ! initialize carma + call carma_init() + + ! Prognostic chemistry. + call chem_init(phys_state,pbuf2d) + + ! Lightning flash frq and NOx prod + call lightning_init( pbuf2d ) + + ! Prescribed tracers + call prescribed_ozone_init() + call prescribed_ghg_init() + call prescribed_aero_init() + call aerodep_flx_init() + call aircraft_emit_init() + call prescribed_volcaero_init() + call prescribed_strataero_init() + + ! co2 cycle + if (co2_transport()) then + call co2_init() + end if + + call gw_init() + + call rayleigh_friction_init() + + call pbl_utils_init(gravit, karman, cpair, rair, zvir) + call vertical_diffusion_init(pbuf2d) + + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then + call waccmx_phys_mspd_init () + ! Initialization of ionosphere module if mode set to ionosphere + if( waccmx_is('ionosphere') ) then + call waccmx_phys_ion_elec_temp_init(pbuf2d) + endif + endif + + call cloud_diagnostics_init() + + call radheat_init(pref_mid) + + call convect_diagnostics_init() + + call cldfrc_init() + call cldfrc2m_init() + + call convect_deep_init(pref_edge) + + if (.not. do_clubb_sgs) call macrop_driver_init(pbuf2d) + ! OSLO_AERO begin + call oslo_aero_microp_init() + ! OSLO_AERO end + call microp_driver_init(pbuf2d) + call conv_water_init + + ! initiate CLUBB within CAM + if (do_clubb_sgs) call clubb_ini_cam(pbuf2d) + + call qbo_init + + call lunar_tides_init() + + call iondrag_init(pref_mid) + ! Geomagnetic module -- after iondrag_init + if (epp_ionization_active) then + call epp_ionization_init() + endif + +#if ( defined OFFLINE_DYN ) + call metdata_phys_init() +#endif + call tropopause_init() + call dadadj_init() + + prec_dp_idx = pbuf_get_index('PREC_DP') + snow_dp_idx = pbuf_get_index('SNOW_DP') + prec_sh_idx = pbuf_get_index('PREC_SH') + snow_sh_idx = pbuf_get_index('SNOW_SH') + + dlfzm_idx = pbuf_get_index('DLFZM', ierr) + cmfmczm_idx = pbuf_get_index('CMFMC_DP', ierr) + + ! OSLO_AERO begin + prog_modal_aero = .true. + ! OSLO_AERO end + + ! Initialize Nudging Parameters + !-------------------------------- + if(Nudge_Model) call nudging_init + + if (clim_modal_aero) then + + ! If climate calculations are affected by prescribed modal aerosols, the + ! initialization routine for the dry mode radius calculation is called + ! here. For prognostic MAM the initialization is called from + ! modal_aero_initialize + if (.not. prog_modal_aero) then + call modal_aero_calcsize_init(pbuf2d) + endif + + call modal_aero_wateruptake_init(pbuf2d) + + end if + + ! Initialize CAM CCPP constituent properties array + ! for use in CCPP-ized physics schemes: + call ccpp_const_props_init() + + ! Initialize qneg3 and qneg4 + call qneg_init() + + ! Initialize phys TEM diagnostics + call phys_grid_ctem_init() + + ! Initialize the snapshot capability + call cam_snapshot_init(cam_in, cam_out, pbuf2d, begchunk) + + ! Initialize the budget capability + call cam_budget_init() + + ! addfld calls for U, V tendency budget variables that are output in + ! tphysac, tphysbc + call addfld ( 'UTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by deep convection') + call addfld ( 'VTEND_DCONV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by deep convection') + call register_vector_field ( 'UTEND_DCONV', 'VTEND_DCONV') + call addfld ( 'UTEND_SHCONV', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by shallow convection') + call addfld ( 'VTEND_SHCONV', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by shallow convection') + call register_vector_field ( 'UTEND_SHCONV', 'VTEND_SHCONV') + call addfld ( 'UTEND_MACROP', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by macrophysics') + call addfld ( 'VTEND_MACROP', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by macrophysics') + call register_vector_field ( 'UTEND_MACROP', 'VTEND_MACROP') + call addfld ( 'UTEND_VDIFF', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by vert. diffus.') + call addfld ( 'VTEND_VDIFF', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by vert. diffus.') + call register_vector_field ( 'UTEND_VDIFF', 'VTEND_VDIFF') + call addfld ( 'UTEND_RAYLEIGH', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by Rayleigh Fric.') + call addfld ( 'VTEND_RAYLEIGH', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by Rayleigh Fric.') + call register_vector_field ( 'UTEND_RAYLEIGH', 'VTEND_RAYLEIGH') + call addfld ( 'UTEND_GWDTOT', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by all GWs') + call addfld ( 'VTEND_GWDTOT', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by all GWs') + call register_vector_field ( 'UTEND_GWDTOT', 'VTEND_GWDTOT') + call addfld ( 'UTEND_QBORLX', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by QBO relaxation') + call addfld ( 'VTEND_QBORLX', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by QBO relaxation') + call register_vector_field ( 'UTEND_QBORLX', 'VTEND_QBORLX') + call addfld ( 'UTEND_LUNART', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by lunar tides') + call addfld ( 'VTEND_LUNART', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by lunar tides') + call register_vector_field ( 'UTEND_LUNART', 'VTEND_LUNART') + call addfld ( 'UTEND_IONDRG', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by ion drag') + call addfld ( 'VTEND_IONDRG', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by ion drag') + call register_vector_field ( 'UTEND_IONDRG', 'VTEND_IONDRG') + call addfld ( 'UTEND_NDG', (/ 'lev' /), 'A', 'm/s2', 'Zonal wind tendency by nudging') + call addfld ( 'VTEND_NDG', (/ 'lev' /), 'A', 'm/s2', 'Meridional wind tendency by nudging') + call register_vector_field ( 'UTEND_NDG', 'VTEND_NDG') + call addfld('UTEND_CORE', (/ 'lev' /), 'A', 'm/s2' , 'Zonal wind tendency due to dynamical core') + call addfld('VTEND_CORE', (/ 'lev' /), 'A', 'm/s2' , 'Meridional wind tendency due to dynamical core') + call register_vector_field('UTEND_CORE','VTEND_CORE') + + + call phys_getopts(history_budget_out = history_budget, & + history_budget_histfile_num_out = history_budget_histfile_num) + + if ( history_budget ) then + call add_default ( 'UTEND_DCONV' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_DCONV' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_SHCONV' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_SHCONV' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_MACROP' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_MACROP' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_VDIFF' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_VDIFF' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_RAYLEIGH' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_RAYLEIGH' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_GWDTOT' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_GWDTOT' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_QBORLX' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_QBORLX' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_LUNART' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_LUNART' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_IONDRG' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_IONDRG' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_NDG' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_NDG' , history_budget_histfile_num, ' ') + call add_default ( 'UTEND_CORE' , history_budget_histfile_num, ' ') + call add_default ( 'VTEND_CORE' , history_budget_histfile_num, ' ') + end if + + ducore_idx = pbuf_get_index('DUCORE') + dvcore_idx = pbuf_get_index('DVCORE') + dtcore_idx = pbuf_get_index('DTCORE') + dqcore_idx = pbuf_get_index('DQCORE') + + end subroutine phys_init + + ! + !----------------------------------------------------------------------- + ! + + subroutine phys_run1(phys_state, ztodt, phys_tend, pbuf2d, cam_in, cam_out) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! First part of atmospheric physics package before updating of surface models + ! + !----------------------------------------------------------------------- + use time_manager, only: get_nstep + use cam_diagnostics,only: diag_allocate, diag_physvar_ic + use check_energy, only: check_energy_gmean + use spmd_utils, only: mpicom + use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_allocate +#if (defined BFB_CAM_SCAM_IOP ) + use cam_history, only: outfld +#endif + use cam_abortutils, only: endrun +#if ( defined OFFLINE_DYN ) + use metdata, only: get_met_srf1 +#endif + ! + ! Input arguments + ! + real(r8), intent(in) :: ztodt ! physics time step unless nstep=0 + ! + ! Input/Output arguments + ! + type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state + type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend + + type(physics_buffer_desc), pointer, dimension(:,:) :: pbuf2d + type(cam_in_t), dimension(begchunk:endchunk) :: cam_in + type(cam_out_t), dimension(begchunk:endchunk) :: cam_out + !----------------------------------------------------------------------- + ! + !---------------------------Local workspace----------------------------- + ! + integer :: c ! indices + integer :: ncol ! number of columns + integer :: nstep ! current timestep number + type(physics_buffer_desc), pointer :: phys_buffer_chunk(:) + + call t_startf ('physpkg_st1') + nstep = get_nstep() + +#if ( defined OFFLINE_DYN ) + ! + ! if offline mode set SNOWH and TS for micro-phys + ! + call get_met_srf1( cam_in ) +#endif + + ! The following initialization depends on the import state (cam_in) + ! being initialized. This isn't true when cam_init is called, so need + ! to postpone this initialization to here. + if (nstep == 0 .and. phys_do_flux_avg()) call flux_avg_init(cam_in, pbuf2d) + + ! Compute total energy of input state and previous output state + call t_startf ('chk_en_gmean') + call check_energy_gmean(phys_state, pbuf2d, ztodt, nstep) + call t_stopf ('chk_en_gmean') + + call pbuf_allocate(pbuf2d, 'physpkg') + call diag_allocate() + + !----------------------------------------------------------------------- + ! Advance time information + !----------------------------------------------------------------------- + + call phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) + + call t_stopf ('physpkg_st1') + +#ifdef TRACER_CHECK + call gmean_mass ('before tphysbc DRY', phys_state) +#endif + + + !----------------------------------------------------------------------- + ! Tendency physics before flux coupler invocation + !----------------------------------------------------------------------- + ! + +#if (defined BFB_CAM_SCAM_IOP ) + do c=begchunk, endchunk + call outfld('Tg',cam_in(c)%ts,pcols ,c ) + end do +#endif + + call t_barrierf('sync_bc_physics', mpicom) + call t_startf ('bc_physics') + call t_adj_detailf(+1) + +!$OMP PARALLEL DO PRIVATE (C, phys_buffer_chunk) + do c=begchunk, endchunk + ! + ! Output physics terms to IC file + ! + phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c) + + call t_startf ('diag_physvar_ic') + call diag_physvar_ic ( c, phys_buffer_chunk, cam_out(c), cam_in(c) ) + call t_stopf ('diag_physvar_ic') + + call tphysbc (ztodt, phys_state(c), & + phys_tend(c), phys_buffer_chunk, & + cam_out(c), cam_in(c) ) + end do + + call t_adj_detailf(-1) + call t_stopf ('bc_physics') + + ! Don't call the rest in CRM mode + if(single_column.and.scm_crm_mode) return + +#ifdef TRACER_CHECK + call gmean_mass ('between DRY', phys_state) +#endif + + end subroutine phys_run1 + + ! + !----------------------------------------------------------------------- + ! + + subroutine phys_run2(phys_state, ztodt, phys_tend, pbuf2d, cam_out, & + cam_in ) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Second part of atmospheric physics package after updating of surface models + ! + !----------------------------------------------------------------------- + use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_deallocate, pbuf_update_tim_idx + use mo_lightning, only: lightning_no_prod + use cam_diagnostics, only: diag_deallocate, diag_surf + use carma_intr, only: carma_accumulate_stats + use spmd_utils, only: mpicom + use iop_forcing, only: scam_use_iop_srf +#if ( defined OFFLINE_DYN ) + use metdata, only: get_met_srf2 +#endif + use hemco_interface, only: HCOI_Chunk_Run + ! + ! Input arguments + ! + real(r8), intent(in) :: ztodt ! physics time step unless nstep=0 + ! + ! Input/Output arguments + ! + type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state + type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend + type(physics_buffer_desc),pointer, dimension(:,:) :: pbuf2d + + type(cam_out_t), intent(inout), dimension(begchunk:endchunk) :: cam_out + type(cam_in_t), intent(inout), dimension(begchunk:endchunk) :: cam_in + ! + !----------------------------------------------------------------------- + !---------------------------Local workspace----------------------------- + ! + integer :: c ! chunk index + integer :: ncol ! number of columns + type(physics_buffer_desc),pointer, dimension(:) :: phys_buffer_chunk + ! + ! If exit condition just return + ! + + if(single_column.and.scm_crm_mode) then + call diag_deallocate() + return + end if + !----------------------------------------------------------------------- + ! if using IOP values for surface fluxes overwrite here after surface components run + !----------------------------------------------------------------------- + if (single_column) call scam_use_iop_srf(cam_in) + + if(use_hemco) then + !---------------------------------------------------------- + ! run hemco (phase 2 before chemistry) + ! only phase 2 is used currently for HEMCO-CESM + !---------------------------------------------------------- + call HCOI_Chunk_Run(cam_in, phys_state, pbuf2d, phase=2) + endif + + !----------------------------------------------------------------------- + ! Tendency physics after coupler + ! Not necessary at terminal timestep. + !----------------------------------------------------------------------- + ! +#if ( defined OFFLINE_DYN ) + ! + ! if offline mode set SHFLX QFLX TAUX TAUY for vert diffusion + ! + call get_met_srf2( cam_in ) +#endif + ! lightning flash freq and prod rate of NOx + call t_startf ('lightning_no_prod') + call lightning_no_prod( phys_state, pbuf2d, cam_in ) + call t_stopf ('lightning_no_prod') + + call t_barrierf('sync_ac_physics', mpicom) + call t_startf ('ac_physics') + call t_adj_detailf(+1) + +!$OMP PARALLEL DO PRIVATE (C, NCOL, phys_buffer_chunk) + + do c=begchunk,endchunk + ncol = get_ncols_p(c) + phys_buffer_chunk => pbuf_get_chunk(pbuf2d, c) + ! + ! surface diagnostics for history files + ! + call t_startf('diag_surf') + call diag_surf(cam_in(c), cam_out(c), phys_state(c), phys_buffer_chunk) + call t_stopf('diag_surf') + + call tphysac(ztodt, cam_in(c), & + cam_out(c), & + phys_state(c), phys_tend(c), phys_buffer_chunk) + end do ! Chunk loop + + call t_adj_detailf(-1) + call t_stopf('ac_physics') + +#ifdef TRACER_CHECK + call gmean_mass ('after tphysac FV:WET)', phys_state) +#endif + + call t_startf ('carma_accumulate_stats') + call carma_accumulate_stats() + call t_stopf ('carma_accumulate_stats') + + call t_startf ('physpkg_st2') + call pbuf_deallocate(pbuf2d, 'physpkg') + + call pbuf_update_tim_idx() + call diag_deallocate() + call t_stopf ('physpkg_st2') + + end subroutine phys_run2 + + ! + !----------------------------------------------------------------------- + ! + + subroutine phys_final( phys_state, phys_tend, pbuf2d ) + use physics_buffer, only: physics_buffer_desc, pbuf_deallocate + use chemistry, only: chem_final + use carma_intr, only: carma_final + use wv_saturation, only: wv_sat_final + use microp_aero, only: microp_aero_final + use phys_grid_ctem, only: phys_grid_ctem_final + use nudging, only: Nudge_Model, nudging_final + use hemco_interface, only: HCOI_Chunk_Final + + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Finalization of physics package + ! + !----------------------------------------------------------------------- + ! Input/output arguments + type(physics_state), pointer :: phys_state(:) + type(physics_tend ), pointer :: phys_tend(:) + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + if(associated(pbuf2d)) then + call pbuf_deallocate(pbuf2d,'global') + deallocate(pbuf2d) + end if + deallocate(phys_state) + deallocate(phys_tend) + call chem_final + call carma_final + call wv_sat_final + ! OSLO_AERO begin + ! microp_aero_final() not called + ! OSLO_AERO end + call phys_grid_ctem_final() + if(Nudge_Model) call nudging_final() + + if(use_hemco) then + ! cleanup hemco + call HCOI_Chunk_Final + endif + + end subroutine phys_final + + + subroutine tphysac (ztodt, cam_in, & + cam_out, state, tend, pbuf) + !----------------------------------------------------------------------- + ! + ! Tendency physics after coupling to land, sea, and ice models. + ! + ! Computes the following: + ! + ! o Aerosol Emission at Surface + ! o Stratiform Macro-Microphysics + ! o Wet Scavenging of Aerosol + ! o Radiation + ! o Source-Sink for Advected Tracers + ! o Symmetric Turbulence Scheme - Vertical Diffusion + ! o Rayleigh Friction + ! o Dry Deposition of Aerosol + ! o Enforce Charge Neutrality ( Only for WACCM ) + ! o Gravity Wave Drag + ! o QBO Relaxation ( Only for WACCM ) + ! o Ion Drag ( Only for WACCM ) + ! o Scale Dry Mass Energy + !----------------------------------------------------------------------- + use physics_buffer, only: physics_buffer_desc, pbuf_set_field, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx + use chemistry, only: chem_is_active, chem_timestep_tend, chem_emissions + use cam_diagnostics, only: diag_phys_tend_writeout + use gw_drag, only: gw_tend + use vertical_diffusion, only: vertical_diffusion_tend + use rayleigh_friction, only: rayleigh_friction_tend + use physics_types, only: physics_dme_adjust, set_dry_to_wet, physics_state_check, & + dyn_te_idx + use waccmx_phys_intr, only: waccmx_phys_mspd_tend ! WACCM-X major diffusion + use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_tend ! WACCM-X + use aoa_tracers, only: aoa_tracers_timestep_tend + use physconst, only: rhoh2o + use aero_model, only: aero_model_drydep + use check_energy, only: check_energy_chng, tot_energy_phys + use check_energy, only: check_tracers_data, check_tracers_init, check_tracers_chng + use time_manager, only: get_nstep + use cam_abortutils, only: endrun + use dycore, only: dycore_is + use cam_control_mod, only: aqua_planet + use mo_gas_phase_chemdr,only: map2chm + use clybry_fam, only: clybry_fam_set + use charge_neutrality, only: charge_balance + use qbo, only: qbo_relax + use iondrag, only: iondrag_calc, do_waccm_ions + use perf_mod + use flux_avg, only: flux_avg_run + use cam_history, only: hist_fld_active, outfld + use qneg_module, only: qneg4 + use co2_cycle, only: co2_cycle_set_ptend + use nudging, only: Nudge_Model,Nudge_ON,nudging_timestep_tend + use cam_snapshot, only: cam_snapshot_all_outfld_tphysac + use cam_snapshot_common,only: cam_snapshot_ptend_outfld + use lunar_tides, only: lunar_tides_tend + use ssatcontrail, only: ssatcontrail_d0 + use physics_types, only: physics_ptend_init, physics_ptend_sum, physics_ptend_scale + use microp_driver, only: microp_driver_tend + use microp_aero, only: microp_aero_run + ! OSLO_AERO begin + use oslo_aero_microp, only: oslo_aero_microp_run + use oslo_aero_share + ! OSLO_AERO end + use clubb_intr, only: clubb_tend_cam, clubb_emissions_cam + use subcol, only: subcol_gen, subcol_ptend_avg + use subcol_utils, only: subcol_ptend_copy, is_subcol_on + use subcol_SILHS, only: subcol_SILHS_var_covar_driver, init_state_subcol + use subcol_SILHS, only: subcol_SILHS_fill_holes_conserv + use subcol_SILHS, only: subcol_SILHS_hydromet_conc_tend_lim + use micro_pumas_cam, only: massless_droplet_destroyer + use convect_deep, only: convect_deep_tend_2, deep_scheme_does_scav_trans + use cloud_diagnostics, only: cloud_diagnostics_calc + use radiation, only: radiation_tend + use tropopause, only: tropopause_output + use cam_diagnostics, only: diag_phys_writeout, diag_conv, diag_clip_tend_writeout + use aero_model, only: aero_model_wetdep + use physics_buffer, only: col_type_subcol + use check_energy, only: check_energy_timestep_init + use carma_intr, only: carma_wetdep_tend, carma_timestep_tend, carma_emission_tend + use carma_flags_mod, only: carma_do_aerosol, carma_do_emission, carma_do_detrain + use carma_flags_mod, only: carma_do_cldice, carma_do_cldliq, carma_do_wetdep + use dyn_tests_utils, only: vc_dycore + use cam_thermo, only: cam_thermo_water_update + 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 + ! + ! Arguments + ! + real(r8), intent(in) :: ztodt ! Two times model timestep (2 delta-t) + + type(cam_in_t), intent(inout) :: cam_in + type(cam_out_t), intent(inout) :: cam_out + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + type(physics_buffer_desc), pointer :: pbuf(:) + + + type(check_tracers_data):: tracerint ! tracer mass integrals and cummulative boundary fluxes + + ! + !---------------------------Local workspace----------------------------- + ! + type(physics_ptend) :: ptend ! indivdual parameterization tendencies + type(physics_ptend) :: ptend_macp_all ! sum of macrophysics tendencies (e.g. CLUBB) over substeps + type(physics_state) :: state_sc ! state for sub-columns + type(physics_ptend) :: ptend_sc ! ptend for sub-columns + type(physics_ptend) :: ptend_aero ! ptend for microp_aero + type(physics_ptend) :: ptend_aero_sc ! ptend for microp_aero on sub-columns + type(physics_tend) :: tend_sc ! tend for sub-columns + + integer :: nstep ! current timestep number + real(r8) :: zero(pcols) ! array of zeros + + integer :: lchnk ! chunk identifier + integer :: ncol ! number of atmospheric columns + integer i,k,m ! Longitude, level indices + integer :: yr, mon, day, tod ! components of a date + integer :: ixq, ixcldice, ixcldliq ! constituent indices for vapor, cloud liquid and ice water. + + ! for macro/micro co-substepping + integer :: macmic_it ! iteration variables + real(r8) :: cld_macmic_ztodt ! modified timestep + + real(r8) :: net_flx(pcols) + + real(r8) :: cmfmc(pcols,pverp) ! Convective mass flux--m sub c + + real(r8) dlf(pcols,pver) ! Detraining cld H20 from shallow + deep convections + real(r8) rtdt ! 1./ztodt + + real(r8) :: rliq(pcols) ! vertical integral of liquid not yet in q(ixcldliq) + real(r8) :: det_s (pcols) ! vertical integral of detrained static energy from ice + real(r8) :: det_ice(pcols) ! vertical integral of detrained ice + real(r8) :: flx_cnd(pcols) + + real(r8) :: zero_sc(pcols*psubcols) ! array of zeros + real(r8) :: zero_tracers(pcols,pcnst) + + real(r8), pointer :: dlfzm(:,:) ! ZM detrained convective cloud water mixing ratio. + real(r8), pointer :: cmfmczm(:,:) ! ZM convective mass fluxes + real(r8), pointer :: rliqbc(:) ! tphysbc reserve liquid + + ! stratiform precipitation variables + real(r8),pointer :: prec_str(:) ! sfc flux of precip from stratiform (m/s) + real(r8),pointer :: snow_str(:) ! sfc flux of snow from stratiform (m/s) + real(r8),pointer :: prec_str_sc(:) ! sfc flux of precip from stratiform (m/s) -- for subcolumns + real(r8),pointer :: snow_str_sc(:) ! sfc flux of snow from stratiform (m/s) -- for subcolumns + real(r8),pointer :: prec_pcw(:) ! total precip from prognostic cloud scheme + real(r8),pointer :: snow_pcw(:) ! snow from prognostic cloud scheme + real(r8),pointer :: prec_sed(:) ! total precip from cloud sedimentation + real(r8),pointer :: snow_sed(:) ! snow from cloud ice sedimentation + + ! Local copies for substepping + real(r8) :: prec_pcw_macmic(pcols) + real(r8) :: snow_pcw_macmic(pcols) + real(r8) :: prec_sed_macmic(pcols) + real(r8) :: snow_sed_macmic(pcols) + + ! carma precipitation variables + real(r8) :: prec_sed_carma(pcols) ! total precip from cloud sedimentation (CARMA) + real(r8) :: snow_sed_carma(pcols) ! snow from cloud ice sedimentation (CARMA) + + logical :: labort ! abort flag + + real(r8) tvm(pcols,pver) ! virtual temperature + real(r8) prect(pcols) ! total precipitation + real(r8) surfric(pcols) ! surface friction velocity + real(r8) obklen(pcols) ! Obukhov length + real(r8) :: fh2o(pcols) ! h2o flux to balance source from methane chemistry + real(r8) :: flx_heat(pcols) ! Heat flux for check_energy_chng. + real(r8) :: tmp_trac (pcols,pver,pcnst) ! tmp space + real(r8) :: tmp_pdel (pcols,pver) ! tmp space + real(r8) :: tmp_ps (pcols) ! tmp space + real(r8) :: scaling(pcols,pver) + logical :: moist_mixing_ratio_dycore + + ! physics buffer fields for total energy and mass adjustment + integer itim_old, ifld + + real(r8), pointer, dimension(:,:) :: cld + real(r8), pointer, dimension(:,:) :: qini + real(r8), pointer, dimension(:,:) :: cldliqini + real(r8), pointer, dimension(:,:) :: cldiceini + real(r8), pointer, dimension(:,:) :: totliqini + real(r8), pointer, dimension(:,:) :: toticeini + real(r8), pointer, dimension(:,:) :: dtcore + real(r8), pointer, dimension(:,:) :: dqcore + real(r8), pointer, dimension(:,:) :: ducore + real(r8), pointer, dimension(:,:) :: dvcore + real(r8), pointer, dimension(:,:) :: ast ! relative humidity cloud fraction + + !----------------------------------------------------------------------- + lchnk = state%lchnk + ncol = state%ncol + + nstep = get_nstep() + rtdt = 1._r8/ztodt + + ! Adjust the surface fluxes to reduce instabilities in near sfc layer + if (phys_do_flux_avg()) then + call flux_avg_run(state, cam_in, pbuf, nstep, ztodt) + endif + + ! Validate the physics state. + if (state_debug_checks) then + call physics_state_check(state, name="before tphysac") + end if + + call t_startf('tphysac_init') + ! Associate pointers with physics buffer fields + itim_old = pbuf_old_tim_idx() + + call pbuf_get_field(pbuf, dtcore_idx, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, dqcore_idx, dqcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, ducore_idx, ducore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, dvcore_idx, dvcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + + call pbuf_get_field(pbuf, qini_idx, qini) + call pbuf_get_field(pbuf, cldliqini_idx, cldliqini) + call pbuf_get_field(pbuf, cldiceini_idx, cldiceini) + call pbuf_get_field(pbuf, totliqini_idx, totliqini) + call pbuf_get_field(pbuf, toticeini_idx, toticeini) + + ifld = pbuf_get_index('CLD') + call pbuf_get_field(pbuf, ifld, cld, start=(/1,1,itim_old/),kount=(/pcols,pver,1/)) + + ifld = pbuf_get_index('AST') + call pbuf_get_field(pbuf, ifld, ast, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + + call cnst_get_ind('Q', ixq) + call cnst_get_ind('CLDLIQ', ixcldliq) + call cnst_get_ind('CLDICE', ixcldice) + + call pbuf_get_field(pbuf, prec_str_idx, prec_str ) + call pbuf_get_field(pbuf, snow_str_idx, snow_str ) + call pbuf_get_field(pbuf, prec_sed_idx, prec_sed ) + call pbuf_get_field(pbuf, snow_sed_idx, snow_sed ) + call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw ) + call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw ) + + if (is_subcol_on()) then + call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol) + call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol) + end if + + if (dlfzm_idx > 0) then + call pbuf_get_field(pbuf, dlfzm_idx, dlfzm) + dlf(:ncol,:) = dlfzm(:ncol,:) + else + dlf(:,:) = 0._r8 + end if + + if (cmfmczm_idx > 0) then + call pbuf_get_field(pbuf, cmfmczm_idx, cmfmczm) + cmfmc(:ncol,:) = cmfmczm(:ncol,:) + else + cmfmc(:ncol,:) = 0._r8 + end if + + call pbuf_get_field(pbuf, rliqbc_idx, rliqbc) + rliq(:ncol) = rliqbc(:ncol) + + ! + ! accumulate fluxes into net flux array for spectral dycores + ! jrm Include latent heat of fusion for snow + ! + do i=1,ncol + tend%flx_net(i) = tend%flx_net(i) + cam_in%shf(i) + (cam_out%precc(i) & + + cam_out%precl(i))*latvap*rhoh2o & + + (cam_out%precsc(i) + cam_out%precsl(i))*latice*rhoh2o + end do + + ! emissions of aerosols and gas-phase chemistry constituents at surface + + if (trim(cam_take_snapshot_before) == "chem_emissions") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call chem_emissions( state, cam_in, pbuf ) + if (trim(cam_take_snapshot_after) == "chem_emissions") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + if (carma_do_emission) then + ! carma emissions + call carma_emission_tend (state, ptend, cam_in, ztodt) + call physics_update(state, ptend, ztodt, tend) + end if + + ! get nstep and zero array for energy checker + zero = 0._r8 + zero_sc(:) = 0._r8 + zero_tracers(:,:) = 0._r8 + nstep = get_nstep() + call check_tracers_init(state, tracerint) + + ! Check if latent heat flux exceeds the total moisture content of the + ! lowest model layer, thereby creating negative moisture. + + call qneg4('TPHYSAC', lchnk, ncol, ztodt , & + state%q(1,pver,1), state%rpdel(1,pver), & + cam_in%shf, cam_in%lhf, cam_in%cflx) + + call t_stopf('tphysac_init') + + !=================================================== + ! Apply tracer surface fluxes to lowest model layer + !=================================================== + call t_startf('clubb_emissions_tend') + + call clubb_emissions_cam(state, cam_in, ptend) + + call physics_update(state, ptend, ztodt, tend) + + call check_energy_chng(state, tend, "clubb_emissions_tend", nstep, ztodt, zero, zero, zero, zero) + + call t_stopf('clubb_emissions_tend') + + !=================================================== + ! Calculate tendencies from CARMA bin microphysics. + !=================================================== + ! + ! If CARMA is doing detrainment, then on output, rliq no longer represents + ! water reserved + ! for detrainment, but instead represents potential snow fall. The mass and + ! number of the + ! snow are stored in the physics buffer and will be incorporated by the MG + ! microphysics. + ! + ! Currently CARMA cloud microphysics is only supported with the MG + ! microphysics. + call t_startf('carma_timestep_tend') + + if (carma_do_cldice .or. carma_do_cldliq) then + call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, dlf=dlf, rliq=rliq, & + prec_str=prec_str, snow_str=snow_str, prec_sed=prec_sed_carma, snow_sed=snow_sed_carma) + call physics_update(state, ptend, ztodt, tend) + + ! Before the detrainment, the reserved condensate is all liquid, but if + ! CARMA is doing + ! detrainment, then the reserved condensate is snow. + if (carma_do_detrain) then + call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str+rliq, snow_str+rliq, zero) + else + call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, prec_str, snow_str, zero) + end if + end if + + call t_stopf('carma_timestep_tend') + + if( microp_scheme == 'MG' ) then + ! Start co-substepping of macrophysics and microphysics + cld_macmic_ztodt = ztodt/cld_macmic_num_steps + + ! Clear precip fields that should accumulate. + prec_sed_macmic = 0._r8 + snow_sed_macmic = 0._r8 + prec_pcw_macmic = 0._r8 + snow_pcw_macmic = 0._r8 + + ! contrail parameterization + ! see Chen et al., 2012: Global contrail coverage simulated + ! by CAM5 with the inventory of 2006 global aircraft emissions, JAMES + ! https://doi.org/10.1029/2011MS000105 + call ssatcontrail_d0(state, pbuf, ztodt, ptend) + call physics_update(state, ptend, ztodt, tend) + + ! initialize ptend structures where macro and microphysics tendencies are + ! accumulated over macmic substeps + call physics_ptend_init(ptend_macp_all,state%psetcols,'macrophysics',lu=.true.,lv=.true.) + + do macmic_it = 1, cld_macmic_num_steps + + !=================================================== + ! Calculate macrophysical tendency (sedimentation, detrain, cloud fraction) + !=================================================== + + call t_startf('macrop_tend') + + ! ===================================================== + ! CLUBB call (PBL, shallow convection, macrophysics) + ! ===================================================== + + if (trim(cam_take_snapshot_before) == "clubb_tend_cam") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call clubb_tend_cam(state, ptend, pbuf, cld_macmic_ztodt,& + cmfmc, cam_in, macmic_it, cld_macmic_num_steps, & + dlf, det_s, det_ice) + + ! Since we "added" the reserved liquid back in this routine, we need + ! to account for it in the energy checker + flx_cnd(:ncol) = -1._r8*rliq(:ncol) + flx_heat(:ncol) = cam_in%shf(:ncol) + det_s(:ncol) + + ! Unfortunately, physics_update does not know what time period + ! "tend" is supposed to cover, and therefore can't update it + ! with substeps correctly. For now, work around this by scaling + ! ptend down by the number of substeps, then applying it for + ! the full time (ztodt). + call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) + + ! Update physics tendencies and copy state to state_eq, because that is + ! input for microphysics + if ( (trim(cam_take_snapshot_after) == "clubb_tend_cam") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + call physics_ptend_sum(ptend,ptend_macp_all,ncol) + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "clubb_tend_cam") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + ! Use actual qflux (not lhf/latvap) for consistency with surface fluxes and revised code + call check_energy_chng(state, tend, "clubb_tend", nstep, ztodt, & + cam_in%cflx(:ncol,1)/cld_macmic_num_steps, & + flx_cnd(:ncol)/cld_macmic_num_steps, & + det_ice(:ncol)/cld_macmic_num_steps, & + flx_heat(:ncol)/cld_macmic_num_steps) + + call t_stopf('macrop_tend') + + !=================================================== + ! Calculate cloud microphysics + !=================================================== + + if (is_subcol_on() .neqv. use_subcol_microp ) then + call endrun("Error calculating cloud microphysics: is_subcol_on() != use_subcol_microp") + end if + + if (is_subcol_on()) then + ! Allocate sub-column structures. + call physics_state_alloc(state_sc, lchnk, psubcols*pcols) + call physics_tend_alloc(tend_sc, psubcols*pcols) + + ! Generate sub-columns using the requested scheme + if (trim(subcol_scheme) == 'SILHS') call init_state_subcol(state, tend, state_sc, tend_sc) + call subcol_gen(state, tend, state_sc, tend_sc, pbuf) + + !Initialize check energy for subcolumns + call check_energy_timestep_init(state_sc, tend_sc, pbuf, col_type_subcol) + end if + + if (trim(cam_take_snapshot_before) == "microp_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + ! OSLO_AERO begin + call t_startf('oslo_aero_microp_run') + call oslo_aero_microp_run(state, ptend_aero, cld_macmic_ztodt, pbuf) + call t_stopf('oslo_aero_microp_run') + ! OSLO_AERO end + + call t_startf('microp_tend') + + if (use_subcol_microp) then + + if (trim(cam_take_snapshot_before) == "microp_driver_tend_subcol") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state_sc, tend_sc, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call microp_driver_tend(state_sc, ptend_sc, cld_macmic_ztodt, pbuf) + ! Parameterize subcolumn effects on covariances, if enabled + if (trim(subcol_scheme) == 'SILHS') & + call subcol_SILHS_var_covar_driver( cld_macmic_ztodt, state_sc, ptend_sc, pbuf ) + + ! Average the sub-column ptend for use in gridded update - will not contain ptend_aero + call subcol_ptend_avg(ptend_sc, state_sc%ngrdcol, lchnk, ptend) + + ! Call the conservative hole filler. + ! Hole filling is only necessary when using subcolumns. + ! Note: this needs to be called after subcol_ptend_avg but before + ! physics_ptend_scale. + if (trim(subcol_scheme) == 'SILHS') & + call subcol_SILHS_fill_holes_conserv( state, cld_macmic_ztodt, & + ptend, pbuf ) + + ! Destroy massless droplets - Note this routine returns with no change unless + ! micro_do_massless_droplet_destroyer has been set to true + call massless_droplet_destroyer( cld_macmic_ztodt, state, & ! Intent(in) + ptend ) ! Intent(inout) + + ! Limit the value of hydrometeor concentrations in order to place + ! reasonable limits on hydrometeor drop size and keep them from + ! becoming too large. + ! Note: this needs to be called after hydrometeor mixing ratio + ! tendencies are adjusted by subcol_SILHS_fill_holes_conserv + ! and after massless drop concentrations are removed by the + ! subcol_SILHS_massless_droplet_destroyer, but before the + ! call to physics_ptend_scale. + if (trim(subcol_scheme) == 'SILHS') & + call subcol_SILHS_hydromet_conc_tend_lim( state, cld_macmic_ztodt, ptend ) + + ! Copy ptend_aero field to one dimensioned by sub-columns before summing with ptend + call subcol_ptend_copy(ptend_aero, state_sc, ptend_aero_sc) + call physics_ptend_sum(ptend_aero_sc, ptend_sc, state_sc%ncol) + call physics_ptend_dealloc(ptend_aero_sc) + + ! Have to scale and apply for full timestep to get tend right + ! (see above note for macrophysics). + call physics_ptend_scale(ptend_sc, 1._r8/cld_macmic_num_steps, ncol) + + if ( (trim(cam_take_snapshot_after) == "microp_driver_tend_subcol") .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_sc, ptend_sc, ztodt, tend_sc) + + if (trim(cam_take_snapshot_after) == "microp_driver_tend_subcol") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state_sc, tend_sc, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call check_energy_chng(state_sc, tend_sc, "microp_tend_subcol", & + nstep, ztodt, zero_sc, & + prec_str_sc(:state_sc%ncol)/cld_macmic_num_steps, & + snow_str_sc(:state_sc%ncol)/cld_macmic_num_steps, zero_sc) + + call physics_state_dealloc(state_sc) + call physics_tend_dealloc(tend_sc) + call physics_ptend_dealloc(ptend_sc) + else + call microp_driver_tend(state, ptend, cld_macmic_ztodt, pbuf) + end if + ! combine aero and micro tendencies for the grid + call physics_ptend_sum(ptend_aero, ptend, ncol) + call physics_ptend_dealloc(ptend_aero) + + ! Have to scale and apply for full timestep to get tend right + ! (see above note for macrophysics). + call physics_ptend_scale(ptend, 1._r8/cld_macmic_num_steps, ncol) + + call diag_clip_tend_writeout(state, ptend, ncol, lchnk, ixcldliq, ixcldice, ixq, ztodt, rtdt) + + if ( (trim(cam_take_snapshot_after) == "microp_section") .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) == "microp_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call check_energy_chng(state, tend, "microp_tend", nstep, ztodt, & + zero, prec_str(:ncol)/cld_macmic_num_steps, & + snow_str(:ncol)/cld_macmic_num_steps, zero) + + call t_stopf('microp_tend') + + prec_sed_macmic(:ncol) = prec_sed_macmic(:ncol) + prec_sed(:ncol) + snow_sed_macmic(:ncol) = snow_sed_macmic(:ncol) + snow_sed(:ncol) + prec_pcw_macmic(:ncol) = prec_pcw_macmic(:ncol) + prec_pcw(:ncol) + snow_pcw_macmic(:ncol) = snow_pcw_macmic(:ncol) + snow_pcw(:ncol) + + end do ! end substepping over macrophysics/microphysics + + call outfld( 'UTEND_MACROP', ptend_macp_all%u, pcols, lchnk) + call outfld( 'VTEND_MACROP', ptend_macp_all%v, pcols, lchnk) + call physics_ptend_dealloc(ptend_macp_all) + + prec_sed(:ncol) = prec_sed_macmic(:ncol)/cld_macmic_num_steps + snow_sed(:ncol) = snow_sed_macmic(:ncol)/cld_macmic_num_steps + prec_pcw(:ncol) = prec_pcw_macmic(:ncol)/cld_macmic_num_steps + snow_pcw(:ncol) = snow_pcw_macmic(:ncol)/cld_macmic_num_steps + prec_str(:ncol) = prec_pcw(:ncol) + prec_sed(:ncol) + snow_str(:ncol) = snow_pcw(:ncol) + snow_sed(:ncol) + + endif + + ! Add the precipitation from CARMA to the precipitation from stratiform. + if (carma_do_cldice .or. carma_do_cldliq) then + prec_sed(:ncol) = prec_sed(:ncol) + prec_sed_carma(:ncol) + snow_sed(:ncol) = snow_sed(:ncol) + snow_sed_carma(:ncol) + end if + + if ( .not. deep_scheme_does_scav_trans() ) then + + ! ------------------------------------------------------------------------------- + ! 1. Wet Scavenging of Aerosols by Convective and Stratiform Precipitation. + ! 2. Convective Transport of Non-Water Aerosol Species. + ! + ! . Aerosol wet chemistry determines scavenging fractions, and transformations + ! . Then do convective transport of all trace species except qv,ql,qi. + ! . We needed to do the scavenging first to determine the interstitial fraction. + ! . When UNICON is used as unified convection, we should still perform + ! wet scavenging but not 'convect_deep_tend2'. + ! ------------------------------------------------------------------------------- + + call t_startf('bc_aerosols') + if (clim_modal_aero .and. .not. prog_modal_aero) then + call modal_aero_calcsize_diag(state, pbuf) + call modal_aero_wateruptake_dr(state, pbuf) + endif + + if (trim(cam_take_snapshot_before) == "aero_model_wetdep") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call aero_model_wetdep( state, ztodt, dlf, cam_out, ptend, pbuf) + if ( (trim(cam_take_snapshot_after) == "aero_model_wetdep") .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) == "aero_model_wetdep") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + if (carma_do_wetdep) then + ! CARMA wet deposition + ! + ! NOTE: It needs to follow aero_model_wetdep, so that + ! cam_out%xxxwetxxx + ! fields have already been set for CAM aerosols and cam_out can be + ! added + ! to for CARMA aerosols. + call t_startf ('carma_wetdep_tend') + call carma_wetdep_tend(state, ptend, ztodt, pbuf, dlf, cam_out) + call physics_update(state, ptend, ztodt, tend) + call t_stopf ('carma_wetdep_tend') + end if + + call t_startf ('convect_deep_tend2') + call convect_deep_tend_2( state, ptend, ztodt, pbuf ) + call physics_update(state, ptend, ztodt, tend) + call t_stopf ('convect_deep_tend2') + + ! check tracer integrals + call check_tracers_chng(state, tracerint, "cmfmca", nstep, ztodt, zero_tracers) + + call t_stopf('bc_aerosols') + + endif + + !=================================================== + ! Moist physical parameteriztions complete: + ! send dynamical variables, and derived variables to history file + !=================================================== + + call t_startf('bc_history_write') + call diag_phys_writeout(state, pbuf) + call diag_conv(state, ztodt, pbuf) + + call t_stopf('bc_history_write') + + !=================================================== + ! Write cloud diagnostics on history file + !=================================================== + + call t_startf('bc_cld_diag_history_write') + + call cloud_diagnostics_calc(state, pbuf) + + call t_stopf('bc_cld_diag_history_write') + + !=================================================== + ! Radiation computations + !=================================================== + call t_startf('radiation') + + if (trim(cam_take_snapshot_before) == "radiation_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call radiation_tend( & + state, ptend, pbuf, cam_out, cam_in, net_flx) + + ! Set net flux used by spectral dycores + do i=1,ncol + tend%flx_net(i) = net_flx(i) + end do + + if ( (trim(cam_take_snapshot_after) == "radiation_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) == "radiation_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call check_energy_chng(state, tend, "radheat", nstep, ztodt, zero, zero, zero, net_flx) + + call t_stopf('radiation') + + ! Diagnose the location of the tropopause and its location to the history file(s). + call t_startf('tropopause') + call tropopause_output(state) + call t_stopf('tropopause') + + !=================================================== + ! Source/sink terms for advected tracers. + !=================================================== + call t_startf('adv_tracer_src_snk') + ! Test tracers + + if (trim(cam_take_snapshot_before) == "aoa_tracers_timestep_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call aoa_tracers_timestep_tend(state, ptend, cam_in%cflx, cam_in%landfrac, ztodt) + if ( (trim(cam_take_snapshot_after) == "aoa_tracers_timestep_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) == "aoa_tracers_timestep_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call check_tracers_chng(state, tracerint, "aoa_tracers_timestep_tend", nstep, ztodt, & + cam_in%cflx) + + if (trim(cam_take_snapshot_before) == "co2_cycle_set_ptend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call co2_cycle_set_ptend(state, pbuf, ptend) + if ( (trim(cam_take_snapshot_after) == "co2_cycle_set_ptend") .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) == "co2_cycle_set_ptend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + !=================================================== + ! Chemistry and MAM calculation + ! MAM core aerosol conversion process is performed in the below 'chem_timestep_tend'. + ! In addition, surface flux of aerosol species other than 'dust' and 'sea salt', and + ! elevated emission of aerosol species are treated in 'chem_timestep_tend' before + ! Gas chemistry and MAM core aerosol conversion. + ! Note that surface flux is not added into the atmosphere, but elevated emission is + ! added into the atmosphere as tendency. + !=================================================== + if (chem_is_active()) then + + if (trim(cam_take_snapshot_before) == "chem_timestep_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call chem_timestep_tend(state, ptend, cam_in, cam_out, ztodt, & + pbuf, fh2o=fh2o) + + + if ( (trim(cam_take_snapshot_after) == "chem_timestep_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) == "chem_timestep_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call check_energy_chng(state, tend, "chem", nstep, ztodt, fh2o, zero, zero, zero) + call check_tracers_chng(state, tracerint, "chem_timestep_tend", nstep, ztodt, & + cam_in%cflx) + end if + call t_stopf('adv_tracer_src_snk') + + !=================================================== + ! Vertical diffusion/pbl calculation + ! Call vertical diffusion (apply tracer emissions, molecular diffusion and pbl form drag) + !=================================================== + + call t_startf('vertical_diffusion_tend') + + if (trim(cam_take_snapshot_before) == "vertical_diffusion_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call vertical_diffusion_tend (ztodt ,state , cam_in, & + surfric ,obklen ,ptend ,ast ,pbuf ) + + !------------------------------------------ + ! Call major diffusion for extended model + !------------------------------------------ + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then + call waccmx_phys_mspd_tend (ztodt ,state ,ptend) + endif + + if ( (trim(cam_take_snapshot_after) == "vertical_diffusion_section") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + if ( ptend%lu ) then + call outfld( 'UTEND_VDIFF', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_VDIFF', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "vertical_diffusion_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call t_stopf ('vertical_diffusion_tend') + + !=================================================== + ! Rayleigh friction calculation + !=================================================== + call t_startf('rayleigh_friction') + call rayleigh_friction_tend( ztodt, state, ptend) + if ( ptend%lu ) then + call outfld( 'UTEND_RAYLEIGH', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_RAYLEIGH', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + call t_stopf('rayleigh_friction') + + if (do_clubb_sgs) then + call check_energy_chng(state, tend, "vdiff", nstep, ztodt, zero, zero, zero, zero) + else + call check_energy_chng(state, tend, "vdiff", nstep, ztodt, cam_in%cflx(:,1), zero, & + zero, cam_in%shf) + endif + + call check_tracers_chng(state, tracerint, "vdiff", nstep, ztodt, cam_in%cflx) + + ! aerosol dry deposition processes + call t_startf('aero_drydep') + + if (trim(cam_take_snapshot_before) == "aero_model_drydep") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call aero_model_drydep( state, pbuf, obklen, surfric, cam_in, ztodt, cam_out, ptend ) + if ( (trim(cam_take_snapshot_after) == "aero_model_drydep") .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) == "aero_model_drydep") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call t_stopf('aero_drydep') + + ! CARMA microphysics + ! + ! NOTE: This does both the timestep_tend for CARMA aerosols as well as doing + ! the dry + ! deposition for CARMA aerosols. It needs to follow vertical_diffusion_tend, + ! so that + ! obklen and surfric have been calculated. It needs to follow + ! aero_model_drydep, so + ! that cam_out%xxxdryxxx fields have already been set for CAM aerosols and + ! cam_out + ! can be added to for CARMA aerosols. + if (carma_do_aerosol) then + call t_startf('carma_timestep_tend') + call carma_timestep_tend(state, cam_in, cam_out, ptend, ztodt, pbuf, obklen=obklen, ustar=surfric) + call physics_update(state, ptend, ztodt, tend) + + call check_energy_chng(state, tend, "carma_tend", nstep, ztodt, zero, zero, zero, zero) + call t_stopf('carma_timestep_tend') + end if + + !--------------------------------------------------------------------------------- + ! ... enforce charge neutrality + !--------------------------------------------------------------------------------- + call charge_balance(state, pbuf) + + !=================================================== + ! Gravity wave drag + !=================================================== + call t_startf('gw_tend') + + if (trim(cam_take_snapshot_before) == "gw_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call gw_tend(state, pbuf, ztodt, ptend, cam_in, flx_heat) + + if ( (trim(cam_take_snapshot_after) == "gw_tend") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + if ( ptend%lu ) then + call outfld( 'UTEND_GWDTOT', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_GWDTOT', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "gw_tend") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + ! Check energy integrals + call check_energy_chng(state, tend, "gwdrag", nstep, ztodt, zero, & + zero, zero, flx_heat) + call t_stopf('gw_tend') + + ! QBO relaxation + + if (trim(cam_take_snapshot_before) == "qbo_relax") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call qbo_relax(state, pbuf, ptend) + if ( (trim(cam_take_snapshot_after) == "qbo_relax") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + if ( ptend%lu ) then + call outfld( 'UTEND_QBORLX', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_QBORLX', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "qbo_relax") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + ! Check energy integrals + call check_energy_chng(state, tend, "qborelax", nstep, ztodt, zero, zero, zero, zero) + + ! Lunar tides + call lunar_tides_tend( state, ptend ) + if ( ptend%lu ) then + call outfld( 'UTEND_LUNART', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_LUNART', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + ! Check energy integrals + call check_energy_chng(state, tend, "lunar_tides", nstep, ztodt, zero, zero, zero, zero) + + ! Ion drag calculation + call t_startf ( 'iondrag' ) + + if (trim(cam_take_snapshot_before) == "iondrag_calc_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + if ( do_waccm_ions ) then + call iondrag_calc( lchnk, ncol, state, ptend, pbuf, ztodt ) + else + call iondrag_calc( lchnk, ncol, state, ptend) + endif + !---------------------------------------------------------------------------- + ! Call ionosphere routines for extended model if mode is set to ionosphere + !---------------------------------------------------------------------------- + if( waccmx_is('ionosphere') ) then + call waccmx_phys_ion_elec_temp_tend(state, ptend, pbuf, ztodt) + endif + + if ( (trim(cam_take_snapshot_after) == "iondrag_calc_section") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + if ( ptend%lu ) then + call outfld( 'UTEND_IONDRG', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_IONDRG', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "iondrag_calc_section") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call tot_energy_phys(state, 'phAP') + call tot_energy_phys(state, 'dyAP',vc=vc_dycore) + + !--------------------------------------------------------------------------------- + ! Enforce charge neutrality after O+ change from ionos_tend + !--------------------------------------------------------------------------------- + if( waccmx_is('ionosphere') ) then + call charge_balance(state, pbuf) + endif + + ! Check energy integrals + call check_energy_chng(state, tend, "iondrag", nstep, ztodt, zero, zero, zero, zero) + + call t_stopf ( 'iondrag' ) + + ! Update Nudging values, if needed + !---------------------------------- + if((Nudge_Model).and.(Nudge_ON)) then + call nudging_timestep_tend(state,ptend) + if ( ptend%lu ) then + call outfld( 'UTEND_NDG', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_NDG', ptend%v, pcols, lchnk) + end if + call physics_update(state,ptend,ztodt,tend) + call check_energy_chng(state, tend, "nudging", nstep, ztodt, zero, zero, zero, zero) + endif + + !-------------- Energy budget checks vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv + ! Save total energy for global fixer in next timestep + ! + ! This call must be after the last parameterization and call to physics_update + ! + call pbuf_set_field(pbuf, teout_idx, state%te_cur(:,dyn_te_idx), (/1,itim_old/),(/pcols,1/)) + ! + ! 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') + ! + ! update cp/cv for energy computation based in updated water variables + ! + call cam_thermo_water_update(state%q(:ncol,:,:), lchnk, ncol, vc_dycore,& + to_dry_factor=state%pdel(:ncol,:)/state%pdeldry(:ncol,:)) + + ! for dry mixing ratio dycore, physics_dme_adjust is called for energy diagnostic purposes only. + ! So, save off tracers + if (.not.moist_mixing_ratio_dycore) then + ! + ! for dry-mixing ratio based dycores dme_adjust takes place in the dynamical core + ! + ! only compute dme_adjust for diagnostics purposes + ! + if (thermo_budget_history) then + tmp_trac(:ncol,:pver,:pcnst) = state%q(:ncol,:pver,:pcnst) + tmp_pdel(:ncol,:pver) = state%pdel(:ncol,:pver) + tmp_ps(:ncol) = state%ps(:ncol) + call physics_dme_adjust(state, tend, qini, totliqini, toticeini, ztodt) + call tot_energy_phys(state, 'phAM') + call tot_energy_phys(state, 'dyAM', vc=vc_dycore) + ! Restore pre-"physics_dme_adjust" tracers + state%q(:ncol,:pver,:pcnst) = tmp_trac(:ncol,:pver,:pcnst) + state%pdel(:ncol,:pver) = tmp_pdel(:ncol,:pver) + state%ps(:ncol) = tmp_ps(:ncol) + end if + 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) + + if (trim(cam_take_snapshot_before) == "physics_dme_adjust") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + call physics_dme_adjust(state, tend, qini, totliqini, toticeini, ztodt) + if (trim(cam_take_snapshot_after) == "physics_dme_adjust") then + call cam_snapshot_all_outfld_tphysac(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf,& + fh2o, surfric, obklen, flx_heat, cmfmc, dlf, det_s, det_ice, net_flx) + end if + + call tot_energy_phys(state, 'phAM') + call tot_energy_phys(state, 'dyAM', vc=vc_dycore) + endif + + 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 + ! are computed correctly; e.g. dtcore) + ! + scaling(1:ncol,:) = cpairv(:ncol,:,lchnk)/cp_or_cv_dycore(:ncol,:,lchnk) + state%T(1:ncol,:) = state%temp_ini(1:ncol,:)+& + scaling(1:ncol,:)*(state%T(1:ncol,:)-state%temp_ini(1:ncol,:)) + tend%dtdt(:ncol,:) = scaling(:ncol,:)*tend%dtdt(:ncol,:) + ! + ! else: do nothing for dycores with energy consistent with CAM physics + ! + end if + + + ! store T, U, and V in buffer for use in computing dynamics T-tendency in next timestep + do k = 1,pver + dtcore(:ncol,k) = state%t(:ncol,k) + dqcore(:ncol,k) = state%q(:ncol,k,ixq) + ducore(:ncol,k) = state%u(:ncol,k) + dvcore(:ncol,k) = state%v(:ncol,k) + end do + + !-------------- Energy budget checks ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + if (aqua_planet) then + labort = .false. + do i=1,ncol + if (cam_in%ocnfrac(i) /= 1._r8) then + labort = .true. + if (masterproc) write(iulog,*) 'oceanfrac(',i,')=',cam_in%ocnfrac(i) + end if + end do + if (labort) then + call endrun ('TPHYSAC error: in aquaplanet mode, but grid contains non-ocean point') + endif + endif + + call diag_phys_tend_writeout (state, pbuf, tend, ztodt, qini, cldliqini, cldiceini) + + call clybry_fam_set( ncol, lchnk, map2chm, state%q, pbuf ) + + end subroutine tphysac + + subroutine tphysbc (ztodt, state, & + tend, pbuf, & + cam_out, cam_in ) + !----------------------------------------------------------------------- + ! + ! Purpose: + ! Evaluate and apply physical processes that are calculated BEFORE + ! coupling to land, sea, and ice models. + ! + ! Processes currently included are: + ! + ! o Resetting Negative Tracers to Positive + ! o Global Mean Total Energy Fixer + ! o Dry Adjustment + ! o Asymmetric Turbulence Scheme - Deep Convection & Shallow Convection + ! + ! Method: + ! + ! Each parameterization should be implemented with this sequence of calls: + ! 1) Call physics interface + ! 2) Check energy + ! 3) Call physics_update + ! See Interface to Column Physics and Chemistry Packages + ! http://www.ccsm.ucar.edu/models/atm-cam/docs/phys-interface/index.html + ! + !----------------------------------------------------------------------- + + use physics_buffer, only: physics_buffer_desc, pbuf_get_field + use physics_buffer, only: pbuf_get_index, pbuf_old_tim_idx + use physics_buffer, only: col_type_subcol, dyn_time_lvls + + use dadadj_cam, only: dadadj_tend + use physics_types, only: physics_update, & + physics_state_check, & + dyn_te_idx + use cam_diagnostics, only: diag_conv_tend_ini, diag_conv, diag_export, diag_state_b4_phys_write + use cam_history, only: outfld + use constituents, only: qmin + use air_composition, only: thermodynamic_active_species_liq_num,thermodynamic_active_species_liq_idx + use air_composition, only: thermodynamic_active_species_ice_num,thermodynamic_active_species_ice_idx + use convect_deep, only: convect_deep_tend + use time_manager, only: is_first_step, get_nstep + use convect_diagnostics,only: convect_diagnostics_calc + use check_energy, only: check_energy_chng, check_energy_fix + use check_energy, only: check_tracers_data, check_tracers_init + use check_energy, only: tot_energy_phys + use dycore, only: dycore_is + use radiation, only: radiation_tend + use perf_mod + use mo_gas_phase_chemdr,only: map2chm + use clybry_fam, only: clybry_fam_adj + use cam_abortutils, only: endrun + use subcol_utils, only: is_subcol_on + use qneg_module, only: qneg3 + use cam_snapshot, only: cam_snapshot_all_outfld_tphysbc + use cam_snapshot_common, only: cam_snapshot_ptend_outfld + use dyn_tests_utils, only: vc_dycore + + ! Arguments + + real(r8), intent(in) :: ztodt ! 2 delta t (model time increment) + + type(physics_state), intent(inout) :: state + type(physics_tend ), intent(inout) :: tend + type(physics_buffer_desc), pointer :: pbuf(:) + + type(cam_out_t), intent(inout) :: cam_out + type(cam_in_t), intent(in) :: cam_in + + + ! + !---------------------------Local workspace----------------------------- + ! + + type(physics_ptend) :: ptend ! indivdual parameterization tendencies + + integer :: nstep ! current timestep number + + real(r8) :: net_flx(pcols) + + real(r8) :: zdu(pcols,pver) ! detraining mass flux from deep convection + real(r8) :: cmfmc(pcols,pverp) ! Convective mass flux--m sub c + + real(r8) cmfcme(pcols,pver) ! cmf condensation - evaporation + + real(r8) dlf(pcols,pver) ! Detraining cld H20 from shallow + deep convections + real(r8) dlf2(pcols,pver) ! Detraining cld H20 from shallow convections + real(r8) rtdt ! 1./ztodt + + integer lchnk ! chunk identifier + integer ncol ! number of atmospheric columns + + integer :: i ! column indicex + integer :: ixcldice, ixcldliq, ixq ! constituent indices for cloud liquid and ice water. + integer :: m, m_cnst + + ! physics buffer fields to compute tendencies for stratiform package + integer itim_old, ifld + real(r8), pointer, dimension(:,:) :: cld ! cloud fraction + + ! physics buffer fields for total energy and mass adjustment + real(r8), pointer, dimension(: ) :: teout + real(r8), pointer, dimension(:,:) :: qini + real(r8), pointer, dimension(:,:) :: cldliqini + real(r8), pointer, dimension(:,:) :: cldiceini + real(r8), pointer, dimension(:,:) :: totliqini + real(r8), pointer, dimension(:,:) :: toticeini + real(r8), pointer, dimension(:,:) :: dtcore + real(r8), pointer, dimension(:,:) :: dqcore + real(r8), pointer, dimension(:,:) :: ducore + real(r8), pointer, dimension(:,:) :: dvcore + + real(r8), pointer, dimension(:,:,:) :: fracis ! fraction of transported species that are insoluble + + real(r8), pointer :: dlfzm(:,:) ! ZM detrained convective cloud water mixing ratio. + real(r8), pointer :: rliqbc(:) ! tphysbc reserve liquid + + ! convective precipitation variables + real(r8),pointer :: prec_dp(:) ! total precipitation from ZM convection + real(r8),pointer :: snow_dp(:) ! snow from ZM convection + real(r8),pointer :: prec_sh(:) ! total precipitation from Hack convection + real(r8),pointer :: snow_sh(:) ! snow from Hack convection + + ! stratiform precipitation variables + real(r8),pointer :: prec_str(:) ! sfc flux of precip from stratiform (m/s) + real(r8),pointer :: snow_str(:) ! sfc flux of snow from stratiform (m/s) + real(r8),pointer :: prec_str_sc(:) ! sfc flux of precip from stratiform (m/s) -- for subcolumns + real(r8),pointer :: snow_str_sc(:) ! sfc flux of snow from stratiform (m/s) -- for subcolumns + real(r8),pointer :: prec_pcw(:) ! total precip from prognostic cloud scheme + real(r8),pointer :: snow_pcw(:) ! snow from prognostic cloud scheme + real(r8),pointer :: prec_sed(:) ! total precip from cloud sedimentation + real(r8),pointer :: snow_sed(:) ! snow from cloud ice sedimentation + + ! energy checking variables + real(r8) :: zero(pcols) ! array of zeros + real(r8) :: zero_sc(pcols*psubcols) ! array of zeros + real(r8) :: rliq(pcols) ! vertical integral of liquid not yet in q(ixcldliq) + real(r8) :: rice(pcols) ! vertical integral of ice not yet in q(ixcldice) + real(r8) :: rliq2(pcols) ! vertical integral of liquid from shallow scheme + real(r8) :: flx_cnd(pcols) + real(r8) :: flx_heat(pcols) + type(check_tracers_data):: tracerint ! energy integrals and cummulative boundary fluxes + real(r8) :: zero_tracers(pcols,pcnst) + + logical :: lq(pcnst) + + !----------------------------------------------------------------------- + + call t_startf('bc_init') + + zero = 0._r8 + zero_tracers(:,:) = 0._r8 + zero_sc(:) = 0._r8 + + lchnk = state%lchnk + ncol = state%ncol + + rtdt = 1._r8/ztodt + + nstep = get_nstep() + + ! Associate pointers with physics buffer fields + itim_old = pbuf_old_tim_idx() + ifld = pbuf_get_index('CLD') + call pbuf_get_field(pbuf, ifld, cld, (/1,1,itim_old/),(/pcols,pver,1/)) + + call pbuf_get_field(pbuf, teout_idx, teout, (/1,itim_old/), (/pcols,1/)) + + call pbuf_get_field(pbuf, qini_idx, qini) + call pbuf_get_field(pbuf, cldliqini_idx, cldliqini) + call pbuf_get_field(pbuf, cldiceini_idx, cldiceini) + call pbuf_get_field(pbuf, totliqini_idx, totliqini) + call pbuf_get_field(pbuf, toticeini_idx, toticeini) + + call pbuf_get_field(pbuf, dtcore_idx, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, dqcore_idx, dqcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, ducore_idx, ducore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + call pbuf_get_field(pbuf, dvcore_idx, dvcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + + ifld = pbuf_get_index('FRACIS') + call pbuf_get_field(pbuf, ifld, fracis, start=(/1,1,1/), kount=(/pcols, pver, pcnst/) ) + fracis (:ncol,:,1:pcnst) = 1._r8 + + ! Set physics tendencies to 0 + tend%dTdt(:ncol,:pver) = 0._r8 + tend%dudt(:ncol,:pver) = 0._r8 + tend%dvdt(:ncol,:pver) = 0._r8 + + ! Verify state coming from the dynamics + if (state_debug_checks) then + call physics_state_check(state, name="before tphysbc (dycore?)") + end if + + call clybry_fam_adj( ncol, lchnk, map2chm, state%q, pbuf ) + + ! Since clybry_fam_adj operates directly on the tracers, and has no + ! physics_update call, re-run qneg3. + call qneg3('TPHYSBCc',lchnk ,ncol ,pcols ,pver , & + 1, pcnst, qmin ,state%q ) + + ! Validate output of clybry_fam_adj. + if (state_debug_checks) then + call physics_state_check(state, name="clybry_fam_adj") + end if + ! + ! Dump out "before physics" state + ! + call diag_state_b4_phys_write (state) + + ! compute mass integrals of input tracers state + call check_tracers_init(state, tracerint) + + call t_stopf('bc_init') + + !=================================================== + ! Global mean total energy fixer + !=================================================== + call t_startf('energy_fixer') + + call tot_energy_phys(state, 'phBF') + call tot_energy_phys(state, 'dyBF',vc=vc_dycore) + + if (.not.dycore_is('EUL')) then + call check_energy_fix(state, ptend, nstep, flx_heat) + call physics_update(state, ptend, ztodt, tend) + call check_energy_chng(state, tend, "chkengyfix", nstep, ztodt, zero, zero, zero, flx_heat) + call outfld( 'EFIX', flx_heat , pcols, lchnk ) + end if + + call tot_energy_phys(state, 'phBP') + call tot_energy_phys(state, 'dyBP',vc=vc_dycore) + ! Save state for convective tendency calculations. + call diag_conv_tend_ini(state, pbuf) + + call cnst_get_ind('Q', ixq) + call cnst_get_ind('CLDLIQ', ixcldliq) + call cnst_get_ind('CLDICE', ixcldice) + qini (:ncol,:pver) = state%q(:ncol,:pver, ixq) + cldliqini(:ncol,:pver) = state%q(:ncol,:pver,ixcldliq) + cldiceini(:ncol,:pver) = state%q(:ncol,:pver,ixcldice) + + totliqini(:ncol,:pver) = 0.0_r8 + do m_cnst=1,thermodynamic_active_species_liq_num + m = thermodynamic_active_species_liq_idx(m_cnst) + totliqini(:ncol,:pver) = totliqini(:ncol,:pver)+state%q(:ncol,:pver,m) + end do + toticeini(:ncol,:pver) = 0.0_r8 + do m_cnst=1,thermodynamic_active_species_ice_num + m = thermodynamic_active_species_ice_idx(m_cnst) + toticeini(:ncol,:pver) = toticeini(:ncol,:pver)+state%q(:ncol,:pver,m) + end do + + + call outfld('TEOUT', teout , pcols, lchnk ) + call outfld('TEINP', state%te_ini(:,dyn_te_idx), pcols, lchnk ) + call outfld('TEFIX', state%te_cur(:,dyn_te_idx), pcols, lchnk ) + + ! T, U, V tendency due to dynamics + if ( nstep > dyn_time_lvls-1 ) then + dtcore(:ncol,:pver) = (state%t(:ncol,:pver) - dtcore(:ncol,:pver))/ztodt + dqcore(:ncol,:pver) = (state%q(:ncol,:pver,ixq) - dqcore(:ncol,:pver))/ztodt + ducore(:ncol,:pver) = (state%u(:ncol,:pver) - ducore(:ncol,:pver))/ztodt + dvcore(:ncol,:pver) = (state%v(:ncol,:pver) - dvcore(:ncol,:pver))/ztodt + call outfld( 'DTCORE', dtcore, pcols, lchnk ) + call outfld( 'DQCORE', dqcore, pcols, lchnk ) + call outfld( 'UTEND_CORE', ducore, pcols, lchnk ) + call outfld( 'VTEND_CORE', dvcore, pcols, lchnk ) + end if + + call t_stopf('energy_fixer') + ! + !=================================================== + ! Dry adjustment + !=================================================== + call t_startf('dry_adjustment') + + if (trim(cam_take_snapshot_before) == "dadadj_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call dadadj_tend(ztodt, state, ptend) + + if ( (trim(cam_take_snapshot_after) == "dadadj_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) == "dadadj_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call t_stopf('dry_adjustment') + + !=================================================== + ! Moist convection + !=================================================== + call t_startf('moist_convection') + + call t_startf ('convect_deep_tend') + + if (trim(cam_take_snapshot_before) == "convect_deep_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call convect_deep_tend( & + cmfmc, cmfcme, & + zdu, & + rliq, rice, & + ztodt, & + state, ptend, cam_in%landfrac, pbuf) + + if ( (trim(cam_take_snapshot_after) == "convect_deep_tend") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + + if ( ptend%lu ) then + call outfld( 'UTEND_DCONV', ptend%u, pcols, lchnk) + end if + if ( ptend%lv ) then + call outfld( 'VTEND_DCONV', ptend%v, pcols, lchnk) + end if + call physics_update(state, ptend, ztodt, tend) + + if (trim(cam_take_snapshot_after) == "convect_deep_tend") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_after_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + + call t_stopf('convect_deep_tend') + + call pbuf_get_field(pbuf, prec_dp_idx, prec_dp ) + call pbuf_get_field(pbuf, snow_dp_idx, snow_dp ) + call pbuf_get_field(pbuf, prec_sh_idx, prec_sh ) + call pbuf_get_field(pbuf, snow_sh_idx, snow_sh ) + + call pbuf_get_field(pbuf, prec_str_idx, prec_str ) + call pbuf_get_field(pbuf, snow_str_idx, snow_str ) + call pbuf_get_field(pbuf, prec_sed_idx, prec_sed ) + call pbuf_get_field(pbuf, snow_sed_idx, snow_sed ) + call pbuf_get_field(pbuf, prec_pcw_idx, prec_pcw ) + call pbuf_get_field(pbuf, snow_pcw_idx, snow_pcw ) + + if (use_subcol_microp) then + call pbuf_get_field(pbuf, prec_str_idx, prec_str_sc, col_type=col_type_subcol) + call pbuf_get_field(pbuf, snow_str_idx, snow_str_sc, col_type=col_type_subcol) + end if + + ! Check energy integrals, including "reserved liquid" + flx_cnd(:ncol) = prec_dp(:ncol) + rliq(:ncol) + snow_dp(:ncol) = snow_dp(:ncol) + rice(:ncol) + call check_energy_chng(state, tend, "convect_deep", nstep, ztodt, zero, flx_cnd, snow_dp, zero) + snow_dp(:ncol) = snow_dp(:ncol) - rice(:ncol) + + !=================================================== + ! Compute convect diagnostics + !=================================================== + + if (dlfzm_idx > 0) then + call pbuf_get_field(pbuf, dlfzm_idx, dlfzm) + dlf(:ncol,:) = dlfzm(:ncol,:) + else + dlf(:,:) = 0._r8 + end if + + if (trim(cam_take_snapshot_before) == "convect_diagnostics_calc") then + call cam_snapshot_all_outfld_tphysbc(cam_snapshot_before_num, state, tend, cam_in, cam_out, pbuf, & + cmfmc, cmfcme, zdu, rliq, rice, dlf, dlf2, rliq2, net_flx) + end if + call convect_diagnostics_calc (ztodt , cmfmc, & + dlf , dlf2 , rliq , rliq2, & + state , pbuf) + if ( (trim(cam_take_snapshot_after) == "convect_diagnostics_calc") .and. & + (trim(cam_take_snapshot_before) == trim(cam_take_snapshot_after))) then + call cam_snapshot_ptend_outfld(ptend, lchnk) + end if + + ! add reserve liquid to pbuf + call pbuf_get_field(pbuf, rliqbc_idx, rliqbc) + rliqbc(:ncol) = rliq(:ncol) + + call t_stopf('moist_convection') + + if (is_first_step()) then + + !initiailize sedimentation arrays + prec_pcw = 0._r8 + snow_pcw = 0._r8 + prec_sed = 0._r8 + snow_sed = 0._r8 + prec_str = 0._r8 + snow_str = 0._r8 + + if (is_subcol_on()) then + prec_str_sc = 0._r8 + snow_str_sc = 0._r8 + end if + + ! OSLO_AERO begin + !=================================================== + ! Run wet deposition routines to intialize aerosols + ! NOT CALLED IN OSLO AERO + !=================================================== + ! OSLO_AERO end + + !=================================================== + ! Radiation computations + ! initialize fluxes only, do not update state + !=================================================== + + call radiation_tend( & + state, ptend, pbuf, cam_out, cam_in, net_flx) + + end if + + ! Save atmospheric fields to force surface models + call t_startf('cam_export') + call cam_export (state,cam_out,pbuf) + call t_stopf('cam_export') + + ! Write export state to history file + call t_startf('diag_export') + call diag_export(cam_out) + call t_stopf('diag_export') + + end subroutine tphysbc + +subroutine phys_timestep_init(phys_state, cam_in, cam_out, pbuf2d) +!----------------------------------------------------------------------------------- +! +! Purpose: The place for parameterizations to call per timestep initializations. +! Generally this is used to update time interpolated fields from boundary +! datasets. +! +!----------------------------------------------------------------------------------- + use chemistry, only: chem_timestep_init + use chem_surfvals, only: chem_surfvals_set + use physics_types, only: physics_state + use physics_buffer, only: physics_buffer_desc + use carma_intr, only: carma_timestep_init + use ghg_data, only: ghg_data_timestep_init + use aoa_tracers, only: aoa_tracers_timestep_init + use vertical_diffusion, only: vertical_diffusion_ts_init + use radheat, only: radheat_timestep_init + use solar_data, only: solar_data_advance + use qbo, only: qbo_timestep_init + use iondrag, only: do_waccm_ions, iondrag_timestep_init + use perf_mod + + use prescribed_ozone, only: prescribed_ozone_adv + use prescribed_ghg, only: prescribed_ghg_adv + use prescribed_aero, only: prescribed_aero_adv + use aerodep_flx, only: aerodep_flx_adv + use aircraft_emit, only: aircraft_emit_adv + use prescribed_volcaero, only: prescribed_volcaero_adv + use prescribed_strataero,only: prescribed_strataero_adv + use mo_apex, only: mo_apex_init + use epp_ionization, only: epp_ionization_active + use iop_forcing, only: scam_use_iop_srf + use nudging, only: Nudge_Model, nudging_timestep_init + use waccmx_phys_intr, only: waccmx_phys_ion_elec_temp_timestep_init + use phys_grid_ctem, only: phys_grid_ctem_diags + ! OSLO_AERO begin + use oslo_aero_ocean, only: oslo_aero_ocean_adv + ! OSLO_AERO end + + implicit none + + type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state + type(cam_in_t), intent(inout), dimension(begchunk:endchunk) :: cam_in + type(cam_out_t), intent(inout), dimension(begchunk:endchunk) :: cam_out + + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + !----------------------------------------------------------------------------- + + if (single_column) call scam_use_iop_srf(cam_in) + + ! update geomagnetic coordinates + if (epp_ionization_active .or. do_waccm_ions) then + call mo_apex_init(phys_state) + endif + + ! Chemistry surface values + call chem_surfvals_set() + + ! Solar irradiance + call solar_data_advance() + + ! Time interpolate for chemistry. + call chem_timestep_init(phys_state, pbuf2d) + + if( waccmx_is('ionosphere') ) then + call waccmx_phys_ion_elec_temp_timestep_init(phys_state,pbuf2d) + endif + + ! Prescribed tracers + call prescribed_ozone_adv(phys_state, pbuf2d) + call prescribed_ghg_adv(phys_state, pbuf2d) + call prescribed_aero_adv(phys_state, pbuf2d) + call aircraft_emit_adv(phys_state, pbuf2d) + call prescribed_volcaero_adv(phys_state, pbuf2d) + call prescribed_strataero_adv(phys_state, pbuf2d) + ! OSLO_AERO begin + call oslo_aero_ocean_adv(phys_state, pbuf2d) + ! OSLO_AERO end + + ! prescribed aerosol deposition fluxes + call aerodep_flx_adv(phys_state, pbuf2d, cam_out) + + ! Time interpolate data models of gasses in pbuf2d + call ghg_data_timestep_init(pbuf2d, phys_state) + + ! Upper atmosphere radiative processes + call radheat_timestep_init(phys_state, pbuf2d) + + ! Time interpolate for vertical diffusion upper boundary condition + call vertical_diffusion_ts_init(pbuf2d, phys_state) + + !---------------------------------------------------------------------- + ! update QBO data for this time step + !---------------------------------------------------------------------- + call qbo_timestep_init + + call iondrag_timestep_init() + + call carma_timestep_init() + + ! age of air tracers + call aoa_tracers_timestep_init(phys_state) + + ! Update Nudging values, if needed + !---------------------------------- + if(Nudge_Model) call nudging_timestep_init(phys_state) + + ! Update TEM diagnostics + call phys_grid_ctem_diags(phys_state) + +end subroutine phys_timestep_init + +end module physpkg diff --git a/src_cam/radiation.F90 b/src_cam/radiation.F90 index 143cb67..4f8b8d7 100644 --- a/src_cam/radiation.F90 +++ b/src_cam/radiation.F90 @@ -46,9 +46,8 @@ module radiation use prescribed_volcaero, only: has_prescribed_volcaero use oslo_aero_optical_params, only: oslo_aero_optical_params_calc use oslo_aero_share, only: nmodes_oslo=>nmodes -#ifdef AEROCOM +use oslo_aero_control, only: use_aerocom use oslo_aero_aerocom, only: dod440, dod550, dod870, abs550, abs550alt -#endif ! OSLO_AERO end implicit none @@ -1317,13 +1316,13 @@ subroutine radiation_tend( & call outfld('FSNS_DRF',fsns , pcols, lchnk) call outfld('FSNTCDRF',rd%fsntc , pcols, lchnk) call outfld('FSNSCDRF',rd%fsnsc , pcols, lchnk) -#ifdef AEROCOM - call outfld('FSUTADRF',rd%fsutoa(:) , pcols, lchnk) - call outfld('FSDS_DRF',fsds(:) , pcols, lchnk) - ftem_1d(1:ncol) = fsds(1:ncol)-fsns(1:ncol) - call outfld('FSUS_DRF',ftem_1d , pcols, lchnk) - call outfld('FSDSCDRF',rd%fsdsc(:) , pcols, lchnk) -#endif + if (use_aerocom) then + call outfld('FSUTADRF',rd%fsutoa(:) , pcols, lchnk) + call outfld('FSDS_DRF',fsds(:) , pcols, lchnk) + ftem_1d(1:ncol) = fsds(1:ncol)-fsns(1:ncol) + call outfld('FSUS_DRF',ftem_1d , pcols, lchnk) + call outfld('FSDSCDRF',rd%fsdsc(:) , pcols, lchnk) + end if call rad_rrtmg_sw( & lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, &