| 1 |
module radiation_mod |
| 2 |
|
| 3 |
! ================================================================================== |
| 4 |
! ================================================================================== |
| 5 |
|
| 6 |
! use fms_mod, only: open_file, check_nml_error, & |
| 7 |
! mpp_pe, close_file |
| 8 |
|
| 9 |
use gcm_params_mod, only: gcm_LEN_MBUF, gcm_SQZ_R, gcm_stdMsgUnit |
| 10 |
use constants_mod, only: stefan, cp_air, grav, pstd_mks |
| 11 |
|
| 12 |
! use diag_manager_mod, only: register_diag_field, send_data |
| 13 |
|
| 14 |
! use time_manager_mod, only: time_type, & |
| 15 |
! operator(+), operator(-), operator(/=) |
| 16 |
|
| 17 |
!================================================================================== |
| 18 |
implicit none |
| 19 |
private |
| 20 |
!================================================================================== |
| 21 |
|
| 22 |
! version information |
| 23 |
|
| 24 |
character(len=128) :: version='$Id: radiation_mod.F90,v 1.3 2012/10/16 15:21:39 jmc Exp $' |
| 25 |
character(len=128) :: tag='homemade' |
| 26 |
|
| 27 |
!================================================================================== |
| 28 |
|
| 29 |
! public interfaces |
| 30 |
|
| 31 |
public :: radiation_init, radiation_down, radiation_up, radiation_end |
| 32 |
!================================================================================== |
| 33 |
|
| 34 |
! module variables |
| 35 |
logical :: initialized =.false. |
| 36 |
|
| 37 |
real :: solar_constant = 1360.0 |
| 38 |
real :: del_sol = 1.4 |
| 39 |
! modif omp: winter/summer hemisphere |
| 40 |
real :: del_sw = 0.0 |
| 41 |
real :: ir_tau_eq = 6.0 |
| 42 |
real :: ir_tau_pole = 1.5 |
| 43 |
real :: linear_tau = 0.1 |
| 44 |
real :: ir_tau_co2 = 0.8678 |
| 45 |
real :: ir_tau_wv = 1.9979e+3 |
| 46 |
real :: atm_abs = 0.0 |
| 47 |
real :: sw_diff = 0.0 |
| 48 |
real :: albedo_value = 0.06 |
| 49 |
real :: window = 0.0 ! spectral window transparent to LW |
| 50 |
real :: wv_exponent = 4.0 |
| 51 |
real :: solar_exponent = 4.0 |
| 52 |
|
| 53 |
real, save :: pi, deg_to_rad , rad_to_deg |
| 54 |
|
| 55 |
namelist/radiation_nml/ solar_constant, del_sol, & |
| 56 |
ir_tau_eq, ir_tau_pole, linear_tau, ir_tau_co2, ir_tau_wv, & |
| 57 |
atm_abs, sw_diff, del_sw, albedo_value, window, wv_exponent, & |
| 58 |
solar_exponent |
| 59 |
|
| 60 |
!================================================================================== |
| 61 |
!-------------------- diagnostics fields ------------------------------- |
| 62 |
|
| 63 |
!integer :: id_olr, id_swdn_sfc, id_swdn_toa, id_lwdn_sfc, id_lwup_sfc, & |
| 64 |
! id_tdt_rad, id_flux_rad, id_flux_lw, id_flux_sw, id_entrop_rad, id_tdt_sw |
| 65 |
|
| 66 |
character(len=10), parameter :: mod_name = 'two_stream' |
| 67 |
|
| 68 |
real :: missing_value = -999. |
| 69 |
|
| 70 |
contains |
| 71 |
|
| 72 |
! ================================================================================== |
| 73 |
! ================================================================================== |
| 74 |
|
| 75 |
!subroutine radiation_init(is, ie, js, je, num_levels, axes, Time) |
| 76 |
subroutine radiation_init(is, ie, js, je, num_levels, nSx,nSy, axes, Time, myThid) |
| 77 |
|
| 78 |
!------------------------------------------------------------------------------------- |
| 79 |
integer, intent(in), dimension(4) :: axes |
| 80 |
!type(time_type), intent(in) :: Time |
| 81 |
real, intent(in) :: Time |
| 82 |
integer, intent(in) :: is, ie, js, je, num_levels |
| 83 |
integer, intent(in) :: nSx, nSy |
| 84 |
integer, intent(in) :: myThid |
| 85 |
!------------------------------------------------------------------------------------- |
| 86 |
integer, dimension(3) :: half = (/1,2,4/) |
| 87 |
!integer :: ierr, io |
| 88 |
integer :: iUnit |
| 89 |
CHARACTER*(gcm_LEN_MBUF) :: msgBuf |
| 90 |
!----------------------------------------------------------------------------------------- |
| 91 |
! read namelist and copy to logfile |
| 92 |
|
| 93 |
! _BARRIER |
| 94 |
! _BEGIN_MASTER(myThid) |
| 95 |
CALL BARRIER(myThid) |
| 96 |
IF ( myThid.EQ.1 ) THEN |
| 97 |
|
| 98 |
WRITE(msgBuf,'(A)') 'RADIATION_INIT: opening data.atm_gray' |
| 99 |
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
| 100 |
CALL OPEN_COPY_DATA_FILE( & |
| 101 |
'data.atm_gray', 'RADIATION_INIT', & |
| 102 |
iUnit, & |
| 103 |
myThid ) |
| 104 |
! Read parameters from open data file |
| 105 |
READ(UNIT=iUnit,NML=radiation_nml) |
| 106 |
WRITE(msgBuf,'(A)') & |
| 107 |
'RADIATION_INIT: finished reading data.atm_gray' |
| 108 |
CALL PRINT_MESSAGE( msgBuf, gcm_stdMsgUnit, gcm_SQZ_R, myThid ) |
| 109 |
! Close the open data file |
| 110 |
CLOSE(iUnit) |
| 111 |
|
| 112 |
pi = 4.0*atan(1.) |
| 113 |
deg_to_rad = 2.*pi/360. |
| 114 |
rad_to_deg = 360.0/2./pi |
| 115 |
|
| 116 |
initialized = .true. |
| 117 |
|
| 118 |
ENDIF |
| 119 |
CALL BARRIER(myThid) |
| 120 |
|
| 121 |
!----------------------------------------------------------------------- |
| 122 |
!------------ initialize diagnostic fields --------------- |
| 123 |
|
| 124 |
! id_olr = & |
| 125 |
! register_diag_field ( mod_name, 'olr', axes(1:2), Time, & |
| 126 |
! 'outgoing longwave radiation', & |
| 127 |
! 'watts/m2', missing_value=missing_value ) |
| 128 |
! id_swdn_sfc = & |
| 129 |
! register_diag_field ( mod_name, 'swdn_sfc', axes(1:2), Time, & |
| 130 |
! 'SW flux down at surface', & |
| 131 |
! 'watts/m2', missing_value=missing_value ) |
| 132 |
! id_swdn_toa = & |
| 133 |
! register_diag_field ( mod_name, 'swdn_toa', axes(1:2), Time, & |
| 134 |
! 'SW flux down at TOA', & |
| 135 |
! 'watts/m2', missing_value=missing_value ) |
| 136 |
! id_lwup_sfc = & |
| 137 |
! register_diag_field ( mod_name, 'lwup_sfc', axes(1:2), Time, & |
| 138 |
! 'LW flux up at surface', & |
| 139 |
! 'watts/m2', missing_value=missing_value ) |
| 140 |
! id_lwdn_sfc = & |
| 141 |
! register_diag_field ( mod_name, 'lwdn_sfc', axes(1:2), Time, & |
| 142 |
! 'LW flux down at surface', & |
| 143 |
! 'watts/m2', missing_value=missing_value ) |
| 144 |
! id_tdt_rad = & |
| 145 |
! register_diag_field ( mod_name, 'tdt_rad', axes(1:3), Time, & |
| 146 |
! 'Temperature tendency due to radiation', & |
| 147 |
! 'K/s', missing_value=missing_value ) |
| 148 |
! id_tdt_sw = & |
| 149 |
! register_diag_field ( mod_name, 'tdt_sw', axes(1:3), Time, & |
| 150 |
! 'Temperature tendency due to SW radiation', & |
| 151 |
! 'K/s', missing_value=missing_value ) |
| 152 |
! id_flux_rad = & |
| 153 |
! register_diag_field ( mod_name, 'flux_rad', axes(half), Time, & |
| 154 |
! 'Total radiative flux (positive up)', & |
| 155 |
! 'W/m^2', missing_value=missing_value ) |
| 156 |
! id_flux_lw = & |
| 157 |
! register_diag_field ( mod_name, 'flux_lw', axes(half), Time, & |
| 158 |
! 'Net longwave radiative flux (positive up)', & |
| 159 |
! 'W/m^2', missing_value=missing_value ) |
| 160 |
! id_flux_sw = & |
| 161 |
! register_diag_field ( mod_name, 'flux_sw', axes(half), Time, & |
| 162 |
! 'Net shortwave radiative flux (positive up)', & |
| 163 |
! 'W/m^2', missing_value=missing_value ) |
| 164 |
! id_entrop_rad = & |
| 165 |
! register_diag_field ( mod_name, 'entrop_rad', axes(1:3), Time, & |
| 166 |
! 'Entropy production by radiation', & |
| 167 |
! '1/s', missing_value=missing_value ) |
| 168 |
|
| 169 |
return |
| 170 |
end subroutine radiation_init |
| 171 |
|
| 172 |
! ================================================================================== |
| 173 |
|
| 174 |
subroutine radiation_down (is, js, Time_diag, lat, p_half, t, q, & |
| 175 |
! net_surf_sw_down, surf_lw_down) |
| 176 |
net_surf_sw_down, surf_lw_down, & |
| 177 |
albedo, dtrans, b, down, solar_down, & |
| 178 |
myThid ) |
| 179 |
|
| 180 |
! Begin the radiation calculation by computing downward fluxes. |
| 181 |
! This part of the calculation does not depend on the surface temperature. |
| 182 |
|
| 183 |
integer, intent(in) :: is, js |
| 184 |
!type(time_type), intent(in) :: Time_diag |
| 185 |
real, intent(in) :: Time_diag |
| 186 |
real, intent(in) , dimension(:,:) :: lat |
| 187 |
real, intent(out) , dimension(:,:) :: net_surf_sw_down |
| 188 |
real, intent(out) , dimension(:,:) :: surf_lw_down |
| 189 |
real, intent(in) , dimension(:,:,:) :: t, q, p_half |
| 190 |
real, intent(out), dimension(:,:) :: albedo |
| 191 |
real, intent(out), dimension(:,:,:) :: dtrans |
| 192 |
real, intent(out), dimension(:,:,:) :: b |
| 193 |
real, intent(out), dimension(:,:,:) :: down |
| 194 |
real, intent(out), dimension(:,:,:) :: solar_down |
| 195 |
integer, intent(in) :: myThid |
| 196 |
|
| 197 |
integer :: i, j, k, n |
| 198 |
integer :: im, jm |
| 199 |
|
| 200 |
logical :: used |
| 201 |
|
| 202 |
! ------------------------------------------------------------------------- |
| 203 |
!real, allocatable, dimension(:,:) :: swin |
| 204 |
real, allocatable, dimension(:,:) :: ss, solar, tau_0, solar_tau_0, p2 |
| 205 |
real, allocatable, dimension(:,:,:) :: solar_tau |
| 206 |
real, allocatable, dimension(:,:) :: del_tau, tau_km, tau_kp |
| 207 |
! ------------------------------------------------------------------------- |
| 208 |
|
| 209 |
n = size(t,3) |
| 210 |
im = size(t,1) |
| 211 |
jm = size(t,2) |
| 212 |
|
| 213 |
! ------------------------------------------------------------------------- |
| 214 |
!allocate (swin (im, jm)) |
| 215 |
allocate (ss (im, jm)) |
| 216 |
allocate (solar (im, jm)) |
| 217 |
allocate (tau_0 (im, jm)) |
| 218 |
allocate (solar_tau_0 (im, jm)) |
| 219 |
allocate (p2 (im, jm)) |
| 220 |
allocate (solar_tau (im, jm, n+1)) |
| 221 |
if ( wv_exponent .eq. 0. ) then |
| 222 |
allocate (del_tau (im, jm)) |
| 223 |
else |
| 224 |
allocate (tau_km (im, jm)) |
| 225 |
allocate (tau_kp (im, jm)) |
| 226 |
endif |
| 227 |
! ------------------------------------------------------------------------- |
| 228 |
|
| 229 |
ss = sin(lat) |
| 230 |
p2 = (1. - 3.*ss*ss)/4. |
| 231 |
|
| 232 |
solar = 0.25*solar_constant*(1.0 + del_sol*p2 + del_sw * ss) |
| 233 |
|
| 234 |
tau_0 = ir_tau_eq + (ir_tau_pole - ir_tau_eq)*ss*ss |
| 235 |
|
| 236 |
solar_tau_0 = (1.0 - sw_diff*ss*ss)*atm_abs |
| 237 |
|
| 238 |
! set a constant albedo for testing |
| 239 |
albedo(:,:) = albedo_value |
| 240 |
|
| 241 |
do k = 1, n+1 |
| 242 |
solar_tau(:,:,k) = solar_tau_0(:,:) & |
| 243 |
*(p_half(:,:,k)/p_half(:,:,n+1))**solar_exponent |
| 244 |
end do |
| 245 |
|
| 246 |
if ( wv_exponent .eq. 0. ) then |
| 247 |
! longwave optical thickness function of specific humidity (M.Byrne & P.O'Gorman): |
| 248 |
do k = 1, n |
| 249 |
del_tau = ( ir_tau_co2 + ir_tau_wv * q(:,:,k) ) & |
| 250 |
* ( p_half(:,:,k+1)-p_half(:,:,k) ) / p_half(:,:,n+1) |
| 251 |
dtrans(:,:,k) = exp( - del_tau ) |
| 252 |
end do |
| 253 |
|
| 254 |
else |
| 255 |
! longwave optical thickness function of latitude and pressure |
| 256 |
k = 1 |
| 257 |
tau_kp = tau_0(:,:) * ( & |
| 258 |
linear_tau * p_half(:,:,k)/p_half(:,:,n+1) & |
| 259 |
+ (1.0 - linear_tau)*(p_half(:,:,k)/p_half(:,:,n+1))**wv_exponent & |
| 260 |
) |
| 261 |
do k = 1, n |
| 262 |
tau_km = tau_kp |
| 263 |
tau_kp = tau_0(:,:) * ( & |
| 264 |
linear_tau * p_half(:,:,k+1)/p_half(:,:,n+1) & |
| 265 |
+ (1.0 - linear_tau)*(p_half(:,:,k+1)/p_half(:,:,n+1))**wv_exponent & |
| 266 |
) |
| 267 |
dtrans(:,:,k) = exp( -(tau_kp - tau_km) ) |
| 268 |
end do |
| 269 |
|
| 270 |
endif |
| 271 |
|
| 272 |
! no radiation from spectral window |
| 273 |
b = (1.0-window)*stefan*t*t*t*t |
| 274 |
|
| 275 |
down(:,:,1) = 0.0 |
| 276 |
do k = 1,n |
| 277 |
down(:,:,k+1) = down(:,:,k)*dtrans(:,:,k) + b(:,:,k)*(1.0 - dtrans(:,:,k)) |
| 278 |
end do |
| 279 |
|
| 280 |
do k = 1,n+1 |
| 281 |
solar_down(:,:,k) = solar(:,:)*exp(-solar_tau(:,:,k)) |
| 282 |
end do |
| 283 |
|
| 284 |
surf_lw_down = down(:,:,n+1) |
| 285 |
net_surf_sw_down = solar_down(:,:,n+1)*(1. - albedo(:,:)) |
| 286 |
!swin = solar_down(:,:,1) |
| 287 |
|
| 288 |
!------- downward sw flux surface ------- |
| 289 |
! if ( id_swdn_sfc > 0 ) then |
| 290 |
! used = send_data ( id_swdn_sfc, net_surf_sw_down, Time_diag) |
| 291 |
! endif |
| 292 |
!------- incoming sw flux toa ------- |
| 293 |
! if ( id_swdn_toa > 0 ) then |
| 294 |
! used = send_data ( id_swdn_toa, swin, Time_diag) |
| 295 |
! endif |
| 296 |
!------- downward lw flux surface ------- |
| 297 |
! if ( id_lwdn_sfc > 0 ) then |
| 298 |
! used = send_data ( id_lwdn_sfc, surf_lw_down, Time_diag) |
| 299 |
! endif |
| 300 |
|
| 301 |
! ------------------------------------------------------------------------- |
| 302 |
!deallocate (swin) |
| 303 |
deallocate (solar_tau) |
| 304 |
deallocate (ss, solar, tau_0, solar_tau_0, p2) |
| 305 |
if ( wv_exponent .eq. 0. ) then |
| 306 |
deallocate (del_tau) |
| 307 |
else |
| 308 |
deallocate (tau_km, tau_kp) |
| 309 |
endif |
| 310 |
! ------------------------------------------------------------------------- |
| 311 |
|
| 312 |
return |
| 313 |
end subroutine radiation_down |
| 314 |
|
| 315 |
! ================================================================================== |
| 316 |
|
| 317 |
!subroutine radiation_up (is, js, Time_diag, lat, p_half, t_surf, t, tdt) |
| 318 |
subroutine radiation_up ( is, js, Time_diag, lat, p_half, t_surf, t, tdt, olr, & |
| 319 |
albedo, dtrans, b, down, solar_down, & |
| 320 |
myThid ) |
| 321 |
|
| 322 |
! Now complete the radiation calculation by computing the upward and net fluxes. |
| 323 |
|
| 324 |
integer, intent(in) :: is, js |
| 325 |
!type(time_type), intent(in) :: Time_diag |
| 326 |
real, intent(in) :: Time_diag |
| 327 |
real, intent(in) , dimension(:,:) :: lat |
| 328 |
real, intent(in) , dimension(:,:) :: t_surf |
| 329 |
real, intent(in) , dimension(:,:,:) :: t, p_half |
| 330 |
real, intent(inout), dimension(:,:,:) :: tdt |
| 331 |
real, intent(out), dimension(:,:) :: olr |
| 332 |
real, intent(in), dimension(:,:) :: albedo |
| 333 |
real, intent(in), dimension(:,:,:) :: dtrans |
| 334 |
real, intent(in), dimension(:,:,:) :: b |
| 335 |
real, intent(in), dimension(:,:,:) :: down |
| 336 |
real, intent(in), dimension(:,:,:) :: solar_down |
| 337 |
integer, intent(in) :: myThid |
| 338 |
|
| 339 |
integer :: i, j, k, n |
| 340 |
integer :: im, jm |
| 341 |
|
| 342 |
logical :: used |
| 343 |
|
| 344 |
! ------------------------------------------------------------------------- |
| 345 |
real, allocatable, dimension(:,:) :: b_surf |
| 346 |
real, allocatable, dimension(:,:,:) :: tdt_rad, entrop_rad, tdt_sw |
| 347 |
real, allocatable, dimension(:,:,:) :: up, net, flux_rad, flux_sw |
| 348 |
! ------------------------------------------------------------------------- |
| 349 |
|
| 350 |
n = size(t,3) |
| 351 |
im = size(t,1) |
| 352 |
jm = size(t,2) |
| 353 |
|
| 354 |
! ------------------------------------------------------------------------- |
| 355 |
allocate (tdt_rad (im, jm, n)) |
| 356 |
allocate (tdt_sw (im, jm, n)) |
| 357 |
allocate (entrop_rad (im, jm, n)) |
| 358 |
allocate (up (im, jm, n+1)) |
| 359 |
allocate (net (im, jm, n+1)) |
| 360 |
allocate (flux_rad (im, jm, n+1)) |
| 361 |
allocate (flux_sw (im, jm, n+1)) |
| 362 |
allocate (b_surf (im, jm)) |
| 363 |
! ------------------------------------------------------------------------- |
| 364 |
|
| 365 |
! total flux from surface |
| 366 |
b_surf = stefan*t_surf*t_surf*t_surf*t_surf |
| 367 |
|
| 368 |
! first deal with non-window upward flux |
| 369 |
up(:,:,n+1) = b_surf*(1.0-window) |
| 370 |
do k = n,1,-1 |
| 371 |
up(:,:,k) = up(:,:,k+1)*dtrans(:,:,k) + b(:,:,k)*(1.0 - dtrans(:,:,k)) |
| 372 |
end do |
| 373 |
|
| 374 |
! add upward flux in spectral window |
| 375 |
do k = 1,n+1 |
| 376 |
up(:,:,k) = up(:,:,k) + b_surf(:,:)*window |
| 377 |
end do |
| 378 |
|
| 379 |
do k = 1,n+1 |
| 380 |
net(:,:,k) = up(:,:,k)-down(:,:,k) |
| 381 |
flux_sw(:,:,k) = albedo(:,:)*solar_down(:,:,n+1) - solar_down(:,:,k) |
| 382 |
flux_rad(:,:,k) = net(:,:,k) + flux_sw(:,:,k) |
| 383 |
end do |
| 384 |
|
| 385 |
do k = 1,n |
| 386 |
tdt_rad(:,:,k) = (net(:,:,k+1) - net(:,:,k) - solar_down(:,:,k+1) + solar_down(:,:,k)) & |
| 387 |
*grav/(cp_air*(p_half(:,:,k+1)-p_half(:,:,k))) |
| 388 |
tdt_sw(:,:,k) = (- solar_down(:,:,k+1) + solar_down(:,:,k)) & |
| 389 |
*grav/(cp_air*(p_half(:,:,k+1)-p_half(:,:,k))) |
| 390 |
tdt(:,:,k) = tdt(:,:,k) + tdt_rad(:,:,k) |
| 391 |
end do |
| 392 |
|
| 393 |
olr = up(:,:,1) |
| 394 |
|
| 395 |
!------- outgoing lw flux toa (olr) ------- |
| 396 |
! if ( id_olr > 0 ) then |
| 397 |
! used = send_data ( id_olr, olr, Time_diag) |
| 398 |
! endif |
| 399 |
!------- upward lw flux surface ------- |
| 400 |
! if ( id_lwup_sfc > 0 ) then |
| 401 |
! used = send_data ( id_lwup_sfc, b_surf, Time_diag) |
| 402 |
! endif |
| 403 |
!------- temperature tendency due to radiation ------------ |
| 404 |
! if ( id_tdt_rad > 0 ) then |
| 405 |
! used = send_data ( id_tdt_rad, tdt_rad, Time_diag) |
| 406 |
! endif |
| 407 |
! if ( id_tdt_sw > 0 ) then |
| 408 |
! used = send_data ( id_tdt_sw, tdt_sw, Time_diag) |
| 409 |
! endif |
| 410 |
!------- total radiative flux (at half levels) ----------- |
| 411 |
! if ( id_flux_rad > 0 ) then |
| 412 |
! used = send_data ( id_flux_rad, flux_rad, Time_diag) |
| 413 |
! endif |
| 414 |
!------- longwave radiative flux (at half levels) -------- |
| 415 |
! if ( id_flux_lw > 0 ) then |
| 416 |
! used = send_data ( id_flux_lw, net, Time_diag) |
| 417 |
! endif |
| 418 |
! if ( id_flux_sw > 0 ) then |
| 419 |
! used = send_data ( id_flux_sw, flux_sw, Time_diag) |
| 420 |
! endif |
| 421 |
! if ( id_entrop_rad > 0 ) then |
| 422 |
! do k=1,n |
| 423 |
! entrop_rad(:,:,k) =tdt_rad(:,:,k)/t(:,:,k)*p_half(:,:,n+1)/1.e5 |
| 424 |
! end do |
| 425 |
! used = send_data ( id_entrop_rad, entrop_rad, Time_diag) |
| 426 |
! endif |
| 427 |
|
| 428 |
! ------------------------------------------------------------------------- |
| 429 |
deallocate (tdt_rad, tdt_sw, entrop_rad) |
| 430 |
deallocate (up, net, flux_rad, flux_sw) |
| 431 |
deallocate (b_surf) |
| 432 |
! ------------------------------------------------------------------------- |
| 433 |
|
| 434 |
return |
| 435 |
end subroutine radiation_up |
| 436 |
|
| 437 |
! ================================================================================== |
| 438 |
|
| 439 |
subroutine radiation_end |
| 440 |
|
| 441 |
!deallocate (b, tdt_rad, tdt_sw, entrop_rad) |
| 442 |
!!deallocate (up, down, net, solar_down, flux_rad, flux_sw) |
| 443 |
!deallocate (up, net, flux_rad, flux_sw) |
| 444 |
!deallocate (b_surf, olr, swin, albedo) |
| 445 |
!!deallocate (dtrans) |
| 446 |
!deallocate (tau, solar_tau) |
| 447 |
!deallocate (ss, solar, tau_0, solar_tau_0, p2) |
| 448 |
|
| 449 |
end subroutine radiation_end |
| 450 |
|
| 451 |
! ================================================================================== |
| 452 |
|
| 453 |
end module radiation_mod |