pkg/atm_phys/radiation_mod.F90

b2ea1d2979 Jean*0001 module radiation_mod
                0002
                0003 ! ==================================================================================
                0004 ! ==================================================================================
                0005
                0006 !  use fms_mod,               only: open_file, check_nml_error, &
                0007 !                                   mpp_pe, close_file
                0008
                0009    use    gcm_params_mod,     only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit
                0010    use    constants_mod,      only: stefan, cp_air, grav, pstd_mks
                0011
                0012 !  use    diag_manager_mod,   only: register_diag_field, send_data
                0013
                0014 !  use    time_manager_mod,   only: time_type, &
                0015 !                                   operator(+), operator(-), operator(/=)
                0016
                0017 !==================================================================================
                0018 implicit none
                0019 private
                0020 !==================================================================================
                0021
                0022 ! version information
                0023
15388b052e Jean*0024 character(len=128) :: version='$Id: radiation_mod.F90,v 1.7 2017/08/11 20:48:51 jmc Exp $'
b2ea1d2979 Jean*0025 character(len=128) :: tag='homemade'
                0026
                0027 !==================================================================================
                0028
                0029 ! public interfaces
                0030
                0031 public :: radiation_init, radiation_down, radiation_up, radiation_end
                0032 !==================================================================================
                0033
                0034 ! module variables
9f7d507fb9 Jean*0035 ! select_incSW :: select expression for Incoming SW radiation @ the top
                0036 !              :: =0 : no season ; =1 : circular orbit planet (obliquity only)
88391fb671 jm-c 0037 ! ozone_in_SW  :: =1 : account for Ozone in downward SW absorption ; =0: ignore it
                0038 ! two_stream_SW :: =1 : also accont for Ozone in upward SW absorption ; =0: ignore
                0039 !                 Note: requires ozone_in_SW=1 to use two_stream_SW=1
9f7d507fb9 Jean*0040 ! yearLength   :: length of solar year in seconds
                0041 ! yearPhase    :: phase in solar year [0-1] relative to NH winter solstice
                0042 ! obliquity    :: obliquity of Earth rotation axis in degre
b2ea1d2979 Jean*0043 logical :: initialized =.false.
9f7d507fb9 Jean*0044 integer :: select_incSW    = 0
88391fb671 jm-c 0045 ! MK added switches for O3 scheme (default off):
                0046 integer :: ozone_in_SW     = 0
                0047 integer :: two_stream_SW   = 0
b2ea1d2979 Jean*0048
                0049 real    :: solar_constant  = 1360.0
                0050 real    :: del_sol         = 1.4
                0051 ! modif omp: winter/summer hemisphere
                0052 real    :: del_sw          = 0.0
9f7d507fb9 Jean*0053 real    :: yearLength = 86400.*360.
                0054 real    :: yearPhase  = 10./365.    ! winter solstice = 22.Dec.h00
                0055 real    :: obliquity  = 23.45
b2ea1d2979 Jean*0056 real    :: atm_abs         = 0.0
                0057 real    :: sw_diff         = 0.0
                0058 real    :: albedo_value    = 0.06
                0059 real    :: solar_exponent  = 4.0
cf7369fae0 Jean*0060 real    :: sw_co2          = 0.0596
                0061 real    :: wv_exponent     = 4.0
                0062 ! only used with wv_exponent > 0:
                0063 real    :: ir_tau_eq       = 6.0
                0064 real    :: ir_tau_pole     = 1.5
                0065 real    :: linear_tau      = 0.1
                0066 ! only used with wv_exponent=-1 :
                0067 real    :: ir_tau_co2_win  = 0.2150
                0068 real    :: ir_tau_wv_win1  = 147.11
                0069 real    :: ir_tau_wv_win2  = 1.0814e4
                0070 ! default value set later (according to wv_exponent):
                0071 real    :: ir_tau_co2      = -999.
                0072 real    :: ir_tau_wv       = -999.
fcc5c8b844 jm-c 0073 real    :: ir_tau_wv2      = -999.
cf7369fae0 Jean*0074 real    :: window          = -999.
e148c170bb Jean*0075 !RG: added carbon dioxide concentration to namelist
                0076 real    :: carbon_conc     = 360.0
b2ea1d2979 Jean*0077
                0078 real, save :: pi, deg_to_rad , rad_to_deg
                0079
9f7d507fb9 Jean*0080 namelist/radiation_nml/ select_incSW, solar_constant, del_sol, &
b2ea1d2979 Jean*0081            ir_tau_eq, ir_tau_pole, linear_tau, ir_tau_co2, ir_tau_wv,   &
fcc5c8b844 jm-c 0082            ir_tau_wv2, atm_abs, sw_diff, del_sw, albedo_value, window,  &
                0083            wv_exponent, solar_exponent, yearLength, yearPhase, obliquity, &
88391fb671 jm-c 0084            sw_co2, ir_tau_co2_win, ir_tau_wv_win1, ir_tau_wv_win2, &
                0085            carbon_conc, ozone_in_SW, two_stream_SW
b2ea1d2979 Jean*0086
                0087 !==================================================================================
                0088 !-------------------- diagnostics fields -------------------------------
                0089
                0090 !integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_lwdn_sfc, id_lwup_sfc, &
                0091 !           id_tdt_rad, id_flux_rad, id_flux_lw, id_flux_sw, id_entrop_rad, id_tdt_sw
                0092
                0093 character(len=10), parameter :: mod_name = 'two_stream'
                0094
                0095 real :: missing_value = -999.
                0096
                0097 contains
                0098
                0099 ! ==================================================================================
                0100 ! ==================================================================================
                0101
                0102 !subroutine radiation_init(is, ie, js, je, num_levels, axes, Time)
c9694dc201 Jean*0103 subroutine radiation_init( is, ie, js, je, num_levels, nSx,nSy, axes,   &
                0104                            Time, cst_albedo, myThid )
b2ea1d2979 Jean*0105
                0106 !-------------------------------------------------------------------------------------
c9694dc201 Jean*0107 integer, intent(in)               :: is, ie, js, je, num_levels
                0108 integer, intent(in)               :: nSx, nSy
b2ea1d2979 Jean*0109 integer, intent(in), dimension(4) :: axes
                0110 !type(time_type), intent(in)       :: Time
                0111 real, intent(in)                  :: Time
