diff --git a/oslo_aero_share.F90 b/oslo_aero_share.F90 deleted file mode 100644 index d8974a0..0000000 --- a/oslo_aero_share.F90 +++ /dev/null @@ -1,1075 +0,0 @@ -module oslo_aero_share - - !--------------------------------------------------------------------------------- - ! Module to set up register aerosols indexes, number of gas and particle - ! species and their scavenging rates. Tables for humidity growth - !--------------------------------------------------------------------------------- - - use shr_kind_mod, only: r8 => shr_kind_r8 - use ppgrid, only: pcols, pver - use constituents, only: pcnst, cnst_name, cnst_get_ind - use mo_tracname, only: solsym - use cam_abortutils, only: endrun - use physics_buffer, only: physics_buffer_desc, pbuf_get_field - use physconst, only: pi - ! - implicit none - public ! Make default type private to the module - - !--------------------------- - ! Public interfaces - !--------------------------- - - public :: aero_register ! register consituents - public :: is_process_mode ! Check is an aerosol specie is a process mode - public :: isAerosol ! Check is specie is aerosol (i.e. gases get .FALSE. here) - public :: getTracerIndex - public :: getNumberOfTracersInMode - public :: getNumberOfBackgroundTracersInMode - public :: getCloudTracerIndex - public :: getCloudTracerIndexDirect - public :: getCloudTracerName - public :: chemistryIndex - public :: physicsIndex - public :: getDryDensity - public :: getConstituentFraction - public :: isTracerInMode - public :: fillAerosolTracerList - public :: getNumberOfAerosolTracers - public :: fillInverseAerosolTracerList - public :: qqcw_get_field - ! - public :: calculateNumberConcentration - public :: calculateNumberMedianRadius - public :: calculateEquivalentDensityOfFractalMode - public :: calculatedNdLogR - public :: calculateLognormalCDF - ! - public :: init_interp_constants - - !--------------------------- - ! Public parameters - !--------------------------- - - ! Define lognormal size parameters for each size mode (dry, at point of emission/production) - integer, parameter :: nmodes = 14 - integer, parameter :: nbmodes = 10 - integer, parameter :: nbands = 14 ! number of aerosol spectral bands in SW - integer, parameter :: nlwbands = 16 ! number of aerosol spectral bands in LW - integer, parameter :: nbmp1 = 11 ! number of first non-background mode - integer, parameter :: max_tracers_per_mode = 7 - - integer, parameter :: numberOfExternallyMixedModes = 4 !Modes 0;11-14 (13 is not used in lifecycle) - integer, parameter :: numberOfInternallyMIxedMOdes = 9 !Modes 1-10 (3 is not used in lifecycle) - integer, parameter :: numberOfProcessModeTracers = 6 - - integer, parameter :: nBinsTab = 44 ![nbr] number of tabulated bins - - real(r8), parameter :: rTabMin = 1.e-9_r8 ![m] smallest lookup table size - real(r8), parameter :: rTabMax = 20.e-6_r8 ![m] largest lookup table size - real(r8), parameter :: rMinAquousChemistry = 0.05e-6_r8 !Smallest particle which can receive aquous chemistry mass - real(r8), parameter :: sq2pi = 1._r8/sqrt(2.0_r8*pi) - real(r8), parameter :: smallConcentration = 1.e-100_r8 !duplicate, sync with smallNumber in Const - real(r8), parameter :: smallNumber = 1.e-100_r8 - - integer, parameter :: n_tracers_in_mode(0:nmodes) = (/ 1, 4, 3, 0, 5, 7, 7, 7, 7, 7, 7, 0, 1, 0, 2 /) - integer, parameter :: n_background_tracers_in_mode(0:nmodes) = (/ 1,2,1,0,2,1,1,1,1,1,1,0,1,0,2 /) - - ! Define aerosol types and their properties.. - integer, parameter :: N_AEROSOL_TYPES = 5 - integer, parameter :: AEROSOL_TYPE_SULFATE = 1 - integer, parameter :: AEROSOL_TYPE_BC = 2 - integer, parameter :: AEROSOL_TYPE_OM = 3 - integer, parameter :: AEROSOL_TYPE_DUST = 4 - integer, parameter :: AEROSOL_TYPE_SALT = 5 - - ! Define aerosol modes - integer, parameter :: MODE_IDX_BC_EXT_AC = 0 !Externally mixed BC accumulation mode - integer, parameter :: MODE_IDX_SO4SOA_AIT = 1 !SO4 and SOA in aitken mode, Created from 11 by growth (condensation) of SO4 - integer, parameter :: MODE_IDX_BC_AIT = 2 !Created from 12 by growth (condensation) SO4 - integer, parameter :: MODE_IDX_NOT_USED = 3 !Not used - integer, parameter :: MODE_IDX_OMBC_INTMIX_COAT_AIT = 4 !Created from 14 by growth (condensation) of SO4 and from cloud processing/wet-phas - integer, parameter :: MODE_IDX_SO4_AC = 5 !Accumulation mode SO4 (mode will have other comps added) - integer, parameter :: MODE_IDX_DST_A2 = 6 !Accumulation mode dust (mode will have other comps added) - integer, parameter :: MODE_IDX_DST_A3 = 7 !Coarse mode dust (mode will have other comps added) - integer, parameter :: MODE_IDX_SS_A1 = 8 !Fine mode sea-salt (mode will have other comps added) - integer, parameter :: MODE_IDX_SS_A2 = 9 !Accumulation mode sea-salt (mode will have other comps added) - integer, parameter :: MODE_IDX_SS_A3 = 10 !Coarse mode sea-salt (mode will have other comps added) - integer, parameter :: MODE_IDX_SO4SOA_NUC = 11 !SO4 and SOA nucleation mode - integer, parameter :: MODE_IDX_BC_NUC = 12 !BC nucleation mode - integer, parameter :: MODE_IDX_LUMPED_ORGANICS = 13 !not used in lifecycle, but some extra mass goes here when max. allowed LUT conc. are too small - integer, parameter :: MODE_IDX_OMBC_INTMIX_AIT = 14 !mix quickly formed in fire-plumes - - ! Number median radius of background emissions THESE DO NOT ASSUME IMPLICIT GROWTH!! - real(r8), parameter :: originalNumberMedianRadius(0:nmodes) = & - 1.e-6_r8* (/ 0.0626_r8, & !0 - 0.0118_r8, 0.024_r8, 0.04_r8, 0.04_r8, 0.075_r8, & !1-5 - 0.22_r8, 0.63_r8, 0.0475_r8, 0.30_r8, 0.75_r8, & !6-10 ! SS: Salter et al. (2015) - 0.0118_r8, 0.024_r8, 0.04_r8, 0.04_r8 /) !11-14 - - ! sigma of background aerosols - real(r8), parameter :: originalSigma(0:nmodes) = & - (/1.6_r8, & !0 - 1.8_r8, 1.8_r8, 1.8_r8, 1.8_r8, 1.59_r8, & !1-5 - 1.59_r8, 2.0_r8, 2.1_r8, 1.72_r8, 1.60_r8, & !6-10 ! SS: Salter et al. (2015) - 1.8_r8, 1.8_r8, 1.8_r8, 1.8_r8 /) !11-14 - - !Radius used for the modes in the lifeCycle MAY ASSUME SOME GROWTH ALREADY HAPPENED - real(r8), parameter :: lifeCycleNumberMedianRadius(0:nmodes) = & - 1.e-6_r8*(/ 0.0626_r8, 0.025_r8, 0.025_r8 , 0.04_r8, 0.06_r8, 0.075_r8, & - 0.22_r8, 0.63_r8, 0.0475_r8, 0.30_r8, 0.75_r8, & ! Salter et al. (2015) - 0.0118_r8, 0.024_r8, 0.04_r8 , 0.04_r8 /) - - !Sigma based on original lifecycle code (taken from "sigmak" used previously in lifecycle code) - real(r8), parameter :: lifeCycleSigma(0:nmodes) = & - (/1.6_r8, 1.8_r8, 1.8_r8, 1.8_r8, 1.8_r8, & !0-4 - 1.59_r8, 1.59_r8, 2.0_r8, & !5,6,7 (SO4+dust) - 2.1_r8, 1.72_r8, 1.6_r8, & !8-10 (SS) ! Salter et al. (2015) - 1.8_r8, 1.8_r8, 1.8_r8, 1.8_r8/) !11-14 - - !Below cloud scavenging coefficients for modes which have an actual size - real(r8), parameter :: belowCloudScavengingCoefficient(0:nmodes) = & - (/ 0.01_r8 , 0.02_r8 , 0.02_r8 , 0.0_r8 , 0.02_r8, 0.01_r8, & !(0-5) - 0.02_r8 , 0.2_r8 , 0.02_r8 , 0.02_r8, 0.5_r8, & !6-10 (DUST+SS) - 0.04_r8 , 0.08_r8 , 0.0_r8 , 0.02_r8 /) ! SO4_n, bc_n, N/A og bc/oc - - ! Treatment of process-modes! - ! The tracers indices can not be set here since they are not known on compile time - ! tracerInProcessMode = (/l_so4_a1, l_so4_a2, l_so4_ac, l_om_ac, l_bc_ac, l_soa_a1 /) - - !The process modes need an "efficient size" (Why does A1 have a different size than the others??) - real(r8), parameter :: processModeNumberMedianRadius(numberOfProcessModeTracers) = & - (/ 0.04e-6_r8, 0.1e-6_r8, 0.1e-6_r8, 0.1e-6_r8, 0.1e-6_r8, 0.04e-6_r8 /) - - !The process modes need an "efficient sigma" - real(r8), parameter :: processModeSigma(numberOfProcessModeTracers) = & - (/ 1.8_r8, 1.59_r8, 1.59_r8, 1.59_r8, 1.59_r8, 1.8_r8 /) - - real(r8), parameter :: belowCloudScavengingCoefficientProcessModes(numberOfProcessModeTracers) = & - (/0.02_r8, 0.01_r8, 0.02_r8, 0.02_r8, 0.02_r8, 0.02_r8 /) - - real(r8), parameter :: aThird = 1.0_r8/3.0_r8 - - real(r8), parameter :: e=2.718281828_r8 - real(r8), parameter :: eps=1.0e-30_r8 - - !--------------------------- - ! Public constants - !--------------------------- - - real(r8) :: nk(0:nmodes,nbinsTab) !dN/dlogr for modes - real(r8) :: normnk(0:nmodes,nbinsTab) !dN for modes (sums to one over size range) - real(r8) :: rBinEdge(nBinsTab+1) - real(r8) :: rBinMidpoint(nBinsTab) - real(r8) :: volumeToNumber(0:nmodes) !m3 ==> # - real(r8) :: numberToSurface(0:nmodes) !# ==> m2 - - ! Constants used in interpolation - ! Internal mixtures of process-tagged mass - ! cate : total added mass (µg/m3 per particle per cm3) from condensation - ! and wet phase chemistry/cloud processing, for kcomp = 1-2. - ! cate should be scaled up/down whenever the modal parameters (modal - ! radius and width) are increased/decreased a lot. - ! cat : total added mass (µg/m3 per particle per cm-3) from coagulation, condensation - ! and wet phase chemistry/cloud processing, for kcomp = 5-10. - ! cat should be scaled up/down whenever the modal parameters (modal - ! radius and width) are increased/decreased a lot. - ! fac : mass fraction of cat or cate from coagulating carbonaceous aerosols (BC+OM). - ! The remaining mass cate*(1-fac) or cat*(1-fac) is SO4. - ! fbc : mass fraction of BC from coagulating carbonaceous aerosols, BC/(BC+OM). - ! faq : mass fraction of sulfate which is produced in wet-phase, SO4aq/SO4. - ! The remaining SO4 mass, SO4*(1-faq), is from condensation. - - real(r8) :: rh(10) - real(r8) :: fombg(6), fbcbg(6), fac(6), fbc(6), faq(6) - real(r8) :: cate(4,16) - real(r8) :: cat(5:10,6) - - ! relative humidity (RH, as integer for output variable names) for use in AeroCom code - integer :: RF(6) - - ! AeroCom specific RH input variables for use in opticsAtConstRh.F90 - integer :: irhrf1(6) - real(r8) :: xrhrf(6) - - integer :: tracerInProcessMode(numberOfProcessModeTracers) - integer :: processModeMap(pcnst) - - ! NUMBERS BELOW ARE ESSENTIAL TO CALCULATE HYGROSCOPICITY AND THEREFORE INDIRECT EFFECT! - ! These numbers define the "hygroscopicity parameter" Numbers are selected so that they give reasonable hygroscipity - ! note that changing numbers individually changes the hygroscopicity! - ! Hygroscopicity is defined in Abdul-Razzak and S. Ghan: (B in their eqn 4) - ! A parameterization of aerosol activation 2. Multiple aerosol types, JGR, vol 105, noD5, pp 6837 - ! http://onlinelibrary.wiley.com/doi/10.1029/1999JD901161/abstract - ! - ! Further note that changing any of these numbers without changing aerotab will lead to - ! inconsistencies in the simulation since Aerotab tabulates hygroscopical growth! - ! - ! Main reference for numbers chosen: Ghan et al MIRAGE paper (JRG, vol 106, D6, pp 5295), 2001 References: - ! SULFATE : Using same numbers as MIRAGE paper (ammonium sulfate) - ! BC : Does not really matter as long as soluble mass fraction is small - ! However, numbers below reproduces values from MIRAGE paper - ! New mass density (October 2016) is based on Bond and Bergstrom (2007): Light Absorption - ! by Carbonaceous Particles: An Investigative Review, Aerosol Science and Technology, 40:27���67. - ! OM : Soluble mass fraction tuned to give B of MIRAGE Paper - ! DUST : The numbers give B of ~ 0.07 (high end of Kohler, Kreidenweis et al, GRL, vol 36, 2009. - ! (10% as soluble mass fraction seems reasonable) - ! (see also Osada et al, Atmospheric Research, vol 124, 2013, pp 101 - ! SEA SALT: Soluble mass fraction tuned to give consistent values for (r/r0) at 99% when using the parametrization in - ! Koepke, Hess, Schult and Shettle: Max-Plack-Institut fur Meteorolgie, report No. 243 "GLOBAL AEROSOL DATA SET" - ! These values give "B" of 1.20 instead of 1.16 in MIRAGE paper. - - character(len=8) :: aerosol_type_name(N_AEROSOL_TYPES) = & - (/"SULFATE ", "BC ","OM ", "DUST ", "SALT " /) - real(r8) :: aerosol_type_density(N_AEROSOL_TYPES) = & - (/1769.0_r8, 1800.0_r8, 1500.0_r8, 2600.0_r8, 2200.0_r8 /) !kg/m3 - real(r8) :: aerosol_type_molecular_weight(N_AEROSOL_TYPES) = & - (/132.0_r8, 12.0_r8, 168.2_r8, 135.0_r8, 58.44_r8 /) !kg/kmol - real(r8) :: aerosol_type_osmotic_coefficient(N_AEROSOL_TYPES) = & - (/0.7_r8, 1.111_r8, 1.0_r8, 1.0_r8, 1.0_r8 /) ![-] - real(r8) :: aerosol_type_soluble_mass_fraction(N_AEROSOL_TYPES) = & - (/1.0_r8, 1.67e-7_r8, 0.8725_r8, 0.1_r8, 0.885_r8 /) ![-] - real(r8) :: aerosol_type_number_of_ions(N_AEROSOL_TYPES) = & - (/3.0_r8, 1.0_r8, 1.0_r8, 2.0_r8, 2.0_r8 /) ![-] - - real(r8) :: rhopart(pcnst) - real(r8) :: sgpart(pcnst) - real(r8) :: osmoticCoefficient(pcnst) - real(r8) :: numberOfIons(pcnst) - real(r8) :: solubleMassFraction(pcnst) - integer :: aerosolType(pcnst) - real(r8) :: numberFractionAvailableAqChem(nbmodes) - real(r8) :: invrhopart(pcnst) - - !These tables describe how the tracers behave chemically - integer, dimension(numberOfExternallyMixedModes) :: externallyMixedMode = & - (/MODE_IDX_BC_EXT_AC, & - MODE_IDX_SO4SOA_NUC, & - MODE_IDX_BC_NUC, & - MODE_IDX_OMBC_INTMIX_AIT /) - - integer, dimension(numberOfInternallyMixedMOdes) :: internallyMixedMode = & - (/MODE_IDX_SO4SOA_AIT, & - MODE_IDX_BC_AIT, & - MODE_IDX_OMBC_INTMIX_COAT_AIT, & - MODE_IDX_SO4_AC, & - MODE_IDX_DST_A2, & - MODE_IDX_DST_A3, & - MODE_IDX_SS_A1, & - MODE_IDX_SS_A2, & - MODE_IDX_SS_A3 /) - - ! species indices for individual camuio species - integer :: l_so4_na, l_so4_a1, l_so4_a2, l_so4_ac - integer :: l_bc_n, l_bc_ax, l_bc_ni, l_bc_a, l_bc_ai,l_bc_ac - integer :: l_om_ni, l_om_ai, l_om_ac - integer :: l_so4_pr - integer :: l_dst_a2, l_dst_a3 - integer :: l_ss_a1, l_ss_a2, l_ss_a3, l_h2so4 - integer :: l_soa_na, l_soa_a1, l_soa_lv, l_soa_sv - - integer :: n_aerosol_tracers !number of aerosol tracers - integer :: imozart - - ! Number of transported tracers in each mode - integer :: tracer_in_mode(0:nmodes, max_tracers_per_mode) - - ! Growth of aerosols, duplicated in oslo_aero_sw_tables - real(r8) :: rdivr0(10,pcnst) - - real(r8) :: rhtab(10) = (/ 0.0_r8, 0.37_r8, 0.47_r8, 0.65_r8, 0.75_r8, 0.80_r8, 0.85_r8, 0.90_r8, 0.95_r8, 0.98_r8/) - - integer :: cloudTracerIndex(pcnst) - character(len=20) :: cloudTracerName(pcnst) - - integer, private :: qqcw(pcnst)=-1 ! Remaps modal_aero indices into pbuf - -!=============================================================================== -contains -!=============================================================================== - - function is_process_mode(l_index_in, isChemistry) result(answer) - - !----------------------------------------------------------------------- - ! For a tracer in an aerosol mode, check if this isactually a real - ! tracer or a process mode - !----------------------------------------------------------------------- - - integer, intent(in) :: l_index_in - logical, intent(in) :: isChemistry !true if called from chemistry - - integer :: l_index_phys - logical :: answer - - l_index_phys = l_index_in - if (isChemistry) then - l_index_phys = l_index_phys + iMozart - 1 - endif - - ! return true if tracer is a "process mode" - answer = .false. - if(l_index_phys .eq. l_so4_a1 .or. & - l_index_phys .eq. l_so4_a2 .or. & - l_index_phys .eq. l_so4_ac .or. & - l_index_phys .eq. l_bc_ac .or. & - l_index_phys .eq. l_om_ac .or. & - l_index_phys .eq. l_soa_a1 ) then - answer = .true. - endif - end function is_process_mode - - !=============================================================================== - - subroutine aero_register - - !----------------------------------------------------------------------- - ! Register aerosol modes and indices, should be changed to read in values - ! instead of hard-coding it. - !----------------------------------------------------------------------- - - use mpishorthand - use physics_buffer, only: pbuf_add_field, dtype_r8 - use ppgrid, only: pcols, pver, pverp - - integer :: idx_dum, l,m,mm - logical :: isAlreadyCounted(pcnst) - - ! register the species - - call cnst_get_ind('SO4_NA' ,l_so4_na, abort=.true.) !Aitken mode sulfate (growth from so4_n) - call cnst_get_ind('SO4_A1' ,l_so4_a1, abort=.true.) !sulfate condensate (gas phase production) - call cnst_get_ind('SO4_A2' ,l_so4_a2, abort=.true.) !sulfate produced in aq. chemistry - call cnst_get_ind('SO4_AC' ,l_so4_ac, abort=.true.) !sulfate from coagulation processes - call cnst_get_ind('SO4_PR' ,l_so4_pr, abort=.true.) !sulfate emitted as primary - - call cnst_get_ind('BC_N' ,l_bc_n, abort=.true.) !emissions (mainly industry) lost through coagulation - call cnst_get_ind('BC_AX' ,l_bc_ax, abort=.true.) !externally mixed (fluffy and impossible to activate) - call cnst_get_ind('BC_NI' ,l_bc_ni, abort=.true.) !mixed with oc (mainly biomass), externally mixed otherwise (before condensation etc) - call cnst_get_ind('BC_A' ,l_bc_a, abort=.true.) !formed when bc_n grows by condensation - call cnst_get_ind('BC_AI' ,l_bc_ai, abort=.true.) !formed when bc_ni grows by condensation - call cnst_get_ind('BC_AC' ,l_bc_ac, abort=.true.) !bc from coagulation processes - - call cnst_get_ind('OM_NI' ,l_om_ni, abort=.true.) !om (mainly from biomass), emitted - call cnst_get_ind('OM_AI' ,l_om_ai, abort=.true.) !om formed when condensation growth of om_ni - call cnst_get_ind('OM_AC' ,l_om_ac, abort=.true.) !om from coagulation processes - - call cnst_get_ind('DST_A2' ,l_dst_a2, abort=.true.) !Dust accumulation mode - call cnst_get_ind('DST_A3' ,l_dst_a3, abort=.true.) !Dust coarse mode - - call cnst_get_ind('SS_A1' ,l_ss_a1, abort=.true.) !Sea salt fine mode - call cnst_get_ind('SS_A2' ,l_ss_a2, abort=.true.) !Sea salt accumulation mode - call cnst_get_ind('SS_A3' ,l_ss_a3, abort=.true.) !Sea salt coarse mode - - ! register SOA species - call cnst_get_ind('SOA_NA' ,l_soa_na, abort=.true.) !Aitken mode SOA with SO4 and SOA condensate - call cnst_get_ind('SOA_A1' ,l_soa_a1, abort=.true.) !SOA condensate - call cnst_get_ind('SOA_LV' ,l_soa_lv, abort=.true.) !Gas phase low volatile SOA - call cnst_get_ind('SOA_SV' ,l_soa_sv, abort=.true.) !Gas phase semi volatile SOA - - ! gas phase h2so4 - call cnst_get_ind('H2SO4' ,l_h2so4, abort=.true.) - - ! Register the tracers in modes - call registerTracersInMode() - - ! Set the aerosol types - aerosolType(:) = -99 - aerosolType(l_so4_na) = AEROSOL_TYPE_SULFATE - aerosolType(l_so4_a1) = AEROSOL_TYPE_SULFATE - aerosolType(l_so4_a2) = AEROSOL_TYPE_SULFATE - aerosolType(l_so4_ac) = AEROSOL_TYPE_SULFATE - aerosolType(l_so4_pr) = AEROSOL_TYPE_SULFATE - aerosolType(l_bc_n) = AEROSOL_TYPE_BC - aerosolType(l_bc_ax) = AEROSOL_TYPE_BC - aerosolType(l_bc_ni) = AEROSOL_TYPE_BC - aerosolType(l_bc_a) = AEROSOL_TYPE_BC - aerosolType(l_bc_ai) = AEROSOL_TYPE_BC - aerosolType(l_bc_ac) = AEROSOL_TYPE_BC - aerosolType(l_om_ni) = AEROSOL_TYPE_OM - aerosolType(l_om_ai) = AEROSOL_TYPE_OM - aerosolType(l_om_ac) = AEROSOL_TYPE_OM - aerosolType(l_dst_a2) = AEROSOL_TYPE_DUST - aerosolType(l_dst_a3) = AEROSOL_TYPE_DUST - aerosolType(l_ss_a1) = AEROSOL_TYPE_SALT - aerosolType(l_ss_a2) = AEROSOL_TYPE_SALT - aerosolType(l_ss_a3) = AEROSOL_TYPE_SALT - aerosolType(l_soa_na) = AEROSOL_TYPE_OM - aerosolType(l_soa_a1) = AEROSOL_TYPE_OM - - rhopart(:)= 1000.0_r8 - - ! assign values based on aerosol type - do m=0,nmodes - do l=1,n_tracers_in_mode(m) - mm= getTracerIndex(m,l,.false.) - osmoticCoefficient(mm) = aerosol_type_osmotic_coefficient(aerosolType(mm)) - rhopart(mm) = aerosol_type_density(aerosolType(mm)) - solubleMassFraction(mm) = aerosol_type_soluble_mass_fraction(aerosolType(mm)) - numberOfIons(mm) = aerosol_type_number_of_ions(aerosolType(mm)) - end do - end do - - !SPECIAL CASES OF AEROSOL PROPERTIES: - !Density of bc_ax is rewritten later (calculated from fractal dimension) - !so4_a2 is different since it is ammonium sulfate and not sulf. acid. - rhopart(l_so4_a2) = 1769.0_r8 - - !These are not really particles, but set densities for the condenseable vapours - !used by condtend - rhopart(l_h2so4)= 1841.0_r8 - rhopart(l_soa_lv) = aerosol_type_density(AEROSOL_TYPE_OM) - rhopart(l_soa_sv) = aerosol_type_density(AEROSOL_TYPE_OM) - - ! Inverse calculated to avoid unneeded divisions in loop - invrhopart(:)=1._r8/rhopart(:) - - !Set process mode sizes - tracerInProcessMode = (/l_so4_a1, l_so4_a2, l_so4_ac, l_om_ac, l_bc_ac, l_soa_a1 /) - processModeMap(:)=-99 !Force error if using unset values - do l =1,pcnst - do m=1,numberOfProcessModeTracers - if(tracerInProcessMode(m) .eq. l)then - processModeMap(l)=m - end if - end do - end do - - ! Find out first mozart tracers (fxm: short lived species might mess up this!) - call cnst_get_ind(trim(solsym(1)), imozart, abort=.true.) - - !Add the cloud-tracers - isAlreadyCounted(:) = .false. - cloudTracerIndex(:) = -1 - do m=1,nmodes - do l=1,n_tracers_in_mode(m) - mm= getTracerIndex(m,l,.false.) - if(.not. isAlreadyCounted(mm))then - cloudTracerName(mm) = trim(cnst_name(mm))//"_OCW" - call pbuf_add_field(trim(cloudTracerName(mm)), 'global', dtype_r8, (/pcols,pver/), idx_dum) - ! Set the module variable qqcw(mm) to be set to idx_dum - call qqcw_set_ptr(mm,idx_dum) - cloudTracerIndex(mm) = idx_dum - isAlreadyCounted(mm) = .true. - endif - end do - end do - - !Find out how many aerosol-tracers we carry - isAlreadyCounted(:) = .false. - n_aerosol_tracers = 0 - do m=1,nmodes - do l=1,n_tracers_in_mode(m) - mm=getTracerIndex(m,l,.false.) - if(.not. isAlreadyCounted(mm))then - n_aerosol_tracers = n_aerosol_tracers + 1 - isAlreadyCounted(mm)=.true. - endif - end do - end do - - !Tabulated rh-growth for all species - call inittabrh - - end subroutine aero_register - - !============================================================================= - function getNumberOfAerosolTracers()RESULT(numberOfTracers) - integer :: numberOfTracers - numberOfTracers = n_aerosol_tracers - end function getNumberOfAerosolTracers - - !============================================================================= - function chemistryIndex(phys_index) RESULT (chemistryIndexOut) - integer, intent(in) :: phys_index - integer :: chemistryIndexOut - chemistryIndexOut = phys_index - imozart + 1 - end function chemistryIndex - - !============================================================================= - function physicsIndex(chem_index) RESULT(physIndexOut) - integer, intent(in) :: chem_index - integer :: physIndexOut - physIndexOut = chem_index + imozart - 1 - end function physicsIndex - - !============================================================================= - function isAerosol(phys_index) RESULT(answer) - integer, intent(in) :: phys_index - logical answer - answer=.FALSE. - if(aerosolType(phys_index) .gt. 0)then - answer = .TRUE. - endif - end function isAerosol - - !============================================================================= - function getNumberOfTracersInMode(modeIndex) RESULT(numberOfSpecies) - integer, intent(in) :: modeIndex - integer numberOfSpecies - numberOfSpecies = n_tracers_in_mode(modeIndex) - end function getNumberOfTracersInMode - - !============================================================================= - function getNumberOfBackgroundTracersInMode(modeIndex) RESULT (numberOfBackgroundSpecies) - integer, intent(in) :: modeIndex - integer numberOfBackgroundSpecies - numberOfBackgroundSpecies = n_background_tracers_in_mode(modeIndex) - end function getNumberOfBackgroundTracersInMode - - !============================================================================= - function getTracerIndex(modeIndex, componentIndex, isChemistry) RESULT(tracerIndex) - !purpose: Ask for an index in mode - !The index is the index in the q-array - !Some tracers may exist in several modes (is that a problem??) - integer, intent(in) :: modeIndex - integer, intent(in) :: componentIndex - logical, intent(in) :: isChemistry - integer tracerIndex - if(isChemistry)then - !This is tracer index in physics array - tracerIndex = tracer_in_mode(modeIndex,componentIndex)-imozart+1 - else - tracerIndex = tracer_in_mode(modeIndex,componentIndex) - endif - end function getTracerIndex - - !=============================================================================== - function getCloudTracerIndex(modeIndex, componentIndex) RESULT(cloud_tracer_index) - - ! Obtain an index in the physics-buffer for a component in the lifecycle scheme - - integer, intent(in) :: modeIndex - integer, intent(in) :: componentIndex - - integer :: tracerIndex - integer :: cloud_tracer_index - - if(componentIndex == 0)then - !Special key for number concentration of a mode - call endrun("error no such species") - else if (componentIndex > 0)then - !Lifecycle specie in a mode - tracerIndex = getTracerIndex(modeIndex,componentIndex,.false.) - cloud_tracer_index = cloudTracerIndex(tracerIndex) !ak: Index in phys-buffer - else - call endrun("negative componentindex in getCloudTracerIndex") - endif - end function getCloudTracerIndex - - !=============================================================================== - function getCloudTracerIndexDirect(tracerIndex) RESULT(cloudTracerIndexOut) - !returns index in pbuf for the corresponding cloud tracer with physics index "tracerIndex" - !returns "-1" if the tracer does not have any corresponding cloud tracer - integer, intent(in) :: tracerIndex - integer :: cloudTracerIndexOut - cloudTracerIndexOut = cloudTracerIndex(tracerIndex) - end function getCloudTracerIndexDirect - - !=============================================================================== - function getDryDensity(m,l) RESULT(density) - integer, intent(in) :: m !mode index - integer, intent(in) :: l !tracer index - real(r8) :: density - density = rhopart(tracer_in_mode(m,l)) - end function getDryDensity - - !=============================================================================== - function getCloudTracerName(tracerIndex) RESULT(cloudTracerNameOut) - integer, intent(in) :: tracerIndex - character(len=20) :: cloudTracerNameOut - cloudTracerNameOut = trim(cloudTracerName(tracerIndex)) - end function getCloudTracerName - - !=============================================================================== - subroutine fillAerosolTracerList(aerosolTracerList) - integer, dimension (:), intent(out) :: aerosolTracerList - logical, dimension(pcnst) :: alreadyFound - integer :: m,l,mm,nTracer - alreadyFound(:) = .FALSE. - nTracer = 0 - do m=1,nmodes - do l=1,n_tracers_in_mode(m) - mm=getTracerIndex(m,l,.FALSE.) - if(.NOT.alreadyFound(mm))then - nTracer = nTracer + 1 - alreadyFound(mm) = .TRUE. - aerosolTracerList(nTracer) = mm - end if - end do - end do - end subroutine fillAerosolTracerList - - !=============================================================================== - subroutine fillInverseAerosolTracerList(aerosolTracerList, inverseAerosolTracerList, n_aerosol_tracers) - integer, dimension(:), intent(in) :: aerosolTracerList - integer, intent(in) :: n_aerosol_tracers - integer, dimension(pcnst), intent(out) :: inverseAerosolTracerList - integer :: i - - inverseAerosolTracerList(:) = -99 - do i=1,n_aerosol_tracers - inverseAerosolTracerList(aerosolTracerList(i)) = i - end do - end subroutine fillInverseAerosolTracerList - - !=============================================================================== - subroutine registerTracersInMode() - !Register tracer index in modes - tracer_in_mode(:,:) = -1 !undefined - - !externally mixed bc - tracer_in_mode(MODE_IDX_BC_EXT_AC, 1:n_tracers_in_mode(MODE_IDX_BC_EXT_AC)) = (/l_bc_ax/) - - !sulphate + soa, sulfate condensate. - tracer_in_mode(MODE_IDX_SO4SOA_AIT, 1:n_tracers_in_mode(MODE_IDX_SO4SOA_AIT) ) = & - (/l_so4_na, l_soa_na, l_so4_a1, l_soa_a1/) - - !bc + sulfate condensate - tracer_in_mode(MODE_IDX_BC_AIT,1:n_tracers_in_mode(MODE_IDX_BC_AIT)) = & - (/l_bc_a, l_so4_a1, l_soa_a1/) - - !index not used - !tracer_in_mode(MODE_IDX_NOT_USED, 1:n_tracers_in_mode(MODE_IDX_NOT_USED)) = (/-1/) - - !om / bc internally mixed with sulfate condensate and aquous phase sulfate - tracer_in_mode(MODE_IDX_OMBC_INTMIX_COAT_AIT, 1:n_tracers_in_mode(MODE_IDX_OMBC_INTMIX_COAT_AIT))= & - (/l_bc_ai, l_om_ai, l_so4_a1, l_so4_a2, l_soa_a1 /) - - !accumulation mode sulfate with coagulate, condensate and aquous phase sulfate - tracer_in_mode(MODE_IDX_SO4_AC, 1:n_tracers_in_mode(MODE_IDX_SO4_AC)) = & - (/l_so4_pr, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !ac-mode dust with sulfate coagulate, condensate sulfate and wet-phase sulfate - tracer_in_mode(MODE_IDX_DST_A2, 1:n_tracers_in_mode(MODE_IDX_DST_A2)) = & - (/l_dst_a2, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !coarse mode dust with sulfate coagulate, condensate sulfate and wet-phase sulfate - tracer_in_mode(MODE_IDX_DST_A3, 1:n_tracers_in_mode(MODE_IDX_DST_A3)) = & - (/l_dst_a3, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !at-mode ss with sulfate coagulate, condensate sulfate and wet-phase sulfate - tracer_in_mode(MODE_IDX_SS_A1, 1:n_tracers_in_mode(MODE_IDX_SS_A1)) = & - (/l_ss_a1, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !ac mode ss with sulfate coagulate, condensate sulfate and wet-phase sulfate - tracer_in_mode(MODE_IDX_SS_A2, 1:n_tracers_in_mode(MODE_IDX_SS_A2)) = & - (/l_ss_a2, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !coarse mode ss sulfate coagulate, condensate sulfate and wet-phase sulfate - tracer_in_mode(MODE_IDX_SS_A3, 1:n_tracers_in_mode(MODE_IDX_SS_A3)) = & - (/l_ss_a3, l_bc_ac, l_om_ac, l_so4_a1, l_so4_ac, l_so4_a2, l_soa_a1 /) - - !sulfate + soa nucleation mode (mode no longer used) - !tracer_in_mode(MODE_IDX_SO4SOA_NUC, 1:n_tracers_in_mode(MODE_IDX_SO4SOA_NUC)) = (/ -1 /) - - !bc in nucleation mode - tracer_in_mode(MODE_IDX_BC_NUC, 1:n_tracers_in_mode(MODE_IDX_BC_NUC)) = (/l_bc_n/) - - !lumped organics - !tracer_in_mode(MODE_IDX_LUMPED_ORGANICS, 1:n_tracers_in_mode(MODE_IDX_LUMPED_ORGANICS)) = (/-1/) - - !intermal mixture bc/oc coated - tracer_in_mode(MODE_IDX_OMBC_INTMIX_AIT, 1:n_tracers_in_mode(MODE_IDX_OMBC_INTMIX_AIT)) = (/l_bc_ni, l_om_ni/) - end subroutine registerTracersInMode - - !=============================================================================== - function isTracerInMode(modeIndex, constituentIndex)RESULT(answer) - integer, intent(in) :: modeIndex - integer, intent(in) :: constituentIndex - integer :: i - logical :: answer - answer = .FALSE. - do i=1,n_tracers_in_mode(modeIndex) - if(tracer_in_mode(modeIndex,i) == constituentIndex)then - answer = .TRUE. - endif - enddo - end function isTracerInMode - - !=============================================================================== - function getConstituentFraction(CProcessModes, f_c, f_bc, f_aq, f_so4_cond, f_soa & - ,Cam, f_acm, f_bcm, f_aqm, f_so4_condm,f_soam, constituentIndex,debugPrint ) RESULT(fraction) ! mass fraction - - real(r8), intent(in) :: CProcessModes - real(r8), intent(in) :: f_c - real(r8), intent(in) :: f_bc - real(r8), intent(in) :: f_aq - real(r8), intent(in) :: f_so4_cond - real(r8), intent(in) :: f_soa - real(r8), intent(in) :: cam - real(r8), intent(in) :: f_aqm - real(r8), intent(in) :: f_bcm - real(r8), intent(in) :: f_acm - real(r8), intent(in) :: f_so4_condm - real(r8), intent(in) :: f_soam - integer, intent(in) :: constituentIndex - logical, optional, intent(in) :: debugPrint - - logical :: doPrint = .false. - real(r8) :: fraction - - if(present(debugPrint))then - if(debugPrint) then - doPrint=.true. - endif - endif - - fraction = 1.0_r8 ! fraction = 1 for all tracers, except special cases (process modes) below - - !This fraction is the mass of a certain tracer in a specific size-mode divided by the total - !mass of the same tracer for (i.e. summed up over) all size-modes. This total mass is what - !is transported in the model, in the life cycle scheme. The word size-mode is here used for a mode in the - !aerosol size-distribution, which is assumed to be log-normal prior to growth. - if((l_so4_a1 .eq. constituentIndex))then !so4 condensation - fraction= (cam & - *(1.0_r8-f_acm) & !sulfate fraction - *(1.0_r8-f_aqm) & !fraction not from aq phase - *(f_so4_condm) & !fraction being condensate - ) & - / & - (CProcessModes*(1.0_r8-f_c)*(1.0_r8-f_aq)*f_so4_cond+smallConcentration) !total so4 condensate - - if (doPrint) then - print*, " " - print*, "conc ==>", CProcessmodes, cam - print*, "modefrc ==>", f_acm, f_aqm, f_so4_condm - print*, "totfrc ==>", f_c, f_aq, f_so4_cond - print*, "fraction ==>", cam/(CProcessModes+smallConcentration)*100.0, fraction*100 , "%" - endif - - else if(l_so4_ac .eq. constituentIndex) then ! so4 coagulation - fraction = (cam & - * (1.0_r8 - f_acm) & !sulfate fraction - * (1.0_r8 - f_aqm) & !fraction not from aq phase - * (1.0_r8 - f_so4_condm) & !fraction not being condensate - ) & - / & - (CProcessModes*(1.0_r8-f_c)*(1.0_r8-f_aq)*(1.0_r8-f_so4_cond) & !total non-aq sulf - + smallConcentration) - - else if(l_so4_a2 .eq. constituentIndex) then !so4 wet phase - fraction = (cam & - *(1.0_r8-f_acm) & !sulfate fraction - *f_aqm) & !aq phase fraction of sulfate - / & - (CProcessModes*(1.0_r8-f_c)*(f_aq)+smallConcentration) - - else if(l_bc_ac .eq. constituentIndex)then !bc coagulated - fraction = (cam & - *f_acm & ! carbonaceous fraction - *f_bcm) & ! bc fraction of carbonaceous - / & - (CProcessModes*f_c*f_bc+smallConcentration) - - else if(l_om_ac .eq. constituentIndex ) then !oc coagulated - fraction = (cam & - *f_acm & ! carbonaceous fraction - *(1.0_r8-f_bcm) & ! oc fraction of carbonaceous - *(1.0_r8-f_soam))& ! oc fraction which is soa - / & - (CProcessModes*f_c*(1.0_r8-f_bc)*(1.0_r8-f_soa)+smallConcentration) - - else if (l_soa_a1 .eq. constituentIndex) then !SOA condensate - fraction = cam & - *f_acm & !carbonaceous fraction - *(1.0_r8 -f_bcm) & !om fraction - *(f_soam) & !fraction of OM is SOA - / & - (CProcessModes * f_c* (1.0_r8 -f_bc)*f_soa + smallConcentration) - end if - - if (fraction .gt. 1.0_r8)then - fraction = 1.0_r8 - endif - end function getConstituentFraction - - !=============================================================================== - subroutine inittabrh() - - ! Tables for hygroscopic growth - - integer :: i - real(r8) :: rr0ss(10),rr0so4(10),rr0bcoc(10) - - data rr0ss / 1.00_r8, 1.00_r8, 1.02_r8, 1.57_r8, 1.88_r8, 1.97_r8, 2.12_r8, 2.35_r8, 2.88_r8, 3.62_r8 / - data rr0so4 / 1.00_r8, 1.34_r8, 1.39_r8, 1.52_r8, 1.62_r8, 1.69_r8, 1.78_r8, 1.92_r8, 2.22_r8, 2.79_r8 / - data rr0bcoc / 1.00_r8, 1.02_r8, 1.03_r8, 1.12_r8, 1.17_r8, 1.20_r8, 1.25_r8, 1.31_r8, 1.46_r8, 1.71_r8 / - - rdivr0(:,:)=1._r8 - - do i=1,10 - rdivr0(i,l_so4_na)=rr0so4(i) - rdivr0(i,l_so4_a1)=rr0so4(i) - rdivr0(i,l_so4_a2)=rr0so4(i) - rdivr0(i,l_so4_ac)=rr0so4(i) - rdivr0(i,l_so4_pr)=rr0so4(i) - - rdivr0(i,l_bc_a)=rr0bcoc(i) - - rdivr0(i,l_bc_ni)=rr0bcoc(i) - rdivr0(i,l_bc_ai)=rr0bcoc(i) - rdivr0(i,l_bc_ac)=rr0bcoc(i) - - rdivr0(i,l_om_ni)=rr0bcoc(i) - rdivr0(i,l_om_ai)=rr0bcoc(i) - rdivr0(i,l_om_ac)=rr0bcoc(i) - - rdivr0(i,l_ss_a1)=rr0ss(i) - rdivr0(i,l_ss_a2)=rr0ss(i) - rdivr0(i,l_ss_a3)=rr0ss(i) - - rdivr0(i,l_soa_na)=rr0bcoc(i) - end do - end subroutine inittabrh - - !=============================================================================== - subroutine qqcw_set_ptr(index, iptr) - integer, intent(in) :: index, iptr - if(index>0 .and. index <= pcnst ) then - qqcw(index)=iptr - else - call endrun('qqcw_set_ptr: attempting to set qqcw pointer already defined') - end if - end subroutine qqcw_set_ptr - - !=============================================================================== - function qqcw_get_field(pbuf, index) - - type(physics_buffer_desc), pointer :: pbuf(:) - integer, intent(in) :: index - - real(r8), pointer :: qqcw_get_field(:,:) - - nullify(qqcw_get_field) - if (index>0 .and. index <= pcnst) then - if (qqcw(index)>0) then - call pbuf_get_field(pbuf, qqcw(index), qqcw_get_field) - endif - end if - end function qqcw_get_field - - !=============================================================================== - subroutine calculateNumberConcentration(ncol, q, rho_air, numberConcentration) - - ! arguments - integer , intent(in) :: ncol !number of columns used - real(r8) , intent(in) :: q(pcols,pver,pcnst) ![kg/kg] mass mixing ratios - real(r8) , intent(in) :: rho_air(pcols,pver) ![kg/m3] air density - real(r8) , intent(out) :: numberConcentration(pcols,pver,0:nmodes) ![#/m3] number concentration - - ! local variables - integer :: m, l, mm, k - - numberConcentration(:,:,:) = 0.0_r8 - do m = 0, nmodes - do l=1,getNumberOfBackgroundTracersInMode(m) - mm = getTracerIndex(m,l,.false.) - do k=1,pver - numberConcentration(:ncol,k,m) = numberConcentration(:ncol,k,m) & - + ( q(:ncol,k,mm) / getDryDensity(m,l)) !Volume of this tracer - end do - end do - end do - - ! until now, the variable "numberConcentration" actually contained "volume mixing ratio" - ! the next couple of lines fixes this! - do m= 0, nmodes - do k=1,pver - numberConcentration(:ncol,k,m) = numberConcentration(:ncol,k,m) * rho_air(:ncol,k) * volumeToNumber(m) - end do - end do - - end subroutine calculateNumberConcentration - - !=================================================== - subroutine calculateNumberMedianRadius(numberConcentration, volumeConcentration, lnSigma, & - numberMedianRadius, ncol) - - !Note the "nmodes" here - real(r8) , intent(in) :: numberConcentration(pcols,pver,0:nmodes) ![#/m3] number concentration - real(r8) , intent(in) :: volumeConcentration(pcols,pver,nmodes) ![kg/kg] mass mixing ratios - real(r8) , intent(in) :: lnSigma(pcols,pver,nmodes) ![kg/m3] air density - integer , intent(in) :: ncol !number of columns used - real(r8) , intent(out) :: numberMedianRadius(pcols,pver,nmodes) ![m] - - integer :: n,k - - do n=1,nmodes - do k=1,pver - where(volumeConcentration(:ncol,k,n) .gt. 1.e-20_r8) - numberMedianRadius(:ncol, k, n) = 0.5_r8 & !diameter ==> radius - * (volumeConcentration(:ncol,k,n) & !conversion formula - * 6.0_r8/pi/numberConcentration(:ncol,k,n) & - *DEXP(-4.5_r8*lnsigma(:ncol,k,n)*lnsigma(:ncol,k,n)))**aThird - elsewhere - numberMedianRadius(:ncol,k,n) = originalNumberMedianRadius(n) - end where - end do - end do - - end subroutine calculateNumberMedianRadius - - !=================================================== - function calculateEquivalentDensityOfFractalMode( & - emissionDensity, emissionRadius, fractalDimension, modeNumberMedianRadius, modeStandardDeviation) & - result (equivalentDensityOfFractal) - - ! output equivalent density of a fractal mode - - ! arguments - real(r8), intent(in) :: emissionDensity ![kg/m3] density at point of emission - real(r8), intent(in) :: emissionRadius ![kg/m3] radius at point of emission - real(r8), intent(in) :: fractalDimension ![kg/m3] fractal dimension of mode - real(r8), intent(in) :: modeNumberMedianRadius ![m] number median radius of mode - real(r8), intent(in) :: modeStandardDeviation ![m] standard deviation of mode - real(r8) :: equivalentDensityOfFractal ! Output - - ! local variables - real(r8) :: sumVolume - real(r8) :: sumMass - real(r8) :: dN, dNdLogR, dLogR - real(r8) :: densityBin - integer :: i - - sumVolume = 0.0_r8 - sumMass = 0.0_r8 - do i=1, nbinsTab - dLogR = log(rBinEdge(i+1)/rBinEdge(i)) - dNdLogR = calculatedNdLogR(rBinMidPoint(i), modeNumberMedianRadius, modeStandardDeviation) - - !Equivalent density (decreases with size since larger particles are long "hair like" threads..) - if (rBinMidPoint(i) < emissionRadius)then - densityBin = emissionDensity - else - densityBin = emissionDensity*(emissionRadius/rBinMidPoint(i))**(3.0 - fractalDimension) - endif - - !number concentration in this bin - dN = dNdLogR * dLogR - - !sum up volume and mass (factor of 4*pi/3 omitted since in both numerator and nominator) - sumVolume = sumVolume + dN * (rBinMidPoint(i)**3) - sumMass = sumMass + dN * densityBin * (rBinMidPoint(i)**3) - end do - - ! Equivalent density is mass by volume - equivalentDensityOfFractal = sumMass / sumVolume - - end function calculateEquivalentDensityOfFractalMode - - !=================================================== - function calculatedNdLogR(actualRadius, numberMedianRadius, sigma) result (dNdLogR) - - real(r8), intent(in) :: actualRadius - real(r8), intent(in) :: numberMedianRadius - real(r8), intent(in) :: sigma - - real(r8) :: logSigma - real(r8) :: dNdLogR - - logSigma = log(sigma) - - ! This is the formula for the lognormal distribution - dNdLogR = 1.0_r8/(sqrt(2.0_r8*pi)*log(sigma)) & - * DEXP(-0.5_r8*(log(actualRadius/numberMedianRadius))**2/(logSigma**2)) - - end function calculatedNdLogR - - !=================================================== - function calculateLognormalCDF(actualRadius, numberMedianRadius, sigma) result(CDF) - - !http://en.wikipedia.org/wiki/Log-normal_distribution#Cumulative_distribution_function - real(r8), intent(in) :: actualRadius - real(r8), intent(in) :: numberMedianRadius - real(r8), intent(in) :: sigma - real(r8) :: CDF ! output - - real(r8) :: argument - - argument = -1.0_r8*(log(actualRadius/numberMedianRadius) / log(sigma) / sqrt(2.0_r8)) - CDF = 0.5_r8 * erfc(argument) - - end function calculateLognormalCDF - - !=================================================== - subroutine init_interp_constants() - - !--------------------------------------------------------------- - ! set module variables for interpolation of tables - !--------------------------------------------------------------- - - ! Local variables - integer :: irf, irelh, kcomp, i - !----------------------------------------------------------- - - ! Defining array bounds for tabulated optical parameters (and r and sigma) - ! relative humidity (only 0 value used for r and sigma tables): - rh = (/ 0.0_r8, 0.37_r8, 0.47_r8, 0.65_r8, 0.75_r8, 0.8_r8, 0.85_r8, 0.9_r8, 0.95_r8, 0.995_r8 /) - - ! relative humidity (RH, as integer for output variable names) for use in AeroCom code - RF = (/0, 40, 55, 65, 75, 85 /) - - ! AeroCom specific RH input variables for use in opticsAtConstRh.F90 - do irf=1,6 - xrhrf(irf) = real(RF(irf))*0.01_r8 - enddo - do irelh=1,9 - do irf=1,6 - if(xrhrf(irf)>=rh(irelh).and.xrhrf(irf)<=rh(irelh+1)) then - irhrf1(irf)=irelh - endif - end do - end do - - ! mass fractions internal mixtures in background (fombg and fbcbg) and mass added to the - ! background modes (fac, faq, faq) - fombg = (/ 0.0_r8, 0.2_r8, 0.4_r8, 0.6_r8, 0.8_r8, 1.0_r8 /) - fac = (/ 0.0_r8, 0.2_r8, 0.4_r8, 0.6_r8, 0.8_r8, 1.0_r8 /) - faq = (/ 0.0_r8, 0.2_r8, 0.4_r8, 0.6_r8, 0.8_r8, 1.0_r8 /) - - ! with more weight on low fractions (thus a logaritmic f axis) for BC, - ! which is less ambundant than sulfate and OC, and the first value - ! corresponding to a clean background mode: - ! and most weight on small concentrations for added mass onto the background: - - fbcbg(1)=1.e-10_r8 - fbc(1)=1.e-10_r8 - do i=2,6 - fbcbg(i)=10**((i-1)/4.0_r8-1.25_r8) - fbc(i)=fbcbg(i) - end do - - do kcomp=1,4 - cate(kcomp,1)=1.e-10_r8 - do i=2,16 - if(kcomp.eq.1.or.kcomp.eq.2) then - cate(kcomp,i)=10.0_r8**((i-1)/3.0_r8-6.222_r8) - elseif(kcomp.eq.3) then - cate(kcomp,i)=1.0e-10_r8 ! not used - else - cate(kcomp,i)=10.0_r8**((i-1)/3.0_r8-4.301_r8) - endif - end do - end do - do kcomp=5,10 - cat(kcomp,1) =1.e-10_r8 - do i=2,6 - if(kcomp.eq.5) then - cat(kcomp,i)=10.0_r8**((i-1)-3.824_r8) - elseif(kcomp.eq.6) then - cat(kcomp,i)=10.0_r8**((i-1)-3.523_r8) - elseif(kcomp.eq.7) then - cat(kcomp,i)=10.0_r8**((i-1)-3.699_r8) - elseif(kcomp.eq.8) then - cat(kcomp,i)=10.0_r8**((i-1)-4.921_r8) - elseif(kcomp.eq.9) then - cat(kcomp,i)=10.0_r8**((i-1)-3.301_r8) - else - cat(kcomp,i)=10.0_r8**((i-1)-3.699_r8) - endif - end do - end do - - end subroutine init_interp_constants - -end module oslo_aero_share diff --git a/src/aero_model.F90 b/src/aero_model.F90 index a252e8c..ca44150 100644 --- a/src/aero_model.F90 +++ b/src/aero_model.F90 @@ -43,8 +43,6 @@ 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_share, only: getNumberofTracersInMode, getCloudTracerIndexDirect, getCloudTracerName - use oslo_aero_share, only: getCloudTracerName, getTracerIndex, aero_register use oslo_aero_control, only: oslo_aero_ctl_readnl 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 @@ -54,6 +52,8 @@ module aero_model use oslo_aero_seasalt, only: oslo_aero_seasalt_init, oslo_aero_seasalt_emis, seasalt_active use oslo_aero_dust, only: oslo_aero_dust_init, oslo_aero_dust_emis, dust_active use oslo_aero_ocean, only: oslo_aero_ocean_init, oslo_aero_dms_emis + use oslo_aero_share, only: getNumberofTracersInMode, getCloudTracerIndexDirect, getCloudTracerName + use oslo_aero_share, only: getCloudTracerName, getTracerIndex, aero_register use oslo_aero_sox_cldaero, only: sox_cldaero_init use oslo_aero_microp, only: oslo_aero_microp_readnl use oslo_aero_sw_tables, only: initopt @@ -213,24 +213,24 @@ subroutine aero_model_init( pbuf2d ) endif enddo - call addfld ('NUCLRATE' ,(/'lev'/), 'A', '#/cm3/s','Nucleation rate') - call addfld ('FORMRATE' ,(/'lev'/), 'A', '#/cm3/s','Formation rate of 12nm particles') - call addfld ('COAGNUCL' ,(/'lev'/), 'A', '/s' ,'Coagulation sink for nucleating particles') - call addfld ('GRH2SO4' ,(/'lev'/), 'A', 'nm/hour','Growth rate H2SO4') - call addfld ('GRSOA' ,(/'lev'/), 'A', 'nm/hour','Growth rate SOA') - call addfld ('GR' ,(/'lev'/), 'A', 'nm/hour','Growth rate, H2SO4+SOA') - call addfld ('NUCLSOA' ,(/'lev'/), 'A', 'kg/kg' ,'SOA nucleate') - call addfld ('ORGNUCL' ,(/'lev'/), 'A', 'kg/kg' ,'Organic gas available for nucleation') + call addfld ('NUCLRATE',(/'lev'/), 'A','#/cm3/s','Nucleation rate') + call addfld ('FORMRATE',(/'lev'/), 'A','#/cm3/s','Formation rate of 12nm particles') + call addfld ('COAGNUCL',(/'lev'/), 'A', '/s','Coagulation sink for nucleating particles') + call addfld ('GRH2SO4',(/'lev'/), 'A', 'nm/hour','Growth rate H2SO4') + call addfld ('GRSOA',(/'lev'/),'A','nm/hour','Growth rate SOA') + call addfld ('GR',(/'lev'/), 'A', 'nm/hour','Growth rate, H2SO4+SOA') + call addfld ('NUCLSOA',(/'lev'/),'A','kg/kg','SOA nucleate') + call addfld ('ORGNUCL',(/'lev'/),'A','kg/kg','Organic gas available for nucleation') if(history_aerosol)then call add_default ('NUCLRATE', 1, ' ') call add_default ('FORMRATE', 1, ' ') call add_default ('COAGNUCL', 1, ' ') - call add_default ('GRH2SO4', 1, ' ') - call add_default ('GRSOA', 1, ' ') - call add_default ('GR', 1, ' ') - call add_default ('NUCLSOA', 1, ' ') - call add_default ('ORGNUCL', 1, ' ') + call add_default ('GRH2SO4', 1, ' ') + call add_default ('GRSOA', 1, ' ') + call add_default ('GR', 1, ' ') + call add_default ('NUCLSOA', 1, ' ') + call add_default ('ORGNUCL', 1, ' ') end if call addfld( 'XPH_LWC', (/ 'lev' /), 'A','kg/kg', 'pH value multiplied by lwc') @@ -238,9 +238,9 @@ subroutine aero_model_init( pbuf2d ) call addfld ('AQSO4_O3', horiz_only, 'A','kg/m2/s', 'SO4 aqueous phase chemistry due to O3') if ( history_aerosol ) then - call add_default ('XPH_LWC', 1, ' ') + call add_default ('XPH_LWC', 1, ' ') call add_default ('AQSO4_H2O2', 1, ' ') - call add_default ('AQSO4_O3', 1, ' ') + call add_default ('AQSO4_O3', 1, ' ') endif end subroutine aero_model_init @@ -668,8 +668,8 @@ end subroutine vmr2qqcw !============================================================================= subroutine aero_model_constants() - ! - ! A number of constants used in the emission and size-calculation in CAM-Oslo Jan 2011. + + ! A number of constants used in the emission and size-calculation in CAM-Oslo ! local variables integer :: kcomp,i diff --git a/src/oslo_aero_coag.F90 b/src/oslo_aero_coag.F90 index 2e8c6b0..b943a0b 100644 --- a/src/oslo_aero_coag.F90 +++ b/src/oslo_aero_coag.F90 @@ -99,7 +99,7 @@ module oslo_aero_coag real(r8), parameter :: temperatureLookupTables = 293.15_r8 !Temperature used in look up tables real(r8), parameter :: mfpAir = 63.3e-9_r8 ![m] mean free path air real(r8), parameter :: viscosityAir = 1.983e-5_r8 ![Pa s] viscosity of air - real(r8), parameter :: rhoh2o = 1000._r8 ! Density of water + real(r8), parameter :: rhoh2o = 1000._r8 !Density of water !================================================================ contains @@ -202,14 +202,14 @@ subroutine initializeCoagulation(rhob,rk) call calculateCoagulationCoefficient(CoagulationCoefficient & !O [m3/s] coagulation coefficient , rk(modeIndexCoagulator) & !I [m] radius of coagulator , rhob(modeIndexCoagulator) & !I [kg/m3] density of coagulator - , rhob(modeIndexReceiver) ) !I [kg/m3] density of receiver + , rhob(modeIndexReceiver) ) !I [kg/m3] density of receiver !Save values CoagCoeffModeAdd(iReceiverMode,:) = CoagulationCoefficient(:) end do !receiver modes - ! Onl one receivermode for cloud coagulation (water) + ! Only one receivermode for cloud coagulation (water) do iCoagulatingMode = 1,numberOfCoagulatingModes !Index of the coagulating mode (0-14), see list above @@ -221,7 +221,7 @@ subroutine initializeCoagulation(rhob,rk) call calculateCoagulationCoefficient(CoagulationCoefficient & !O [m3/s] coagulation coefficient , rk(modeIndexCoagulator) & !I [m] radius of coagulator , rhob(modeIndexCoagulator) & !I [kg/m3] density of coagulator - , rhoh2o ) !I [kg/m3] density of receiver + , rhoh2o ) !I [kg/m3] density of receiver !Save values K12Cl(iCoagulatingMode,:) = CoagulationCoefficient(:) @@ -361,8 +361,8 @@ subroutine calculateCoagulationCoefficient(CoagulationCoefficient, modeRadius, m real(r8) :: mfv1 ![m] mean free path particle real(r8) :: mfv2 ![m] mean free path particle - ! coagulation coefficient for SO4 (Brownian, Fuchs form) - !Loop through indexes in look-up table + ! coagulation coefficient for SO4 (Brownian, Fuchs form) + ! Loop through indexes in look-up table do i=1,nBinsTab c1=calculateThermalVelocity(rBinMidPoint(i), receiverDensity) !receiving size c2=calculateThermalVelocity(modeRadius, modeDensity) !coagulating aerosol @@ -379,14 +379,13 @@ subroutine calculateCoagulationCoefficient(CoagulationCoefficient, modeRadius, m g12=sqrt(g1**2+g2**2) - !Coagulation coefficient of receiver size "i" with the coagulating - !mode "kcomp" + !Coagulation coefficient of receiver size "i" with the coagulating mode "kcomp" CoagulationCoefficient(i) = & 4.0_r8*pi*(rBinMidPoint(i)+modeRadius)*(diff1+diff2) & /((rBinMidPoint(i)+modeRadius)/(rBinMidPoint(i)+modeRadius+g12) & +(4.0_r8/c12)*(diff1+diff2)/(modeRadius+rBinMidPoint(i))) - enddo ! loop on imax + enddo end subroutine calculateCoagulationCoefficient !================================================================ @@ -495,9 +494,9 @@ subroutine coagtend(q, pmid, pdel, temperature, delt_inverse, ncol, lchnk) !process modes don't change mode except so4 condensate which becomes coagulate instead !assumed to have same sink as MODE_IDX_OMBC_INTMIX_AIT - if( .NOT. is_process_mode(l_index_donor,.true.) & - .OR. ( (l_index_donor.eq.chemistryIndex(l_so4_a1)) & - .AND. modeIndexCoagulator .eq. MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then + if( .NOT. is_process_mode(l_index_donor,.true.) .or. & + ( (l_index_donor == chemistryIndex(l_so4_a1)) .and. & + modeIndexCoagulator == MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then !Done summing total loss of this coagulating specie totalLoss(i,k,l_index_donor) = coagulationSink & !loss rate for a mode in [1/s] summed over all receivers @@ -515,16 +514,16 @@ subroutine coagtend(q, pmid, pdel, temperature, delt_inverse, ncol, lchnk) end do ! k - !UPDATE THE TRACERS AND DO DIAGNOSTICS + ! UPDATE THE TRACERS AND DO DIAGNOSTICS do iCoagulator = 1, numberOfCoagulatingModes do ispecie = 1, getNumberOfTracersInMode(coagulatingMode(iCoagulator)) l_index_donor = getTracerIndex(coagulatingMode(iCoagulator) , ispecie ,.true.) !so4_a1 is a process mode (condensate), but is still lost in coagulation - if( .NOT. is_process_mode(l_index_donor, .true.) & - .OR. ( (l_index_donor.eq.chemistryIndex(l_so4_a1)) .AND. & - coagulatingMode(iCoagulator) .eq. MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then + if( .NOT. is_process_mode(l_index_donor, .true.) .or. & + ( (l_index_donor == chemistryIndex(l_so4_a1)) .and. & + coagulatingMode(iCoagulator) == MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then l_index_donor = getTracerIndex(coagulatingMode(iCoagulator) , ispecie,.true. ) @@ -639,13 +638,14 @@ subroutine clcoag(q, pmid, pdel, temperature, cldnum, cldfrc, delt_inverse, ncol !process modes don't change mode except so4 condensate which becomes coagulate instead !assumed to have same sink as MODE_IDX_OMBC_INTMIX_AIT - if( .NOT. is_process_mode(l_index_donor,.true.) & - .OR. ( (l_index_donor.eq.chemistryIndex(l_so4_a1)) .AND. modeIndexCoagulator .eq. MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then + if( .NOT. is_process_mode(l_index_donor,.true.) .or. & + ( (l_index_donor == chemistryIndex(l_so4_a1)) .and. & + modeIndexCoagulator == MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then !Done summing total loss of this coagulating specie - cloudLoss(i,k,l_index_donor) = coagulationSink & !loss rate for a mode in [1/s] summed over all receivers - * cldfrc(i,k)*q(i,k,l_index_donor) & !* mixing ratio ==> MMR/s - / delt_inverse ! seconds ==> MMR + cloudLoss(i,k,l_index_donor) = coagulationSink & !loss rate for a mode in [1/s] summed over all receivers + * cldfrc(i,k)*q(i,k,l_index_donor) & !* mixing ratio ==> MMR/s + / delt_inverse !/ seconds ==> MMR !Can not loose more than we have ! At present day assumed lost within the cloud @@ -659,15 +659,15 @@ subroutine clcoag(q, pmid, pdel, temperature, cldnum, cldfrc, delt_inverse, ncol end do ! i end do ! k - !UPDATE THE TRACERS AND DO DIAGNOSTICS + ! UPDATE THE TRACERS AND DO DIAGNOSTICS do iCoagulator = 1, numberOfCoagulatingModes do ispecie = 1, getNumberOfTracersInMode(coagulatingMode(iCoagulator)) l_index_donor = getTracerIndex(coagulatingMode(iCoagulator) , ispecie ,.true.) !so4_a1 is a process mode (condensate), but is still lost in coagulation - if( .NOT. is_process_mode(l_index_donor, .true.) & - .OR. ( (l_index_donor.eq.chemistryIndex(l_so4_a1)) .AND. coagulatingMode(iCoagulator) & - .eq. MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then + if( .NOT. is_process_mode(l_index_donor, .true.) .or. & + ( (l_index_donor == chemistryIndex(l_so4_a1)) .and. & + coagulatingMode(iCoagulator) == MODE_IDX_OMBC_INTMIX_COAT_AIT) ) then l_index_donor = getTracerIndex(coagulatingMode(iCoagulator), ispecie, .true.) diff --git a/src/oslo_aero_nucleate_ice.F90 b/src/oslo_aero_nucleate_ice.F90 index e5eba9d..0640b5c 100644 --- a/src/oslo_aero_nucleate_ice.F90 +++ b/src/oslo_aero_nucleate_ice.F90 @@ -39,7 +39,6 @@ module oslo_aero_nucleate_ice use tropopause, only: tropopause_findChemTrop use cam_logfile, only: iulog use cam_abortutils, only: endrun - use time_manager, only: is_first_step use nucleate_ice, only: nucleati_init, nucleati ! portable module ! use oslo_aero_share, only: l_dst_a2, l_dst_a3, MODE_IDX_DST_A2, MODE_IDX_DST_A3, rhopart, qqcw_get_field @@ -164,13 +163,6 @@ subroutine nucleate_ice_oslo_init(mincld_in, bulk_scale_in) mincld = mincld_in bulk_scale = bulk_scale_in - ! TODO: IS this necessary? - ! if (is_first_step()) then - ! call pbuf_set_field(pbuf, naai_idx, 0.0_r8) - ! end if - - ! Initialize naai. - if( masterproc ) then write(iulog,*) 'nucleate_ice parameters:' write(iulog,*) ' mincld = ', mincld_in @@ -299,16 +291,14 @@ subroutine nucleate_ice_oslo_calc( state, wsubi, pbuf, dtime, ptend, numberConce real(r8), pointer :: qi(:,:) ! cloud ice mixing ratio (kg/kg) real(r8), pointer :: ni(:,:) ! cloud ice number conc (1/kg) real(r8), pointer :: pmid(:,:) ! pressure at layer midpoints (pa) - + real(r8), pointer :: cld_dst_a2(:,:) ! mmr cld dst a2 + real(r8), pointer :: cld_dst_a3(:,:) ! mass m.r. of coarse dust real(r8), pointer :: ast(:,:) - real(r8) :: icecldf(pcols,pver) ! ice cloud fraction real(r8), pointer :: qsatfac(:,:) ! Subgrid cloud water saturation scaling factor. - real(r8), pointer :: cld_dst_a2(:,:) ! mmr cld dst a2 - real(r8), pointer :: cld_dst_a3(:,:) ! mass m.r. of coarse dust + real(r8) :: icecldf(pcols,pver) ! ice cloud fraction real(r8) :: rho(pcols,pver) ! air density (kg m-3) - real(r8) :: qs(pcols) ! liquid-ice weighted sat mixing rat (kg/kg) real(r8) :: es(pcols) ! liquid-ice weighted sat vapor press (pa) real(r8) :: gammas(pcols) ! parameter for cond/evap of cloud water diff --git a/src/oslo_aero_ocean.F90 b/src/oslo_aero_ocean.F90 index dd99fbb..07410fb 100644 --- a/src/oslo_aero_ocean.F90 +++ b/src/oslo_aero_ocean.F90 @@ -40,7 +40,7 @@ module oslo_aero_ocean ! Public interfaces public :: oslo_aero_ocean_init ! initializing, reading file - public :: oslo_aero_ocean_time ! time interpolation + public :: oslo_aero_ocean_adv ! spatial and time interpolation public :: oslo_aero_dms_emis ! calculate dms surface emissions public :: oslo_aero_dms_inq ! logical function which tells mo_srf_emis what to do public :: oslo_aero_opom_emis ! calculate opom surface emissions @@ -58,7 +58,7 @@ module oslo_aero_ocean integer :: dms_cycle_yr = 0 !will be collected from NAMELIST !will be collected from NAMELIST - but can be overwritten by atm_import_export if the ocean is sending DMS to the atm - character(len=20), public :: dms_source = 'emission_file' + character(len=20), public :: dms_source = 'emission_file' ! character(len=16) :: opomo_fld_name = 'chlor_a' !not set from namelist, hard coded, name of nc var character(len=16) :: opomn_fld_name = 'poc' !not set from namelist, hard coded, name of nc var @@ -189,7 +189,7 @@ subroutine oslo_aero_ocean_init() endsubroutine oslo_aero_ocean_init !=============================================================================== - subroutine oslo_aero_ocean_time(state, pbuf2d) + subroutine oslo_aero_ocean_adv(state, pbuf2d) ! Interpolate ocean_species in space and time to model grid and time @@ -204,7 +204,7 @@ subroutine oslo_aero_ocean_time(state, pbuf2d) call advance_trcdata( oceanspcs(m)%fields, oceanspcs(m)%file, state, pbuf2d ) end do - endsubroutine oslo_aero_ocean_time + endsubroutine oslo_aero_ocean_adv !=============================================================================== subroutine oslo_aero_dms_emis(ncol, lchnk, u, v, zm, ocnfrc, icefrc, sst, fdms, cflx) @@ -219,7 +219,7 @@ subroutine oslo_aero_dms_emis(ncol, lchnk, u, v, zm, ocnfrc, icefrc, sst, fdms, real(r8) , intent(in) :: icefrc(pcols) real(r8) , intent(in) :: sst(pcols) real(r8) , intent(in) :: fdms(pcols) - real(r8) , intent(inout) :: cflx(pcols,pcnst) + real(r8) , intent(inout) :: cflx(pcols,pcnst) ! local variables real(r8) :: u10m(pcols) ! [m/s] diff --git a/src/oslo_aero_optical_params.F90 b/src/oslo_aero_optical_params.F90 index 1d9e4ed..b49ddf1 100644 --- a/src/oslo_aero_optical_params.F90 +++ b/src/oslo_aero_optical_params.F90 @@ -31,7 +31,7 @@ module oslo_aero_optical_params !=============================================================================== subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, & - coszrs, state, t, cld, qm1, Nnatk, & + coszrs, state, t, cld, qm1, & per_tau, per_tau_w, per_tau_w_g, per_tau_w_f, per_lw_abs, & volc_ext_sun, volc_omega_sun, volc_g_sun, volc_ext_earth, volc_omega_earth, & aodvis, absvis) @@ -52,9 +52,6 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, real(r8), intent(in) :: volc_omega_earth(pcols,pver,nlwbands) ! volcanic aerosol SSA for terrestrial bands, CMIP6 type(physics_state), intent(in), target :: state - ! Input-output arguments - real(r8), intent(inout) :: Nnatk(pcols,pver,0:nmodes) ! aerosol mode number concentration - ! Output arguments ! AOD and absorptive AOD for visible wavelength closest to 0.55 um (0.442-0.625) ! Note that aodvis and absvis output should be divided by dayfoc to give physical (A)AOD values @@ -69,13 +66,14 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, ! Local variables integer :: i, k, ib, icol, mplus10 integer :: iloop - logical :: daylight(pcols) ! SW calculations also at (polar) night in interpol* if daylight=.true. - real(r8) :: aodvisvolc(pcols) ! AOD vis for CMIP6 volcanic aerosol - real(r8) :: absvisvolc(pcols) ! AAOD vis for CMIP6 volcanic aerosol - real(r8) :: bevisvolc(pcols,pver) ! Extinction in vis wavelength band for CMIP6 volcanic aerosol - real(r8) :: rhum(pcols,pver) ! (trimmed) relative humidity for the aerosol calculations - real(r8) :: deltah_km(pcols,pver) ! Layer thickness, unit km - real(r8) :: deltah, airmassl(pcols,pver), airmass(pcols) !akc6 + real(r8) :: Nnatk(pcols,pver,0:nmodes) ! aerosol mode number concentration + logical :: daylight(pcols) ! SW calculations also at (polar) night in interpol* if daylight=.true. + real(r8) :: aodvisvolc(pcols) ! AOD vis for CMIP6 volcanic aerosol + real(r8) :: absvisvolc(pcols) ! AAOD vis for CMIP6 volcanic aerosol + real(r8) :: bevisvolc(pcols,pver) ! Extinction in vis wavelength band for CMIP6 volcanic aerosol + real(r8) :: rhum(pcols,pver) ! (trimmed) relative humidity for the aerosol calculations + real(r8) :: deltah_km(pcols,pver) ! Layer thickness, unit km + real(r8) :: deltah, airmassl(pcols,pver), airmass(pcols) real(r8) :: Ca(pcols,pver), f_c(pcols,pver), f_bc(pcols,pver), f_aq(pcols,pver) real(r8) :: fnbc(pcols,pver), faitbc(pcols,pver), f_so4_cond(pcols,pver) real(r8) :: f_soa(pcols,pver),f_soana(pcols,pver) @@ -419,7 +417,6 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, per_tau_w_f(i,k,ib)=per_tau_w_g(i,k,ib)*asymtot(i,k,ib) end do end do ! ncol - !------------------------------------------------------------------------------------------------ ! LW Optical properties of total aerosol: do ib=1,nlwbands @@ -521,7 +518,7 @@ subroutine oslo_aero_optical_params_calc(lchnk, ncol, pint, pmid, ! Layer thickness, unit km, and layer airmass, unit kg/m2 deltah=deltah_km(icol,k) airmassl(icol,k)=1.e3_r8*deltah*rhoda(icol,k) - airmass(icol)=airmass(icol)+airmassl(icol,k) !akc6 + airmass(icol)=airmass(icol)+airmassl(icol,k) ! Optical depths at ca. 550 nm (0.442-0.625um) all aerosols aodvis(icol)=aodvis(icol)+betotvis(icol,k)*deltah diff --git a/src_cam/radiation.F90 b/src_cam/radiation.F90 index 642e924..9f5bd33 100644 --- a/src_cam/radiation.F90 +++ b/src_cam/radiation.F90 @@ -1,1867 +1,1728 @@ module radiation - !--------------------------------------------------------------------------------- - ! - ! CAM interface to RRTMG radiation parameterization - ! - !--------------------------------------------------------------------------------- - - use shr_kind_mod, only: r8=>shr_kind_r8 - use spmd_utils, only: masterproc - use ppgrid, only: pcols, pver, pverp, begchunk, endchunk - use physics_types, only: physics_state, physics_ptend - use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx - use camsrfexch, only: cam_out_t, cam_in_t - use physconst, only: cappa, cpair - - use time_manager, only: get_nstep, is_first_restart_step, & - get_curr_calday, get_step_size - - use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_info, & - rad_cnst_get_gas, rad_cnst_out, oldcldoptics, & - liqcldoptics, icecldoptics - - use radconstants, only: nswbands, nlwbands, rrtmg_sw_cloudsim_band, rrtmg_lw_cloudsim_band, & - idx_sw_diag - - use cospsimulator_intr, only: docosp, cospsimulator_intr_init, & - cospsimulator_intr_run, cosp_nradsteps - - use scamMod, only: scm_crm_mode, single_column, have_cld, cldobs - - use cam_history, only: addfld, add_default, horiz_only, outfld, hist_fld_active - use cam_history_support, only: fillvalue - - use pio, only: file_desc_t, var_desc_t, & - pio_int, pio_double, pio_noerr, & - pio_seterrorhandling, pio_bcast_error, & - pio_inq_varid, pio_def_var, & - pio_put_var, pio_get_var, pio_put_att - - use cam_abortutils, only: endrun - use error_messages, only: handle_err - use perf_mod, only: t_startf, t_stopf - use cam_logfile, only: iulog -#ifdef OSLO_AERO - 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 +!--------------------------------------------------------------------------------- +! +! CAM interface to RRTMG radiation parameterization +! +!--------------------------------------------------------------------------------- + +use shr_kind_mod, only: r8=>shr_kind_r8 +use spmd_utils, only: masterproc +use ppgrid, only: pcols, pver, pverp, begchunk, endchunk +use physics_types, only: physics_state, physics_ptend +use physics_buffer, only: physics_buffer_desc, pbuf_get_field, pbuf_old_tim_idx +use camsrfexch, only: cam_out_t, cam_in_t +use physconst, only: cappa, cpair + +use time_manager, only: get_nstep, is_first_restart_step, & + get_curr_calday, get_step_size + +use rad_constituents, only: N_DIAG, rad_cnst_get_call_list, rad_cnst_get_info, & + rad_cnst_get_gas, rad_cnst_out, oldcldoptics, & + liqcldoptics, icecldoptics + +use radconstants, only: nswbands, nlwbands, rrtmg_sw_cloudsim_band, rrtmg_lw_cloudsim_band, & + idx_sw_diag + +use cospsimulator_intr, only: docosp, cospsimulator_intr_init, & + cospsimulator_intr_run, cosp_nradsteps + +use scamMod, only: scm_crm_mode, single_column, have_cld, cldobs + +use cam_history, only: addfld, add_default, horiz_only, outfld, hist_fld_active +use cam_history_support, only: fillvalue + +use pio, only: file_desc_t, var_desc_t, & + pio_int, pio_double, pio_noerr, & + pio_seterrorhandling, pio_bcast_error, & + pio_inq_varid, pio_def_var, & + pio_put_var, pio_get_var, pio_put_att + +use cam_abortutils, only: endrun +use error_messages, only: handle_err +use perf_mod, only: t_startf, t_stopf +use cam_logfile, only: iulog +! OSLO_AERO begin +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_aerocom, only: dod440, dod550, dod870, abs550, abs550alt -#endif +use oslo_aero_aerocom, only: dod440, dod550, dod870, abs550, abs550alt #endif - - implicit none - private - - public :: & - radiation_readnl, &! read namelist variables - radiation_register, &! registers radiation physics buffer fields - radiation_nextsw_cday, &! calendar day of next radiation calculation - radiation_do, &! query which radiation calcs are done this timestep - radiation_init, &! initialization - radiation_define_restart, &! define variables for restart - radiation_write_restart, &! write variables to restart - radiation_read_restart, &! read variables from restart - radiation_tend, &! compute heating rates and fluxes - rad_out_t ! type for diagnostic outputs - - integer,public, allocatable :: cosp_cnt(:) ! counter for cosp - integer,public :: cosp_cnt_init = 0 !initial value for cosp counter - real(r8), public, protected :: nextsw_cday ! future radiation calday for surface models - - type rad_out_t - real(r8) :: solin(pcols) ! Solar incident flux - - real(r8) :: qrsc(pcols,pver) - - real(r8) :: fsntc(pcols) ! Clear sky total column abs solar flux - real(r8) :: fsntoa(pcols) ! Net solar flux at TOA - real(r8) :: fsntoac(pcols) ! Clear sky net solar flux at TOA - real(r8) :: fsutoa(pcols) ! upwelling solar flux at TOA - - real(r8) :: fsnirt(pcols) ! Near-IR flux absorbed at toa - real(r8) :: fsnrtc(pcols) ! Clear sky near-IR flux absorbed at toa - real(r8) :: fsnirtsq(pcols) ! Near-IR flux absorbed at toa >= 0.7 microns - - real(r8) :: fsn200(pcols) ! fns interpolated to 200 mb - real(r8) :: fsn200c(pcols) ! fcns interpolated to 200 mb - real(r8) :: fsnr(pcols) ! fns interpolated to tropopause - - real(r8) :: fsnsc(pcols) ! Clear sky surface abs solar flux - real(r8) :: fsdsc(pcols) ! Clear sky surface downwelling solar flux - - real(r8) :: qrlc(pcols,pver) - - real(r8) :: flntc(pcols) ! Clear sky lw flux at model top - real(r8) :: flut(pcols) ! Upward flux at top of model - real(r8) :: flutc(pcols) ! Upward Clear Sky flux at top of model - real(r8) :: lwcf(pcols) ! longwave cloud forcing - - real(r8) :: fln200(pcols) ! net longwave flux interpolated to 200 mb - real(r8) :: fln200c(pcols) ! net clearsky longwave flux interpolated to 200 mb - real(r8) :: flnr(pcols) ! net longwave flux interpolated to tropopause - - real(r8) :: flnsc(pcols) ! Clear sky lw flux at srf (up-down) - real(r8) :: fldsc(pcols) ! Clear sky lw flux at srf (down) - - real(r8) :: tot_cld_vistau(pcols,pver) ! gbx water+ice cloud optical depth (only during day, night = fillvalue) - real(r8) :: tot_icld_vistau(pcols,pver) ! in-cld water+ice cloud optical depth (only during day, night = fillvalue) - real(r8) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth (only during day, night = fillvalue) - real(r8) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth (only during day, night = fillvalue) - real(r8) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth for output on history files - real(r8) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth for output on history files - - real(r8) :: cld_tau_cloudsim(pcols,pver) - real(r8) :: aer_tau400(pcols,0:pver) - real(r8) :: aer_tau550(pcols,0:pver) - real(r8) :: aer_tau700(pcols,0:pver) - - end type rad_out_t - - ! Namelist variables - - integer :: iradsw = -1 ! freq. of shortwave radiation calc in time steps (positive) - ! or hours (negative). - integer :: iradlw = -1 ! frequency of longwave rad. calc. in time steps (positive) - ! or hours (negative). - - integer :: irad_always = 0 ! Specifies length of time in timesteps (positive) - ! or hours (negative) SW/LW radiation will be - ! run continuously from the start of an - ! initial or restart run - logical :: use_rad_dt_cosz = .false. ! if true, use radiation dt for all cosz calculations - logical :: spectralflux = .false. ! calculate fluxes (up and down) per band. - logical :: graupel_in_rad = .false. ! graupel in radiation code - logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation - - - ! Physics buffer indices - integer :: qrs_idx = 0 - integer :: qrl_idx = 0 - integer :: su_idx = 0 - integer :: sd_idx = 0 - integer :: lu_idx = 0 - integer :: ld_idx = 0 - integer :: fsds_idx = 0 - integer :: fsns_idx = 0 - integer :: fsnt_idx = 0 - integer :: flns_idx = 0 - integer :: flnt_idx = 0 - integer :: cldfsnow_idx = 0 - integer :: cld_idx = 0 - integer :: cldfgrau_idx = 0 - integer :: volc_idx = 0 - - character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ','_d6 ','_d7 ','_d8 ','_d9 ','_d10'/) - - ! averaging time interval for zenith angle - real(r8) :: dt_avg = 0._r8 - - real(r8) :: rad_uniform_angle = -99._r8 - - ! PIO descriptors (for restarts) - type(var_desc_t) :: cospcnt_desc - type(var_desc_t) :: nextsw_cday_desc - +! OSLO_AERO end + +implicit none +private +save + +public :: & + radiation_readnl, &! read namelist variables + radiation_register, &! registers radiation physics buffer fields + radiation_nextsw_cday, &! calendar day of next radiation calculation + radiation_do, &! query which radiation calcs are done this timestep + radiation_init, &! initialization + radiation_define_restart, &! define variables for restart + radiation_write_restart, &! write variables to restart + radiation_read_restart, &! read variables from restart + radiation_tend, &! compute heating rates and fluxes + rad_out_t ! type for diagnostic outputs + +integer,public, allocatable :: cosp_cnt(:) ! counter for cosp +integer,public :: cosp_cnt_init = 0 !initial value for cosp counter +real(r8), public, protected :: nextsw_cday ! future radiation calday for surface models + +type rad_out_t + real(r8) :: solin(pcols) ! Solar incident flux + + real(r8) :: qrsc(pcols,pver) + + real(r8) :: fsntc(pcols) ! Clear sky total column abs solar flux + real(r8) :: fsntoa(pcols) ! Net solar flux at TOA + real(r8) :: fsntoac(pcols) ! Clear sky net solar flux at TOA + real(r8) :: fsutoa(pcols) ! upwelling solar flux at TOA + + real(r8) :: fsnirt(pcols) ! Near-IR flux absorbed at toa + real(r8) :: fsnrtc(pcols) ! Clear sky near-IR flux absorbed at toa + real(r8) :: fsnirtsq(pcols) ! Near-IR flux absorbed at toa >= 0.7 microns + + real(r8) :: fsn200(pcols) ! fns interpolated to 200 mb + real(r8) :: fsn200c(pcols) ! fcns interpolated to 200 mb + real(r8) :: fsnr(pcols) ! fns interpolated to tropopause + + real(r8) :: fsnsc(pcols) ! Clear sky surface abs solar flux + real(r8) :: fsdsc(pcols) ! Clear sky surface downwelling solar flux + + real(r8) :: qrlc(pcols,pver) + + real(r8) :: flntc(pcols) ! Clear sky lw flux at model top + real(r8) :: flut(pcols) ! Upward flux at top of model + real(r8) :: flutc(pcols) ! Upward Clear Sky flux at top of model + real(r8) :: lwcf(pcols) ! longwave cloud forcing + + real(r8) :: fln200(pcols) ! net longwave flux interpolated to 200 mb + real(r8) :: fln200c(pcols) ! net clearsky longwave flux interpolated to 200 mb + real(r8) :: flnr(pcols) ! net longwave flux interpolated to tropopause + + real(r8) :: flnsc(pcols) ! Clear sky lw flux at srf (up-down) + real(r8) :: fldsc(pcols) ! Clear sky lw flux at srf (down) + + real(r8) :: tot_cld_vistau(pcols,pver) ! gbx water+ice cloud optical depth (only during day, night = fillvalue) + real(r8) :: tot_icld_vistau(pcols,pver) ! in-cld water+ice cloud optical depth (only during day, night = fillvalue) + real(r8) :: liq_icld_vistau(pcols,pver) ! in-cld liq cloud optical depth (only during day, night = fillvalue) + real(r8) :: ice_icld_vistau(pcols,pver) ! in-cld ice cloud optical depth (only during day, night = fillvalue) + real(r8) :: snow_icld_vistau(pcols,pver) ! snow in-cloud visible sw optical depth for output on history files + real(r8) :: grau_icld_vistau(pcols,pver) ! Graupel in-cloud visible sw optical depth for output on history files + + real(r8) :: cld_tau_cloudsim(pcols,pver) + real(r8) :: aer_tau400(pcols,0:pver) + real(r8) :: aer_tau550(pcols,0:pver) + real(r8) :: aer_tau700(pcols,0:pver) + +end type rad_out_t + +! Namelist variables + +integer :: iradsw = -1 ! freq. of shortwave radiation calc in time steps (positive) + ! or hours (negative). +integer :: iradlw = -1 ! frequency of longwave rad. calc. in time steps (positive) + ! or hours (negative). + +integer :: irad_always = 0 ! Specifies length of time in timesteps (positive) + ! or hours (negative) SW/LW radiation will be + ! run continuously from the start of an + ! initial or restart run +logical :: use_rad_dt_cosz = .false. ! if true, use radiation dt for all cosz calculations +logical :: spectralflux = .false. ! calculate fluxes (up and down) per band. +logical :: graupel_in_rad = .false. ! graupel in radiation code +logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation + + +! Physics buffer indices +integer :: qrs_idx = 0 +integer :: qrl_idx = 0 +integer :: su_idx = 0 +integer :: sd_idx = 0 +integer :: lu_idx = 0 +integer :: ld_idx = 0 +integer :: fsds_idx = 0 +integer :: fsns_idx = 0 +integer :: fsnt_idx = 0 +integer :: flns_idx = 0 +integer :: flnt_idx = 0 +integer :: cldfsnow_idx = 0 +integer :: cld_idx = 0 +integer :: cldfgrau_idx = 0 +integer :: volc_idx = 0 + +character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ','_d6 ','_d7 ','_d8 ','_d9 ','_d10'/) + +! averaging time interval for zenith angle +real(r8) :: dt_avg = 0._r8 + +real(r8) :: rad_uniform_angle = -99._r8 + +! PIO descriptors (for restarts) +type(var_desc_t) :: cospcnt_desc +type(var_desc_t) :: nextsw_cday_desc !=============================================================================== contains !=============================================================================== - subroutine radiation_readnl(nlfile) +subroutine radiation_readnl(nlfile) - ! Read radiation_nl namelist group. + ! Read radiation_nl namelist group. - use namelist_utils, only: find_group_name - use units, only: getunit, freeunit - use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, mpi_real8 + use namelist_utils, only: find_group_name + use units, only: getunit, freeunit + use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_integer, mpi_logical, mpi_real8 - character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input - ! Local variables - integer :: unitn, ierr - integer :: dtime ! timestep size - character(len=*), parameter :: sub = 'radiation_readnl' + ! Local variables + integer :: unitn, ierr + integer :: dtime ! timestep size + character(len=*), parameter :: sub = 'radiation_readnl' - namelist /radiation_nl/ iradsw, iradlw, irad_always, & - use_rad_dt_cosz, spectralflux,graupel_in_rad, use_rad_uniform_angle, rad_uniform_angle + namelist /radiation_nl/ iradsw, iradlw, irad_always, & + use_rad_dt_cosz, spectralflux,graupel_in_rad, use_rad_uniform_angle, rad_uniform_angle - !----------------------------------------------------------------------------- + !----------------------------------------------------------------------------- - if (masterproc) then - unitn = getunit() - open( unitn, file=trim(nlfile), status='old' ) - call find_group_name(unitn, 'radiation_nl', status=ierr) - if (ierr == 0) then - read(unitn, radiation_nl, iostat=ierr) - if (ierr /= 0) then - call endrun(sub // ':: ERROR reading namelist') - end if + if (masterproc) then + unitn = getunit() + open( unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'radiation_nl', status=ierr) + if (ierr == 0) then + read(unitn, radiation_nl, iostat=ierr) + if (ierr /= 0) then + call endrun(sub // ':: ERROR reading namelist') end if - close(unitn) - call freeunit(unitn) end if - - ! Broadcast namelist variables - call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw") - call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw") - call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always") - call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz") - call mpi_bcast(spectralflux, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: spectralflux") - call mpi_bcast(graupel_in_rad, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: graupel_in_rad") - call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle") - call mpi_bcast(rad_uniform_angle, 1, mpi_real8, mstrid, mpicom, ierr) - if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle") - - if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then - call endrun(sub // ' ERROR - use_rad_uniform_angle is set to .true, but rad_uniform_angle is not set ') + close(unitn) + call freeunit(unitn) + end if + + ! Broadcast namelist variables + call mpi_bcast(iradsw, 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradsw") + call mpi_bcast(iradlw, 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: iradlw") + call mpi_bcast(irad_always, 1, mpi_integer, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: irad_always") + call mpi_bcast(use_rad_dt_cosz, 1, mpi_logical, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_dt_cosz") + call mpi_bcast(spectralflux, 1, mpi_logical, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: spectralflux") + call mpi_bcast(graupel_in_rad, 1, mpi_logical, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: graupel_in_rad") + call mpi_bcast(use_rad_uniform_angle, 1, mpi_logical, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: use_rad_uniform_angle") + call mpi_bcast(rad_uniform_angle, 1, mpi_real8, mstrid, mpicom, ierr) + if (ierr /= 0) call endrun(sub//": FATAL: mpi_bcast: rad_uniform_angle") + + if (use_rad_uniform_angle .and. rad_uniform_angle == -99._r8) then + call endrun(sub // ' ERROR - use_rad_uniform_angle is set to .true, but rad_uniform_angle is not set ') + end if + + ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary + dtime = get_step_size() + if (iradsw < 0) iradsw = nint((-iradsw *3600._r8)/dtime) + if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime) + if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime) + + !----------------------------------------------------------------------- + ! Print runtime options to log. + !----------------------------------------------------------------------- + + if (masterproc) then + write(iulog,*) 'RRTMG radiation scheme parameters:' + write(iulog,10) iradsw, iradlw, irad_always, use_rad_dt_cosz, spectralflux, graupel_in_rad + end if + +10 format(' Frequency (timesteps) of Shortwave Radiation calc: ',i5/, & + ' Frequency (timesteps) of Longwave Radiation calc: ',i5/, & + ' SW/LW calc done every timestep for first N steps. N=',i5/, & + ' Use average zenith angle: ',l5/, & + ' Output spectrally resolved fluxes: ',l5/, & + ' Graupel in Radiation Code: ',l5/) + +end subroutine radiation_readnl + +!================================================================================================ + +subroutine radiation_register + + ! Register radiation fields in the physics buffer + + use physics_buffer, only: pbuf_add_field, dtype_r8 + use radiation_data, only: rad_data_register + + call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate + call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate + + call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux + call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux + call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux + + call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux + call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux + + ! If the namelist has been configured for preserving the spectral fluxes, then create + ! physics buffer variables to store the results. + if (spectralflux) then + call pbuf_add_field('SU' , 'global',dtype_r8,(/pcols,pverp,nswbands/), su_idx) ! shortwave upward flux (per band) + call pbuf_add_field('SD' , 'global',dtype_r8,(/pcols,pverp,nswbands/), sd_idx) ! shortwave downward flux (per band) + call pbuf_add_field('LU' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), lu_idx) ! longwave upward flux (per band) + call pbuf_add_field('LD' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), ld_idx) ! longwave downward flux (per band) + end if + + call rad_data_register() + +end subroutine radiation_register + +!================================================================================================ + +function radiation_do(op, timestep) + + ! Return true if the specified operation is done this timestep. + + character(len=*), intent(in) :: op ! name of operation + integer, intent(in), optional:: timestep + logical :: radiation_do ! return value + + ! Local variables + integer :: nstep ! current timestep number + !----------------------------------------------------------------------- + + if (present(timestep)) then + nstep = timestep + else + nstep = get_nstep() + end if + + select case (op) + + case ('sw') ! do a shortwave heating calc this timestep? + radiation_do = nstep == 0 .or. iradsw == 1 & + .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + + case ('lw') ! do a longwave heating calc this timestep? + radiation_do = nstep == 0 .or. iradlw == 1 & + .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & + .or. nstep <= irad_always + + case default + call endrun('radiation_do: unknown operation:'//op) + + end select +end function radiation_do + +!================================================================================================ + +real(r8) function radiation_nextsw_cday() + + ! Return calendar day of next sw radiation calculation + + ! Local variables + integer :: nstep ! timestep counter + logical :: dosw ! true => do shosrtwave calc + integer :: offset ! offset for calendar day calculation + integer :: dtime ! integer timestep size + real(r8):: calday ! calendar day of + real(r8):: caldayp1 ! calendar day of next time-step + !----------------------------------------------------------------------- + + radiation_nextsw_cday = -1._r8 + dosw = .false. + nstep = get_nstep() + dtime = get_step_size() + offset = 0 + do while (.not. dosw) + nstep = nstep + 1 + offset = offset + dtime + if (radiation_do('sw', nstep)) then + radiation_nextsw_cday = get_curr_calday(offset=offset) + dosw = .true. end if - - ! Convert iradsw, iradlw and irad_always from hours to timesteps if necessary + end do + if(radiation_nextsw_cday == -1._r8) then + call endrun('error in radiation_nextsw_cday') + end if + + ! determine if next radiation time-step not equal to next time-step + if (get_nstep() >= 1) then + caldayp1 = get_curr_calday(offset=int(dtime)) + if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8 + end if + +end function radiation_nextsw_cday + +!================================================================================================ + +subroutine radiation_init(pbuf2d) + + ! Initialize the radiation parameterization, add fields to the history buffer + + use physics_buffer, only: pbuf_get_index, pbuf_set_field + use phys_control, only: phys_getopts + use radsw, only: radsw_init + use radlw, only: radlw_init + use rad_solar_var, only: rad_solar_var_init + use radiation_data, only: rad_data_init + use cloud_rad_props, only: cloud_rad_props_init + use modal_aer_opt, only: modal_aer_opt_init + use rrtmg_state, only: rrtmg_state_init + use time_manager, only: is_first_step + + + ! arguments + type(physics_buffer_desc), pointer :: pbuf2d(:,:) + + ! local variables + integer :: icall, nmodes + logical :: active_calls(0:N_DIAG) + integer :: nstep ! current timestep number + logical :: history_amwg ! output the variables used by the AMWG diag package + logical :: history_vdiag ! output the variables used by the AMWG variability diag package + logical :: history_budget ! output tendencies and state variables for CAM4 + ! temperature, water vapor, cloud ice and cloud + ! liquid budgets. + integer :: history_budget_histfile_num ! output history file number for budget fields + integer :: err + + integer :: dtime + !----------------------------------------------------------------------- + + call rad_solar_var_init() + call rrtmg_state_init() + call rad_data_init(pbuf2d) ! initialize output fields for offline driver + call radsw_init() + call radlw_init() + call cloud_rad_props_init() + + cld_idx = pbuf_get_index('CLD') + cldfsnow_idx = pbuf_get_index('CLDFSNOW',errcode=err) + cldfgrau_idx = pbuf_get_index('CLDFGRAU',errcode=err) + + if (is_first_step()) then + call pbuf_set_field(pbuf2d, qrl_idx, 0._r8) + end if + + ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation + if (use_rad_dt_cosz) then dtime = get_step_size() - if (iradsw < 0) iradsw = nint((-iradsw *3600._r8)/dtime) - if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime) - if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime) - - !----------------------------------------------------------------------- - ! Print runtime options to log. - !----------------------------------------------------------------------- + dt_avg = iradsw*dtime + end if + + ! Surface components to get radiation computed today + if (.not. is_first_restart_step()) then + nextsw_cday = get_curr_calday() + end if + + call phys_getopts(history_amwg_out = history_amwg, & + history_vdiag_out = history_vdiag, & + history_budget_out = history_budget, & + history_budget_histfile_num_out = history_budget_histfile_num) + + ! Determine whether modal aerosols are affecting the climate, and if so + ! then initialize the modal aerosol optics module + call rad_cnst_get_info(0, nmodes=nmodes) + if (nmodes > 0) call modal_aer_opt_init() + + ! "irad_always" is number of time steps to execute radiation continuously from start of + ! initial OR restart run + nstep = get_nstep() + if (irad_always > 0) then + nstep = get_nstep() + irad_always = irad_always + nstep + end if + + if (docosp) call cospsimulator_intr_init + + allocate(cosp_cnt(begchunk:endchunk)) + if (is_first_restart_step()) then + cosp_cnt(begchunk:endchunk) = cosp_cnt_init + else + cosp_cnt(begchunk:endchunk) = 0 + end if + + call addfld('O3colAbove', horiz_only, 'A', 'DU', 'Column O3 above model top', sampling_seq='rad_lwsw') + + call addfld('TOT_CLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total gbx cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total in-cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Liquid in-cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Ice in-cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + + if (cldfsnow_idx > 0) then + call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Snow in-cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + endif + if (cldfgrau_idx > 0 .and. graupel_in_rad) then + call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Graupel in-cloud extinction visible sw optical depth', & + sampling_seq='rad_lwsw', flag_xyfill=.true.) + endif + + ! get list of active radiation calls + call rad_cnst_get_call_list(active_calls) + + ! Add shortwave radiation fields to history master field list. + + do icall = 0, N_DIAG + + if (active_calls(icall)) then + + call addfld('SOLIN'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar insolation', sampling_seq='rad_lwsw') + + call addfld('QRS'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Solar heating rate', sampling_seq='rad_lwsw') + call addfld('QRSC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky solar heating rate', & + sampling_seq='rad_lwsw') + call addfld('FSNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('FSNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('FSNTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of atmosphere', & + sampling_seq='rad_lwsw') + call addfld('FSNTOAC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of atmosphere', & + sampling_seq='rad_lwsw') + call addfld('SWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Shortwave cloud forcing', & + sampling_seq='rad_lwsw') + call addfld('FSUTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling solar flux at top of atmosphere', & + sampling_seq='rad_lwsw') + call addfld('FSNIRTOA'//diag(icall), horiz_only, 'A', 'W/m2', & + 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw') + call addfld('FSNRTOAC'//diag(icall), horiz_only, 'A', 'W/m2', & + 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw') + call addfld('FSNRTOAS'//diag(icall), horiz_only, 'A', 'W/m2', & + 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw') + + call addfld('FSN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net shortwave flux at 200 mb', & + sampling_seq='rad_lwsw') + call addfld('FSN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net shortwave flux at 200 mb', & + sampling_seq='rad_lwsw') + + call addfld('FSNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at tropopause', & + sampling_seq='rad_lwsw') + + call addfld('SOLL'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared direct to surface', & + sampling_seq='rad_lwsw') + call addfld('SOLS'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible direct to surface', & + sampling_seq='rad_lwsw') + call addfld('SOLLD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared diffuse to surface', & + sampling_seq='rad_lwsw') + call addfld('SOLSD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible diffuse to surface', & + sampling_seq='rad_lwsw') + call addfld('FSNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FSNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at surface', & + sampling_seq='rad_lwsw') + + call addfld('FSDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling solar flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FSDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky downwelling solar flux at surface', & + sampling_seq='rad_lwsw') + + call addfld('FUS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave upward flux') + call addfld('FDS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave downward flux') + call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky upward flux') + call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky downward flux') + + if (history_amwg) then + call add_default('SOLIN'//diag(icall), 1, ' ') + call add_default('QRS'//diag(icall), 1, ' ') + call add_default('FSNT'//diag(icall), 1, ' ') + call add_default('FSNTC'//diag(icall), 1, ' ') + call add_default('FSNTOA'//diag(icall), 1, ' ') + call add_default('FSNTOAC'//diag(icall), 1, ' ') + call add_default('SWCF'//diag(icall), 1, ' ') + call add_default('FSNS'//diag(icall), 1, ' ') + call add_default('FSNSC'//diag(icall), 1, ' ') + call add_default('FSUTOA'//diag(icall), 1, ' ') + call add_default('FSDSC'//diag(icall), 1, ' ') + call add_default('FSDS'//diag(icall), 1, ' ') + endif - if (masterproc) then - write(iulog,*) 'RRTMG radiation scheme parameters:' - write(iulog,10) iradsw, iradlw, irad_always, use_rad_dt_cosz, spectralflux, graupel_in_rad end if + end do + ! OSLO_AERO begin + call addfld('FDSCDRF', (/ 'ilev' /), 'A', 'W/m2', 'Shortwave clear-sky downward flux') + call addfld('FUSCDRF', (/ 'ilev' /), 'A', 'W/m2', 'Shortwave clear-sky upward flux') + ! OSLO_AERO end + + if (scm_crm_mode) then + call add_default('FUS ', 1, ' ') + call add_default('FUSC ', 1, ' ') + call add_default('FDS ', 1, ' ') + call add_default('FDSC ', 1, ' ') + endif + + ! Add longwave radiation fields to history master field list. + + do icall = 0, N_DIAG + + if (active_calls(icall)) then + + call addfld('QRL'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Longwave heating rate', sampling_seq='rad_lwsw') + call addfld('QRLC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky longwave heating rate', & + sampling_seq='rad_lwsw') + call addfld('FLNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('FLNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('FLNTCLR'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky ONLY points net longwave flux at top of model',& + sampling_seq='rad_lwsw') + call addfld('FREQCLR'//diag(icall), horiz_only, 'A', 'Frac', 'Frequency of Occurrence of Clearsky', & + sampling_seq='rad_lwsw') + call addfld('FLUT'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling longwave flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('FLUTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky upwelling longwave flux at top of model', & + sampling_seq='rad_lwsw') + call addfld('LWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Longwave cloud forcing', sampling_seq='rad_lwsw') + + call addfld('FLN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at 200 mb', & + sampling_seq='rad_lwsw') + call addfld('FLN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at 200 mb', & + sampling_seq='rad_lwsw') + call addfld('FLNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at tropopause', & + sampling_seq='rad_lwsw') + + call addfld('FLNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FLNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FLDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling longwave flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FLDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky Downwelling longwave flux at surface', & + sampling_seq='rad_lwsw') + call addfld('FUL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave upward flux') + call addfld('FDL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave downward flux') + call addfld('FULC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky upward flux') + call addfld('FDLC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky downward flux') + + if (history_amwg) then + call add_default('QRL'//diag(icall), 1, ' ') + call add_default('FLNT'//diag(icall), 1, ' ') + call add_default('FLNTC'//diag(icall), 1, ' ') + call add_default('FLNTCLR'//diag(icall), 1, ' ') + call add_default('FREQCLR'//diag(icall), 1, ' ') + call add_default('FLUT'//diag(icall), 1, ' ') + call add_default('FLUTC'//diag(icall), 1, ' ') + call add_default('LWCF'//diag(icall), 1, ' ') + call add_default('FLNS'//diag(icall), 1, ' ') + call add_default('FLNSC'//diag(icall), 1, ' ') + call add_default('FLDS'//diag(icall), 1, ' ') + endif -10 format(' Frequency (timesteps) of Shortwave Radiation calc: ',i5/, & - ' Frequency (timesteps) of Longwave Radiation calc: ',i5/, & - ' SW/LW calc done every timestep for first N steps. N=',i5/, & - ' Use average zenith angle: ',l5/, & - ' Output spectrally resolved fluxes: ',l5/, & - ' Graupel in Radiation Code: ',l5/) - - end subroutine radiation_readnl - - !================================================================================================ - - subroutine radiation_register - - ! Register radiation fields in the physics buffer - - use physics_buffer, only: pbuf_add_field, dtype_r8 - use radiation_data, only: rad_data_register + end if + end do - call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate - call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate + call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity') - call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux - call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux - call pbuf_add_field('FSNT' , 'global',dtype_r8,(/pcols/), fsnt_idx) ! Top-of-model net shortwave flux + if (scm_crm_mode) then + call add_default ('FUL ', 1, ' ') + call add_default ('FULC ', 1, ' ') + call add_default ('FDL ', 1, ' ') + call add_default ('FDLC ', 1, ' ') + endif - call pbuf_add_field('FLNS' , 'global',dtype_r8,(/pcols/), flns_idx) ! Surface net longwave flux - call pbuf_add_field('FLNT' , 'global',dtype_r8,(/pcols/), flnt_idx) ! Top-of-model net longwave flux + ! Heating rate needed for d(theta)/dt computation + call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation') - ! If the namelist has been configured for preserving the spectral fluxes, then create - ! physics buffer variables to store the results. - if (spectralflux) then - call pbuf_add_field('SU' , 'global',dtype_r8,(/pcols,pverp,nswbands/), su_idx) ! shortwave upward flux (per band) - call pbuf_add_field('SD' , 'global',dtype_r8,(/pcols,pverp,nswbands/), sd_idx) ! shortwave downward flux (per band) - call pbuf_add_field('LU' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), lu_idx) ! longwave upward flux (per band) - call pbuf_add_field('LD' , 'global',dtype_r8,(/pcols,pverp,nlwbands/), ld_idx) ! longwave downward flux (per band) - end if + if ( history_budget .and. history_budget_histfile_num > 1 ) then + call add_default ('QRL ', history_budget_histfile_num, ' ') + call add_default ('QRS ', history_budget_histfile_num, ' ') + end if - call rad_data_register() + if (history_vdiag) then + call add_default('FLUT', 2, ' ') + call add_default('FLUT', 3, ' ') + end if - end subroutine radiation_register +end subroutine radiation_init - !================================================================================================ +!=============================================================================== - function radiation_do(op, timestep) +subroutine radiation_define_restart(file) - ! Return true if the specified operation is done this timestep. + ! define variables to be written to restart file - character(len=*), intent(in) :: op ! name of operation - integer, intent(in), optional:: timestep - logical :: radiation_do ! return value + ! arguments + type(file_desc_t), intent(inout) :: file - ! Local variables - integer :: nstep ! current timestep number - !----------------------------------------------------------------------- + ! local variables + integer :: ierr + !---------------------------------------------------------------------------- - if (present(timestep)) then - nstep = timestep - else - nstep = get_nstep() - end if + call pio_seterrorhandling(File, PIO_BCAST_ERROR) - select case (op) + ierr = pio_def_var(File, 'nextsw_cday', pio_double, nextsw_cday_desc) + ierr = pio_put_att(File, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models') + if (docosp) then + ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc) + end if - case ('sw') ! do a shortwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradsw == 1 & - .or. (mod(nstep-1,iradsw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always +end subroutine radiation_define_restart - case ('lw') ! do a longwave heating calc this timestep? - radiation_do = nstep == 0 .or. iradlw == 1 & - .or. (mod(nstep-1,iradlw) == 0 .and. nstep /= 1) & - .or. nstep <= irad_always +!=============================================================================== - case default - call endrun('radiation_do: unknown operation:'//op) +subroutine radiation_write_restart(file) - end select - end function radiation_do + ! write variables to restart file - !================================================================================================ + ! arguments + type(file_desc_t), intent(inout) :: file - real(r8) function radiation_nextsw_cday() + ! local variables + integer :: ierr + !---------------------------------------------------------------------------- - ! Return calendar day of next sw radiation calculation + ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /)) + if (docosp) then + ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/)) + end if - ! Local variables - integer :: nstep ! timestep counter - logical :: dosw ! true => do shosrtwave calc - integer :: offset ! offset for calendar day calculation - integer :: dtime ! integer timestep size - real(r8):: calday ! calendar day of - real(r8):: caldayp1 ! calendar day of next time-step - !----------------------------------------------------------------------- +end subroutine radiation_write_restart - radiation_nextsw_cday = -1._r8 - dosw = .false. - nstep = get_nstep() - dtime = get_step_size() - offset = 0 - do while (.not. dosw) - nstep = nstep + 1 - offset = offset + dtime - if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) - dosw = .true. - end if - end do - if(radiation_nextsw_cday == -1._r8) then - call endrun('error in radiation_nextsw_cday') - end if +!=============================================================================== - ! determine if next radiation time-step not equal to next time-step - if (get_nstep() >= 1) then - caldayp1 = get_curr_calday(offset=int(dtime)) - if (caldayp1 /= radiation_nextsw_cday) radiation_nextsw_cday = -1._r8 - end if +subroutine radiation_read_restart(file) - end function radiation_nextsw_cday - - !================================================================================================ - - subroutine radiation_init(pbuf2d) - - ! Initialize the radiation parameterization, add fields to the history buffer - - use physics_buffer, only: pbuf_get_index, pbuf_set_field - use phys_control, only: phys_getopts - use radsw, only: radsw_init - use radlw, only: radlw_init - use rad_solar_var, only: rad_solar_var_init - use radiation_data, only: rad_data_init - use cloud_rad_props, only: cloud_rad_props_init - use modal_aer_opt, only: modal_aer_opt_init - use rrtmg_state, only: rrtmg_state_init - use time_manager, only: is_first_step - - - ! arguments - type(physics_buffer_desc), pointer :: pbuf2d(:,:) - - ! local variables - integer :: icall, nmodes - logical :: active_calls(0:N_DIAG) - integer :: nstep ! current timestep number - logical :: history_amwg ! output the variables used by the AMWG diag package - logical :: history_vdiag ! output the variables used by the AMWG variability diag package - logical :: history_budget ! output tendencies and state variables for CAM4 - ! temperature, water vapor, cloud ice and cloud - ! liquid budgets. - integer :: history_budget_histfile_num ! output history file number for budget fields - integer :: err - - integer :: dtime - !----------------------------------------------------------------------- - - call rad_solar_var_init() - call rrtmg_state_init() - call rad_data_init(pbuf2d) ! initialize output fields for offline driver - call radsw_init() - call radlw_init() - call cloud_rad_props_init() - - cld_idx = pbuf_get_index('CLD') - cldfsnow_idx = pbuf_get_index('CLDFSNOW',errcode=err) - cldfgrau_idx = pbuf_get_index('CLDFGRAU',errcode=err) - - if (is_first_step()) then - call pbuf_set_field(pbuf2d, qrl_idx, 0._r8) - end if + ! read variables from restart file - ! Set the radiation timestep for cosz calculations if requested using the adjusted iradsw value from radiation - if (use_rad_dt_cosz) then - dtime = get_step_size() - dt_avg = iradsw*dtime - end if + ! arguments + type(file_desc_t), intent(inout) :: file - ! Surface components to get radiation computed today - if (.not. is_first_restart_step()) then - nextsw_cday = get_curr_calday() - end if + ! local variables - call phys_getopts(history_amwg_out = history_amwg, & - history_vdiag_out = history_vdiag, & - history_budget_out = history_budget, & - history_budget_histfile_num_out = history_budget_histfile_num) + integer :: err_handling + integer :: ierr + real(r8) :: temp_var - ! Determine whether modal aerosols are affecting the climate, and if so - ! then initialize the modal aerosol optics module - call rad_cnst_get_info(0, nmodes=nmodes) - if (nmodes > 0) call modal_aer_opt_init() + type(var_desc_t) :: vardesc + !---------------------------------------------------------------------------- - ! "irad_always" is number of time steps to execute radiation continuously from start of - ! initial OR restart run - nstep = get_nstep() - if (irad_always > 0) then - nstep = get_nstep() - irad_always = irad_always + nstep - end if - - if (docosp) call cospsimulator_intr_init - - allocate(cosp_cnt(begchunk:endchunk)) - if (is_first_restart_step()) then - cosp_cnt(begchunk:endchunk) = cosp_cnt_init + if (docosp) then + call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling) + ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc) + call pio_seterrorhandling(File, err_handling) + if (ierr /= PIO_NOERR) then + cosp_cnt_init = 0 else - cosp_cnt(begchunk:endchunk) = 0 + ierr = pio_get_var(File, vardesc, cosp_cnt_init) end if + end if - call addfld('O3colAbove', horiz_only, 'A', 'DU', 'Column O3 above model top', sampling_seq='rad_lwsw') - - call addfld('TOT_CLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total gbx cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) - call addfld('TOT_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Total in-cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) - call addfld('LIQ_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Liquid in-cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) - call addfld('ICE_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Ice in-cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) + ierr = pio_inq_varid(File, 'nextsw_cday', vardesc) + ierr = pio_get_var(File, vardesc, temp_var) + nextsw_cday = temp_var - if (cldfsnow_idx > 0) then - call addfld('SNOW_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Snow in-cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) - endif - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - call addfld('GRAU_ICLD_VISTAU', (/ 'lev' /), 'A', '1', 'Graupel in-cloud extinction visible sw optical depth', & - sampling_seq='rad_lwsw', flag_xyfill=.true.) - endif - - ! get list of active radiation calls - call rad_cnst_get_call_list(active_calls) - - ! Add shortwave radiation fields to history master field list. - - do icall = 0, N_DIAG - - if (active_calls(icall)) then - - call addfld('SOLIN'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar insolation', sampling_seq='rad_lwsw') - - call addfld('QRS'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Solar heating rate', sampling_seq='rad_lwsw') - call addfld('QRSC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky solar heating rate', & - sampling_seq='rad_lwsw') - call addfld('FSNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('FSNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('FSNTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at top of atmosphere', & - sampling_seq='rad_lwsw') - call addfld('FSNTOAC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at top of atmosphere', & - sampling_seq='rad_lwsw') - call addfld('SWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Shortwave cloud forcing', & - sampling_seq='rad_lwsw') - call addfld('FSUTOA'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling solar flux at top of atmosphere', & - sampling_seq='rad_lwsw') - call addfld('FSNIRTOA'//diag(icall), horiz_only, 'A', 'W/m2', & - 'Net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw') - call addfld('FSNRTOAC'//diag(icall), horiz_only, 'A', 'W/m2', & - 'Clearsky net near-infrared flux (Nimbus-7 WFOV) at top of atmosphere', sampling_seq='rad_lwsw') - call addfld('FSNRTOAS'//diag(icall), horiz_only, 'A', 'W/m2', & - 'Net near-infrared flux (>= 0.7 microns) at top of atmosphere', sampling_seq='rad_lwsw') - - call addfld('FSN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net shortwave flux at 200 mb', & - sampling_seq='rad_lwsw') - call addfld('FSN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net shortwave flux at 200 mb', & - sampling_seq='rad_lwsw') - - call addfld('FSNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at tropopause', & - sampling_seq='rad_lwsw') - - call addfld('SOLL'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared direct to surface', & - sampling_seq='rad_lwsw') - call addfld('SOLS'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible direct to surface', & - sampling_seq='rad_lwsw') - call addfld('SOLLD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward near infrared diffuse to surface', & - sampling_seq='rad_lwsw') - call addfld('SOLSD'//diag(icall), horiz_only, 'A', 'W/m2', 'Solar downward visible diffuse to surface', & - sampling_seq='rad_lwsw') - call addfld('FSNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net solar flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FSNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net solar flux at surface', & - sampling_seq='rad_lwsw') - - call addfld('FSDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling solar flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FSDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky downwelling solar flux at surface', & - sampling_seq='rad_lwsw') - - call addfld('FUS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave upward flux') - call addfld('FDS'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave downward flux') - call addfld('FUSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky upward flux') - call addfld('FDSC'//diag(icall), (/ 'ilev' /), 'I', 'W/m2', 'Shortwave clear-sky downward flux') - - if (history_amwg) then - call add_default('SOLIN'//diag(icall), 1, ' ') - call add_default('QRS'//diag(icall), 1, ' ') - call add_default('FSNT'//diag(icall), 1, ' ') - call add_default('FSNTC'//diag(icall), 1, ' ') - call add_default('FSNTOA'//diag(icall), 1, ' ') - call add_default('FSNTOAC'//diag(icall), 1, ' ') - call add_default('SWCF'//diag(icall), 1, ' ') - call add_default('FSNS'//diag(icall), 1, ' ') - call add_default('FSNSC'//diag(icall), 1, ' ') - call add_default('FSUTOA'//diag(icall), 1, ' ') - call add_default('FSDSC'//diag(icall), 1, ' ') - call add_default('FSDS'//diag(icall), 1, ' ') - endif - - end if - end do +end subroutine radiation_read_restart -#ifdef OSLO_AERO - call addfld('FDSCDRF', (/ 'ilev' /), 'A', 'W/m2', 'Shortwave clear-sky downward flux') - call addfld('FUSCDRF', (/ 'ilev' /), 'A', 'W/m2', 'Shortwave clear-sky upward flux') -#endif +!=============================================================================== - if (scm_crm_mode) then - call add_default('FUS ', 1, ' ') - call add_default('FUSC ', 1, ' ') - call add_default('FDS ', 1, ' ') - call add_default('FDSC ', 1, ' ') - endif - - ! Add longwave radiation fields to history master field list. - - do icall = 0, N_DIAG - - if (active_calls(icall)) then - - call addfld('QRL'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Longwave heating rate', sampling_seq='rad_lwsw') - call addfld('QRLC'//diag(icall), (/ 'lev' /), 'A', 'K/s', 'Clearsky longwave heating rate', & - sampling_seq='rad_lwsw') - call addfld('FLNT'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('FLNTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('FLNTCLR'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky ONLY points net longwave flux at top of model',& - sampling_seq='rad_lwsw') - call addfld('FREQCLR'//diag(icall), horiz_only, 'A', 'Frac', 'Frequency of Occurrence of Clearsky', & - sampling_seq='rad_lwsw') - call addfld('FLUT'//diag(icall), horiz_only, 'A', 'W/m2', 'Upwelling longwave flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('FLUTC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky upwelling longwave flux at top of model', & - sampling_seq='rad_lwsw') - call addfld('LWCF'//diag(icall), horiz_only, 'A', 'W/m2', 'Longwave cloud forcing', sampling_seq='rad_lwsw') - - call addfld('FLN200'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at 200 mb', & - sampling_seq='rad_lwsw') - call addfld('FLN200C'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at 200 mb', & - sampling_seq='rad_lwsw') - call addfld('FLNR'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at tropopause', & - sampling_seq='rad_lwsw') - - call addfld('FLNS'//diag(icall), horiz_only, 'A', 'W/m2', 'Net longwave flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FLNSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky net longwave flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FLDS'//diag(icall), horiz_only, 'A', 'W/m2', 'Downwelling longwave flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FLDSC'//diag(icall), horiz_only, 'A', 'W/m2', 'Clearsky Downwelling longwave flux at surface', & - sampling_seq='rad_lwsw') - call addfld('FUL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave upward flux') - call addfld('FDL'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave downward flux') - call addfld('FULC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky upward flux') - call addfld('FDLC'//diag(icall), (/ 'ilev' /),'I', 'W/m2', 'Longwave clear-sky downward flux') - - if (history_amwg) then - call add_default('QRL'//diag(icall), 1, ' ') - call add_default('FLNT'//diag(icall), 1, ' ') - call add_default('FLNTC'//diag(icall), 1, ' ') - call add_default('FLNTCLR'//diag(icall), 1, ' ') - call add_default('FREQCLR'//diag(icall), 1, ' ') - call add_default('FLUT'//diag(icall), 1, ' ') - call add_default('FLUTC'//diag(icall), 1, ' ') - call add_default('LWCF'//diag(icall), 1, ' ') - call add_default('FLNS'//diag(icall), 1, ' ') - call add_default('FLNSC'//diag(icall), 1, ' ') - call add_default('FLDS'//diag(icall), 1, ' ') - endif +subroutine radiation_tend( & + state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out) - end if + !----------------------------------------------------------------------- + ! + ! Driver for radiation computation. + ! + ! Revision history: + ! 2007-11-05 M. Iacono Install rrtmg_lw and sw as radiation model. + ! 2007-12-27 M. Iacono Modify to use CAM cloud optical properties with rrtmg. + !----------------------------------------------------------------------- + + use phys_grid, only: get_rlat_all_p, get_rlon_all_p + use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr + use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz + + use cloud_rad_props, only: get_ice_optics_sw, get_liquid_optics_sw, liquid_cloud_get_rad_props_lw, & + ice_cloud_get_rad_props_lw, cloud_rad_props_get_lw, & + grau_cloud_get_rad_props_lw, get_grau_optics_sw, & + snow_cloud_get_rad_props_lw, get_snow_optics_sw + use slingo, only: slingo_liq_get_rad_props_lw, slingo_liq_optics_sw + use ebert_curry, only: ec_ice_optics_sw, ec_ice_get_rad_props_lw + + use rad_solar_var, only: get_variability + use radsw, only: rad_rrtmg_sw + use radlw, only: rad_rrtmg_lw + use radheat, only: radheat_tend + + use radiation_data, only: rad_data_write + use rrtmg_state, only: rrtmg_state_create, rrtmg_state_update, rrtmg_state_destroy, rrtmg_state_t, & + num_rrtmg_levs + + use interpolate_data, only: vertinterp + use tropopause, only: tropopause_find, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE + + use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps + ! OSLO_AERO begin + use physics_buffer, only: pbuf_get_index, pbuf_get_field + use oslo_aero_control, only: oslo_aero_getopts + use oslo_aero_share + ! OSLO_AERO end + + ! Arguments + type(physics_state), intent(in), target :: state + type(physics_ptend), intent(out) :: ptend + + type(physics_buffer_desc), pointer :: pbuf(:) + type(cam_out_t), intent(inout) :: cam_out + type(cam_in_t), intent(in) :: cam_in + real(r8), intent(out) :: net_flx(pcols) + + type(rad_out_t), target, optional, intent(out) :: rd_out + + ! Local variables + real(r8) :: volc_fraction_coarse ! Fraction of volcanic aerosols going to coarse mode + real(r8), pointer :: rvolcmmr(:,:) ! stratospheric volcanoes aerosol mmr + integer :: band + logical :: idrf + type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object + ! if the argument is not present + logical :: write_output + + integer :: i, k + integer :: lchnk, ncol + logical :: dosw, dolw + real(r8) :: calday ! current calendar day + real(r8) :: delta ! Solar declination angle in radians + real(r8) :: eccf ! Earth orbit eccentricity factor + real(r8) :: clat(pcols) ! current latitudes(radians) + real(r8) :: clon(pcols) ! current longitudes(radians) + real(r8) :: coszrs(pcols) ! Cosine solar zenith angle + + ! Gathered indices of day and night columns + ! chunk_column_index = IdxDay(daylight_column_index) + integer :: Nday ! Number of daylight columns + integer :: Nnite ! Number of night columns + integer :: IdxDay(pcols) ! Indices of daylight columns + integer :: IdxNite(pcols) ! Indices of night columns + + integer :: itim_old + + real(r8), pointer :: cld(:,:) ! cloud fraction + real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds- whatever they are" + real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "snow clouds- whatever they are" + real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate + real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate + real(r8), pointer :: fsds(:) ! Surface solar down flux + real(r8), pointer :: fsns(:) ! Surface solar absorbed flux + real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top + real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux + real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top + + real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up + real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down + real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up + real(r8), pointer, dimension(:,:,:) :: ld => NULL() ! longwave spectral flux down + + ! tropopause diagnostic + integer :: troplev(pcols) + real(r8) :: p_trop(pcols) + + type(rrtmg_state_t) :: r_state ! contains the atm concentrations in layers needed for RRTMG + + ! cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: ice_tau (nswbands,pcols,pver) ! ice extinction optical depth + real(r8) :: ice_tau_w (nswbands,pcols,pver) ! ice single scattering albedo * tau + real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice assymetry parameter * tau * w + real(r8) :: ice_tau_w_f(nswbands,pcols,pver) ! ice forward scattered fraction * tau * w + real(r8) :: ice_lw_abs (nlwbands,pcols,pver) ! ice absorption optics depth (LW) + + ! cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth + real(r8) :: liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau + real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid assymetry parameter * tau * w + real(r8) :: liq_tau_w_f(nswbands,pcols,pver) ! liquid forward scattered fraction * tau * w + real(r8) :: liq_lw_abs (nlwbands,pcols,pver) ! liquid absorption optics depth (LW) + + ! cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: cld_tau (nswbands,pcols,pver) ! cloud extinction optical depth + real(r8) :: cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau + real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud assymetry parameter * w * tau + real(r8) :: cld_tau_w_f(nswbands,pcols,pver) ! cloud forward scattered fraction * w * tau + real(r8) :: cld_lw_abs (nlwbands,pcols,pver) ! cloud absorption optics depth (LW) + + ! cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth + real(r8) :: snow_tau_w (nswbands,pcols,pver) ! snow single scattering albedo * tau + real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow assymetry parameter * tau * w + real(r8) :: snow_tau_w_f(nswbands,pcols,pver) ! snow forward scattered fraction * tau * w + real(r8) :: snow_lw_abs (nlwbands,pcols,pver)! snow absorption optics depth (LW) + + ! Add graupel as another snow species. + ! cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: grau_tau (nswbands,pcols,pver) ! graupel extinction optical depth + real(r8) :: grau_tau_w (nswbands,pcols,pver) ! graupel single scattering albedo * tau + real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel assymetry parameter * tau * w + real(r8) :: grau_tau_w_f(nswbands,pcols,pver) ! graupel forward scattered fraction * tau * w + real(r8) :: grau_lw_abs (nlwbands,pcols,pver)! graupel absorption optics depth (LW) + + ! combined cloud radiative parameters are "in cloud" not "in cell" + real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction (snow plus regular) + real(r8) :: c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth + real(r8) :: c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau + real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud assymetry parameter * w * tau + real(r8) :: c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau + real(r8) :: c_cld_lw_abs (nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW) + + real(r8) :: sfac(1:nswbands) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band + + integer :: icall ! index through climate/diagnostic radiation calls + logical :: active_calls(0:N_DIAG) + + ! OSLO_AERO begin + ! Local variables used for calculating aerosol optics and direct and indirect forcings. + ! aodvis and absvis are AOD and absorptive AOD for visible wavelength close to 0.55 um (0.35-0.64) + ! Note that aodvis and absvis output should be devided by dayfoc to give physical (A)AOD values + real(r8) :: qdirind(pcols,pver,pcnst) ! Common tracers for indirect and direct calculations + real(r8) :: aodvis(pcols) ! AOD vis + real(r8) :: absvis(pcols) ! absorptive AOD vis + real(r8) :: clearodvis(pcols) + real(r8) :: clearabsvis(pcols) + real(r8) :: cloudfree(pcols) + real(r8) :: cloudfreemax(pcols) + real(r8) :: clearod440(pcols) + real(r8) :: clearod550(pcols) + real(r8) :: clearod870(pcols) + real(r8) :: clearabs550(pcols) + real(r8) :: clearabs550alt(pcols) + real(r8) :: ftem_1d(pcols) ! work-array to avoid NAN and pcols/ncol confusion + real(r8) :: volc_ext_sun(pcols,pver,nswbands) ! volcanic aerosol extinction for solar bands, CMIP6 + real(r8) :: volc_omega_sun(pcols,pver,nswbands) ! volcanic aerosol SSA for solar bands, CMIP6 + real(r8) :: volc_g_sun(pcols,pver,nswbands) ! volcanic aerosol g for solar bands, CMIP6 + real(r8) :: volc_ext_earth(pcols,pver,nlwbands) ! volcanic aerosol extinction for terrestrial bands, CMIP6 + real(r8) :: volc_omega_earth(pcols,pver,nlwbands) ! volcanic aerosol SSA for terrestrial bands, CMIP6 + ! OSLO_AERO end + + ! Aerosol radiative properties + real(r8) :: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth + real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau + real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau + real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau + real(r8) :: aer_lw_abs (pcols,pver,nlwbands) ! aerosol absorption optics depth (LW) + + real(r8) :: fns(pcols,pverp) ! net shortwave flux + real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux + real(r8) :: fnl(pcols,pverp) ! net longwave flux + real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux + + ! for COSP + real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity + real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau + real(r8) :: gb_snow_lw(pcols,pver) ! grid-box mean LW snow optical depth + + real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables + + real(r8) :: freqclr(pcols) ! Frequency of occurrence of clear sky columns + real(r8) :: flntclr(pcols) ! Clearsky only columns (zero if cloudy) + + character(*), parameter :: name = 'radiation_tend' + !-------------------------------------------------------------------------------------- + + ! OSLO_AERO begin + aer_lw_abs(:,:,:) =0._r8 + aer_tau(:,:,:) =0._r8 + aer_tau_w(:,:,:) =0._r8 + aer_tau_w_g(:,:,:) =0._r8 + aer_tau_w_f(:,:,:) =0._r8 + ! OSLO_AERO end + + lchnk = state%lchnk + ncol = state%ncol + + if (present(rd_out)) then + rd => rd_out + write_output = .false. + else + allocate(rd) + write_output=.true. + end if + + dosw = radiation_do('sw') ! do shortwave heating calc this timestep? + dolw = radiation_do('lw') ! do longwave heating calc this timestep? + + ! Cosine solar zenith angle for current time step + calday = get_curr_calday() + call get_rlat_all_p(lchnk, ncol, clat) + call get_rlon_all_p(lchnk, ncol, clon) + + call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, & + delta, eccf) + + if (use_rad_uniform_angle) then + do i = 1, ncol + coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, uniform_angle=rad_uniform_angle) end do - - call addfld('EMIS', (/ 'lev' /), 'A', '1', 'Cloud longwave emissivity') - - if (scm_crm_mode) then - call add_default ('FUL ', 1, ' ') - call add_default ('FULC ', 1, ' ') - call add_default ('FDL ', 1, ' ') - call add_default ('FDLC ', 1, ' ') - endif - - ! Heating rate needed for d(theta)/dt computation - call addfld ('HR',(/ 'lev' /), 'A','K/s','Heating rate needed for d(theta)/dt computation') - - if ( history_budget .and. history_budget_histfile_num > 1 ) then - call add_default ('QRL ', history_budget_histfile_num, ' ') - call add_default ('QRS ', history_budget_histfile_num, ' ') - end if - - if (history_vdiag) then - call add_default('FLUT', 2, ' ') - call add_default('FLUT', 3, ' ') - end if - - end subroutine radiation_init - - !=============================================================================== - - subroutine radiation_define_restart(file) - - ! define variables to be written to restart file - - ! arguments - type(file_desc_t), intent(inout) :: file - - ! local variables - integer :: ierr - !---------------------------------------------------------------------------- - - call pio_seterrorhandling(File, PIO_BCAST_ERROR) - - ierr = pio_def_var(File, 'nextsw_cday', pio_double, nextsw_cday_desc) - ierr = pio_put_att(File, nextsw_cday_desc, 'long_name', 'future radiation calday for surface models') - if (docosp) then - ierr = pio_def_var(File, 'cosp_cnt_init', pio_int, cospcnt_desc) - end if - - end subroutine radiation_define_restart - - !=============================================================================== - - subroutine radiation_write_restart(file) - - ! write variables to restart file - - ! arguments - type(file_desc_t), intent(inout) :: file - - ! local variables - integer :: ierr - !---------------------------------------------------------------------------- - - ierr = pio_put_var(File, nextsw_cday_desc, (/ nextsw_cday /)) - if (docosp) then - ierr = pio_put_var(File, cospcnt_desc, (/cosp_cnt(begchunk)/)) - end if - - end subroutine radiation_write_restart - - !=============================================================================== - - subroutine radiation_read_restart(file) - - ! read variables from restart file - - ! arguments - type(file_desc_t), intent(inout) :: file - - ! local variables - - integer :: err_handling - integer :: ierr - real(r8) :: temp_var - - type(var_desc_t) :: vardesc - !---------------------------------------------------------------------------- - - if (docosp) then - call pio_seterrorhandling(File, PIO_BCAST_ERROR, err_handling) - ierr = pio_inq_varid(File, 'cosp_cnt_init', vardesc) - call pio_seterrorhandling(File, err_handling) - if (ierr /= PIO_NOERR) then - cosp_cnt_init = 0 - else - ierr = pio_get_var(File, vardesc, cosp_cnt_init) - end if + else + do i = 1, ncol + coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg) + end do + end if + + ! Gather night/day column indices. + Nday = 0 + Nnite = 0 + do i = 1, ncol + if ( coszrs(i) > 0.0_r8 ) then + Nday = Nday + 1 + IdxDay(Nday) = i + else + Nnite = Nnite + 1 + IdxNite(Nnite) = i end if + end do + + ! Associate pointers to physics buffer fields + itim_old = pbuf_old_tim_idx() + if (cldfsnow_idx > 0) then + call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + endif + if (cldfgrau_idx > 0 .and. graupel_in_rad) then + call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + endif + + call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + + call pbuf_get_field(pbuf, qrs_idx, qrs) + call pbuf_get_field(pbuf, qrl_idx, qrl) + + call pbuf_get_field(pbuf, fsnt_idx, fsnt) + call pbuf_get_field(pbuf, fsds_idx, fsds) + call pbuf_get_field(pbuf, fsns_idx, fsns) + call pbuf_get_field(pbuf, flns_idx, flns) + call pbuf_get_field(pbuf, flnt_idx, flnt) + + if (spectralflux) then + call pbuf_get_field(pbuf, su_idx, su) + call pbuf_get_field(pbuf, sd_idx, sd) + call pbuf_get_field(pbuf, lu_idx, lu) + call pbuf_get_field(pbuf, ld_idx, ld) + end if + + ! For CRM, make cloud equal to input observations: + if (scm_crm_mode .and. have_cld) then + do k = 1, pver + cld(:ncol,k)= cldobs(k) + end do + end if - ierr = pio_inq_varid(File, 'nextsw_cday', vardesc) - ierr = pio_get_var(File, vardesc, temp_var) - nextsw_cday = temp_var - - end subroutine radiation_read_restart - - !=============================================================================== - - subroutine radiation_tend( & - state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out) - - !----------------------------------------------------------------------- - ! - ! Driver for radiation computation. - ! - ! Revision history: - ! 2007-11-05 M. Iacono Install rrtmg_lw and sw as radiation model. - ! 2007-12-27 M. Iacono Modify to use CAM cloud optical properties with rrtmg. - !----------------------------------------------------------------------- - - use phys_grid, only: get_rlat_all_p, get_rlon_all_p - use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr - use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz - - use aer_rad_props, only: aer_rad_props_sw, aer_rad_props_lw - - use cloud_rad_props, only: get_ice_optics_sw, get_liquid_optics_sw, liquid_cloud_get_rad_props_lw, & - ice_cloud_get_rad_props_lw, cloud_rad_props_get_lw, & - grau_cloud_get_rad_props_lw, get_grau_optics_sw, & - snow_cloud_get_rad_props_lw, get_snow_optics_sw - use slingo, only: slingo_liq_get_rad_props_lw, slingo_liq_optics_sw - use ebert_curry, only: ec_ice_optics_sw, ec_ice_get_rad_props_lw - - use rad_solar_var, only: get_variability - use radsw, only: rad_rrtmg_sw - use radlw, only: rad_rrtmg_lw - use radheat, only: radheat_tend - - use radiation_data, only: rad_data_write - use rrtmg_state, only: rrtmg_state_create, rrtmg_state_update, rrtmg_state_destroy, rrtmg_state_t, & - num_rrtmg_levs - - use interpolate_data, only: vertinterp - use tropopause, only: tropopause_find, TROP_ALG_HYBSTOB, TROP_ALG_CLIMATE - - use cospsimulator_intr, only: docosp, cospsimulator_intr_run, cosp_nradsteps -#ifdef OSLO_AERO - use physics_buffer, only: pbuf_get_index, pbuf_get_field - use oslo_aero_control, only: oslo_aero_getopts - use oslo_aero_share -#endif - - ! Arguments - type(physics_state), intent(in), target :: state - type(physics_ptend), intent(out) :: ptend - type(physics_buffer_desc), pointer :: pbuf(:) - type(cam_out_t), intent(inout) :: cam_out - type(cam_in_t), intent(in) :: cam_in - real(r8), intent(out) :: net_flx(pcols) - type(rad_out_t), target, optional, intent(out) :: rd_out - - ! Local variables - real(r8) :: volc_fraction_coarse ! Fraction of volcanic aerosols going to coarse mode - logical :: has_prescribed_volcaero_cmip6 = .false. ! TODO: the following is hardwired for now - real(r8), pointer :: rvolcmmr(:,:) ! stratospheric volcanoes aerosol mmr - real(r8), pointer :: volcopt(:,:) ! Read in stratospheric volcano SW optical parameter (CMIP6) - integer :: band - character(len=3) :: c3 - logical :: idrf - type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object if argument is not present - logical :: write_output - ! - integer :: i, k - integer :: lchnk, ncol - logical :: dosw, dolw - real(r8) :: calday ! current calendar day - real(r8) :: delta ! Solar declination angle in radians - real(r8) :: eccf ! Earth orbit eccentricity factor - real(r8) :: clat(pcols) ! current latitudes(radians) - real(r8) :: clon(pcols) ! current longitudes(radians) - real(r8) :: coszrs(pcols) ! Cosine solar zenith angle - - ! Gathered indices of day and night columns - ! chunk_column_index = IdxDay(daylight_column_index) - integer :: Nday ! Number of daylight columns - integer :: Nnite ! Number of night columns - integer :: IdxDay(pcols) ! Indices of daylight columns - integer :: IdxNite(pcols) ! Indices of night columns - - integer :: itim_old - - real(r8), pointer :: cld(:,:) ! cloud fraction - real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds- whatever they are" - real(r8), pointer :: cldfgrau(:,:) ! cloud fraction of just "snow clouds- whatever they are" - real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate - real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate - real(r8), pointer :: fsds(:) ! Surface solar down flux - real(r8), pointer :: fsns(:) ! Surface solar absorbed flux - real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top - real(r8), pointer :: flns(:) ! Srf longwave cooling (up-down) flux - real(r8), pointer :: flnt(:) ! Net outgoing lw flux at model top - - real(r8), pointer, dimension(:,:,:) :: su => NULL() ! shortwave spectral flux up - real(r8), pointer, dimension(:,:,:) :: sd => NULL() ! shortwave spectral flux down - real(r8), pointer, dimension(:,:,:) :: lu => NULL() ! longwave spectral flux up - real(r8), pointer, dimension(:,:,:) :: ld => NULL() ! longwave spectral flux down - - ! tropopause diagnostic - integer :: troplev(pcols) - real(r8) :: p_trop(pcols) - - type(rrtmg_state_t) :: r_state ! contains the atm concentrations in layers needed for RRTMG - - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: ice_tau (nswbands,pcols,pver) ! ice extinction optical depth - real(r8) :: ice_tau_w (nswbands,pcols,pver) ! ice single scattering albedo * tau - real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! ice assymetry parameter * tau * w - real(r8) :: ice_tau_w_f(nswbands,pcols,pver) ! ice forward scattered fraction * tau * w - real(r8) :: ice_lw_abs (nlwbands,pcols,pver) ! ice absorption optics depth (LW) - - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth - real(r8) :: liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau - real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! liquid assymetry parameter * tau * w - real(r8) :: liq_tau_w_f(nswbands,pcols,pver) ! liquid forward scattered fraction * tau * w - real(r8) :: liq_lw_abs (nlwbands,pcols,pver) ! liquid absorption optics depth (LW) - - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: cld_tau (nswbands,pcols,pver) ! cloud extinction optical depth - real(r8) :: cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau - real(r8) :: cld_tau_w_g(nswbands,pcols,pver) ! cloud assymetry parameter * w * tau - real(r8) :: cld_tau_w_f(nswbands,pcols,pver) ! cloud forward scattered fraction * w * tau - real(r8) :: cld_lw_abs (nlwbands,pcols,pver) ! cloud absorption optics depth (LW) - - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth - real(r8) :: snow_tau_w (nswbands,pcols,pver) ! snow single scattering albedo * tau - real(r8) :: snow_tau_w_g(nswbands,pcols,pver) ! snow assymetry parameter * tau * w - real(r8) :: snow_tau_w_f(nswbands,pcols,pver) ! snow forward scattered fraction * tau * w - real(r8) :: snow_lw_abs (nlwbands,pcols,pver)! snow absorption optics depth (LW) - - ! Add graupel as another snow species. - ! cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: grau_tau (nswbands,pcols,pver) ! graupel extinction optical depth - real(r8) :: grau_tau_w (nswbands,pcols,pver) ! graupel single scattering albedo * tau - real(r8) :: grau_tau_w_g(nswbands,pcols,pver) ! graupel assymetry parameter * tau * w - real(r8) :: grau_tau_w_f(nswbands,pcols,pver) ! graupel forward scattered fraction * tau * w - real(r8) :: grau_lw_abs (nlwbands,pcols,pver)! graupel absorption optics depth (LW) - - ! combined cloud radiative parameters are "in cloud" not "in cell" - real(r8) :: cldfprime(pcols,pver) ! combined cloud fraction (snow plus regular) - real(r8) :: c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth - real(r8) :: c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau - real(r8) :: c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud assymetry parameter * w * tau - real(r8) :: c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau - real(r8) :: c_cld_lw_abs (nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW) - - real(r8) :: sfac(1:nswbands) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band - - integer :: icall ! index through climate/diagnostic radiation calls - logical :: active_calls(0:N_DIAG) - - ! Aerosol radiative properties - real(r8) :: aer_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth - real(r8) :: aer_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau - real(r8) :: aer_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau - real(r8) :: aer_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau - real(r8) :: aer_lw_abs (pcols,pver,nlwbands) ! aerosol absorption optics depth (LW) - -#ifdef OSLO_AERO - ! Local variables used for calculating aerosol optics and direct and indirect forcings. - ! aodvis and absvis are AOD and absorptive AOD for visible wavelength close to 0.55 um (0.35-0.64) - ! Note that aodvis and absvis output should be devided by dayfoc to give physical (A)AOD values - real(r8) :: qdirind(pcols,pver,pcnst) ! Common tracers for indirect and direct calculations - real(r8) :: aodvis(pcols) ! AOD vis - real(r8) :: absvis(pcols) ! absorptive AOD vis - real(r8) :: clearodvis(pcols), clearabsvis(pcols), cloudfree(pcols), cloudfreemax(pcols) - real(r8) :: clearod440(pcols),clearod550(pcols),clearod870(pcols),clearabs550(pcols),clearabs550alt(pcols) - real(r8) :: ftem_1d(pcols) ! work-array to avoid NAN and pcols/ncol confusion - real(r8) :: Nnatk(pcols,pver,0:nmodes_oslo) ! Modal aerosol number concentration - real(r8) :: per_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth - real(r8) :: per_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau - real(r8) :: per_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau - real(r8) :: per_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau - real(r8) :: per_lw_abs (pcols,pver,nlwbands) ! aerosol absorption optics depth (LW) - real(r8) :: volc_ext_sun(pcols,pver,nswbands) ! volcanic aerosol extinction for solar bands, CMIP6 - real(r8) :: volc_omega_sun(pcols,pver,nswbands) ! volcanic aerosol SSA for solar bands, CMIP6 - real(r8) :: volc_g_sun(pcols,pver,nswbands) ! volcanic aerosol g for solar bands, CMIP6 - real(r8) :: volc_ext_earth(pcols,pver,nlwbands) ! volcanic aerosol extinction for terrestrial bands, CMIP6 - real(r8) :: volc_omega_earth(pcols,pver,nlwbands) ! volcanic aerosol SSA for terrestrial bands, CMIP6 - integer :: ns ! spectral loop index -#endif - real(r8) :: fns(pcols,pverp) ! net shortwave flux - real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux - real(r8) :: fnl(pcols,pverp) ! net longwave flux - real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux - - ! for COSP - real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity - real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau - real(r8) :: gb_snow_lw(pcols,pver) ! grid-box mean LW snow optical depth - - real(r8) :: ftem(pcols,pver) ! Temporary workspace for outfld variables - - real(r8) :: freqclr(pcols) ! Frequency of occurrence of clear sky columns - real(r8) :: flntclr(pcols) ! Clearsky only columns (zero if cloudy) - - character(*), parameter :: name = 'radiation_tend' - !-------------------------------------------------------------------------------------- - - lchnk = state%lchnk - ncol = state%ncol - - per_lw_abs(:,:,:) =0._r8 - per_tau(:,:,:) =0._r8 - per_tau_w(:,:,:) =0._r8 - per_tau_w_g(:,:,:) =0._r8 - per_tau_w_f(:,:,:) =0._r8 + ! Find tropopause height if needed for diagnostic output + if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then + call tropopause_find(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE) + endif - if (present(rd_out)) then - rd => rd_out - write_output = .false. - else - allocate(rd) - write_output=.true. - end if + ! Get time of next radiation calculation - albedos will need to be + ! calculated by each surface model at this time + nextsw_cday = radiation_nextsw_cday() - dosw = radiation_do('sw') ! do shortwave heating calc this timestep? - dolw = radiation_do('lw') ! do longwave heating calc this timestep? + if (dosw .or. dolw) then - ! Cosine solar zenith angle for current time step - calday = get_curr_calday() - call get_rlat_all_p(lchnk, ncol, clat) - call get_rlon_all_p(lchnk, ncol, clon) + ! construct an RRTMG state object + r_state = rrtmg_state_create( state, cam_in ) - call shr_orb_decl(calday, eccen, mvelpp, lambm0, obliqr, & - delta, eccf) + call t_startf('cldoptics') - if (use_rad_uniform_angle) then - do i = 1, ncol - coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg, uniform_angle=rad_uniform_angle) + if (cldfsnow_idx > 0) then + do k = 1, pver + do i = 1, ncol + cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k)) + end do end do else - do i = 1, ncol - coszrs(i) = shr_orb_cosz(calday, clat(i), clon(i), delta, dt_avg) - end do + cldfprime(:ncol,:) = cld(:ncol,:) end if - ! Gather night/day column indices. - Nday = 0 - Nnite = 0 - do i = 1, ncol - if ( coszrs(i) > 0.0_r8 ) then - Nday = Nday + 1 - IdxDay(Nday) = i - else - Nnite = Nnite + 1 - IdxNite(Nnite) = i - end if - end do - - ! Associate pointers to physics buffer fields - itim_old = pbuf_old_tim_idx() - if (cldfsnow_idx > 0) then - call pbuf_get_field(pbuf, cldfsnow_idx, cldfsnow, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - endif if (cldfgrau_idx > 0 .and. graupel_in_rad) then - call pbuf_get_field(pbuf, cldfgrau_idx, cldfgrau, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - endif - - call pbuf_get_field(pbuf, cld_idx, cld, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - - call pbuf_get_field(pbuf, qrs_idx, qrs) - call pbuf_get_field(pbuf, qrl_idx, qrl) - - call pbuf_get_field(pbuf, fsnt_idx, fsnt) - call pbuf_get_field(pbuf, fsds_idx, fsds) - call pbuf_get_field(pbuf, fsns_idx, fsns) - call pbuf_get_field(pbuf, flns_idx, flns) - call pbuf_get_field(pbuf, flnt_idx, flnt) - - if (spectralflux) then - call pbuf_get_field(pbuf, su_idx, su) - call pbuf_get_field(pbuf, sd_idx, sd) - call pbuf_get_field(pbuf, lu_idx, lu) - call pbuf_get_field(pbuf, ld_idx, ld) - end if - - ! For CRM, make cloud equal to input observations: - if (scm_crm_mode .and. have_cld) then do k = 1, pver - cld(:ncol,k)= cldobs(k) + do i = 1, ncol + cldfprime(i,k) = max(cld(i,k), cldfgrau(i,k)) + end do end do end if -#ifdef OSLO_AERO - qdirind(:ncol,:,:) = state%q(:ncol,:,:) - if (has_prescribed_volcaero) then - call oslo_aero_getopts(volc_fraction_coarse_out=volc_fraction_coarse) - call pbuf_get_field(pbuf, volc_idx, rvolcmmr, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - qdirind(:ncol,:,l_so4_pr) = qdirind(:ncol,:,l_so4_pr) + (1.0_r8 - volc_fraction_coarse)*rvolcmmr(:ncol,:) - qdirind(:ncol,:,l_ss_a3) = qdirind(:ncol,:,l_ss_a3) + volc_fraction_coarse*rvolcmmr(:ncol,:) - end if -#endif + if (dosw) then + + if (oldcldoptics) then + call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.false.) + call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.false.) + else + select case (icecldoptics) + case ('ebertcurry') + call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.true.) + case ('mitchell') + call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f) + case default + call endrun('iccldoptics must be one either ebertcurry or mitchell') + end select + + select case (liqcldoptics) + case ('slingo') + call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.true.) + case ('gammadist') + call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) + case default + call endrun('liqcldoptics must be either slingo or gammadist') + end select + end if - ! Find tropopause height if needed for diagnostic output - if (hist_fld_active('FSNR') .or. hist_fld_active('FLNR')) then - call tropopause_find(state, troplev, tropP=p_trop, primary=TROP_ALG_HYBSTOB, backup=TROP_ALG_CLIMATE) - endif + cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:) + cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:) + cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:) + cld_tau_w_f(:,:ncol,:) = liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:) - ! Get time of next radiation calculation - albedos will need to be - ! calculated by each surface model at this time - nextsw_cday = radiation_nextsw_cday() + if (cldfsnow_idx > 0) then + ! add in snow + call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, snow_tau_w_f) + do i = 1, ncol + do k = 1, pver - if (dosw .or. dolw) then + if (cldfprime(i,k) > 0._r8) then - ! construct an RRTMG state object - r_state = rrtmg_state_create( state, cam_in ) + c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) & + + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k) - call t_startf('cldoptics') + c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) & + + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k) - if (cldfsnow_idx > 0) then - do k = 1, pver - do i = 1, ncol - cldfprime(i,k) = max(cld(i,k), cldfsnow(i,k)) + c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) & + + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k) + + c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) & + + cld(i,k)*cld_tau_w_f(:,i,k) )/cldfprime(i,k) + else + c_cld_tau(:,i,k) = 0._r8 + c_cld_tau_w(:,i,k) = 0._r8 + c_cld_tau_w_g(:,i,k) = 0._r8 + c_cld_tau_w_f(:,i,k) = 0._r8 + end if end do end do else - cldfprime(:ncol,:) = cld(:ncol,:) + c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:) + c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:) + c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:) + c_cld_tau_w_f(:,:ncol,:) = cld_tau_w_f(:,:ncol,:) end if if (cldfgrau_idx > 0 .and. graupel_in_rad) then - do k = 1, pver - do i = 1, ncol - cldfprime(i,k) = max(cld(i,k), cldfgrau(i,k)) - end do - end do - end if + ! add in graupel + call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, grau_tau_w_f) + do i = 1, ncol + do k = 1, pver - if (dosw) then - - if (oldcldoptics) then - call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.false.) - call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.false.) - else - select case (icecldoptics) - case ('ebertcurry') - call ec_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f, oldicewp=.true.) - case ('mitchell') - call get_ice_optics_sw(state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f) - case default - call endrun('iccldoptics must be one either ebertcurry or mitchell') - end select - - select case (liqcldoptics) - case ('slingo') - call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.true.) - case ('gammadist') - call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) - case default - call endrun('liqcldoptics must be either slingo or gammadist') - end select - end if + if (cldfprime(i,k) > 0._r8) then - cld_tau(:,:ncol,:) = liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:) - cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:) - cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:) - cld_tau_w_f(:,:ncol,:) = liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:) + c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) & + + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k) - if (cldfsnow_idx > 0) then - ! add in snow - call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, snow_tau_w_f) - do i = 1, ncol - do k = 1, pver + c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) & + + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k) - if (cldfprime(i,k) > 0._r8) then + c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) & + + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k) - c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) & - + cld(i,k)*cld_tau(:,i,k) )/cldfprime(i,k) + c_cld_tau_w_f(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_f(:,i,k) & + + cld(i,k)*c_cld_tau_w_f(:,i,k) )/cldfprime(i,k) + else + c_cld_tau(:,i,k) = 0._r8 + c_cld_tau_w(:,i,k) = 0._r8 + c_cld_tau_w_g(:,i,k) = 0._r8 + c_cld_tau_w_f(:,i,k) = 0._r8 + end if + end do + end do + end if - c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) & - + cld(i,k)*cld_tau_w(:,i,k) )/cldfprime(i,k) + ! Output cloud optical depth fields for the visible band + rd%tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:) + rd%liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:) + rd%ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:) - c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) & - + cld(i,k)*cld_tau_w_g(:,i,k) )/cldfprime(i,k) + if (cldfsnow_idx > 0) then + rd%snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:) + endif - c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) & - + cld(i,k)*cld_tau_w_f(:,i,k) )/cldfprime(i,k) - else - c_cld_tau(:,i,k) = 0._r8 - c_cld_tau_w(:,i,k) = 0._r8 - c_cld_tau_w_g(:,i,k) = 0._r8 - c_cld_tau_w_f(:,i,k) = 0._r8 - end if - end do - end do - else - c_cld_tau(:,:ncol,:) = cld_tau(:,:ncol,:) - c_cld_tau_w(:,:ncol,:) = cld_tau_w(:,:ncol,:) - c_cld_tau_w_g(:,:ncol,:) = cld_tau_w_g(:,:ncol,:) - c_cld_tau_w_f(:,:ncol,:) = cld_tau_w_f(:,:ncol,:) + if (cldfgrau_idx > 0 .and. graupel_in_rad) then + rd%grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:) + endif + + ! multiply by total cloud fraction to get gridbox value + rd%tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:) + + ! add fillvalue for night columns + do i = 1, Nnite + rd%tot_cld_vistau(IdxNite(i),:) = fillvalue + rd%tot_icld_vistau(IdxNite(i),:) = fillvalue + rd%liq_icld_vistau(IdxNite(i),:) = fillvalue + rd%ice_icld_vistau(IdxNite(i),:) = fillvalue + if (cldfsnow_idx > 0) then + rd%snow_icld_vistau(IdxNite(i),:) = fillvalue end if - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - ! add in graupel - call get_grau_optics_sw(state, pbuf, grau_tau, grau_tau_w, grau_tau_w_g, grau_tau_w_f) - do i = 1, ncol - do k = 1, pver - - if (cldfprime(i,k) > 0._r8) then - - c_cld_tau(:,i,k) = ( cldfgrau(i,k)*grau_tau(:,i,k) & - + cld(i,k)*c_cld_tau(:,i,k) )/cldfprime(i,k) - - c_cld_tau_w(:,i,k) = ( cldfgrau(i,k)*grau_tau_w(:,i,k) & - + cld(i,k)*c_cld_tau_w(:,i,k) )/cldfprime(i,k) - - c_cld_tau_w_g(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_g(:,i,k) & - + cld(i,k)*c_cld_tau_w_g(:,i,k) )/cldfprime(i,k) - - c_cld_tau_w_f(:,i,k) = ( cldfgrau(i,k)*grau_tau_w_f(:,i,k) & - + cld(i,k)*c_cld_tau_w_f(:,i,k) )/cldfprime(i,k) - else - c_cld_tau(:,i,k) = 0._r8 - c_cld_tau_w(:,i,k) = 0._r8 - c_cld_tau_w_g(:,i,k) = 0._r8 - c_cld_tau_w_f(:,i,k) = 0._r8 - end if - end do - end do + rd%grau_icld_vistau(IdxNite(i),:) = fillvalue end if + end do - ! Output cloud optical depth fields for the visible band - rd%tot_icld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:) - rd%liq_icld_vistau(:ncol,:) = liq_tau(idx_sw_diag,:ncol,:) - rd%ice_icld_vistau(:ncol,:) = ice_tau(idx_sw_diag,:ncol,:) + if (write_output) call radiation_output_cld(lchnk, ncol, rd) - if (cldfsnow_idx > 0) then - rd%snow_icld_vistau(:ncol,:) = snow_tau(idx_sw_diag,:ncol,:) - endif + end if ! if (dosw) - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - rd%grau_icld_vistau(:ncol,:) = grau_tau(idx_sw_diag,:ncol,:) - endif - - ! multiply by total cloud fraction to get gridbox value - rd%tot_cld_vistau(:ncol,:) = c_cld_tau(idx_sw_diag,:ncol,:)*cldfprime(:ncol,:) - - ! add fillvalue for night columns - do i = 1, Nnite - rd%tot_cld_vistau(IdxNite(i),:) = fillvalue - rd%tot_icld_vistau(IdxNite(i),:) = fillvalue - rd%liq_icld_vistau(IdxNite(i),:) = fillvalue - rd%ice_icld_vistau(IdxNite(i),:) = fillvalue - if (cldfsnow_idx > 0) then - rd%snow_icld_vistau(IdxNite(i),:) = fillvalue - end if - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - rd%grau_icld_vistau(IdxNite(i),:) = fillvalue - end if - end do + if (dolw) then - if (write_output) call radiation_output_cld(lchnk, ncol, rd) + if (oldcldoptics) then + call cloud_rad_props_get_lw(state, pbuf, cld_lw_abs, oldcloud=.true.) + else + select case (icecldoptics) + case ('ebertcurry') + call ec_ice_get_rad_props_lw(state, pbuf, ice_lw_abs, oldicewp=.true.) + case ('mitchell') + call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs) + case default + call endrun('iccldoptics must be one either ebertcurry or mitchell') + end select + + select case (liqcldoptics) + case ('slingo') + call slingo_liq_get_rad_props_lw(state, pbuf, liq_lw_abs, oldliqwp=.true.) + case ('gammadist') + call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs) + case default + call endrun('liqcldoptics must be either slingo or gammadist') + end select + + cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:) - end if ! if (dosw) + end if - if (dolw) then + if (cldfsnow_idx > 0) then - if (oldcldoptics) then - call cloud_rad_props_get_lw(state, pbuf, cld_lw_abs, oldcloud=.true.) - else - select case (icecldoptics) - case ('ebertcurry') - call ec_ice_get_rad_props_lw(state, pbuf, ice_lw_abs, oldicewp=.true.) - case ('mitchell') - call ice_cloud_get_rad_props_lw(state, pbuf, ice_lw_abs) - case default - call endrun('iccldoptics must be one either ebertcurry or mitchell') - end select + ! add in snow + call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs) - select case (liqcldoptics) - case ('slingo') - call slingo_liq_get_rad_props_lw(state, pbuf, liq_lw_abs, oldliqwp=.true.) - case ('gammadist') - call liquid_cloud_get_rad_props_lw(state, pbuf, liq_lw_abs) - case default - call endrun('liqcldoptics must be either slingo or gammadist') - end select + do i = 1, ncol + do k = 1, pver + if (cldfprime(i,k) > 0._r8) then + c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) & + + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k) + else + c_cld_lw_abs(:,i,k) = 0._r8 + end if + end do + end do + else + c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:) + end if - cld_lw_abs(:,:ncol,:) = liq_lw_abs(:,:ncol,:) + ice_lw_abs(:,:ncol,:) + if (cldfgrau_idx > 0 .and. graupel_in_rad) then - end if + ! add in graupel + call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs) - if (cldfsnow_idx > 0) then + do i = 1, ncol + do k = 1, pver + if (cldfprime(i,k) > 0._r8) then + c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) & + + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k) + else + c_cld_lw_abs(:,i,k) = 0._r8 + end if + end do + end do + end if - ! add in snow - call snow_cloud_get_rad_props_lw(state, pbuf, snow_lw_abs) + end if ! if (dolw) - do i = 1, ncol - do k = 1, pver - if (cldfprime(i,k) > 0._r8) then - c_cld_lw_abs(:,i,k) = ( cldfsnow(i,k)*snow_lw_abs(:,i,k) & - + cld(i,k)*cld_lw_abs(:,i,k) )/cldfprime(i,k) - else - c_cld_lw_abs(:,i,k) = 0._r8 - end if - end do - end do - else - c_cld_lw_abs(:,:ncol,:) = cld_lw_abs(:,:ncol,:) - end if + call t_stopf('cldoptics') - if (cldfgrau_idx > 0 .and. graupel_in_rad) then + ! Solar radiation computation - ! add in graupel - call grau_cloud_get_rad_props_lw(state, pbuf, grau_lw_abs) + ! OSLO_AERO begin + if (dosw) then - do i = 1, ncol - do k = 1, pver - if (cldfprime(i,k) > 0._r8) then - c_cld_lw_abs(:,i,k) = ( cldfgrau(i,k)*grau_lw_abs(:,i,k) & - + cld(i,k)*c_cld_lw_abs(:,i,k) )/cldfprime(i,k) - else - c_cld_lw_abs(:,i,k) = 0._r8 - end if - end do - end do - end if + qdirind(:ncol,:,:) = state%q(:ncol,:,:) + if (has_prescribed_volcaero) then + call oslo_aero_getopts(volc_fraction_coarse_out = volc_fraction_coarse) + call pbuf_get_field(pbuf, volc_idx, rvolcmmr, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) + qdirind(:ncol,:,l_so4_pr) = qdirind(:ncol,:,l_so4_pr) + (1.0_r8 - volc_fraction_coarse)*rvolcmmr(:ncol,:) + qdirind(:ncol,:,l_ss_a3) = qdirind(:ncol,:,l_ss_a3) + volc_fraction_coarse*rvolcmmr(:ncol,:) + end if - end if ! if (dolw) - - call t_stopf('cldoptics') - - ! ------------------------------------------ - ! Solar radiation computation - ! ------------------------------------------ - - ! OSLO_AERO aerosol - -#ifdef OSLO_AERO - if (dosw) then - - ! Volcanic optics for solar (SW) bands - do band=1,nswbands - volc_ext_sun(1:ncol,1:pver,band)=0.0_r8 - volc_omega_sun(1:ncol,1:pver,band)=0.999_r8 - volc_g_sun(1:ncol,1:pver,band)=0.5_r8 - enddo - if (has_prescribed_volcaero_cmip6) then - do band=1,nlwbands - write(c3,'(i3)') band - volc_idx = pbuf_get_index('ext_sun'//trim(adjustl(c3))) - call pbuf_get_field(pbuf, volc_idx, volcopt, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - volc_ext_sun(1:ncol,1:pver,band)=volcopt(1:ncol,1:pver) - volc_idx = pbuf_get_index('omega_sun'//trim(adjustl(c3))) - call pbuf_get_field(pbuf, volc_idx, volcopt, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - volc_omega_sun(1:ncol,1:pver,band)=volcopt(1:ncol,1:pver) - volc_idx = pbuf_get_index('g_sun'//trim(adjustl(c3))) - call pbuf_get_field(pbuf, volc_idx, volcopt, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - volc_g_sun(1:ncol,1:pver,band)=volcopt(1:ncol,1:pver) - enddo - endif - - ! Volcanic optics for terrestrial (LW) bands (g is not used here) - do band=1,nlwbands - volc_ext_earth(1:ncol,1:pver,band)=0.0_r8 - volc_omega_earth(1:ncol,1:pver,band)=0.999_r8 - enddo - if (has_prescribed_volcaero_cmip6) then - do band=1,nlwbands - write(c3,'(i3)') band - volc_idx = pbuf_get_index('ext_earth'//trim(adjustl(c3))) - call pbuf_get_field(pbuf, volc_idx, volcopt, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - volc_ext_earth(1:ncol,1:pver,band)=volcopt(1:ncol,1:pver) - - volc_idx = pbuf_get_index('omega_earth'//trim(adjustl(c3))) - call pbuf_get_field(pbuf, volc_idx, volcopt, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) - volc_omega_earth(1:ncol,1:pver,band)=volcopt(1:ncol,1:pver) - enddo - endif - - ! No aerocom variables passed for now - ! dod440, dod550, dod870, abs550, abs550alt - call oslo_aero_optical_params_calc(lchnk, ncol, 10.0_r8*state%pint, state%pmid, & - coszrs, state, state%t, cld, qdirind, Nnatk, & - per_tau, per_tau_w, per_tau_w_g, per_tau_w_f, per_lw_abs, & - volc_ext_sun, volc_omega_sun, volc_g_sun, volc_ext_earth, volc_omega_earth, & - aodvis, absvis) - - call get_variability(sfac) - - ! Get the active climate/diagnostic shortwave calculations - call rad_cnst_get_call_list(active_calls) - - ! The climate (icall==0) calculation must occur last. - do icall = N_DIAG, 0, -1 - - if (active_calls(icall)) then - - ! update the concentrations in the RRTMG state object - call rrtmg_state_update(state, pbuf, icall, r_state) - - ! A first call with Oslo aerosols set to zero for radiative forcing diagnostics - ! follwoing the Ghan (2013) method: - ! for calculation of direct radiative forcing, not necessarily "offline" as such anymore - ! (just nudged), but with an extra call with 0 aerosol extiction. - ! - idrf = .true. - call rad_rrtmg_sw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - cldfprime, aer_tau*0.0_r8, aer_tau_w, aer_tau_w_g, aer_tau_w_f, & - eccf, coszrs, rd%solin, sfac, cam_in%asdir, & - cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, & - fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, & - rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & - rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & - cam_out%solld, fns, fcns, Nday, Nnite, & - IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & - E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & - E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false., idrf=idrf) - - ! - ! Dump shortwave radiation information to history tape buffer (diagnostics) - ! Note that DRF fields are now from the aer_tau=0 call (clean), no longer with - ! aer_tau from oslo_aero_optical_params_calc - - ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair - call outfld('QRS_DRF ',ftem ,pcols,lchnk) - - ftem(:ncol,:pver) = rd%qrsc(:ncol,:pver)/cpair - call outfld('QRSC_DRF',ftem ,pcols,lchnk) - - call outfld('FSNT_DRF',fsnt , pcols, lchnk) - call outfld('FSNS_DRF',fsns , pcols, lchnk) - call outfld('FSNTCDRF',rd%fsntc , pcols, lchnk) - call outfld('FSNSCDRF',rd%fsnsc , pcols, lchnk) + ! Volcanic optics for solar (SW) bands + do band=1,nswbands + volc_ext_sun(1:ncol,1:pver,band)=0.0_r8 + volc_omega_sun(1:ncol,1:pver,band)=0.999_r8 + volc_g_sun(1:ncol,1:pver,band)=0.5_r8 + enddo + ! Volcanic optics for terrestrial (LW) bands + do band=1,nlwbands + volc_ext_earth(1:ncol,1:pver,band) = 0.0_r8 + volc_omega_earth(1:ncol,1:pver,band) = 0.999_r8 + enddo + + ! No aerocom variables passed for now + ! dod440, dod550, dod870, abs550, abs550alt + call oslo_aero_optical_params_calc(lchnk, ncol, 10.0_r8*state%pint, state%pmid, & + coszrs, state, state%t, cld, qdirind, & + aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, aer_lw_abs, & + volc_ext_sun, volc_omega_sun, volc_g_sun, volc_ext_earth, volc_omega_earth, & + aodvis, absvis) + + call get_variability(sfac) + + ! Get the active climate/diagnostic shortwave calculations + call rad_cnst_get_call_list(active_calls) + + ! The climate (icall==0) calculation must occur last. + do icall = N_DIAG, 0, -1 + + if (active_calls(icall)) then + + ! update the concentrations in the RRTMG state object + call rrtmg_state_update(state, pbuf, icall, r_state) + + ! A first call with Oslo aerosols set to zero for radiative forcing diagnostics + ! follwoing the Ghan (2013) method: + ! for calculation of direct radiative forcing, not necessarily "offline" as such anymore + ! (just nudged), but with an extra call with 0 aerosol extiction. + + call rad_rrtmg_sw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + cldfprime, aer_tau*0.0_r8, aer_tau_w, aer_tau_w_g, aer_tau_w_f, & + eccf, coszrs, rd%solin, sfac, cam_in%asdir, & + cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, & + fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, & + rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & + rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & + cam_out%solld, fns, fcns, Nday, Nnite, & + IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & + E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & + E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false., idrf=.true.) + + ! Dump shortwave radiation information to history tape buffer (diagnostics) + ! Note that DRF fields are now from the aer_tau=0 call (clean), no longer with + ! aer_tau from oslo_aero_optical_params_calc + + ftem(:ncol,:pver) = qrs(:ncol,:pver)/cpair + call outfld('QRS_DRF ',ftem ,pcols,lchnk) + + ftem(:ncol,:pver) = rd%qrsc(:ncol,:pver)/cpair + call outfld('QRSC_DRF',ftem ,pcols,lchnk) + + call outfld('FSNT_DRF',fsnt , pcols, lchnk) + 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) + 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 - idrf = .false. - call rad_rrtmg_sw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - cldfprime, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, & - eccf, coszrs, rd%solin, sfac, cam_in%asdir, & - cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, & - fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, & - rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & - rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & - cam_out%solld, fns, fcns, Nday, Nnite, & - IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & - E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & - E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false., idrf=idrf) - - ! Output net fluxes at 200 mb - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c) - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200) - if (hist_fld_active('FSNR')) then - do i = 1,ncol - call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i)) - end do - end if - - if (write_output) call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out) + call rad_rrtmg_sw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + cldfprime, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, & + eccf, coszrs, rd%solin, sfac, cam_in%asdir, & + cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, & + fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, & + rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & + rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & + cam_out%solld, fns, fcns, Nday, Nnite, & + IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & + E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & + E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false., idrf=.false.) + + ! Output net fluxes at 200 mb + call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c) + call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200) + if (hist_fld_active('FSNR')) then + do i = 1,ncol + call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i)) + end do end if - end do - end if - ! Calculate cloud-free fraction assuming random overlap - ! (kind of duplicated from cloud_cover_diags::cldsav) - ! (note this duplicated code and may not be consistent with cldtot calculated elsewhere) - cloudfree(1:ncol) = 1.0_r8 - cloudfreemax(1:ncol) = 1.0_r8 - do k = 1, pver - do i=1,ncol - cloudfree(i) = cloudfree(i) * cloudfreemax(i) - cloudfreemax(i) = min(cloudfreemax(i),1.0_r8-cld(i,k)) - end do - end do + if (write_output) call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out) - ! Calculate AOD (visible) for cloud free - do i = 1, ncol - clearodvis(i) = cloudfree(i)*aodvis(i) - clearabsvis(i) = cloudfree(i)*absvis(i) + end if end do + end if - ! clear-sky AOD and absorptive AOD for visible wavelength close to 0.55 um (0.35-0.64) - ! Note that caodvis and cabsvis output should be devided by dayfoc*cloudfree to give physical (A)AOD values - call outfld('CAODVIS ',clearodvis ,pcols,lchnk) - call outfld('CABSVIS ',clearabsvis ,pcols,lchnk) - call outfld('CLDFREE ',cloudfree ,pcols,lchnk) -#ifdef AEROCOM - do i = 1, ncol - clearod440(i) = cloudfree(i)*dod440(i) - clearod550(i) = cloudfree(i)*dod550(i) - clearod870(i) = cloudfree(i)*dod870(i) - clearabs550(i) = cloudfree(i)*abs550(i) - clearabs550alt(i) = cloudfree(i)*abs550alt(i) + ! Calculate cloud-free fraction assuming random overlap + ! (kind of duplicated from cloud_cover_diags::cldsav) + ! (note this duplicated code and may not be consistent with cldtot calculated elsewhere) + cloudfree(1:ncol) = 1.0_r8 + cloudfreemax(1:ncol) = 1.0_r8 + do k = 1, pver + do i=1,ncol + cloudfree(i) = cloudfree(i) * cloudfreemax(i) + cloudfreemax(i) = min(cloudfreemax(i),1.0_r8-cld(i,k)) end do - call outfld('CDOD440 ',clearod440 ,pcols,lchnk) - call outfld('CDOD550 ',clearod550 ,pcols,lchnk) - call outfld('CDOD870 ',clearod870 ,pcols,lchnk) - call outfld('CABS550 ',clearabs550 ,pcols,lchnk) - call outfld('CABS550A',clearabs550alt ,pcols,lchnk) -#endif -#endif - - ! CAM Aerosol - -#ifndef OSLO_AERO - if (dosw) then - - call get_variability(sfac) - - ! Get the active climate/diagnostic shortwave calculations - call rad_cnst_get_call_list(active_calls) - - ! The climate (icall==0) calculation must occur last. - do icall = N_DIAG, 0, -1 - - if (active_calls(icall)) then - - ! update the concentrations in the RRTMG state object - call rrtmg_state_update(state, pbuf, icall, r_state) - - call aer_rad_props_sw(icall, state, pbuf, nnite, idxnite, & - aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f) - - rd%cld_tau_cloudsim(:ncol,:) = cld_tau(rrtmg_sw_cloudsim_band,:ncol,:) - rd%aer_tau550(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag) - rd%aer_tau400(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag+1) - rd%aer_tau700(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag-1) - - call rad_rrtmg_sw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - cldfprime, aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f, & - eccf, coszrs, rd%solin, sfac, cam_in%asdir, & - cam_in%asdif, cam_in%aldir, cam_in%aldif, qrs, rd%qrsc, & - fsnt, rd%fsntc, rd%fsntoa, rd%fsutoa, rd%fsntoac, & - rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & - rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & - cam_out%solld, fns, fcns, Nday, Nnite, & - IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & - E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & - E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false.) - - ! Output net fluxes at 200 mb - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcns, rd%fsn200c) - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fns, rd%fsn200) - if (hist_fld_active('FSNR')) then - do i = 1,ncol - call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fns(i,:), rd%fsnr(i)) - end do - end if + end do - if (write_output) call radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out) + ! Calculate AOD (visible) for cloud free + do i = 1, ncol + clearodvis(i) = cloudfree(i)*aodvis(i) + clearabsvis(i) = cloudfree(i)*absvis(i) + end do - end if - end do - end if ! end of (if do_sw) + ! clear-sky AOD and absorptive AOD for visible wavelength close to 0.55 um (0.35-0.64) + ! Note that caodvis and cabsvis output should be devided by dayfoc*cloudfree to give physical (A)AOD values + call outfld('CAODVIS ',clearodvis ,pcols,lchnk) + call outfld('CABSVIS ',clearabsvis ,pcols,lchnk) + call outfld('CLDFREE ',cloudfree ,pcols,lchnk) +#ifdef AEROCOM + do i = 1, ncol + clearod440(i) = cloudfree(i)*dod440(i) + clearod550(i) = cloudfree(i)*dod550(i) + clearod870(i) = cloudfree(i)*dod870(i) + clearabs550(i) = cloudfree(i)*abs550(i) + clearabs550alt(i) = cloudfree(i)*abs550alt(i) + end do + call outfld('CDOD440 ',clearod440 ,pcols,lchnk) + call outfld('CDOD550 ',clearod550 ,pcols,lchnk) + call outfld('CDOD870 ',clearod870 ,pcols,lchnk) + call outfld('CABS550 ',clearabs550 ,pcols,lchnk) + call outfld('CABS550A',clearabs550alt ,pcols,lchnk) #endif - ! Output aerosol mmr - call rad_cnst_out(0, state, pbuf) - - ! ------------------------------------------ - ! Longwave radiation computation - ! ------------------------------------------ + ! Output aerosol mmr + call rad_cnst_out(0, state, pbuf) -#ifdef OSLO_AERO - if (dolw) then + ! Longwave radiation computation - call rad_cnst_get_call_list(active_calls) + if (dolw) then - ! The climate (icall==0) calculation must occur last. - do icall = N_DIAG, 0, -1 + call rad_cnst_get_call_list(active_calls) - if (active_calls(icall)) then + ! The climate (icall==0) calculation must occur last. + do icall = N_DIAG, 0, -1 - ! update the conctrations in the RRTMG state object - call rrtmg_state_update( state, pbuf, icall, r_state) + if (active_calls(icall)) then - call rad_rrtmg_lw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - aer_lw_abs*0.0_r8, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & - flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, & - rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, & - lu, ld) + ! update the concentrations in the RRTMG state object + call rrtmg_state_update( state, pbuf, icall, r_state) - call outfld('FLNT_DRF',flnt(:) , pcols, lchnk) - call outfld('FLNTCDRF',rd%flntc(:), pcols, lchnk) + call rad_rrtmg_lw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + aer_lw_abs*0.0_r8, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & + flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, & + rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, & + lu, ld) - call rad_rrtmg_lw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & - flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, & - rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, & - lu, ld) + call outfld('FLNT_DRF',flnt(:) , pcols, lchnk) + call outfld('FLNTCDRF',rd%flntc(:), pcols, lchnk) - ! FLNT_ORG is just for temporary testing vs. FLNT - ftem_1d(1:ncol) = cam_out%flwds(1:ncol) - flns(1:ncol) - call outfld('FLUS ',ftem_1d, pcols, lchnk) + call rad_rrtmg_lw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & + flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, & + rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, & + lu, ld) - ! Output fluxes at 200 mb - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200) - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c) - if (hist_fld_active('FLNR')) then - do i = 1,ncol - call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i)) - end do - end if + ! FLNT_ORG is just for temporary testing vs. FLNT + ftem_1d(1:ncol) = cam_out%flwds(1:ncol) - flns(1:ncol) + call outfld('FLUS ',ftem_1d ,pcols,lchnk) - flntclr(:) = 0._r8 - freqclr(:) = 0._r8 - do i = 1, ncol - if (maxval(cldfprime(i,:)) <= 0.1_r8) then - freqclr(i) = 1._r8 - flntclr(i) = rd%flntc(i) - end if + ! Output fluxes at 200 mb + call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200) + call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c) + if (hist_fld_active('FLNR')) then + do i = 1,ncol + call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i)) end do - - if (write_output) call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr) end if - end do - end if -#endif - -#ifndef OSLO_AERO - if (dolw) then - - call rad_cnst_get_call_list(active_calls) - - ! The climate (icall==0) calculation must occur last. - do icall = N_DIAG, 0, -1 - - if (active_calls(icall)) then - ! update the conctrations in the RRTMG state object - call rrtmg_state_update( state, pbuf, icall, r_state) - - call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs) - - call rad_rrtmg_lw( & - lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & - aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & - flns, flnt, rd%flnsc, rd%flntc, cam_out%flwds, & - rd%flut, rd%flutc, fnl, fcnl, rd%fldsc, & - lu, ld) - - ! Output fluxes at 200 mb - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fnl, rd%fln200) - call vertinterp(ncol, pcols, pverp, state%pint, 20000._r8, fcnl, rd%fln200c) - if (hist_fld_active('FLNR')) then - do i = 1,ncol - call vertinterp(1, 1, pverp, state%pint(i,:), p_trop(i), fnl(i,:), rd%flnr(i)) - end do + flntclr(:) = 0._r8 + freqclr(:) = 0._r8 + do i = 1, ncol + if (maxval(cldfprime(i,:)) <= 0.1_r8) then + freqclr(i) = 1._r8 + flntclr(i) = rd%flntc(i) end if + end do - flntclr(:) = 0._r8 - freqclr(:) = 0._r8 - do i = 1, ncol - if (maxval(cldfprime(i,:)) <= 0.1_r8) then - freqclr(i) = 1._r8 - flntclr(i) = rd%flntc(i) - end if - end do + if (write_output) call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr) - if (write_output) call radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr) + end if + end do - end if - end do - end if -#endif + end if + ! OSLO_AERO end - ! deconstruct the RRTMG state object - call rrtmg_state_destroy(r_state) + ! deconstruct the RRTMG state object + call rrtmg_state_destroy(r_state) - if (docosp) then + if (docosp) then - ! initialize and calculate emis - emis(:,:) = 0._r8 - emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(rrtmg_lw_cloudsim_band,:ncol,:)) - call outfld('EMIS', emis, pcols, lchnk) + ! initialize and calculate emis + emis(:,:) = 0._r8 + emis(:ncol,:) = 1._r8 - exp(-cld_lw_abs(rrtmg_lw_cloudsim_band,:ncol,:)) + call outfld('EMIS', emis, pcols, lchnk) - ! compute grid-box mean SW and LW snow optical depth for use by COSP - gb_snow_tau(:,:) = 0._r8 - gb_snow_lw(:,:) = 0._r8 - if (cldfsnow_idx > 0) then - do i = 1, ncol - do k = 1, pver - if (cldfsnow(i,k) > 0._r8) then - - ! Add graupel to snow tau for cosp - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) + & - grau_tau(rrtmg_sw_cloudsim_band,i,k)*cldfgrau(i,k) - gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) + & - grau_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfgrau(i,k) - else - gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) - gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) - end if + ! compute grid-box mean SW and LW snow optical depth for use by COSP + gb_snow_tau(:,:) = 0._r8 + gb_snow_lw(:,:) = 0._r8 + if (cldfsnow_idx > 0) then + do i = 1, ncol + do k = 1, pver + if (cldfsnow(i,k) > 0._r8) then + + ! Add graupel to snow tau for cosp + if (cldfgrau_idx > 0 .and. graupel_in_rad) then + gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) + & + grau_tau(rrtmg_sw_cloudsim_band,i,k)*cldfgrau(i,k) + gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) + & + grau_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfgrau(i,k) + else + gb_snow_tau(i,k) = snow_tau(rrtmg_sw_cloudsim_band,i,k)*cldfsnow(i,k) + gb_snow_lw(i,k) = snow_lw_abs(rrtmg_lw_cloudsim_band,i,k)*cldfsnow(i,k) end if - end do + end if end do - end if + end do + end if - ! advance counter for this timestep (chunk dimension required for thread safety) - cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1 + ! advance counter for this timestep (chunk dimension required for thread safety) + cosp_cnt(lchnk) = cosp_cnt(lchnk) + 1 - ! if counter is the same as cosp_nradsteps, run cosp and reset counter - if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then + ! if counter is the same as cosp_nradsteps, run cosp and reset counter + if (cosp_nradsteps .eq. cosp_cnt(lchnk)) then - ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave - ! optical depths are passed. - call cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs, & - cld_swtau_in=cld_tau(rrtmg_sw_cloudsim_band,:,:),& - snow_tau_in=gb_snow_tau, snow_emis_in=gb_snow_lw) - cosp_cnt(lchnk) = 0 - end if - end if ! end of if (docosp) + ! N.B.: For snow optical properties, the GRID-BOX MEAN shortwave and longwave + ! optical depths are passed. + call cospsimulator_intr_run(state, pbuf, cam_in, emis, coszrs, & + cld_swtau_in=cld_tau(rrtmg_sw_cloudsim_band,:,:),& + snow_tau_in=gb_snow_tau, snow_emis_in=gb_snow_lw) + cosp_cnt(lchnk) = 0 + end if + end if - else ! if (dosw .or. dolw) then + else ! if (dosw .or. dolw) then - ! convert radiative heating rates from Q*dp to Q for energy conservation - do k =1 , pver - do i = 1, ncol - qrs(i,k) = qrs(i,k)/state%pdel(i,k) - qrl(i,k) = qrl(i,k)/state%pdel(i,k) - end do + ! convert radiative heating rates from Q*dp to Q for energy conservation + do k =1 , pver + do i = 1, ncol + qrs(i,k) = qrs(i,k)/state%pdel(i,k) + qrl(i,k) = qrl(i,k)/state%pdel(i,k) end do + end do - end if ! if (dosw .or. dolw) then - - ! output rad inputs and resulting heating rates - call rad_data_write( pbuf, state, cam_in, coszrs ) + end if ! if (dosw .or. dolw) then - ! Compute net radiative heating tendency - call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & - fsnt, flns, flnt, cam_in%asdir, net_flx) + ! output rad inputs and resulting heating rates + call rad_data_write( pbuf, state, cam_in, coszrs ) - if (write_output) then - ! Compute heating rate for dtheta/dt - do k = 1, pver - do i = 1, ncol - ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa - end do - end do - call outfld('HR', ftem, pcols, lchnk) - end if + ! Compute net radiative heating tendency + call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & + fsnt, flns, flnt, cam_in%asdir, net_flx) - ! convert radiative heating rates to Q*dp for energy conservation + if (write_output) then + ! Compute heating rate for dtheta/dt do k = 1, pver do i = 1, ncol - qrs(i,k) = qrs(i,k)*state%pdel(i,k) - qrl(i,k) = qrl(i,k)*state%pdel(i,k) + ftem(i,k) = (qrs(i,k) + qrl(i,k))/cpair * (1.e5_r8/state%pmid(i,k))**cappa end do end do + call outfld('HR', ftem, pcols, lchnk) + end if - cam_out%netsw(:ncol) = fsns(:ncol) + ! convert radiative heating rates to Q*dp for energy conservation + do k = 1, pver + do i = 1, ncol + qrs(i,k) = qrs(i,k)*state%pdel(i,k) + qrl(i,k) = qrl(i,k)*state%pdel(i,k) + end do + end do - if (.not. present(rd_out)) then - deallocate(rd) - end if + cam_out%netsw(:ncol) = fsns(:ncol) - end subroutine radiation_tend + if (.not. present(rd_out)) then + deallocate(rd) + end if - !=============================================================================== +end subroutine radiation_tend - subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out) +!=============================================================================== - ! Dump shortwave radiation information to history buffer. +subroutine radiation_output_sw(lchnk, ncol, icall, rd, pbuf, cam_out) - integer , intent(in) :: lchnk - integer, intent(in) :: ncol - integer, intent(in) :: icall - type(rad_out_t), intent(in) :: rd - type(physics_buffer_desc), pointer :: pbuf(:) - type(cam_out_t), intent(in) :: cam_out + ! Dump shortwave radiation information to history buffer. - ! local variables - real(r8), pointer :: qrs(:,:) - real(r8), pointer :: fsnt(:) - real(r8), pointer :: fsns(:) - real(r8), pointer :: fsds(:) + integer , intent(in) :: lchnk + integer, intent(in) :: ncol + integer, intent(in) :: icall + type(rad_out_t), intent(in) :: rd + type(physics_buffer_desc), pointer :: pbuf(:) + type(cam_out_t), intent(in) :: cam_out - real(r8) :: ftem(pcols) - !---------------------------------------------------------------------------- + ! local variables + real(r8), pointer :: qrs(:,:) + real(r8), pointer :: fsnt(:) + real(r8), pointer :: fsns(:) + real(r8), pointer :: fsds(:) - call pbuf_get_field(pbuf, qrs_idx, qrs) - call pbuf_get_field(pbuf, fsnt_idx, fsnt) - call pbuf_get_field(pbuf, fsns_idx, fsns) - call pbuf_get_field(pbuf, fsds_idx, fsds) + real(r8) :: ftem(pcols) + !---------------------------------------------------------------------------- - call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk) + call pbuf_get_field(pbuf, qrs_idx, qrs) + call pbuf_get_field(pbuf, fsnt_idx, fsnt) + call pbuf_get_field(pbuf, fsns_idx, fsns) + call pbuf_get_field(pbuf, fsds_idx, fsds) - call outfld('QRS'//diag(icall), qrs/cpair, pcols, lchnk) - call outfld('QRSC'//diag(icall), rd%qrsc/cpair, pcols, lchnk) + call outfld('SOLIN'//diag(icall), rd%solin, pcols, lchnk) - call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk) - call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk) - call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk) - call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk) + call outfld('QRS'//diag(icall), qrs(:ncol,:)/cpair, ncol, lchnk) + call outfld('QRSC'//diag(icall), rd%qrsc(:ncol,:)/cpair, ncol, lchnk) - ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol) - call outfld('SWCF'//diag(icall), ftem, pcols, lchnk) + call outfld('FSNT'//diag(icall), fsnt, pcols, lchnk) + call outfld('FSNTC'//diag(icall), rd%fsntc, pcols, lchnk) + call outfld('FSNTOA'//diag(icall), rd%fsntoa, pcols, lchnk) + call outfld('FSNTOAC'//diag(icall), rd%fsntoac, pcols, lchnk) - call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk) + ftem(:ncol) = rd%fsntoa(:ncol) - rd%fsntoac(:ncol) + call outfld('SWCF'//diag(icall), ftem, pcols, lchnk) - call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk) - call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk) - call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk) + call outfld('FSUTOA'//diag(icall), rd%fsutoa, pcols, lchnk) - call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk) - call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk) + call outfld('FSNIRTOA'//diag(icall), rd%fsnirt, pcols, lchnk) + call outfld('FSNRTOAC'//diag(icall), rd%fsnrtc, pcols, lchnk) + call outfld('FSNRTOAS'//diag(icall), rd%fsnirtsq, pcols, lchnk) - call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk) + call outfld('FSN200'//diag(icall), rd%fsn200, pcols, lchnk) + call outfld('FSN200C'//diag(icall), rd%fsn200c, pcols, lchnk) - call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk) - call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk) - call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk) - call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk) + call outfld('FSNR'//diag(icall), rd%fsnr, pcols, lchnk) - call outfld('FSNS'//diag(icall), fsns, pcols, lchnk) - call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk) + call outfld('SOLS'//diag(icall), cam_out%sols, pcols, lchnk) + call outfld('SOLL'//diag(icall), cam_out%soll, pcols, lchnk) + call outfld('SOLSD'//diag(icall), cam_out%solsd, pcols, lchnk) + call outfld('SOLLD'//diag(icall), cam_out%solld, pcols, lchnk) - call outfld('FSDS'//diag(icall), fsds, pcols, lchnk) - call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk) + call outfld('FSNS'//diag(icall), fsns, pcols, lchnk) + call outfld('FSNSC'//diag(icall), rd%fsnsc, pcols, lchnk) - end subroutine radiation_output_sw + call outfld('FSDS'//diag(icall), fsds, pcols, lchnk) + call outfld('FSDSC'//diag(icall), rd%fsdsc, pcols, lchnk) +end subroutine radiation_output_sw - !=============================================================================== - subroutine radiation_output_cld(lchnk, ncol, rd) +!=============================================================================== - ! Dump shortwave cloud optics information to history buffer. +subroutine radiation_output_cld(lchnk, ncol, rd) - integer , intent(in) :: lchnk - integer, intent(in) :: ncol - type(rad_out_t), intent(in) :: rd - !---------------------------------------------------------------------------- + ! Dump shortwave cloud optics information to history buffer. - call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk) - call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk) - call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk) - call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk) - if (cldfsnow_idx > 0) then - call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk) - endif - if (cldfgrau_idx > 0 .and. graupel_in_rad) then - call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau, pcols, lchnk) - endif - end subroutine radiation_output_cld + integer , intent(in) :: lchnk + integer, intent(in) :: ncol + type(rad_out_t), intent(in) :: rd + !---------------------------------------------------------------------------- - !=============================================================================== + call outfld('TOT_CLD_VISTAU', rd%tot_cld_vistau, pcols, lchnk) + call outfld('TOT_ICLD_VISTAU', rd%tot_icld_vistau, pcols, lchnk) + call outfld('LIQ_ICLD_VISTAU', rd%liq_icld_vistau, pcols, lchnk) + call outfld('ICE_ICLD_VISTAU', rd%ice_icld_vistau, pcols, lchnk) + if (cldfsnow_idx > 0) then + call outfld('SNOW_ICLD_VISTAU', rd%snow_icld_vistau, pcols, lchnk) + endif + if (cldfgrau_idx > 0 .and. graupel_in_rad) then + call outfld('GRAU_ICLD_VISTAU', rd%grau_icld_vistau, pcols, lchnk) + endif +end subroutine radiation_output_cld - subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr) +!=============================================================================== - ! Dump longwave radiation information to history buffer +subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, flntclr) - integer, intent(in) :: lchnk - integer, intent(in) :: ncol - integer, intent(in) :: icall ! icall=0 for climate diagnostics - type(rad_out_t), intent(in) :: rd - type(physics_buffer_desc), pointer :: pbuf(:) - type(cam_out_t), intent(in) :: cam_out - real(r8), intent(in) :: freqclr(pcols) - real(r8), intent(in) :: flntclr(pcols) + ! Dump longwave radiation information to history buffer - ! local variables - real(r8), pointer :: qrl(:,:) - real(r8), pointer :: flnt(:) - real(r8), pointer :: flns(:) + integer, intent(in) :: lchnk + integer, intent(in) :: ncol + integer, intent(in) :: icall ! icall=0 for climate diagnostics + type(rad_out_t), intent(in) :: rd + type(physics_buffer_desc), pointer :: pbuf(:) + type(cam_out_t), intent(in) :: cam_out + real(r8), intent(in) :: freqclr(pcols) + real(r8), intent(in) :: flntclr(pcols) - real(r8) :: ftem(pcols) - !---------------------------------------------------------------------------- + ! local variables + real(r8), pointer :: qrl(:,:) + real(r8), pointer :: flnt(:) + real(r8), pointer :: flns(:) - call pbuf_get_field(pbuf, qrl_idx, qrl) - call pbuf_get_field(pbuf, flnt_idx, flnt) - call pbuf_get_field(pbuf, flns_idx, flns) + real(r8) :: ftem(pcols) + !---------------------------------------------------------------------------- - call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk) - call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk) + call pbuf_get_field(pbuf, qrl_idx, qrl) + call pbuf_get_field(pbuf, flnt_idx, flnt) + call pbuf_get_field(pbuf, flns_idx, flns) - call outfld('FLNT'//diag(icall), flnt, pcols, lchnk) - call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk) + call outfld('QRL'//diag(icall), qrl(:ncol,:)/cpair, ncol, lchnk) + call outfld('QRLC'//diag(icall), rd%qrlc(:ncol,:)/cpair, ncol, lchnk) - call outfld('FREQCLR'//diag(icall), freqclr, pcols, lchnk) - call outfld('FLNTCLR'//diag(icall), flntclr, pcols, lchnk) + call outfld('FLNT'//diag(icall), flnt, pcols, lchnk) + call outfld('FLNTC'//diag(icall), rd%flntc, pcols, lchnk) - call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk) - call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk) + call outfld('FREQCLR'//diag(icall), freqclr, pcols, lchnk) + call outfld('FLNTCLR'//diag(icall), flntclr, pcols, lchnk) - ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol) - call outfld('LWCF'//diag(icall), ftem, pcols, lchnk) + call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk) + call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk) - call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk) - call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk) + ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol) + call outfld('LWCF'//diag(icall), ftem, pcols, lchnk) - call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk) + call outfld('FLN200'//diag(icall), rd%fln200, pcols, lchnk) + call outfld('FLN200C'//diag(icall), rd%fln200c, pcols, lchnk) - call outfld('FLNS'//diag(icall), flns, pcols, lchnk) - call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk) + call outfld('FLNR'//diag(icall), rd%flnr, pcols, lchnk) - call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk) - call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk) + call outfld('FLNS'//diag(icall), flns, pcols, lchnk) + call outfld('FLNSC'//diag(icall), rd%flnsc, pcols, lchnk) - end subroutine radiation_output_lw + call outfld('FLDS'//diag(icall), cam_out%flwds, pcols, lchnk) + call outfld('FLDSC'//diag(icall), rd%fldsc, pcols, lchnk) - !=============================================================================== +end subroutine radiation_output_lw - subroutine calc_col_mean(state, mmr_pointer, mean_value) +!=============================================================================== - ! Compute the column mean mass mixing ratio. +subroutine calc_col_mean(state, mmr_pointer, mean_value) - type(physics_state), intent(in) :: state - real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev) - real(r8), dimension(pcols), intent(out) :: mean_value ! column mean mmr + ! Compute the column mean mass mixing ratio. - integer :: i, k, ncol - real(r8) :: ptot(pcols) - !----------------------------------------------------------------------- + type(physics_state), intent(in) :: state + real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev) + real(r8), dimension(pcols), intent(out) :: mean_value ! column mean mmr - ncol = state%ncol - mean_value = 0.0_r8 - ptot = 0.0_r8 + integer :: i, k, ncol + real(r8) :: ptot(pcols) + !----------------------------------------------------------------------- - do k=1,pver - do i=1,ncol - mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k) - ptot(i) = ptot(i) + state%pdeldry(i,k) - end do - end do + ncol = state%ncol + mean_value = 0.0_r8 + ptot = 0.0_r8 + + do k=1,pver do i=1,ncol - mean_value(i) = mean_value(i) / ptot(i) + mean_value(i) = mean_value(i) + mmr_pointer(i,k)*state%pdeldry(i,k) + ptot(i) = ptot(i) + state%pdeldry(i,k) end do + end do + do i=1,ncol + mean_value(i) = mean_value(i) / ptot(i) + end do - end subroutine calc_col_mean +end subroutine calc_col_mean - !=============================================================================== +!=============================================================================== end module radiation diff --git a/src_cam/radlw.F90 b/src_cam/radlw.F90 index 7e2d80d..6864039 100644 --- a/src_cam/radlw.F90 +++ b/src_cam/radlw.F90 @@ -1,7 +1,7 @@ module radlw -!----------------------------------------------------------------------- -! +!----------------------------------------------------------------------- +! ! Purpose: Longwave radiation calculations. ! !----------------------------------------------------------------------- @@ -27,14 +27,14 @@ module radlw public ::& radlw_init, &! initialize constants rad_rrtmg_lw ! driver for longwave radiation code - + ! Private data integer :: ntoplw ! top level to solve for longwave cooling ! Flag for cloud overlap method ! 0=clear, 1=random, 2=maximum/random, 3=maximum integer, parameter :: icld = 2 - + !=============================================================================== CONTAINS @@ -91,7 +91,7 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & real(r8), pointer, dimension(:,:,:) :: lu ! longwave spectral flux up real(r8), pointer, dimension(:,:,:) :: ld ! longwave spectral flux down - + ! !---------------------------Local variables----------------------------- ! @@ -109,7 +109,7 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & real(r8), parameter :: dps = 1._r8/86400._r8 ! Inverse of seconds per day - ! Cloud arrays for McICA + ! Cloud arrays for McICA integer, parameter :: nsubclw = ngptlw ! rrtmg_lw g-point (quadrature point) dimension integer :: permuteseed ! permute seed for sub-column generator @@ -139,10 +139,10 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & ! mji/rrtmg ! Calculate cloud optical properties here if using CAM method, or if using one of the - ! methods in RRTMG_LW, then pass in cloud physical properties and zero out cloud optical + ! methods in RRTMG_LW, then pass in cloud physical properties and zero out cloud optical ! properties here - - ! Zero optional cloud optical depth input array tauc_lw, + + ! Zero optional cloud optical depth input array tauc_lw, ! if inputting cloud physical properties into RRTMG_LW ! tauc_lw(:,:,:) = 0. ! Or, pass in CAM cloud longwave optical depth to RRTMG_LW @@ -155,7 +155,7 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & ! Call sub-column generator for McICA in radiation call t_startf('mcica_subcol_lw') - ! Set permute seed (must be offset between LW and SW by at least 140 to insure + ! Set permute seed (must be offset between LW and SW by at least 140 to insure ! effective randomization) permuteseed = 150 @@ -171,7 +171,7 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & call t_stopf('mcica_subcol_lw') - + call t_startf('rrtmg_lw') ! Convert incoming water amounts from specific humidity to vmr as needed; @@ -204,7 +204,7 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & ! extra layer above model top with vertical indexing from bottom to top. ! Heating units are in K/d on output from RRTMG and contain output for ! extra layer above model top with vertical indexing from bottom to top. - ! Heating units are converted to J/kg/s below for use in CAM. + ! Heating units are converted to J/kg/s below for use in CAM. flwds(:ncol) = dflx (:ncol,1) fldsc(:ncol) = dflxc(:ncol,1) @@ -227,17 +227,13 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & fsul(:ncol,pverp-rrtmg_levs+1:pverp)=uflxc(:ncol,rrtmg_levs:1:-1) fsdl(:ncol,pverp-rrtmg_levs+1:pverp)=dflxc(:ncol,rrtmg_levs:1:-1) -#ifndef OSLO_AERO - if (single_column.and.scm_crm_mode) then -#endif - call outfld('FUL ',ful,pcols,lchnk) - call outfld('FDL ',fdl,pcols,lchnk) - call outfld('FULC ',fsul,pcols,lchnk) - call outfld('FDLC ',fsdl,pcols,lchnk) -#ifndef OSLO_AERO - endif -#endif - + ! OSLO_AERO begin + call outfld('FUL ',ful,pcols,lchnk) + call outfld('FDL ',fdl,pcols,lchnk) + call outfld('FULC ',fsul,pcols,lchnk) + call outfld('FDLC ',fsdl,pcols,lchnk) + ! OSLO_AERO end + fnl(:ncol,:) = ful(:ncol,:) - fdl(:ncol,:) ! mji/ cam excluded this? fcnl(:ncol,:) = fsul(:ncol,:) - fsdl(:ncol,:) @@ -260,12 +256,12 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & lu(:ncol,pverp-rrtmg_levs+1:pverp,:) = reshape(lwuflxs(:,:ncol,rrtmg_levs:1:-1), & (/ncol,rrtmg_levs,nbndlw/), order=(/3,1,2/)) end if - + if (associated(ld)) then ld(:ncol,pverp-rrtmg_levs+1:pverp,:) = reshape(lwdflxs(:,:ncol,rrtmg_levs:1:-1), & (/ncol,rrtmg_levs,nbndlw/), order=(/3,1,2/)) end if - + call t_stopf('rrtmg_lw') end subroutine rad_rrtmg_lw @@ -273,9 +269,9 @@ end subroutine rad_rrtmg_lw !------------------------------------------------------------------------------- subroutine radlw_init() -!----------------------------------------------------------------------- -! -! Purpose: +!----------------------------------------------------------------------- +! +! Purpose: ! Initialize various constants for radiation scheme. ! !----------------------------------------------------------------------- diff --git a/src_cam/radsw.F90 b/src_cam/radsw.F90 index 0ec8be5..b76648a 100644 --- a/src_cam/radsw.F90 +++ b/src_cam/radsw.F90 @@ -1,7 +1,7 @@ module radsw -!----------------------------------------------------------------------- -! +!----------------------------------------------------------------------- +! ! Purpose: Solar radiation calculations. ! !----------------------------------------------------------------------- @@ -57,33 +57,33 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & !----------------------------------------------------------------------- -! -! Purpose: +! +! Purpose: ! Solar radiation code -! -! Method: +! +! Method: ! mji/rrtmg ! RRTMG, two-stream, with McICA -! +! ! Divides solar spectrum into 14 intervals from 0.2-12.2 micro-meters. ! solar flux fractions specified for each interval. allows for ! seasonally and diurnally varying solar input. Includes molecular, -! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud, +! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud, ! and surface absorption. Computes delta-eddington reflections and -! transmissions assuming homogeneously mixed layers. Adds the layers -! assuming scattering between layers to be isotropic, and distinguishes +! transmissions assuming homogeneously mixed layers. Adds the layers +! assuming scattering between layers to be isotropic, and distinguishes ! direct solar beam from scattered radiation. -! +! ! Longitude loops are broken into 1 or 2 sections, so that only daylight ! (i.e. coszrs > 0) computations are done. -! +! ! Note that an extra layer above the model top layer is added. -! +! ! mks units are used. -! +! ! Special diagnostic calculation of the clear sky surface and total column ! absorbed flux is also done for cloud forcing diagnostics. -! +! !----------------------------------------------------------------------- use cmparray_mod, only: CmpDayNite, ExpDayNite @@ -91,8 +91,8 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & use mcica_subcol_gen_sw, only: mcica_subcol_sw use physconst, only: cpair use rrtmg_state, only: rrtmg_state_t - - ! Minimum cloud amount (as a fraction of the grid-box area) to + + ! Minimum cloud amount (as a fraction of the grid-box area) to ! distinguish from clear sky real(r8), parameter :: cldmin = 1.0e-80_r8 @@ -126,12 +126,12 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8), intent(in) :: E_aldir(pcols) ! 0.7-5.0 micro-meter srfc alb: direct rad real(r8), intent(in) :: E_asdif(pcols) ! 0.2-0.7 micro-meter srfc alb: diffuse rad real(r8), intent(in) :: E_aldif(pcols) ! 0.7-5.0 micro-meter srfc alb: diffuse rad - real(r8), intent(in) :: sfac(nbndsw) ! factor to account for solar variability in each band + real(r8), intent(in) :: sfac(nbndsw) ! factor to account for solar variability in each band real(r8), optional, intent(in) :: E_cld_tau (nbndsw, pcols, pver) ! cloud optical depth - real(r8), optional, intent(in) :: E_cld_tau_w (nbndsw, pcols, pver) ! cloud optical - real(r8), optional, intent(in) :: E_cld_tau_w_g(nbndsw, pcols, pver) ! cloud optical - real(r8), optional, intent(in) :: E_cld_tau_w_f(nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w (nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w_g(nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w_f(nbndsw, pcols, pver) ! cloud optical logical , optional, intent(in) :: old_convert logical , optional, intent(in) :: idrf @@ -184,10 +184,10 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: h2ovmr(pcols,rrtmg_levs) ! h2o volume mixing ratio real(r8) :: o3vmr(pcols,rrtmg_levs) ! o3 volume mixing ratio - real(r8) :: co2vmr(pcols,rrtmg_levs) ! co2 volume mixing ratio - real(r8) :: ch4vmr(pcols,rrtmg_levs) ! ch4 volume mixing ratio - real(r8) :: o2vmr(pcols,rrtmg_levs) ! o2 volume mixing ratio - real(r8) :: n2ovmr(pcols,rrtmg_levs) ! n2o volume mixing ratio + real(r8) :: co2vmr(pcols,rrtmg_levs) ! co2 volume mixing ratio + real(r8) :: ch4vmr(pcols,rrtmg_levs) ! ch4 volume mixing ratio + real(r8) :: o2vmr(pcols,rrtmg_levs) ! o2 volume mixing ratio + real(r8) :: n2ovmr(pcols,rrtmg_levs) ! n2o volume mixing ratio real(r8) :: tsfc(pcols) ! surface temperature @@ -210,7 +210,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: asm_aer_sw(pcols, rrtmg_levs-1, nbndsw) ! aer asymmetry parameter real(r8) :: cld_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud fraction - real(r8) :: rei_stosw(pcols, rrtmg_levs-1) ! stochastic ice particle size + real(r8) :: rei_stosw(pcols, rrtmg_levs-1) ! stochastic ice particle size real(r8) :: rel_stosw(pcols, rrtmg_levs-1) ! stochastic liquid particle size real(r8) :: cicewp_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud ice water path real(r8) :: cliqwp_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud liquid wter path @@ -220,7 +220,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: fsfc_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud forward scattering fraction (optional) real(r8), parameter :: dps = 1._r8/86400._r8 ! Inverse of seconds per day - + real(r8) :: swuflx(pcols,rrtmg_levs+1) ! Total sky shortwave upward flux (W/m2) real(r8) :: swdflx(pcols,rrtmg_levs+1) ! Total sky shortwave downward flux (W/m2) real(r8) :: swhr(pcols,rrtmg_levs) ! Total sky shortwave radiative heating rate (K/d) @@ -305,16 +305,12 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & qrsc(1:ncol,1:pver) = 0.0_r8 fns(1:ncol,1:pverp) = 0.0_r8 fcns(1:ncol,1:pverp) = 0.0_r8 -#ifndef OSLO_AERO - if (single_column.and.scm_crm_mode) then -#endif - fus(1:ncol,1:pverp) = 0.0_r8 - fds(1:ncol,1:pverp) = 0.0_r8 - fusc(:ncol,:pverp) = 0.0_r8 - fdsc(:ncol,:pverp) = 0.0_r8 -#ifndef OSLO_AERO - endif -#endif + ! OSLO_AERO begin + fus(1:ncol,1:pverp) = 0.0_r8 + fds(1:ncol,1:pverp) = 0.0_r8 + fusc(:ncol,:pverp) = 0.0_r8 + fdsc(:ncol,:pverp) = 0.0_r8 + ! OSLO_AERO end if (associated(su)) su(1:ncol,:,:) = 0.0_r8 if (associated(sd)) sd(1:ncol,:,:) = 0.0_r8 @@ -394,10 +390,10 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & end do ! Calculate cloud optical properties here if using CAM method, or if using one of the - ! methods in RRTMG_SW, then pass in cloud physical properties and zero out cloud optical + ! methods in RRTMG_SW, then pass in cloud physical properties and zero out cloud optical ! properties here - ! Zero optional cloud optical property input arrays tauc_sw, ssac_sw, asmc_sw, + ! Zero optional cloud optical property input arrays tauc_sw, ssac_sw, asmc_sw, ! if inputting cloud physical properties to RRTMG_SW !tauc_sw(:,:,:) = 0.0_r8 !ssac_sw(:,:,:) = 1.0_r8 @@ -420,7 +416,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsfc_sw(ns,i,k) = 0._r8 asmc_sw(ns,i,k) = 0._r8 endif - + tauc_sw(ns,i,k)=E_cld_tau(ns,IdxDay(i),kk) if (tauc_sw(ns,i,k) > 0._r8) then ssac_sw(ns,i,k)=E_cld_tau_w(ns,IdxDay(i),kk)/tauc_sw(ns,i,k) @@ -446,7 +442,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsfc_sw(ns,i,k) = 0._r8 asmc_sw(ns,i,k) = 0._r8 endif - + tauc_sw(ns,i,k)=E_cld_tau(ns,IdxDay(i),kk) if (tauc_sw(ns,i,k) > 0._r8) then ssac_sw(ns,i,k)=max(E_cld_tau_w(ns,IdxDay(i),kk),1.e-80_r8)/max(tauc_sw(ns,i,k),1.e-80_r8) @@ -492,7 +488,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Call sub-column generator for McICA in radiation call t_startf('mcica_subcol_sw') - ! Set permute seed (must be offset between LW and SW by at least 140 to insure + ! Set permute seed (must be offset between LW and SW by at least 140 to insure ! effective randomization) permuteseed = 1 @@ -531,8 +527,8 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Flux units are in W/m2 on output from rrtmg_sw and contain output for ! extra layer above model top with vertical indexing from bottom to top. ! - ! Heating units are in J/kg/s on output from rrtmg_sw and contain output - ! for extra layer above model top with vertical indexing from bottom to top. + ! Heating units are in J/kg/s on output from rrtmg_sw and contain output + ! for extra layer above model top with vertical indexing from bottom to top. ! ! Reverse vertical indexing to go from top to bottom for CAM output. @@ -550,7 +546,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsnt(1:Nday) = swdflx(1:Nday,rrtmg_levs) - swuflx(1:Nday,rrtmg_levs) fsntc(1:Nday) = swdflxc(1:Nday,rrtmg_levs) - swuflxc(1:Nday,rrtmg_levs) - ! Set the downwelling flux at the surface + ! Set the downwelling flux at the surface fsds(1:Nday) = swdflx(1:Nday,1) fsdsc(1:Nday) = swdflxc(1:Nday,1) @@ -627,21 +623,17 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & end if ! these outfld calls don't work for spmd only outfield in scm mode (nonspmd) -#ifndef OSLO_AERO - if (single_column .and. scm_crm_mode) then -#endif - ! Following outputs added for CRM - call ExpDayNite(fus,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) - call ExpDayNite(fds,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) - call ExpDayNite(fusc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) - call ExpDayNite(fdsc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) - call outfld('FUS ', fus(:ncol,:), ncol, lchnk) - call outfld('FDS ', fds(:ncol,:), ncol, lchnk) - call outfld('FUSC ', fusc(:ncol,:), ncol, lchnk) - call outfld('FDSC ', fdsc(:ncol,:), ncol, lchnk) -#ifndef OSLO_AERO - end if -#endif + ! OSLO_AERO begin + ! Following outputs added for CRM + call ExpDayNite(fus,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) + call ExpDayNite(fds,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) + call ExpDayNite(fusc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) + call ExpDayNite(fdsc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) + call outfld('FUS ', fus(:ncol,:), ncol, lchnk) + call outfld('FDS ', fds(:ncol,:), ncol, lchnk) + call outfld('FUSC ', fusc(:ncol,:), ncol, lchnk) + call outfld('FDSC ', fdsc(:ncol,:), ncol, lchnk) + ! OSLO_AERO end if (present(idrf)) then if (idrf) then call outfld('FUSCDRF ', fusc(:ncol,:), ncol, lchnk) @@ -654,9 +646,9 @@ end subroutine rad_rrtmg_sw !------------------------------------------------------------------------------- subroutine radsw_init() -!----------------------------------------------------------------------- -! -! Purpose: +!----------------------------------------------------------------------- +! +! Purpose: ! Initialize various constants for radiation scheme. ! !----------------------------------------------------------------------- @@ -669,7 +661,7 @@ subroutine radsw_init() ! Initialize rrtmg_sw call rrtmg_sw_ini - + end subroutine radsw_init diff --git a/src_cam/vertical_diffusion.F90 b/src_cam/vertical_diffusion.F90 index 6ad60e4..68478a5 100644 --- a/src_cam/vertical_diffusion.F90 +++ b/src_cam/vertical_diffusion.F90 @@ -72,9 +72,9 @@ module vertical_diffusion use ref_pres, only : do_molec_diff, nbot_molec use phys_control, only : phys_getopts use time_manager, only : is_first_step -#ifdef OSLO_AERO +! OSLO_AERO begin use oslo_aero_share, only: getNumberOfAerosolTracers, fillAerosolTracerList -#endif +! OSLO_AERO end implicit none private @@ -324,36 +324,12 @@ subroutine vertical_diffusion_init(pbuf2d) ! aerosols. ! N.B. - This implementation assumes that the prognostic modal aerosols are ! impacting the climate calculation (i.e., can get info from list 0). -#ifdef OSLO_AERO + ! OSLO_AERO begin prog_modal_aero = .true. pmam_ncnst = getNumberOfAerosolTracers() allocate(pmam_cnst_idx(pmam_ncnst)) call fillAerosolTracerList(pmam_cnst_idx) -#else - call phys_getopts(prog_modal_aero_out=prog_modal_aero) - if (prog_modal_aero) then - ! First need total number of mam constituents - call rad_cnst_get_info(0, nmodes=nmodes) - do m = 1, nmodes - call rad_cnst_get_info(0, m, nspec=nspec) - pmam_ncnst = pmam_ncnst + 1 + nspec - end do - - allocate(pmam_cnst_idx(pmam_ncnst)) - - ! Get the constituent indicies - im = 1 - do m = 1, nmodes - call rad_cnst_get_mode_num_idx(m, pmam_cnst_idx(im)) - im = im + 1 - call rad_cnst_get_info(0, m, nspec=nspec) - do l = 1, nspec - call rad_cnst_get_mam_mmr_idx(m, l, pmam_cnst_idx(im)) - im = im + 1 - end do - end do - end if -#endif + ! OSLO_AERO end ! Initialize upper boundary condition module @@ -577,10 +553,10 @@ subroutine vertical_diffusion_init(pbuf2d) if( history_budget ) then call add_default( vdiffnam(ixcldliq), history_budget_histfile_num, ' ' ) call add_default( vdiffnam(ixcldice), history_budget_histfile_num, ' ' ) -#ifdef OSLO_AERO + ! OSLO_AERO begin call add_default( vdiffnam(ixnumliq), history_budget_histfile_num, ' ' ) call add_default( vdiffnam(ixnumice), history_budget_histfile_num, ' ' ) -#endif + ! OSLO_AERO end if( history_budget_histfile_num > 1 ) then call add_default( vdiffnam(1), history_budget_histfile_num, ' ' ) call add_default( 'DTV' , history_budget_histfile_num, ' ' ) @@ -1181,22 +1157,9 @@ subroutine vertical_diffusion_tend( & end if -#ifdef OSLO_AERO ! Do nothing if OSLO_AERO ! Oslo aero adds emissions together with dry deposition - so do not add the explicit ! surface fluxes to the lowest layer -#else - if (prog_modal_aero) then - ! Modal aerosol species not diffused, so just add the explicit surface fluxes to the - ! lowest layer. **NOTE** This code assumes wet mmr. - - tmp1(:ncol) = ztodt * gravit * state%rpdel(:ncol,pver) - do m = 1, pmam_ncnst - l = pmam_cnst_idx(m) - q_tmp(:ncol,pver,l) = q_tmp(:ncol,pver,l) + tmp1(:ncol) * cflux(:ncol,l) - enddo - end if -#endif ! -------------------------------------------------------- ! ! Diagnostics and output writing after applying PBL scheme !