Line data Source code
1 : #ifdef AIX
2 : #define USE_ESSL
3 : #endif
4 : #define USE_BDE
5 :
6 : module mo_jlong
7 :
8 : use shr_kind_mod, only : r4 => shr_kind_r4
9 : use shr_kind_mod, only : r8 => shr_kind_r8
10 : use cam_logfile, only : iulog
11 : use cam_abortutils, only : endrun
12 : #ifdef SPMD
13 : use mpishorthand, only : mpicom,mpiint,mpir8, mpilog, mpir4
14 : #endif
15 : use spmd_utils, only : masterproc
16 :
17 : implicit none
18 :
19 : interface jlong
20 : module procedure jlong_photo
21 : module procedure jlong_hrates
22 : end interface
23 :
24 : private
25 : public :: jlong_init
26 : public :: jlong_timestep_init
27 : public :: jlong
28 : public :: numj
29 :
30 : save
31 :
32 : real(r8), parameter :: hc = 6.62608e-34_r8 * 2.9979e8_r8 / 1.e-9_r8
33 : real(r8), parameter :: wc_o2_b = 242.37_r8 ! (nm)
34 : real(r8), parameter :: wc_o3_a = 310.32_r8 ! (nm)
35 : real(r8), parameter :: wc_o3_b = 1179.87_r8 ! (nm)
36 :
37 : integer :: nw ! wavelengths >200nm
38 : integer :: nt ! number of temperatures in xsection table
39 : integer :: np_xs ! number of pressure levels in xsection table
40 : integer :: numj ! number of photorates in xsqy, rsf
41 : integer :: nump ! number of altitudes in rsf
42 : integer :: numsza ! number of zen angles in rsf
43 : integer :: numalb ! number of albedos in rsf
44 : integer :: numcolo3 ! number of o3 columns in rsf
45 : real(r4), allocatable :: xsqy(:,:,:,:)
46 : real(r8), allocatable :: wc(:)
47 : real(r8), allocatable :: we(:)
48 : real(r8), allocatable :: wlintv(:)
49 : real(r8), allocatable :: etfphot(:)
50 : real(r8), allocatable :: bde_o2_b(:)
51 : real(r8), allocatable :: bde_o3_a(:)
52 : real(r8), allocatable :: bde_o3_b(:)
53 : real(r8), allocatable :: xs_o2b(:,:,:)
54 : real(r8), allocatable :: xs_o3a(:,:,:)
55 : real(r8), allocatable :: xs_o3b(:,:,:)
56 : real(r8), allocatable :: p(:)
57 : real(r8), allocatable :: del_p(:)
58 : real(r8), allocatable :: prs(:)
59 : real(r8), allocatable :: dprs(:)
60 : real(r8), allocatable :: sza(:)
61 : real(r8), allocatable :: del_sza(:)
62 : real(r8), allocatable :: alb(:)
63 : real(r8), allocatable :: del_alb(:)
64 : real(r8), allocatable :: o3rat(:)
65 : real(r8), allocatable :: del_o3rat(:)
66 : real(r8), allocatable :: colo3(:)
67 : real(r4), allocatable :: rsf_tab(:,:,:,:,:)
68 : logical :: jlong_used = .false.
69 :
70 : contains
71 :
72 0 : subroutine jlong_init( xs_long_file, rsf_file, lng_indexer )
73 :
74 : use ppgrid, only : pver
75 : use mo_util, only : rebin
76 : use solar_irrad_data,only : data_nw => nbins, data_we => we, data_etf => sol_etf
77 :
78 : implicit none
79 :
80 : !------------------------------------------------------------------------------
81 : ! ... dummy arguments
82 : !------------------------------------------------------------------------------
83 : integer, intent(inout) :: lng_indexer(:)
84 : character(len=*), intent(in) :: xs_long_file, rsf_file
85 :
86 : !------------------------------------------------------------------------------
87 : ! ... read Cross Section * QY NetCDF file
88 : ! find temperature index for given altitude
89 : ! derive cross*QY results, returns xsqy(nj,nz,nw)
90 : !------------------------------------------------------------------------------
91 0 : call get_xsqy( xs_long_file, lng_indexer )
92 :
93 : !------------------------------------------------------------------------------
94 : ! ... read radiative source function NetCDF file
95 : !------------------------------------------------------------------------------
96 0 : if(masterproc) write(iulog,*) 'jlong_init: before get_rsf'
97 0 : call get_rsf(rsf_file)
98 0 : if(masterproc) write(iulog,*) 'jlong_init: after get_rsf'
99 :
100 0 : we(:nw) = wc(:nw) - .5_r8*wlintv(:nw)
101 0 : we(nw+1) = wc(nw) + .5_r8*wlintv(nw)
102 :
103 0 : if (masterproc) then
104 0 : write(iulog,*) ' '
105 0 : write(iulog,*) '--------------------------------------------------'
106 : endif
107 0 : call rebin( data_nw, nw, data_we, we, data_etf, etfphot )
108 0 : if (masterproc) then
109 0 : write(iulog,*) 'jlong_init: etfphot after data rebin'
110 0 : write(iulog,'(1p,5g15.7)') etfphot(:)
111 0 : write(iulog,*) '--------------------------------------------------'
112 0 : write(iulog,*) ' '
113 : endif
114 :
115 0 : jlong_used = .true.
116 :
117 0 : end subroutine jlong_init
118 :
119 0 : subroutine get_xsqy( xs_long_file, lng_indexer )
120 : !=============================================================================!
121 : ! PURPOSE: !
122 : ! Reads a NetCDF file that contains: !
123 : ! cross section * QY temperature dependence, >200nm !
124 : ! !
125 : !=============================================================================!
126 : ! PARAMETERS: !
127 : ! Input: !
128 : ! filepath.... NetCDF filepath that contains the "cross sections" !
129 : !=============================================================================!
130 : ! EDIT HISTORY: !
131 : ! Created by Doug Kinnison, 3/14/2002 !
132 : !=============================================================================!
133 :
134 0 : use ioFileMod, only : getfil
135 : use error_messages, only : alloc_err
136 : use chem_mods, only : phtcnt, pht_alias_lst, rxt_tag_lst
137 : use netcdf, only: &
138 : nf90_open, &
139 : nf90_nowrite, &
140 : nf90_inq_dimid, &
141 : nf90_inquire_dimension, &
142 : nf90_inq_varid, &
143 : nf90_noerr, &
144 : nf90_get_var, &
145 : nf90_close
146 :
147 : !------------------------------------------------------------------------------
148 : ! ... Dummy arguments
149 : !------------------------------------------------------------------------------
150 : integer, intent(inout) :: lng_indexer(:)
151 : character(len=*) :: xs_long_file
152 :
153 : !------------------------------------------------------------------------------
154 : ! ... Local variables
155 : !------------------------------------------------------------------------------
156 : integer :: varid, dimid, ndx
157 : integer :: ncid
158 : integer :: iret
159 : integer :: i, k, m, n
160 : integer :: wrk_ndx(phtcnt)
161 : character(len=256) :: locfn
162 :
163 0 : Masterproc_only : if( masterproc ) then
164 : !------------------------------------------------------------------------------
165 : ! ... open NetCDF File
166 : !------------------------------------------------------------------------------
167 0 : call getfil(xs_long_file, locfn, 0)
168 0 : iret = nf90_open(trim(locfn), NF90_NOWRITE, ncid)
169 :
170 : !------------------------------------------------------------------------------
171 : ! ... get dimensions
172 : !------------------------------------------------------------------------------
173 0 : iret = nf90_inq_dimid( ncid, 'numtemp', dimid )
174 0 : iret = nf90_inquire_dimension( ncid, dimid,len= nt )
175 0 : iret = nf90_inq_dimid( ncid, 'numwl', dimid )
176 0 : iret = nf90_inquire_dimension( ncid, dimid,len= nw )
177 0 : iret = nf90_inq_dimid( ncid, 'numprs', dimid )
178 0 : iret = nf90_inquire_dimension( ncid, dimid,len= np_xs )
179 : !------------------------------------------------------------------------------
180 : ! ... check for cross section in dataset
181 : !------------------------------------------------------------------------------
182 : do m = 1,phtcnt
183 : if( pht_alias_lst(m,2) == ' ' ) then
184 : iret = nf90_inq_varid( ncid, rxt_tag_lst(m), varid )
185 : if( iret == nf90_noerr ) then
186 : lng_indexer(m) = varid
187 : end if
188 : else if( pht_alias_lst(m,2) == 'userdefined' ) then
189 : lng_indexer(m) = -1
190 : else
191 : iret = nf90_inq_varid( ncid, pht_alias_lst(m,2), varid )
192 : if( iret == nf90_noerr ) then
193 : lng_indexer(m) = varid
194 : else
195 : write(iulog,*) 'get_xsqy : ',rxt_tag_lst(m)(:len_trim(rxt_tag_lst(m))),' alias ', &
196 : pht_alias_lst(m,2)(:len_trim(pht_alias_lst(m,2))),' not in dataset'
197 : call endrun
198 : end if
199 : end if
200 : end do
201 0 : numj = 0
202 : do m = 1,phtcnt
203 : if( lng_indexer(m) > 0 ) then
204 : if( any( lng_indexer(:m-1) == lng_indexer(m) ) ) then
205 : cycle
206 : end if
207 : numj = numj + 1
208 : end if
209 : end do
210 :
211 : !------------------------------------------------------------------------------
212 : ! ... allocate arrays
213 : !------------------------------------------------------------------------------
214 :
215 0 : allocate( xsqy(numj,nw,nt,np_xs),stat=iret )
216 0 : if( iret /= 0 ) then
217 0 : call alloc_err( iret, 'get_xsqy', 'xsqy', numj*nt*np_xs*nw )
218 : end if
219 0 : allocate( prs(np_xs),dprs(np_xs-1),stat=iret )
220 0 : if( iret /= 0 ) then
221 0 : call alloc_err( iret, 'get_xsqy', 'prs,dprs', np_xs )
222 : end if
223 0 : allocate( xs_o2b(nw,nt,np_xs),xs_o3a(nw,nt,np_xs),xs_o3b(nw,nt,np_xs),stat=iret )
224 0 : if( iret /= 0 ) then
225 0 : call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs )
226 : end if
227 : !------------------------------------------------------------------------------
228 : ! ... read cross sections
229 : !------------------------------------------------------------------------------
230 : ndx = 0
231 : do m = 1,phtcnt
232 : if( lng_indexer(m) > 0 ) then
233 : if( any( lng_indexer(:m-1) == lng_indexer(m) ) ) then
234 : cycle
235 : end if
236 : ndx = ndx + 1
237 : iret = nf90_get_var( ncid, lng_indexer(m), xsqy(ndx,:,:,:) )
238 : end if
239 : end do
240 0 : if( ndx /= numj ) then
241 0 : write(iulog,*) 'get_xsqy : ndx count /= cross section count'
242 0 : call endrun
243 : end if
244 0 : iret = nf90_inq_varid( ncid, 'jo2_b', varid )
245 0 : iret = nf90_get_var( ncid, varid, xs_o2b )
246 0 : iret = nf90_inq_varid( ncid, 'jo3_a', varid )
247 0 : iret = nf90_get_var( ncid, varid, xs_o3a )
248 0 : iret = nf90_inq_varid( ncid, 'jo3_b', varid )
249 0 : iret = nf90_get_var( ncid, varid, xs_o3b )
250 : !------------------------------------------------------------------------------
251 : ! ... setup final lng_indexer
252 : !------------------------------------------------------------------------------
253 0 : ndx = 0
254 0 : wrk_ndx(:) = lng_indexer(:)
255 : do m = 1,phtcnt
256 : if( wrk_ndx(m) > 0 ) then
257 : ndx = ndx + 1
258 : i = wrk_ndx(m)
259 : where( wrk_ndx(:) == i )
260 : lng_indexer(:) = ndx
261 : wrk_ndx(:) = -100000
262 : end where
263 : end if
264 : end do
265 :
266 0 : iret = nf90_inq_varid( ncid, 'pressure', varid )
267 0 : iret = nf90_get_var( ncid, varid, prs )
268 0 : iret = nf90_close( ncid )
269 : end if Masterproc_only
270 :
271 : #ifdef SPMD
272 : ! call mpibarrier( mpicom )
273 0 : call mpibcast( numj, 1, mpiint, 0, mpicom )
274 0 : call mpibcast( nt, 1, mpiint, 0, mpicom )
275 0 : call mpibcast( nw, 1, mpiint, 0, mpicom )
276 0 : call mpibcast( np_xs, 1, mpiint, 0, mpicom )
277 0 : call mpibcast( lng_indexer, phtcnt, mpiint, 0, mpicom )
278 : #endif
279 0 : if( .not. masterproc ) then
280 : !------------------------------------------------------------------------------
281 : ! ... allocate arrays
282 : !------------------------------------------------------------------------------
283 0 : allocate( xsqy(numj,nw,nt,np_xs),stat=iret )
284 0 : if( iret /= nf90_noerr) then
285 0 : write(iulog,*) 'get_xsqy : failed to allocate xsqy ; error = ',iret
286 0 : call endrun
287 : end if
288 0 : allocate( prs(np_xs),dprs(np_xs-1),stat=iret )
289 0 : if( iret /= nf90_noerr) then
290 0 : write(iulog,*) 'get_xsqy : failed to allocate prs,dprs ; error = ',iret
291 0 : call endrun
292 : end if
293 0 : allocate( xs_o2b(nw,nt,np_xs),xs_o3a(nw,nt,np_xs),xs_o3b(nw,nt,np_xs),stat=iret )
294 0 : if( iret /= 0 ) then
295 0 : call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs )
296 : end if
297 : end if
298 : #ifdef SPMD
299 : ! call mpibarrier( mpicom )
300 0 : call mpibcast( prs, np_xs, mpir8, 0, mpicom )
301 0 : call mpibcast( xsqy, numj*nt*np_xs*nw, mpir4, 0, mpicom )
302 0 : call mpibcast( xs_o2b, nw*nt*np_xs, mpir8, 0, mpicom )
303 0 : call mpibcast( xs_o3a, nw*nt*np_xs, mpir8, 0, mpicom )
304 0 : call mpibcast( xs_o3b, nw*nt*np_xs, mpir8, 0, mpicom )
305 : #endif
306 0 : dprs(:np_xs-1) = 1._r8/(prs(1:np_xs-1) - prs(2:np_xs))
307 :
308 0 : end subroutine get_xsqy
309 :
310 0 : subroutine get_rsf(rsf_file)
311 : !=============================================================================!
312 : ! PURPOSE: !
313 : ! Reads a NetCDF file that contains:
314 : ! Radiative Souce function !
315 : !=============================================================================!
316 : ! PARAMETERS: !
317 : ! Input: !
318 : ! filepath.... NetCDF file that contains the RSF !
319 : ! !
320 : ! Output: !
321 : ! rsf ........ Radiative Source Function (quanta cm-2 sec-1 !
322 : ! !
323 : ! EDIT HISTORY: !
324 : ! Created by Doug Kinnison, 3/14/2002 !
325 : !=============================================================================!
326 :
327 : use ioFileMod, only : getfil
328 : use error_messages, only : alloc_err
329 : use netcdf, only: &
330 : nf90_open, &
331 : nf90_nowrite, &
332 : nf90_inq_dimid, &
333 : nf90_inquire_dimension, &
334 : nf90_inq_varid, &
335 : nf90_noerr, &
336 : nf90_get_var, &
337 : nf90_close
338 :
339 : !------------------------------------------------------------------------------
340 : ! ... dummy arguments
341 : !------------------------------------------------------------------------------
342 : character(len=*) :: rsf_file
343 :
344 : !------------------------------------------------------------------------------
345 : ! ... local variables
346 : !------------------------------------------------------------------------------
347 : integer :: varid, dimid
348 : integer :: ncid
349 : integer :: i, j, k, l, w
350 : integer :: iret
351 : integer :: count(5)
352 : integer :: start(5)
353 : real(r8) :: wrk
354 : character(len=256) :: locfn
355 :
356 0 : Masterproc_only : if( masterproc ) then
357 : !------------------------------------------------------------------------------
358 : ! ... open NetCDF File
359 : !------------------------------------------------------------------------------
360 0 : call getfil(rsf_file, locfn, 0)
361 0 : iret = nf90_open(trim(locfn), NF90_NOWRITE, ncid)
362 :
363 : !------------------------------------------------------------------------------
364 : ! ... get dimensions
365 : !------------------------------------------------------------------------------
366 0 : iret = nf90_inq_dimid( ncid, 'numz', dimid )
367 0 : iret = nf90_inquire_dimension( ncid, dimid,len= nump )
368 0 : iret = nf90_inq_dimid( ncid, 'numsza', dimid )
369 0 : iret = nf90_inquire_dimension( ncid, dimid,len= numsza )
370 0 : iret = nf90_inq_dimid( ncid, 'numalb', dimid )
371 0 : iret = nf90_inquire_dimension( ncid, dimid,len= numalb )
372 0 : iret = nf90_inq_dimid( ncid, 'numcolo3fact', dimid )
373 0 : iret = nf90_inquire_dimension( ncid, dimid,len= numcolo3 )
374 :
375 : end if Masterproc_only
376 : #ifdef SPMD
377 : ! call mpibarrier( mpicom )
378 0 : call mpibcast( nump, 1, mpiint, 0, mpicom )
379 0 : call mpibcast( numsza, 1, mpiint, 0, mpicom )
380 0 : call mpibcast( numalb, 1, mpiint, 0, mpicom )
381 0 : call mpibcast( numcolo3, 1, mpiint, 0, mpicom )
382 : #endif
383 : !------------------------------------------------------------------------------
384 : ! ... allocate arrays
385 : !------------------------------------------------------------------------------
386 0 : allocate( wc(nw),stat=iret )
387 0 : if( iret /= 0 ) then
388 0 : call alloc_err( iret, 'get_rsf', 'wc', nw )
389 : end if
390 0 : allocate( wlintv(nw),we(nw+1),etfphot(nw),stat=iret )
391 0 : if( iret /= 0 ) then
392 0 : call alloc_err( iret, 'get_rsf', 'wlintv,etfphot', nw )
393 : end if
394 0 : allocate( bde_o2_b(nw),bde_o3_a(nw),bde_o3_b(nw),stat=iret )
395 0 : if( iret /= 0 ) then
396 0 : call alloc_err( iret, 'get_rsf', 'bde', nw )
397 : end if
398 0 : allocate( p(nump),del_p(nump-1),stat=iret )
399 0 : if( iret /= 0 ) then
400 0 : call alloc_err( iret, 'get_rsf', 'p,delp', nump )
401 : end if
402 0 : allocate( sza(numsza),del_sza(numsza-1),stat=iret )
403 0 : if( iret /= 0 ) then
404 0 : call alloc_err( iret, 'get_rsf', 'sza,del_sza', numsza )
405 : end if
406 0 : allocate( alb(numalb),del_alb(numalb-1),stat=iret )
407 0 : if( iret /= 0 ) then
408 0 : call alloc_err( iret, 'get_rsf', 'alb,del_alb', numalb )
409 : end if
410 0 : allocate( o3rat(numcolo3),del_o3rat(numcolo3-1),stat=iret )
411 0 : if( iret /= 0 ) then
412 0 : call alloc_err( iret, 'get_rsf', 'o3rat,del_o3rat', numcolo3 )
413 : end if
414 0 : allocate( colo3(nump),stat=iret )
415 0 : if( iret /= 0 ) then
416 0 : call alloc_err( iret, 'get_rsf', 'colo3', nump )
417 : end if
418 0 : allocate( rsf_tab(nw,nump,numsza,numcolo3,numalb),stat=iret )
419 0 : if( iret /= 0 ) then
420 0 : write(iulog,*) 'get_rsf : dimensions = ',nw,nump,numsza,numcolo3,numalb
421 0 : call alloc_err( iret, 'get_rsf', 'rsf_tab', numalb*numcolo3*numsza*nump )
422 : end if
423 :
424 0 : Masterproc_only2 : if( masterproc ) then
425 : !------------------------------------------------------------------------------
426 : ! ... read variables
427 : !------------------------------------------------------------------------------
428 0 : iret = nf90_inq_varid( ncid, 'wc', varid )
429 0 : iret = nf90_get_var( ncid, varid, wc )
430 0 : iret = nf90_inq_varid( ncid, 'wlintv', varid )
431 0 : iret = nf90_get_var( ncid, varid, wlintv )
432 0 : iret = nf90_inq_varid( ncid, 'pm', varid )
433 0 : iret = nf90_get_var( ncid, varid, p )
434 0 : iret = nf90_inq_varid( ncid, 'sza', varid )
435 0 : iret = nf90_get_var( ncid, varid, sza )
436 0 : iret = nf90_inq_varid( ncid, 'alb', varid )
437 0 : iret = nf90_get_var( ncid, varid, alb )
438 0 : iret = nf90_inq_varid( ncid, 'colo3fact', varid )
439 0 : iret = nf90_get_var( ncid, varid, o3rat )
440 0 : iret = nf90_inq_varid( ncid, 'colo3', varid )
441 0 : iret = nf90_get_var( ncid, varid, colo3 )
442 :
443 0 : iret = nf90_inq_varid( ncid, 'RSF', varid )
444 :
445 0 : if (masterproc) then
446 0 : write(iulog,*) ' '
447 0 : write(iulog,*) '----------------------------------------------'
448 0 : write(iulog,*) 'get_rsf: numalb, numcolo3, numsza, nump = ', &
449 0 : numalb, numcolo3, numsza, nump
450 0 : write(iulog,*) 'get_rsf: size of rsf_tab = ',size(rsf_tab,dim=1),size(rsf_tab,dim=2), &
451 0 : size(rsf_tab,dim=3),size(rsf_tab,dim=4),size(rsf_tab,dim=5)
452 0 : write(iulog,*) '----------------------------------------------'
453 0 : write(iulog,*) ' '
454 : endif
455 :
456 0 : iret = nf90_get_var( ncid, varid, rsf_tab )
457 0 : iret = nf90_close( ncid )
458 :
459 0 : do w = 1,nw
460 0 : wrk = wlintv(w)
461 0 : rsf_tab(w,:,:,:,:) = wrk*rsf_tab(w,:,:,:,:)
462 : enddo
463 :
464 : end if Masterproc_only2
465 : #ifdef SPMD
466 0 : call mpibcast( wc, nw, mpir8, 0, mpicom )
467 0 : call mpibcast( wlintv, nw, mpir8, 0, mpicom )
468 0 : call mpibcast( p, nump, mpir8, 0, mpicom )
469 0 : call mpibcast( sza, numsza, mpir8, 0, mpicom )
470 0 : call mpibcast( alb, numalb, mpir8, 0, mpicom )
471 0 : call mpibcast( o3rat, numcolo3, mpir8, 0, mpicom )
472 0 : call mpibcast( colo3, nump, mpir8, 0, mpicom )
473 0 : do w = 1,nw
474 0 : call mpibcast( rsf_tab(w,:,:,:,:), numalb*numcolo3*numsza*nump, mpir4, 0, mpicom )
475 : enddo
476 : #endif
477 : #ifdef USE_BDE
478 0 : if (masterproc) write(iulog,*) 'Jlong using bdes'
479 0 : bde_o2_b(:) = max( 0._r8, hc*(wc_o2_b - wc(:))/(wc_o2_b*wc(:)) )
480 0 : bde_o3_a(:) = max( 0._r8, hc*(wc_o3_a - wc(:))/(wc_o3_a*wc(:)) )
481 0 : bde_o3_b(:) = max( 0._r8, hc*(wc_o3_b - wc(:))/(wc_o3_b*wc(:)) )
482 : #else
483 : if (masterproc) write(iulog,*) 'Jlong not using bdes'
484 : bde_o2_b(:) = hc/wc(:)
485 : bde_o3_a(:) = hc/wc(:)
486 : bde_o3_b(:) = hc/wc(:)
487 : #endif
488 :
489 0 : del_p(:nump-1) = 1._r8/abs(p(1:nump-1)- p(2:nump))
490 0 : del_sza(:numsza-1) = 1._r8/(sza(2:numsza) - sza(1:numsza-1))
491 0 : del_alb(:numalb-1) = 1._r8/(alb(2:numalb) - alb(1:numalb-1))
492 0 : del_o3rat(:numcolo3-1) = 1._r8/(o3rat(2:numcolo3) - o3rat(1:numcolo3-1))
493 :
494 0 : end subroutine get_rsf
495 :
496 0 : subroutine jlong_timestep_init
497 : !---------------------------------------------------------------
498 : ! ... set etfphot if required
499 : !---------------------------------------------------------------
500 :
501 : use mo_util, only : rebin
502 :
503 : use solar_irrad_data,only : data_nw => nbins, data_we => we, data_etf => sol_etf
504 :
505 : implicit none
506 :
507 0 : if (.not.jlong_used) return
508 :
509 0 : call rebin( data_nw, nw, data_we, we, data_etf, etfphot )
510 :
511 0 : end subroutine jlong_timestep_init
512 :
513 0 : subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, &
514 0 : mw, o2_vmr, o3_vmr, colo3_in, qrl_col, &
515 0 : cparg, kbot )
516 : !==============================================================================
517 : ! Purpose:
518 : ! To calculate the thermal heating rates longward of 200nm.
519 : !==============================================================================
520 : ! Approach:
521 : ! 1) Reads the Cross Section*QY NetCDF file
522 : ! 2) Given a temperature profile, derives the appropriate XS*QY
523 : !
524 : ! 3) Reads the Radiative Source function (RSF) NetCDF file
525 : ! Units = quanta cm-2 sec-1
526 : !
527 : ! 4) Indices are supplied to select a RSF that is consistent with
528 : ! the reference atmosphere in TUV (for direct comparision of J's).
529 : ! This approach will be replaced in the global model. Here colo3, zenith
530 : ! angle, and altitude will be inputed and the correct entry in the table
531 : ! will be derived.
532 : !==============================================================================
533 :
534 0 : use physconst, only : avogad
535 : use error_messages, only : alloc_err
536 :
537 : implicit none
538 :
539 : !------------------------------------------------------------------------------
540 : ! ... dummy arguments
541 : !------------------------------------------------------------------------------
542 : integer, intent (in) :: nlev ! number vertical levels
543 : integer, intent (in) :: kbot ! heating levels
544 : real(r8), intent(in) :: o2_vmr(nlev) ! o2 conc (mol/mol)
545 : real(r8), intent(in) :: o3_vmr(nlev) ! o3 conc (mol/mol)
546 : real(r8), intent(in) :: sza_in ! solar zenith angle (degrees)
547 : real(r8), intent(in) :: alb_in(nlev) ! albedo
548 : real(r8), intent(in) :: p_in(nlev) ! midpoint pressure (hPa)
549 : real(r8), intent(in) :: t_in(nlev) ! Temperature profile (K)
550 : real(r8), intent(in) :: colo3_in(nlev) ! o3 column density (molecules/cm^3)
551 : real(r8), intent(in) :: mw(nlev) ! atms molecular weight
552 : real(r8), intent(in) :: cparg(nlev) ! specific heat capacity
553 : real(r8), intent(inout) :: qrl_col(:,:) ! heating rates
554 :
555 : !----------------------------------------------------------------------
556 : ! ... local variables
557 : !----------------------------------------------------------------------
558 : integer :: astat
559 : integer :: k, km, m
560 : integer :: t_index ! Temperature index
561 : integer :: pndx
562 : real(r8) :: ptarget
563 : real(r8) :: delp
564 : real(r8) :: hfactor
565 0 : real(r8), allocatable :: rsf(:,:) ! Radiative source function
566 0 : real(r8), allocatable :: xswk(:,:) ! working xsection array
567 0 : real(r8), allocatable :: wrk(:) ! work vector
568 :
569 : !----------------------------------------------------------------------
570 : ! ... allocate variables
571 : !----------------------------------------------------------------------
572 0 : allocate( rsf(nw,kbot),stat=astat )
573 0 : if( astat /= 0 ) then
574 0 : call alloc_err( astat, 'jlong_hrates', 'rsf', nw*nlev )
575 : end if
576 0 : allocate( xswk(nw,3),wrk(nw),stat=astat )
577 0 : if( astat /= 0 ) then
578 0 : call alloc_err( astat, 'jlong_hrates', 'xswk,wrk', 3*nw )
579 : end if
580 :
581 : !----------------------------------------------------------------------
582 : ! ... interpolate table rsf to model variables
583 : !----------------------------------------------------------------------
584 0 : call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf )
585 :
586 : !------------------------------------------------------------------------------
587 : ! ... calculate thermal heating rates for wavelengths >200nm
588 : !------------------------------------------------------------------------------
589 : !------------------------------------------------------------------------------
590 : ! LLNL LUT approach to finding temperature index...
591 : ! Calculate the temperature index into the cross section
592 : ! data which lists coss sections for temperatures from
593 : ! 150 to 350 degrees K. Make sure the index is a value
594 : ! between 1 and 201.
595 : !------------------------------------------------------------------------------
596 0 : qrl_col(kbot+1:nlev,:) = 0._r8
597 : level_loop_1 : &
598 0 : do k = 1,kbot
599 : !----------------------------------------------------------------------
600 : ! ... get index into xsqy
601 : !----------------------------------------------------------------------
602 0 : t_index = t_in(k) - 148.5_r8
603 0 : t_index = min( 201,max( t_index,1) )
604 : !----------------------------------------------------------------------
605 : ! ... find pressure level
606 : !----------------------------------------------------------------------
607 0 : ptarget = p_in(k)
608 0 : if( ptarget >= prs(1) ) then
609 0 : xswk(:,1) = xs_o2b(:,t_index,1)
610 0 : xswk(:,2) = xs_o3a(:,t_index,1)
611 0 : xswk(:,3) = xs_o3b(:,t_index,1)
612 0 : else if( ptarget <= prs(np_xs) ) then
613 0 : xswk(:,1) = xs_o2b(:,t_index,np_xs)
614 0 : xswk(:,2) = xs_o3a(:,t_index,np_xs)
615 0 : xswk(:,3) = xs_o3b(:,t_index,np_xs)
616 : else
617 0 : do km = 2,np_xs
618 0 : if( ptarget >= prs(km) ) then
619 0 : pndx = km - 1
620 0 : delp = (prs(pndx) - ptarget)*dprs(pndx)
621 0 : exit
622 : end if
623 : end do
624 0 : xswk(:,1) = xs_o2b(:,t_index,pndx) &
625 0 : + delp*(xs_o2b(:,t_index,pndx+1) - xs_o2b(:,t_index,pndx))
626 0 : xswk(:,2) = xs_o3a(:,t_index,pndx) &
627 0 : + delp*(xs_o3a(:,t_index,pndx+1) - xs_o3a(:,t_index,pndx))
628 0 : xswk(:,3) = xs_o3b(:,t_index,pndx) &
629 0 : + delp*(xs_o3b(:,t_index,pndx+1) - xs_o3b(:,t_index,pndx))
630 : end if
631 0 : hfactor = avogad/(cparg(k)*mw(k))
632 0 : wrk(:) = xswk(:,1)*bde_o2_b(:)
633 0 : qrl_col(k,2) = dot_product( wrk(:),rsf(:,k) ) * o2_vmr(k) * hfactor
634 0 : wrk(:) = xswk(:,2)*bde_o3_a(:)
635 0 : qrl_col(k,3) = dot_product( wrk(:),rsf(:,k) ) * o3_vmr(k) * hfactor
636 0 : wrk(:) = xswk(:,3)*bde_o3_b(:)
637 0 : qrl_col(k,4) = dot_product( wrk(:),rsf(:,k) ) * o3_vmr(k) * hfactor
638 : end do level_loop_1
639 :
640 0 : deallocate( rsf, xswk, wrk )
641 :
642 0 : end subroutine jlong_hrates
643 :
644 0 : subroutine jlong_photo( nlev, sza_in, alb_in, p_in, t_in, &
645 0 : colo3_in, j_long )
646 : !==============================================================================
647 : ! Purpose:
648 : ! To calculate the total J for selective species longward of 200nm.
649 : !==============================================================================
650 : ! Approach:
651 : ! 1) Reads the Cross Section*QY NetCDF file
652 : ! 2) Given a temperature profile, derives the appropriate XS*QY
653 : !
654 : ! 3) Reads the Radiative Source function (RSF) NetCDF file
655 : ! Units = quanta cm-2 sec-1
656 : !
657 : ! 4) Indices are supplied to select a RSF that is consistent with
658 : ! the reference atmosphere in TUV (for direct comparision of J's).
659 : ! This approach will be replaced in the global model. Here colo3, zenith
660 : ! angle, and altitude will be inputed and the correct entry in the table
661 : ! will be derived.
662 : !==============================================================================
663 :
664 : use spmd_utils, only : masterproc
665 : use error_messages, only : alloc_err
666 :
667 : implicit none
668 :
669 : !------------------------------------------------------------------------------
670 : ! ... dummy arguments
671 : !------------------------------------------------------------------------------
672 : integer, intent (in) :: nlev ! number vertical levels
673 : real(r8), intent(in) :: sza_in ! solar zenith angle (degrees)
674 : real(r8), intent(in) :: alb_in(nlev) ! albedo
675 : real(r8), intent(in) :: p_in(nlev) ! midpoint pressure (hPa)
676 : real(r8), intent(in) :: t_in(nlev) ! Temperature profile (K)
677 : real(r8), intent(in) :: colo3_in(nlev) ! o3 column density (molecules/cm^3)
678 : real(r8), intent(out) :: j_long(:,:) ! photo rates (1/s)
679 :
680 : !----------------------------------------------------------------------
681 : ! ... local variables
682 : !----------------------------------------------------------------------
683 : integer :: astat
684 : integer :: k, km, m
685 : integer :: wn
686 : integer :: t_index ! Temperature index
687 : integer :: pndx
688 : real(r8) :: ptarget
689 : real(r8) :: delp
690 : real(r8) :: hfactor
691 0 : real(r8), allocatable :: rsf(:,:) ! Radiative source function
692 0 : real(r8), allocatable :: xswk(:,:) ! working xsection array
693 :
694 : !----------------------------------------------------------------------
695 : ! ... allocate variables
696 : !----------------------------------------------------------------------
697 0 : allocate( rsf(nw,nlev),stat=astat )
698 0 : if( astat /= 0 ) then
699 0 : call alloc_err( astat, 'jlong_photo', 'rsf', nw*nlev )
700 : end if
701 0 : allocate( xswk(numj,nw),stat=astat )
702 0 : if( astat /= 0 ) then
703 0 : call alloc_err( astat, 'jlong_photo', 'xswk', numj*nw )
704 : end if
705 :
706 : !----------------------------------------------------------------------
707 : ! ... interpolate table rsf to model variables
708 : !----------------------------------------------------------------------
709 0 : call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, nlev, rsf )
710 :
711 : !------------------------------------------------------------------------------
712 : ! ... calculate total Jlong for wavelengths >200nm
713 : !------------------------------------------------------------------------------
714 : !------------------------------------------------------------------------------
715 : ! LLNL LUT approach to finding temperature index...
716 : ! Calculate the temperature index into the cross section
717 : ! data which lists coss sections for temperatures from
718 : ! 150 to 350 degrees K. Make sure the index is a value
719 : ! between 1 and 201.
720 : !------------------------------------------------------------------------------
721 : level_loop_1 : &
722 0 : do k = 1,nlev
723 : !----------------------------------------------------------------------
724 : ! ... get index into xsqy
725 : !----------------------------------------------------------------------
726 0 : t_index = t_in(k) - 148.5_r8
727 0 : t_index = min( 201,max( t_index,1) )
728 : !----------------------------------------------------------------------
729 : ! ... find pressure level
730 : !----------------------------------------------------------------------
731 0 : ptarget = p_in(k)
732 0 : if( ptarget >= prs(1) ) then
733 0 : do wn = 1,nw
734 0 : xswk(:,wn) = xsqy(:,wn,t_index,1)
735 : end do
736 0 : else if( ptarget <= prs(np_xs) ) then
737 0 : do wn = 1,nw
738 0 : xswk(:,wn) = xsqy(:,wn,t_index,np_xs)
739 : end do
740 : else
741 0 : do km = 2,np_xs
742 0 : if( ptarget >= prs(km) ) then
743 0 : pndx = km - 1
744 0 : delp = (prs(pndx) - ptarget)*dprs(pndx)
745 0 : exit
746 : end if
747 : end do
748 0 : do wn = 1,nw
749 0 : xswk(:,wn) = xsqy(:,wn,t_index,pndx) &
750 0 : + delp*(xsqy(:,wn,t_index,pndx+1) - xsqy(:,wn,t_index,pndx))
751 : end do
752 : end if
753 : #ifdef USE_ESSL
754 : call dgemm( 'N', 'N', numj, 1, nw, &
755 : 1._r8, xswk, numj, rsf(1,k), nw, &
756 : 0._r8, j_long(1,k), numj )
757 : #else
758 0 : j_long(:,k) = matmul( xswk(:,:),rsf(:,k) )
759 : #endif
760 : end do level_loop_1
761 :
762 0 : deallocate( rsf, xswk )
763 :
764 0 : end subroutine jlong_photo
765 :
766 : !----------------------------------------------------------------------
767 : !----------------------------------------------------------------------
768 : ! ... interpolate table rsf to model variables
769 : !----------------------------------------------------------------------
770 : !----------------------------------------------------------------------
771 0 : subroutine interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf )
772 :
773 : use error_messages, only : alloc_err
774 :
775 : implicit none
776 :
777 : !------------------------------------------------------------------------------
778 : ! ... dummy arguments
779 : !------------------------------------------------------------------------------
780 : real(r8), intent(in) :: alb_in(:) ! albedo
781 : real(r8), intent(in) :: sza_in ! solar zenith angle (degrees)
782 : integer, intent(in) :: kbot ! heating levels
783 : real(r8), intent(in) :: p_in(:) ! midpoint pressure (hPa)
784 : real(r8), intent(in) :: colo3_in(:) ! o3 column density (molecules/cm^3)
785 : real(r8), intent(out) :: rsf(:,:)
786 :
787 : !----------------------------------------------------------------------
788 : ! ... local variables
789 : !----------------------------------------------------------------------
790 : integer :: astat
791 : integer :: is, iv, ial
792 : integer :: isp1, ivp1, ialp1
793 : real(r8), dimension(3) :: dels
794 : real(r8), dimension(0:1,0:1,0:1) :: wghtl, wghtu
795 : real(r8) :: psum_u
796 0 : real(r8), allocatable :: psum_l(:)
797 : real(r8) :: v3ratl, v3ratu
798 : integer :: pind, albind
799 : real(r8) :: wrk0, wrk1, wght1
800 : integer :: iz, k, m
801 : integer :: izl
802 : integer :: ratindl, ratindu
803 : integer :: wn
804 :
805 : !----------------------------------------------------------------------
806 : ! ... allocate variables
807 : !----------------------------------------------------------------------
808 0 : allocate( psum_l(nw),stat=astat )
809 0 : if( astat /= 0 ) then
810 0 : call alloc_err( astat, 'jlong_hrates', 'psum_l', nw )
811 : end if
812 :
813 : !----------------------------------------------------------------------
814 : ! ... find the zenith angle index ( same for all levels )
815 : !----------------------------------------------------------------------
816 0 : do is = 1,numsza
817 0 : if( sza(is) > sza_in ) then
818 : exit
819 : end if
820 : end do
821 0 : is = max( min( is,numsza ) - 1,1 )
822 0 : isp1 = is + 1
823 0 : dels(1) = max( 0._r8,min( 1._r8,(sza_in - sza(is)) * del_sza(is) ) )
824 0 : wrk0 = 1._r8 - dels(1)
825 :
826 0 : izl = 2
827 : Level_loop : &
828 0 : do k = kbot,1,-1
829 : !----------------------------------------------------------------------
830 : ! ... find albedo indicies
831 : !----------------------------------------------------------------------
832 0 : do ial = 1,numalb
833 0 : if( alb(ial) > alb_in(k) ) then
834 : exit
835 : end if
836 : end do
837 0 : albind = max( min( ial,numalb ) - 1,1 )
838 : !----------------------------------------------------------------------
839 : ! ... find pressure level indicies
840 : !----------------------------------------------------------------------
841 0 : if( p_in(k) > p(1) ) then
842 : pind = 2
843 : wght1 = 1._r8
844 0 : else if( p_in(k) <= p(nump) ) then
845 : pind = nump
846 : wght1 = 0._r8
847 : else
848 0 : do iz = izl,nump
849 0 : if( p(iz) < p_in(k) ) then
850 : izl = iz
851 : exit
852 : end if
853 : end do
854 0 : pind = max( min( iz,nump ),2 )
855 0 : wght1 = max( 0._r8,min( 1._r8,(p_in(k) - p(pind)) * del_p(pind-1) ) )
856 : end if
857 : !----------------------------------------------------------------------
858 : ! ... find "o3 ratios"
859 : !----------------------------------------------------------------------
860 0 : v3ratu = colo3_in(k) / colo3(pind-1)
861 0 : do iv = 1,numcolo3
862 0 : if( o3rat(iv) > v3ratu ) then
863 : exit
864 : end if
865 : end do
866 0 : ratindu = max( min( iv,numcolo3 ) - 1,1 )
867 :
868 0 : if( colo3(pind) /= 0._r8 ) then
869 0 : v3ratl = colo3_in(k) / colo3(pind)
870 0 : do iv = 1,numcolo3
871 0 : if( o3rat(iv) > v3ratl ) then
872 : exit
873 : end if
874 : end do
875 0 : ratindl = max( min( iv,numcolo3 ) - 1,1 )
876 : else
877 0 : ratindl = ratindu
878 0 : v3ratl = o3rat(ratindu)
879 : end if
880 :
881 : !----------------------------------------------------------------------
882 : ! ... compute the weigths
883 : !----------------------------------------------------------------------
884 0 : ial = albind
885 0 : ialp1 = ial + 1
886 0 : iv = ratindl
887 :
888 0 : dels(2) = max( 0._r8,min( 1._r8,(v3ratl - o3rat(iv)) * del_o3rat(iv) ) )
889 0 : dels(3) = max( 0._r8,min( 1._r8,(alb_in(k) - alb(ial)) * del_alb(ial) ) )
890 :
891 0 : wrk1 = (1._r8 - dels(2))*(1._r8 - dels(3))
892 0 : wghtl(0,0,0) = wrk0*wrk1
893 0 : wghtl(1,0,0) = dels(1)*wrk1
894 0 : wrk1 = (1._r8 - dels(2))*dels(3)
895 0 : wghtl(0,0,1) = wrk0*wrk1
896 0 : wghtl(1,0,1) = dels(1)*wrk1
897 0 : wrk1 = dels(2)*(1._r8 - dels(3))
898 0 : wghtl(0,1,0) = wrk0*wrk1
899 0 : wghtl(1,1,0) = dels(1)*wrk1
900 0 : wrk1 = dels(2)*dels(3)
901 0 : wghtl(0,1,1) = wrk0*wrk1
902 0 : wghtl(1,1,1) = dels(1)*wrk1
903 :
904 0 : iv = ratindu
905 0 : dels(2) = max( 0._r8,min( 1._r8,(v3ratu - o3rat(iv)) * del_o3rat(iv) ) )
906 :
907 0 : wrk1 = (1._r8 - dels(2))*(1._r8 - dels(3))
908 0 : wghtu(0,0,0) = wrk0*wrk1
909 0 : wghtu(1,0,0) = dels(1)*wrk1
910 0 : wrk1 = (1._r8 - dels(2))*dels(3)
911 0 : wghtu(0,0,1) = wrk0*wrk1
912 0 : wghtu(1,0,1) = dels(1)*wrk1
913 0 : wrk1 = dels(2)*(1._r8 - dels(3))
914 0 : wghtu(0,1,0) = wrk0*wrk1
915 0 : wghtu(1,1,0) = dels(1)*wrk1
916 0 : wrk1 = dels(2)*dels(3)
917 0 : wghtu(0,1,1) = wrk0*wrk1
918 0 : wghtu(1,1,1) = dels(1)*wrk1
919 :
920 0 : iz = pind
921 0 : iv = ratindl
922 0 : ivp1 = iv + 1
923 0 : do wn = 1,nw
924 0 : psum_l(wn) = wghtl(0,0,0) * rsf_tab(wn,iz,is,iv,ial) &
925 0 : + wghtl(0,0,1) * rsf_tab(wn,iz,is,iv,ialp1) &
926 0 : + wghtl(0,1,0) * rsf_tab(wn,iz,is,ivp1,ial) &
927 : + wghtl(0,1,1) * rsf_tab(wn,iz,is,ivp1,ialp1) &
928 0 : + wghtl(1,0,0) * rsf_tab(wn,iz,isp1,iv,ial) &
929 : + wghtl(1,0,1) * rsf_tab(wn,iz,isp1,iv,ialp1) &
930 : + wghtl(1,1,0) * rsf_tab(wn,iz,isp1,ivp1,ial) &
931 0 : + wghtl(1,1,1) * rsf_tab(wn,iz,isp1,ivp1,ialp1)
932 : end do
933 :
934 0 : iz = iz - 1
935 0 : iv = ratindu
936 0 : ivp1 = iv + 1
937 0 : do wn = 1,nw
938 0 : psum_u = wghtu(0,0,0) * rsf_tab(wn,iz,is,iv,ial) &
939 0 : + wghtu(0,0,1) * rsf_tab(wn,iz,is,iv,ialp1) &
940 0 : + wghtu(0,1,0) * rsf_tab(wn,iz,is,ivp1,ial) &
941 : + wghtu(0,1,1) * rsf_tab(wn,iz,is,ivp1,ialp1) &
942 0 : + wghtu(1,0,0) * rsf_tab(wn,iz,isp1,iv,ial) &
943 : + wghtu(1,0,1) * rsf_tab(wn,iz,isp1,iv,ialp1) &
944 : + wghtu(1,1,0) * rsf_tab(wn,iz,isp1,ivp1,ial) &
945 0 : + wghtu(1,1,1) * rsf_tab(wn,iz,isp1,ivp1,ialp1)
946 0 : rsf(wn,k) = (psum_l(wn) + wght1*(psum_u - psum_l(wn)))
947 : end do
948 : !------------------------------------------------------------------------------
949 : ! etfphot comes in as photons/cm^2/sec/nm (rsf includes the wlintv factor -- nm)
950 : ! ... --> convert to photons/cm^2/s
951 : !------------------------------------------------------------------------------
952 0 : rsf(:,k) = etfphot(:) * rsf(:,k)
953 :
954 : end do Level_loop
955 :
956 0 : deallocate( psum_l )
957 :
958 0 : end subroutine interpolate_rsf
959 :
960 :
961 : end module mo_jlong
|