c9694dc201 Jean*0112 real, intent(out)                 :: cst_albedo
b2ea1d2979 Jean*0113 integer, intent(in)               :: myThid
                0114 !-------------------------------------------------------------------------------------
31f8711b60 Jean*0115 !integer, dimension(3) :: half = (/1,2,4/)
b2ea1d2979 Jean*0116 !integer :: ierr, io
                0117 integer         :: iUnit
                0118 CHARACTER*(gcm_LEN_MBUF) :: msgBuf
                0119 !-----------------------------------------------------------------------------------------
                0120 ! read namelist and copy to logfile
                0121
                0122 !    _BARRIER
                0123 !    _BEGIN_MASTER(myThid)
                0124      CALL BARRIER(myThid)
                0125      IF ( myThid.EQ.1 ) THEN
                0126
                0127      WRITE(msgBuf,'(A)') 'RADIATION_INIT: opening data.atm_gray'
                0128      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0129      CALL OPEN_COPY_DATA_FILE(                                      &
                0130                            'data.atm_gray', 'RADIATION_INIT',       &
                0131                            iUnit,                                   &
                0132                            myThid )
                0133 !    Read parameters from open data file
                0134      READ(UNIT=iUnit,NML=radiation_nml)
                0135      WRITE(msgBuf,'(A)')                                            &
                0136           'RADIATION_INIT: finished reading data.atm_gray'
                0137      CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid )
                0138 !    Close the open data file
15388b052e Jean*0139 #ifdef SINGLE_DISK_IO
b2ea1d2979 Jean*0140      CLOSE(iUnit)
15388b052e Jean*0141 #else
                0142      CLOSE(iUnit,STATUS='DELETE')
                0143 #endif /* SINGLE_DISK_IO */
b2ea1d2979 Jean*0144
                0145 pi    = 4.0*atan(1.)
                0146 deg_to_rad = 2.*pi/360.
                0147 rad_to_deg = 360.0/2./pi
                0148
cf7369fae0 Jean*0149 if ( wv_exponent .eq. -1. ) then
                0150 ! default value for wv_exponent=-1:
fcc5c8b844 jm-c 0151    if ( ir_tau_co2.eq. -999. ) ir_tau_co2 =   3.14
                0152    if ( ir_tau_wv .eq. -999. ) ir_tau_wv  = 199.25
                0153    if ( ir_tau_wv2.eq. -999. ) ir_tau_wv2 =  14.78
cf7369fae0 Jean*0154    if ( window    .eq. -999. ) window     = 0.3732 ! spectral window for Water Vapour
                0155 else
                0156 ! default value for wv_exponent= 0:
                0157    if ( ir_tau_co2.eq. -999. ) ir_tau_co2 = 0.8678
                0158    if ( ir_tau_wv .eq. -999. ) ir_tau_wv  = 1.9979e+3
                0159    if ( window    .eq. -999. ) window     = 0.0    ! spectral window transparent to LW
                0160 endif
                0161
b2ea1d2979 Jean*0162 initialized = .true.
                0163
                0164      ENDIF
                0165      CALL BARRIER(myThid)
c9694dc201 Jean*0166      cst_albedo = albedo_value
b2ea1d2979 Jean*0167
                0168 !-----------------------------------------------------------------------
                0169 !------------ initialize diagnostic fields ---------------
                0170
                0171 !   id_olr = &
                0172 !   register_diag_field ( mod_name, 'olr', axes(1:2), Time, &
                0173 !              'outgoing longwave radiation', &
                0174 !              'watts/m2', missing_value=missing_value               )
                0175 !   id_swdn_sfc = &
                0176 !   register_diag_field ( mod_name, 'swdn_sfc', axes(1:2), Time, &
                0177 !              'SW flux down at surface', &
                0178 !              'watts/m2', missing_value=missing_value               )
                0179 !   id_swdn_toa = &
                0180 !   register_diag_field ( mod_name, 'swdn_toa', axes(1:2), Time, &
                0181 !              'SW flux down at TOA', &
                0182 !              'watts/m2', missing_value=missing_value               )
                0183 !   id_lwup_sfc = &
                0184 !   register_diag_field ( mod_name, 'lwup_sfc', axes(1:2), Time, &
                0185 !              'LW flux up at surface', &
                0186 !              'watts/m2', missing_value=missing_value               )
                0187 !   id_lwdn_sfc = &
                0188 !   register_diag_field ( mod_name, 'lwdn_sfc', axes(1:2), Time, &
                0189 !              'LW flux down at surface', &
                0190 !              'watts/m2', missing_value=missing_value               )
                0191 !   id_tdt_rad = &
                0192 !       register_diag_field ( mod_name, 'tdt_rad', axes(1:3), Time, &
                0193 !              'Temperature tendency due to radiation', &
                0194 !              'K/s', missing_value=missing_value               )
                0195 !   id_tdt_sw  = &
                0196 !       register_diag_field ( mod_name, 'tdt_sw', axes(1:3), Time, &
                0197 !              'Temperature tendency due to SW radiation', &
                0198 !              'K/s', missing_value=missing_value               )
                0199 !   id_flux_rad = &
                0200 !       register_diag_field ( mod_name, 'flux_rad', axes(half), Time, &
                0201 !              'Total radiative flux (positive up)', &
                0202 !              'W/m^2', missing_value=missing_value               )
                0203 !   id_flux_lw = &
                0204 !       register_diag_field ( mod_name, 'flux_lw', axes(half), Time, &
                0205 !              'Net longwave radiative flux (positive up)', &
                0206 !              'W/m^2', missing_value=missing_value               )
                0207 !   id_flux_sw = &
                0208 !       register_diag_field ( mod_name, 'flux_sw', axes(half), Time, &
                0209 !              'Net shortwave radiative flux (positive up)', &
                0210 !              'W/m^2', missing_value=missing_value               )
                0211 !   id_entrop_rad = &
                0212 !           register_diag_field ( mod_name, 'entrop_rad', axes(1:3), Time, &
                0213 !              'Entropy production by radiation', &
                0214 !              '1/s', missing_value=missing_value               )
                0215
                0216 return
                0217 end subroutine radiation_init
                0218
                0219 ! ==================================================================================
                0220
                0221 subroutine radiation_down (is, js, Time_diag, lat, p_half, t, q,      &
88391fb671 jm-c 0222                            albedo, o3conc, ozone_column,              &
                0223                            net_surf_sw_down, surf_lw_down,            &
c9694dc201 Jean*0224                            dtrans, dtrans_win, b, b_win,              &
cf7369fae0 Jean*0225                            down, solar_down,                          &
b2ea1d2979 Jean*0226                            myThid )
                0227
                0228 ! Begin the radiation calculation by computing downward fluxes.
                0229 ! This part of the calculation does not depend on the surface temperature.
                0230
                0231 integer, intent(in)                 :: is, js
                0232 !type(time_type), intent(in)         :: Time_diag
                0233 real, intent(in)                    :: Time_diag
                0234 real, intent(in) , dimension(:,:)   :: lat
