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