/[MITgcm]/MITgcm_contrib/verification_other/atm_gray/code_gray/radiation_mod.F90
ViewVC logotype

Contents of /MITgcm_contrib/verification_other/atm_gray/code_gray/radiation_mod.F90

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.4 - (show annotations) (download)
Wed May 8 23:55:45 2013 UTC (13 years, 6 months ago) by jmc
Branch: MAIN
CVS Tags: HEAD
Changes since 1.3: +1 -1 lines
FILE REMOVED
updated to use new pkg "atm_phys" (instead of "mypackage" + code in local
  dir "code_gray")

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

  ViewVC Help
Powered by ViewVC 1.1.22