cf7369fae0 Jean*0235 real, intent(out), dimension(:,:)   :: net_surf_sw_down
                0236 real, intent(out), dimension(:,:)   :: surf_lw_down
b2ea1d2979 Jean*0237 real, intent(in) , dimension(:,:,:) :: t, q, p_half
c9694dc201 Jean*0238 real, intent(in) , dimension(:,:)   :: albedo
88391fb671 jm-c 0239 real, intent(in) , dimension(:,:,:) :: o3conc
                0240 real, intent(out), dimension(:,:,:) :: ozone_column
b2ea1d2979 Jean*0241 real, intent(out), dimension(:,:,:) :: dtrans
cf7369fae0 Jean*0242 real, intent(out), dimension(:,:,:) :: dtrans_win
b2ea1d2979 Jean*0243 real, intent(out), dimension(:,:,:) :: b
cf7369fae0 Jean*0244 real, intent(out), dimension(:,:,:) :: b_win
b2ea1d2979 Jean*0245 real, intent(out), dimension(:,:,:) :: down
                0246 real, intent(out), dimension(:,:,:) :: solar_down
                0247 integer, intent(in)                 :: myThid
                0248
31f8711b60 Jean*0249 !integer :: i, j
                0250 integer :: k, n
b2ea1d2979 Jean*0251 integer :: im, jm
9f7d507fb9 Jean*0252 ! Variables for seasonal Incoming SW
                0253 real    :: tYear, largeTan, cDecl, sDecl, tanDecl
e148c170bb Jean*0254 real    :: lgCO2conc
b2ea1d2979 Jean*0255
31f8711b60 Jean*0256 !logical :: used
b2ea1d2979 Jean*0257
                0258 ! -------------------------------------------------------------------------
                0259 !real, allocatable, dimension(:,:)   :: swin
9f7d507fb9 Jean*0260 real, allocatable, dimension(:,:)   :: ss, solar, solar_tau_0, p2
cf7369fae0 Jean*0261 real, allocatable, dimension(:,:)   :: solar_tau_k, sw_wv, del_sol_tau, mag_fac
                0262 real, allocatable, dimension(:,:,:) :: solar_tau, dtrans_sol, down_win
                0263 real, allocatable, dimension(:,:)   :: del_tau, tau_0, tau_km, tau_kp, del_tau_win
88391fb671 jm-c 0264 ! MK: Variables for Lacis & Hansen ozone SW scheme
                0265 real, allocatable, dimension(:,:,:) :: ozone_mag, ozone_dF0_down
                0266 real, allocatable, dimension(:,:)   :: abs_uv_LH, abs_uv_LH_FS
                0267 real, allocatable, dimension(:,:)   :: abs_vis_LH, abs_vis_LH_FS
9f7d507fb9 Jean*0268 ! Variables for seasonal Incoming SW
                0269 real, allocatable, dimension(:,:)   :: cLat, sLat, cos_H, HourAng
b2ea1d2979 Jean*0270 ! -------------------------------------------------------------------------
                0271
                0272 n = size(t,3)
                0273 im = size(t,1)
                0274 jm = size(t,2)
                0275
                0276 ! -------------------------------------------------------------------------
                0277 !allocate (swin             (im, jm))
                0278 allocate (ss               (im, jm))
                0279 allocate (solar            (im, jm))
9f7d507fb9 Jean*0280 if ( select_incSW .eq. 0 ) then
                0281   allocate (p2             (im, jm))
88391fb671 jm-c 0282 else
9f7d507fb9 Jean*0283 ! Variables for seasonal Incoming SW
                0284   allocate (cLat           (im, jm))
                0285   allocate (sLat           (im, jm))
                0286   allocate (cos_H          (im, jm))
                0287   allocate (HourAng        (im, jm))
                0288 endif
88391fb671 jm-c 0289 if ( ozone_in_SW .eq. 1 ) then
                0290   allocate (ozone_mag      (im, jm, n+1))
                0291   allocate (ozone_dF0_down (im, jm, n+1))
                0292   allocate (abs_uv_LH      (im, jm))
                0293   allocate (abs_vis_LH     (im, jm))
                0294   allocate (abs_uv_LH_FS   (im, jm))
                0295   allocate (abs_vis_LH_FS  (im, jm))
                0296 endif
cf7369fae0 Jean*0297 if ( solar_exponent .eq. 0. ) then
                0298   allocate (solar_tau_k    (im, jm))
                0299   allocate (sw_wv          (im, jm))
                0300   allocate (del_sol_tau    (im, jm))
                0301   allocate (mag_fac        (im, jm))
                0302   allocate (dtrans_sol     (im, jm, n))
                0303 else
                0304   allocate (solar_tau_0    (im, jm))
                0305   allocate (solar_tau      (im, jm, n+1))
                0306 endif
                0307 if ( wv_exponent .eq. -1. ) then
                0308   allocate (del_tau        (im, jm))
                0309   allocate (del_tau_win    (im, jm))
                0310   allocate (down_win       (im, jm, n+1))
                0311 elseif ( wv_exponent .eq. 0. ) then
9f7d507fb9 Jean*0312   allocate (del_tau        (im, jm))
b2ea1d2979 Jean*0313 else
9f7d507fb9 Jean*0314   allocate (tau_0          (im, jm))
                0315   allocate (tau_km         (im, jm))
                0316   allocate (tau_kp         (im, jm))
b2ea1d2979 Jean*0317 endif
                0318 ! -------------------------------------------------------------------------
                0319
                0320 ss  = sin(lat)
                0321
