From af7f2ccf8edb8dcc3823b7e631c50b63eadb23c0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Alf=20Kirkev=C3=A5g?= Date: Mon, 15 Jun 2020 11:24:47 +0200 Subject: [PATCH] #12: Added CMIP6 SpAer code changes (as well as some not used code for IRF-RFMIP) --- src/NorESM/cam_diagnostics.F90 | 178 ++-- src/NorESM/cloud_rad_props.F90 | 848 ++++++++++++++++++ src/NorESM/micro_mg2_0.F90 | 43 +- src/NorESM/micro_mg_cam.F90 | 179 +++- src/NorESM/physpkg.F90 | 39 +- src/physics/cam_oslo/mo_simple_plumes_v1.F90 | 399 ++++++++ src/physics/cam_oslo/pmxsub.F90 | 19 +- src/physics/cam_oslo/radconstants.F90 | 281 ++++++ src/physics/cam_oslo/radiation.F90 | 513 +++++------ src/physics/cam_oslo/radlw.F90 | 24 +- src/physics/cam_oslo/radsw.F90 | 47 +- .../cam_oslo/simple_plumes_interface.F90 | 242 +++++ 12 files changed, 2394 insertions(+), 418 deletions(-) create mode 100644 src/NorESM/cloud_rad_props.F90 create mode 100644 src/physics/cam_oslo/mo_simple_plumes_v1.F90 create mode 100644 src/physics/cam_oslo/radconstants.F90 create mode 100644 src/physics/cam_oslo/simple_plumes_interface.F90 diff --git a/src/NorESM/cam_diagnostics.F90 b/src/NorESM/cam_diagnostics.F90 index e8dff42def..07e54b61f2 100644 --- a/src/NorESM/cam_diagnostics.F90 +++ b/src/NorESM/cam_diagnostics.F90 @@ -186,9 +186,9 @@ subroutine diag_init_dry(pbuf2d) !+ #ifdef AEROCOM use commondefinitions, only: nbmodes -!#ifdef RFMIPIRF -! use radconstants, only: nswbands, nlwbands -!#endif +#ifdef RFMIPIRF + use radconstants, only: nswbands, nlwbands +#endif #endif !- @@ -206,10 +206,10 @@ subroutine diag_init_dry(pbuf2d) character(len=10) :: modeString character(len=20) :: varname integer :: i, irh -!#ifdef RFMIPIRF -! character(len=2) :: c2 -! integer :: ib -!#endif +#ifdef RFMIPIRF + character(len=2) :: c2 + integer :: ib +#endif #endif !- @@ -221,7 +221,6 @@ subroutine diag_init_dry(pbuf2d) call addfld ('PS', horiz_only, 'A', 'Pa', 'Surface pressure') call addfld ('T', (/ 'lev' /), 'A', 'K', 'Temperature') call addfld ('U', (/ 'lev' /), 'A', 'm/s', 'Zonal wind') - call addfld ('UA010', horiz_only, 'A', 'm/s', 'Zonal wind U at 10 mbar pressure surface') call addfld ('V', (/ 'lev' /), 'A', 'm/s', 'Meridional wind') call register_vector_field('U','V') @@ -259,12 +258,10 @@ subroutine diag_init_dry(pbuf2d) call addfld ('Z200', horiz_only, 'A', 'm', 'Geopotential Z at 200 mbar pressure surface') call addfld ('Z100', horiz_only, 'A', 'm', 'Geopotential Z at 100 mbar pressure surface') call addfld ('Z050', horiz_only, 'A', 'm', 'Geopotential Z at 50 mbar pressure surface') - call addfld ('Z010', horiz_only, 'A', 'm', 'Geopotential Z at 10 mbar pressure surface') call addfld ('ZZ', (/ 'lev' /), 'A', 'm2', 'Eddy height variance' ) call addfld ('VZ', (/ 'lev' /), 'A', 'm2/s', 'Meridional transport of geopotential height') call addfld ('VT', (/ 'lev' /), 'A', 'K m/s ', 'Meridional heat transport') - call addfld ('VT100', horiz_only, 'A', 'K m/s ', 'Meridional heat transport at 100 mbar pressure level') call addfld ('VU', (/ 'lev' /), 'A', 'm2/s2', 'Meridional flux of zonal momentum' ) call addfld ('VV', (/ 'lev' /), 'A', 'm2/s2', 'Meridional velocity squared' ) call addfld ('OMEGAV', (/ 'lev' /), 'A', 'm Pa/s2 ', 'Vertical flux of meridional momentum' ) @@ -327,7 +324,7 @@ subroutine diag_init_dry(pbuf2d) call addfld ('ATMEINT', horiz_only, 'A', 'J/m2','Vertically integrated total atmospheric energy ') -!akc6+ CNVCLD is zero... +!akc6+ CNVCLD is actually zero... ! call addfld ('CNVCLD', horiz_only, 'A', 'fraction', 'Vertically integrated convective cloud cover') !akc6- @@ -352,13 +349,13 @@ subroutine diag_init_dry(pbuf2d) !akc6+ call addfld ('BVISVOLC ',(/'lev'/), 'A','1/km ','CMIP6 volcanic aerosol extinction at 0.442-0.625um') !akc6- -!#ifdef SPAERO -! call addfld ('AODVISSP',horiz_only, 'A','unitless' ,'Simple plumes aerosol optical depth at 0.35-0.64um') -! call addfld ('ABSVISSP',horiz_only, 'A','unitless' ,'Simple plumes aerosol absorptive optical depth at 0.35-0.64um') -! call addfld ('XCDNC_SP',horiz_only, 'A','unitless' ,'CDNC modification factor for simple plume aerosols') -! call addfld ('AODV3DSP',(/'lev'/), 'A','unitless','Simple plumes 3D aerosol optical depth at 0.35-0.64um') -! call addfld ('ABSV3DSP',(/'lev'/), 'A','unitless','Simple plumes 3D absorption AOD at 0.35-0.64um') -!#endif +#ifdef SPAERO + call addfld ('AODVISSP',horiz_only, 'A','unitless' ,'Simple plumes aerosol optical depth at 0.35-0.64um') + call addfld ('ABSVISSP',horiz_only, 'A','unitless' ,'Simple plumes aerosol absorptive optical depth at 0.35-0.64um') + call addfld ('XCDNC_SP',horiz_only, 'A','unitless' ,'CDNC modification factor for simple plume aerosols') + call addfld ('AODV3DSP',(/'lev'/), 'A','unitless','Simple plumes 3D aerosol optical depth at 0.35-0.64um') + call addfld ('ABSV3DSP',(/'lev'/), 'A','unitless','Simple plumes 3D absorption AOD at 0.35-0.64um') +#endif #ifdef COLTST4INTCONS ! optical depth for each mode/mixture: call addfld ('TAUKC0 ',horiz_only, 'A','unitless','Aerosol optical depth at 0.442-0.625um for kcomp 0') @@ -426,7 +423,6 @@ subroutine diag_init_dry(pbuf2d) call addfld ('PM2P5 ',(/'lev'/), 'A','ug/m3 ','3D aerosol PM2.5') call addfld ('MMRPM2P5',(/'lev'/), 'A','kg/kg ','3D aerosol PM2.5 mass mixing ratio') call addfld ('MMRPM1 ',(/'lev'/), 'A','kg/kg ','3D aerosol PM1.0 mass mixing ratio') - call addfld ('MMRPM2P5_SRF',horiz_only, 'A','kg/kg ','Aerosol PM2.5 mass mixing ratio in bottom layer') !akc6- call addfld ('GRIDAREA',horiz_only, 'A','m2 ','Grid area for 1.9x2.5 horizontal resolution') call addfld ('DAERH2O ',horiz_only, 'A', 'mg/m2 ','Aerosol water load') @@ -588,22 +584,22 @@ subroutine diag_init_dry(pbuf2d) if(i.ne.3) call addfld(varName, horiz_only, 'A', 'unitless', 'relative exessive added mass column for mode'//modeString) enddo -!#ifdef RFMIPIRF -! do ib=1,nswbands -! write(c2,'(I2)') ib -! call addfld('AERTAUBND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol extinction optical depth for wavelength band '//trim(adjustl(c2))) -! call addfld('AERSSABND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol single scattering albedo for wavelength band '//c2) -! call addfld('AERASYBND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol asymmetry parameter for wavelength band '//c2) -! -! call addfld('SDBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'shortwave spectral flux down for wavelength band '//c2) -! call addfld('SUBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'shortwave spectral flux up for wavelength band '//c2) -! enddo -! do ib=1,nlwbands -! write(c2,'(I2)') ib -! call addfld('LDBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'longwave spectral flux down for wavelength band '//c2) -! call addfld('LUBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'longwave spectral flux up for wavelength band '//c2) -! enddo -!#endif +#ifdef RFMIPIRF + do ib=1,nswbands + write(c2,'(I2)') ib + call addfld('AERTAUBND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol extinction optical depth for wavelength band '//trim(adjustl(c2))) + call addfld('AERSSABND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol single scattering albedo for wavelength band '//c2) + call addfld('AERASYBND'//trim(adjustl(c2)), (/'lev'/),'A', 'unitless', 'aerosol asymmetry parameter for wavelength band '//c2) + + call addfld('SDBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'shortwave spectral flux down for wavelength band '//c2) + call addfld('SUBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'shortwave spectral flux up for wavelength band '//c2) + enddo + do ib=1,nlwbands + write(c2,'(I2)') ib + call addfld('LDBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'longwave spectral flux down for wavelength band '//c2) + call addfld('LUBND'//trim(adjustl(c2)), (/'ilev'/),'A', 'W/m^2', 'longwave spectral flux up for wavelength band '//c2) + enddo +#endif #ifdef AEROCOM_INSITU ! Note that this code has not yet been updated to CESM2 standard @@ -832,13 +828,13 @@ subroutine diag_init_dry(pbuf2d) !akc6+ call add_default ('BVISVOLC', 1, ' ') !akc6- -!#ifdef SPAERO -! call add_default ('AODVISSP', 1, ' ') -! call add_default ('ABSVISSP', 1, ' ') -! call add_default ('XCDNC_SP', 1, ' ') -! call add_default ('AODV3DSP', 1, ' ') -! call add_default ('ABSV3DSP', 1, ' ') -!#endif +#ifdef SPAERO + call add_default ('AODVISSP', 1, ' ') + call add_default ('ABSVISSP', 1, ' ') + call add_default ('XCDNC_SP', 1, ' ') + call add_default ('AODV3DSP', 1, ' ') + call add_default ('ABSV3DSP', 1, ' ') +#endif #ifdef AEROFFL call add_default ('FSNT_DRF', 1, ' ') call add_default ('FSNTCDRF', 1, ' ') @@ -1020,23 +1016,22 @@ subroutine diag_init_dry(pbuf2d) enddo !++ -!#ifdef RFMIPIRF -! do i=1,nbands -! do ib=1,nswbands -! write(c2,'(I2)') ib -! call add_default('AERTAUBND'//trim(adjustl(c2)), 1, ' ') -! call add_default('AERSSABND'//trim(adjustl(c2)), 1, ' ') -! call add_default('AERASYBND'//trim(adjustl(c2)), 1, ' ') -! -! call add_default('SDBND'//trim(adjustl(c2)), 1, ' ') -! call add_default('SUBND'//trim(adjustl(c2)), 1, ' ') -! enddo -! do ib=1,nlwbands -! write(c2,'(I2)') ib -! call add_default('LDBND'//trim(adjustl(c2)), 1, ' ') -! call add_default('LUBND'//trim(adjustl(c2)), 1, ' ') -! enddo -!#endif +#ifdef RFMIPIRF + do ib=1,nswbands + write(c2,'(I2)') ib + call add_default('AERTAUBND'//trim(adjustl(c2)), 1, ' ') + call add_default('AERSSABND'//trim(adjustl(c2)), 1, ' ') + call add_default('AERASYBND'//trim(adjustl(c2)), 1, ' ') + + call add_default('SDBND'//trim(adjustl(c2)), 1, ' ') + call add_default('SUBND'//trim(adjustl(c2)), 1, ' ') + enddo + do ib=1,nlwbands + write(c2,'(I2)') ib + call add_default('LDBND'//trim(adjustl(c2)), 1, ' ') + call add_default('LUBND'//trim(adjustl(c2)), 1, ' ') + enddo +#endif #ifdef AEROCOM_INSITU @@ -1071,32 +1066,32 @@ subroutine diag_init_dry(pbuf2d) #endif ! aerocom #endif ! dirind -!#ifdef SPAERO -! call addfld ('FSNT_SP ', horiz_only, 'A','W/m^2','Total column absorbed solar flux (without SP aerosols)') -! call addfld ('FSNTC_SP', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (without SP aerosols)') -! call addfld ('FSNS_SP ', horiz_only, 'A','W/m^2','Surface absorbed solar flux (without SP aerosols)') -! call addfld ('FSNSC_SP', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (without SP aerosols)') -! call addfld ('FSNT_SP2', horiz_only, 'A','W/m^2','Total column absorbed solar flux (SP aerosols for DRF only)') -! call addfld ('FSNTCSP2', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (SP aerosols for DRF only)') -! call addfld ('FSNS_SP2', horiz_only, 'A','W/m^2','Surface absorbed solar flux (SP aerosols for DRF only)') -! call addfld ('FSNSCSP2', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (SP aerosols for DRF only)') -! call addfld ('FSNT_SP3', horiz_only, 'A','W/m^2','Total column absorbed solar flux (SP aerosols)') -! call addfld ('FSNTCSP3', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (SP aerosols)') -! call addfld ('FSNS_SP3', horiz_only, 'A','W/m^2','Surface absorbed solar flux (SP aerosols)') -! call addfld ('FSNSCSP3', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (SP aerosols)') -! call add_default ('FSNT_SP' , 1, ' ') -! call add_default ('FSNTC_SP', 1, ' ') -! call add_default ('FSNS_SP' , 1, ' ') -! call add_default ('FSNSC_SP', 1, ' ') -! call add_default ('FSNT_SP2', 1, ' ') -! call add_default ('FSNTCSP2', 1, ' ') -! call add_default ('FSNS_SP2', 1, ' ') -! call add_default ('FSNSCSP2', 1, ' ') -! call add_default ('FSNT_SP3', 1, ' ') -! call add_default ('FSNTCSP3', 1, ' ') -! call add_default ('FSNS_SP3', 1, ' ') -! call add_default ('FSNSCSP3', 1, ' ') -!#endif +#ifdef SPAERO + call addfld ('FSNT_SP ', horiz_only, 'A','W/m^2','Total column absorbed solar flux (without SP aerosols)') + call addfld ('FSNTC_SP', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (without SP aerosols)') + call addfld ('FSNS_SP ', horiz_only, 'A','W/m^2','Surface absorbed solar flux (without SP aerosols)') + call addfld ('FSNSC_SP', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (without SP aerosols)') + call addfld ('FSNT_SP2', horiz_only, 'A','W/m^2','Total column absorbed solar flux (SP aerosols for DRF only)') + call addfld ('FSNTCSP2', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (SP aerosols for DRF only)') + call addfld ('FSNS_SP2', horiz_only, 'A','W/m^2','Surface absorbed solar flux (SP aerosols for DRF only)') + call addfld ('FSNSCSP2', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (SP aerosols for DRF only)') + call addfld ('FSNT_SP3', horiz_only, 'A','W/m^2','Total column absorbed solar flux (SP aerosols)') + call addfld ('FSNTCSP3', horiz_only, 'A','W/m^2','Clear sky total column absorbed solar flux (SP aerosols)') + call addfld ('FSNS_SP3', horiz_only, 'A','W/m^2','Surface absorbed solar flux (SP aerosols)') + call addfld ('FSNSCSP3', horiz_only, 'A','W/m^2','Clear sky surface absorbed solar flux (SP aerosols)') + call add_default ('FSNT_SP' , 1, ' ') + call add_default ('FSNTC_SP', 1, ' ') + call add_default ('FSNS_SP' , 1, ' ') + call add_default ('FSNSC_SP', 1, ' ') + call add_default ('FSNT_SP2', 1, ' ') + call add_default ('FSNTCSP2', 1, ' ') + call add_default ('FSNS_SP2', 1, ' ') + call add_default ('FSNSCSP2', 1, ' ') + call add_default ('FSNT_SP3', 1, ' ') + call add_default ('FSNTCSP3', 1, ' ') + call add_default ('FSNS_SP3', 1, ' ') + call add_default ('FSNSCSP3', 1, ' ') +#endif end subroutine diag_init_dry @@ -1714,14 +1709,6 @@ subroutine diag_phys_writeout_dry(state, pbuf, p_surf_t) call vertinterp(ncol, pcols, pver, state%pmid, 5000._r8, z3, p_surf, ln_interp=.true.) call outfld('Z050 ', p_surf, pcols, lchnk) end if - if (hist_fld_active('Z010')) then - call vertinterp(ncol, pcols, pver, state%pmid, 1000._r8, z3, p_surf, ln_interp=.true.) - call outfld('Z010 ', p_surf, pcols, lchnk) - end if - if (hist_fld_active('UA010')) then - call vertinterp(ncol, pcols, pver, state%pmid, 1000._r8, state%u, p_surf, ln_interp=.true.) - call outfld('UA010 ', p_surf, pcols, lchnk) - end if ! ! Quadratic height fiels Z3*Z3 ! @@ -1736,11 +1723,6 @@ subroutine diag_phys_writeout_dry(state, pbuf, p_surf_t) ftem(:ncol,:) = state%v(:ncol,:)*state%t(:ncol,:) call outfld ('VT ',ftem ,pcols ,lchnk ) - if (hist_fld_active('VT100')) then - call vertinterp(ncol, pcols, pver, state%pmid, 10000._r8, ftem, p_surf, ln_interp=.true.) - call outfld('VT100 ', p_surf, pcols, lchnk) - end if - ftem(:ncol,:) = state%v(:ncol,:)**2 call outfld ('VV ',ftem ,pcols ,lchnk ) diff --git a/src/NorESM/cloud_rad_props.F90 b/src/NorESM/cloud_rad_props.F90 new file mode 100644 index 0000000000..d821f81371 --- /dev/null +++ b/src/NorESM/cloud_rad_props.F90 @@ -0,0 +1,848 @@ +module cloud_rad_props + +!------------------------------------------------------------------------------------------------ +!------------------------------------------------------------------------------------------------ +! May 2019 A. Kirkevåg Added MACv2 SP-aerosol code wrt. Twomey effect (SPAERO) +!------------------------------------------------------------------------------------------------ +!ak+ +#include +!ak- + +use shr_kind_mod, only: r8 => shr_kind_r8 +use ppgrid, only: pcols, pver, pverp +use physics_types, only: physics_state +use physics_buffer, only: physics_buffer_desc, pbuf_get_index, pbuf_get_field, pbuf_old_tim_idx +use radconstants, only: nswbands, nlwbands, idx_sw_diag, ot_length, idx_lw_diag +use cam_abortutils, only: endrun +use rad_constituents, only: iceopticsfile, liqopticsfile +use oldcloud, only: oldcloud_lw, old_liq_get_rad_props_lw, old_ice_get_rad_props_lw, oldcloud_init + +use ebert_curry, only: scalefactor +use cam_logfile, only: iulog + +use interpolate_data, only: interp_type, lininterp_init, lininterp, & + extrap_method_bndry, lininterp_finish + +implicit none +private +save + +public :: & + cloud_rad_props_init, & + cloud_rad_props_get_sw, & ! return SW optical props of total bulk aerosols + cloud_rad_props_get_lw, & ! return LW optical props of total bulk aerosols + get_ice_optics_sw, & ! return Mitchell SW ice radiative properties + ice_cloud_get_rad_props_lw, & ! Mitchell LW ice rad props + get_liquid_optics_sw, & ! return Conley SW rad props + liquid_cloud_get_rad_props_lw, & ! return Conley LW rad props + snow_cloud_get_rad_props_lw, & + get_snow_optics_sw + +integer :: nmu, nlambda +real(r8), allocatable :: g_mu(:) ! mu samples on grid +real(r8), allocatable :: g_lambda(:,:) ! lambda scale samples on grid +real(r8), allocatable :: ext_sw_liq(:,:,:) +real(r8), allocatable :: ssa_sw_liq(:,:,:) +real(r8), allocatable :: asm_sw_liq(:,:,:) +real(r8), allocatable :: abs_lw_liq(:,:,:) + +integer :: n_g_d +real(r8), allocatable :: g_d_eff(:) ! radiative effective diameter samples on grid +real(r8), allocatable :: ext_sw_ice(:,:) +real(r8), allocatable :: ssa_sw_ice(:,:) +real(r8), allocatable :: asm_sw_ice(:,:) +real(r8), allocatable :: abs_lw_ice(:,:) + +! +! indexes into pbuf for optical parameters of MG clouds +! + integer :: i_dei, i_mu, i_lambda, i_iciwp, i_iclwp, i_des, i_icswp +#ifdef SPAERO + integer :: i_sp_mu, i_sp_lambda +#endif + +! indexes into constituents for old optics + integer :: & + ixcldice, & ! cloud ice water index + ixcldliq ! cloud liquid water index + + +!============================================================================== +contains +!============================================================================== + +subroutine cloud_rad_props_init() + + use netcdf + use spmd_utils, only: masterproc + use ioFileMod, only: getfil + use error_messages, only: handle_ncerr +#if ( defined SPMD ) + use mpishorthand +#endif + use constituents, only: cnst_get_ind + use slingo, only: slingo_rad_props_init + use ebert_curry, only: ec_rad_props_init, scalefactor + + character(len=256) :: liquidfile + character(len=256) :: icefile + character(len=256) :: locfn + + integer :: ncid, dimid, f_nlwbands, f_nswbands, ierr + integer :: vdimids(NF90_MAX_VAR_DIMS), ndims, templen + ! liquid clouds + integer :: mudimid, lambdadimid + integer :: mu_id, lambda_id, ext_sw_liq_id, ssa_sw_liq_id, asm_sw_liq_id, abs_lw_liq_id + + ! ice clouds + integer :: d_dimid ! diameters + integer :: d_id, ext_sw_ice_id, ssa_sw_ice_id, asm_sw_ice_id, abs_lw_ice_id + + integer :: err + + liquidfile = liqopticsfile + icefile = iceopticsfile + + call slingo_rad_props_init + call ec_rad_props_init + call oldcloud_init + + i_dei = pbuf_get_index('DEI',errcode=err) + i_mu = pbuf_get_index('MU',errcode=err) + i_lambda = pbuf_get_index('LAMBDAC',errcode=err) + i_iciwp = pbuf_get_index('ICIWP',errcode=err) + i_iclwp = pbuf_get_index('ICLWP',errcode=err) + i_des = pbuf_get_index('DES',errcode=err) + i_icswp = pbuf_get_index('ICSWP',errcode=err) +#ifdef SPAERO + i_sp_mu = pbuf_get_index('SP_MU',errcode=err) + i_sp_lambda = pbuf_get_index('SP_LAMBDAC',errcode=err) +!tst4 i_sp_mu = i_mu +!tst4 i_sp_lambda = i_lambda +#endif + + ! old optics + call cnst_get_ind('CLDICE', ixcldice) + call cnst_get_ind('CLDLIQ', ixcldliq) + + ! read liquid cloud optics + if(masterproc) then + call getfil( trim(liquidfile), locfn, 0) + call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'liquid optics file missing') + write(iulog,*)' reading liquid cloud optics from file ',locfn + + call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim') + call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands') + if (f_nlwbands /= nlwbands) call endrun('number of lw bands does not match') + + call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim') + call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands') + if (f_nswbands /= nswbands) call endrun('number of sw bands does not match') + + call handle_ncerr(nf90_inq_dimid( ncid, 'mu', mudimid), 'getting mu dim') + call handle_ncerr(nf90_inquire_dimension( ncid, mudimid, len=nmu), 'getting n mu samples') + + call handle_ncerr(nf90_inq_dimid( ncid, 'lambda_scale', lambdadimid), 'getting lambda dim') + call handle_ncerr(nf90_inquire_dimension( ncid, lambdadimid, len=nlambda), 'getting n lambda samples') + endif ! if (masterproc) + +#if ( defined SPMD ) + call mpibcast(nmu, 1, mpiint, 0, mpicom, ierr) + call mpibcast(nlambda, 1, mpiint, 0, mpicom, ierr) +#endif + + allocate(g_mu(nmu)) + allocate(g_lambda(nmu,nlambda)) + allocate(ext_sw_liq(nmu,nlambda,nswbands) ) + allocate(ssa_sw_liq(nmu,nlambda,nswbands)) + allocate(asm_sw_liq(nmu,nlambda,nswbands)) + allocate(abs_lw_liq(nmu,nlambda,nlwbands)) + + if(masterproc) then + call handle_ncerr( nf90_inq_varid(ncid, 'mu', mu_id),& + 'cloud optics mu get') + call handle_ncerr( nf90_get_var(ncid, mu_id, g_mu),& + 'read cloud optics mu values') + + call handle_ncerr( nf90_inq_varid(ncid, 'lambda', lambda_id),& + 'cloud optics lambda get') + call handle_ncerr( nf90_get_var(ncid, lambda_id, g_lambda),& + 'read cloud optics lambda values') + + call handle_ncerr( nf90_inq_varid(ncid, 'k_ext_sw', ext_sw_liq_id),& + 'cloud optics ext_sw_liq get') + call handle_ncerr( nf90_get_var(ncid, ext_sw_liq_id, ext_sw_liq),& + 'read cloud optics ext_sw_liq values') + + call handle_ncerr( nf90_inq_varid(ncid, 'ssa_sw', ssa_sw_liq_id),& + 'cloud optics ssa_sw_liq get') + call handle_ncerr( nf90_get_var(ncid, ssa_sw_liq_id, ssa_sw_liq),& + 'read cloud optics ssa_sw_liq values') + + call handle_ncerr( nf90_inq_varid(ncid, 'asm_sw', asm_sw_liq_id),& + 'cloud optics asm_sw_liq get') + call handle_ncerr( nf90_get_var(ncid, asm_sw_liq_id, asm_sw_liq),& + 'read cloud optics asm_sw_liq values') + + call handle_ncerr( nf90_inq_varid(ncid, 'k_abs_lw', abs_lw_liq_id),& + 'cloud optics abs_lw_liq get') + call handle_ncerr( nf90_get_var(ncid, abs_lw_liq_id, abs_lw_liq),& + 'read cloud optics abs_lw_liq values') + + call handle_ncerr( nf90_close(ncid), 'liquid optics file missing') + endif ! if masterproc + +#if ( defined SPMD ) + call mpibcast(g_mu, nmu, mpir8, 0, mpicom, ierr) + call mpibcast(g_lambda, nmu*nlambda, mpir8, 0, mpicom, ierr) + call mpibcast(ext_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(ssa_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(asm_sw_liq, nmu*nlambda*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(abs_lw_liq, nmu*nlambda*nlwbands, mpir8, 0, mpicom, ierr) +#endif + ! I forgot to convert kext from m^2/Volume to m^2/Kg + ext_sw_liq = ext_sw_liq / 0.9970449e3_r8 + abs_lw_liq = abs_lw_liq / 0.9970449e3_r8 + + ! read ice cloud optics + if(masterproc) then + call getfil( trim(icefile), locfn, 0) + call handle_ncerr( nf90_open(locfn, NF90_NOWRITE, ncid), 'ice optics file missing') + write(iulog,*)' reading ice cloud optics from file ',locfn + + call handle_ncerr(nf90_inq_dimid( ncid, 'lw_band', dimid), 'getting lw_band dim') + call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nlwbands), 'getting n lw bands') + if (f_nlwbands /= nlwbands) call endrun('number of lw bands does not match') + + call handle_ncerr(nf90_inq_dimid( ncid, 'sw_band', dimid), 'getting sw_band_dim') + call handle_ncerr(nf90_inquire_dimension( ncid, dimid, len=f_nswbands), 'getting n sw bands') + if (f_nswbands /= nswbands) call endrun('number of sw bands does not match') + + call handle_ncerr(nf90_inq_dimid( ncid, 'd_eff', d_dimid), 'getting deff dim') + call handle_ncerr(nf90_inquire_dimension( ncid, d_dimid, len=n_g_d), 'getting n deff samples') + + endif ! if (masterproc) + +#if ( defined SPMD ) + call mpibcast(n_g_d, 1, mpiint, 0, mpicom, ierr) +! call mpibcast(nswbands, 1, mpiint, 0, mpicom, ierr) +! call mpibcast(nlwbands, 1, mpiint, 0, mpicom, ierr) +#endif + + allocate(g_d_eff(n_g_d)) + allocate(ext_sw_ice(n_g_d,nswbands)) + allocate(ssa_sw_ice(n_g_d,nswbands)) + allocate(asm_sw_ice(n_g_d,nswbands)) + allocate(abs_lw_ice(n_g_d,nlwbands)) + + if(masterproc) then + call handle_ncerr( nf90_inq_varid(ncid, 'd_eff', d_id),& + 'cloud optics deff get') + call handle_ncerr( nf90_get_var(ncid, d_id, g_d_eff),& + 'read cloud optics deff values') + + call handle_ncerr( nf90_inq_varid(ncid, 'sw_ext', ext_sw_ice_id),& + 'cloud optics ext_sw_ice get') + call handle_ncerr(nf90_inquire_variable ( ncid, ext_sw_ice_id, ndims=ndims, dimids=vdimids),& + 'checking dimensions of ext_sw_ice') + call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(1), len=templen),& + 'getting first dimension sw_ext') + !write(iulog,*) 'expected length',n_g_d,'actual len',templen + call handle_ncerr(nf90_inquire_dimension( ncid, vdimids(2), len=templen),& + 'getting first dimension sw_ext') + !write(iulog,*) 'expected length',nswbands,'actual len',templen + call handle_ncerr( nf90_get_var(ncid, ext_sw_ice_id, ext_sw_ice),& + 'read cloud optics ext_sw_ice values') + + call handle_ncerr( nf90_inq_varid(ncid, 'sw_ssa', ssa_sw_ice_id),& + 'cloud optics ssa_sw_ice get') + call handle_ncerr( nf90_get_var(ncid, ssa_sw_ice_id, ssa_sw_ice),& + 'read cloud optics ssa_sw_ice values') + + call handle_ncerr( nf90_inq_varid(ncid, 'sw_asm', asm_sw_ice_id),& + 'cloud optics asm_sw_ice get') + call handle_ncerr( nf90_get_var(ncid, asm_sw_ice_id, asm_sw_ice),& + 'read cloud optics asm_sw_ice values') + + call handle_ncerr( nf90_inq_varid(ncid, 'lw_abs', abs_lw_ice_id),& + 'cloud optics abs_lw_ice get') + call handle_ncerr( nf90_get_var(ncid, abs_lw_ice_id, abs_lw_ice),& + 'read cloud optics abs_lw_ice values') + + call handle_ncerr( nf90_close(ncid), 'ice optics file missing') + + endif ! if masterproc +#if ( defined SPMD ) + call mpibcast(g_d_eff, n_g_d, mpir8, 0, mpicom, ierr) + call mpibcast(ext_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(ssa_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(asm_sw_ice, n_g_d*nswbands, mpir8, 0, mpicom, ierr) + call mpibcast(abs_lw_ice, n_g_d*nlwbands, mpir8, 0, mpicom, ierr) +#endif + + return + +end subroutine cloud_rad_props_init + +!============================================================================== + +!ak SPAERO : hvor kalles denne fra???? Finner ikke kallet noe sted... + +subroutine cloud_rad_props_get_sw(state, pbuf, & + tau, tau_w, tau_w_g, tau_w_f,& + diagnosticindex, oldliq, oldice) + +! return totaled (across all species) layer tau, omega, g, f +! for all spectral interval for aerosols affecting the climate + + ! Arguments + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer :: pbuf(:) + integer, optional, intent(in) :: diagnosticindex ! index (if present) to radiation diagnostic information + + real(r8), intent(out) :: tau (nswbands,pcols,pver) ! aerosol extinction optical depth + real(r8), intent(out) :: tau_w (nswbands,pcols,pver) ! aerosol single scattering albedo * tau + real(r8), intent(out) :: tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w + real(r8), intent(out) :: tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w + + logical, optional, intent(in) :: oldliq,oldice + + ! Local variables + + integer :: ncol + integer :: lchnk + integer :: k, i ! lev and daycolumn indices + integer :: iswband ! sw band indices + + ! optical props for each aerosol + real(r8), pointer :: h_ext(:,:) + real(r8), pointer :: h_ssa(:,:) + real(r8), pointer :: h_asm(:,:) + real(r8), pointer :: n_ext(:) + real(r8), pointer :: n_ssa(:) + real(r8), pointer :: n_asm(:) + + ! rad properties for liquid clouds + real(r8) :: liq_tau (nswbands,pcols,pver) ! aerosol extinction optical depth + real(r8) :: liq_tau_w (nswbands,pcols,pver) ! aerosol single scattering albedo * tau + real(r8) :: liq_tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w + real(r8) :: liq_tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w + +#ifdef SPAERO + real(r8) :: xcdnc(pcols) ! CDNC multiplication factor for SP aerosols +#endif + + ! rad properties for ice clouds + real(r8) :: ice_tau (nswbands,pcols,pver) ! aerosol extinction optical depth + real(r8) :: ice_tau_w (nswbands,pcols,pver) ! aerosol single scattering albedo * tau + real(r8) :: ice_tau_w_g(nswbands,pcols,pver) ! aerosol assymetry parameter * tau * w + real(r8) :: ice_tau_w_f(nswbands,pcols,pver) ! aerosol forward scattered fraction * tau * w + + !----------------------------------------------------------------------------- + + ncol = state%ncol + lchnk = state%lchnk + + ! initialize to conditions that would cause failure + tau (:,:,:) = -100._r8 + tau_w (:,:,:) = -100._r8 + tau_w_g (:,:,:) = -100._r8 + tau_w_f (:,:,:) = -100._r8 + + ! initialize layers to accumulate od's + tau (:,1:ncol,:) = 0._r8 + tau_w (:,1:ncol,:) = 0._r8 + tau_w_g(:,1:ncol,:) = 0._r8 + tau_w_f(:,1:ncol,:) = 0._r8 + + +#ifdef SPAERO +! Denne sybrutinen (cloud_rad_props_get_sw) blir aldri kalt (i vår versjon)!!! +! så endringene nedenfor er bare for ordens og kompileringens skyld + xcdnc(1:ncol) = 1.0_r8 + call get_liquid_optics_sw(state, pbuf, xcdnc, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) +#else + call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) +#endif + + call get_ice_optics_sw (state, pbuf, ice_tau, ice_tau_w, ice_tau_w_g, ice_tau_w_f) + + tau (:,1:ncol,:) = liq_tau (:,1:ncol,:) + ice_tau (:,1:ncol,:) + tau_w (:,1:ncol,:) = liq_tau_w (:,1:ncol,:) + ice_tau_w (:,1:ncol,:) + tau_w_g(:,1:ncol,:) = liq_tau_w_g(:,1:ncol,:) + ice_tau_w_g(:,1:ncol,:) + tau_w_f(:,1:ncol,:) = liq_tau_w_f(:,1:ncol,:) + ice_tau_w_f(:,1:ncol,:) + +end subroutine cloud_rad_props_get_sw +!============================================================================== + +subroutine cloud_rad_props_get_lw(state, pbuf, cld_abs_od, diagnosticindex, oldliq, oldice, oldcloud) + +! Purpose: Compute cloud longwave absorption optical depth +! cloud_rad_props_get_lw() is called by radlw() + + ! Arguments + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer:: pbuf(:) + real(r8), intent(out) :: cld_abs_od(nlwbands,pcols,pver) ! [fraction] absorption optical depth, per layer + integer, optional, intent(in) :: diagnosticindex + logical, optional, intent(in) :: oldliq ! use old liquid optics + logical, optional, intent(in) :: oldice ! use old ice optics + logical, optional, intent(in) :: oldcloud ! use old optics for both (b4b) + + ! Local variables + + integer :: bnd_idx ! LW band index + integer :: i ! column index + integer :: k ! lev index + integer :: ncol ! number of columns + integer :: lchnk + + ! rad properties for liquid clouds + real(r8) :: liq_tau_abs_od(nlwbands,pcols,pver) ! liquid cloud absorption optical depth + + ! rad properties for ice clouds + real(r8) :: ice_tau_abs_od(nlwbands,pcols,pver) ! ice cloud absorption optical depth + + !----------------------------------------------------------------------------- + + ncol = state%ncol + lchnk = state%lchnk + + ! compute optical depths cld_absod + cld_abs_od = 0._r8 + + if(present(oldcloud))then + if(oldcloud) then + ! make diagnostic calls to these first to output ice and liq OD's + !call old_liq_get_rad_props_lw(state, pbuf, liq_tau_abs_od, oldliqwp=.false.) + !call old_ice_get_rad_props_lw(state, pbuf, ice_tau_abs_od, oldicewp=.false.) + ! This affects climate (cld_abs_od) + call oldcloud_lw(state,pbuf,cld_abs_od,oldwp=.false.) + return + endif + endif + + if(present(oldliq))then + if(oldliq) then + call old_liq_get_rad_props_lw(state, pbuf, liq_tau_abs_od, oldliqwp=.false.) + else + call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od) + endif + else + call liquid_cloud_get_rad_props_lw(state, pbuf, liq_tau_abs_od) + endif + + if(present(oldice))then + if(oldice) then + call old_ice_get_rad_props_lw(state, pbuf, ice_tau_abs_od, oldicewp=.false.) + else + call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od) + endif + else + call ice_cloud_get_rad_props_lw(state, pbuf, ice_tau_abs_od) + endif + + cld_abs_od(:,1:ncol,:) = liq_tau_abs_od(:,1:ncol,:) + ice_tau_abs_od(:,1:ncol,:) + +end subroutine cloud_rad_props_get_lw + +!============================================================================== + +subroutine get_snow_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f) + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer :: pbuf(:) + + real(r8),intent(out) :: tau (nswbands,pcols,pver) ! extinction optical depth + real(r8),intent(out) :: tau_w (nswbands,pcols,pver) ! single scattering albedo * tau + real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! assymetry parameter * tau * w + real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w + + real(r8), pointer :: icswpth(:,:), des(:,:) + + ! This does the same thing as get_ice_optics_sw, except with a different + ! water path and effective diameter. + call pbuf_get_field(pbuf, i_icswp, icswpth) + call pbuf_get_field(pbuf, i_des, des) + + call interpolate_ice_optics_sw(state%ncol, icswpth, des, tau, tau_w, & + tau_w_g, tau_w_f) + +end subroutine get_snow_optics_sw + +!============================================================================== +! Private methods +!============================================================================== + +subroutine get_ice_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f) + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer :: pbuf(:) + + real(r8),intent(out) :: tau (nswbands,pcols,pver) ! extinction optical depth + real(r8),intent(out) :: tau_w (nswbands,pcols,pver) ! single scattering albedo * tau + real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! assymetry parameter * tau * w + real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w + + real(r8), pointer :: iciwpth(:,:), dei(:,:) + + ! Get relevant pbuf fields, and interpolate optical properties from + ! the lookup tables. + call pbuf_get_field(pbuf, i_iciwp, iciwpth) + call pbuf_get_field(pbuf, i_dei, dei) + + call interpolate_ice_optics_sw(state%ncol, iciwpth, dei, tau, tau_w, & + tau_w_g, tau_w_f) + +end subroutine get_ice_optics_sw + +!============================================================================== + +subroutine interpolate_ice_optics_sw(ncol, iciwpth, dei, tau, tau_w, & + tau_w_g, tau_w_f) + + integer, intent(in) :: ncol + real(r8), intent(in) :: iciwpth(pcols,pver) + real(r8), intent(in) :: dei(pcols,pver) + + real(r8),intent(out) :: tau (nswbands,pcols,pver) ! extinction optical depth + real(r8),intent(out) :: tau_w (nswbands,pcols,pver) ! single scattering albedo * tau + real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! assymetry parameter * tau * w + real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w + + type(interp_type) :: dei_wgts + + integer :: i, k, swband + real(r8) :: ext(nswbands), ssa(nswbands), asm(nswbands) + + do k = 1,pver + do i = 1,ncol + if( iciwpth(i,k) < 1.e-80_r8 .or. dei(i,k) == 0._r8) then + ! if ice water path is too small, OD := 0 + tau (:,i,k) = 0._r8 + tau_w (:,i,k) = 0._r8 + tau_w_g(:,i,k) = 0._r8 + tau_w_f(:,i,k) = 0._r8 + else + ! for each cell interpolate to find weights in g_d_eff grid. + call lininterp_init(g_d_eff, n_g_d, dei(i:i,k), 1, & + extrap_method_bndry, dei_wgts) + ! interpolate into grid and extract radiative properties + do swband = 1, nswbands + call lininterp(ext_sw_ice(:,swband), n_g_d, & + ext(swband:swband), 1, dei_wgts) + call lininterp(ssa_sw_ice(:,swband), n_g_d, & + ssa(swband:swband), 1, dei_wgts) + call lininterp(asm_sw_ice(:,swband), n_g_d, & + asm(swband:swband), 1, dei_wgts) + end do + tau (:,i,k) = iciwpth(i,k) * ext + tau_w (:,i,k) = tau(:,i,k) * ssa + tau_w_g(:,i,k) = tau_w(:,i,k) * asm + tau_w_f(:,i,k) = tau_w_g(:,i,k) * asm + call lininterp_finish(dei_wgts) + endif + enddo + enddo + +end subroutine interpolate_ice_optics_sw + +!============================================================================== + +#ifdef SPAERO ! this seems to be called from radiation_tend only: +subroutine get_liquid_optics_sw(state, pbuf, xcdnc, tau, tau_w, tau_w_g, tau_w_f) +#else +subroutine get_liquid_optics_sw(state, pbuf, tau, tau_w, tau_w_g, tau_w_f) +#endif + + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer :: pbuf(:) + +#ifdef SPAERO + real(r8), intent (in) :: xcdnc(pcols) ! CDNC multiplication factor for SP aerosols +#endif + + real(r8),intent(out) :: tau (nswbands,pcols,pver) ! extinction optical depth + real(r8),intent(out) :: tau_w (nswbands,pcols,pver) ! single scattering albedo * tau + real(r8),intent(out) :: tau_w_g(nswbands,pcols,pver) ! asymetry parameter * tau * w + real(r8),intent(out) :: tau_w_f(nswbands,pcols,pver) ! forward scattered fraction * tau * w + + real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth +#ifdef SPAERO + real(r8), pointer, dimension(:,:) :: sp_lamc, sp_pgam +#endif + real(r8), dimension(pcols,pver) :: kext + integer i,k,swband,lchnk,ncol + + lchnk = state%lchnk + ncol = state%ncol + + + call pbuf_get_field(pbuf, i_lambda, lamc) + call pbuf_get_field(pbuf, i_mu, pgam) + call pbuf_get_field(pbuf, i_iclwp, iclwpth) +#ifdef SPAERO + call pbuf_get_field(pbuf, i_sp_lambda, sp_lamc) + call pbuf_get_field(pbuf, i_sp_mu, sp_pgam) +! do i=10,10 ! TEST -> ser ok ut, men mange 0-verdier (i starten, bare?) +! write(*,*) 'lamc(i,25), sp_lamc(i,25) = ', lamc(i,25), sp_lamc(i,25) +! write(*,*) 'pgam(i,25), sp_pgam(i,25) = ', pgam(i,25), sp_pgam(i,25) +! end do +#endif + + do k = 1,pver + do i = 1,ncol +#ifdef SPAERO + if(sp_lamc(i,k) > 0._r8) then ! This seems to be clue from microphysics of no cloud + if(xcdnc(i).gt.1._r8) then + call gam_liquid_sw(iclwpth(i,k), sp_lamc(i,k), sp_pgam(i,k), & + tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k)) + else + call gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), & + tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k)) + endif +#else + if(lamc(i,k) > 0._r8) then ! This seems to be clue from microphysics of no cloud +!#ifdef SPAERO +! if(xcdnc(i).gt.1._r8) then +! call gam_liquid_sw(iclwpth(i,k), sp_lamc(i,k), sp_pgam(i,k), & +! tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k)) +! else +! call gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), & +! tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k)) +! endif +!#else + call gam_liquid_sw(iclwpth(i,k), lamc(i,k), pgam(i,k), & + tau(1:nswbands,i,k), tau_w(1:nswbands,i,k), tau_w_g(1:nswbands,i,k), tau_w_f(1:nswbands,i,k)) +#endif + else + tau(1:nswbands,i,k) = 0._r8 + tau_w(1:nswbands,i,k) = 0._r8 + tau_w_g(1:nswbands,i,k) = 0._r8 + tau_w_f(1:nswbands,i,k) = 0._r8 + endif + enddo + enddo +!test +! write(*,*) 're_mult = ', 1, 1, re_mult(1,1) +! write(*,*) 're_mult = ', ncol, 1, re_mult(ncol,1) +! write(*,*) 're_mult = ', 1, pver, re_mult(1,pver) +! write(*,*) 'tau 1,1,1 = ', 1, 1, 1, tau(1,1,1) +! write(*,*) 'tau 1,1,p= ', 1, 1, pver, tau(1,1,pver) +! write(*,*) 'tau n,1,p= ', nswbands, 1, pver, tau(nswbands,1,pver) +!test + +end subroutine get_liquid_optics_sw + +!============================================================================== + +subroutine liquid_cloud_get_rad_props_lw(state, pbuf, abs_od) + type(physics_state), intent(in) :: state + type(physics_buffer_desc),pointer :: pbuf(:) + real(r8), intent(out) :: abs_od(nlwbands,pcols,pver) + + integer :: lchnk, ncol + real(r8), pointer, dimension(:,:) :: lamc, pgam, iclwpth + + integer lwband, i, k + + abs_od = 0._r8 + + lchnk = state%lchnk + ncol = state%ncol + + call pbuf_get_field(pbuf, i_lambda, lamc) + call pbuf_get_field(pbuf, i_mu, pgam) + call pbuf_get_field(pbuf, i_iclwp, iclwpth) + + do k = 1,pver + do i = 1,ncol + if(lamc(i,k) > 0._r8) then ! This seems to be the clue for no cloud from microphysics formulation + call gam_liquid_lw(iclwpth(i,k), lamc(i,k), pgam(i,k), abs_od(1:nlwbands,i,k)) + else + abs_od(1:nlwbands,i,k) = 0._r8 + endif + enddo + enddo + +end subroutine liquid_cloud_get_rad_props_lw +!============================================================================== + +subroutine snow_cloud_get_rad_props_lw(state, pbuf, abs_od) + type(physics_state), intent(in) :: state + type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), intent(out) :: abs_od(nlwbands,pcols,pver) + + real(r8), pointer :: icswpth(:,:), des(:,:) + + ! This does the same thing as ice_cloud_get_rad_props_lw, except with a + ! different water path and effective diameter. + call pbuf_get_field(pbuf, i_icswp, icswpth) + call pbuf_get_field(pbuf, i_des, des) + + call interpolate_ice_optics_lw(state%ncol,icswpth, des, abs_od) + +end subroutine snow_cloud_get_rad_props_lw + +!============================================================================== + +subroutine ice_cloud_get_rad_props_lw(state, pbuf, abs_od) + type(physics_state), intent(in) :: state + type(physics_buffer_desc), pointer :: pbuf(:) + real(r8), intent(out) :: abs_od(nlwbands,pcols,pver) + + real(r8), pointer :: iciwpth(:,:), dei(:,:) + + ! Get relevant pbuf fields, and interpolate optical properties from + ! the lookup tables. + call pbuf_get_field(pbuf, i_iciwp, iciwpth) + call pbuf_get_field(pbuf, i_dei, dei) + + call interpolate_ice_optics_lw(state%ncol,iciwpth, dei, abs_od) + +end subroutine ice_cloud_get_rad_props_lw + +!============================================================================== + +subroutine interpolate_ice_optics_lw(ncol, iciwpth, dei, abs_od) + + integer, intent(in) :: ncol + real(r8), intent(in) :: iciwpth(pcols,pver) + real(r8), intent(in) :: dei(pcols,pver) + + real(r8),intent(out) :: abs_od(nlwbands,pcols,pver) + + type(interp_type) :: dei_wgts + + integer :: i, k, lwband + real(r8) :: absor(nlwbands) + + do k = 1,pver + do i = 1,ncol + ! if ice water path is too small, OD := 0 + if( iciwpth(i,k) < 1.e-80_r8 .or. dei(i,k) == 0._r8) then + abs_od (:,i,k) = 0._r8 + else + ! for each cell interpolate to find weights in g_d_eff grid. + call lininterp_init(g_d_eff, n_g_d, dei(i:i,k), 1, & + extrap_method_bndry, dei_wgts) + ! interpolate into grid and extract radiative properties + do lwband = 1, nlwbands + call lininterp(abs_lw_ice(:,lwband), n_g_d, & + absor(lwband:lwband), 1, dei_wgts) + enddo + abs_od(:,i,k) = iciwpth(i,k) * absor + where(abs_od(:,i,k) > 50.0_r8) abs_od(:,i,k) = 50.0_r8 + call lininterp_finish(dei_wgts) + endif + enddo + enddo + +end subroutine interpolate_ice_optics_lw + +!============================================================================== + +subroutine gam_liquid_lw(clwptn, lamc, pgam, abs_od) + real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)? + real(r8), intent(in) :: lamc ! prognosed value of lambda for cloud + real(r8), intent(in) :: pgam ! prognosed value of mu for cloud + real(r8), intent(out) :: abs_od(1:nlwbands) + + integer :: lwband ! sw band index + + type(interp_type) :: mu_wgts + type(interp_type) :: lambda_wgts + + if (clwptn < 1.e-80_r8) then + abs_od = 0._r8 + return + endif + + call get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts) + + do lwband = 1, nlwbands + call lininterp(abs_lw_liq(:,:,lwband), nmu, nlambda, & + abs_od(lwband:lwband), 1, mu_wgts, lambda_wgts) + enddo + + abs_od = clwptn * abs_od + + call lininterp_finish(mu_wgts) + call lininterp_finish(lambda_wgts) + +end subroutine gam_liquid_lw + +!============================================================================== + +subroutine gam_liquid_sw(clwptn, lamc, pgam, tau, tau_w, tau_w_g, tau_w_f) + real(r8), intent(in) :: clwptn ! cloud water liquid path new (in cloud) (in g/m^2)? + real(r8), intent(in) :: lamc ! prognosed value of lambda for cloud + real(r8), intent(in) :: pgam ! prognosed value of mu for cloud + real(r8), intent(out) :: tau(1:nswbands), tau_w(1:nswbands), tau_w_f(1:nswbands), tau_w_g(1:nswbands) + + integer :: swband ! sw band index + + real(r8) :: ext(nswbands), ssa(nswbands), asm(nswbands) + + type(interp_type) :: mu_wgts + type(interp_type) :: lambda_wgts + + if (clwptn < 1.e-80_r8) then + tau = 0._r8 + tau_w = 0._r8 + tau_w_g = 0._r8 + tau_w_f = 0._r8 + return + endif + + call get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts) + + do swband = 1, nswbands + call lininterp(ext_sw_liq(:,:,swband), nmu, nlambda, & + ext(swband:swband), 1, mu_wgts, lambda_wgts) + call lininterp(ssa_sw_liq(:,:,swband), nmu, nlambda, & + ssa(swband:swband), 1, mu_wgts, lambda_wgts) + call lininterp(asm_sw_liq(:,:,swband), nmu, nlambda, & + asm(swband:swband), 1, mu_wgts, lambda_wgts) + enddo + + ! compute radiative properties + tau = clwptn * ext + tau_w = tau * ssa + tau_w_g = tau_w * asm + tau_w_f = tau_w_g * asm + + call lininterp_finish(mu_wgts) + call lininterp_finish(lambda_wgts) + +end subroutine gam_liquid_sw + +!============================================================================== + +subroutine get_mu_lambda_weights(lamc, pgam, mu_wgts, lambda_wgts) + real(r8), intent(in) :: lamc ! prognosed value of lambda for cloud + real(r8), intent(in) :: pgam ! prognosed value of mu for cloud + ! Output interpolation weights. Caller is responsible for freeing these. + type(interp_type), intent(out) :: mu_wgts + type(interp_type), intent(out) :: lambda_wgts + + integer :: ilambda + real(r8) :: g_lambda_interp(nlambda) + + ! Make interpolation weights for mu. + ! (Put pgam in a temporary array for this purpose.) + call lininterp_init(g_mu, nmu, [pgam], 1, extrap_method_bndry, mu_wgts) + + ! Use mu weights to interpolate to a row in the lambda table. + do ilambda = 1, nlambda + call lininterp(g_lambda(:,ilambda), nmu, & + g_lambda_interp(ilambda:ilambda), 1, mu_wgts) + end do + + ! Make interpolation weights for lambda. + call lininterp_init(g_lambda_interp, nlambda, [lamc], 1, & + extrap_method_bndry, lambda_wgts) + +end subroutine get_mu_lambda_weights + +!============================================================================== + +end module cloud_rad_props diff --git a/src/NorESM/micro_mg2_0.F90 b/src/NorESM/micro_mg2_0.F90 index 97567f491e..3db85fb331 100644 --- a/src/NorESM/micro_mg2_0.F90 +++ b/src/NorESM/micro_mg2_0.F90 @@ -99,13 +99,19 @@ module micro_mg2_0 ! ......................................................................... ! NOTE: List of all inputs/outputs passed through the call/subroutine statement ! for micro_mg_tend is given below at the start of subroutine micro_mg_tend. -!--------------------------------------------------------------------------------- +!----------------------------------------------------------------------------------- +! Sept. 2019 A. Kirkevåg Added MACv2 SP-aerosol code wrt. Twomey effect (SPAERO) +!----------------------------------------------------------------------------------- ! Procedures required: ! 1) An implementation of the gamma function (if not intrinsic). ! 2) saturation vapor pressure and specific humidity over water ! 3) svp over ice +!ak+ +#include +!ak- + #ifndef HAVE_GAMMA_INTRINSICS use shr_spfn_mod, only: gamma => shr_spfn_gamma #endif @@ -377,6 +383,10 @@ subroutine micro_mg_tend ( & prain, prodsnow, & cmeout, deffi, & pgamrad, lamcrad, & +#ifdef SPAERO + sp_xcdnc, & + sp_pgamrad, sp_lamcrad, & +#endif qsout, dsout, & lflx, iflx, & rflx, sflx, qrout, & @@ -525,6 +535,10 @@ subroutine micro_mg_tend ( & real(r8), intent(out) :: deffi(mgncol,nlev) ! ice effective diameter for optics (radiation) (micron) real(r8), intent(out) :: pgamrad(mgncol,nlev) ! ice gamma parameter for optics (radiation) (no units) real(r8), intent(out) :: lamcrad(mgncol,nlev) ! slope of droplet distribution for optics (radiation) (1/m) +#ifdef SPAERO + real(r8), intent(out) :: sp_pgamrad(mgncol,nlev) ! ice gamma parameter for SP optics (radiation) (no units) + real(r8), intent(out) :: sp_lamcrad(mgncol,nlev) ! slope of droplet distribution for SP optics (radiation) (1/m) +#endif real(r8), intent(out) :: qsout(mgncol,nlev) ! snow mixing ratio (kg/kg) real(r8), intent(out) :: dsout(mgncol,nlev) ! snow diameter (m) real(r8), intent(out) :: lflx(mgncol,nlev+1) ! grid-box average liquid condensate flux (kg m^-2 s^-1) @@ -852,6 +866,10 @@ subroutine micro_mg_tend ( & ! dummies for in-cloud variables real(r8) :: dumc(mgncol,nlev) ! qc real(r8) :: dumnc(mgncol,nlev) ! nc +#ifdef SPAERO + real(r8) :: sp_dumnc(mgncol,nlev) + real(r8), intent(in) :: sp_xcdnc(mgncol) +#endif real(r8) :: dumi(mgncol,nlev) ! qi real(r8) :: dumni(mgncol,nlev) ! ni real(r8) :: dumr(mgncol,nlev) ! rain mixing ratio @@ -1130,6 +1148,10 @@ subroutine micro_mg_tend ( & effc = 10._r8 lamcrad = 0._r8 pgamrad = 0._r8 +#ifdef SPAERO + sp_lamcrad = 0._r8 + sp_pgamrad = 0._r8 +#endif effc_fn = 10._r8 effi = 25._r8 sadice = 0._r8 @@ -2188,7 +2210,7 @@ subroutine micro_mg_tend ( & ! switch for specification of droplet and crystal number if (nccons) then dumnc(i,k)=ncnst/rho(i,k) - end if + end if ! switch for specification of cloud ice number if (nicons) then @@ -2881,6 +2903,9 @@ subroutine micro_mg_tend ( & dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k) dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k) dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k) +#ifdef SPAERO + sp_dumnc(i,k) = dumnc(i,k) * sp_xcdnc(i) +#endif dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k) dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat,0._r8)/precip_frac(i,k) @@ -2973,6 +2998,13 @@ subroutine micro_mg_tend ( & dum = dumnc(i,k) +#ifdef SPAERO +! Extra call with SP aerosol to obtain new droplet size distrobution for +! radiative transfer (only). Does not change the results when xcdnc=1.0_r8. +! NB: dumnc er en in-out parameter i kallet under, mens dumc ikke er det ! + call size_dist_param_liq(mg_liq_props, dumc(i,k), sp_dumnc(i,k), rho(i,k), & + sp_pgamrad(i,k), sp_lamcrad(i,k)) +#endif call size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), & pgam(i,k), lamc(i,k)) @@ -2986,12 +3018,15 @@ subroutine micro_mg_tend ( & nctend(i,k)=(dumnc(i,k)*lcldm(i,k)-nc(i,k))/deltat end if - effc(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8 +#ifdef SPAERO + effc(i,k) = (sp_pgamrad(i,k)+3._r8)/sp_lamcrad(i,k)/2._r8*1.e6_r8 +#else + effc(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8 +#endif !assign output fields for shape here lamcrad(i,k)=lamc(i,k) pgamrad(i,k)=pgam(i,k) - ! recalculate effective radius for constant number, in order to separate ! first and second indirect effects !====================================== diff --git a/src/NorESM/micro_mg_cam.F90 b/src/NorESM/micro_mg_cam.F90 index bbae68e378..bc9f5908f1 100644 --- a/src/NorESM/micro_mg_cam.F90 +++ b/src/NorESM/micro_mg_cam.F90 @@ -61,7 +61,13 @@ module micro_mg_cam ! "post_proc%process_and_unpack", which sets every single field that was ! added with post_proc%add_field. ! -!--------------------------------------------------------------------------------- +!----------------------------------------------------------------------------------- +! Sept. 2019 A. Kirkevåg Added MACv2 SP-aerosol code wrt. Twomey effect (SPAERO) +!----------------------------------------------------------------------------------- + +!ak+ +#include +!ak- use shr_kind_mod, only: r8=>shr_kind_r8 use spmd_utils, only: masterproc @@ -162,8 +168,13 @@ module micro_mg_cam rel_idx, & dei_idx, & mu_idx, & - prer_evap_idx, & + prer_evap_idx, & lambdac_idx, & +#ifdef SPAERO + sp_mu_idx, & + sp_lambdac_idx, & + sp_rel_idx, & +#endif iciwpst_idx, & iclwpst_idx, & des_idx, & @@ -214,6 +225,9 @@ module micro_mg_cam cld_idx = -1, & concld_idx = -1, & qsatfac_idx = -1 +#ifdef SPAERO + integer :: sp_xcdnc_idx =-1 +#endif ! Pbuf fields needed for subcol_SILHS integer :: & @@ -452,6 +466,9 @@ subroutine micro_mg_cam_register call pbuf_add_field('SADICE', 'physpkg',dtype_r8,(/pcols,pver/), sadice_idx) call pbuf_add_field('SADSNOW', 'physpkg',dtype_r8,(/pcols,pver/), sadsnow_idx) call pbuf_add_field('REL', 'physpkg',dtype_r8,(/pcols,pver/), rel_idx) +#ifdef SPAERO + call pbuf_add_field('SPREL', 'physpkg',dtype_r8,(/pcols,pver/), sp_rel_idx) +#endif ! Mitchell ice effective diameter for radiation call pbuf_add_field('DEI', 'physpkg',dtype_r8,(/pcols,pver/), dei_idx) @@ -459,7 +476,12 @@ subroutine micro_mg_cam_register call pbuf_add_field('MU', 'physpkg',dtype_r8,(/pcols,pver/), mu_idx) ! Size distribution shape parameter for radiation call pbuf_add_field('LAMBDAC', 'physpkg',dtype_r8,(/pcols,pver/), lambdac_idx) - +#ifdef SPAERO + ! Size distribution shape parameter for radiation + call pbuf_add_field('SP_MU', 'physpkg',dtype_r8,(/pcols,pver/), sp_mu_idx) + ! Size distribution shape parameter for radiation + call pbuf_add_field('SP_LAMBDAC', 'physpkg',dtype_r8,(/pcols,pver/), sp_lambdac_idx) +#endif ! Stratiform only in cloud ice water path for radiation call pbuf_add_field('ICIWPST', 'physpkg',dtype_r8,(/pcols,pver/), iciwpst_idx) ! Stratiform in cloud liquid water path for radiation @@ -530,6 +552,9 @@ subroutine micro_mg_cam_register call pbuf_register_subcol('SADICE', 'micro_mg_cam_register', sadice_idx) call pbuf_register_subcol('SADSNOW', 'micro_mg_cam_register', sadsnow_idx) call pbuf_register_subcol('REL', 'micro_mg_cam_register', rel_idx) +#ifdef SPAERO + call pbuf_register_subcol('SPREL', 'micro_mg_cam_register', sp_rel_idx) +#endif ! Mitchell ice effective diameter for radiation call pbuf_register_subcol('DEI', 'micro_mg_cam_register', dei_idx) @@ -537,7 +562,12 @@ subroutine micro_mg_cam_register call pbuf_register_subcol('MU', 'micro_mg_cam_register', mu_idx) ! Size distribution shape parameter for radiation call pbuf_register_subcol('LAMBDAC', 'micro_mg_cam_register', lambdac_idx) - +#ifdef SPAERO + ! Size distribution shape parameter for radiation + call pbuf_register_subcol('SP_MU', 'micro_mg_cam_register', sp_mu_idx) + ! Size distribution shape parameter for radiation + call pbuf_register_subcol('SP_LAMBDAC', 'micro_mg_cam_register', sp_lambdac_idx) +#endif ! Stratiform only in cloud ice water path for radiation call pbuf_register_subcol('ICIWPST', 'micro_mg_cam_register', iciwpst_idx) ! Stratiform in cloud liquid water path for radiation @@ -873,6 +903,9 @@ subroutine micro_mg_cam_init(pbuf2d) call addfld ('AWNC', (/ 'lev' /), 'A', 'm-3', 'Average cloud water number conc' ) call addfld ('AWNI', (/ 'lev' /), 'A', 'm-3', 'Average cloud ice number conc' ) call addfld ('AREL', (/ 'lev' /), 'A', 'Micron', 'Average droplet effective radius' ) +#ifdef SPAERO + call addfld ('SPAREL', (/ 'lev' /), 'A', 'Micron', 'Average droplet effective radius SP' ) +#endif call addfld ('AREI', (/ 'lev' /), 'A', 'Micron', 'Average ice effective radius' ) ! Frequency arrays for above call addfld ('FREQL', (/ 'lev' /), 'A', 'fraction', 'Fractional occurrence of liquid' ) @@ -888,8 +921,7 @@ subroutine micro_mg_cam_init(pbuf2d) call addfld ('FCTI', horiz_only, 'A', 'fraction', 'Fractional occurrence of cloud top ice' ) !++IH !For comparing to Bernartz CDNC concentrations -!akc6 call addfld ('ACTNL_B ', horiz_only, 'A', 'Micron ', 'Average Cloud Top droplet number (Bennartz)' ) - call addfld ('ACTNL_B ', horiz_only, 'A', 'm-3', 'Average Cloud Top droplet number (Bennartz)' ) + call addfld ('ACTNL_B ', horiz_only, 'A', 'm-3', 'Average Cloud Top droplet number (Bennartz)' ) call addfld ('FCTL_B ', horiz_only, 'A', 'fraction', 'Fractional occurrence of cloud top liquid (Bennartz)' ) !ak6+ call addfld ('CCN_B ', horiz_only, 'A', 'm-3', 'Average Cloud Top liquid CCN (Bennartz)' ) @@ -984,6 +1016,9 @@ subroutine micro_mg_cam_init(pbuf2d) call add_default ('ADSNOW ', 1, ' ') call add_default ('AREI ', 1, ' ') call add_default ('AREL ', 1, ' ') +#ifdef SPAERO + call add_default ('SPAREL ', 1, ' ') +#endif call add_default ('AWNC ', 1, ' ') call add_default ('AWNI ', 1, ' ') call add_default ('CDNUMC ', 1, ' ') @@ -1142,6 +1177,10 @@ subroutine micro_mg_cam_init(pbuf2d) nrain_idx = pbuf_get_index('NRAIN', ierr) nsnow_idx = pbuf_get_index('NSNOW', ierr) +#ifdef SPAERO + sp_xcdnc_idx = pbuf_get_index('SP_XCDNC', ierr) +#endif + ! fields for heterogeneous freezing frzimm_idx = pbuf_get_index('FRZIMM', ierr) frzcnt_idx = pbuf_get_index('FRZCNT', ierr) @@ -1281,6 +1320,11 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle real(r8), pointer :: dei(:,:) ! Ice effective diameter (meters) (AG: microns?) real(r8), pointer :: mu(:,:) ! Size distribution shape parameter for radiation real(r8), pointer :: lambdac(:,:) ! Size distribution slope parameter for radiation +#ifdef SPAERO + real(r8), pointer :: sp_mu(:,:) ! Size distribution shape parameter for radiation + real(r8), pointer :: sp_lambdac(:,:) ! Size distribution slope parameter for radiation + real(r8), pointer :: sp_xcdnc(:) ! CDNC modification factor from MACv2SP +#endif real(r8), pointer :: des(:,:) ! Snow effective diameter (m) real(r8) :: rho(state%psetcols,pver) @@ -1526,11 +1570,19 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle real(r8), target :: packed_qcrat(mgncol,nlev) real(r8), target :: packed_rel(mgncol,nlev) +#ifdef SPAERO + real(r8), target :: packed_sp_rel(mgncol,nlev) +#endif real(r8), target :: packed_rei(mgncol,nlev) real(r8), target :: packed_sadice(mgncol,nlev) real(r8), target :: packed_sadsnow(mgncol,nlev) real(r8), target :: packed_lambdac(mgncol,nlev) real(r8), target :: packed_mu(mgncol,nlev) +#ifdef SPAERO + real(r8), target :: packed_sp_lambdac(mgncol,nlev) + real(r8), target :: packed_sp_mu(mgncol,nlev) + real(r8), target :: packed_sp_xcdnc(mgncol) +#endif real(r8), target :: packed_des(mgncol,nlev) real(r8), target :: packed_dei(mgncol,nlev) @@ -1625,6 +1677,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle real(r8) :: liqcldf(state%psetcols,pver) ! Liquid cloud fraction (combined into cloud) real(r8), pointer :: rel(:,:) ! Liquid effective drop radius (microns) +#ifdef SPAERO + real(r8), pointer :: sp_rel(:,:) ! Liquid effective drop radius (microns) +#endif real(r8), pointer :: rei(:,:) ! Ice effective drop size (microns) real(r8), pointer :: sadice(:,:) ! Ice surface area density (cm2/cm3) real(r8), pointer :: sadsnow(:,:) ! Snow surface area density (cm2/cm3) @@ -1650,6 +1705,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle ! Averaging arrays for effective radius and number.... real(r8) :: efiout_grid(pcols,pver) real(r8) :: efcout_grid(pcols,pver) +#ifdef SPAERO + real(r8) :: spefcout_grid(pcols,pver) +#endif real(r8) :: ncout_grid(pcols,pver) real(r8) :: niout_grid(pcols,pver) real(r8) :: freqi_grid(pcols,pver) @@ -1715,6 +1773,11 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle real(r8), pointer :: lambdac_grid(:,:) real(r8), pointer :: mu_grid(:,:) +#ifdef SPAERO + real(r8), pointer :: sp_lambdac_grid(:,:) + real(r8), pointer :: sp_mu_grid(:,:) + real(r8), pointer :: sp_rel_grid(:,:) +#endif real(r8), pointer :: rel_grid(:,:) real(r8), pointer :: rei_grid(:,:) real(r8), pointer :: sadice_grid(:,:) @@ -1727,6 +1790,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle real(r8) :: liqcldf_grid(pcols,pver) real(r8) :: qsout_grid(pcols,pver) real(r8) :: ncic_grid(pcols,pver) +#ifdef SPAERO + real(r8) :: sp_ncic_grid(pcols,pver) +#endif real(r8) :: niic_grid(pcols,pver) real(r8) :: rel_fn_grid(pcols,pver) ! Ice effective drop size at fixed number (indirect effect) (microns) - on grid real(r8) :: qrout_grid(pcols,pver) @@ -1872,6 +1938,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle if (qsatfac_idx > 0) call pbuf_get_field(pbuf, qsatfac_idx, qsatfac, col_type=col_type, copy_if_needed=use_subcol_microp) +#ifdef SPAERO + call pbuf_get_field(pbuf, sp_xcdnc_idx, sp_xcdnc, col_type=col_type) +#endif + !----------------------- ! These physics buffer fields are calculated and set in this parameterization ! If subcolumns is turned on, then these fields will be calculated on a subcolumn grid, otherwise they will be a normal grid @@ -1888,6 +1958,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle call pbuf_get_field(pbuf, dei_idx, dei, col_type=col_type) call pbuf_get_field(pbuf, mu_idx, mu, col_type=col_type) call pbuf_get_field(pbuf, lambdac_idx, lambdac, col_type=col_type) +#ifdef SPAERO + call pbuf_get_field(pbuf, sp_mu_idx, sp_mu, col_type=col_type) + call pbuf_get_field(pbuf, sp_lambdac_idx, sp_lambdac, col_type=col_type) +#endif call pbuf_get_field(pbuf, des_idx, des, col_type=col_type) call pbuf_get_field(pbuf, ls_flxprc_idx, mgflxprc, col_type=col_type) call pbuf_get_field(pbuf, ls_flxsnw_idx, mgflxsnw, col_type=col_type) @@ -1899,6 +1973,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle call pbuf_get_field(pbuf, iclwpst_idx, iclwpst, col_type=col_type) call pbuf_get_field(pbuf, icswp_idx, icswp, col_type=col_type) call pbuf_get_field(pbuf, rel_idx, rel, col_type=col_type) +#ifdef SPAERO + call pbuf_get_field(pbuf, sp_rel_idx, sp_rel, col_type=col_type) +#endif call pbuf_get_field(pbuf, rei_idx, rei, col_type=col_type) call pbuf_get_field(pbuf, sadice_idx, sadice, col_type=col_type) call pbuf_get_field(pbuf, sadsnow_idx, sadsnow, col_type=col_type) @@ -1940,6 +2017,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle call pbuf_get_field(pbuf, dei_idx, dei_grid) call pbuf_get_field(pbuf, mu_idx, mu_grid) call pbuf_get_field(pbuf, lambdac_idx, lambdac_grid) +#ifdef SPAERO + call pbuf_get_field(pbuf, sp_mu_idx, sp_mu_grid) + call pbuf_get_field(pbuf, sp_lambdac_idx, sp_lambdac_grid) +#endif call pbuf_get_field(pbuf, des_idx, des_grid) call pbuf_get_field(pbuf, ls_flxprc_idx, mgflxprc_grid) call pbuf_get_field(pbuf, ls_flxsnw_idx, mgflxsnw_grid) @@ -1951,6 +2032,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle call pbuf_get_field(pbuf, iclwpst_idx, iclwpst_grid) call pbuf_get_field(pbuf, icswp_idx, icswp_grid) call pbuf_get_field(pbuf, rel_idx, rel_grid) +#ifdef SPAERO + call pbuf_get_field(pbuf, sp_rel_idx, sp_rel_grid) +#endif call pbuf_get_field(pbuf, rei_idx, rei_grid) call pbuf_get_field(pbuf, sadice_idx, sadice_grid) call pbuf_get_field(pbuf, sadsnow_idx, sadsnow_grid) @@ -2200,6 +2284,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle ! the value from the last substep, which is what "accum_null" does. call post_proc%add_field(p(rel), p(packed_rel), & fillvalue=10._r8, accum_method=accum_null) +#ifdef SPAERO + call post_proc%add_field(p(sp_rel), p(packed_sp_rel), & + fillvalue=10._r8, accum_method=accum_null) +#endif call post_proc%add_field(p(rei), p(packed_rei), & fillvalue=25._r8, accum_method=accum_null) call post_proc%add_field(p(sadice), p(packed_sadice), & @@ -2210,6 +2298,12 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle accum_method=accum_null) call post_proc%add_field(p(mu), p(packed_mu), & accum_method=accum_null) +#ifdef SPAERO + call post_proc%add_field(p(sp_lambdac), p(packed_sp_lambdac), & + accum_method=accum_null) + call post_proc%add_field(p(sp_mu), p(packed_sp_mu), & + accum_method=accum_null) +#endif call post_proc%add_field(p(des), p(packed_des), & accum_method=accum_null) call post_proc%add_field(p(dei), p(packed_dei), & @@ -2236,6 +2330,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle packed_naai = packer%pack(naai) packed_npccn = packer%pack(npccn) +#ifdef SPAERO + packed_sp_xcdnc = packer%pack(sp_xcdnc) +#endif + allocate(packed_rndst(mgncol,nlev,size(rndst, 3))) packed_rndst = packer%pack(rndst) @@ -2341,6 +2439,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle packed_prain, packed_prodsnow, & packed_cmeout, packed_dei, & packed_mu, packed_lambdac, & +#ifdef SPAERO + packed_sp_xcdnc, & + packed_sp_mu, packed_sp_lambdac, & +#endif packed_qsout, packed_des, & packed_cflx, packed_iflx, & packed_rflx, packed_sflx, packed_qrout, & @@ -2648,6 +2750,11 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle ! as they are reset before being used, so it would be a needless calculation lambdac_grid => lambdac mu_grid => mu +#ifdef SPAERO + sp_lambdac_grid => sp_lambdac + sp_mu_grid => sp_mu + sp_rel_grid => sp_rel +#endif rel_grid => rel rei_grid => rei sadice_grid => sadice @@ -2758,10 +2865,10 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle end if ! Effective radius for cloud liquid, fixed number. + mu_grid = 0._r8 lambdac_grid = 0._r8 rel_fn_grid = 10._r8 - ncic_grid = 1.e8_r8 call size_dist_param_liq(mg_liq_props, icwmrst_grid(:ngrdcol,top_lev:), & @@ -2774,6 +2881,55 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle lambdac_grid(:ngrdcol,top_lev:)/2._r8 * 1.e6_r8 end where +!------------------------------------------------------------------------- +#ifdef SPAERO + ! Calculating cloud droplet size distribution and effective radii for + ! the MACv2SP aerosols + + ! Effective radius for cloud liquid, and size parameters + ! sp_mu_grid and sp_lambdac_grid. + sp_mu_grid = 0._r8 + sp_lambdac_grid = 0._r8 + sp_rel_grid = 10._r8 + + ! Calculate sp_ncic on the grid + sp_ncic_grid = 1.e8_r8 + ncic_grid(:ngrdcol,top_lev:) = nc_grid(:ngrdcol,top_lev:) / & + max(mincld,liqcldf_grid(:ngrdcol,top_lev:)) + + do k = top_lev, pver + sp_ncic_grid(:,k) = sp_xcdnc(:) * ncic_grid(:,k) + end do + + call size_dist_param_liq(mg_liq_props, icwmrst_grid(:ngrdcol,top_lev:), & + sp_ncic_grid(:ngrdcol,top_lev:), rho_grid(:ngrdcol,top_lev:), & + sp_mu_grid(:ngrdcol,top_lev:), sp_lambdac_grid(:ngrdcol,top_lev:)) + + where (icwmrst_grid(:ngrdcol,top_lev:) >= qsmall) + sp_rel_grid(:ngrdcol,top_lev:) = & + (sp_mu_grid(:ngrdcol,top_lev:) + 3._r8) / & + sp_lambdac_grid(:ngrdcol,top_lev:)/2._r8 * 1.e6_r8 + elsewhere + ! Deal with the fact that size_dist_param_liq sets mu_grid to -100 + ! wherever there is no cloud. + sp_mu_grid(:ngrdcol,top_lev:) = 0._r8 + end where + + ! Limiters for low cloud fraction. + do k = top_lev, pver + do i = 1, ngrdcol + ! Convert snow effective diameter to microns + if ( ast_grid(i,k) < 1.e-4_r8 ) then + sp_mu_grid(i,k) = mucon + sp_lambdac_grid(i,k) = (mucon + 1._r8)/dcon + end if + end do + end do + +#endif +!------------------------------------------------------------------------- + + ! Effective radius for cloud liquid, and size parameters ! mu_grid and lambdac_grid. mu_grid = 0._r8 @@ -3005,6 +3161,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle ! Averaging for new output fields efcout_grid = 0._r8 +#ifdef SPAERO + spefcout_grid = 0._r8 +#endif efiout_grid = 0._r8 ncout_grid = 0._r8 niout_grid = 0._r8 @@ -3020,6 +3179,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle do i = 1, ngrdcol if ( liqcldf_grid(i,k) > 0.01_r8 .and. icwmrst_grid(i,k) > 5.e-5_r8 ) then efcout_grid(i,k) = rel_grid(i,k) * liqcldf_grid(i,k) +#ifdef SPAERO + spefcout_grid(i,k) = sp_rel_grid(i,k) * liqcldf_grid(i,k) +#endif ncout_grid(i,k) = icwnc_grid(i,k) * liqcldf_grid(i,k) freql_grid(i,k) = liqcldf_grid(i,k) icwmrst_grid_out(i,k) = icwmrst_grid(i,k) @@ -3298,6 +3460,9 @@ subroutine micro_mg_cam_tend_pack(state, ptend, dtime, pbuf, mgncol, mgcols, nle call outfld('VPRCO', vprco_grid, pcols, lchnk) call outfld('RACAU', racau_grid, pcols, lchnk) call outfld('AREL', efcout_grid, pcols, lchnk) +#ifdef SPAERO + call outfld('SPAREL', spefcout_grid, pcols, lchnk) +#endif call outfld('AREI', efiout_grid, pcols, lchnk) call outfld('AWNC' , ncout_grid, pcols, lchnk) call outfld('AWNI' , niout_grid, pcols, lchnk) diff --git a/src/NorESM/physpkg.F90 b/src/NorESM/physpkg.F90 index 7bc17bdba5..7d79ff3e9b 100644 --- a/src/NorESM/physpkg.F90 +++ b/src/NorESM/physpkg.F90 @@ -9,6 +9,8 @@ module physpkg ! 2005-10-17 B. Eaton Add contents of inti.F90 to phys_init(). Add ! initialization of grid info in phys_state. ! Nov 2010 A. Gettelman Put micro/macro physics into separate routines + ! Sept 2019 A. Kirkevåg Added MACv2 SP-aerosol code wrt. Twomey effect + ! (SPAERO) !----------------------------------------------------------------------- #ifdef OSLO_AERO @@ -89,6 +91,10 @@ module physpkg integer :: snow_sh_idx = 0 integer :: dlfzm_idx = 0 ! detrained convective cloud water mixing ratio. +#ifdef SPAERO + integer :: sp_xcdnc_idx = 0 +#endif + !======================================================================= contains !======================================================================= @@ -199,6 +205,11 @@ subroutine phys_register call pbuf_add_field('CLDNCINI', 'physpkg', dtype_r8, (/pcols,pver/), cldncini_idx) call pbuf_add_field('CLDNIINI', 'physpkg', dtype_r8, (/pcols,pver/), cldniini_idx) !AL + +#ifdef SPAERO + call pbuf_add_field('SP_XCDNC', 'physpkg', dtype_r8, (/pcols/), sp_xcdnc_idx) +#endif + ! check energy package call check_energy_register @@ -1785,6 +1796,7 @@ subroutine tphysbc (ztodt, state, & integer kcomp ! mode number (1-14) #endif + ! physics buffer fields to compute tendencies for stratiform package integer itim_old, ifld real(r8), pointer, dimension(:,:) :: cld ! cloud fraction @@ -1800,6 +1812,11 @@ subroutine tphysbc (ztodt, state, & real(r8), pointer, dimension(:,:) :: cldncini real(r8), pointer, dimension(:,:) :: cldniini !AL + +#ifdef SPAERO + real(r8), pointer, dimension(:) :: sp_xcdnc ! CDNC modification factor from MACv2SP +#endif + real(r8), pointer, dimension(:,:,:) :: fracis ! fraction of transported species that are insoluble real(r8), pointer :: dlfzm(:,:) ! ZM detrained convective cloud water mixing ratio. @@ -1919,6 +1936,11 @@ subroutine tphysbc (ztodt, state, & call pbuf_get_field(pbuf, cldniini_idx, cldniini) !AL +#ifdef SPAERO + ifld = pbuf_get_index('SP_XCDNC') + call pbuf_get_field(pbuf, ifld, sp_xcdnc, start=(/1/), kount=(/pcols/) ) +#endif + ifld = pbuf_get_index('DTCORE') call pbuf_get_field(pbuf, ifld, dtcore, start=(/1,1,itim_old/), kount=(/pcols,pver,1/) ) @@ -2403,9 +2425,9 @@ subroutine tphysbc (ztodt, state, & ! call outfld('RNEWD8 ',rnewdry8,pcols,lchnk) ! call outfld('RNEWD9 ',rnewdry9,pcols,lchnk) ! call outfld('RNEWD10 ',rnewdry10,pcols,lchnk) -!! call outfld('RNEWD11 ',rnewdry11,pcols,lchnk) ! always = 0.0118 -!! call outfld('RNEWD13 ',rnewdry13,pcols,lchnk) ! always = 0.04 -!! call outfld('RNEWD14 ',rnewdry14,pcols,lchnk) ! always = 0.04 +!! call outfld('RNEWD11 ',rnewdry11,pcols,lchnk) ! constant +!! call outfld('RNEWD13 ',rnewdry13,pcols,lchnk) ! constant +!! call outfld('RNEWD14 ',rnewdry14,pcols,lchnk) ! constant ! call outfld('RNEW1 ',rnew1,pcols,lchnk) ! call outfld('RNEW2 ',rnew2,pcols,lchnk) ! call outfld('RNEW4 ',rnew4,pcols,lchnk) @@ -2427,9 +2449,9 @@ subroutine tphysbc (ztodt, state, & ! call outfld('LOGSIG8 ',logsig8,pcols,lchnk) ! call outfld('LOGSIG9 ',logsig9,pcols,lchnk) ! call outfld('LOGSIG10',logsig10,pcols,lchnk) -!! call outfld('LOGSIG11',logsig11,pcols,lchnk) ! always = 0.2553 -!! call outfld('LOGSIG13',logsig13,pcols,lchnk) ! always = 0.2553 -!! call outfld('LOGSIG14',logsig14,pcols,lchnk) ! always = 0.2553 +!! call outfld('LOGSIG11',logsig11,pcols,lchnk) ! constant +!! call outfld('LOGSIG13',logsig13,pcols,lchnk) ! constant +!! call outfld('LOGSIG14',logsig14,pcols,lchnk) ! constant #endif ! aerocom #endif ! dirind @@ -2485,7 +2507,12 @@ subroutine tphysbc (ztodt, state, & call radiation_tend( & + +#ifdef SPAERO + state, ptend, pbuf, cam_out, cam_in, net_flx, sp_xcdnc) +#else state, ptend, pbuf, cam_out, cam_in, net_flx) +#endif ! Set net flux used by spectral dycores do i=1,ncol diff --git a/src/physics/cam_oslo/mo_simple_plumes_v1.F90 b/src/physics/cam_oslo/mo_simple_plumes_v1.F90 new file mode 100644 index 0000000000..34719997d8 --- /dev/null +++ b/src/physics/cam_oslo/mo_simple_plumes_v1.F90 @@ -0,0 +1,399 @@ +!> +!! +!! @brief Module MO_SIMPLE_PLUMES: provides anthropogenic aerosol optical properties as a function of lat, lon +!! height, time, and wavelength +!! +!! @remarks +!! +!! @author Bjorn Stevens, Stephanie Fiedler and Karsten Peters MPI-Met, Hamburg (v1 release 2016-11-10) +!! +!! @change-log: +!! - 2016-12-05: beta release (BS, SF and KP, MPI-Met) +!! - 2016-09-28: revised representation of Twomey effect (SF, MPI-Met) +!! - 2015-09-28: bug fixes (SF, MPI-Met) +!! - 2016-10-12: revised maximum longitudinal extent of European plume (KP, SF, MPI-Met) +!! $ID: n/a$ +!! +!! @par Origin +!! Based on code originally developed at the MPI-Met by Karsten Peters, Bjorn Stevens, Stephanie Fiedler +!! and Stefan Kinne with input from Thorsten Mauritsen and Robert Pincus +!! +!! @par Copyright +!! +!!! +!! Adapted for use in NorESM (P. Räisänen, FMI, 15 Feb 2017) +!! +!! Adapted for use in NorESM2 (A. Kirkevåg, September 2019) +!! + +MODULE MO_SIMPLE_PLUMES + + USE netcdf + USE shr_kind_mod, only: r8=>shr_kind_r8 + USE ppgrid, only: pcols +!ak+ + USE oslo_control, only: oslo_getopts, dir_string_length +!ak- + + IMPLICIT NONE + + INTEGER, PARAMETER :: & + nplumes = 9 ,& !< Number of plumes + nfeatures = 2 ,& !< Number of features per plume + ntimes = 52 ,& !< Number of times resolved per year (52 => weekly resolution) + nyears = 251 !< Number of years of available forcing + + LOGICAL, SAVE :: & + sp_initialized = .FALSE. !< parameter determining whether input needs to be read + + REAL(r8) :: & + plume_lat (nplumes) ,& !< latitude of plume center (AOD maximum) + plume_lon (nplumes) ,& !< longitude of plume center (AOD maximum) + beta_a (nplumes) ,& !< parameter a for beta function vertical profile + beta_b (nplumes) ,& !< parameter b for beta function vertical profile + aod_spmx (nplumes) ,& !< anthropogenic AOD maximum at 550 for plumes + aod_fmbg (nplumes) ,& !< anthropogenic AOD at 550 for fine-mode natural background (idealized to mimic Twomey effect) + asy550 (nplumes) ,& !< asymmetry parameter at 550nm for plume + ssa550 (nplumes) ,& !< single scattering albedo at 550nm for plume + angstrom (nplumes) ,& !< Angstrom parameter for plume + sig_lon_E (nfeatures,nplumes) ,& !< Eastward extent of plume feature + sig_lon_W (nfeatures,nplumes) ,& !< Westward extent of plume feature + sig_lat_E (nfeatures,nplumes) ,& !< Southward extent of plume feature + sig_lat_W (nfeatures,nplumes) ,& !< Northward extent of plume feature + theta (nfeatures,nplumes) ,& !< Rotation angle of plume feature + ftr_weight (nfeatures,nplumes) ,& !< Feature weights + time_weight (nfeatures,nplumes) ,& !< Time weights + time_weight_bg (nfeatures,nplumes) ,& !< as time_weight but for natural background in Twomey effect + year_weight (nyears,nplumes) ,& !< Yearly weight for plume + ann_cycle (nfeatures,ntimes,nplumes) !< annual cycle for plume feature + + PUBLIC sp_aop_profile + +CONTAINS + ! + ! ------------------------------------------------------------------------------------------------------------------------ + ! SP_SETUP: This subroutine should be called at initialization to read the netcdf data that describes the simple plume + ! climatology. The information needs to be either read by each processor or distributed to processors. + ! + SUBROUTINE sp_setup + ! + ! ---------- + ! + INTEGER :: iret, ncid, DimID, VarID, xdmy +!ak CHARACTER(LEN=100) :: infile + ! + ! ---------- + ! + ! So far, hard-coded: +!ak+ +! infile = "/stornext/field/users/raisanen/RECIA/MACv2.0-SP_v1.nc" +! iret = nf90_open(infile, NF90_NOWRITE, ncid) + character(len=dir_string_length) :: aerotab_table_dir + call oslo_getopts(aerotab_table_dir_out = aerotab_table_dir) + iret = nf90_open(trim(aerotab_table_dir)//'/MACv2.0-SP_v1.nc', NF90_NOWRITE, ncid) +!ak- + + IF (iret /= NF90_NOERR) STOP 'NetCDF File not opened' + ! + ! read dimensions and make sure file conforms to expected size + ! + iret = nf90_inq_dimid(ncid, "plume_number" , DimId) + iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) + IF (xdmy /= nplumes) STOP 'NetCDF improperly dimensioned -- plume_number' + + iret = nf90_inq_dimid(ncid, "plume_feature", DimId) + iret = nf90_inquire_dimension(ncid, DimId, len = xdmy) + IF (xdmy /= nfeatures) STOP 'NetCDF improperly dimensioned -- plume_feature' + + iret = nf90_inq_dimid(ncid, "year_fr" , DimId) + iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) + IF (xdmy /= ntimes) STOP 'NetCDF improperly dimensioned -- year_fr' + + iret = nf90_inq_dimid(ncid, "years" , DimId) + iret = nf90_inquire_dimension(ncid, DimID, len = xdmy) + IF (xdmy /= nyears) STOP 'NetCDF improperly dimensioned -- years' + ! + ! read variables that define the simple plume climatology + ! + iret = nf90_inq_varid(ncid, "plume_lat", VarId) + iret = nf90_get_var(ncid, VarID, plume_lat(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat' + iret = nf90_inq_varid(ncid, "plume_lon", VarId) + iret = nf90_get_var(ncid, VarID, plume_lon(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lon' + iret = nf90_inq_varid(ncid, "beta_a" , VarId) + iret = nf90_get_var(ncid, VarID, beta_a(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_a' + iret = nf90_inq_varid(ncid, "beta_b" , VarId) + iret = nf90_get_var(ncid, VarID, beta_b(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading beta_b' + iret = nf90_inq_varid(ncid, "aod_spmx" , VarId) + iret = nf90_get_var(ncid, VarID, aod_spmx(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_spmx' + iret = nf90_inq_varid(ncid, "aod_fmbg" , VarId) + iret = nf90_get_var(ncid, VarID, aod_fmbg(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading aod_fmbg' + iret = nf90_inq_varid(ncid, "ssa550" , VarId) + iret = nf90_get_var(ncid, VarID, ssa550(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ssa550' + iret = nf90_inq_varid(ncid, "asy550" , VarId) + iret = nf90_get_var(ncid, VarID, asy550(:) , start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading asy550' + iret = nf90_inq_varid(ncid, "angstrom" , VarId) + iret = nf90_get_var(ncid, VarID, angstrom(:), start=(/1/),count=(/nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading angstrom' + + iret = nf90_inq_varid(ncid, "sig_lat_W" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lat_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_W' + iret = nf90_inq_varid(ncid, "sig_lat_E" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lat_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lat_E' + iret = nf90_inq_varid(ncid, "sig_lon_E" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lon_E(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_E' + iret = nf90_inq_varid(ncid, "sig_lon_W" , VarId) + iret = nf90_get_var(ncid, VarID, sig_lon_W(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading sig_lon_W' + iret = nf90_inq_varid(ncid, "theta" , VarId) + iret = nf90_get_var(ncid, VarID, theta(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading theta' + iret = nf90_inq_varid(ncid, "ftr_weight" , VarId) + iret = nf90_get_var(ncid, VarID, ftr_weight(:,:) , start=(/1,1/),count=(/nfeatures,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading plume_lat' + iret = nf90_inq_varid(ncid, "year_weight" , VarId) + iret = nf90_get_var(ncid, VarID, year_weight(:,:) , start=(/1,1/),count=(/nyears,nplumes /)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading year_weight' + iret = nf90_inq_varid(ncid, "ann_cycle" , VarId) + iret = nf90_get_var(ncid, VarID, ann_cycle(:,:,:) , start=(/1,1,1/),count=(/nfeatures,ntimes,nplumes/)) + IF (iret /= NF90_NOERR) STOP 'NetCDF Error reading ann_cycle' + + iret = nf90_close(ncid) + + sp_initialized = .TRUE. + + RETURN + END SUBROUTINE sp_setup + ! + ! ------------------------------------------------------------------------------------------------------------------------ + ! SET_TIME_WEIGHT: The simple plume model assumes that meteorology constrains plume shape and that only source strength + ! influences the amplitude of a plume associated with a given source region. This routine retrieves the temporal weights + ! for the plumes. Each plume feature has its own temporal weights which varies yearly. The annual cycle is indexed by + ! week in the year and superimposed on the yearly mean value of the weight. + ! + SUBROUTINE set_time_weight(year_fr) + ! + ! ---------- + ! + REAL(r8), INTENT(IN) :: & + year_fr !< Fractional Year (1850.0 - 2100.99) + + INTEGER :: & + iyear ,& !< Integer year values between 1 and 156 (1850-2100) + iweek ,& !< Integer index (between 1 and ntimes); for ntimes=52 this corresponds to weeks (roughly) + iplume ! plume number + ! + ! ---------- + ! + iyear = FLOOR(year_fr) - 1849 + iweek = FLOOR((year_fr - FLOOR(year_fr)) * ntimes) + 1 + + IF ((iweek > ntimes) .OR. (iweek < 1) .OR. (iyear > nyears) .OR. (iyear < 1)) STOP 'Time out of bounds in set_time_weight' + DO iplume=1,nplumes + time_weight(1,iplume) = year_weight(iyear,iplume) * ann_cycle(1,iweek,iplume) + time_weight(2,iplume) = year_weight(iyear,iplume) * ann_cycle(2,iweek,iplume) + time_weight_bg(1,iplume) = ann_cycle(1,iweek,iplume) + time_weight_bg(2,iplume) = ann_cycle(2,iweek,iplume) + END DO + + RETURN + END SUBROUTINE set_time_weight + ! + ! ------------------------------------------------------------------------------------------------------------------------ + ! SP_AOP_PROFILE: This subroutine calculates the simple plume aerosol and cloud active optical properties based on the + ! the simple plume fit to the MPI Aerosol Climatology (Version 2). It sums over nplumes to provide a profile of aerosol + ! optical properties on a host models vertical grid. + ! + SUBROUTINE sp_aop_profile ( & + nlevels ,ncol ,lambda ,oro ,lon ,lat , & + year_fr ,z ,dz ,dNovrN ,aod_prof ,ssa_prof , & + asy_prof ) + ! + ! ---------- + ! + INTEGER, INTENT(IN) :: & + nlevels, & !< number of levels + ncol !< number of columns + + REAL(r8), INTENT(IN) :: & + lambda, & !< wavelength + year_fr, & !< Fractional Year (1903.0 is the 0Z on the first of January 1903, Gregorian) + oro(pcols), & !< orographic height (m) + lon(pcols), & !< longitude + lat(pcols), & !< latitude + z (pcols,nlevels), & !< height above sea-level (m) + dz(pcols,nlevels) !< level thickness (difference between half levels) (m) + + REAL(r8), INTENT(OUT) :: & + dNovrN(pcols) , & !< anthropogenic increase in cloud drop number concentration (factor) + aod_prof(pcols,nlevels), & !< profile of aerosol optical depth + ssa_prof(pcols,nlevels), & !< profile of single scattering albedo + asy_prof(pcols,nlevels) !< profile of asymmetry parameter + + INTEGER :: iplume, icol, k + + REAL(r8) :: & + eta(pcols,nlevels), & !< normalized height (by 15 km) + z_beta(pcols,nlevels), & !< profile for scaling column optical depth + prof(pcols,nlevels), & !< scaled profile (by beta function) + beta_sum(pcols), & !< vertical sum of beta function + ssa(pcols), & !< single scattering albedo + asy(pcols), & !< asymmetry parameter + cw_an(pcols), & !< column weight for simple plume (anthropogenic) AOD at 550 nm + cw_bg(pcols), & !< column weight for fine-mode natural background AOD at 550 nm + caod_sp(pcols), & !< column simple plume anthropogenic AOD at 550 nm + caod_bg(pcols), & !< column fine-mode natural background AOD at 550 nm + a_plume1, & !< gaussian longitude factor for feature 1 + a_plume2, & !< gaussian longitude factor for feature 2 + b_plume1, & !< gaussian latitude factor for feature 1 + b_plume2, & !< gaussian latitude factor for feature 2 + delta_lat, & !< latitude offset + delta_lon, & !< longitude offset + delta_lon_t, & !< threshold for maximum longitudinal plume extent used in transition from 360 to 0 degrees + lon1, & !< rotated longitude for feature 1 + lat1, & !< rotated latitude for feature 2 + lon2, & !< rotated longitude for feature 1 + lat2, & !< rotated latitude for feature 2 + f1, & !< contribution from feature 1 + f2, & !< contribution from feature 2 + f3, & !< contribution from feature 1 in natural background of Twomey effect + f4, & !< contribution from feature 2 in natural background of Twomey effect + aod_550, & !< aerosol optical depth at 550nm + aod_lmd, & !< aerosol optical depth at input wavelength + lfactor !< factor to compute wavelength dependence of optical properties + ! + ! ---------- + ! + ! initialize input data (by calling setup at first instance) + ! + IF (.NOT.sp_initialized) CALL sp_setup + ! + ! get time weights + ! + CALL set_time_weight(year_fr) + ! + ! initialize variables, including output + ! + DO k=1,nlevels + DO icol=1,ncol + aod_prof(icol,k) = 0.0 + ssa_prof(icol,k) = 0.0 + asy_prof(icol,k) = 0.0 + z_beta(icol,k) = MERGE(1.0, 0.0, z(icol,k) >= oro(icol)) + eta(icol,k) = MAX(0.0,MIN(1.0,z(icol,k)/15000.)) + END DO + END DO + DO icol=1,ncol + dNovrN(icol) = 1.0 + caod_sp(icol) = 0.0 + caod_bg(icol) = 0.02 + END DO + ! + ! sum contribution from plumes to construct composite profiles of aerosol optical properties + ! + DO iplume=1,nplumes + ! + ! calculate vertical distribution function from parameters of beta distribution + ! + DO icol=1,ncol + beta_sum(icol) = 0. + END DO + DO k=1,nlevels + DO icol=1,ncol + prof(icol,k) = (eta(icol,k)**(beta_a(iplume)-1.) * (1.-eta(icol,k))**(beta_b(iplume)-1.)) * dz(icol,k) + beta_sum(icol) = beta_sum(icol) + prof(icol,k) + END DO + END DO + DO k=1,nlevels + DO icol=1,ncol + prof(icol,k) = ( prof(icol,k) / beta_sum(icol) ) * z_beta(icol,k) + END DO + END DO + ! + ! calculate plume weights + ! + DO icol=1,ncol + ! + ! get plume-center relative spatial parameters for specifying amplitude of plume at given lat and lon + ! + delta_lat = lat(icol) - plume_lat(iplume) + delta_lon = lon(icol) - plume_lon(iplume) + delta_lon_t = MERGE (260., 180., iplume == 1) + delta_lon = MERGE ( delta_lon-SIGN(360.,delta_lon) , delta_lon , ABS(delta_lon) > delta_lon_t) + + a_plume1 = 0.5 / (MERGE(sig_lon_E(1,iplume), sig_lon_W(1,iplume), delta_lon > 0)**2) + b_plume1 = 0.5 / (MERGE(sig_lat_E(1,iplume), sig_lat_W(1,iplume), delta_lon > 0)**2) + a_plume2 = 0.5 / (MERGE(sig_lon_E(2,iplume), sig_lon_W(2,iplume), delta_lon > 0)**2) + b_plume2 = 0.5 / (MERGE(sig_lat_E(2,iplume), sig_lat_W(2,iplume), delta_lon > 0)**2) + ! + ! adjust for a plume specific rotation which helps match plume state to climatology. + ! + lon1 = COS(theta(1,iplume))*(delta_lon) + SIN(theta(1,iplume))*(delta_lat) + lat1 = - SIN(theta(1,iplume))*(delta_lon) + COS(theta(1,iplume))*(delta_lat) + lon2 = COS(theta(2,iplume))*(delta_lon) + SIN(theta(2,iplume))*(delta_lat) + lat2 = - SIN(theta(2,iplume))*(delta_lon) + COS(theta(2,iplume))*(delta_lat) + ! + ! calculate contribution to plume from its different features, to get a column weight for the anthropogenic + ! (cw_an) and the fine-mode natural background aerosol (cw_bg) + ! + f1 = time_weight(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) + f2 = time_weight(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) + f3 = time_weight_bg(1,iplume) * ftr_weight(1,iplume) * EXP(-1.* (a_plume1 * ((lon1)**2) + (b_plume1 * ((lat1)**2)))) + f4 = time_weight_bg(2,iplume) * ftr_weight(2,iplume) * EXP(-1.* (a_plume2 * ((lon2)**2) + (b_plume2 * ((lat2)**2)))) + + cw_an(icol) = f1 * aod_spmx(iplume) + f2 * aod_spmx(iplume) + cw_bg(icol) = f3 * aod_fmbg(iplume) + f4 * aod_fmbg(iplume) + ! + ! calculate wavelength-dependent scattering properties + ! + lfactor = MIN(1.0,700.0/lambda) + ssa(icol) = (ssa550(iplume) * lfactor**4) / ((ssa550(iplume) * lfactor**4) + ((1-ssa550(iplume)) * lfactor)) + asy(icol) = asy550(iplume) * SQRT(lfactor) + END DO + ! + ! distribute plume optical properties across its vertical profile weighting by optical depth and scaling for + ! wavelength using the angstrom parameter. + ! + lfactor = EXP(-angstrom(iplume) * LOG(lambda/550.0)) + DO k=1,nlevels + DO icol = 1,ncol + aod_550 = prof(icol,k) * cw_an(icol) + aod_lmd = aod_550 * lfactor + caod_sp(icol) = caod_sp(icol) + aod_550 + caod_bg(icol) = caod_bg(icol) + prof(icol,k) * cw_bg(icol) + asy_prof(icol,k) = asy_prof(icol,k) + aod_lmd * ssa(icol) * asy(icol) + ssa_prof(icol,k) = ssa_prof(icol,k) + aod_lmd * ssa(icol) + aod_prof(icol,k) = aod_prof(icol,k) + aod_lmd + END DO + END DO + END DO + ! + ! complete optical depth weighting + ! + DO k=1,nlevels + DO icol = 1,ncol + asy_prof(icol,k) = MERGE(asy_prof(icol,k)/ssa_prof(icol,k), 0.0, ssa_prof(icol,k) > TINY(1.)) + ssa_prof(icol,k) = MERGE(ssa_prof(icol,k)/aod_prof(icol,k), 1.0, aod_prof(icol,k) > TINY(1.)) + END DO + END DO + ! + ! calculate effective radius normalization (divisor) factor + ! + DO icol=1,ncol + dNovrN(icol) = LOG((1000.0 * (caod_sp(icol) + caod_bg(icol))) + 1.0)/LOG((1000.0 * caod_bg(icol)) + 1.0) + END DO + + RETURN + END SUBROUTINE sp_aop_profile + +END MODULE MO_SIMPLE_PLUMES diff --git a/src/physics/cam_oslo/pmxsub.F90 b/src/physics/cam_oslo/pmxsub.F90 index 58687750cb..1900124891 100644 --- a/src/physics/cam_oslo/pmxsub.F90 +++ b/src/physics/cam_oslo/pmxsub.F90 @@ -53,7 +53,7 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & integer, intent(in) :: lchnk ! chunk identifier integer, intent(in) :: ncol ! number of atmospheric columns real(r8), intent(in) :: coszrs(pcols) ! Cosine solar zenith angle - real(r8), intent(in) :: pint(pcols,pverp) ! Model interface pressures (10*Pa) + real(r8), intent(in) :: pint(pcols,pverp) ! Model interface pressures (10*Pa) !!!!!!!! real(r8), intent(in) :: pmid(pcols,pver) ! Model level pressures (Pa) real(r8), intent(in) :: t(pcols,pver) ! Model level temperatures (K) real(r8), intent(in) :: cld(pcols,pver) ! cloud fraction @@ -85,7 +85,7 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & !---------------------------Local variables----------------------------- ! integer i, k, ib, icol, mplus10 - integer iloop + integer iloop ! Loop index for testing of CPU time of parts of the code logical daylight(pcols) ! SW calculations also at (polar) night in interpol* if daylight=.true. real(r8) aodvisvolc(pcols) ! AOD vis for CMIP6 volcanic aerosol @@ -97,7 +97,6 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & !tst ! real(r8) aodvis3d(pcols,pver) ! 3D AOD in VIS !tst - real(r8) deltah_km(pcols,pver) ! Layer thickness, unit km !akc6 real(r8) deltah, airmass(pcols,pver) @@ -971,13 +970,11 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & per_tau_w(i,k,ib)=per_tau(i,k,ib)*max(min(ssatot(i,k,14-ib),0.999999_r8),1.e-6_r8) per_tau_w_g(i,k,ib)=per_tau_w(i,k,ib)*asymtot(i,k,14-ib) per_tau_w_f(i,k,ib)=per_tau_w_g(i,k,ib)*asymtot(i,k,14-ib) -!tst ! if(ib.eq.4.and.k.eq.pver.and.i.eq.1) then ! write(*,*) 'per_tau =', per_tau(i,k,ib) ! write(*,*) 'per_tau_w =', per_tau_w(i,k,ib) ! write(*,*) 'per_tau_w_g =', per_tau_w_g(i,k,ib) ! endif -!tst end do end do ib=14 @@ -1009,7 +1006,7 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & enddo enddo -! Adding also the volcanic contribution (CMIP6), which is also using +! Adding the volcanic contribution (CMIP6), which is also using ! AeroTab band numbering, so that a remapping is required here do ib=1,nlwbands volc_balw(1:ncol,1:pver,ib) = volc_ext_earth(:ncol,1:pver,ib) & @@ -1179,9 +1176,10 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & do k=1,pver ! Layer thickness, unit km, and layer airmass, unit kg/m2 deltah=deltah_km(icol,k) -!akc6 airmass(icol,k)=1.e3_r8*deltah*rhoda(icol,k) +!akc6+ airmass(icol,k)=1.e3_r8*deltah*rhoda(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) +!akc6- ! Optical depths at ca. 550 nm (0.442-0.625um) all aerosols !tst ! aodvis3d(icol,k)=betotvis(icol,k)*deltah @@ -1976,6 +1974,7 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & !tst ! call outfld('DOD5503D',dod5503d,pcols,lchnk) !tst +!- call outfld('DOD5503D',dod5503d,pcols,lchnk) !- call outfld('ABS5503D',abs5503d,pcols,lchnk) !- call outfld('D443_SS ',dod4403d_ss ,pcols,lchnk) !- call outfld('D443_DU ',dod4403d_dust,pcols,lchnk) @@ -2392,7 +2391,6 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & call outfld('PM2P5 ',c_pm25 ,pcols,lchnk) call outfld('MMRPM2P5',mmr_pm25,pcols,lchnk) call outfld('MMRPM1 ',mmr_pm1 ,pcols,lchnk) - call outfld('MMRPM2P5_SRF',mmr_pm25(:pcols,pver),pcols,lchnk) !akc6- ! total (all r) dry concentrations (ug/m3) and loadings (mg/m2) call outfld('DLOAD_MI',dload_mi,pcols,lchnk) @@ -2433,7 +2431,6 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & !akc6 call outfld('AIRMASSL',airmassl,pcols,lchnk) call outfld('AIRMASSL',airmassl,pcols,lchnk) call outfld('AIRMASS ',airmass,pcols,lchnk) !akc6 - !c_er3d ! effective dry radii (um) in each layer ! call outfld('ERLT053D',erlt053d,pcols,lchnk) @@ -2448,7 +2445,7 @@ subroutine pmxsub(lchnk, ncol, pint, pmid, coszrs, state, t, cld, qm1, Nnatk, & !000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000 -! Extra AeroCom diagnostics requiring table look-ups with RH = constant +! Extra AeroCom diagnostics requiring table look-ups with RH = constant. #ifdef AEROCOM_INSITU irfmax=6 diff --git a/src/physics/cam_oslo/radconstants.F90 b/src/physics/cam_oslo/radconstants.F90 new file mode 100644 index 0000000000..45c735dfaa --- /dev/null +++ b/src/physics/cam_oslo/radconstants.F90 @@ -0,0 +1,281 @@ +module radconstants + +! This module contains constants that are specific to the radiative transfer +! code used in the RRTMG model. +! +! Alf Kirkevåg Sept. 2019: Added "wav_sp" for using the MACv2SP "simple plume" +! aerosol climatology, based on modified NorESM1 code from P. Räisänen (2016). + +#ifdef OSLO_AERO +#include +#endif + +use shr_kind_mod, only: r8 => shr_kind_r8 +use cam_abortutils, only: endrun + +implicit none +private +save + +! SHORTWAVE DATA + +! number of shorwave spectral intervals +integer, parameter, public :: nswbands = 14 +integer, parameter, public :: nbndsw = 14 + +! Wavenumbers of band boundaries +! +! Note: Currently rad_solar_var extends the lowest band down to +! 100 cm^-1 if it is too high to cover the far-IR. Any changes meant +! to affect IR solar variability should take note of this. + +real(r8), parameter :: wavenum_low(nbndsw) = & ! in cm^-1 + (/2600._r8, 3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, & + 8050._r8,12850._r8,16000._r8,22650._r8,29000._r8,38000._r8, 820._r8/) +real(r8), parameter :: wavenum_high(nbndsw) = & ! in cm^-1 + (/3250._r8, 4000._r8, 4650._r8, 5150._r8, 6150._r8, 7700._r8, 8050._r8, & + 12850._r8,16000._r8,22650._r8,29000._r8,38000._r8,50000._r8, 2600._r8/) + +#ifdef SPAERO + ! Band "representative wavenumber" for using MACv2-SP aerosol climatology, + ! assumed to be the midband wavenumbers, (wavenum_low+wavenum_high)/2.0: +real(r8), public, parameter :: wav_sp(nbndsw) = & ! in cm^-1 + (/2925._r8, 3650._r8, 4325._r8, 4900._r8, 5650._r8, 6925._r8, 7875._r8, & + 10450._r8,14425._r8,19325._r8,25825._r8,33500._r8,44000._r8, 1710._r8/) +integer iw +#endif + +! Solar irradiance at 1 A.U. in W/m^2 assumed by radiation code +! Rescaled so that sum is precisely 1368.22 and fractional amounts sum to 1.0 +real(r8), parameter :: solar_ref_band_irradiance(nbndsw) = & + (/ & + 12.11_r8, 20.3600000000001_r8, 23.73_r8, & + 22.43_r8, 55.63_r8, 102.93_r8, 24.29_r8, & + 345.74_r8, 218.19_r8, 347.20_r8, & + 129.49_r8, 50.15_r8, 3.08_r8, 12.89_r8 & + /) + +! None of the following comment appears to be the case any more? This +! should be reevalutated and/or removed. + +! rrtmg (coarse) reference solar flux in rrtmg is initialized as the following +! reference data inside rrtmg seems to indicate 1366.44 instead +! This data references 1366.442114152342 +!real(r8), parameter :: solar_ref_band_irradiance(nbndsw) = & +! (/ & +! 12.10956827000000_r8, 20.36508467999999_r8, 23.72973826333333_r8, & +! 22.42769644333333_r8, 55.62661262000000_r8, 102.9314315544444_r8, 24.29361887666667_r8, & +! 345.7425138000000_r8, 218.1870300666667_r8, 347.1923147000001_r8, & +! 129.4950181200000_r8, 48.37217043000000_r8, 3.079938997898001_r8, 12.88937733000000_r8 & +! /) +! Kurucz (fine) reference would seem to imply the following but the above values are from rrtmg_sw_init +! (/12.109559, 20.365097, 23.729752, 22.427697, 55.626622, 102.93142, 24.293593, & +! 345.73655, 218.18416, 347.18406, 129.49407, 50.147238, 3.1197130, 12.793834 /) + +! These are indices to the band for diagnostic output +integer, parameter, public :: idx_sw_diag = 10 ! index to sw visible band +integer, parameter, public :: idx_nir_diag = 8 ! index to sw near infrared (778-1240 nm) band +integer, parameter, public :: idx_uv_diag = 11 ! index to sw uv (345-441 nm) band + +integer, parameter, public :: rrtmg_sw_cloudsim_band = 9 ! rrtmg band for .67 micron + +! Number of evenly spaced intervals in rh +! The globality of this mesh may not be necessary +! Perhaps it could be specific to the aerosol +! But it is difficult to see how refined it must be +! for lookup. This value was found to be sufficient +! for Sulfate and probably necessary to resolve the +! high variation near rh = 1. Alternative methods +! were found to be too slow. +! Optimal approach would be for cam to specify size of aerosol +! based on each aerosol's characteristics. Radiation +! should know nothing about hygroscopic growth! +integer, parameter, public :: nrh = 1000 + +! LONGWAVE DATA + +! These are indices to the band for diagnostic output +integer, parameter, public :: idx_lw_diag = 7 ! index to (H20 window) LW band + +integer, parameter, public :: rrtmg_lw_cloudsim_band = 6 ! rrtmg band for 10.5 micron + +! number of lw bands +integer, parameter, public :: nlwbands = 16 +integer, parameter, public :: nbndlw = 16 + +real(r8), parameter :: wavenumber1_longwave(nlwbands) = &! Longwave spectral band limits (cm-1) + (/ 10._r8, 350._r8, 500._r8, 630._r8, 700._r8, 820._r8, 980._r8, 1080._r8, & + 1180._r8, 1390._r8, 1480._r8, 1800._r8, 2080._r8, 2250._r8, 2390._r8, 2600._r8 /) + +real(r8), parameter :: wavenumber2_longwave(nlwbands) = &! Longwave spectral band limits (cm-1) + (/ 350._r8, 500._r8, 630._r8, 700._r8, 820._r8, 980._r8, 1080._r8, 1180._r8, & + 1390._r8, 1480._r8, 1800._r8, 2080._r8, 2250._r8, 2390._r8, 2600._r8, 3250._r8 /) + +!These can go away when old camrt disappears +! Index of volc. abs., H2O non-window +integer, public, parameter :: idx_LW_H2O_NONWND=1 +! Index of volc. abs., H2O window +integer, public, parameter :: idx_LW_H2O_WINDOW=2 +! Index of volc. cnt. abs. 0500--0650 cm-1 +integer, public, parameter :: idx_LW_0500_0650=3 +! Index of volc. cnt. abs. 0650--0800 cm-1 +integer, public, parameter :: idx_LW_0650_0800=4 +! Index of volc. cnt. abs. 0800--1000 cm-1 +integer, public, parameter :: idx_LW_0800_1000=5 +! Index of volc. cnt. abs. 1000--1200 cm-1 +integer, public, parameter :: idx_LW_1000_1200=6 +! Index of volc. cnt. abs. 1200--2000 cm-1 +integer, public, parameter :: idx_LW_1200_2000=7 + +! GASES TREATED BY RADIATION (line spectrae) + +! gasses required by radiation +integer, public, parameter :: gasnamelength = 5 +integer, public, parameter :: nradgas = 8 +character(len=gasnamelength), public, parameter :: gaslist(nradgas) & + = (/'H2O ','O3 ', 'O2 ', 'CO2 ', 'N2O ', 'CH4 ', 'CFC11', 'CFC12'/) + +! what is the minimum mass mixing ratio that can be supported by radiation implementation? +real(r8), public, parameter :: minmmr(nradgas) & + = epsilon(1._r8) + +! Length of "optics type" string specified in optics files. +integer, parameter, public :: ot_length = 32 + +public :: rad_gas_index + +public :: get_number_sw_bands, & + get_sw_spectral_boundaries, & + get_lw_spectral_boundaries, & + get_ref_solar_band_irrad, & + get_ref_total_solar_irrad, & + get_solar_band_fraction_irrad + +!#ifdef SPAERO + ! This would be a safer way to define wav_sp, but would have to be moved... + ! Band "representative wavenumber" for using MACv2-SP aerosol climatology, + ! assumed to be the midband wavenumbers, (wavenum_low+wavenum_high)/2.0: +!real(r8), public :: wav_sp(nbndsw) ! in cm^-1 +! do iw = 1, nbndsw +! wav_sp(iw) = (wavenum_low(iw)+wavenum_high(iw))/2._r8 +! end do +!#endif + +contains +!------------------------------------------------------------------------------ +subroutine get_solar_band_fraction_irrad(fractional_irradiance) + ! provide Solar Irradiance for each band in RRTMG + + ! fraction of solar irradiance in each band + real(r8), intent(out) :: fractional_irradiance(1:nswbands) + real(r8) :: tsi ! total solar irradiance + + tsi = sum(solar_ref_band_irradiance) + fractional_irradiance = solar_ref_band_irradiance / tsi + +end subroutine get_solar_band_fraction_irrad +!------------------------------------------------------------------------------ +subroutine get_ref_total_solar_irrad(tsi) + ! provide Total Solar Irradiance assumed by RRTMG + + real(r8), intent(out) :: tsi + + tsi = sum(solar_ref_band_irradiance) + +end subroutine get_ref_total_solar_irrad +!------------------------------------------------------------------------------ +subroutine get_ref_solar_band_irrad( band_irrad ) + + ! solar irradiance in each band (W/m^2) + real(r8), intent(out) :: band_irrad(nswbands) + + band_irrad = solar_ref_band_irradiance + +end subroutine get_ref_solar_band_irrad +!------------------------------------------------------------------------------ +subroutine get_number_sw_bands(number_of_bands) + + ! number of solar (shortwave) bands in the rrtmg code + integer, intent(out) :: number_of_bands + + number_of_bands = nswbands + +end subroutine get_number_sw_bands + +!------------------------------------------------------------------------------ +subroutine get_lw_spectral_boundaries(low_boundaries, high_boundaries, units) + ! provide spectral boundaries of each longwave band + + real(r8), intent(out) :: low_boundaries(nlwbands), high_boundaries(nlwbands) + character(*), intent(in) :: units ! requested units + + select case (units) + case ('inv_cm','cm^-1','cm-1') + low_boundaries = wavenumber1_longwave + high_boundaries = wavenumber2_longwave + case('m','meter','meters') + low_boundaries = 1.e-2_r8/wavenumber2_longwave + high_boundaries = 1.e-2_r8/wavenumber1_longwave + case('nm','nanometer','nanometers') + low_boundaries = 1.e7_r8/wavenumber2_longwave + high_boundaries = 1.e7_r8/wavenumber1_longwave + case('um','micrometer','micrometers','micron','microns') + low_boundaries = 1.e4_r8/wavenumber2_longwave + high_boundaries = 1.e4_r8/wavenumber1_longwave + case('cm','centimeter','centimeters') + low_boundaries = 1._r8/wavenumber2_longwave + high_boundaries = 1._r8/wavenumber1_longwave + case default + call endrun('get_lw_spectral_boundaries: spectral units not acceptable'//units) + end select + +end subroutine get_lw_spectral_boundaries + +!------------------------------------------------------------------------------ +subroutine get_sw_spectral_boundaries(low_boundaries, high_boundaries, units) + ! provide spectral boundaries of each shortwave band + + real(r8), intent(out) :: low_boundaries(nswbands), high_boundaries(nswbands) + character(*), intent(in) :: units ! requested units + + select case (units) + case ('inv_cm','cm^-1','cm-1') + low_boundaries = wavenum_low + high_boundaries = wavenum_high + case('m','meter','meters') + low_boundaries = 1.e-2_r8/wavenum_high + high_boundaries = 1.e-2_r8/wavenum_low + case('nm','nanometer','nanometers') + low_boundaries = 1.e7_r8/wavenum_high + high_boundaries = 1.e7_r8/wavenum_low + case('um','micrometer','micrometers','micron','microns') + low_boundaries = 1.e4_r8/wavenum_high + high_boundaries = 1.e4_r8/wavenum_low + case('cm','centimeter','centimeters') + low_boundaries = 1._r8/wavenum_high + high_boundaries = 1._r8/wavenum_low + case default + call endrun('rad_constants.F90: spectral units not acceptable'//units) + end select + +end subroutine get_sw_spectral_boundaries + +!------------------------------------------------------------------------------ +integer function rad_gas_index(gasname) + + ! return the index in the gaslist array of the specified gasname + + character(len=*),intent(in) :: gasname + integer :: igas + + rad_gas_index = -1 + do igas = 1, nradgas + if (trim(gaslist(igas)).eq.trim(gasname)) then + rad_gas_index = igas + return + endif + enddo + call endrun ("rad_gas_index: can not find gas with name "//gasname) +end function rad_gas_index + +end module radconstants diff --git a/src/physics/cam_oslo/radiation.F90 b/src/physics/cam_oslo/radiation.F90 index 71e5697419..b413e2fd71 100644 --- a/src/physics/cam_oslo/radiation.F90 +++ b/src/physics/cam_oslo/radiation.F90 @@ -130,11 +130,11 @@ module radiation ! 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. -!#ifdef RFMIPIRF -! logical :: spectralflux = .true. ! calculate fluxes (up and down) per band. -!#else +#ifdef RFMIPIRF + logical :: spectralflux = .true. ! calculate fluxes (up and down) per band. +#else logical :: spectralflux = .false. ! calculate fluxes (up and down) per band. -!#endif +#endif ! Physics buffer indices integer :: qrs_idx = 0 @@ -184,9 +184,9 @@ subroutine radiation_readnl(nlfile) namelist /radiation_nl/ iradsw, iradlw, irad_always, & use_rad_dt_cosz, spectralflux -!#ifdef RFMIPIRF -! spectralflux = .true. ! calculate fluxes (up and down) per band. -!#endif +#ifdef RFMIPIRF + spectralflux = .true. ! calculate fluxes (up and down) per band. +#endif !----------------------------------------------------------------------------- @@ -261,12 +261,16 @@ subroutine radiation_register ! If the namelist has been configured for preserving the spectral fluxes, then create ! physics buffer variables to store the results. ! legg til #ifndef RFMIPIRF her ogsaa?! +#ifndef RFMIPIRF if (spectralflux) then +#endif 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) +#ifndef RFMIPIRF end if +#endif call rad_data_register() @@ -506,10 +510,6 @@ subroutine radiation_init(pbuf2d) 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') -!#ifdef AEROFFL -! 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 (history_amwg) then call add_default('SOLIN'//diag(icall), 1, ' ') @@ -698,11 +698,11 @@ end subroutine radiation_read_restart !=============================================================================== subroutine radiation_tend( & -!#ifdef SPAERO -! state, ptend, pbuf, cam_out, cam_in, net_flx, xcdnc, rd_out) -!#else +#ifdef SPAERO + state, ptend, pbuf, cam_out, cam_in, net_flx, xcdnc, rd_out) +#else state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out) -!#endif +#endif !----------------------------------------------------------------------- ! @@ -713,8 +713,10 @@ subroutine radiation_tend( & ! 2007-12-27 M. Iacono Modify to use CAM cloud optical properties with rrtmg. ! ! - ! 2019-05-06 A. Kirkevåg: Changes for testing the - ! "simple plumes" aerosols, based on NorESM1 code P. Räisänen. + ! Spet. 2019 A. Kirkevåg: The code has been modified to facilitate use of the + ! MACv2-SP (simple plumes) prescribed aerosol optics and CDNC scaling (for use in RFMIP) + ! (SPAERO), based on NorESM1 code from P. Räisänen. + ! !----------------------------------------------------------------------- use phys_grid, only: get_rlat_all_p, get_rlon_all_p @@ -750,10 +752,10 @@ subroutine radiation_tend( & use constituents, only: pcnst use oslo_control, only: oslo_getopts use physics_buffer, only: pbuf_get_index -!#ifdef SPAERO -! use time_manager, only: get_curr_date -! use physconst, only: rair -!#endif +#ifdef SPAERO + use time_manager, only: get_curr_date + use physconst, only: rair +#endif #endif #ifdef DIRIND @@ -849,18 +851,18 @@ subroutine radiation_tend( & 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) -!#ifdef SPAERO - ! cloud radiative parameters are "in cloud" not "in cell" with SP aerosols -! real(r8) :: sp_liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth -! real(r8) :: sp_liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau -! real(r8) :: sp_liq_tau_w_g(nswbands,pcols,pver) ! liquid assymetry parameter * tau * w -! real(r8) :: sp_liq_tau_w_f(nswbands,pcols,pver) ! liquid forward scattered fraction * tau * w -! real(r8) :: sp_cld_tau (nswbands,pcols,pver) ! liquid extinction optical depth -! real(r8) :: sp_cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau -! real(r8) :: sp_cld_tau_w_g(nswbands,pcols,pver) ! cloud assymetry parameter * w * tau -! real(r8) :: sp_cld_tau_w_f(nswbands,pcols,pver) ! cloud forward scattered fraction * w * tau -! real(r8) :: sp_cld_lw_abs (nlwbands,pcols,pver) ! cloud absorption optics depth (LW) -!#endif +#ifdef SPAERO + ! cloud radiative parameters with SP aerosols are "in cloud" not "in cell" + real(r8) :: sp_liq_tau (nswbands,pcols,pver) ! liquid extinction optical depth + real(r8) :: sp_liq_tau_w (nswbands,pcols,pver) ! liquid single scattering albedo * tau + real(r8) :: sp_liq_tau_w_g(nswbands,pcols,pver) ! liquid assymetry parameter * tau * w + real(r8) :: sp_liq_tau_w_f(nswbands,pcols,pver) ! liquid forward scattered fraction * tau * w + real(r8) :: sp_cld_tau (nswbands,pcols,pver) ! liquid extinction optical depth + real(r8) :: sp_cld_tau_w (nswbands,pcols,pver) ! cloud single scattering albedo * tau + real(r8) :: sp_cld_tau_w_g(nswbands,pcols,pver) ! cloud assymetry parameter * w * tau + real(r8) :: sp_cld_tau_w_f(nswbands,pcols,pver) ! cloud forward scattered fraction * w * tau + real(r8) :: sp_cld_lw_abs (nlwbands,pcols,pver) ! cloud absorption optics depth (LW) +#endif ! cloud radiative parameters are "in cloud" not "in cell" real(r8) :: snow_tau (nswbands,pcols,pver) ! snow extinction optical depth @@ -876,12 +878,12 @@ subroutine radiation_tend( & 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) -!#ifdef SPAERO ! and for SP aerosols (only for SW) -! real(r8) :: sp_c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth -! real(r8) :: sp_c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau -! real(r8) :: sp_c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud assymetry parameter * w * tau -! real(r8) :: sp_c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau -!#endif +#ifdef SPAERO ! and for SP aerosols (only for SW) + real(r8) :: sp_c_cld_tau (nswbands,pcols,pver) ! combined cloud extinction optical depth + real(r8) :: sp_c_cld_tau_w (nswbands,pcols,pver) ! combined cloud single scattering albedo * tau + real(r8) :: sp_c_cld_tau_w_g(nswbands,pcols,pver) ! combined cloud assymetry parameter * w * tau + real(r8) :: sp_c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau +#endif real(r8) :: sfac(1:nswbands) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band @@ -897,34 +899,35 @@ subroutine radiation_tend( & #ifdef DIRIND -!#ifdef SPAERO +#ifdef SPAERO ! Aerosol optical properties, simple plume aerosols -! real(r8) :: sp_tau (pcols,pver,nswbands) ! aerosol extinction optical depth, simple plumes -! real(r8) :: sp_ssa (pcols,pver,nswbands) ! aerosol single scattering albedo, simple plumes -! real(r8) :: sp_asy (pcols,pver,nswbands) ! aerosol assymetry parameter, simple plumes -! + real(r8) :: sp_tau (pcols,pver,nswbands) ! aerosol extinction optical depth, simple plumes + real(r8) :: sp_ssa (pcols,pver,nswbands) ! aerosol single scattering albedo, simple plumes + real(r8) :: sp_asy (pcols,pver,nswbands) ! aerosol assymetry parameter, simple plumes + ! Aerosol optical properties, sum of NorESM + simple plume aerosols -! real(r8) :: sp_per_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth -! real(r8) :: sp_per_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau -! real(r8) :: sp_per_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau -! real(r8) :: sp_per_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau -! -! real(r8), intent(out) :: xcdnc(pcols) ! CDNC modification factor -! real(r8) :: re_mult(pcols,pver) ! Multiplication factor of liquid cloud effective radius, simple plumes -! real(r8) :: re_multeq1(pcols,pver) ! Dummy multiplication factor (=1.0) of liquid cloud effective radius, simple plumes -! real(r8) :: xcdnceq1(pcols) ! Dummy xcdnc -! -!!ak real(r8) :: sp_relca(pcols,pver) ! Modified liquid cloud effective radius, simple plumes -! -! real(r8) :: year_fr ! Fractional year (1903.0 is the 0Z on the first of January 1903, Gregorian) -! -! integer :: yr, mon, day, tod ! date components -!#endif -!#ifdef RFMIPIRF -! real(r8) :: per_aod (pcols,pver,nswbands) ! aerosol single scattering albedo -! real(r8) :: per_ssa (pcols,pver,nswbands) ! aerosol single scattering albedo -! real(r8) :: per_asy (pcols,pver,nswbands) ! aerosol assymetry parameter -!#endif + real(r8) :: sp_per_tau (pcols,0:pver,nswbands) ! aerosol extinction optical depth + real(r8) :: sp_per_tau_w (pcols,0:pver,nswbands) ! aerosol single scattering albedo * tau + real(r8) :: sp_per_tau_w_g(pcols,0:pver,nswbands) ! aerosol assymetry parameter * w * tau + real(r8) :: sp_per_tau_w_f(pcols,0:pver,nswbands) ! aerosol forward scattered fraction * w * tau + + real(r8), intent(out) :: xcdnc(pcols) ! CDNC modification factor + real(r8) :: re_mult(pcols,pver) ! Multiplication factor of liquid cloud effective radius, simple plumes + real(r8) :: re_multeq1(pcols,pver) ! Dummy multiplication factor (=1.0) of liquid cloud effective radius, simple plumes + real(r8) :: xcdnceq1(pcols) ! Dummy xcdnc + +!ak real(r8) :: sp_relca(pcols,pver) ! Modified liquid cloud effective radius, simple plumes + + real(r8) :: year_fr ! Fractional year (1903.0 is the 0Z on the first of January 1903, Gregorian) + + integer :: yr, mon, day, tod ! date components +#endif + +#ifdef RFMIPIRF + real(r8) :: per_aod (pcols,pver,nswbands) ! aerosol single scattering albedo + real(r8) :: per_ssa (pcols,pver,nswbands) ! aerosol single scattering albedo + real(r8) :: per_asy (pcols,pver,nswbands) ! aerosol assymetry parameter +#endif ! 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) @@ -936,9 +939,9 @@ subroutine radiation_tend( & #ifdef AEROCOM real(r8) dod440(pcols),dod550(pcols),dod870(pcols),abs550(pcols),abs550alt(pcols) real(r8) clearod440(pcols),clearod550(pcols),clearod870(pcols),clearabs550(pcols),clearabs550alt(pcols) -!#ifdef RFMIPIRF -! character(len=2) :: c2 -!#endif +#ifdef RFMIPIRF + character(len=2) :: c2 +#endif #endif ! AEROCOM real(r8) ftem_1d(pcols) ! work-array to avoid NAN and pcols/ncol confusion real(r8) Nnatk(pcols,pver,0:nmodes) ! Modal aerosol number concentration @@ -958,9 +961,9 @@ subroutine radiation_tend( & 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 -!#ifdef SPAERO -! real(r8) deltah_km(pcols,pver) ! Layer thickness, unit km -!#endif +#ifdef SPAERO + real(r8) deltah_km(pcols,pver) ! Layer thickness, unit km +#endif #endif real(r8) :: fns(pcols,pverp) ! net shortwave flux @@ -1053,7 +1056,7 @@ subroutine radiation_tend( & #ifdef DIRIND qdirind(:ncol,:,:) = state%q(:ncol,:,:) - if (has_prescribed_volcaero) then + if (has_prescribed_volcaero) then ! old (CMIP5) method call oslo_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,:) @@ -1085,47 +1088,63 @@ subroutine radiation_tend( & if (dosw) then -!#ifdef SPAERO +#ifdef SPAERO !*********************************** SPAERO + ********************************************* ! Define anthrop. aerosol optical properties and "cdnc" for the "simple plumes" climatology -! CALL get_curr_date(yr, mon, day, tod) + CALL get_curr_date(yr, mon, day, tod) -! Petri used a hard-coded year for BACCHUS: either 1850, 1975 or 2005 +!ak yr=2005 + yr=2014 -! yr=2005 + year_fr = yr + (calday-1.0_r8) / 365.0_r8 -! year_fr = yr + (calday-1.0_r8) / 365.0_r8 +!test+ +! write(*,*) 'test yr = ', yr +! write(*,*) 'test year_fr = ', year_fr +!test- !ak+ Need deltah_km before pmxsub is called, due to cloud optics (-> input to pmxsub later) -! do k=1,pver -!! NB have to multiply with 10 to get the same values as in pmxsub, due to different p units! -! rhoda(1:ncol,k) = state%pmid(1:ncol,k)/(rair*state%t(1:ncol,k)) ! unit kg/m^3 -! deltah_km(1:ncol,k)=10._r8*1.e-4_r8*(state%pint(1:ncol,k+1)-state%pint(1:ncol,k))/(rhoda(1:ncol,k)*9.8_r8) -! end do + do k=1,pver +! NB have to multiply with 10 to get the same values as in pmxsub, due to different p units! + rhoda(1:ncol,k) = state%pmid(1:ncol,k)/(rair*state%t(1:ncol,k)) ! unit kg/m^3 + deltah_km(1:ncol,k)=10._r8*1.e-4_r8*(state%pint(1:ncol,k+1)-state%pint(1:ncol,k))/(rhoda(1:ncol,k)*9.8_r8) + end do ! initialization -! re_mult(1:ncol,1:pver) = 1._r8 -! xcdnc(1:ncol) = 1._r8 + re_mult(1:ncol,1:pver) = 1._r8 + xcdnc(1:ncol) = 1._r8 ! for use in calls without the effect of SP aerosols -! re_multeq1(1:ncol,1:pver) = 1._r8 -! xcdnceq1(1:ncol) = 1._r8 + re_multeq1(1:ncol,1:pver) = 1._r8 + xcdnceq1(1:ncol) = 1._r8 !ak- -! CALL simple_plumes_interface(lchnk, ncol, nswbands,state%phis, & -! deltah_km, clon, clat, year_fr, & -! sp_tau, sp_ssa, sp_asy, re_mult, xcdnc) + CALL simple_plumes_interface(lchnk, ncol, nswbands,state%phis, & + deltah_km, clon, clat, year_fr, & + sp_tau, sp_ssa, sp_asy, re_mult, xcdnc) ! When using year 1850 for the MACv2-SP aerosols, switch them off entirely -! IF (yr==1850) THEN -! sp_tau(1:ncol,1:pver,1:nswbands)=0._r8 -! sp_ssa(1:ncol,1:pver,1:nswbands)=0._r8 -! sp_asy(1:ncol,1:pver,1:nswbands)=0._r8 -! re_mult(1:ncol,1:pver) = 1._r8 -! xcdnc(1:ncol) = 1._r8 -! END IF + IF (yr==1850) THEN + sp_tau(1:ncol,1:pver,1:nswbands)=0._r8 + sp_ssa(1:ncol,1:pver,1:nswbands)=0._r8 + sp_asy(1:ncol,1:pver,1:nswbands)=0._r8 + re_mult(1:ncol,1:pver) = 1._r8 + xcdnc(1:ncol) = 1._r8 + END IF + +!!!!!!!!!!!!!!!!!!!!!!!!!! TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! +! To verify/test that we get the same results as without the added +! SPAERO when the SP contribution is zero, activate the following +! (testet to be OK September 2019): +! sp_tau(1:ncol,1:pver,1:nswbands)=0._r8 +! sp_ssa(1:ncol,1:pver,1:nswbands)=0._r8 +! sp_asy(1:ncol,1:pver,1:nswbands)=0._r8 +! re_mult(1:ncol,1:pver) = 1._r8 +! xcdnc(1:ncol) = 1._r8 +!!!!!!!!!!!!!!!!!!!!!!!!!! TEST !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + !*********************************** SPAERO - ********************************************* -!#endif +#endif 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.) @@ -1145,13 +1164,13 @@ subroutine radiation_tend( & call slingo_liq_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f, oldliqwp=.true.) case ('gammadist') -!#ifdef SPAERO +#ifdef SPAERO ! The order of the two calls below has been tested not to make any difference -! call get_liquid_optics_sw(state, pbuf, xcdnceq1, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) -! call get_liquid_optics_sw(state, pbuf, xcdnc, sp_liq_tau, sp_liq_tau_w, sp_liq_tau_w_g, sp_liq_tau_w_f) -!#else + call get_liquid_optics_sw(state, pbuf, xcdnceq1, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) + call get_liquid_optics_sw(state, pbuf, xcdnc, sp_liq_tau, sp_liq_tau_w, sp_liq_tau_w_g, sp_liq_tau_w_f) +#else call get_liquid_optics_sw(state, pbuf, liq_tau, liq_tau_w, liq_tau_w_g, liq_tau_w_f) -!#endif +#endif case default call endrun('liqcldoptics must be either slingo or gammadist') end select @@ -1161,12 +1180,12 @@ subroutine radiation_tend( & 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,:) -!#ifdef SPAERO -! sp_cld_tau(:,:ncol,:) = sp_liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:) -! sp_cld_tau_w(:,:ncol,:) = sp_liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:) -! sp_cld_tau_w_g(:,:ncol,:) = sp_liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:) -! sp_cld_tau_w_f(:,:ncol,:) = sp_liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:) -!#endif +#ifdef SPAERO + sp_cld_tau(:,:ncol,:) = sp_liq_tau(:,:ncol,:) + ice_tau(:,:ncol,:) + sp_cld_tau_w(:,:ncol,:) = sp_liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:) + sp_cld_tau_w_g(:,:ncol,:) = sp_liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:) + sp_cld_tau_w_f(:,:ncol,:) = sp_liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:) +#endif if (cldfsnow_idx > 0) then ! add in snow @@ -1195,38 +1214,38 @@ subroutine radiation_tend( & end if end do end do -!#ifdef SPAERO -! do i = 1, ncol -! do k = 1, pver -! if (cldfprime(i,k) > 0._r8) then -! sp_c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) & -! + cld(i,k)*sp_cld_tau(:,i,k) )/cldfprime(i,k) -! sp_c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) & -! + cld(i,k)*sp_cld_tau_w(:,i,k) )/cldfprime(i,k) -! sp_c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) & -! + cld(i,k)*sp_cld_tau_w_g(:,i,k) )/cldfprime(i,k) -! sp_c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) & -! + cld(i,k)*sp_cld_tau_w_f(:,i,k) )/cldfprime(i,k) -! else -! sp_c_cld_tau(:,i,k) = 0._r8 -! sp_c_cld_tau_w(:,i,k) = 0._r8 -! sp_c_cld_tau_w_g(:,i,k) = 0._r8 -! sp_c_cld_tau_w_f(:,i,k) = 0._r8 -! end if -! end do -! end do -!#endif +#ifdef SPAERO + do i = 1, ncol + do k = 1, pver + if (cldfprime(i,k) > 0._r8) then + sp_c_cld_tau(:,i,k) = ( cldfsnow(i,k)*snow_tau(:,i,k) & + + cld(i,k)*sp_cld_tau(:,i,k) )/cldfprime(i,k) + sp_c_cld_tau_w(:,i,k) = ( cldfsnow(i,k)*snow_tau_w(:,i,k) & + + cld(i,k)*sp_cld_tau_w(:,i,k) )/cldfprime(i,k) + sp_c_cld_tau_w_g(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_g(:,i,k) & + + cld(i,k)*sp_cld_tau_w_g(:,i,k) )/cldfprime(i,k) + sp_c_cld_tau_w_f(:,i,k) = ( cldfsnow(i,k)*snow_tau_w_f(:,i,k) & + + cld(i,k)*sp_cld_tau_w_f(:,i,k) )/cldfprime(i,k) + else + sp_c_cld_tau(:,i,k) = 0._r8 + sp_c_cld_tau_w(:,i,k) = 0._r8 + sp_c_cld_tau_w_g(:,i,k) = 0._r8 + sp_c_cld_tau_w_f(:,i,k) = 0._r8 + end if + end do + end do +#endif 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,:) -!#ifdef SPAERO -! sp_c_cld_tau(:,:ncol,:) = sp_cld_tau(:,:ncol,:) -! sp_c_cld_tau_w(:,:ncol,:) = sp_cld_tau_w(:,:ncol,:) -! sp_c_cld_tau_w_g(:,:ncol,:) = sp_cld_tau_w_g(:,:ncol,:) -! sp_c_cld_tau_w_f(:,:ncol,:) = sp_cld_tau_w_f(:,:ncol,:) -!#endif +#ifdef SPAERO + sp_c_cld_tau(:,:ncol,:) = sp_cld_tau(:,:ncol,:) + sp_c_cld_tau_w(:,:ncol,:) = sp_cld_tau_w(:,:ncol,:) + sp_c_cld_tau_w_g(:,:ncol,:) = sp_cld_tau_w_g(:,:ncol,:) + sp_c_cld_tau_w_f(:,:ncol,:) = sp_cld_tau_w_f(:,:ncol,:) +#endif end if ! Output cloud optical depth fields for the visible band @@ -1317,8 +1336,6 @@ subroutine radiation_tend( & ! qdirind(:ncol,:,l_so4_a1) = 0.0_r8 ! qdirind(:ncol,:,l_so4_na) = 0.0_r8 !TEST -!cak+ Calculate CAM5-Oslo/NorESM2 aerosol optical parameters -! (move to aer_rad_props.F90? No, then it cannot be called for night-time calculations...) ! ! Volcanic optics for solar (SW) bands do band=1, solar_bands @@ -1370,54 +1387,53 @@ subroutine radiation_tend( & aodvis, absvis) #endif -!#ifdef RFMIPIRF -!! Extra RFMIP-IRF diagnostics for each SW wave-length/number band -! per_aod(:,:,:)=0._r8 -! per_ssa(:,:,:)=0._r8 -! per_asy(:,:,:)=0._r8 -! DO i=1,ncol -!! DO k=0,pver -! DO k=1,pver -! DO ns=1,nswbands -! per_aod(i,k,ns)=per_tau(i,k,ns) -! per_ssa(i,k,ns)=min(per_tau_w(i,k,ns)/(per_tau(i,k,ns)+eps),1._r8) -! per_asy(i,k,ns)=min(per_tau_w_g(i,k,ns)/(per_tau_w(i,k,ns)+eps),1._r8) -! ENDDO -! ENDDO -! ENDDO -! do ns=1,nswbands -! write(c2,'(I2)') ns -! call outfld('AERTAUBND'//trim(adjustl(c2)),per_aod(:,:,ns),pcols,lchnk) -! call outfld('AERSSABND'//trim(adjustl(c2)),per_ssa(:,:,ns),pcols,lchnk) -! call outfld('AERASYBND'//trim(adjustl(c2)),per_asy(:,:,ns),pcols,lchnk) -! enddo -!#endif +#ifdef RFMIPIRF +! Extra RFMIP-IRF diagnostics for each SW wave-length/number band + per_aod(:,:,:)=0._r8 + per_ssa(:,:,:)=0._r8 + per_asy(:,:,:)=0._r8 + DO i=1,ncol + DO k=1,pver + DO ns=1,nswbands + per_aod(i,k,ns)=per_tau(i,k,ns) + per_ssa(i,k,ns)=min(per_tau_w(i,k,ns)/(per_tau(i,k,ns)+eps),1._r8) + per_asy(i,k,ns)=min(per_tau_w_g(i,k,ns)/(per_tau_w(i,k,ns)+eps),1._r8) + ENDDO + ENDDO + ENDDO + do ns=1,nswbands + write(c2,'(I2)') ns + call outfld('AERTAUBND'//trim(adjustl(c2)),per_aod(:,:,ns),pcols,lchnk) + call outfld('AERSSABND'//trim(adjustl(c2)),per_ssa(:,:,ns),pcols,lchnk) + call outfld('AERASYBND'//trim(adjustl(c2)),per_asy(:,:,ns),pcols,lchnk) + enddo +#endif -!#ifdef SPAERO +#ifdef SPAERO !*********************************** SPAERO + ********************************************* ! Use the anthropogenic aerosol optical properties for the "simple plumes" climatology ! Add "simple plumes" to NorESM2 aerosols (which should be defined using year 1850 emissions) ! Set aerosol optical properties zero for the top layer (0) -! sp_per_tau(1:ncol,0,1:nswbands) = 0._r8 -! sp_per_tau_w(1:ncol,0,1:nswbands) = 0._r8 -! sp_per_tau_w_g(1:ncol,0,1:nswbands) = 0._r8 -! sp_per_tau_w_f(1:ncol,0,1:nswbands) = 0._r8 + sp_per_tau(1:ncol,0,1:nswbands) = 0._r8 + sp_per_tau_w(1:ncol,0,1:nswbands) = 0._r8 + sp_per_tau_w_g(1:ncol,0,1:nswbands) = 0._r8 + sp_per_tau_w_f(1:ncol,0,1:nswbands) = 0._r8 ! Other layers -! DO i=1,ncol -! DO k=1,pver -! DO ns=1,nswbands -! sp_per_tau(i,k,ns)=per_tau(i,k,ns) + sp_tau(i,k,ns) -! sp_per_tau_w(i,k,ns)=per_tau_w(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns) -! sp_per_tau_w_g(i,k,ns)=per_tau_w_g(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns)*sp_asy(i,k,ns) -! sp_per_tau_w_f(i,k,ns)=per_tau_w_f(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns)*sp_asy(i,k,ns)**2 -! ENDDO -! ENDDO -! ENDDO + DO i=1,ncol + DO k=1,pver + DO ns=1,nswbands + sp_per_tau(i,k,ns)=per_tau(i,k,ns) + sp_tau(i,k,ns) + sp_per_tau_w(i,k,ns)=per_tau_w(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns) + sp_per_tau_w_g(i,k,ns)=per_tau_w_g(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns)*sp_asy(i,k,ns) + sp_per_tau_w_f(i,k,ns)=per_tau_w_f(i,k,ns) + sp_tau(i,k,ns)*sp_ssa(i,k,ns)*sp_asy(i,k,ns)**2 + ENDDO + ENDDO + ENDDO !*********************************** SPAERO - ********************************************* -!#endif +#endif #endif ! DIRIND @@ -1457,12 +1473,7 @@ subroutine radiation_tend( & rd%fsnirt, rd%fsnrtc, rd%fsnirtsq, fsns, rd%fsnsc, & rd%fsdsc, fsds, cam_out%sols, cam_out%soll, cam_out%solsd, & !akc6+ -!#ifdef AEROFFL -! cam_out%solld, fns, fcns, fds, fdsc, Nday, Nnite, & cam_out%solld, fns, fcns, idrf, Nday, Nnite, & -!#else -! cam_out%solld, fns, fcns, Nday, Nnite, & -!#endif !akc6- 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, & @@ -1525,102 +1536,100 @@ subroutine radiation_tend( & E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false.) -!#ifdef SPAERO +#ifdef SPAERO ! ! Dump shortwave radiation information to history tape buffer (diagnostics) ! ! Added, as for BACCHUS by P. Räisänen -! call outfld('FSNT_SP ',fsnt ,pcols,lchnk) -! call outfld('FSNS_SP ',fsns ,pcols,lchnk) -! call outfld('FSNTC_SP',rd%fsntc ,pcols,lchnk) -! call outfld('FSNSC_SP',rd%fsnsc ,pcols,lchnk) -!#endif + call outfld('FSNT_SP ',fsnt ,pcols,lchnk) + call outfld('FSNS_SP ',fsns ,pcols,lchnk) + call outfld('FSNTC_SP',rd%fsntc ,pcols,lchnk) + call outfld('FSNSC_SP',rd%fsnsc ,pcols,lchnk) +#endif -!#ifdef SPAERO +#ifdef SPAERO !*********************************** SPAERO + ********************************************* ! THIRD CALL INCLUDING SIMPLE PLUME AEROSOLS FOR ONLY THE DIRECT EFFECT (SW ERF ari) -! call rad_rrtmg_sw( & -! lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & -! cldfprime, & + call rad_rrtmg_sw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + cldfprime, & !ak+ per_tau, per_tau_w, per_tau_w_g, per_tau_w_f, & -! sp_per_tau, sp_per_tau_w, sp_per_tau_w_g, sp_per_tau_w_f, & + sp_per_tau, sp_per_tau_w, sp_per_tau_w_g, sp_per_tau_w_f, & !ak- -! 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, & + 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, & !akc6+ -!#ifdef AEROFFL -! cam_out%solld, fns, fcns, idrf, Nday, Nnite, & -!#else -! cam_out%solld, fns, fcns, Nday, Nnite, & -!#endif +#ifdef AEROFFL + cam_out%solld, fns, fcns, idrf, Nday, Nnite, & +#else + cam_out%solld, fns, fcns, Nday, Nnite, & +#endif !akc6- -! 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.) + 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.) ! ! Dump shortwave radiation information to history tape buffer (diagnostics) ! -! Added, as for BACCHUS by P. Räisänen -! call outfld('FSNT_SP2',fsnt ,pcols,lchnk) -! call outfld('FSNS_SP2',fsns ,pcols,lchnk) -! call outfld('FSNTCSP2',rd%fsntc ,pcols,lchnk) -! call outfld('FSNSCSP2',rd%fsnsc ,pcols,lchnk) + call outfld('FSNT_SP2',fsnt ,pcols,lchnk) + call outfld('FSNS_SP2',fsns ,pcols,lchnk) + call outfld('FSNTCSP2',rd%fsntc ,pcols,lchnk) + call outfld('FSNSCSP2',rd%fsnsc ,pcols,lchnk) ! FOURTH CALL INCLUDING SIMPLE PLUME AEROSOLS FOR BOTH THE DIRECT AND THE 1. INDIRECT EFFECT -! call rad_rrtmg_sw( & -! lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & -! cldfprime, & + call rad_rrtmg_sw( & + lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & + cldfprime, & !ak+ per_tau, per_tau_w, per_tau_w_g, per_tau_w_f, & -! sp_per_tau, sp_per_tau_w, sp_per_tau_w_g, sp_per_tau_w_f, & + sp_per_tau, sp_per_tau_w, sp_per_tau_w_g, sp_per_tau_w_f, & !ak- -! 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, & + 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, & !akc6+ -!#ifdef AEROFFL -! cam_out%solld, fns, fcns, idrf, Nday, Nnite, & -!#else -! cam_out%solld, fns, fcns, Nday, Nnite, & -!#endif +#ifdef AEROFFL + cam_out%solld, fns, fcns, idrf, Nday, Nnite, & +#else + cam_out%solld, fns, fcns, Nday, Nnite, & +#endif !akc6- !ak+ IdxDay, IdxNite, su, sd, E_cld_tau=c_cld_tau, & !ak+ E_cld_tau_w=c_cld_tau_w, E_cld_tau_w_g=c_cld_tau_w_g, & !ak+ E_cld_tau_w_f=c_cld_tau_w_f, old_convert=.false.) -! IdxDay, IdxNite, su, sd, E_cld_tau=sp_c_cld_tau, & -! E_cld_tau_w=sp_c_cld_tau_w, E_cld_tau_w_g=sp_c_cld_tau_w_g,& -! E_cld_tau_w_f=sp_c_cld_tau_w_f, old_convert=.false.) + IdxDay, IdxNite, su, sd, E_cld_tau=sp_c_cld_tau, & + E_cld_tau_w=sp_c_cld_tau_w, E_cld_tau_w_g=sp_c_cld_tau_w_g,& + E_cld_tau_w_f=sp_c_cld_tau_w_f, old_convert=.false.) !ak- ! ! Dump shortwave radiation information to history tape buffer (diagnostics) ! -! Added, as for BACCHUS by P. Räisänen -! call outfld('FSNT_SP3',fsnt ,pcols,lchnk) -! call outfld('FSNS_SP3',fsns ,pcols,lchnk) -! call outfld('FSNTCSP3',rd%fsntc ,pcols,lchnk) -! call outfld('FSNSCSP3',rd%fsnsc ,pcols,lchnk) + call outfld('FSNT_SP3',fsnt ,pcols,lchnk) + call outfld('FSNS_SP3',fsns ,pcols,lchnk) + call outfld('FSNTCSP3',rd%fsntc ,pcols,lchnk) + call outfld('FSNSCSP3',rd%fsnsc ,pcols,lchnk) !*********************************** SPAERO + ********************************************* -!#endif +#endif -!#ifdef RFMIPIRF +#ifdef RFMIPIRF ! Extra RFMIP-IRF diagnostics for each SW wave-length/number band -! do ns=1,nswbands -! write(c2,'(I2)') ns -! call outfld('SDBND'//trim(adjustl(c2)),sd(:,:,ns),pcols,lchnk) -! call outfld('SUBND'//trim(adjustl(c2)),su(:,:,ns),pcols,lchnk) -! enddo -!#endif + do ns=1,nswbands + write(c2,'(I2)') ns + call outfld('SDBND'//trim(adjustl(c2)),sd(:,:,ns),pcols,lchnk) + call outfld('SUBND'//trim(adjustl(c2)),su(:,:,ns),pcols,lchnk) + enddo +#endif !ak+ Has been moved from above to after the last rad_rrtmg_sw call... ! Output net fluxes at 200 mb @@ -1724,14 +1733,12 @@ subroutine radiation_tend( & lu, ld) #ifdef DIRIND -#ifdef AEROFFL ! FLNT_ORG is just for temporary testing vs. FLNT -!#ifdef AEROCOM -! call outfld('FLNT_ORG',flnt(:) ,pcols,lchnk) +#ifdef AEROFFL ftem_1d(1:ncol) = cam_out%flwds(1:ncol) - flns(1:ncol) call outfld('FLUS ',ftem_1d ,pcols,lchnk) -!#endif ! AEROCOM #endif ! AEROFFL -!#ifdef AEROCOM +!#ifdef AEROCOM +! only used for testing of internal consistency ! do i=1,ncol ! do k=1,pver ! aerlwabs01(i,k) = aer_lw_abs(i,k,16) @@ -1740,14 +1747,14 @@ subroutine radiation_tend( & ! call outfld('AERLWA01',aerlwabs01,pcols,lchnk) !#endif -!#ifdef RFMIPIRF +#ifdef RFMIPIRF ! Extra RFMIP-IRF diagnostics for each LW wave-length/number band -! do ns=1,nlwbands -! write(c2,'(I2)') ns -! call outfld('LDBND'//trim(adjustl(c2)),ld(:,:,ns),pcols,lchnk) -! call outfld('LUBND'//trim(adjustl(c2)),lu(:,:,ns),pcols,lchnk) -! enddo -!#endif + do ns=1,nlwbands + write(c2,'(I2)') ns + call outfld('LDBND'//trim(adjustl(c2)),ld(:,:,ns),pcols,lchnk) + call outfld('LUBND'//trim(adjustl(c2)),lu(:,:,ns),pcols,lchnk) + enddo +#endif #endif ! DIRIND @@ -1827,8 +1834,8 @@ subroutine radiation_tend( & end if ! if (dosw .or. dolw) then - ! output rad inputs and resulting heating rates - call rad_data_write( pbuf, state, cam_in, coszrs ) + ! output rad inputs and resulting heating rates + call rad_data_write( pbuf, state, cam_in, coszrs ) ! Compute net radiative heating tendency call radheat_tend(state, pbuf, ptend, qrl, qrs, fsns, & diff --git a/src/physics/cam_oslo/radlw.F90 b/src/physics/cam_oslo/radlw.F90 index df8dd0c4b5..8ad7fed365 100644 --- a/src/physics/cam_oslo/radlw.F90 +++ b/src/physics/cam_oslo/radlw.F90 @@ -193,10 +193,10 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & if (associated(lu)) lu(1:ncol,:,:) = 0.0_r8 if (associated(ld)) ld(1:ncol,:,:) = 0.0_r8 -!#ifdef RFMIPIRF -! lu(1:ncol,:,:) = 0.0_r8 -! ld(1:ncol,:,:) = 0.0_r8 -!#endif +#ifdef RFMIPIRF + lu(1:ncol,:,:) = 0.0_r8 + ld(1:ncol,:,:) = 0.0_r8 +#endif call rrtmg_lw(lchnk ,ncol ,rrtmg_levs ,icld , & r_state%pmidmb ,r_state%pintmb ,r_state%tlay ,r_state%tlev ,tsfc ,r_state%h2ovmr, & @@ -266,23 +266,23 @@ subroutine rad_rrtmg_lw(lchnk ,ncol ,rrtmg_levs,r_state, & ! Pass spectral fluxes, reverse layering ! order=(/3,1,2/) maps the first index of lwuflxs to the third index of lu. -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(lu)) then -!#endif +#endif lu(:ncol,pverp-rrtmg_levs+1:pverp,:) = reshape(lwuflxs(:,:ncol,rrtmg_levs:1:-1), & (/ncol,rrtmg_levs,nbndlw/), order=(/3,1,2/)) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(ld)) then -!#endif +#endif ld(:ncol,pverp-rrtmg_levs+1:pverp,:) = reshape(lwdflxs(:,:ncol,rrtmg_levs:1:-1), & (/ncol,rrtmg_levs,nbndlw/), order=(/3,1,2/)) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif call t_stopf('rrtmg_lw') diff --git a/src/physics/cam_oslo/radsw.F90 b/src/physics/cam_oslo/radsw.F90 index 24a3b865fd..8e677b9246 100644 --- a/src/physics/cam_oslo/radsw.F90 +++ b/src/physics/cam_oslo/radsw.F90 @@ -276,11 +276,6 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: fdsc(pcols,pverp) ! Downward clear-sky flux (added for CRM) #ifdef AEROFFL -! real(r8), intent(out) :: fds(pcols,pverp) ! Downward flux (added for CRM) -! real(r8), intent(out) :: fdsc(pcols,pverp) ! Downward clear-sky flux (added for CRM) -!#else -! real(r8) :: fds(pcols,pverp) ! Downward flux (added for CRM) -! real(r8) :: fdsc(pcols,pverp) ! Downward clear-sky flux (added for CRM) logical, intent(in) :: idrf #endif @@ -337,10 +332,10 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & if (associated(su)) su(1:ncol,:,:) = 0.0_r8 if (associated(sd)) sd(1:ncol,:,:) = 0.0_r8 -!#ifdef RFMIPIRF -! su(1:ncol,:,:) = 0.0_r8 -! sd(1:ncol,:,:) = 0.0_r8 -!#endif +#ifdef RFMIPIRF + su(1:ncol,:,:) = 0.0_r8 + sd(1:ncol,:,:) = 0.0_r8 +#endif ! If night everywhere, return: if ( Nday == 0 ) then @@ -603,23 +598,23 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Set spectral fluxes, reverse layering ! order=(/3,1,2/) maps the first index of swuflxs to the third index of su. -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(su)) then -!#endif +#endif su(1:Nday,pverp-rrtmg_levs+1:pverp,:) = reshape(swuflxs(:,1:Nday,rrtmg_levs:1:-1), & (/Nday,rrtmg_levs,nbndsw/), order=(/3,1,2/)) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(sd)) then -!#endif +#endif sd(1:Nday,pverp-rrtmg_levs+1:pverp,:) = reshape(swdflxs(:,1:Nday,rrtmg_levs:1:-1), & (/Nday,rrtmg_levs,nbndsw/), order=(/3,1,2/)) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif call t_stopf('rrtmg_sw') @@ -649,21 +644,21 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & call ExpDayNite(fsnrtoac, Nday, IdxDay, Nnite, IdxNite, 1, pcols) call ExpDayNite(fsnrtoaq, Nday, IdxDay, Nnite, IdxNite, 1, pcols) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(su)) then -!#endif +#endif call ExpDayNite(su, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp, 1, nbndsw) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif -!#ifndef RFMIPIRF +#ifndef RFMIPIRF if (associated(sd)) then -!#endif +#endif call ExpDayNite(sd, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp, 1, nbndsw) -!#ifndef RFMIPIRF +#ifndef RFMIPIRF end if -!#endif +#endif ! these outfld calls don't work for spmd only outfield in scm mode (nonspmd) #ifndef OSLO_AERO @@ -684,8 +679,6 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & #ifdef AEROFFL if(idrf) then -! call ExpDayNite(fusc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) -! call ExpDayNite(fdsc,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) call outfld('FUSCDRF ', fusc, pcols, lchnk) call outfld('FDSCDRF ', fdsc, pcols, lchnk) endif diff --git a/src/physics/cam_oslo/simple_plumes_interface.F90 b/src/physics/cam_oslo/simple_plumes_interface.F90 new file mode 100644 index 0000000000..7c1c833323 --- /dev/null +++ b/src/physics/cam_oslo/simple_plumes_interface.F90 @@ -0,0 +1,242 @@ + +SUBROUTINE simple_plumes_interface (lchnk,ncol, nswbands,surf_geopot, & +!ak+ deltah_km, cdnc, clon, clat, year_fr, & + deltah_km, clon, clat, year_fr, & + sp_tau, sp_ssa, sp_asy, re_mult, xcdnc) +!ak- + +!--------------------------------------------------------------------------- +! Purpose: +! - interface between the MAC-SP ("simple plumes") aerosol climatology +! and CAM/NorESM radiation scheme +! * called from "radiation_tend" in radiation.F90 +! * calls "sp_aop_profile" in mo_simple_plumes.F90 +! +! May 2016, P. Räisänen : initial version +! Sept 2019, A. Kirkevåg : modified for CAM6-Oslo / NorESM2 (+ extra output) +!--------------------------------------------------------------------------- + +#include + + USE ppgrid, only: pcols, pver + USE shr_kind_mod, only: r8=>shr_kind_r8 + USE physconst, ONLY: pi, rga +#ifdef SPAERO + USE radconstants, ONLY: wav_sp +#endif + USE CAM_HISTORY, ONLY: outfld + USE mo_simple_plumes, ONLY: sp_aop_profile +!ak+ + use cam_logfile, only: iulog +!ak- + + IMPLICIT NONE + +!ak (NorESM1-M): INTEGER, PARAMETER :: nextra = 10 ! Number of layers added between the + INTEGER, PARAMETER :: nextra = 13 ! Number of layers added between the + ! sea level and model surface for computing + ! the simple plumes aerosols + +! Arguments + INTEGER, INTENT(in) :: lchnk ! chunk identifier + INTEGER, INTENT(in) :: ncol ! number of columns + INTEGER, INTENT(in) :: nswbands ! number of shortwave bands + + REAL(r8), INTENT(in) :: surf_geopot(pcols) ! surface geopotential + REAL(r8), INTENT(in) :: deltah_km(pcols,pver) ! layer thicknesses (1=uppermost layer) +!ak REAL(r8), INTENT(in) :: cdnc(pcols,pver) ! in-cloud CDNC [cm-3?] (1=uppermost layer) + REAL(r8), INTENT(in) :: clon(pcols) ! longitude (in radians) + REAL(r8), INTENT(in) :: clat(pcols) ! latitude (in radians) + REAL(r8), INTENT(in) :: year_fr ! Fractional year (1903.0 is the 0Z on + !the first of January 1903, Gregorian) + + REAL(r8), INTENT(out) :: sp_tau (pcols,pver,nswbands) ! aerosol extinction optical depth (1=uppermost layer) + REAL(r8), INTENT(out) :: sp_ssa (pcols,pver,nswbands) ! aerosol single scattering albedo (1=uppermost layer) + REAL(r8), INTENT(out) :: sp_asy (pcols,pver,nswbands) ! aerosol assymetry parameter (1=uppermost layer) + REAL(r8), INTENT(out) :: re_mult(pcols,pver) ! Multiplication factor of liquid cloud effective radius + REAL(r8), INTENT(out) :: xcdnc(pcols) ! CDNC modification factor + +! Local variables: input for SP_AOP_PROFILE + + REAL(r8) :: col_lon(pcols) ! longitudes in degrees + REAL(r8) :: col_lat(pcols) ! latitudes in degrees + REAL(r8) :: oro(pcols) ! surface orography (height above sea level in [m]) + REAL(r8) :: z_sp(pcols,pver+nextra) ! layer mid-point height above sea level [m], (1=lowermost layer) + REAL(r8) :: dz_sp(pcols,pver+nextra) ! layer thickness [m] (1=lowermost layer) + + REAL(r8) :: lambda ! wavelength in [nm] + +! Local variables: output from SP_AOP_PROFILE + + REAL(r8) :: aod_prof(pcols,pver+nextra) ! aerosol optical depth (1=lowermost layer) + REAL(r8) :: ssa_prof(pcols,pver+nextra) ! aerosol single-scattering albedo (1=lowermost layer) + REAL(r8) :: asy_prof(pcols,pver+nextra) ! aerosol asymmetry parameter (1=lowermost layer) + +! Local variables: fields written to model output. + REAL(r8) :: aodvis_sp(pcols) ! AOD at 0.35-0.64 µm + REAL(r8) :: absvis_sp(pcols) ! Absorption AOD at 0.35-0.64 µm +!ak+ + REAL(r8) :: aodv3d_sp(pcols,pver) ! 3D AOD at 0.35-0.64 µm + REAL(r8) :: absv3d_sp(pcols,pver) ! 3D absorption AOD at 0.35-0.64 µm +!ak- +! Other local variables + REAL(r8) :: zhalf(pcols,pver+1) ! layer interface height above sea level [m] (1=uppermost) + REAL(r8) :: zhalf_sp(pcols,pver+nextra+1) ! layer interface height above sea level [m] + ! grid used for simple plume climatology (1=surface) + + REAL(r8), PARAMETER :: rad2deg = 180._r8 / pi + +!ak+ REAL(r8) :: zmin, zmax, deltaz, rv_mult, cdnc_sp, epsilon, epsilon_sp, & +!ak beta, beta_sp, cdnc_test + REAL(r8) :: zmin, zmax, deltaz +!ak- + + INTEGER :: i,j,k,ns + +!**************************************** +! Prepare input data for SP_AOP_PROFILE +!**************************************** + +! Longitude and latitude in degrees + col_lon(:) = rad2deg * clon(:) + col_lat(:) = rad2deg * clat(:) +! Surface elevation + oro(:) = surf_geopot(:) * rga + +!test + do i=1,pcols +! write(*,*) 'i, surf_geopot(i) = ', i, surf_geopot(i) + if(surf_geopot(i).gt.100000._r8) oro(i)=0._r8 + enddo +!test + +! Define the vertical grids +! 1) Model half-levels, layers indices decreasing upwards + + zhalf(:,pver+1) = oro(:) + DO k=pver,1,-1 + zhalf(:,k) = zhalf(:,k+1) + 1000._r8 * deltah_km(:,k) + ENDDO + +! 2) Layer mid-point and half-level elevations for calling the simple plume +! climatology. A separate grid mist be used for this, since the grid must +! start from the sea level (0 m), not from the model orographic height. +! Otherwise, elevated terrain does not result in reduced total AOD, as it +! should. +! +! To accomplish this, NEXTRA equally thick layers are added between the +! sea-level and the model orographic height. The numbering of the simple +! plume layers goes upwards, lowest level = 1. + + zhalf_sp(:,1) = 0._r8 + zhalf_sp(:,nextra+1) = MAX(zhalf(:,pver+1),0.001_r8) + + DO j=2,nextra + zhalf_sp(:,j) = REAL(j-1)/REAL(nextra) * zhalf_sp(:,nextra+1) + ENDDO + +! Above the model surface, use the model grid for computing the simple +! plumes aerosols. Indexing: +! - lowest model half-level plev+1 = half-level nextra+1 for simple plumes +! - uppermost model half-level 1 = half-level pver+nextra+1 for simple plumes + + DO j=nextra+2,pver+nextra+1 + k = pver+nextra+2-j + zhalf_sp(:,j) = zhalf(:,k) + ENDDO + + DO j=1,pver+nextra + z_sp(:,j) = 0.5_r8*(zhalf_sp(:,j)+zhalf_sp(:,j+1)) + dz_sp(:,j) = zhalf_sp(:,j+1)-zhalf_sp(:,j) + ENDDO + +!*********************************************************************** +! Loop over spectral bands (SW only) + + DO ns=1,nswbands +! Wavelength in nanometres +!ak lambda = 1000._r8*wav_sp(ns) ! wav_sp is now a wave number (cm-1) +#ifdef SPAERO + lambda = 1.e7_r8/wav_sp(ns) +#else ! Hack to make the model compile without SPAERO defined + lambda = 1._r8 +#endif + +! Calculate aerosol profiles from the "simple plumes" climatology + + CALL sp_aop_profile(pver+nextra, ncol, lambda, oro, col_lon, col_lat, & + year_fr, z_sp, dz_sp, & + xcdnc, aod_prof, ssa_prof, asy_prof) + +! Mapping to model vertical grid +! model layer 1 <=> simple plumes layer pver+nextra +! model layer pver <=> simple plumes layer nextra+1 + + DO i=1,ncol + DO k=nextra+1,nextra+pver + j = pver+nextra+1-k + sp_tau(i,j,ns) = aod_prof(i,k) + sp_ssa(i,j,ns) = ssa_prof(i,k) + sp_asy(i,j,ns) = asy_prof(i,k) +! Security settings (are these needed?) + sp_ssa(i,j,ns) = MIN(MAX(sp_ssa(i,j,ns),1.E-6_r8),0.999999_r8) + sp_asy(i,j,ns) = MIN(MAX(sp_asy(i,j,ns),0._r8),0.99_r8) + ENDDO + ENDDO + ENDDO !1,nswbands + +! Conversion from CDNC modification factor to effective radius multiplier +! This assumes that XCDNC = total CDNC / CDNC without simple-plume aerosols + + DO i=1,ncol + xcdnc(i) = MIN(MAX(xcdnc(i),1._r8),100._r8) +! First, multiplication factor for volume-mean radius (r_v is proportional to CDNC**(-1/3)) +!ak+ rv_mult = xcdnc(i)**(-0.333333_r8) +!ak- +! Include dispersion effect following Rotstayn & Liu (GRL 2009, 36, L10801), +! eq. 2 (OLDBETA), similar to "cldwat.f90". +! beta = effective radius / volume-mean radius + + DO j=1,pver +!ak cdnc_sp = cdnc(i,j) * xcdnc(i) +!ak epsilon = 1._r8-0.7_r8*EXP(-0.003*cdnc(i,j)) +!ak epsilon_sp = 1._r8-0.7_r8*EXP(-0.003*cdnc_sp) +!ak beta = (1._r8+2._r8*epsilon**2 )**0.666667_r8 / (1._r8+epsilon**2 )**0.333333_r8 +!ak beta_sp = (1._r8+2._r8*epsilon_sp**2)**0.666667_r8 / (1._r8+epsilon_sp**2)**0.333333_r8 +! +! Multiplication factor for effective radius, including the dispersion effect +! re_mult(i,j) = beta_sp/beta * rv_mult +! For BACCHUS, without the dispersion effect (P. Räisänen, 14 March 2017) +!ak re_mult(i,:) = xcdnc(i)**(-0.333333_r8) + re_mult(i,j) = xcdnc(i)**(-0.333333_r8) + ENDDO + ENDDO + +! Write fields to model output + + aodvis_sp(:) = 0._r8 + absvis_sp(:) = 0._r8 +!ak+ + aodv3d_sp(:,:) = 0._r8 + absv3d_sp(:,:) = 0._r8 +!ak- + + DO j=1,pver +!ak aodvis_sp(:) = aodvis_sp(:)+sp_tau(:,j,8) +!ak absvis_sp(:) = absvis_sp(:)+sp_tau(:,j,8)*(1._r8-sp_ssa(:,j,8)) + aodv3d_sp(:,j) = sp_tau(:,j,10) + absv3d_sp(:,j) = sp_tau(:,j,10)*(1._r8-sp_ssa(:,j,10)) + aodvis_sp(:) = aodvis_sp(:)+aodv3d_sp(:,j) + absvis_sp(:) = absvis_sp(:)+absv3d_sp(:,j) + ENDDO + +!ak+ + call outfld('AODVISSP',aodvis_sp ,pcols,lchnk) + call outfld('ABSVISSP',absvis_sp ,pcols,lchnk) + call outfld('XCDNC_SP',xcdnc ,pcols,lchnk) + call outfld('AODV3DSP',aodv3d_sp ,pcols,lchnk) + call outfld('ABSV3DSP',absv3d_sp ,pcols,lchnk) +!ak- + + RETURN +END SUBROUTINE simple_plumes_interface