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

Annotation 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 - (hide 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 jmc 1.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 jmc 1.4 character(len=128) :: version='$Id: radiation_mod.F90,v 1.3 2012/10/16 15:21:39 jmc Exp $'
25 jmc 1.1 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 jmc 1.3 real :: linear_tau = 0.1
44     real :: ir_tau_co2 = 0.8678
45     real :: ir_tau_wv = 1.9979e+3
46 jmc 1.1 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 jmc 1.3 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 jmc 1.1 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 jmc 1.3 subroutine radiation_down (is, js, Time_diag, lat, p_half, t, q, &
175 jmc 1.1 ! 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 jmc 1.3 real, intent(in) , dimension(:,:,:) :: t, q, p_half
190 jmc 1.1 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 jmc 1.3 !real, allocatable, dimension(:,:) :: swin
204 jmc 1.1 real, allocatable, dimension(:,:) :: ss, solar, tau_0, solar_tau_0, p2
205 jmc 1.3 real, allocatable, dimension(:,:,:) :: solar_tau
206     real, allocatable, dimension(:,:) :: del_tau, tau_km, tau_kp
207 jmc 1.1 ! -------------------------------------------------------------------------
208    
209     n = size(t,3)
210     im = size(t,1)
211     jm = size(t,2)
212    
213     ! -------------------------------------------------------------------------
214 jmc 1.3 !allocate (swin (im, jm))
215 jmc 1.1 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 jmc 1.3 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 jmc 1.1 ! -------------------------------------------------------------------------
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 jmc 1.3 end do
245 jmc 1.1
246 jmc 1.3 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 jmc 1.1
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 jmc 1.3 !swin = solar_down(:,:,1)
287 jmc 1.1
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 jmc 1.3 !deallocate (swin)
303     deallocate (solar_tau)
304 jmc 1.1 deallocate (ss, solar, tau_0, solar_tau_0, p2)
305 jmc 1.3 if ( wv_exponent .eq. 0. ) then
306     deallocate (del_tau)
307     else
308     deallocate (tau_km, tau_kp)
309     endif
310 jmc 1.1 ! -------------------------------------------------------------------------
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 jmc 1.2 subroutine radiation_up ( is, js, Time_diag, lat, p_half, t_surf, t, tdt, olr, &
319     albedo, dtrans, b, down, solar_down, &
320 jmc 1.1 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 jmc 1.2 real, intent(out), dimension(:,:) :: olr
332 jmc 1.1 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 jmc 1.2 deallocate (up, net, flux_rad, flux_sw)
431 jmc 1.1 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