88391fb671 jm-c 0322 if ( select_incSW .eq. 1 .or. ozone_in_SW .eq. 1 ) then
fcc5c8b844 jm-c 0323 ! daily-mean Incoming SW with simple seasonal cycle accounting
9f7d507fb9 Jean*0324 !  only for obliquity (i.e., circular orbit planet)
                0325
                0326 ! tYear = time in the solar year, in [0-1]
                0327    tYear = MOD( Time_diag/yearLength + yearPhase , 1.0 )
                0328
                0329 ! Compute the declination angle
                0330 ! a) approximate estimate of declination angle: relative error is less
                0331 !    than 0.03 for current obliq but as large as 0.21 for obliq=60^o
88391fb671 jm-c 0332 !    xDecl = - obliquity*deg_to_rad * cos(2.*pi*tYear)
fcc5c8b844 jm-c 0333 ! b) unapproximate expression:
88391fb671 jm-c 0334    sDecl = -sin( obliquity*deg_to_rad ) * cos(2.*pi*tYear)
9f7d507fb9 Jean*0335    cDecl =  cos( asin( sDecl ) )
88391fb671 jm-c 0336 endif
                0337
                0338 if ( select_incSW .eq. 0 ) then
                0339 ! Original Incoming SW (no seasonal cycle):
                0340    p2 = (1. - 3.*ss*ss)/4.
                0341    solar = 0.25*solar_constant*(1.0 + del_sol*p2 + del_sw * ss)
                0342
                0343 elseif ( select_incSW .eq. 1 ) then
                0344
                0345    largeTan = 1.e+16
9f7d507fb9 Jean*0346    if ( cDecl.EQ.0. ) then
                0347     tanDecl = sign( largeTan, sDecl )
                0348    else
                0349     tanDecl = sDecl / cDecl
                0350    endif
                0351
                0352 ! Compute Insolation
                0353    cLat = cos(lat)
                0354 !  sLat = sin(lat)
                0355    sLat = ss
                0356    cos_H = sign( largeTan, sLat )
                0357    where ( cLat .ne. 0. )
                0358      cos_H = sLat/cLat
                0359    endwhere
                0360    cos_H = max( min( - cos_H * tanDecl, 1. ) , -1. )
                0361    HourAng = acos( cos_H )
                0362    solar = (solar_constant/pi)                             &
                0363          *( HourAng*sLat*sDecl + cLat*cDecl*sin(HourAng) )
                0364 else
                0365   stop 'invalid select_incSW'
                0366 endif
b2ea1d2979 Jean*0367
e148c170bb Jean*0368 ! set the log of CO2 concentration
                0369 lgCO2conc = log(carbon_conc/360.)
                0370
b2ea1d2979 Jean*0371 ! set a constant albedo for testing
c9694dc201 Jean*0372 !albedo(:,:) = albedo_value
b2ea1d2979 Jean*0373
cf7369fae0 Jean*0374 if ( solar_exponent .eq. 0. ) then
                0375 ! (RG) scheme based on dtau/dsigma = bq + amu where b is a function of solar_tau
                0376 ! ref: Ruth Geen etal, GRL 2016 (supp. information).
e148c170bb Jean*0377 ! RG: added carbon dioxide log function
cf7369fae0 Jean*0378   solar_tau_k = 0.
                0379   mag_fac(:,:) = 1. ! 35.0 / sqrt(1224 * (cos(lat(:,:))) ** 2 + 1)
                0380
                0381   do k = 1, n
                0382 !   sw_wv = exp( 0.01887 / (solar_tau_k + 0.009522)                            &
                0383 !              + 1.603 / ( (solar_tau_k + 0.5194)**2 ) )
                0384     sw_wv = solar_tau_k + 0.5194
                0385     sw_wv = exp( 0.01887 / (solar_tau_k + 0.009522)                            &
                0386                  + 1.603 / ( sw_wv*sw_wv ) )
e148c170bb Jean*0387     del_sol_tau(:,:) = ( sw_co2 + 0.0029 * lgCO2conc                           &
                0388                                 + sw_wv(:,:) * q(:,:,k) ) * mag_fac(:,:)       &
cf7369fae0 Jean*0389                      * ( p_half(:,:,k+1) - p_half(:,:,k) ) / p_half(:,:,n+1)
                0390     dtrans_sol(:,:,k) = exp( - del_sol_tau(:,:) )
                0391     solar_tau_k = solar_tau_k + del_sol_tau(:,:)
                0392   end do
                0393
                0394   solar_down(:,:,1) = solar(:,:)
                0395   do k = 1,n
                0396     solar_down(:,:,k+1) = solar_down(:,:,k)*dtrans_sol(:,:,k)
                0397   end do
                0398
                0399 else
                0400 ! tau proportional to p^solar_exponent, scaled by atm_abs
                0401
                0402   solar_tau_0 = (1.0 - sw_diff*ss*ss)*atm_abs
                0403
                0404   do k = 1, n+1
                0405     solar_tau(:,:,k) = solar_tau_0(:,:)                                        &
                0406                      *(p_half(:,:,k)/p_half(:,:,n+1))**solar_exponent
                0407   end do
                0408
                0409   do k = 1,n+1
                0410     solar_down(:,:,k) = solar(:,:)*exp(-solar_tau(:,:,k))
                0411   end do
                0412
                0413 endif
                0414
