module_phot_fastj.F
References to this file elsewhere.
1 !**********************************************************************************
2 ! This computer software was prepared by Battelle Memorial Institute, hereinafter
3 ! the Contractor, under Contract No. DE-AC05-76RL0 1830 with the Department of
4 ! Energy (DOE). NEITHER THE GOVERNMENT NOR THE CONTRACTOR MAKES ANY WARRANTY,
5 ! EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
6 !
7 ! Photolysis Option: Fast-J
8 ! * Primary investigators: Elaine G. Chapman and James C. Barnard
9 ! * Co-investigators: Jerome D. Fast, William I. Gustafson Jr.
10 ! Last update: September 2005
11 !
12 ! Contact:
13 ! Jerome D. Fast, PhD
14 ! Staff Scientist
15 ! Pacific Northwest National Laboratory
16 ! P.O. Box 999, MSIN K9-30
17 ! Richland, WA, 99352
18 ! Phone: (509) 372-6116
19 ! Email: Jerome.Fast@pnl.gov
20 !
21 ! The original Fast-J code was provided by Oliver Wild (Univ. of Calif. Irvine);
22 ! however, substantial modifications were necessary to make it compatible with
23 ! WRF-chem and to include the effect of prognostic aerosols on photolysis rates.
24 !
25 ! Please report any bugs or problems to Jerome Fast, the WRF-chem implmentation
26 ! team leader for PNNL
27 !
28 ! References:
29 ! * Wild, O., X. Zhu, M.J. Prather, (2000), Accurate simulation of in- and below
30 ! cloud photolysis in tropospheric chemical models, J. Atmos. Chem., 37, 245-282.
31 ! * Barnard, J.C., E.G. Chapman, J.D. Fast, J.R. Schmelzer, J.R. Schulsser, and
32 ! R.E. Shetter (2004), An evaluation of the FAST-J photolysis model for
33 ! predicting nitrogen dioxide photolysis rates under clear and cloudy sky
34 ! conditions, Atmos. Environ., 38, 3393-3403.
35 ! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. Barnard, E.G.
36 ! Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution of ozone, particulates,
37 ! and aerosol direct radiative forcing in the vicinity of Houston using a fully-
38 ! coupled meteorology-chemistry-aerosol model, Submitted to J. Geophys. Res.
39 !
40 ! Contact Jerome Fast for updates on the status of manuscripts under review.
41 !
42 ! Additional information:
43 ! * www.pnl.gov/atmos_sciences/Jdf/wrfchem.html
44 !
45 ! Support:
46 ! Funding for adapting Fast-J was provided by the U.S. Department of Energy
47 ! under the auspices of Atmospheric Science Program of the Office of Biological
48 ! and Environmental Research the PNNL Laboratory Research and Directed Research
49 ! and Development program.
50 !**********************************************************************************
51
52 !WRF:MODEL_LAYER:CHEMISTRY
53 !
54 module module_phot_fastj
55 integer, parameter :: lunerr = -1
56
57 contains
58 !***********************************************************************
59 subroutine fastj_driver(id,ktau,dtstep,config_flags, &
60 gmt,julday,t_phy,moist,p8w,p_phy, &
61 chem,rho_phy,dz8w,xlat,xlong,z_at_w, &
62 ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
63 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
64 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
65 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
66 ph_n2o5, &
67 tauaer1,tauaer2,tauaer3,tauaer4, &
68 gaer1,gaer2,gaer3,gaer4, &
69 waer1,waer2,waer3,waer4, &
70 ids,ide, jds,jde, kds,kde, &
71 ims,ime, jms,jme, kms,kme, &
72 its,ite, jts,jte, kts,kte )
73
74
75 USE module_configure
76 USE module_state_description
77 USE module_data_mosaic_therm, only: nbin_a, nbin_a_maxd
78 USE module_data_mosaic_asect
79 USE module_data_mosaic_other
80 USE module_fastj_mie
81 USE module_mosaic_therm, only: aerosol_optical_properties
82 USE module_mosaic_driver, only: mapaer_tofrom_host
83 USE module_fastj_data, only: nb, nc
84
85 implicit none
86 !jdf
87 ! integer, parameter :: iprint = 0
88 integer, parameter :: single = 4 !compiler dependent value real*4
89 integer, parameter :: double = 8 !compiler dependent value real*8
90 integer,parameter :: ipar_fastj=1,jpar=1
91 integer,parameter :: jppj=14 !Number of photolytic reactions supplied
92 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
93 integer,save :: lpar !Number of levels in CTM
94 integer,save :: jpnl !Number of levels requiring chemistry
95 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
96 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
97 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
98 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
99 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
100 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
101 integer month_fastj ! Number of month (1-12)
102 integer iday_fastj ! Day of year
103 integer nspint ! Num of spectral intervals across solar
104 parameter ( nspint = 4 ) ! spectrum for FAST-J
105 real, dimension (nspint),save :: wavmid !cm
106 real, dimension (nspint, kmaxd+1),save :: sizeaer,extaer,waer,gaer,tauaer
107 real, dimension (nspint, kmaxd+1),save :: l2,l3,l4,l5,l6,l7
108 data wavmid &
109 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
110 integer nphoto_fastj
111 parameter (nphoto_fastj = 14)
112 integer &
113 lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, &
114 lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, &
115 lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, &
116 lfastj_hno4
117 parameter( lfastj_no2 = 1 )
118 parameter( lfastj_o3a = 2 )
119 parameter( lfastj_o3b = 3 )
120 parameter( lfastj_h2o2 = 4 )
121 parameter( lfastj_hchoa = 5 )
122 parameter( lfastj_hchob = 6 )
123 parameter( lfastj_ch3ooh= 7 )
124 parameter( lfastj_no3x = 8 )
125 parameter( lfastj_no3l = 9 )
126 parameter( lfastj_hono = 10 )
127 parameter( lfastj_n2o5 = 11 )
128 parameter( lfastj_hno3 = 12 )
129 parameter( lfastj_hno4 = 13 )
130 !jdf
131 INTEGER, INTENT(IN ) :: id,julday, &
132 ids,ide, jds,jde, kds,kde, &
133 ims,ime, jms,jme, kms,kme, &
134 its,ite, jts,jte, kts,kte
135 INTEGER, INTENT(IN ) :: &
136 ktau
137 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
138 INTENT(IN ) :: moist
139 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
140 INTENT(INOUT ) :: &
141 ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
142 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
143 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
144 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
145 ph_n2o5, &
146 tauaer1,tauaer2,tauaer3,tauaer4, &
147 gaer1,gaer2,gaer3,gaer4, &
148 waer1,waer2,waer3,waer4
149 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
150 INTENT(INOUT ) :: chem
151 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
152 INTENT(IN ) :: &
153 t_phy, &
154 p_phy, &
155 dz8w, &
156 p8w, &
157 rho_phy, &
158 z_at_w
159 REAL, DIMENSION( ims:ime , jms:jme ) , &
160 INTENT(IN ) :: xlat, &
161 xlong
162 REAL, INTENT(IN ) :: &
163 dtstep,gmt
164
165 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
166
167
168 !local variables
169
170 integer iclm, jclm
171
172 real number_bin(nbin_a_maxd,kmaxd) !These have to be the full kmaxd to
173 real radius_wet(nbin_a_maxd,kmaxd) !match arrays in MOSAIC routines.
174 complex refindx(nbin_a_maxd,kmaxd)
175
176 integer i,j, k, nsub, isecfrm0, ixhour
177 real xtime, xhour, xmin, gmtp, tmidh
178 real sla, slo
179 real psfc
180
181 real cos_sza
182 real, dimension(kts:kte-1) :: temp, ozone, dz
183 real, dimension(0:kte-1) :: pbnd
184 real, dimension(kts:kte-1) :: cloudmr, airdensity, relhum
185 real, dimension(kts:kte) :: zatw
186
187 real valuej(kte-1,nphoto_fastj)
188
189 logical processingAerosols
190
191 ! set "pegasus" grid size variables
192 ! itot = ite
193 ! jtot = jte
194 lpar = kte - 1
195 jpnl = kte - 1
196 nb = lpar + 1 !for module module_fastj_data
197 nc = 2*nb !ditto, and don't confuse this with nc in module_fastj_mie
198 nsubareas = 1
199 if ((kte-1 > kmaxd) .or. (lpar <= 0)) then
200 write( wrf_err_message, '(a,4i5)' ) &
201 '*** subr fastj_driver -- ' // &
202 'lpar, kmaxd, kts, kte', lpar, kmaxd, kts, kte
203 call wrf_message( trim(wrf_err_message) )
204 wrf_err_message = '*** subr fastj_driver -- ' // &
205 'kte-1>kmaxd OR lpar<=0'
206 call wrf_error_fatal( wrf_err_message )
207 end if
208
209 ! Determine if aerosol data is provided in the chem array. Currently,
210 ! only MOSAIC will work. The Mie routine does not know how to handle
211 ! SORGAM aerosols.
212 select case (config_flags%chem_opt)
213 case (CBMZ_MOSAIC_AA,CBMZ_MOSAIC_BB)
214 processingAerosols = .true.
215 case default
216 processingAerosols = .false.
217 end select
218
219 ! Set nbin_a = ntot_amode and check nbin_a <= nbin_a_maxd.
220 ! This duplicates the same assignment and check in module_mosaic_therm.F,
221 ! but photolysis is called before aeorosols so this must be set too.
222 !
223 ! rce 2004-dec-07 -- nbin_a is initialized elsewhere
224 !!$ nbin_a = ntot_amode
225 !!$ if ((nbin_a .gt. nbin_a_maxd) .or. (nbin_a .le. 0)) then
226 !!$ write( wrf_err_message, '(a,2(1x,i4))' ) &
227 !!$ '*** subr fastj_driver -- nbin_a, nbin_a_maxd =', &
228 !!$ nbin_a, nbin_a_maxd
229 !!$ call wrf_message( wrf_err_message )
230 !!$ call wrf_error_fatal( &
231 !!$ '*** subr fastj_driver -- BAD VALUE for nbin_a' )
232 !!$ end if
233
234 ! determine current time of day in Greenwich Mean Time at middle
235 ! of current time step, tmidh. do this by computing the number of minutes
236 ! from beginning of run to middle of current time step
237 xtime=(ktau-1)*dtstep/60. + dtstep/120.
238 ixhour = ifix(gmt + 0.01) + ifix(xtime/60.)
239 xhour=float(ixhour) !current hour
240 xmin = 60.*gmt + (xtime-xhour*60)
241 gmtp=mod(xhour,24.)
242 tmidh= gmtp + xmin/60.
243 isecfrm0 = ifix ( (ktau-1) * dtstep )
244
245 ! execute for each i,j column and each nsub subarea
246 do nsub = 1, nsubareas
247 do jclm = jts, jte
248 do iclm = its, ite
249
250 do k = kts, lpar
251 dz(k) = dz8w(iclm, k, jclm) ! cell depth (m)
252 end do
253
254 if( processingAerosols ) then
255 ! take chem data and extract 1 column to create rsub(l,k,m) array
256 call mapaer_tofrom_host( 0, &
257 ims,ime, jms,jme, kms,kme, &
258 its,ite, jts,jte, kts,kte, &
259 iclm, jclm, kts,lpar, &
260 num_moist, num_chem, moist, chem, &
261 t_phy, p_phy, rho_phy )
262 ! generate aerosol optical properties for cells in this column
263 ! subroutine is located in file module_mosaic_therm
264 call aerosol_optical_properties(iclm, jclm, lpar, refindx, &
265 radius_wet, number_bin)
266 ! execute mie code , located in file module_fastj_mie
267 CALL wrf_debug(150,'fastj_driver: calling mieaer')
268 call mieaer( &
269 iclm, jclm, nbin_a, &
270 number_bin, radius_wet, refindx, &
271 dz, isecfrm0, lpar, &
272 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
273 end if
274
275 ! set lat, long
276 sla = xlat(iclm,jclm)
277 slo = xlong(iclm,jclm)
278 ! set column pressures, temperature, and ozone
279 psfc = p8w(iclm,1,jclm) * 10. ! convert pascals to dynes/cm2
280 do k = kts, lpar
281 pbnd(k) = p8w(iclm,k+1,jclm) *10. ! convert pascals to dynes/cm2
282 temp(k) = t_phy(iclm,k,jclm)
283 ozone(k) = chem(iclm,k,jclm,p_o3) / 1.0e6 ! ppm->mol/mol air
284 cloudmr(k) = moist(iclm,k,jclm,p_qc)/0.622
285 airdensity(k) = rho_phy(iclm,k,jclm)/28.966e3
286 relhum(k) = MIN( .95, moist(iclm,k,jclm,p_qv) / &
287 (3.80*exp(17.27*(t_phy(iclm,k,jclm)-273.)/ &
288 (t_phy(iclm,k,jclm)-36.))/(.01*p_phy(iclm,k,jclm))))
289 relhum(k) = MAX(.001,relhum(k))
290 zatw(k)=z_at_w(iclm,k,jclm)
291 end do
292 zatw(lpar+1)=z_at_w(iclm,lpar+1,jclm)
293 ! call interface_fastj
294 CALL wrf_debug(150,'fastj_driver: calling interface_fastj')
295 call interface_fastj(tmidh,sla,slo,julday, &
296 pbnd, psfc, temp, ozone, &
297 dz, cloudmr, airdensity, relhum, zatw, &
298 iclm, jclm, lpar, jpnl, &
299 isecfrm0, valuej, cos_sza, processingAerosols, &
300 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
301 ! put column photolysis rates (valuej) into wrf photolysis (i,k,j) arrays
302 CALL wrf_debug(150,'fastj_driver: calling mapJrates_tofrom_host')
303 call mapJrates_tofrom_host( 0, &
304 ims,ime, jms,jme, kms,kme, &
305 its,ite, jts,jte, kts,kte, &
306 iclm, jclm, kts,lpar, &
307 valuej, &
308 ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
309 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
310 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
311 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
312 ph_n2o5 )
313 ! put the aerosol optical properties into the wrf arrays (this is hard-
314 ! coded to 4 spectral bins, nspint)
315 do k=kts,kte-1
316 tauaer1(iclm,k,jclm) = tauaer(1,k)
317 tauaer2(iclm,k,jclm) = tauaer(2,k)
318 tauaer3(iclm,k,jclm) = tauaer(3,k)
319 tauaer4(iclm,k,jclm) = tauaer(4,k)
320 gaer1(iclm,k,jclm) = gaer(1,k)
321 gaer2(iclm,k,jclm) = gaer(2,k)
322 gaer3(iclm,k,jclm) = gaer(3,k)
323 gaer4(iclm,k,jclm) = gaer(4,k)
324 waer1(iclm,k,jclm) = waer(1,k)
325 waer2(iclm,k,jclm) = waer(2,k)
326 waer3(iclm,k,jclm) = waer(3,k)
327 waer4(iclm,k,jclm) = waer(4,k)
328 end do
329
330 end do
331 end do
332 end do
333
334 return
335 end subroutine fastj_driver
336
337 !-----------------------------------------------------------------------
338 subroutine mapJrates_tofrom_host( iflag, &
339 ims,ime, jms,jme, kms,kme, &
340 its,ite, jts,jte, kts,kte, &
341 iclm, jclm, ktmaps,ktmape, &
342 valuej, &
343 ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
344 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
345 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
346 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
347 ph_n2o5 )
348
349 USE module_data_cbmz
350
351 implicit none
352 !jdf
353 integer nphoto_fastj
354 parameter (nphoto_fastj = 14)
355 integer &
356 lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, &
357 lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, &
358 lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, &
359 lfastj_hno4
360 parameter( lfastj_no2 = 1 )
361 parameter( lfastj_o3a = 2 )
362 parameter( lfastj_o3b = 3 )
363 parameter( lfastj_h2o2 = 4 )
364 parameter( lfastj_hchoa = 5 )
365 parameter( lfastj_hchob = 6 )
366 parameter( lfastj_ch3ooh= 7 )
367 parameter( lfastj_no3x = 8 )
368 parameter( lfastj_no3l = 9 )
369 parameter( lfastj_hono = 10 )
370 parameter( lfastj_n2o5 = 11 )
371 parameter( lfastj_hno3 = 12 )
372 parameter( lfastj_hno4 = 13 )
373 !jdf
374 INTEGER, INTENT(IN ) :: iflag, &
375 ims,ime, jms,jme, kms,kme, &
376 its,ite, jts,jte, kts,kte, &
377 iclm, jclm, ktmaps, ktmape
378 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
379 INTENT(INOUT ) :: &
380 ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
381 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
382 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
383 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
384 ph_n2o5
385
386 REAL, DIMENSION( kte-1,nphoto_fastj ), INTENT(INOUT) :: valuej
387
388 ! local variables
389 real ft
390 integer kt
391
392 ft = 60.
393
394 if (iflag .gt. 0) go to 2000
395 ! flag is <=0, put pegasus column J rates (in 1/sec) into WRF arrays (in 1/min)
396 do kt = ktmaps, ktmape
397 ph_no2(iclm,kt,jclm) = valuej(kt,lfastj_no2) * ft
398 ph_no3o(iclm,kt,jclm) = valuej(kt,lfastj_no3x) * ft
399 ph_no3o2(iclm,kt,jclm) = valuej(kt,lfastj_no3l) * ft
400 ph_o33p(iclm,kt,jclm) = valuej(kt,lfastj_o3a) * ft
401 ph_o31d(iclm,kt,jclm) = valuej(kt,lfastj_o3b) * ft
402 ph_hno2(iclm,kt,jclm) = valuej(kt,lfastj_hono) * ft
403 ph_hno3(iclm,kt,jclm) = valuej(kt,lfastj_hno3) * ft
404 ph_hno4(iclm,kt,jclm) = valuej(kt,lfastj_hno4) * ft
405 ph_h2o2(iclm,kt,jclm) = valuej(kt,lfastj_h2o2) * ft
406 ph_ch3o2h(iclm,kt,jclm) = valuej(kt,lfastj_ch3ooh) * ft
407 ph_ch2or(iclm,kt,jclm) = valuej(kt,lfastj_hchoa) * ft
408 ph_ch2om(iclm,kt,jclm) = valuej(kt,lfastj_hchob) * ft
409 ph_n2o5(iclm,kt,jclm) = valuej(kt,lfastj_n2o5) * ft
410
411 ph_o2(iclm,kt,jclm) = 0.0
412 ph_ch3cho(iclm,kt,jclm) = 0.0
413 ph_ch3coch3(iclm,kt,jclm) = 0.0
414 ph_ch3coc2h5(iclm,kt,jclm) = 0.0
415 ph_hcocho(iclm,kt,jclm) = 0.0
416 ph_ch3cocho(iclm,kt,jclm) = 0.0
417 ph_hcochest(iclm,kt,jclm) = 0.0
418 ph_ch3coo2h(iclm,kt,jclm) = 0.0
419 ph_ch3ono2(iclm,kt,jclm) = 0.0
420 ph_hcochob(iclm,kt,jclm) = 0.0
421
422 end do
423 return !finished peg-> wrf mapping
424
425 2000 continue
426 ! iflag > 0 ; put wrf ph_xxx Jrates (1/min) into pegasus column valuej (1/sec)
427 do kt = ktmaps, ktmape
428 valuej(kt,lfastj_no2) = ph_no2(iclm,kt,jclm) / ft
429 valuej(kt,lfastj_no3x) = ph_no3o(iclm,kt,jclm) / ft
430 valuej(kt,lfastj_no3l) = ph_no3o2(iclm,kt,jclm)/ ft
431 valuej(kt,lfastj_o3a) = ph_o33p(iclm,kt,jclm) / ft
432 valuej(kt,lfastj_o3b) = ph_o31d(iclm,kt,jclm) / ft
433 valuej(kt,lfastj_hono) = ph_hno2(iclm,kt,jclm) / ft
434 valuej(kt,lfastj_hno3) = ph_hno3(iclm,kt,jclm) / ft
435 valuej(kt,lfastj_hno4) = ph_hno4(iclm,kt,jclm) / ft
436 valuej(kt,lfastj_h2o2) = ph_h2o2(iclm,kt,jclm) / ft
437 valuej(kt,lfastj_ch3ooh) = ph_ch3o2h(iclm,kt,jclm)/ft
438 valuej(kt,lfastj_hchoa) = ph_ch2or(iclm,kt,jclm)/ ft
439 valuej(kt,lfastj_hchob) = ph_ch2om(iclm,kt,jclm)/ ft
440 valuej(kt,lfastj_n2o5) = ph_n2o5(iclm,kt,jclm) / ft
441 end do
442
443 return !finished wrf->peg mapping
444
445 end subroutine mapJrates_tofrom_host
446 !-----------------------------------------------------------------------
447
448
449
450 !***********************************************************************
451 subroutine interface_fastj(tmidh,sla,slo,julian_day, &
452 pbnd, psfc, temp, ozone, &
453 dz, cloudmr, airdensity, relhum, zatw, &
454 isvode, jsvode, lpar, jpnl, &
455 isecfrm0, valuej, cos_sza, processingAerosols, &
456 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
457 !-----------------------------------------------------------------
458 ! sets parameters for fastj call.
459 ! inputs
460 ! tmidh -- GMT time in decimal hours at which to calculate
461 ! photolysis rates
462 ! sla -- latitude, decimal degrees in real*8
463 ! slo -- negative of the longitude, decimal degrees in real*8
464 ! julian_day -- day of the year in julian days
465 ! pbnd(0:lpar) = pressure at top boundary of cell k (dynes/cm^2).
466 ! psfc = surface pressure (dynes/cm^2).
467 ! temp(lpar)= mid-cell temperature values (deg K)
468 ! ozone(lpar) = mid-cell ozone mixing ratios
469 ! surface_albedo -- broadband albedo (dimensionless)
470 ! isecfrm0 -- elapsed time from start of simulation in seconds.
471 ! isvode,jsvode -- current column i,j.
472 !
473 ! lpar -- vertical extent of column (from module_fastj_cmnh)
474 !
475 ! outputs
476 ! cos_sza -- cosine of solar zenith angle.
477 ! valuej(lpar,nphoto_fastj-1) -- array of photolysis rates, s-1.
478 !
479 !
480 ! local variables
481 ! surface_pressure_mb -- surface pressure (mb). equal to col_press_mb(1).
482 ! col_press_mb(lpar) -- for the column, grid cell boundary pressures
483 ! (not at cell centers) up until the bottom pressure for the
484 ! top cell (mb).
485 ! col_temp_K(lpar+1) -- for the column, grid cell center temperature (deg K)
486 ! col_ozone(lpar+1) -- for the column, grid cell center ozone mixing
487 ! ratios (dimensionless)
488 ! col_optical_depth(lpar+1) -- for the column, grid cell center cloud
489 ! optical depths (dimensionless).SET TO ZERO IN THIS VERSION
490 ! tauaer_550 -- aerosol optical thickness at 550 nm.
491 ! note: photolysis rates are calculated at centers of model layers
492 ! the pressures are given at the boundaries defining
493 ! the top and bottom of the layers
494 ! so the number of pressure values is equal
495 ! to the (number of layers) + 1 ; the last pressure is set = 0 in fastj code.
496 ! pressures from the surface up to the bottom of the top (lpar<=kmaxd) cell
497 ! ******** pressure 2
498 ! layer 1 - temperature,optical depth, and O3 given here
499 ! ******** pressure 1
500 ! the optical depth is appropriate for the layer depth
501 ! conversion factor: 1 dyne/cm2 = 0.001 mb
502 !-----------------------------------------------------------------
503 USE module_data_mosaic_other, only : kmaxd
504 USE module_peg_util, only: peg_message, peg_error_fatal
505
506 IMPLICIT NONE
507
508 !jdf
509 integer, parameter :: iprint = 0
510 integer, parameter :: single = 4 !compiler dependent value real*4
511 integer, parameter :: double = 8 !compiler dependent value real*8
512 integer,parameter :: ipar_fastj=1,jpar=1
513 integer,parameter :: jppj=14 !Number of photolytic reactions supplied
514 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
515 integer lpar !Number of levels in CTM
516 integer jpnl !Number of levels requiring chemistry
517 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
518 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
519 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
520 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
521 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
522 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
523 integer month_fastj ! Number of month (1-12)
524 integer iday_fastj ! Day of year
525 integer nphoto_fastj
526 parameter (nphoto_fastj = 14)
527 integer &
528 lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, &
529 lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, &
530 lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, &
531 lfastj_hno4
532 parameter( lfastj_no2 = 1 )
533 parameter( lfastj_o3a = 2 )
534 parameter( lfastj_o3b = 3 )
535 parameter( lfastj_h2o2 = 4 )
536 parameter( lfastj_hchoa = 5 )
537 parameter( lfastj_hchob = 6 )
538 parameter( lfastj_ch3ooh= 7 )
539 parameter( lfastj_no3x = 8 )
540 parameter( lfastj_no3l = 9 )
541 parameter( lfastj_hono = 10 )
542 parameter( lfastj_n2o5 = 11 )
543 parameter( lfastj_hno3 = 12 )
544 parameter( lfastj_hno4 = 13 )
545 integer nspint ! Num of spectral intervals across solar
546 parameter ( nspint = 4 ) ! spectrum for FAST-J
547 real, dimension (nspint),save :: wavmid !cm
548 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
549 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
550 data wavmid &
551 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
552 !jdf
553 real pbnd(0:lpar), psfc
554 real temp(lpar), ozone(lpar), surface_albedo
555 real dz(lpar), cloudmr(lpar), airdensity(lpar), relhum(lpar), zatw(lpar+1)
556 integer isecfrm0, isvode, jsvode
557
558 real cos_sza
559 integer,parameter :: lunout=41
560
561 real valuej(lpar,nphoto_fastj)
562
563 real hl,rhl,factor1,part1,part2,cfrac,rhfrac
564 real emziohl(lpar+1),clwp(lpar)
565 !ec material to check output
566 real valuej_no3rate(lpar)
567
568 real*8 lat,lon
569 real*8 jvalue(lpar,nphoto_fastj)
570 real sza
571 real tau1
572 real tmidh, sla, slo
573
574 integer julian_day,iozone1
575 integer,parameter :: nfastj_rxns = 14
576 integer k, l
577
578 real surface_pressure_mb, tauaer_550, &
579 col_press_mb,col_temp_K,col_ozone,col_optical_depth
580 dimension col_press_mb(lpar+2),col_temp_K(lpar+1), &
581 col_ozone(lpar+1),col_optical_depth(lpar+1)
582 character*80 msg
583
584 ! define logical processingAerosols
585 ! if processingAerosols = true, uses values calculated in subroutine
586 ! mieaer for variables & arrays in common block mie.
587 ! if processingAerosols = false, sets all variables & arrays in common
588 ! block mie to zero. (JCB-revised Fast-J requires common block mie info,
589 ! regardless of whether aerosols are present or not. Original Wild Fast-J
590 ! did not use common block mie info.)
591
592 logical processingAerosols
593
594 ! set lat and longitude as real*8 for consistency with fastj code.
595 ! variables lat and lon previously declared as reals
596 lat = sla
597 lon = slo
598 !
599 ! cloud optical depths currently treated by using fractional cloudiness
600 ! based on relative humidity. cloudmr set up to use cloud liquid water
601 ! but hooks into microphysics need to be tested - for now set cloudmr=0
602 !
603 ! parameters to calculate 'typical' liquid cloudwater path values for
604 ! non convective clouds based on approximations in NCAR's CCM2
605 ! 0.18 = reference liquid water concentration (gh2o/m3)
606 ! hl = liquid water scale height (m)
607 !
608 hl=1080.+2000.0*cos(lat*0.017454329)
609 rhl=1.0/hl
610 do k =1, lpar+1
611 emziohl(k)=exp(-zatw(k)*rhl)
612 enddo
613 do k =1, lpar
614 clwp(k)=0.18*hl*(emziohl(k)-emziohl(k+1))
615 enddo
616 ! assume radius of cloud droplets is constant at 10 microns (0.001 cm) and
617 ! that density of water is constant at 1 g2ho/cm3
618 ! factor1=3./2./0.001/1.
619 factor1=1500.0
620 do k =1, lpar
621 col_optical_depth(k) = 0.0
622 cfrac=0.0
623 cloudmr(k)=0.0
624 if(cloudmr(k).gt.0.0) cfrac=1.0
625 ! 18.0*airdensity converts mole h2o/mole air to g h2o/cm3, part1 is in g h2o/cm2
626 part1=cloudmr(k)*cfrac*18.0*airdensity(k)*dz(k)*100.0
627 if(relhum(k).lt.0.8) then
628 rhfrac=0.0
629 elseif(relhum(k).le.1.0.and.relhum(k).ge.0.8) then
630 ! rhfrac=(relhum(k)-0.8)/(1.0-0.8)
631 rhfrac=(relhum(k)-0.8)/0.2
632 else
633 rhfrac=1.0
634 endif
635 if(rhfrac.ge.0.01) then
636 ! factor 1.0e4 converts clwp of g h2o/m2 to g h2o/cm2
637 part2=rhfrac*clwp(k)/1.e4
638 else
639 part2=0.0
640 endif
641 if(cfrac.gt.0) part2=0.0
642 col_optical_depth(k) = factor1*(part1+part2)
643 ! col_optical_depth(k) = 0.0
644 ! if(isvode.eq.33.and.jsvode.eq.34) &
645 ! print*,'jdf opt',isvode,jsvode,k,col_optical_depth(k), &
646 ! cfrac,rhfrac,relhum(k),clwp(k)
647 end do
648 col_optical_depth(lpar+1) = 0.0
649 if (.not.processingAerosols) then
650 ! set common block mie variables to 0 if
651 call set_common_mie(lpar, &
652 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
653 end if ! processingAerosols=false
654
655 ! set pressure, temperature, ozone of each cell in the column
656 ! set iozone1 = lpar to allow replacement of climatological ozone with model
657 ! predicted ozone to top of chemistry column; standard fastj climatological o3
658 ! thereafter.
659 surface_pressure_mb = psfc * 0.001
660 tau1 = tmidh
661 col_press_mb(1) = psfc * 0.001
662 iozone1 = lpar
663 do k =1, lpar
664 col_press_mb(k+1) = pbnd(k) * 0.001
665 col_temp_K(k) = temp(k)
666 col_ozone(k) = ozone(k)
667 end do
668
669 ! surface_albedo=0.055
670 !jdf
671 surface_albedo=0.05
672
673 ! set aerosol parameters needed by Fast-J
674 if (processingAerosols) then
675 tauaer_550 = 0.0 ! needed parameters already calculated by subroutine
676 ! mieaer and passed into proper parts of fastj code
677 ! via module_fastj_cmnmie
678 else
679 tauaer_550 = 0.05 ! no aerosols, assume typical constant aerosol optical thickness
680 end if
681
682 CALL wrf_debug(200,'interface_fastj: calling fastj')
683 call fastj(isvode,jsvode,lat,lon,surface_pressure_mb,surface_albedo, &
684 julian_day, tau1, &
685 col_press_mb, col_temp_K, col_optical_depth, col_ozone, &
686 iozone1,tauaer_550,jvalue,sza,lpar,jpnl, &
687 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
688
689 cos_sza = cosd(sza)
690
691
692 ! array jvalue (real*8) is returned from fastj. array valuej(unspecified
693 ! real, default of real*4) is sent on to
694 ! other chemistry subroutines
695 do k = 1, lpar
696 valuej(k, lfastj_no2) = jvalue(k,lfastj_no2)
697 valuej(k, lfastj_o3a) = jvalue(k,lfastj_o3a)
698 valuej(k, lfastj_o3b) = jvalue(k,lfastj_o3b)
699 valuej(k, lfastj_h2o2) = jvalue(k,lfastj_h2o2)
700 valuej(k, lfastj_hchoa) = jvalue(k,lfastj_hchoa)
701 valuej(k, lfastj_hchob) = jvalue(k,lfastj_hchob)
702 valuej(k, lfastj_ch3ooh) = jvalue(k,lfastj_ch3ooh)
703 valuej(k, lfastj_no3x) = jvalue(k,lfastj_no3x)
704 valuej(k, lfastj_no3l) = jvalue(k,lfastj_no3l)
705 valuej(k, lfastj_hono) = jvalue(k,lfastj_hono)
706 valuej(k, lfastj_n2o5) = jvalue(k,lfastj_n2o5)
707 valuej(k, lfastj_hno3) = jvalue(k,lfastj_hno3)
708 valuej(k, lfastj_hno4) = jvalue(k,lfastj_hno4)
709 end do
710 ! diagnostic output and zeroed value if negative photolysis rates returned
711 do k = 1, lpar
712 valuej(k,nphoto_fastj)=0.0
713 do l = 1, nphoto_fastj-1
714 if (valuej(k,l) .lt. 0) then
715 write( msg, '(a,i8,4i4,1x,e11.4)' ) &
716 'FASTJ negative Jrate ' // &
717 'tsec i j k l J(k,l)', isecfrm0,isvode,jsvode,k,l,valuej(k,l)
718 call peg_message( lunerr, msg )
719 valuej(k,l) = 0.0
720 ! following code used if want run stopped with negative Jrate
721 ! msg = '*** subr interface_fastj -- ' // &
722 ! 'Negative J rate returned from Fast-J'
723 ! call peg_error_fatal( lunerr, msg )
724 end if
725 end do
726 end do
727 ! compute overall no3 photolysis rate
728 ! wig: commented out since it is not used anywhere
729 ! do k = 1, lpar
730 ! valuej_no3rate(k) = &
731 ! valuej(k, lfastj_no3x) + valuej(k,lfastj_no3l)
732 ! end do
733 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
734 !EC FOLLOWING PRINT LOOP SHOULD BE ELIMINATED FROM FINAL WRF CODE
735 ! print outs of aerosol optical properties are performed
736 ! in mie subroutines. now output jrates.
737 ! write(27,909)
738 ! 909 format( &
739 ! 'isecfrm0',2x, 'i', 2x, 'j',2x 'k',3x, 'cos_sza',7x, &
740 ! 'no2',13x,'o3a', 13x,'o3b',13x, &
741 ! 'h2o2' , 12x, 'hchoa',11x,'hchob', 11x,'ch3ooh',10x,'no3', &
742 ! 13x,'hono',12x,'n2o5',12x, 'hno3',12x,'hno4')
743 !
744 ! do k = 1, lpar
745 ! write(27, 911) isecfrm0,isvode,jsvode,k, cos_sza, &
746 ! valuej(k, lfastj_no2), &
747 ! valuej(k, lfastj_o3a) , &
748 ! valuej(k, lfastj_o3b), &
749 ! valuej(k, lfastj_h2o2), &
750 ! valuej(k, lfastj_hchoa), &
751 ! valuej(k, lfastj_hchob) , &
752 ! valuej(k, lfastj_ch3ooh), &
753 ! valuej_no3rate(k), &
754 ! valuej(k, lfastj_hono), &
755 ! valuej(k, lfastj_n2o5), &
756 ! valuej(k, lfastj_hno3) , &
757 ! valuej(k, lfastj_hno4)
758 !
759 !911 format(i7,3(2x,i4),2x,f7.4, 2x, &
760 ! 12(e14.6,2x))
761 ! end do
762 !EC END PRINT LOOP.
763 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
764 return
765 end subroutine interface_fastj
766 !***********************************************************************
767 subroutine set_common_mie(lpar, &
768 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
769 ! for use when aerosols are not included in a model run. sets variables
770 ! in common block mie to zero, except for wavelengths.
771 ! OUTPUT: in module_fastj_cmnmie
772 ! wavmid ! fast-J wavelengths (cm)
773 ! tauaer ! aerosol optical depth
774 ! waer ! aerosol single scattering albedo
775 ! gaer ! aerosol asymmetery factor
776 ! extaer ! aerosol extinction
777 ! l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,......
778 ! sizeaer ! average wet radius
779 ! INPUTS
780 ! lpar = total number of vertical layers in chemistry model. Passed
781 ! via module_fastf_cmnh
782 ! NB = total vertical layers + 1 considered by FastJ (lpar+1=kmaxd+1).
783 ! passed via module_fast j_data
784 !------------------------------------------------------------------------
785
786 USE module_data_mosaic_other, only : kmaxd
787
788 IMPLICIT NONE
789 !jdf
790 integer lpar ! Number of levels in CTM
791 integer nspint ! Num of spectral intervals across solar
792 parameter ( nspint = 4 ) ! spectrum for FAST-J
793 real, dimension (nspint),save :: wavmid !cm
794 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
795 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
796 data wavmid &
797 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
798 !jdf
799
800 ! LOCAL VARIABLES
801 integer klevel ! vertical level index
802 integer ns ! spectral loop index
803
804
805 ! aerosol optical properties: set everything = 0 when no aerosol
806 do 1000 ns=1,nspint
807 do 1000 klevel = 1, lpar
808 tauaer(ns,klevel)=0.
809 waer(ns,klevel)=0.
810 gaer(ns,klevel)=0.
811 sizeaer(ns,klevel)=0.0
812 extaer(ns,klevel)=0.0
813 l2(ns,klevel)=0.0
814 l3(ns,klevel)=0.0
815 l4(ns,klevel)=0.0
816 l5(ns,klevel)=0.0
817 l6(ns,klevel)=0.0
818 l7(ns,klevel)=0.0
819 1000 continue
820
821 return
822 end subroutine set_common_mie
823 !***********************************************************************
824 subroutine fastj(isvode,jsvode,lat,lon,surface_pressure,surface_albedo, &
825 julian_day,tau1, pressure, temperature, optical_depth, my_ozone1, &
826 iozone1,tauaer_550_1,jvalue,SZA_dum,lpar,jpnl, &
827 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
828 ! input:
829 ! lat = latitute; must be real*8
830 ! lon = longitude; must be real*8
831 ! surface_pressure (mb); real*4
832 ! surface_albedo (broadband albedo); real*4
833 ! julian_day; integer
834 ! tau1 = time of calculation (GMT); real*4
835 ! pressure (mb) = vector of pressure values, pressure(NB);
836 ! real*4; NB is the number of model layers;
837 ! pressure (NB+1) is defined as 0 mb in model
838 ! temperature (degree K)= vector of temperature values, temperature(NB);
839 ! real*4
840 ! optical_depth (dimensionless) = vector of cloud optical depths,
841 ! optical_depth(NB); real*4
842 ! my_ozone1 (volume mixing ratio) = ozone at layer center
843 ! ozone(iozone); real*4; if iozone1 <= NB-1, then climatology is
844 ! used in the upper layers
845 ! tauaer_550; real*4 aerosol optical thickness @ 550 nm
846 ! input note: NB is the number of model layers -- photolysis rates are calculated
847 ! at layer centers while pressures are given at the boundaries defining
848 ! the top and bottom of the layers. The number of pressure values =
849 ! (number of layers) + 1 ; see below
850 ! ******** pressure 2
851 ! layer 1 - optical depth, O3, and temperature given here
852 ! ******** pressure 1
853 ! temperature and o3 are defined at the layer center. optical depth is
854 ! appropriate for the layer depth.
855 ! output:
856 ! jvalue = photolysis rates, an array of dimension jvalue(jpnl,jppj) where
857 ! jpnl = # of models level at which photolysis rates are calculated
858 ! note: level 1 = first level of model (adjacent to ground)
859 ! jppj = # of chemical species for which photolysis rates are calculated;
860 ! this is fixed and is not easy to change on the fly
861 ! jpnl land jppl are defined in the common block "cmn_h.f"
862 ! SZA_dum = solar zenith angle
863 !-----------------------------------------------------------------------
864
865 USE module_data_mosaic_other, only : kmaxd
866 USE module_fastj_data
867
868 IMPLICIT NONE
869 !jdf
870 ! Print Fast-J diagnostics if iprint /= 0
871 integer, parameter :: iprint = 0
872 integer, parameter :: single = 4 !compiler dependent value real*4
873 ! integer, parameter :: double = 8 !compiler dependent value real*8
874 ! following specific for fastj
875 ! jppj will be gas phase mechanism dependent
876 integer,parameter :: ipar_fastj=1,jpar=1
877 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
878 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
879 ! The vertical level variables are set in fastj_driver.
880 integer lpar !Number of levels in CTM
881 integer jpnl !Number of levels requiring chemistry
882 ! following should be available from other wrf modules and passed into
883 ! photodriver
884 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
885 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
886 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
887 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
888 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
889 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
890 integer month_fastj ! Number of month (1-12)
891 integer iday_fastj ! Day of year
892 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
893 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
894 real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios
895 real tauaer_550
896 integer iozone
897 integer nslat ! Latitude of current profile point
898 integer nslon ! Longitude of current profile point
899 save :: nslat, nslon
900 integer nspint ! Num of spectral intervals across solar
901 parameter ( nspint = 4 ) ! spectrum for FAST-J
902 real, dimension (nspint),save :: wavmid !cm
903 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
904 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
905 data wavmid &
906 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
907 !jdf
908
909 integer julian_day
910 real surface_pressure,surface_albedo,pressure(lpar+2), &
911 temperature(lpar+1)
912 real optical_depth(lpar+1)
913 real tau1
914 real*8 pi_fastj,lat,lon,timej,jvalue(lpar,jppj)
915 integer isvode,jsvode
916
917 integer iozone1,i
918 real my_ozone1(lpar+1)
919
920 real tauaer_550_1
921 real sza_dum
922
923 integer ientryno_fastj
924 save ientryno_fastj
925 data ientryno_fastj / 0 /
926
927
928
929 ! Just focus on one column
930 nslat = 1
931 nslon = 1
932 pi_fastj=3.141592653589793D0
933 !
934
935 ! JCB - note that pj(NB+1) = p and is defined such elsewhere
936 do i=1,lpar+1
937 pj(i)=pressure(i)
938 T(nslon,nslat,i)=temperature(i)
939 OD(nslon,nslat,i)=optical_depth(i)
940 enddo
941 ! surface albedo
942 SA(nslon,nslat)=surface_albedo
943 !
944 iozone=iozone1
945 do i=1,iozone1
946 my_ozone(i)=my_ozone1(i)
947 enddo
948 !
949 tau_fastj=tau1 ! fix time
950 iday_fastj=julian_day
951 ! fix optical depth for situations where aerosols not considered
952 tauaer_550=tauaer_550_1
953 !
954 month_fastj=int(dble(iday_fastj)*12.d0/365.d0)+1 ! Approximately
955 xgrd(nslon)=lon*pi_fastj/180.d0
956 ygrd(nslat)=lat*pi_fastj/180.d0
957 ydgrd(nslat)=lat
958
959 ! Initial call to Fast-J to set things up--done only once
960 if (ientryno_fastj .eq. 0) then
961 call inphot2
962 ientryno_fastj = 1
963 end if
964 !
965 ! Now call fastj as appropriate
966 timej=0.0 ! manually set offset to zero (JCB: 14 November 2001)
967 call photoj(isvode,jsvode,jvalue,timej,nslat,nslon,iozone,tauaer_550, &
968 my_ozone,p,t,od,sa,lpar,jpnl, &
969 xgrd,ygrd,tau_fastj,month_fastj,iday_fastj,ydgrd, &
970 sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
971 sza_dum=SZA
972
973 return
974 end subroutine fastj
975 !**********************************************************************
976 !---(pphot.f)-------generic CTM shell from UCIrvine (p-code 4.0, 7/99)
977 !---------PPHOT calculates photolysis rates with the Fast-J scheme
978 !---------subroutines: inphot, photoj, Fast-J schemes...
979 !-----------------------------------------------------------------------
980 !
981 subroutine inphot2
982 !-----------------------------------------------------------------------
983 ! Routine to initialise photolysis rate data, called directly from the
984 ! cinit routine in ASAD. Currently use it to read the JPL spectral data
985 ! and standard O3 and T profiles and to set the appropriate reaction index.
986 !-----------------------------------------------------------------------
987 !
988 ! iph Channel number for reading all data files
989 ! rad Radius of Earth (cm)
990 ! zzht Effective scale height above top of atmosphere (cm)
991 ! dtaumax Maximum opt.depth above which a new level should be inserted
992 ! dtausub No. of opt.depths at top of cloud requiring subdivision
993 ! dsubdiv Number of additional levels to add at top of cloud
994 ! szamax Solar zenith angle cut-off, above which to skip calculation
995 !
996 !-----------------------------------------------------------------------
997 !
998
999 USE module_data_mosaic_other, only : kmaxd
1000 USE module_fastj_data
1001
1002 IMPLICIT NONE
1003 !jdf
1004 ! Print Fast-J diagnostics if iprint /= 0
1005 integer, parameter :: iprint = 0
1006 integer, parameter :: single = 4 !compiler dependent value real*4
1007 ! integer, parameter :: double = 8 !compiler dependent value real*8
1008 integer,parameter :: ipar_fastj=1,jpar=1
1009 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1010 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1011 integer lpar !Number of levels in CTM
1012 integer jpnl !Number of levels requiring chemistry
1013 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1014 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1015 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1016 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1017 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1018 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1019 integer month_fastj ! Number of month (1-12)
1020 integer iday_fastj ! Day of year
1021 !jdf
1022 !
1023 ! Set labels of photolysis rates required
1024 !ec032504 CALL RD_JS(iph,path_fastj_ratjd)
1025 ! call rd_js2
1026 !
1027 ! Read in JPL spectral data set
1028 !ec032504 CALL RD_TJPL(iph,path_fastj_jvspec)
1029 call rd_tjpl2
1030 !
1031 ! Read in T & O3 climatology
1032 !ec032504 CALL RD_PROF(iph,path_fastj_jvatms)
1033 ! call rd_prof2
1034 !
1035 ! Select Aerosol/Cloud types to be used
1036 ! call set_aer2
1037
1038 return
1039 end subroutine inphot2
1040 !*************************************************************************
1041
1042 subroutine photoj(isvode,jsvode,zpj,timej,nslat,nslon,iozone,tauaer_550_1, &
1043 my_ozone,p,t,od,sa,lpar,jpnl,xgrd,ygrd,tau_fastj,month_fastj,iday_fastj, &
1044 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1045 !-----------------------------------------------------------------------
1046 !----jv_trop.f: new FAST J-Value code, troposphere only (mjprather 6/96)
1047 !---- uses special wavelength quadrature spectral data (jv_spec.dat)
1048 !--- that includes only 289 nm - 800 nm (later a single 205 nm add-on)
1049 !--- uses special compact Mie code based on Feautrier/Auer/Prather vers.
1050 !-----------------------------------------------------------------------
1051 !
1052 ! zpj External array providing J-values to main CTM code
1053 ! timej Offset in hours from start of timestep to time J-values
1054 ! required for - take as half timestep for mid-step Js.
1055 ! solf Solar distance factor, for scaling; normally given by:
1056 ! 1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.))
1057 !
1058 !----------basic common blocks:-----------------------------------------
1059
1060 USE module_data_mosaic_other, only : kmaxd
1061 USE module_fastj_data
1062
1063 IMPLICIT NONE
1064 !jdf
1065 ! Print Fast-J diagnostics if iprint /= 0
1066 integer, parameter :: iprint = 0
1067 integer, parameter :: single = 4 !compiler dependent value real*4
1068 ! integer, parameter :: double = 8 !compiler dependent value real*8
1069 integer,parameter :: ipar_fastj=1,jpar=1
1070 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1071 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1072 integer lpar !Number of levels in CTM
1073 integer jpnl !Number of levels requiring chemistry
1074 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1075 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1076 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1077 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1078 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1079 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1080 integer month_fastj ! Number of month (1-12)
1081 integer iday_fastj ! Day of year
1082 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
1083 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1084 real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios
1085 real tauaer_550_1
1086 integer iozone
1087 integer nslat ! Latitude of current profile point
1088 integer nslon ! Longitude of current profile point
1089 integer nspint ! Num of spectral intervals across solar
1090 parameter ( nspint = 4 ) ! spectrum for FAST-J
1091 real, dimension (nspint),save :: wavmid !cm
1092 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
1093 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
1094 data wavmid &
1095 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
1096 !jdf
1097 real*8 zpj(lpar,jppj),timej,solf
1098 real*8 pi_fastj
1099
1100 integer i,j
1101 integer isvode,jsvode
1102
1103 !-----------------------------------------------------------------------
1104 !
1105 do i=1,jpnl
1106 do j=1,jppj
1107 zj(i,j)=0.D0
1108 zpj(i,j)=0.D0
1109 enddo
1110 enddo
1111 !
1112 !---Calculate new solar zenith angle
1113 CALL SOLAR2(timej,nslat,nslon, &
1114 xgrd,ygrd,tau_fastj,month_fastj,iday_fastj)
1115 if(SZA.gt.szamax) go to 10
1116 !
1117 !---Set up profiles on model levels
1118 CALL SET_PROF(isvode,jsvode,nslat,nslon,iozone,tauaer_550_1, &
1119 my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd)
1120 !
1121 !---Print out atmosphere
1122 if(iprint.ne.0)CALL PRTATM(3,nslat,nslon,tau_fastj,month_fastj,ydgrd) ! code change jcb
1123 ! call prtatm(0)
1124 !
1125 !-----------------------------------------------------------------------
1126 CALL JVALUE(isvode,jsvode,lpar,jpnl, &
1127 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1128 !-----------------------------------------------------------------------
1129 !---Print solar flux terms
1130 ! WRITE(6,'(A16,I5,20I9)') ' wave (beg/end)',(i,i=1,jpnl)
1131 ! DO j=NW1,NW2
1132 ! WRITE(6,'(2F8.2,20F9.6)') WBIN(j),WBIN(j+1),
1133 ! $ (FFF(j,i)/FL(j),i=1,jpnl)
1134 ! ENDDO
1135 !
1136 !---Include variation in distance to sun
1137 pi_fastj=3.1415926536d0
1138 solf=1.d0-(0.034d0*cos(dble(iday_fastj-186)*2.d0 &
1139 *pi_fastj/365.d0))
1140 if(iprint.ne.0)then
1141 ! write(6,'('' solf = '', f10.5)')solf
1142 write(*,'('' solf = '', f10.5)')solf
1143 endif
1144 ! solf=1.d0 ! code change jcb
1145 !-----------------------------------------------------------------------
1146 CALL JRATET(solf,nslat,nslon,p,t,od,sa,lpar,jpnl)
1147 !-----------------------------------------------------------------------
1148 !
1149
1150 ! "zj" updated in JRATET - pass this back to ASAD as "zpj"
1151 do i=1,jpnl
1152 do j=1,jppj
1153 zpj(i,j)= zj(i,j)
1154 enddo
1155 enddo
1156
1157 !
1158 !---Output selected values
1159 10 if((.not.ldeg45.and.nslon.eq.37.and.nslat.eq.36).or. &
1160 (ldeg45.and.nslon.eq.19.and.nslat.eq.18)) then
1161 i=min(jppj,8)
1162 ! write(6,1000)iday_fastj,tau_fastj+timej,sza,jlabel(i),zpj(1,i)
1163 endif
1164 !
1165 return
1166 ! 1000 format(' Photolysis on day ',i4,' at ',f4.1,' hrs: SZA = ',f7.3, &
1167 ! ' J',a7,'= ',1pE10.3)
1168 end subroutine photoj
1169
1170 !*************************************************************************
1171 subroutine set_prof(isvode,jsvode,nslat,nslon,iozone,tauaer_550, &
1172 my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd)
1173 !-----------------------------------------------------------------------
1174 ! Routine to set up atmospheric profiles required by Fast-J using a
1175 ! doubled version of the level scheme used in the CTM. First pressure
1176 ! and z* altitude are defined, then O3 and T are taken from the supplied
1177 ! climatology and integrated to the CTM levels (may be overwritten with
1178 ! values directly from the CTM, if desired) and then black carbon and
1179 ! aerosol profiles are constructed.
1180 ! Oliver (04/07/99)
1181 !-----------------------------------------------------------------------
1182 !
1183 ! pj Pressure at boundaries of model levels (hPa)
1184 ! z Altitude of boundaries of model levels (cm)
1185 ! odcol Optical depth at each model level
1186 ! masfac Conversion factor for pressure to column density
1187 !
1188 ! TJ Temperature profile on model grid
1189 ! DM Air column for each model level (molecules.cm-2)
1190 ! DO3 Ozone column for each model level (molecules.cm-2)
1191 ! DBC Mass of Black Carbon at each model level (g.cm-3) ! .....!
1192 ! PSTD Approximate pressures of levels for supplied climatology
1193 !
1194 !-----------------------------------------------------------------------
1195
1196 USE module_data_mosaic_other, only : kmaxd
1197 USE module_fastj_data
1198
1199 IMPLICIT NONE
1200 !jdf
1201 ! Print Fast-J diagnostics if iprint /= 0
1202 integer, parameter :: iprint = 0
1203 integer, parameter :: single = 4 !compiler dependent value real*4
1204 ! integer, parameter :: double = 8 !compiler dependent value real*8
1205 integer,parameter :: ipar_fastj=1,jpar=1
1206 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1207 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1208 integer lpar !Number of levels in CTM
1209 integer jpnl !Number of levels requiring chemistry
1210 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1211 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1212 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1213 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1214 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1215 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1216 integer month_fastj ! Number of month (1-12)
1217 integer iday_fastj ! Day of year
1218 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
1219 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1220 real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios
1221 real tauaer_550
1222 integer iozone
1223 integer nslat ! Latitude of current profile point
1224 integer nslon ! Longitude of current profile point
1225 !jdf
1226 real*8 pstd(52),oref2(51),tref2(51),bref2(51)
1227 real*8 odcol(lpar),dlogp,f0,t0,b0,pb,pc,xc,masfac,scaleh
1228 real vis, aerd1, aerd2
1229
1230 integer i, k, l, m
1231 integer isvode,jsvode
1232
1233 pj(NB+1) = 0.d0 ! define top level
1234
1235 !
1236 ! Set up cloud and surface properties
1237 call CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa,lpar,jpnl)
1238
1239 ! Mass factor - delta-Pressure (mbars) to delta-Column (molecules.cm-2)
1240 masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0)
1241 !
1242 ! Set up pressure levels for O3/T climatology - assume that value
1243 ! given for each 2 km z* level applies from 1 km below to 1 km above,
1244 ! so select pressures at these boundaries. Surface level values at
1245 ! 1000 mb are assumed to extend down to the actual P(nslon,nslat).
1246 !
1247 pstd(1) = max(pj(1),1000.d0)
1248 pstd(2) = 1000.d0*10.d0**(-1.d0/16.d0)
1249 dlogp = 10.d0**(-2.d0/16.d0)
1250 do i=3,51
1251 pstd(i) = pstd(i-1)*dlogp
1252 enddo
1253 pstd(52) = 0.d0
1254 !
1255 ! Select appropriate monthly and latitudinal profiles
1256 m = max(1,min(12,month_fastj))
1257 l = max(1,min(18,(int(ydgrd(nslat))+99)/10))
1258 !
1259 ! Temporary arrays for climatology data
1260 do i=1,51
1261 oref2(i)=oref(i,l,m)
1262 tref2(i)=tref(i,l,m)
1263 bref2(i)=bref(i)
1264 enddo
1265 !
1266 ! Apportion O3 and T on supplied climatology z* levels onto CTM levels
1267 ! with mass (pressure) weighting, assuming constant mixing ratio and
1268 ! temperature half a layer on either side of the point supplied.
1269 !
1270 do i = 1,NB
1271 F0 = 0.d0
1272 T0 = 0.d0
1273 B0 = 0.d0
1274 do k = 1,51
1275 PC = min(pj(i),pstd(k))
1276 PB = max(pj(i+1),pstd(k+1))
1277 if(PC.gt.PB) then
1278 XC = (PC-PB)/(pj(i)-pj(i+1))
1279 F0 = F0 + oref2(k)*XC
1280 T0 = T0 + tref2(k)*XC
1281 B0 = B0 + bref2(k)*XC
1282 endif
1283 enddo
1284 TJ(i) = T0
1285 DO3(i)= F0*1.d-6
1286 DBC(i)= B0
1287 enddo
1288 !
1289 ! Insert model values here to replace or supplement climatology.
1290 ! Note that CTM temperature is always used in x-section calculations
1291 ! (see JRATET); TJ is used in actinic flux calculation only.
1292 !
1293 do i=1,lpar ! JCB code change; just use climatlogy for upper levels
1294 if(i.le.iozone)DO3(i) = my_ozone(i) ! Volume Mixing Ratio
1295 ! TJ(i) = T(nslon,nslat,I) ! Kelvin
1296 ! JCB - overwrite climatology
1297 ! TJ(i) = (T(nslon,nslat,i)+T(nslon,nslat,i+1))/2. ! JCB - take midpoint
1298 ! code change - now take temperature as appropriate for midpoint of layer
1299 TJ(i)=T(nslon,nslat,i)
1300 enddo
1301 if(lpar+1.le.iozone)then
1302 DO3(lpar+1) = my_ozone(lpar+1) ! Above top of model (or use climatology)
1303 endif
1304 ! TJ(lpar+1) = my_temp(lpar) ! Above top of model (or use climatology)
1305 !wig 26-Aug-2000: Comment out following line so that climatology is used for
1306 ! above the model top.
1307 ! TJ(lpar+1) = T(nslon,nslat,NB) ! JCB - just use climatology or given temperature
1308 ! JCB read in O3
1309 !
1310 !
1311 ! Calculate effective altitudes using scale height at each level
1312 z(1) = 0.d0
1313 do i=1,lpar
1314 scaleh=1.3806d-19*masfac*TJ(i)
1315 z(i+1) = z(i)-(log(pj(i+1)/pj(i))*scaleh)
1316 enddo
1317 !
1318 ! Add Aerosol Column - include aerosol types here. Currently use soot
1319 ! water and ice; assume black carbon x-section of 10 m2/g, independent
1320 ! of wavelength; assume limiting temperature for ice of -40 deg C.
1321
1322 do i=1,lpar
1323 ! AER(1,i) = DBC(i)*10.d0*(z(i+1)-z(i)) ! DBC must be g/m^3
1324 ! calculate AER(1,i) according to aerosol density - use trap rule
1325 vis=23.0
1326 call aeroden(z(i)/100000.,vis,aerd1) ! convert cm to km
1327 call aeroden(z(i+1)/100000.,vis,aerd2)
1328 ! trap rule used here; convert cm to km; divide by 100000.
1329 AER(1,i)=(z(i+1)-z(i))/100000.*(aerd1+aerd2)/2./4287.55*tauaer_550
1330 ! write(6,*)i,z(i)/100000.,aerd1,aerd2,tauaer_550,AER(1,i)
1331 !
1332 if(T(nslon,nslat,I).gt.233.d0) then
1333 AER(2,i) = odcol(i)
1334 AER(3,i) = 0.d0
1335 else
1336 AER(2,i) = 0.d0
1337 AER(3,i) = odcol(i)
1338 endif
1339 enddo
1340 do k=1,MX
1341 AER(k,lpar+1) = 0.d0 ! just set equal to zero
1342 enddo
1343
1344 AER(1,lpar+1)=2.0*AER(1,lpar) ! kludge
1345 !
1346 ! Calculate column quantities for Fast-J
1347 do i=1,NB
1348 DM(i) = (PJ(i)-PJ(i+1))*masfac
1349 DO3(i) = DO3(i)*DM(i)
1350 enddo
1351 !
1352 end subroutine set_prof
1353
1354 !******************************************************************
1355 ! SUBROUTINE CLDSRF(odcol)
1356 SUBROUTINE CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa, &
1357 lpar,jpnl)
1358 !-----------------------------------------------------------------------
1359 ! Routine to set cloud and surface properties
1360 !-----------------------------------------------------------------------
1361 ! rflect Surface albedo (Lambertian)
1362 ! odmax Maximum allowed optical depth, above which they are scaled
1363 ! odcol Optical depth at each model level
1364 ! odsum Column optical depth
1365 ! nlbatm Level of lower photolysis boundary - usually surface
1366 !-----------------------------------------------------------------------
1367
1368 USE module_data_mosaic_other, only : kmaxd
1369 USE module_fastj_data
1370
1371 IMPLICIT NONE
1372 !jdf
1373 ! Print Fast-J diagnostics if iprint /= 0
1374 integer, parameter :: iprint = 0
1375 integer, parameter :: single = 4 !compiler dependent value real*4
1376 ! integer, parameter :: double = 8 !compiler dependent value real*8
1377 integer,parameter :: ipar_fastj=1,jpar=1
1378 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1379 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1380 integer lpar !Number of levels in CTM
1381 integer jpnl !Number of levels requiring chemistry
1382 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1383 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1384 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1385 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1386 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1387 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1388 integer month_fastj ! Number of month (1-12)
1389 integer iday_fastj ! Day of year
1390 integer nslat ! Latitude of current profile point
1391 integer nslon ! Longitude of current profile point
1392 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
1393 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1394 !jdf
1395 integer i, j, k
1396 integer isvode, jsvode
1397 real*8 odcol(lpar), odsum, odmax, odtot
1398 !
1399 ! Default lower photolysis boundary as bottom of level 1
1400 nlbatm = 1
1401 !
1402 ! Set surface albedo
1403 RFLECT = dble(SA(nslon,nslat))
1404 RFLECT = max(0.d0,min(1.d0,RFLECT))
1405 !
1406 ! Zero aerosol column
1407 do k=1,MX
1408 do i=1,NB
1409 AER(k,i) = 0.d0
1410 enddo
1411 enddo
1412 !
1413 ! Scale optical depths as appropriate - limit column to 'odmax'
1414 odmax = 200.d0
1415 odsum = 0.d0
1416 do i=1,lpar
1417 odcol(i) = dble(OD(nslon,nslat,i))
1418 odsum = odsum + odcol(i)
1419 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf odcol',odcol(i),odsum
1420 enddo
1421 if(odsum.gt.odmax) then
1422 odsum = odmax/odsum
1423 do i=1,lpar
1424 odcol(i) = odcol(i)*odsum
1425 enddo
1426 odsum = odmax
1427 endif
1428 ! Set sub-division switch if appropriate
1429 odtot=0.d0
1430 jadsub(nb)=0
1431 jadsub(nb-1)=0
1432 do i=nb-1,1,-1
1433 k=2*i
1434 jadsub(k)=0
1435 jadsub(k-1)=0
1436 odtot=odtot+odcol(i)
1437 if(odcol(i).gt.0.d0.and.dtausub.gt.0.d0) then
1438 if(odtot.le.dtausub) then
1439 jadsub(k)=1
1440 jadsub(k-1)=1
1441 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf1',i,k,jadsub(k)
1442 elseif(odtot.gt.dtausub) then
1443 jadsub(k)=1
1444 jadsub(k-1)=0
1445 do j=1,2*(i-1)
1446 jadsub(j)=0
1447 enddo
1448 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf2',i,k,jadsub(k)
1449 go to 20
1450 endif
1451 endif
1452 enddo
1453 20 continue
1454 !
1455 return
1456 end SUBROUTINE CLDSRF
1457
1458 !********************************************************************
1459 subroutine solar2(timej,nslat,nslon, &
1460 xgrd,ygrd,tau_fastj,month_fastj,iday_fastj)
1461 !-----------------------------------------------------------------------
1462 ! Routine to set up SZA for given lat, lon and time
1463 !-----------------------------------------------------------------------
1464 ! timej Offset in hours from start of timestep to time J-values
1465 ! required for - take as half timestep for mid-step Js.
1466 !-----------------------------------------------------------------------
1467
1468 USE module_data_mosaic_other, only : kmaxd
1469 USE module_fastj_data
1470 IMPLICIT NONE
1471 !jdf
1472 ! Print Fast-J diagnostics if iprint /= 0
1473 integer, parameter :: iprint = 0
1474 integer, parameter :: single = 4 !compiler dependent value real*4
1475 ! integer, parameter :: double = 8 !compiler dependent value real*8
1476 integer,parameter :: ipar_fastj=1,jpar=1
1477 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1478 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1479 integer lpar !Number of levels in CTM
1480 integer jpnl !Number of levels requiring chemistry
1481 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1482 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1483 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1484 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1485 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1486 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1487 integer month_fastj ! Number of month (1-12)
1488 integer iday_fastj ! Day of year
1489 integer nslat ! Latitude of current profile point
1490 integer nslon ! Longitude of current profile point
1491 !jdf
1492 real*8 pi_fastj, pi180, loct, timej
1493 real*8 sindec, soldek, cosdec, sinlat, sollat, coslat, cosz
1494 !
1495 pi_fastj=3.141592653589793D0
1496 pi180=pi_fastj/180.d0
1497 sindec=0.3978d0*sin(0.9863d0*(dble(iday_fastj)-80.d0)*pi180)
1498 soldek=asin(sindec)
1499 cosdec=cos(soldek)
1500 sinlat=sin(ygrd(nslat))
1501 sollat=asin(sinlat)
1502 coslat=cos(sollat)
1503 !
1504 loct = (((tau_fastj+timej)*15.d0)-180.d0)*pi180 + xgrd(nslon)
1505 cosz = cosdec*coslat*cos(loct) + sindec*sinlat
1506 sza = acos(cosz)/pi180
1507 U0 = cos(SZA*pi180)
1508 !
1509 return
1510 end subroutine solar2
1511
1512 !**********************************************************************
1513
1514
1515 SUBROUTINE JRATET(SOLF,nslat,nslon,p,t,od,sa,lpar,jpnl)
1516 !-----------------------------------------------------------------------
1517 ! Calculate and print J-values. Note that the loop in this routine
1518 ! only covers the jpnl levels actually needed by the CTM.
1519 !-----------------------------------------------------------------------
1520 !
1521 ! FFF Actinic flux at each level for each wavelength bin
1522 ! QQQ Cross sections for species (read in in RD_TJPL)
1523 ! SOLF Solar distance factor, for scaling; normally given by:
1524 ! 1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.))
1525 ! Assumes aphelion day 186, perihelion day 3.
1526 ! TQQ Temperatures at which QQQ cross sections supplied
1527 !
1528 !-----------------------------------------------------------------------
1529
1530 USE module_data_mosaic_other, only : kmaxd
1531 USE module_fastj_data
1532
1533 IMPLICIT NONE
1534 !jdf
1535 ! Print Fast-J diagnostics if iprint /= 0
1536 integer, parameter :: iprint = 0
1537 integer, parameter :: single = 4 !compiler dependent value real*4
1538 ! integer, parameter :: double = 8 !compiler dependent value real*8
1539 integer,parameter :: ipar_fastj=1,jpar=1
1540 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1541 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1542 integer lpar !Number of levels in CTM
1543 integer jpnl !Number of levels requiring chemistry
1544 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1545 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1546 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1547 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1548 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1549 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1550 integer month_fastj ! Number of month (1-12)
1551 integer iday_fastj ! Day of year
1552 integer nslat ! Latitude of current profile point
1553 integer nslon ! Longitude of current profile point
1554 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
1555 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1556 !jdf
1557 integer i, j, k
1558 real*8 qo2tot, qo3tot, qo31d, qo33p, qqqt
1559 ! real*8 xseco2, xseco3, xsec1d, solf, tfact
1560 real*8 solf, tfact
1561 !
1562 do I=1,jpnl
1563 VALJ(1) = 0.d0
1564 VALJ(2) = 0.d0
1565 VALJ(3) = 0.d0
1566 do K=NW1,NW2 ! Using model 'T's here
1567 QO2TOT= xseco2(K,dble(T(nslon,nslat,I)))
1568 VALJ(1) = VALJ(1) + QO2TOT*FFF(K,I)
1569 QO3TOT= xseco3(K,dble(T(nslon,nslat,I)))
1570 QO31D = xsec1d(K,dble(T(nslon,nslat,I)))*QO3TOT
1571 QO33P = QO3TOT - QO31D
1572 VALJ(2) = VALJ(2) + QO33P*FFF(K,I)
1573 VALJ(3) = VALJ(3) + QO31D*FFF(K,I)
1574 enddo
1575 !------Calculate remaining J-values with T-dep X-sections
1576 do J=4,NJVAL
1577 VALJ(J) = 0.d0
1578 TFACT = 0.d0
1579 if(TQQ(2,J).gt.TQQ(1,J)) TFACT = max(0.d0,min(1.d0, &
1580 (T(nslon,nslat,I)-TQQ(1,J))/(TQQ(2,J)-TQQ(1,J)) ))
1581 do K=NW1,NW2
1582 QQQT = QQQ(K,1,J-3) + (QQQ(K,2,J-3) - QQQ(K,1,J-3))*TFACT
1583 VALJ(J) = VALJ(J) + QQQT*FFF(K,I)
1584 !------Additional code for pressure dependencies
1585 ! if(jpdep(J).ne.0) then
1586 ! VALJ(J) = VALJ(J) + QQQT*FFF(K,I)*
1587 ! $ (zpdep(K,L)*(pj(i)+pj(i+1))*0.5d0)
1588 ! endif
1589 enddo
1590 enddo
1591 do j=1,jppj
1592 zj(i,j)=VALJ(jind(j))*jfacta(j)*SOLF
1593 enddo
1594 ! Herzberg bin
1595 do j=1,nhz
1596 zj(i,hzind(j))=hztoa(j)*fhz(i)*SOLF
1597 enddo
1598 enddo
1599 return
1600 end SUBROUTINE JRATET
1601
1602 !*********************************************************************
1603
1604
1605 SUBROUTINE PRTATM(N,nslat,nslon,tau_fastj,month_fastj,ydgrd)
1606 !-----------------------------------------------------------------------
1607 ! Print out the atmosphere and calculate appropriate columns
1608 ! N=1 Print out column totals only
1609 ! N=2 Print out full columns
1610 ! N=3 Print out full columns and climatology
1611 !-----------------------------------------------------------------------
1612
1613 USE module_data_mosaic_other, only : kmaxd
1614 USE module_fastj_data
1615
1616 !jdf
1617 ! Print Fast-J diagnostics if iprint /= 0
1618 integer, parameter :: iprint = 0
1619 integer, parameter :: single = 4 !compiler dependent value real*4
1620 ! integer, parameter :: double = 8 !compiler dependent value real*8
1621 integer,parameter :: ipar_fastj=1,jpar=1
1622 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1623 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1624 integer lpar !Number of levels in CTM
1625 integer jpnl !Number of levels requiring chemistry
1626 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1627 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1628 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1629 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1630 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1631 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1632 integer month_fastj ! Number of month (1-12)
1633 integer iday_fastj ! Day of year
1634 integer nslat ! Latitude of current profile point
1635 integer nslon ! Longitude of current profile point
1636 !jdf
1637 integer n, i, k, l, m
1638 real*8 COLO3(NB),COLO2(NB),COLAX(MX,NB),ZKM,ZSTAR,PJC
1639 real*8 climat(9),masfac,dlogp
1640 if(N.eq.0) return
1641 !---Calculate columns, for diagnostic output only:
1642 COLO3(NB) = DO3(NB)
1643 COLO2(NB) = DM(NB)*0.20948d0
1644 do K=1,MX
1645 COLAX(K,NB) = AER(K,NB)
1646 enddo
1647 do I=NB-1,1,-1
1648 COLO3(i) = COLO3(i+1)+DO3(i)
1649 COLO2(i) = COLO2(i+1)+DM(i)*0.20948d0
1650 do K=1,MX
1651 COLAX(k,i) = COLAX(k,i+1)+AER(k,i)
1652 enddo
1653 enddo
1654 write(*,1200) ' Tau=',tau_fastj,' SZA=',sza
1655 write(*,1200) ' O3-column(DU)=',COLO3(1)/2.687d16, &
1656 ' column aerosol @1000nm=',(COLAX(K,1),K=1,MX)
1657 !---Print out atmosphere
1658 if(N.gt.1) then
1659 write(*,1000) (' AER-X ','col-AER',k=1,mx)
1660 do I=NB,1,-1
1661 PJC = PJ(I)
1662 ZKM =1.d-5*Z(I)
1663 ZSTAR = 16.d0*DLOG10(1013.d0/PJC)
1664 write(*,1100) I,ZKM,ZSTAR,DM(I),DO3(I),1.d6*DO3(I)/DM(I), &
1665 TJ(I),PJC,COLO3(I),COLO2(I),(AER(K,I),COLAX(K,I),K=1,MX)
1666 enddo
1667 endif
1668 !
1669 !---Print out climatology
1670 if(N.gt.2) then
1671 do i=1,9
1672 climat(i)=0.d0
1673 enddo
1674 m = max(1,min(12,month_fastj))
1675 l = max(1,min(18,(int(ydgrd(nslat))+99)/10))
1676 masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0)
1677 write(*,*) 'Specified Climatology'
1678 write(*,1000)
1679 do i=51,1,-1
1680 dlogp=10.d0**(-1.d0/16.d0)
1681 PJC = 1000.d0*dlogp**(2*i-2)
1682 climat(1) = 16.d0*DLOG10(1000.D0/PJC)
1683 climat(2) = climat(1)
1684 climat(3) = PJC*(1.d0/dlogp-dlogp)*masfac
1685 if(i.eq.1) climat(3)=PJC*(1.d0-dlogp)*masfac
1686 climat(4)=climat(3)*oref(i,l,m)*1.d-6
1687 climat(5)=oref(i,l,m)
1688 climat(6)=tref(i,l,m)
1689 climat(7)=PJC
1690 climat(8)=climat(8)+climat(4)
1691 climat(9)=climat(9)+climat(3)*0.20948d0
1692 write(*,1100) I,(climat(k),k=1,9)
1693 enddo
1694 write(*,1200) ' O3-column(DU)=',climat(8)/2.687d16
1695 endif
1696 return
1697 1000 format(5X,'Zkm',3X,'Z*',8X,'M',8X,'O3',6X,'f-O3',5X,'T',7X,'P',6x, &
1698 'col-O3',3X,'col-O2',2X,10(a7,2x))
1699 1100 format(1X,I2,0P,2F6.2,1P,2E10.3,0P,F7.3,F8.2,F10.4,1P,10E9.2)
1700 1200 format(A,F8.1,A,10(1pE10.3))
1701 end SUBROUTINE PRTATM
1702
1703 SUBROUTINE JVALUE(isvode,jsvode,lpar,jpnl, &
1704 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1705 !-----------------------------------------------------------------------
1706 ! Calculate the actinic flux at each level for the current SZA value.
1707 ! quit when SZA > 98.0 deg ==> tangent height = 63 km
1708 ! or 99. 80 km
1709 !-----------------------------------------------------------------------
1710 !
1711 ! AVGF Attenuation of beam at each level for each wavelength
1712 ! FFF Actinic flux at each desired level
1713 ! FHZ Actinic flux in Herzberg bin
1714 ! WAVE Effective wavelength of each wavelength bin
1715 ! XQO2 Absorption cross-section of O2
1716 ! XQO3 Absorption cross-section of O3
1717 !
1718 !-----------------------------------------------------------------------
1719
1720 USE module_data_mosaic_other, only : kmaxd
1721 USE module_fastj_data
1722 USE module_peg_util, only: peg_message
1723 !jdf
1724 ! Print Fast-J diagnostics if iprint /= 0
1725 integer, parameter :: iprint = 0
1726 integer, parameter :: single = 4 !compiler dependent value real*4
1727 ! integer, parameter :: double = 8 !compiler dependent value real*8
1728 integer,parameter :: ipar_fastj=1,jpar=1
1729 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
1730 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
1731 integer lpar !Number of levels in CTM
1732 integer jpnl !Number of levels requiring chemistry
1733 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
1734 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
1735 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
1736 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
1737 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
1738 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
1739 integer month_fastj ! Number of month (1-12)
1740 integer iday_fastj ! Day of year
1741 real(kind=double), dimension(ipar_fastj,jpar) :: P , SA
1742 real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD
1743 integer nspint ! Num of spectral intervals across solar
1744 parameter ( nspint = 4 ) ! spectrum for FAST-J
1745 real, dimension (nspint),save :: wavmid !cm
1746 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
1747 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
1748 data wavmid &
1749 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
1750 !jdf
1751 integer j, k
1752 ! real*8 wave, xseco3, xseco2
1753 real*8 wave
1754 real*8 AVGF(lpar),XQO3(NB),XQO2(NB)
1755 ! diagnostics for error situations
1756 ! integer lunout
1757 ! parameter (lunout = 41)
1758 integer isvode,jsvode
1759 character*80 msg
1760 !
1761 do J=1,jpnl
1762 do K=NW1,NW2
1763 FFF(K,J) = 0.d0
1764 enddo
1765 FHZ(J) = 0.d0
1766 enddo
1767 !
1768 !---SZA check
1769 if(SZA.gt.szamax) GOTO 99
1770 !
1771 !---Calculate spherical weighting functions
1772 CALL SPHERE(lpar)
1773 !
1774 !---Loop over all wavelength bins
1775 do K=NW1,NW2
1776 WAVE = WL(K)
1777 do J=1,NB
1778 XQO3(J) = xseco3(K,dble(TJ(J)))
1779 enddo
1780 do J=1,NB
1781 XQO2(J) = xseco2(K,dble(TJ(J)))
1782 enddo
1783 !-----------------------------------------
1784 CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, &
1785 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1786 !-----------------------------------------
1787 do J=1,jpnl
1788 FFF(K,J) = FFF(K,J) + FL(K)*AVGF(J)
1789 ! diagnostic
1790 if ( FFF(K,J) .lt. 0) then
1791 write( msg, '(a,2i4,e14.6)' ) &
1792 'FASTJ neg actinic flux ' // &
1793 'k j FFF(K,J) ', k, j, fff(k,j)
1794 call peg_message( lunerr, msg )
1795 end if
1796 ! end diagnostic
1797 enddo
1798 enddo
1799 !
1800 !---Herzberg continuum bin above 10 km, if required
1801 if(NHZ.gt.0) then
1802 K=NW2+1
1803 WAVE = 204.d0
1804 do J=1,NB
1805 XQO3(J) = HZO3
1806 XQO2(J) = HZO2
1807 enddo
1808 CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, &
1809 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1810 do J=1,jpnl
1811 if(z(j).gt.1.d6) FHZ(J)=AVGF(J)
1812 enddo
1813 endif
1814 !
1815 99 continue
1816 1000 format(' SZA=',f6.1,' Reflectvty=',f6.3,' OD=',10(1pe10.3))
1817
1818 return
1819 end SUBROUTINE JVALUE
1820
1821 FUNCTION xseco3(K,TTT)
1822 !-----------------------------------------------------------------------
1823 ! Cross-sections for O3 for all processes interpolated across 3 temps
1824 !-----------------------------------------------------------------------
1825
1826 USE module_fastj_data
1827
1828 integer k
1829 ! real*8 ttt, flint, xseco3
1830 real*8 ttt, xseco3
1831 xseco3 = &
1832 flint(TTT,TQQ(1,2),TQQ(2,2),TQQ(3,2),QO3(K,1),QO3(K,2),QO3(K,3))
1833 return
1834 end FUNCTION xseco3
1835
1836 FUNCTION xsec1d(K,TTT)
1837 !-----------------------------------------------------------------------
1838 ! Quantum yields for O3 --> O2 + O(1D) interpolated across 3 temps
1839 !-----------------------------------------------------------------------
1840
1841 USE module_fastj_data
1842
1843 integer k
1844 ! real*8 ttt, flint, xsec1d
1845 real*8 ttt, xsec1d
1846 xsec1d = &
1847 flint(TTT,TQQ(1,3),TQQ(2,3),TQQ(3,3),Q1D(K,1),Q1D(K,2),Q1D(K,3))
1848 return
1849 end FUNCTION xsec1d
1850
1851 FUNCTION xseco2(K,TTT)
1852 !-----------------------------------------------------------------------
1853 ! Cross-sections for O2 interpolated across 3 temps; No S_R Bands yet!
1854 !-----------------------------------------------------------------------
1855
1856 USE module_fastj_data
1857
1858 integer k
1859 ! real*8 ttt, flint, xseco2
1860 real*8 ttt, xseco2
1861 xseco2 = &
1862 flint(TTT,TQQ(1,1),TQQ(2,1),TQQ(3,1),QO2(K,1),QO2(K,2),QO2(K,3))
1863 return
1864 end FUNCTION xseco2
1865
1866 REAL*8 FUNCTION flint (TINT,T1,T2,T3,F1,F2,F3)
1867 !-----------------------------------------------------------------------
1868 ! Three-point linear interpolation function
1869 !-----------------------------------------------------------------------
1870 real*8 TINT,T1,T2,T3,F1,F2,F3
1871 IF (TINT .LE. T2) THEN
1872 IF (TINT .LE. T1) THEN
1873 flint = F1
1874 ELSE
1875 flint = F1 + (F2 - F1)*(TINT -T1)/(T2 -T1)
1876 ENDIF
1877 ELSE
1878 IF (TINT .GE. T3) THEN
1879 flint = F3
1880 ELSE
1881 flint = F2 + (F3 - F2)*(TINT -T2)/(T3 -T2)
1882 ENDIF
1883 ENDIF
1884 return
1885 end FUNCTION flint
1886
1887 SUBROUTINE SPHERE(lpar)
1888 !-----------------------------------------------------------------------
1889 ! Calculation of spherical geometry; derive tangent heights, slant path
1890 ! lengths and air mass factor for each layer. Not called when
1891 ! SZA > 98 degrees. Beyond 90 degrees, include treatment of emergent
1892 ! beam (where tangent height is below altitude J-value desired at).
1893 !-----------------------------------------------------------------------
1894 !
1895 ! GMU MU, cos(solar zenith angle)
1896 ! RZ Distance from centre of Earth to each point (cm)
1897 ! RQ Square of radius ratios
1898 ! TANHT Tangent height for the current SZA
1899 ! XL Slant path between points
1900 ! AMF Air mass factor for slab between level and level above
1901 !
1902 !-----------------------------------------------------------------------
1903
1904 USE module_fastj_data
1905
1906 !jdf
1907 integer lpar
1908 !jdf
1909 integer i, j, k, ii
1910 real*8 airmas, gmu, xmu1, xmu2, xl, diff
1911 REAL*8 Ux,H,RZ(NB),RQ(NB),ZBYR
1912 !
1913 ! Inlined air mass factor function for top of atmosphere
1914 AIRMAS(Ux,H) = (1.0d0+H)/SQRT(Ux*Ux+2.0d0*H*(1.0d0- &
1915 0.6817d0*EXP(-57.3d0*ABS(Ux)/SQRT(1.0d0+5500.d0*H))/ &
1916 (1.0d0+0.625d0*H)))
1917 !
1918 GMU = U0
1919 RZ(1)=RAD+Z(1)
1920 ZBYR = ZZHT/RAD
1921 DO 2 II=2,NB
1922 RZ(II) = RAD + Z(II)
1923 RQ(II-1) = (RZ(II-1)/RZ(II))**2
1924 2 CONTINUE
1925 IF (GMU.LT.0.0D0) THEN
1926 TANHT = RZ(nlbatm)/DSQRT(1.0D0-GMU**2)
1927 ELSE
1928 TANHT = RZ(nlbatm)
1929 ENDIF
1930 !
1931 ! Go up from the surface calculating the slant paths between each level
1932 ! and the level above, and deriving the appropriate Air Mass Factor
1933 DO 16 J=1,NB
1934 DO K=1,NB
1935 AMF(K,J)=0.D0
1936 ENDDO
1937 !
1938 ! Air Mass Factors all zero if below the tangent height
1939 IF (RZ(J).LT.TANHT) GOTO 16
1940 ! Ascend from layer J calculating AMFs
1941 XMU1=ABS(GMU)
1942 DO 12 I=J,lpar
1943 XMU2=DSQRT(1.0D0-RQ(I)*(1.0D0-XMU1**2))
1944 XL=RZ(I+1)*XMU2-RZ(I)*XMU1
1945 AMF(I,J)=XL/(RZ(I+1)-RZ(I))
1946 XMU1=XMU2
1947 12 CONTINUE
1948 ! Use function and scale height to provide AMF above top of model
1949 AMF(NB,J)=AIRMAS(XMU1,ZBYR)
1950 !
1951 ! Twilight case - Emergent Beam
1952 IF (GMU.GE.0.0D0) GOTO 16
1953 XMU1=ABS(GMU)
1954 ! Descend from layer J
1955 DO 14 II=J-1,1,-1
1956 DIFF=RZ(II+1)*DSQRT(1.0D0-XMU1**2)-RZ(II)
1957 if(II.eq.1) DIFF=max(DIFF,0.d0) ! filter
1958 ! Tangent height below current level - beam passes through twice
1959 IF (DIFF.LT.0.0D0) THEN
1960 XMU2=DSQRT(1.0D0-(1.0D0-XMU1**2)/RQ(II))
1961 XL=ABS(RZ(II+1)*XMU1-RZ(II)*XMU2)
1962 AMF(II,J)=2.d0*XL/(RZ(II+1)-RZ(II))
1963 XMU1=XMU2
1964 ! Lowest level intersected by emergent beam
1965 ELSE
1966 XL=RZ(II+1)*XMU1*2.0D0
1967 ! WTING=DIFF/(RZ(II+1)-RZ(II))
1968 ! AMF(II,J)=(1.0D0-WTING)*2.D0**XL/(RZ(II+1)-RZ(II))
1969 AMF(II,J)=XL/(RZ(II+1)-RZ(II))
1970 GOTO 16
1971 ENDIF
1972 14 CONTINUE
1973 !
1974 16 CONTINUE
1975 RETURN
1976 END SUBROUTINE SPHERE
1977
1978
1979 SUBROUTINE OPMIE(isvode,jsvode,KW,WAVEL,XQO2,XQO3,FMEAN,lpar,jpnl, &
1980 ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7)
1981 !-----------------------------------------------------------------------
1982 ! NEW Mie code for J's, only uses 8-term expansion, 4-Gauss pts
1983 ! Currently allow up to NP aerosol phase functions (at all altitudes) to
1984 ! be associated with optical depth AER(1:NC) = aerosol opt.depth @ 1000 nm
1985 !
1986 ! Pick Mie-wavelength with phase function and Qext:
1987 !
1988 ! 01 RAYLE = Rayleigh phase
1989 ! 02 ISOTR = isotropic
1990 ! 03 ABSRB = fully absorbing 'soot', wavelength indep.
1991 ! 04 S_Bkg = backgrnd stratospheric sulfate (n=1.46,log-norm:r=.09um/sigma=.6)
1992 ! 05 S_Vol = volcanic stratospheric sulfate (n=1.46,log-norm:r=.08um/sigma=.8)
1993 ! 06 W_H01 = water haze (H1/Deirm.) (n=1.335, gamma: r-mode=0.1um /alpha=2)
1994 ! 07 W_H04 = water haze (H1/Deirm.) (n=1.335, gamma: r-mode=0.4um /alpha=2)
1995 ! 08 W_C02 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=2.0um /alpha=6)
1996 ! 09 W_C04 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=4.0um /alpha=6)
1997 ! 10 W_C08 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=8.0um /alpha=6)
1998 ! 11 W_C13 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=13.3um /alpha=6)
1999 ! 12 W_L06 = water cloud (Lacis) (n=1.335, r-mode=5.5um / alpha=11/3)
2000 ! 13 Ice-H = hexagonal ice cloud (Mishchenko)
2001 ! 14 Ice-I = irregular ice cloud (Mishchenko)
2002 !
2003 ! Choice of aerosol index MIEDX is made in SET_AER; optical depths are
2004 ! apportioned to the AER array in SET_PROF
2005 !
2006 !-----------------------------------------------------------------------
2007 ! FUNCTION RAYLAY(WAVE)---RAYLEIGH CROSS-SECTION for wave > 170 nm
2008 ! WSQI = 1.E6/(WAVE*WAVE)
2009 ! REFRM1 = 1.0E-6*(64.328+29498.1/(146.-WSQI)+255.4/(41.-WSQI))
2010 ! RAYLAY = 5.40E-21*(REFRM1*WSQI)**2
2011 !-----------------------------------------------------------------------
2012 !
2013 ! DTAUX Local optical depth of each CTM level
2014 ! PIRAY Contribution of Rayleigh scattering to extinction
2015 ! PIAER Contribution of Aerosol scattering to extinction
2016 ! TTAU Optical depth of air vertically above each point (to top of atm)
2017 ! FTAU Attenuation of solar beam
2018 ! POMEGA Scattering phase function
2019 ! FMEAN Mean actinic flux at desired levels
2020 !
2021 !-----------------------------------------------------------------------
2022
2023 USE module_data_mosaic_other, only : kmaxd
2024 USE module_fastj_data
2025 USE module_peg_util, only: peg_message, peg_error_fatal
2026 !jdf
2027 ! Print Fast-J diagnostics if iprint /= 0
2028 integer, parameter :: iprint = 0
2029 integer, parameter :: single = 4 !compiler dependent value real*4
2030 ! integer, parameter :: double = 8 !compiler dependent value real*8
2031 integer,parameter :: ipar_fastj=1,jpar=1
2032 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
2033 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
2034 integer lpar !Number of levels in CTM
2035 integer jpnl !Number of levels requiring chemistry
2036 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
2037 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
2038 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
2039 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
2040 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
2041 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
2042 integer month_fastj ! Number of month (1-12)
2043 integer iday_fastj ! Day of year
2044 integer nspint ! Num of spectral intervals across solar
2045 parameter ( nspint = 4 ) ! spectrum for FAST-J
2046 real, dimension (nspint),save :: wavmid !cm
2047 real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer
2048 real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7
2049 data wavmid &
2050 / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 /
2051 INTEGER NL, N__, M__
2052 PARAMETER (NL=350, N__=2*NL, M__=4)
2053 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2054 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2055 REAL(kind=double), dimension(M__,2*M__) :: PM
2056 REAL(kind=double), dimension(2*M__) :: PM0
2057 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2058 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2059 REAL(kind=double), dimension(M__,M__,N__) :: DD
2060 REAL(kind=double), dimension(M__,N__) :: RR
2061 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2062 INTEGER ND,N,M,MFIT
2063 !jdf
2064 integer jndlev(lpar),jaddlv(nc),jaddto(nc+1)
2065 integer KW,km,i,j,k,l,ix,j1
2066 integer isvode,jsvode
2067 real*8 QXMIE(MX),XLAER(MX),SSALB(MX)
2068 real*8 xlo2,xlo3,xlray,xltau,zk,taudn,tauup,zk2
2069 real*8 WAVEL,XQO2(NB),XQO3(NB),FMEAN(lpar),POMEGAJ(2*M__,NC+1)
2070 real*8 DTAUX(NB),PIRAY(NB),PIAER(MX,NB),TTAU(NC+1),FTAU(NC+1)
2071 real*8 ftaulog,dttau,dpomega(2*M__)
2072 real*8 ftaulog2,dttau2,dpomega2(2*M__)
2073 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2074 real*8 PIAER_MX1(NB)
2075 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2076 character*80 msg
2077 !
2078 !---Pick nearest Mie wavelength, no interpolation--------------
2079 KM=1
2080 if( WAVEL .gt. 355.d0 ) KM=2
2081 if( WAVEL .gt. 500.d0 ) KM=3
2082 ! if( WAVEL .gt. 800.d0 ) KM=4 !drop the 1000 nm wavelength
2083 !
2084 !---For Mie code scale extinction at 1000 nm to wavelength WAVEL (QXMIE)
2085 ! define angstrom exponent
2086 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2087 ang=log(QAA(1,MIEDX(1))/QAA(4,MIEDX(1)))/log(300./999.)
2088 do I=1,MX
2089 ! QAA is extinction efficiency
2090 QXMIE(I) = QAA(KM,MIEDX(I))/QAA(4,MIEDX(I))
2091 ! scale to 550 nm using angstrom relationship
2092 ! note that this gives QXMIE at 550.0 nm = 1.0, aerosol optical thickness
2093 ! is defined at 550 nm
2094 ! convention -- I = 1 is aerosol, I > 1 are clouds
2095 if(I.eq.1) QXMIE(I) = (WAVEL/550.0)**ang
2096 SSALB(I) = SSA(KM,MIEDX(I)) ! single scattering albedo
2097 enddo
2098 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2099 !
2100 !---Reinitialize arrays
2101 do j=1,nc+1
2102 ttau(j)=0.d0
2103 ftau(j)=0.d0
2104 enddo
2105 !
2106 !---Set up total optical depth over each CTM level, DTAUX
2107 J1 = NLBATM
2108 do J=J1,NB
2109 XLO3=DO3(J)*XQO3(J)
2110 XLO2=DM(J)*XQO2(J)*0.20948d0
2111 XLRAY=DM(J)*QRAYL(KW)
2112 ! Zero absorption for testing purposes
2113 ! call NOABS(XLO3,XLO2,XLRAY,AER(1,j),RFLECT)
2114 do I=1,MX
2115 ! I is aerosol type, j is level, AER(I,J)*QXMIE(I) is the layer aerosol optical thickness
2116 ! at 1000 nm , AER(I,J), times extinction efficiency for the layer (normalized to be one at 1000 nm)
2117 ! therefore xlaer(i) is the layer optical depth at the wavelength index KM
2118 XLAER(I)=AER(I,J)*QXMIE(I)
2119 enddo
2120 ! Total optical depth from all elements
2121 DTAUX(J)=XLO3+XLO2+XLRAY
2122 do I=1,MX
2123 DTAUX(J)=DTAUX(J)+XLAER(I)
2124 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf xlaer',&
2125 ! i,j,km,dtaux(j),xlaer(i),aer(i,j),qxmie(i),xlo3,xlo2,xlray
2126 enddo
2127 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2128 ! add in new aerosol information from Mie code
2129 ! layer aerosol optical thickness at wavelength index KM, layer j
2130 ! tauaer(km,j)=0.0
2131 dtaux(j)=dtaux(j)+tauaer(km,j)
2132 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf dtaux',&
2133 ! j,km,dtaux(j),tauaer(km,j)
2134 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2135 ! Fractional extinction for Rayleigh scattering and each aerosol type
2136 PIRAY(J)=XLRAY/DTAUX(J)
2137 do I=1,MX
2138 PIAER(I,J)=SSALB(I)*XLAER(I)/DTAUX(J)
2139 enddo
2140 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2141 ! note the level is now important
2142 PIAER_MX1(J)=waer(km,j)*tauaer(km,j)/DTAUX(J)
2143 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2144 enddo ! of the level "j" loop
2145 !
2146 !---Define the scattering phase fn. with mix of Rayleigh(1) & Mie(MIEDX)
2147 ! No. of quadrature pts fixed at 4 (M__), expansion of phase fn @ 8
2148 N = M__
2149 MFIT = 2*M__
2150 do j=j1,NB ! jcb: layer index
2151 do i=1,MFIT
2152 pomegaj(i,j) = PIRAY(J)*PAA(I,KM,1) ! jcb: paa are the expansion coefficients of the phase function
2153 do k=1,MX ! jcb: mx is # of aerosols
2154 pomegaj(i,j) = pomegaj(i,j) + PIAER(K,J)*PAA(I,KM,MIEDX(K))
2155 enddo
2156 enddo
2157 ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
2158 ! i is the # of coefficients, KM is the wavelength index, j is the level
2159 ! note the level in now relevant because we allow aerosol properties to
2160 ! vary by level
2161 pomegaj(1,j) = pomegaj(1,j) + PIAER_MX1(J)*1.0 ! 1.0 is l0
2162 pomegaj(2,j) = pomegaj(2,j) + PIAER_MX1(J)*gaer(KM,j)*3.0 ! the three converts gear to l1
2163 pomegaj(3,j) = pomegaj(3,j) + PIAER_MX1(J)*l2(KM,j)
2164 pomegaj(4,j) = pomegaj(4,j) + PIAER_MX1(J)*l3(KM,j)
2165 pomegaj(5,j) = pomegaj(5,j) + PIAER_MX1(J)*l4(KM,j)
2166 pomegaj(6,j) = pomegaj(6,j) + PIAER_MX1(J)*l5(KM,j)
2167 pomegaj(7,j) = pomegaj(7,j) + PIAER_MX1(J)*l6(KM,j)
2168 pomegaj(8,j) = pomegaj(7,j) + PIAER_MX1(J)*l7(KM,j)
2169 ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB
2170
2171 enddo
2172 !
2173 !---Calculate attenuated incident beam EXP(-TTAU/U0) and flux on surface
2174 do J=J1,NB
2175 if(AMF(J,J).gt.0.0D0) then
2176 XLTAU=0.0D0
2177 do I=1,NB
2178 XLTAU=XLTAU + DTAUX(I)*AMF(I,J)
2179 enddo
2180 if(XLTAU.gt.450.d0) then ! for compilers with no underflow trapping
2181 FTAU(j)=0.d0
2182 else
2183 FTAU(J)=DEXP(-XLTAU)
2184 endif
2185 else
2186 FTAU(J)=0.0D0
2187 endif
2188 enddo
2189 if(U0.gt.0.D0) then
2190 ZFLUX = U0*FTAU(J1)*RFLECT/(1.d0+RFLECT)
2191 else
2192 ZFLUX = 0.d0
2193 endif
2194 !
2195 !------------------------------------------------------------------------
2196 ! Take optical properties on CTM layers and convert to a photolysis
2197 ! level grid corresponding to layer centres and boundaries. This is
2198 ! required so that J-values can be calculated for the centre of CTM
2199 ! layers; the index of these layers is kept in the jndlev array.
2200 !------------------------------------------------------------------------
2201 !
2202 ! Set lower boundary and levels to calculate J-values at
2203 J1=2*J1-1
2204 do j=1,lpar
2205 jndlev(j)=2*j
2206 enddo
2207 !
2208 ! Calculate column optical depths above each level, TTAU
2209 TTAU(NC+1)=0.0D0
2210 do J=NC,J1,-1
2211 I=(J+1)/2
2212 TTAU(J)=TTAU(J+1) + 0.5d0*DTAUX(I)
2213 jaddlv(j)=int(0.5d0*DTAUX(I)/dtaumax)
2214 ! Subdivide cloud-top levels if required
2215 if(jadsub(j).gt.0) then
2216 jadsub(j)=min(jaddlv(j)+1,nint(dtausub))*(nint(dsubdiv)-1)
2217 jaddlv(j)=jaddlv(j)+jadsub(j)
2218 endif
2219 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in opmie',&
2220 ! j,jadsub(j),jaddlv(j),ttau(j),dtaux(i)
2221 enddo
2222 !
2223 ! Calculate attenuated beam, FTAU, level boundaries then level centres
2224 FTAU(NC+1)=1.0D0
2225 do J=NC-1,J1,-2
2226 I=(J+1)/2
2227 FTAU(J)=FTAU(I)
2228 enddo
2229 do J=NC,J1,-2
2230 FTAU(J)=sqrt(FTAU(J+1)*FTAU(J-1))
2231 enddo
2232 !
2233 ! Calculate scattering properties, level centres then level boundaries
2234 ! using an inverse interpolation to give correctly-weighted values
2235 do j=NC,J1,-2
2236 do i=1,MFIT
2237 pomegaj(i,j) = pomegaj(i,j/2)
2238 enddo
2239 enddo
2240 do j=J1+2,nc,2
2241 taudn = ttau(j-1)-ttau(j)
2242 tauup = ttau(j)-ttau(j+1)
2243 do i=1,MFIT
2244 pomegaj(i,j) = (pomegaj(i,j-1)*taudn + &
2245 pomegaj(i,j+1)*tauup) / (taudn+tauup)
2246 enddo
2247 enddo
2248 ! Define lower and upper boundaries
2249 do i=1,MFIT
2250 pomegaj(i,J1) = pomegaj(i,J1+1)
2251 pomegaj(i,nc+1) = pomegaj(i,nc)
2252 enddo
2253 !
2254 !------------------------------------------------------------------------
2255 ! Calculate cumulative total and define levels we want J-values at.
2256 ! Sum upwards for levels, and then downwards for Mie code readjustments.
2257 !
2258 ! jaddlv(i) Number of new levels to add between (i) and (i+1)
2259 ! jaddto(i) Total number of new levels to add to and above level (i)
2260 ! jndlev(j) Level needed for J-value for CTM layer (j)
2261 !
2262 !------------------------------------------------------------------------
2263 !
2264 ! Reinitialize level arrays
2265 do j=1,nc+1
2266 jaddto(j)=0
2267 enddo
2268 !
2269 jaddto(J1)=jaddlv(J1)
2270 do j=J1+1,nc
2271 jaddto(j)=jaddto(j-1)+jaddlv(j)
2272 enddo
2273 if((jaddto(nc)+nc).gt.nl) then
2274 ! print*,'jdf mie',isvode,jsvode,jaddto(nc),nc,nl
2275 write ( msg, '(a, 2i6)' ) &
2276 'FASTJ Max NL exceeded ' // &
2277 'jaddto(nc)+nc NL', jaddto(nc)+nc,NL
2278 call peg_message( lunerr, msg )
2279 msg = 'FASTJ subr OPMIE error. Max NL exceeded'
2280 call peg_error_fatal( lunerr, msg )
2281 ! write(6,1500) jaddto(nc)+nc, 'NL',NL
2282 ! stop
2283 endif
2284 do i=1,lpar
2285 jndlev(i)=jndlev(i)+jaddto(jndlev(i)-1)
2286 enddo
2287 jaddto(nc)=jaddlv(nc)
2288 do j=nc-1,J1,-1
2289 jaddto(j)=jaddto(j+1)+jaddlv(j)
2290 enddo
2291 !
2292 !---------------------SET UP FOR MIE CODE-------------------------------
2293 !
2294 ! Transpose the ascending TTAU grid to a descending ZTAU grid.
2295 ! Double the resolution - TTAU points become the odd points on the
2296 ! ZTAU grid, even points needed for asymm phase fn soln, contain 'h'.
2297 ! Odd point added at top of grid for unattenuated beam (Z='inf')
2298 !
2299 ! Surface: TTAU(1) now use ZTAU(2*NC+1)
2300 ! Top: TTAU(NC) now use ZTAU(3)
2301 ! Infinity: now use ZTAU(1)
2302 !
2303 ! Mie scattering code only used from surface to level NC
2304 !------------------------------------------------------------------------
2305 !
2306 ! Initialise all Fast-J optical property arrays
2307 do k=1,N__
2308 do i=1,MFIT
2309 pomega(i,k) = 0.d0
2310 enddo
2311 ztau(k) = 0.d0
2312 fz(k) = 0.d0
2313 enddo
2314 !
2315 ! Ascend through atmosphere transposing grid and adding extra points
2316 do j=J1,nc+1
2317 k = 2*(nc+1-j)+2*jaddto(j)+1
2318 ztau(k)= ttau(j)
2319 fz(k) = ftau(j)
2320 do i=1,MFIT
2321 pomega(i,k) = pomegaj(i,j)
2322 enddo
2323 enddo
2324 !
2325 ! Check profiles if desired
2326 ! ND = 2*(NC+jaddto(J1)-J1) + 3
2327 ! if(kw.eq.1) call CH_PROF
2328 !
2329 !------------------------------------------------------------------------
2330 ! Insert new levels, working downwards from the top of the atmosphere
2331 ! to the surface (down in 'j', up in 'k'). This allows ztau and pomega
2332 ! to be incremented linearly (in a +ve sense), and the flux fz to be
2333 ! attenuated top-down (avoiding problems where lower level fluxes are
2334 ! zero).
2335 !
2336 ! zk fractional increment in level
2337 ! dttau change in ttau per increment (linear, positive)
2338 ! dpomega change in pomega per increment (linear)
2339 ! ftaulog change in ftau per increment (exponential, normally < 1)
2340 !
2341 !------------------------------------------------------------------------
2342 !
2343 do j=nc,J1,-1
2344 zk = 0.5d0/(1.d0+dble(jaddlv(j)-jadsub(j)))
2345 dttau = (ttau(j)-ttau(j+1))*zk
2346 do i=1,MFIT
2347 dpomega(i) = (pomegaj(i,j)-pomegaj(i,j+1))*zk
2348 enddo
2349 ! Filter attenuation factor - set minimum at 1.0d-05
2350 if(ftau(j+1).eq.0.d0) then
2351 ftaulog=0.d0
2352 else
2353 ftaulog = ftau(j)/ftau(j+1)
2354 if(ftaulog.lt.1.d-150) then
2355 ftaulog=1.0d-05
2356 else
2357 ftaulog=exp(log(ftaulog)*zk)
2358 endif
2359 endif
2360 k = 2*(nc-j+jaddto(j)-jaddlv(j))+1 ! k at level j+1
2361 l = 0
2362 ! Additional subdivision of first level if required
2363 if(jadsub(j).ne.0) then
2364 l=jadsub(j)/nint(dsubdiv-1)
2365 zk2=1.d0/dsubdiv
2366 dttau2=dttau*zk2
2367 ftaulog2=ftaulog**zk2
2368 do i=1,MFIT
2369 dpomega2(i)=dpomega(i)*zk2
2370 enddo
2371 do ix=1,2*(jadsub(j)+l)
2372 ztau(k+1) = ztau(k) + dttau2
2373 fz(k+1) = fz(k)*ftaulog2
2374 do i=1,MFIT
2375 pomega(i,k+1) = pomega(i,k) + dpomega2(i)
2376 enddo
2377 k = k+1
2378 enddo
2379 endif
2380 l = 2*(jaddlv(j)-jadsub(j)-l)+1
2381 !
2382 ! Add values at all intermediate levels
2383 do ix=1,l
2384 ztau(k+1) = ztau(k) + dttau
2385 fz(k+1) = fz(k)*ftaulog
2386 do i=1,MFIT
2387 pomega(i,k+1) = pomega(i,k) + dpomega(i)
2388 enddo
2389 k = k+1
2390 enddo
2391 !
2392 ! Alternate method to attenuate fluxes, fz, using 2nd-order finite
2393 ! difference scheme - just need to comment in section below
2394 ! ix = 2*(jaddlv(j)-jadsub(j))+1
2395 ! if(l.le.0) then
2396 ! l=k-ix-1
2397 ! else
2398 ! l=k-ix
2399 ! endif
2400 ! call efold(ftau(j+1),ftau(j),ix+1,fz(l))
2401 ! if(jadsub(j).ne.0) then
2402 ! k = 2*(nc-j+jaddto(j)-jaddlv(j))+1 ! k at level j+1
2403 ! ix=2*(jadsub(j)+(jadsub(j)/nint(dsubdiv-1)))
2404 ! call efold(ftau(j+1),fz(k+ix),ix,fz(k))
2405 ! endif
2406 !
2407 enddo
2408 !
2409 !---Update total number of levels and check doesn't exceed N__
2410 ND = 2*(NC+jaddto(J1)-J1) + 3
2411 if(nd.gt.N__) then
2412 write ( msg, '(a, 2i6)' ) &
2413 'FASTJ Max N__ exceeded ' // &
2414 'ND N__', ND, N__
2415 call peg_message( lunerr, msg )
2416 msg = 'FASTJ subr OPMIE error. Max N__ exceeded'
2417 call peg_error_fatal( lunerr, msg )
2418 ! write(6,1500) ND, 'N__',N__
2419 ! stop
2420 endif
2421 !
2422 !---Add boundary/ground layer to ensure no negative J's caused by
2423 !---too large a TTAU-step in the 2nd-order lower b.c.
2424 ZTAU(ND+1) = ZTAU(ND)*1.000005d0
2425 ZTAU(ND+2) = ZTAU(ND)*1.000010d0
2426 zk=max(abs(U0),0.01d0)
2427 zk=dexp(-ZTAU(ND)*5.d-6/zk)
2428 FZ(ND+1) = FZ(ND)*zk
2429 FZ(ND+2) = FZ(ND+1)*zk
2430 do I=1,MFIT
2431 POMEGA(I,ND+1) = POMEGA(I,ND)
2432 POMEGA(I,ND+2) = POMEGA(I,ND)
2433 enddo
2434 ND = ND+2
2435 !
2436 ZU0 = U0
2437 ZREFL = RFLECT
2438 !
2439 !-----------------------------------------
2440 CALL MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2441 ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2442 !-----------------------------------------
2443 ! Accumulate attenuation for selected levels
2444 l=2*(NC+jaddto(J1))+3
2445 do j=1,lpar
2446 k=l-(2*jndlev(j))
2447 if(k.gt.ND-2) then
2448 FMEAN(j) = 0.d0
2449 else
2450 FMEAN(j) = FJ(k)
2451 endif
2452 enddo
2453 !
2454 return
2455 1000 format(1x,i3,3(2x,1pe10.4),1x,i3)
2456 1300 format(1x,50(i3))
2457 1500 format(' Too many levels in photolysis code: need ',i3,' but ',a, &
2458 ' dimensioned as ',i3)
2459 END SUBROUTINE OPMIE
2460
2461 !*********************************************************************
2462 subroutine EFOLD (F0, F1, N, F)
2463 !-----------------------------------------------------------------------
2464 !--- calculate the e-fold between two boundaries, given the value
2465 !--- at both boundaries F0(x=0) = top, F1(x=1) = bottom.
2466 !--- presume that F(x) proportional to exp[-A*x] for x=0 to x=1
2467 !--- d2F/dx2 = A*A*F and thus expect F1 = F0 * exp[-A]
2468 !--- alternatively, could define A = ln[F0/F1]
2469 !--- let X = A*x, d2F/dX2 = F
2470 !--- assume equal spacing (not necessary, but makes this easier)
2471 !--- with N-1 intermediate points (and N layers of thickness dX = A/N)
2472 !---
2473 !--- 2nd-order finite difference: (F(i-1) - 2F(i) + F(i+1)) / dX*dX = F(i)
2474 !--- let D = 1 / dX*dX:
2475 !
2476 ! 1 | 1 0 0 0 0 0 | | F0 |
2477 ! | | | 0 |
2478 ! 2 | -D 2D+1 -D 0 0 0 | | 0 |
2479 ! | | | 0 |
2480 ! 3 | 0 -D 2D+1 -D 0 0 | | 0 |
2481 ! | | | 0 |
2482 ! | 0 0 -D 2D+1 -D 0 | | 0 |
2483 ! | | | 0 |
2484 ! N | 0 0 0 -D 2D+1 -D | | 0 |
2485 ! | | | 0 |
2486 ! N+1 | 0 0 0 0 0 1 | | F1 |
2487 !
2488 !-----------------------------------------------------------------------
2489 ! Advantage of scheme over simple attenuation factor: conserves total
2490 ! number of photons - very useful when using scheme for heating rates.
2491 ! Disadvantage: although reproduces e-folds very well for small flux
2492 ! differences, starts to drift off when many orders of magnitude are
2493 ! involved.
2494 !-----------------------------------------------------------------------
2495 implicit none
2496 real*8 F0,F1,F(250) !F(N+1)
2497 integer N
2498 integer I
2499 real*8 A,DX,D,DSQ,DDP1, B(101),R(101)
2500 !
2501 if(F0.eq.0.d0) then
2502 do I=1,N
2503 F(I)=0.d0
2504 enddo
2505 return
2506 elseif(F1.eq.0.d0) then
2507 A = DLOG(F0/1.d-250)
2508 else
2509 A = DLOG(F0/F1)
2510 endif
2511 !
2512 DX = float(N)/A
2513 D = DX*DX
2514 DSQ = D*D
2515 DDP1 = D+D+1.d0
2516 !
2517 B(2) = DDP1
2518 R(2) = +D*F0
2519 do I=3,N
2520 B(I) = DDP1 - DSQ/B(I-1)
2521 R(I) = +D*R(I-1)/B(I-1)
2522 enddo
2523 F(N+1) = F1
2524 do I=N,2,-1
2525 F(I) = (R(I) + D*F(I+1))/B(I)
2526 enddo
2527 F(1) = F0
2528 return
2529 end subroutine EFOLD
2530
2531
2532 subroutine CH_PROF(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2533 ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2534 !-----------------------------------------------------------------------
2535 ! Check profiles to be passed to MIESCT
2536 !-----------------------------------------------------------------------
2537
2538 USE module_peg_util, only: peg_message
2539
2540 implicit none
2541 !jdf
2542 integer, parameter :: single = 4 !compiler dependent value real*4
2543 integer, parameter :: double = 8 !compiler dependent value real*8
2544 INTEGER NL, N__, M__
2545 PARAMETER (NL=350, N__=2*NL, M__=4)
2546 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2547 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2548 REAL(kind=double), dimension(M__,2*M__) :: PM
2549 REAL(kind=double), dimension(2*M__) :: PM0
2550 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2551 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2552 REAL(kind=double), dimension(M__,M__,N__) :: DD
2553 REAL(kind=double), dimension(M__,N__) :: RR
2554 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2555 INTEGER ND,N,M,MFIT
2556 !jdf
2557 integer i,j
2558 character*80 msg
2559 ! write(6,1100) 'lev','ztau','fz ','pomega( )'
2560 do i=1,ND
2561 if(ztau(i).ne.0.d0) then
2562 write ( msg, '(a, i3, 2(1x,1pe9.3))' ) &
2563 'FASTJ subr CH_PROF ztau ne 0. check pomega. ' // &
2564 'k ztau(k) fz(k) ', i,ztau(i),fz(i)
2565 call peg_message( lunerr, msg )
2566 ! write(6,1200) i,ztau(i),fz(i),(pomega(j,i),j=1,8)
2567 endif
2568 enddo
2569 return
2570 1100 format(1x,a3,4(a9,2x))
2571 1200 format(1x,i3,11(1x,1pe9.3))
2572 end subroutine CH_PROF
2573
2574
2575 SUBROUTINE MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2576 ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2577 ! SUBROUTINE MIESCT
2578 !-----------------------------------------------------------------------
2579 ! This is an adaption of the Prather radiative transfer code, (mjp, 10/95)
2580 ! Prather, 1974, Astrophys. J. 192, 787-792.
2581 ! Sol'n of inhomogeneous Rayleigh scattering atmosphere.
2582 ! (original Rayleigh w/ polarization)
2583 ! Cochran and Trafton, 1978, Ap.J., 219, 756-762.
2584 ! Raman scattering in the atmospheres of the major planets.
2585 ! (first use of anisotropic code)
2586 ! Jacob, Gottlieb and Prather, 1989, J.Geophys.Res., 94, 12975-13002.
2587 ! Chemistry of a polluted cloudy boundary layer,
2588 ! (documentation of extension to anisotropic scattering)
2589 !
2590 ! takes atmospheric structure and source terms from std J-code
2591 ! ALSO limited to 4 Gauss points, only calculates mean field!
2592 !
2593 ! mean rad. field ONLY (M=1)
2594 ! initialize variables FIXED/UNUSED in this special version:
2595 ! FTOP = 1.0 = astrophysical flux (unit of pi) at SZA, -ZU0, use for scaling
2596 ! FBOT = 0.0 = external isotropic flux on lower boundary
2597 ! SISOTP = 0.0 = Specific Intensity of isotropic radiation incident from top
2598 !
2599 ! SUBROUTINES: MIESCT needs 'jv_mie.cmn'
2600 ! BLKSLV needs 'jv_mie.cmn'
2601 ! GEN (ID) needs 'jv_mie.cmn'
2602 ! LEGND0 (X,PL,N)
2603 ! MATIN4 (A)
2604 ! GAUSSP (N,XPT,XWT)
2605 !-----------------------------------------------------------------------
2606
2607 ! INCLUDE 'jv_mie.h'
2608
2609 implicit none
2610 !jdf
2611 integer, parameter :: single = 4 !compiler dependent value real*4
2612 integer, parameter :: double = 8 !compiler dependent value real*8
2613 INTEGER NL, N__, M__
2614 PARAMETER (NL=350, N__=2*NL, M__=4)
2615 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2616 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2617 REAL(kind=double), dimension(M__,2*M__) :: PM
2618 REAL(kind=double), dimension(2*M__) :: PM0
2619 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2620 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2621 REAL(kind=double), dimension(M__,M__,N__) :: DD
2622 REAL(kind=double), dimension(M__,N__) :: RR
2623 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2624 INTEGER ND,N,M,MFIT
2625 !jdf
2626 integer i, id, im
2627 real*8 cmeq1
2628 !-----------------------------------------------------------------------
2629 !---fix scattering to 4 Gauss pts = 8-stream
2630 CALL GAUSSP (N,EMU,WT)
2631 !---solve eqn of R.T. only for first-order M=1
2632 ! ZFLUX = (ZU0*FZ(ND)*ZREFL+FBOT)/(1.0d0+ZREFL)
2633 ZFLUX = (ZU0*FZ(ND)*ZREFL)/(1.0d0+ZREFL)
2634 M=1
2635 DO I=1,N
2636 CALL LEGND0 (EMU(I),PM0,MFIT)
2637 DO IM=M,MFIT
2638 PM(I,IM) = PM0(IM)
2639 ENDDO
2640 ENDDO
2641 !
2642 CMEQ1 = 0.25D0
2643 CALL LEGND0 (-ZU0,PM0,MFIT)
2644 DO IM=M,MFIT
2645 PM0(IM) = CMEQ1*PM0(IM)
2646 ENDDO
2647 !
2648 ! CALL BLKSLV
2649 CALL BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2650 ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2651 !
2652 DO ID=1,ND,2
2653 FJ(ID) = 4.0d0*FJ(ID) + FZ(ID)
2654 ENDDO
2655
2656 RETURN
2657 END SUBROUTINE MIESCT
2658
2659 SUBROUTINE BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2660 ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2661 !-----------------------------------------------------------------------
2662 ! Solves the block tri-diagonal system:
2663 ! A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I)
2664 !-----------------------------------------------------------------------
2665 ! INCLUDE 'jv_mie.h'
2666
2667 implicit none
2668 !jdf
2669 integer, parameter :: single = 4 !compiler dependent value real*4
2670 integer, parameter :: double = 8 !compiler dependent value real*8
2671 INTEGER NL, N__, M__
2672 PARAMETER (NL=350, N__=2*NL, M__=4)
2673 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2674 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2675 REAL(kind=double), dimension(M__,2*M__) :: PM
2676 REAL(kind=double), dimension(2*M__) :: PM0
2677 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2678 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2679 REAL(kind=double), dimension(M__,M__,N__) :: DD
2680 REAL(kind=double), dimension(M__,N__) :: RR
2681 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2682 INTEGER ND,N,M,MFIT
2683 !jdf
2684 integer i, j, k, id
2685 real*8 thesum
2686 !-----------UPPER BOUNDARY ID=1
2687 CALL GEN(1,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2688 ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2689 CALL MATIN4 (B)
2690 DO I=1,N
2691 RR(I,1) = 0.0d0
2692 DO J=1,N
2693 THESUM = 0.0d0
2694 DO K=1,N
2695 THESUM = THESUM - B(I,K)*CC(K,J)
2696 ENDDO
2697 DD(I,J,1) = THESUM
2698 RR(I,1) = RR(I,1) + B(I,J)*H(J)
2699 ENDDO
2700 ENDDO
2701 !----------CONTINUE THROUGH ALL DEPTH POINTS ID=2 TO ID=ND-1
2702 DO ID=2,ND-1
2703 CALL GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2704 ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2705 DO I=1,N
2706 DO J=1,N
2707 B(I,J) = B(I,J) + A(I)*DD(I,J,ID-1)
2708 ENDDO
2709 H(I) = H(I) - A(I)*RR(I,ID-1)
2710 ENDDO
2711 CALL MATIN4 (B)
2712 DO I=1,N
2713 RR(I,ID) = 0.0d0
2714 DO J=1,N
2715 RR(I,ID) = RR(I,ID) + B(I,J)*H(J)
2716 DD(I,J,ID) = - B(I,J)*C1(J)
2717 ENDDO
2718 ENDDO
2719 ENDDO
2720 !---------FINAL DEPTH POINT: ND
2721 CALL GEN(ND,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, &
2722 ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2723 DO I=1,N
2724 DO J=1,N
2725 THESUM = 0.0d0
2726 DO K=1,N
2727 THESUM = THESUM + AA(I,K)*DD(K,J,ND-1)
2728 ENDDO
2729 B(I,J) = B(I,J) + THESUM
2730 H(I) = H(I) - AA(I,J)*RR(J,ND-1)
2731 ENDDO
2732 ENDDO
2733 CALL MATIN4 (B)
2734 DO I=1,N
2735 RR(I,ND) = 0.0d0
2736 DO J=1,N
2737 RR(I,ND) = RR(I,ND) + B(I,J)*H(J)
2738 ENDDO
2739 ENDDO
2740 !-----------BACK SOLUTION
2741 DO ID=ND-1,1,-1
2742 DO I=1,N
2743 DO J=1,N
2744 RR(I,ID) = RR(I,ID) + DD(I,J,ID)*RR(J,ID+1)
2745 ENDDO
2746 ENDDO
2747 ENDDO
2748 !----------MEAN J & H
2749 DO ID=1,ND,2
2750 FJ(ID) = 0.0d0
2751 DO I=1,N
2752 FJ(ID) = FJ(ID) + RR(I,ID)*WT(I)
2753 ENDDO
2754 ENDDO
2755 DO ID=2,ND,2
2756 FJ(ID) = 0.0d0
2757 DO I=1,N
2758 FJ(ID) = FJ(ID) + RR(I,ID)*WT(I)*EMU(I)
2759 ENDDO
2760 ENDDO
2761 ! Output fluxes for testing purposes
2762 ! CALL CH_FLUX
2763 ! CALL CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2764 ! ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2765 !
2766 RETURN
2767 END SUBROUTINE BLKSLV
2768
2769
2770 SUBROUTINE CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2771 ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2772 !-----------------------------------------------------------------------
2773 ! Diagnostic routine to check fluxes at each level - makes most sense
2774 ! when running a conservative atmosphere (zero out absorption in
2775 ! OPMIE by calling the NOABS routine below)
2776 !-----------------------------------------------------------------------
2777
2778 ! INCLUDE 'jv_mie.h'
2779 IMPLICIT NONE
2780 !jdf
2781 integer, parameter :: single = 4 !compiler dependent value real*4
2782 integer, parameter :: double = 8 !compiler dependent value real*8
2783 INTEGER NL, N__, M__
2784 PARAMETER (NL=350, N__=2*NL, M__=4)
2785 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2786 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2787 REAL(kind=double), dimension(M__,2*M__) :: PM
2788 REAL(kind=double), dimension(2*M__) :: PM0
2789 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2790 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2791 REAL(kind=double), dimension(M__,M__,N__) :: DD
2792 REAL(kind=double), dimension(M__,N__) :: RR
2793 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2794 INTEGER ND,N,M,MFIT
2795 !jdf
2796 integer I,ID
2797 real*8 FJCHEK(N__),FZMEAN
2798 !
2799 ! Odd (h) levels held as actinic flux, so recalculate irradiances
2800 DO ID=1,ND,2
2801 FJCHEK(ID) = 0.0d0
2802 DO I=1,N
2803 FJCHEK(ID) = FJCHEK(ID) + RR(I,ID)*WT(I)*EMU(i)
2804 ENDDO
2805 ENDDO
2806 !
2807 ! Even (j) levels are already held as irradiances
2808 DO ID=2,ND,2
2809 DO I=1,N
2810 FJCHEK(ID) = FJ(ID)
2811 ENDDO
2812 ENDDO
2813 !
2814 ! Output Downward and Upward fluxes down through atmosphere
2815 ! WRITE(6,1200)
2816 WRITE(34,1200)
2817 DO ID=2,ND,2
2818 FZMEAN=sqrt(FZ(ID)*FZ(ID-1))
2819 ! WRITE(6,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)), &
2820 WRITE(34,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)), &
2821 2.0*(FJCHEK(id)+FJCHEK(id-1)), &
2822 2.0*(FJCHEK(id)+FJCHEK(id-1))/ &
2823 (ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)))
2824 ENDDO
2825 RETURN
2826 1000 FORMAT(1x,i3,1p,2E12.4,1x,0p,f9.4)
2827 1200 FORMAT(1x,'Lev',3x,'Downward',4x,'Upward',7x,'Ratio')
2828 END SUBROUTINE CH_FLUX
2829
2830 SUBROUTINE NOABS(XLO3,XLO2,XLRAY,BCAER,RFLECT)
2831 !-----------------------------------------------------------------------
2832 ! Zero out absorption terms to check scattering code. Leave a little
2833 ! Rayleigh to provide a minimal optical depth, and set surface albedo
2834 ! to unity.
2835 !-----------------------------------------------------------------------
2836 IMPLICIT NONE
2837 real*8 XLO3,XLO2,XLRAY,BCAER,RFLECT
2838 XLO3=0.d0
2839 XLO2=0.d0
2840 XLRAY=XLRAY*1.d-10
2841 BCAER=0.d0
2842 RFLECT=1.d0
2843 RETURN
2844 END SUBROUTINE NOABS
2845
2846 SUBROUTINE GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, &
2847 ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0)
2848 !-----------------------------------------------------------------------
2849 ! Generates coefficient matrices for the block tri-diagonal system:
2850 ! A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I)
2851 !-----------------------------------------------------------------------
2852
2853 ! INCLUDE 'jv_mie.h'
2854 IMPLICIT NONE
2855 !jdf
2856 integer, parameter :: single = 4 !compiler dependent value real*4
2857 integer, parameter :: double = 8 !compiler dependent value real*8
2858 INTEGER NL, N__, M__
2859 PARAMETER (NL=350, N__=2*NL, M__=4)
2860 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU
2861 REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1
2862 REAL(kind=double), dimension(M__,2*M__) :: PM
2863 REAL(kind=double), dimension(2*M__) :: PM0
2864 REAL(kind=double), dimension(2*M__,N__) :: POMEGA
2865 REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ
2866 REAL(kind=double), dimension(M__,M__,N__) :: DD
2867 REAL(kind=double), dimension(M__,N__) :: RR
2868 REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0
2869 INTEGER ND,N,M,MFIT
2870 !jdf
2871 integer id, id0, id1, im, i, j, k, mstart
2872 real*8 sum0, sum1, sum2, sum3
2873 real*8 deltau, d1, d2, surfac
2874 !---------------------------------------------
2875 IF(ID.EQ.1 .OR. ID.EQ.ND) THEN
2876 !---------calculate generic 2nd-order terms for boundaries
2877 ID0 = ID
2878 ID1 = ID+1
2879 IF(ID.GE.ND) ID1 = ID-1
2880 DO 10 I=1,N
2881 SUM0 = 0.0d0
2882 SUM1 = 0.0d0
2883 SUM2 = 0.0d0
2884 SUM3 = 0.0d0
2885 DO IM=M,MFIT,2
2886 SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM)
2887 SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM)
2888 ENDDO
2889 DO IM=M+1,MFIT,2
2890 SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM)
2891 SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM)
2892 ENDDO
2893 H(I) = 0.5d0*(SUM0*FZ(ID0) + SUM2*FZ(ID1))
2894 A(I) = 0.5d0*(SUM1*FZ(ID0) + SUM3*FZ(ID1))
2895 DO J=1,I
2896 SUM0 = 0.0d0
2897 SUM1 = 0.0d0
2898 SUM2 = 0.0d0
2899 SUM3 = 0.0d0
2900 DO IM=M,MFIT,2
2901 SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM)
2902 SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM)
2903 ENDDO
2904 DO IM=M+1,MFIT,2
2905 SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM)
2906 SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM)
2907 ENDDO
2908 S(I,J) = - SUM2*WT(J)
2909 S(J,I) = - SUM2*WT(I)
2910 W(I,J) = - SUM1*WT(J)
2911 W(J,I) = - SUM1*WT(I)
2912 U1(I,J) = - SUM3*WT(J)
2913 U1(J,I) = - SUM3*WT(I)
2914 SUM0 = 0.5d0*(SUM0 + SUM2)
2915 B(I,J) = - SUM0*WT(J)
2916 B(J,I) = - SUM0*WT(I)
2917 ENDDO
2918 S(I,I) = S(I,I) + 1.0d0
2919 W(I,I) = W(I,I) + 1.0d0
2920 U1(I,I) = U1(I,I) + 1.0d0
2921 B(I,I) = B(I,I) + 1.0d0
2922 10 CONTINUE
2923 DO I=1,N
2924 SUM0 = 0.0d0
2925 DO J=1,N
2926 SUM0 = SUM0 + S(I,J)*A(J)/EMU(J)
2927 ENDDO
2928 C1(I) = SUM0
2929 ENDDO
2930 DO I=1,N
2931 DO J=1,N
2932 SUM0 = 0.0d0
2933 SUM2 = 0.0d0
2934 DO K=1,N
2935 SUM0 = SUM0 + S(J,K)*W(K,I)/EMU(K)
2936 SUM2 = SUM2 + S(J,K)*U1(K,I)/EMU(K)
2937 ENDDO
2938 A(J) = SUM0
2939 V1(J) = SUM2
2940 ENDDO
2941 DO J=1,N
2942 W(J,I) = A(J)
2943 U1(J,I) = V1(J)
2944 ENDDO
2945 ENDDO
2946 IF (ID.EQ.1) THEN
2947 !-------------upper boundary, 2nd-order, C-matrix is full (CC)
2948 DELTAU = ZTAU(2) - ZTAU(1)
2949 D2 = 0.25d0*DELTAU
2950 DO I=1,N
2951 D1 = EMU(I)/DELTAU
2952 DO J=1,N
2953 B(I,J) = B(I,J) + D2*W(I,J)
2954 CC(I,J) = D2*U1(I,J)
2955 ENDDO
2956 B(I,I) = B(I,I) + D1
2957 CC(I,I) = CC(I,I) - D1
2958 ! H(I) = H(I) + 2.0d0*D2*C1(I) + D1*SISOTP
2959 H(I) = H(I) + 2.0d0*D2*C1(I)
2960 A(I) = 0.0d0
2961 ENDDO
2962 ELSE
2963 !-------------lower boundary, 2nd-order, A-matrix is full (AA)
2964 DELTAU = ZTAU(ND) - ZTAU(ND-1)
2965 D2 = 0.25d0*DELTAU
2966 SURFAC = 4.0d0*ZREFL/(1.0d0 + ZREFL)
2967 DO I=1,N
2968 D1 = EMU(I)/DELTAU
2969 H(I) = H(I) - 2.0d0*D2*C1(I)
2970 SUM0 = 0.0d0
2971 DO J=1,N
2972 SUM0 = SUM0 + W(I,J)
2973 ENDDO
2974 SUM0 = D1 + D2*SUM0
2975 SUM1 = SURFAC*SUM0
2976 DO J=1,N
2977 B(I,J) = B(I,J) + D2*W(I,J) - SUM1*EMU(J)*WT(J)
2978 ENDDO
2979 B(I,I) = B(I,I) + D1
2980 H(I) = H(I) + SUM0*ZFLUX
2981 DO J=1,N
2982 AA(I,J) = - D2*U1(I,J)
2983 ENDDO
2984 AA(I,I) = AA(I,I) + D1
2985 C1(I) = 0.0d0
2986 ENDDO
2987 ENDIF
2988 !------------intermediate points: can be even or odd, A & C diagonal
2989 ELSE
2990 DELTAU = ZTAU(ID+1) - ZTAU(ID-1)
2991 MSTART = M + MOD(ID+1,2)
2992 DO I=1,N
2993 A(I) = EMU(I)/DELTAU
2994 C1(I) = -A(I)
2995 SUM0 = 0.0d0
2996 DO IM=MSTART,MFIT,2
2997 SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM0(IM)
2998 ENDDO
2999 H(I) = SUM0*FZ(ID)
3000 DO J=1,I
3001 SUM0 = 0.0d0
3002 DO IM=MSTART,MFIT,2
3003 SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM(J,IM)
3004 ENDDO
3005 B(I,J) = - SUM0*WT(J)
3006 B(J,I) = - SUM0*WT(I)
3007 ENDDO
3008 B(I,I) = B(I,I) + 1.0d0
3009 ENDDO
3010 ENDIF
3011 RETURN
3012 END SUBROUTINE GEN
3013
3014 SUBROUTINE LEGND0 (X,PL,N)
3015 !---Calculates ORDINARY LEGENDRE fns of X (real)
3016 !--- from P[0] = PL(1) = 1, P[1] = X, .... P[N-1] = PL(N)
3017 IMPLICIT NONE
3018 INTEGER N,I
3019 REAL*8 X,PL(N),DEN
3020 !---Always does PL(2) = P[1]
3021 PL(1) = 1.D0
3022 PL(2) = X
3023 DO I=3,N
3024 DEN = (I-1)
3025 PL(I) = PL(I-1)*X*(2.d0-1.D0/DEN) - PL(I-2)*(1.d0-1.D0/DEN)
3026 ENDDO
3027 RETURN
3028 END SUBROUTINE LEGND0
3029
3030 SUBROUTINE MATIN4 (A)
3031 !-----------------------------------------------------------------------
3032 ! invert 4x4 matrix A(4,4) in place with L-U decomposition (mjp, old...)
3033 !-----------------------------------------------------------------------
3034 IMPLICIT NONE
3035 REAL*8 A(4,4)
3036 !---SETUP L AND U
3037 A(2,1) = A(2,1)/A(1,1)
3038 A(2,2) = A(2,2)-A(2,1)*A(1,2)
3039 A(2,3) = A(2,3)-A(2,1)*A(1,3)
3040 A(2,4) = A(2,4)-A(2,1)*A(1,4)
3041 A(3,1) = A(3,1)/A(1,1)
3042 A(3,2) = (A(3,2)-A(3,1)*A(1,2))/A(2,2)
3043 A(3,3) = A(3,3)-A(3,1)*A(1,3)-A(3,2)*A(2,3)
3044 A(3,4) = A(3,4)-A(3,1)*A(1,4)-A(3,2)*A(2,4)
3045 A(4,1) = A(4,1)/A(1,1)
3046 A(4,2) = (A(4,2)-A(4,1)*A(1,2))/A(2,2)
3047 A(4,3) = (A(4,3)-A(4,1)*A(1,3)-A(4,2)*A(2,3))/A(3,3)
3048 A(4,4) = A(4,4)-A(4,1)*A(1,4)-A(4,2)*A(2,4)-A(4,3)*A(3,4)
3049 !---INVERT L
3050 A(4,3) = -A(4,3)
3051 A(4,2) = -A(4,2)-A(4,3)*A(3,2)
3052 A(4,1) = -A(4,1)-A(4,2)*A(2,1)-A(4,3)*A(3,1)
3053 A(3,2) = -A(3,2)
3054 A(3,1) = -A(3,1)-A(3,2)*A(2,1)
3055 A(2,1) = -A(2,1)
3056 !---INVERT U
3057 A(4,4) = 1.D0/A(4,4)
3058 A(3,4) = -A(3,4)*A(4,4)/A(3,3)
3059 A(3,3) = 1.D0/A(3,3)
3060 A(2,4) = -(A(2,3)*A(3,4)+A(2,4)*A(4,4))/A(2,2)
3061 A(2,3) = -A(2,3)*A(3,3)/A(2,2)
3062 A(2,2) = 1.D0/A(2,2)
3063 A(1,4) = -(A(1,2)*A(2,4)+A(1,3)*A(3,4)+A(1,4)*A(4,4))/A(1,1)
3064 A(1,3) = -(A(1,2)*A(2,3)+A(1,3)*A(3,3))/A(1,1)
3065 A(1,2) = -A(1,2)*A(2,2)/A(1,1)
3066 A(1,1) = 1.D0/A(1,1)
3067 !---MULTIPLY (U-INVERSE)*(L-INVERSE)
3068 A(1,1) = A(1,1)+A(1,2)*A(2,1)+A(1,3)*A(3,1)+A(1,4)*A(4,1)
3069 A(1,2) = A(1,2)+A(1,3)*A(3,2)+A(1,4)*A(4,2)
3070 A(1,3) = A(1,3)+A(1,4)*A(4,3)
3071 A(2,1) = A(2,2)*A(2,1)+A(2,3)*A(3,1)+A(2,4)*A(4,1)
3072 A(2,2) = A(2,2)+A(2,3)*A(3,2)+A(2,4)*A(4,2)
3073 A(2,3) = A(2,3)+A(2,4)*A(4,3)
3074 A(3,1) = A(3,3)*A(3,1)+A(3,4)*A(4,1)
3075 A(3,2) = A(3,3)*A(3,2)+A(3,4)*A(4,2)
3076 A(3,3) = A(3,3)+A(3,4)*A(4,3)
3077 A(4,1) = A(4,4)*A(4,1)
3078 A(4,2) = A(4,4)*A(4,2)
3079 A(4,3) = A(4,4)*A(4,3)
3080 RETURN
3081 END SUBROUTINE MATIN4
3082
3083 SUBROUTINE GAUSSP (N,XPT,XWT)
3084 !-----------------------------------------------------------------------
3085 ! Loads in pre-set Gauss points for 4 angles from 0 to +1 in cos(theta)=mu
3086 !-----------------------------------------------------------------------
3087 IMPLICIT NONE
3088 INTEGER N,I
3089 REAL*8 XPT(N),XWT(N)
3090 REAL*8 GPT4(4),GWT4(4)
3091 DATA GPT4/.06943184420297D0,.33000947820757D0,.66999052179243D0, &
3092 .93056815579703D0/
3093 DATA GWT4/.17392742256873D0,.32607257743127D0,.32607257743127D0, &
3094 .17392742256873D0/
3095 N = 4
3096 DO I=1,N
3097 XPT(I) = GPT4(I)
3098 XWT(I) = GWT4(I)
3099 ENDDO
3100 RETURN
3101 END SUBROUTINE GAUSSP
3102 !
3103 subroutine aeroden(zz,v,aerd)
3104
3105 ! purpose: find number density of boundary layer aerosols, aerd,
3106 ! at a given altitude, zz, and for a specified visibility
3107 ! input:
3108 ! zz altitude (km)
3109 ! v visibility for a horizontal surface path (km)
3110 ! output:
3111 ! aerd aerosol density at altitude z
3112
3113 ! the vertical distribution of the boundary layer aerosol density is
3114 ! based on the 5s vertical profile models for 5 and 23 km visibility.
3115 ! above 5 km, the aden05 and aden23 models are the same
3116 ! below 5 km, the models differ as follows;
3117 ! aden05 0.99 km scale height (94% of extinction occurs below 5 km)
3118 ! aden23 1.45 km scale heigth (80% of extinction occurs below 5 km)
3119 !
3120
3121 implicit none
3122 integer mz, nz
3123 parameter (mz=33)
3124 real v,aerd
3125 real*8 zz ! compatability with fastj
3126 real alt, aden05, aden23, aer05,aer23
3127 dimension alt(mz),aden05(mz),aden23(mz)
3128 !jdf dimension zbaer(*),dbaer(*)
3129
3130 real z, f, wth
3131 integer k, kp
3132 save alt,aden05,aden23,nz
3133
3134
3135 data nz/mz/
3136
3137 data alt/ &
3138 0.0, 1.0, 2.0, 3.0, 4.0, &
3139 5.0, 6.0, 7.0, 8.0, 9.0, &
3140 10.0, 11.0, 12.0, 13.0, 14.0, &
3141 15.0, 16.0, 17.0, 18.0, 19.0, &
3142 20.0, 21.0, 22.0, 23.0, 24.0, &
3143 25.0, 30.0, 35.0, 40.0, 45.0, &
3144 50.0, 70.0, 100.0/
3145 data aden05/ &
3146 1.378E+04, 5.030E+03, 1.844E+03, 6.731E+02, 2.453E+02, &
3147 8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, &
3148 5.675E+01, 5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, &
3149 4.744E+01, 4.511E+01, 4.458E+01, 4.314E+01, 3.634E+01, &
3150 2.667E+01, 1.933E+01, 1.455E+01, 1.113E+01, 8.826E+00, &
3151 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01, 4.082E-02, &
3152 1.078E-02, 5.550E-05, 1.969E-08/
3153 data aden23/ &
3154 2.828E+03, 1.244E+03, 5.371E+02, 2.256E+02, 1.192E+02, &
3155 8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, &
3156 5.675E+01, 5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, &
3157 4.744E+01, 4.511E+01, 4.458E+01, 4.314E+01, 3.634E+01, &
3158 2.667E+01, 1.933E+01, 1.455E+01, 1.113E+01, 8.826E+00, &
3159 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01, 4.082E-02, &
3160 1.078E-02, 5.550E-05, 1.969E-08/
3161 !
3162 z=max(0.,min(100.,real(zz)))
3163 aerd=0.
3164 if(z.gt.alt(nz)) return
3165
3166 call locate(alt,nz,z,k)
3167 kp=k+1
3168 f=(z-alt(k))/(alt(kp)-alt(k))
3169
3170 if(min(aden05(k),aden05(kp),aden23(k),aden23(kp)).le.0.) then
3171 aer05=aden05(k)*(1.-f)+aden05(kp)*f
3172 aer23=aden23(k)*(1.-f)+aden23(kp)*f
3173 else
3174 aer05=aden05(k)*(aden05(kp)/aden05(k))**f
3175 aer23=aden23(k)*(aden23(kp)/aden23(k))**f
3176 endif
3177
3178 wth=(1./v-1/5.)/(1./23.-1./5.)
3179 wth=max(0.,min(1.,wth))
3180
3181 aerd=(1.-wth)*aer05+wth*aer23
3182
3183 ! write(*,*) 'aeroden k,kp,z,aer05(k),aer05(kp),f,aerd'
3184 ! write(*,'(2i5,1p5e11.3)') k,kp,z,aden05(k),aden05(kp),f,aerd
3185
3186 return
3187 end subroutine aeroden
3188 !=======================================================================
3189 subroutine locate(xx,n,x,j)
3190 !
3191 ! purpose: given an array xx of length n, and given a value X, returns
3192 ! a value J such that X is between xx(j) and xx(j+1). xx must
3193 ! be monotonic, either increasing of decreasing. this function
3194 ! returns j=1 or j=n-1 if x is out of range.
3195 !c
3196 ! input:
3197 ! xx monitonic table
3198 ! n size of xx
3199 ! x single floating point value perhaps within the range of xx
3200 !
3201 ! output:
3202 ! function returns index value j, such that
3203 !
3204 ! for an increasing table
3205 !
3206 ! xx(j) .lt. x .le. xx(j+1),
3207 ! j=1 for x .lt. xx(1)
3208 ! j=n-1 for x .gt. xx(n)
3209 !
3210 ! for a decreasing table
3211 ! xx(j) .le. x .lt. xx(j+1)
3212 ! j=n-1 for x .lt. xx(n)
3213 ! j=1 for x .gt. xx(1)
3214 !
3215 implicit none
3216 integer j,n
3217 real x,xx(n)
3218 integer jl,jm,ju
3219
3220 ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3221
3222 if(x.eq.xx(1)) then
3223 j=1
3224 return
3225 endif
3226 if(x.eq.xx(n)) then
3227 j=n-1
3228 return
3229 endif
3230 jl=1
3231 ju=n
3232 10 if(ju-jl.gt.1) then
3233 jm=(ju+jl)/2
3234 if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then
3235 jl=jm
3236 else
3237 ju=jm
3238 endif
3239 goto 10
3240 endif
3241 j=jl
3242 return
3243 end subroutine locate
3244 !************************************************************************
3245 subroutine rd_tjpl2
3246 !-----------------------------------------------------------------------
3247 ! set wavelength bins, solar fluxes, Rayleigh parameters, temperature-
3248 ! dependent cross sections and Rayleigh/aerosol scattering phase functions
3249 ! with temperature dependences. Current data originates from JPL 2000
3250 !-----------------------------------------------------------------------
3251 !
3252 ! NJVAL Number of species to calculate J-values for
3253 ! NWWW Number of wavelength bins, from NW1:NW2
3254 ! WBIN Boundaries of wavelength bins
3255 ! WL Centres of wavelength bins - 'effective wavelength'
3256 ! FL Solar flux incident on top of atmosphere (cm-2.s-1)
3257 ! QRAYL Rayleigh parameters (effective cross-section) (cm2)
3258 ! QBC Black Carbon absorption extinct. (specific cross-sect.) (m2/g)
3259 ! QO2 O2 cross-sections
3260 ! QO3 O3 cross-sections
3261 ! Q1D O3 => O(1D) quantum yield
3262 ! TQQ Temperature for supplied cross sections
3263 ! QQQ Supplied cross sections in each wavelength bin (cm2)
3264 ! NAA Number of categories for scattering phase functions
3265 ! QAA Aerosol scattering phase functions
3266 ! NK Number of wavelengths at which functions supplied (set as 4)
3267 ! WAA Wavelengths for the NK supplied phase functions
3268 ! PAA Phase function: first 8 terms of expansion
3269 ! RAA Effective radius associated with aerosol type
3270 ! SSA Single scattering albedo
3271 !
3272 ! npdep Number of pressure dependencies
3273 ! zpdep Pressure dependencies by wavelength bin
3274 ! jpdep Index of cross sections requiring pressure dependence
3275 ! lpdep Label for pressure dependence
3276 !
3277 !-----------------------------------------------------------------------
3278
3279 USE module_data_mosaic_other, only : kmaxd
3280 USE module_fastj_data
3281 USE module_peg_util, only: peg_message, peg_error_fatal
3282
3283 IMPLICIT NONE
3284 !jdf
3285 ! Print Fast-J diagnostics if iprint /= 0
3286 integer, parameter :: iprint = 0
3287 integer, parameter :: single = 4 !compiler dependent value real*4
3288 ! integer, parameter :: double = 8 !compiler dependent value real*8
3289 integer,parameter :: ipar_fastj=1,jpar=1
3290 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied
3291 logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution
3292 integer lpar !Number of levels in CTM
3293 integer jpnl !Number of levels requiring chemistry
3294 real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians)
3295 real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians)
3296 real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees)
3297 real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries
3298 real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries
3299 real(kind=double) :: tau_fastj ! Time of Day (hours, GMT)
3300 integer month_fastj ! Number of month (1-12)
3301 integer iday_fastj ! Day of year
3302 !jdf
3303 integer i, j, k
3304 character*7 lpdep(3)
3305 character*80 msg
3306
3307 if(NJVAL.gt.NS) then
3308 ! fastj input files are not set up for current situation
3309 write ( msg, '(a, 2i6)' ) &
3310 'FASTJ # xsect supplied > max allowed ' // &
3311 'NJVAL NS ', NJVAL, NS
3312 call peg_message( lunerr, msg )
3313 msg = &
3314 'FASTJ Setup Error: # xsect supplied > max allowed. Increase NS'
3315 call peg_error_fatal( lunerr, msg )
3316 ! write(6,300) NJVAL,NS
3317 ! stop
3318 endif
3319
3320 if(NAA.gt.NP) then
3321 write ( msg, '(a, 2i6)' ) &
3322 'FASTJ # aerosol/cloud types > NP ' // &
3323 'NAA NP ', NAA ,NP
3324 call peg_message( lunerr, msg )
3325 msg = &
3326 'FASTJ Setup Error: Too many phase functions supplied. Increase NP'
3327 call peg_error_fatal( lunerr, msg )
3328 ! write(6,350) NAA
3329 ! stop
3330 endif
3331 !---Zero index arrays
3332 do j=1,jppj
3333 jind(j)=0
3334 enddo
3335 do j=1,NJVAL
3336 jpdep(j)=0
3337 enddo
3338 do j=1,nh
3339 hzind(j)=0
3340 enddo
3341 !
3342 !---Set mapping index
3343 do j=1,NJVAL
3344 do k=1,jppj
3345 if(jlabel(k).eq.titlej(1,j)) jind(k)=j
3346 ! write(6,*)k,jind(k) ! jcb
3347 ! write(6,*)jlabel(k),titlej(1,j) ! jcb
3348 enddo
3349 do k=1,npdep
3350 if(lpdep(k).eq.titlej(1,j)) jpdep(j)=k
3351 enddo
3352 enddo
3353 do k=1,jppj
3354 if(jfacta(k).eq.0.d0) then
3355 ! write(6,*) 'Not using photolysis reaction ',k
3356 write ( msg, '(a, i6)' ) &
3357 'FASTJ Not using photolysis reaction ' , k
3358 call peg_message( lunerr, msg )
3359 end if
3360 if(jind(k).eq.0) then
3361 if(jfacta(k).eq.0.d0) then
3362 jind(k)=1
3363 else
3364 write ( msg, '(a, i6)' ) &
3365 'FASTJ Which J-rate for photolysis reaction ' , k
3366 call peg_message( lunerr, msg )
3367 ! write(6,*) 'Which J-rate for photolysis reaction ',k,' ?'
3368 ! stop
3369 msg = 'FASTJ subr rd_tjpl2 Unknown Jrate. Incorrect FASTJ setup'
3370 call peg_error_fatal( lunerr, msg )
3371 endif
3372 endif
3373 enddo
3374 ! Herzberg index
3375 i=0
3376 do j=1,nhz
3377 do k=1,jppj
3378 if(jlabel(k).eq.hzlab(j)) then
3379 i=i+1
3380 hzind(i)=k
3381 hztoa(i)=hztmp(j)*jfacta(k)
3382 endif
3383 enddo
3384 enddo
3385 nhz=i
3386 if(nhz.eq.0) then
3387 if(iprint.ne.0) then
3388 write ( msg, '(a)' ) &
3389 'FASTJ Not using Herzberg bin '
3390 call peg_message( lunerr, msg )
3391 ! write(6,400)
3392 end if
3393 else
3394 if(iprint.ne.0) then
3395 write ( msg, '(a)' ) &
3396 'FASTJ Using Herzberg bin for: '
3397 call peg_message( lunerr, msg )
3398 write( msg, '(a,10a7)' ) &
3399 'FASTJ ', (jlabel(hzind(i)),i=1,nhz)
3400 ! write(6,420) (jlabel(hzind(i)),i=1,nhz)
3401 end if
3402 endif
3403
3404 ! 300 format(' Number of x-sections supplied to Fast-J: ',i3,/, &
3405 ! ' Maximum number allowed (NS) only set to: ',i3, &
3406 ! ' - increase in jv_cmn.h',/, &
3407 ! 'RESULTS WILL BE IN ERROR' )
3408 ! 350 format(' Too many phase functions supplied; increase NP to ',i2, &
3409 ! /,'RESULTS WILL BE IN ERROR' )
3410 ! 400 format(' Not using Herzberg bin')
3411 ! 420 format(' Using Herzberg bin for: ',10a7)
3412
3413
3414 return
3415 end subroutine rd_tjpl2
3416 !********************************************************************
3417
3418
3419 end module module_phot_fastj