88391fb671 jm-c 0415 if ( ozone_in_SW .eq. 1 ) then
                0416 !-------
                0417 ! MK 2017: Add Lacis and Hansen (JAS, 1974) absorption by ozone,
                0418 ! and apply Forster and Shine (JGR, 1997) correction.
                0419 !-------
                0420 ! Computes broadband absorption, as fraction of total solar flux at TOA.
                0421 ! However, O3 band does not significantly overlap with water vapour SW absorption,
                0422 ! so can safely add O3 heating on top of water vapour heating without changing the
                0423 ! radiation seen by the water vapour optical depths in Ruth's scheme above.
                0424
                0425 ! o3conc = Ozone VMR on pressure levels.  Need to convert to centimetre of ozone
                0426 ! to use with Lacis and Hansen scheme.  Make ozone_column array containing
                0427 ! ozone amount in cm from TOA down to top of level k
                0428
                0429 ! ozone_column(:,:,1:n+1) = Column O3 in cm above level 1 <= k <= n+1
                0430
                0431   ozone_column(:,:,1) = 0.
                0432   do k = 1, n
                0433     ozone_column(:,:,k+1) = ozone_column(:,:,k) + (                            &
                0434                           ( 287.1 / (1.38065e-23 * grav) )                     &
                0435                           * ( p_half(:,:,k+1) - p_half(:,:,k) )                &
                0436                           * o3conc(:,:,k) / 2.687e+23                          &
                0437                                                   )
                0438   end do
                0439
                0440 ! Apply magnification factor:
                0441   ozone_mag(:,:,:) = ozone_column(:,:,:) * 35.                                  &
                0442                    / sqrt( 1224.*cDecl*cDecl + 1. )
                0443
                0444 ! Parameterization for fraction of total flux absorbed from LH74:
                0445   do k = 1, n+1
                0446
                0447     abs_vis_LH(:,:)   = ( 0.02118 * ozone_mag(:,:,k) ) / (                     &
                0448                           1. + 0.042 * ozone_mag(:,:,k)                        &
                0449                              + 0.000323*ozone_mag(:,:,k)*ozone_mag(:,:,k)      &
                0450                                                          )
                0451   ! FS97 correction:
                0452     abs_vis_LH_FS(:,:) = ( -0.002894 * ozone_mag(:,:,k) + 1.0663 )             &
                0453                        * abs_vis_LH(:,:)
                0454
                0455     abs_uv_LH(:,:)   = ( 1.082 * ozone_mag(:,:,k) ) / (                        &
                0456                          ( 1. + 138.6 * ozone_mag(:,:,k) )**0.805              &
                0457                                                       )                        &
                0458                      + ( 0.0658 * ozone_mag(:,:,k) ) / (                       &
                0459                          1. + ( 103.6 * ozone_mag(:,:,k) )**3                  &
                0460                                                        )
                0461   ! FS97 correction:
                0462     abs_uv_LH_FS(:,:) = ( -0.01632 * ozone_mag(:,:,k) + 1.08964 )              &
                0463                       * abs_uv_LH(:,:)
                0464
                0465   ! Total downward SW flux absorbed above level k due to ozone heating.
                0466   ! N.B. Flux absorbed in level k = ozone_dF0_down(:,:,k+1) - ozone_dF0_down(:,:,k)
                0467     ozone_dF0_down(:,:,k) = solar(:,:) * ( abs_vis_LH_FS(:,:)                  &
                0468                                          + abs_uv_LH_FS(:,:)                   &
                0469                                          )
                0470
                0471   ! Subtract off ozone absorption from total SW flux:
                0472     solar_down(:,:,k) = solar_down(:,:,k) - ozone_dF0_down(:,:,k)
                0473
                0474   end do
                0475
                0476 endif
                0477
cf7369fae0 Jean*0478 if ( wv_exponent .eq. -1. ) then
                0479 ! split LW in 2 bands: water-vapour window and remaining = non-window:
                0480 ! ref: Ruth Geen etal, GRL 2016 (supp. information).
e148c170bb Jean*0481 ! RG: added carbon dioxide log function
cf7369fae0 Jean*0482   do k = 1, n
e148c170bb Jean*0483     del_tau    = ( ir_tau_co2 + 0.2023 * lgCO2conc                             &
fcc5c8b844 jm-c 0484                     + ir_tau_wv*log(ir_tau_wv2*q(:,:,k) + 1) )                 &
cf7369fae0 Jean*0485                * ( p_half(:,:,k+1)-p_half(:,:,k) ) / p_half(:,:,n+1)
                0486     dtrans(:,:,k) = exp( - del_tau )
e148c170bb Jean*0487     del_tau_win   = ( ir_tau_co2_win + 0.0954 * lgCO2conc                      &
                0488                                      + ir_tau_wv_win1*q(:,:,k)                 &
cf7369fae0 Jean*0489                                      + ir_tau_wv_win2*q(:,:,k)*q(:,:,k) )      &
                0490                   * ( p_half(:,:,k+1)-p_half(:,:,k) ) / p_half(:,:,n+1)
                0491     dtrans_win(:,:,k) = exp( - del_tau_win )
                0492   end do
b2ea1d2979 Jean*0493
cf7369fae0 Jean*0494 elseif ( wv_exponent .eq. 0. ) then
                0495   dtrans_win = 1.
b2ea1d2979 Jean*0496 ! longwave optical thickness function of specific humidity (M.Byrne & P.O'Gorman):
                0497   do k = 1, n
                0498     del_tau    = ( ir_tau_co2 + ir_tau_wv * q(:,:,k) )                         &
                0499                * ( p_half(:,:,k+1)-p_half(:,:,k) ) / p_half(:,:,n+1)
                0500     dtrans(:,:,k) = exp( - del_tau )
                0501   end do
                0502
                0503 else
cf7369fae0 Jean*0504   dtrans_win = 1.
b2ea1d2979 Jean*0505 ! longwave optical thickness function of latitude and pressure
9f7d507fb9 Jean*0506   tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*ss*ss
                0507
b2ea1d2979 Jean*0508   k = 1
                0509   tau_kp =   tau_0(:,:) * (                                                    &
                0510                   linear_tau * p_half(:,:,k)/p_half(:,:,n+1)                   &
                0511          + (1.0 - linear_tau)*(p_half(:,:,k)/p_half(:,:,n+1))**wv_exponent     &
                0512                           )
                0513   do k = 1, n
                0514     tau_km = tau_kp
                0515     tau_kp = tau_0(:,:) * (                                                    &
                0516                   linear_tau * p_half(:,:,k+1)/p_half(:,:,n+1)                 &
                0517          + (1.0 - linear_tau)*(p_half(:,:,k+1)/p_half(:,:,n+1))**wv_exponent   &
                0518                           )
                0519     dtrans(:,:,k) = exp( -(tau_kp - tau_km) )
                0520   end do
                0521
                0522 endif
                0523
                0524 ! no radiation from spectral window
cf7369fae0 Jean*0525 b = stefan*t*t*t*t
                0526 b_win = window*b
                0527 b = (1.0-window)*b
b2ea1d2979 Jean*0528
                0529 down(:,:,1) = 0.0
                0530 do k = 1,n
cf7369fae0 Jean*0531   down(:,:,k+1)     = down(:,:,k)*dtrans(:,:,k)                        &
                0532                     + b(:,:,k)*(1.0 - dtrans(:,:,k))
b2ea1d2979 Jean*0533 end do
                0534
cf7369fae0 Jean*0535 if ( wv_exponent .eq. -1. ) then
                0536   down_win(:,:,1) = 0.0
                0537   do k = 1,n
                0538     down_win(:,:,k+1) = down_win(:,:,k)*dtrans_win(:,:,k)              &
                0539                       + b_win(:,:,k)*(1.0 - dtrans_win(:,:,k))
                0540   end do
                0541   down = down + down_win
                0542 endif
b2ea1d2979 Jean*0543
                0544 surf_lw_down     = down(:,:,n+1)
                0545 net_surf_sw_down = solar_down(:,:,n+1)*(1. - albedo(:,:))
                0546 !swin = solar_down(:,:,1)
                0547
                0548 !------- downward sw flux surface -------
                0549 !     if ( id_swdn_sfc > 0 ) then
                0550 !         used = send_data ( id_swdn_sfc, net_surf_sw_down, Time_diag)
                0551 !     endif
                0552 !------- incoming sw flux toa -------
                0553 !     if ( id_swdn_toa > 0 ) then
                0554 !         used = send_data ( id_swdn_toa, swin, Time_diag)
                0555 !     endif
                0556 !------- downward lw flux surface -------
                0557 !     if ( id_lwdn_sfc > 0 ) then
                0558 !         used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag)
                0559 !     endif
                0560
                0561 ! -------------------------------------------------------------------------
                0562 !deallocate (swin)
cf7369fae0 Jean*0563 deallocate (ss, solar)
9f7d507fb9 Jean*0564 if ( select_incSW .eq. 0 ) then
                0565   deallocate (p2)
88391fb671 jm-c 0566 else
9f7d507fb9 Jean*0567 ! Variables for seasonal Incoming SW
88391fb671 jm-c 0568   deallocate (cLat, sLat)
                0569   deallocate (cos_H, HourAng)
9f7d507fb9 Jean*0570 endif
cf7369fae0 Jean*0571 if ( solar_exponent .eq. 0. ) then
                0572   deallocate (solar_tau_k, sw_wv, del_sol_tau, mag_fac, dtrans_sol)
                0573 else
                0574   deallocate (solar_tau_0, solar_tau)
                0575 endif
88391fb671 jm-c 0576 if ( ozone_in_SW .eq. 1 ) then
                0577   deallocate (ozone_mag, ozone_dF0_down)
                0578   deallocate (abs_uv_LH, abs_vis_LH, abs_uv_LH_FS, abs_vis_LH_FS)
                0579 endif
cf7369fae0 Jean*0580 if ( wv_exponent .eq. -1. ) then
                0581   deallocate (del_tau, del_tau_win, down_win)
                0582 elseif ( wv_exponent .eq. 0. ) then
b2ea1d2979 Jean*0583   deallocate (del_tau)
                0584 else
9f7d507fb9 Jean*0585   deallocate (tau_0, tau_km, tau_kp)
b2ea1d2979 Jean*0586 endif
                0587 ! -------------------------------------------------------------------------
                0588
                0589 return
                0590 end subroutine radiation_down
                0591
                0592 ! ==================================================================================
                0593
                0594 !subroutine radiation_up (is, js, Time_diag, lat, p_half, t_surf, t, tdt)
88391fb671 jm-c 0595 subroutine radiation_up ( is, js, Time_diag, lat, p_half, t_surf, t, tdt,      &
                0596                           flux_lw, flux_sw, albedo, ozone_column,              &
                0597                           dtrans, dtrans_win, b, b_win, down, solar_down,      &
cf7369fae0 Jean*0598                           myThid )
b2ea1d2979 Jean*0599
                0600 ! Now complete the radiation calculation by computing the upward and net fluxes.
                0601
                0602 integer, intent(in)                 :: is, js
                0603 !type(time_type), intent(in)         :: Time_diag
                0604 real, intent(in)                    :: Time_diag
                0605 real, intent(in) , dimension(:,:)   :: lat
                0606 real, intent(in) , dimension(:,:)   :: t_surf
                0607 real, intent(in) , dimension(:,:,:) :: t, p_half
                0608 real, intent(inout), dimension(:,:,:) :: tdt
88391fb671 jm-c 0609 real, intent(out), dimension(:,:,:) :: flux_lw
cf7369fae0 Jean*0610 real, intent(out), dimension(:,:,:) :: flux_sw
b2ea1d2979 Jean*0611 real, intent(in),  dimension(:,:)   :: albedo
88391fb671 jm-c 0612 real, intent(in) , dimension(:,:,:) :: ozone_column
b2ea1d2979 Jean*0613 real, intent(in),  dimension(:,:,:) :: dtrans
cf7369fae0 Jean*0614 real, intent(in),  dimension(:,:,:) :: dtrans_win
b2ea1d2979 Jean*0615 real, intent(in),  dimension(:,:,:) :: b
cf7369fae0 Jean*0616 real, intent(in),  dimension(:,:,:) :: b_win
b2ea1d2979 Jean*0617 real, intent(in),  dimension(:,:,:) :: down
                0618 real, intent(in),  dimension(:,:,:) :: solar_down
                0619 integer, intent(in)                 :: myThid
                0620
31f8711b60 Jean*0621 !integer :: i, j
                0622 integer :: k, n
b2ea1d2979 Jean*0623 integer :: im, jm
88391fb671 jm-c 0624 ! Variables for seasonal Incoming SW
                0625 real    :: tYear, cDecl, sDecl
b2ea1d2979 Jean*0626
31f8711b60 Jean*0627 !logical :: used
b2ea1d2979 Jean*0628
                0629 ! -------------------------------------------------------------------------
88391fb671 jm-c 0630 real, allocatable, dimension(:,:)   :: solar, b_surf
                0631 real, allocatable, dimension(:,:,:) :: tdt_rad, up, up_win, solar_up
                0632 !real, allocatable, dimension(:,:,:) :: tdt_sw, flux_rad   ! not used
                0633 ! MK: Variables for Lacis & Hansen ozone SW scheme
                0634 real, allocatable, dimension(:,:,:) :: ozone_mag_up, ozone_dF0_up
                0635 real, allocatable, dimension(:,:)   :: abs_uv_LH_up, abs_uv_LH_FS_up
                0636 real, allocatable, dimension(:,:)   :: abs_vis_LH_up, abs_vis_LH_FS_up
                0637 real, allocatable, dimension(:,:)   :: r_bar
b2ea1d2979 Jean*0638 ! -------------------------------------------------------------------------
                0639
                0640 n = size(t,3)
                0641 im = size(t,1)
                0642 jm = size(t,2)
                0643
                0644 ! -------------------------------------------------------------------------
88391fb671 jm-c 0645 if( two_stream_SW .eq. 1 ) then
                0646   allocate (solar            (im, jm))
                0647   if ( ozone_in_SW .eq. 1 ) then
                0648     allocate (ozone_mag_up    (im, jm, n+1))
                0649     allocate (ozone_dF0_up    (im, jm, n+1))
                0650     allocate (abs_uv_LH_up    (im, jm))
                0651     allocate (abs_vis_LH_up   (im, jm))
                0652     allocate (abs_uv_LH_FS_up (im, jm))
                0653     allocate (abs_vis_LH_FS_up(im, jm))
                0654     allocate (r_bar           (im, jm))
                0655   endif
                0656 endif
b2ea1d2979 Jean*0657 allocate (tdt_rad          (im, jm, n))
                0658 allocate (up               (im, jm, n+1))
cf7369fae0 Jean*0659 allocate (up_win           (im, jm, n+1))
88391fb671 jm-c 0660 allocate (solar_up         (im, jm, n+1))
                0661 !allocate (tdt_sw           (im, jm, n))    ! not used
                0662 !allocate (flux_rad         (im, jm, n+1))  ! not used
b2ea1d2979 Jean*0663 allocate (b_surf           (im, jm))
                0664 ! -------------------------------------------------------------------------
                0665
88391fb671 jm-c 0666 ! MK add in solar_up variables
                0667 do k = 1,n+1
                0668   solar_up(:,:,k) = albedo(:,:) * solar_down(:,:,n+1)
                0669 end do
                0670
                0671 if ( two_stream_SW .eq. 1 ) then
                0672 ! MK do upward SW for ozone
                0673
                0674 ! Need solar variables again:
                0675   solar = solar_down(:,:,1)
                0676
                0677   if ( ozone_in_SW .eq. 1 ) then
                0678   ! tYear = time in the solar year, in [0-1]
                0679     tYear = MOD( Time_diag/yearLength + yearPhase , 1.0 )
                0680
                0681   ! Compute the declination angle
                0682   ! a) approximate estimate of declination angle: relative error is less
                0683   !    than 0.03 for current obliq but as large as 0.21 for obliq=60^o
                0684   !    xDecl = - obliquity*deg_to_rad * cos(2.*pi*tYear)
                0685   ! b) unapproximate expression:
                0686     sDecl = -sin(  obliquity*deg_to_rad ) * cos(2.*pi*tYear)
                0687     cDecl =  cos( asin( sDecl ) )
                0688
                0689   ! Lacis and Hansen scheme for absorption of upward SW flux by O3:
                0690   ! Apply magnification factor:
                0691
                0692     do k = n+1,1,-1
                0693       ozone_mag_up(:,:,k) = ozone_column(:,:,n+1) * 35. /                      &
                0694                             sqrt( 1224.*cDecl*cDecl + 1. )                     &
                0695                           + 1.9 * ( ozone_column(:,:,n+1)                      &
                0696                                     - ozone_column(:,:,k)                      &
                0697                                   )
                0698     end do
                0699
                0700     r_bar(:,:) = ( 0.291 / ( 1. + 0.816*cDecl ) ) + (                          &
                0701                    ( 1. - ( 0.291 / ( 1. + 0.816*cDecl ) ) )                   &
                0702                    * 0.856 * albedo(:,:) / ( 1. - 0.144*albedo(:,:) )          &
                0703                                                     )
                0704   ! LH74 parameterisation for fraction of total solar flux absorbed:
                0705
                0706     do k = n+1,1,-1
                0707
                0708       abs_vis_LH_up(:,:) = r_bar(:,:) * (0.02118 * ozone_mag_up(:,:,k)) / (    &
                0709                          1. + 0.042 * ozone_mag_up(:,:,k)                      &
                0710                             + 0.000323*ozone_mag_up(:,:,k)*ozone_mag_up(:,:,k) &
                0711                                                                           )
                0712     ! FS97 correction
                0713       abs_vis_LH_FS_up(:,:) = ( -0.002894 * ozone_mag_up(:,:,k) + 1.0663 )     &
                0714                               * abs_vis_LH_up(:,:)
                0715
                0716       abs_uv_LH_up(:,:) = r_bar(:,:) * (1.082 * ozone_mag_up(:,:,k)) / (       &
                0717                             ( 1. + 138.6 * ozone_mag_up(:,:,k) )**0.805        &
                0718                                                                        )       &
                0719                         + ( 0.0658 * ozone_mag_up(:,:,k) ) / (                 &
                0720                             1. + ( 103.6 * ozone_mag_up(:,:,k) )**3            &
                0721                                                              )
                0722     ! FS97 correction
                0723       abs_uv_LH_FS_up(:,:) = ( -0.01632 * ozone_mag_up(:,:,k) + 1.08964 )      &
                0724                              * abs_uv_LH_up(:,:)
                0725
                0726     ! Total SW flux absorbed along path by reflected radiation up till and including level k.
                0727     ! N.B. Flux absorbed in level k = ozone_dF0_up(:,:,k) - ozone_dF0_up(:,:,k+1)
                0728       ozone_dF0_up(:,:,k) = solar(:,:) * ( abs_vis_LH_FS_up(:,:)               &
                0729                                          + abs_uv_LH_FS_up(:,:)                &
                0730                                          )
                0731
                0732     end do
                0733
                0734     do k = n,1,-1
                0735     ! Subtract off ozone absorption from total SW flux:
                0736       solar_up(:,:,k) = solar_up(:,:,k)                                        &
                0737                       - ( ozone_dF0_up(:,:,k) - ozone_dF0_up(:,:,n+1) )
                0738     end do
                0739
                0740   endif   ! ozone_in_SW
                0741
                0742 endif   ! two_stream_SW
                0743
b2ea1d2979 Jean*0744 ! total flux from surface
                0745 b_surf = stefan*t_surf*t_surf*t_surf*t_surf
                0746
                0747 ! first deal with non-window upward flux
                0748 up(:,:,n+1) = b_surf*(1.0-window)
cf7369fae0 Jean*0749 up_win(:,:,n+1) = b_surf*window
b2ea1d2979 Jean*0750 do k = n,1,-1
                0751   up(:,:,k) = up(:,:,k+1)*dtrans(:,:,k) + b(:,:,k)*(1.0 - dtrans(:,:,k))
cf7369fae0 Jean*0752   up_win(:,:,k) = up_win(:,:,k+1)*dtrans_win(:,:,k)            &
                0753                 + b_win(:,:,k)*(1.0 - dtrans_win(:,:,k))
b2ea1d2979 Jean*0754 end do
                0755
                0756 ! add upward flux in spectral window
cf7369fae0 Jean*0757 up = up + up_win
b2ea1d2979 Jean*0758
88391fb671 jm-c 0759 ! MK Compute flux_sw in terms of solar_up, and use flux_sw in tdt calculation
b2ea1d2979 Jean*0760 do k = 1,n+1
88391fb671 jm-c 0761   flux_lw(:,:,k) = up(:,:,k)-down(:,:,k)   ! LW only, upward +ve
                0762   flux_sw(:,:,k) = solar_down(:,:,k) - solar_up(:,:,k)  ! Net downward SW
                0763 ! flux_rad(:,:,k) = flux_lw(:,:,k) + flux_sw(:,:,k)   ! not used
b2ea1d2979 Jean*0764 end do
                0765
88391fb671 jm-c 0766 if ( two_stream_SW .eq. 1 ) then
                0767   do k = 1,n
                0768     tdt_rad(:,:,k) = ( flux_lw(:,:,k+1) - flux_lw(:,:,k)                  &
                0769                      - flux_sw(:,:,k+1) + flux_sw(:,:,k)                  &
                0770                      )                                                    &
                0771                    * grav / ( cp_air*( p_half(:,:,k+1)-p_half(:,:,k) ) )
                0772   end do
                0773 else
                0774   do k = 1,n
                0775     tdt_rad(:,:,k) = ( flux_lw(:,:,k+1) - flux_lw(:,:,k)                  &
                0776                      - solar_down(:,:,k+1) + solar_down(:,:,k) )          &
                0777                    * grav / ( cp_air*( p_half(:,:,k+1)-p_half(:,:,k) ) )
                0778   end do
                0779 endif
                0780
b2ea1d2979 Jean*0781 do k = 1,n
88391fb671 jm-c 0782 ! tdt_sw is not used:
                0783 ! tdt_sw(:,:,k) = ( - flux_sw(:,:,k+1) + flux_sw(:,:,k) )                 &
                0784 !                  * grav / ( cp_air*( p_half(:,:,k+1)-p_half(:,:,k) ) )
b2ea1d2979 Jean*0785   tdt(:,:,k) = tdt(:,:,k) + tdt_rad(:,:,k)
                0786 end do
                0787
                0788 !------- outgoing lw flux toa (olr) -------
                0789 !     if ( id_olr > 0 ) then
                0790 !         used = send_data ( id_olr, olr, Time_diag)
                0791 !     endif
                0792 !------- upward lw flux surface -------
                0793 !     if ( id_lwup_sfc > 0 ) then
                0794 !         used = send_data ( id_lwup_sfc, b_surf, Time_diag)
                0795 !     endif
                0796 !------- temperature tendency due to radiation ------------
                0797 !     if ( id_tdt_rad > 0 ) then
                0798 !        used = send_data ( id_tdt_rad, tdt_rad, Time_diag)
                0799 !     endif
                0800 !     if ( id_tdt_sw > 0 ) then
                0801 !        used = send_data ( id_tdt_sw, tdt_sw, Time_diag)
                0802 !     endif
                0803 !------- total radiative flux (at half levels) -----------
                0804 !     if ( id_flux_rad > 0 ) then
                0805 !        used = send_data ( id_flux_rad, flux_rad, Time_diag)
                0806 !     endif
                0807 !------- longwave radiative flux (at half levels) --------
                0808 !     if ( id_flux_lw > 0 ) then
                0809 !        used = send_data ( id_flux_lw, net, Time_diag)
                0810 !     endif
                0811 !     if ( id_flux_sw > 0 ) then
                0812 !        used = send_data ( id_flux_sw, flux_sw, Time_diag)
                0813 !     endif
                0814 !     if ( id_entrop_rad > 0 ) then
                0815 !        do k=1,n
                0816 !           entrop_rad(:,:,k) =tdt_rad(:,:,k)/t(:,:,k)*p_half(:,:,n+1)/1.e5
                0817 !        end do
                0818 !        used = send_data ( id_entrop_rad, entrop_rad, Time_diag)
                0819 !     endif
                0820
                0821 ! -------------------------------------------------------------------------
88391fb671 jm-c 0822 if ( two_stream_SW .eq. 1 ) then
                0823   deallocate (solar)
                0824   if ( ozone_in_SW .eq. 1 ) then
                0825     deallocate (ozone_mag_up, ozone_dF0_up)
                0826     deallocate (abs_uv_LH_up, abs_vis_LH_up)
                0827     deallocate (abs_uv_LH_FS_up, abs_vis_LH_FS_up, r_bar)
                0828   endif
                0829 endif
                0830 deallocate (tdt_rad, up, up_win, solar_up)
                0831 !deallocate (tdt_sw, flux_rad)  ! not used
b2ea1d2979 Jean*0832 deallocate (b_surf)
                0833 ! -------------------------------------------------------------------------
                0834
                0835 return
                0836 end subroutine radiation_up
                0837
                0838 ! ==================================================================================
                0839
                0840 subroutine radiation_end
                0841
                0842 !deallocate (b, tdt_rad, tdt_sw, entrop_rad)
                0843 !!deallocate (up, down, net, solar_down, flux_rad, flux_sw)
                0844 !deallocate (up, net, flux_rad, flux_sw)
                0845 !deallocate (b_surf, olr, swin, albedo)
                0846 !!deallocate (dtrans)
                0847 !deallocate (tau, solar_tau)
                0848 !deallocate (ss, solar, tau_0, solar_tau_0, p2)
                0849
                0850 end subroutine radiation_end
                0851
                0852 ! ==================================================================================
                0853
                0854 end module radiation_mod