source: trunk/LMDZ.MARS/libf/aeronomars/photolysis_online.F

Last change on this file was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

  • remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
  • turn aeronomars routines to modules, for those which aren't in modules yet.

EM
File size: 52.5 KB
Line 
1      module photolysis_online_mod
2     
3      implicit none
4     
5      contains
6
7!==============================================================================
8
9      subroutine photolysis_online(nlayer, deutchem, nb_phot_max,
10     $           alt, press, temp, mmean,
11     $           i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
12     $           i_h2, i_oh, i_ho2, i_h2o2, i_h2o,               
13     $           i_n, i_n2d, i_no, i_no2, i_n2, nesp, rm,
14     $           tau, sza, dist_sol, v_phot)
15
16      use photolysis_mod
17
18      implicit none
19
20!     input
21
22      logical, intent(in) :: deutchem
23      integer, intent(in) :: nesp                                    ! total number of chemical species
24      integer, intent(in) :: nlayer
25      integer, intent(in) :: nb_phot_max
26      integer, intent(in) :: i_co2, i_co, i_o, i_o1d, i_o2, i_o3, i_h,
27     $                       i_h2, i_oh, i_ho2, i_h2o2, i_h2o,
28     $                       i_n, i_n2d, i_no, i_no2, i_n2
29
30      real, dimension(nlayer), intent(in) :: press, temp, mmean      ! pressure (hpa)/temperature (k)/mean molecular mass (g.mol-1)
31      real, dimension(nlayer), intent(in) :: alt                     ! altitude (km)
32      real, dimension(nlayer,nesp), intent(in) :: rm                 ! mixing ratios
33      real, intent(in) :: tau                                        ! integrated aerosol optical depth at the surface
34      real, intent(in) :: sza                                        ! solar zenith angle (degrees)
35      real, intent(in) :: dist_sol                                   ! solar distance (au)
36
37!     output
38
39      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot       ! photolysis rates (s-1)
40 
41!     solar flux at mars
42
43      real, dimension(nw) :: fmars                                   ! solar flux (w.m-2.nm-1)
44      real :: factor
45
46!     cross-sections
47
48      real, dimension(nlayer,nw,nphot) :: sj                         ! general cross-section array (cm2)
49
50!     atmosphere
51
52      real, dimension(nlayer+1) :: zpress, zalt, ztemp, zmmean       ! pressure (hpa)/altitude (km)/temperature (k)/mean molecular mass (g.mol-1)
53
54      real, dimension(nlayer+1) :: colinc                            ! air column increment (molecule.cm-2)
55      real, dimension(nlayer+1,nw) :: dtrl                           ! rayleigh optical depth
56      real, dimension(nlayer+1,nw) :: dtaer                          ! aerosol optical depth
57      real, dimension(nlayer+1,nw) :: omaer                          ! aerosol single scattering albedo
58      real, dimension(nlayer+1,nw) :: gaer                           ! aerosol asymmetry parameter
59      real, dimension(nlayer+1,nw) :: dtcld                          ! cloud optical depth
60      real, dimension(nlayer+1,nw) :: omcld                          ! cloud single scattering albedo
61      real, dimension(nlayer+1,nw) :: gcld                           ! cloud asymmetry parameter
62      real, dimension(nlayer+1,nw,nabs) :: dtgas                     ! optical depth for each gas
63      real, dimension(nlayer+1,nw) :: dagas                          ! total gas optical depth
64      real, dimension(nlayer+1) :: edir, edn, eup                    ! normalised irradiances
65      real, dimension(nlayer+1) :: fdir, fdn, fup                    ! normalised actinic fluxes
66      real, dimension(nlayer+1) :: saflux                            ! total actinic flux
67
68      integer, dimension(0:nlayer+1) :: nid
69      real, dimension(0:nlayer+1,nlayer+1) :: dsdh
70
71      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d, j_o3_o,
72     $           j_h2o, j_h2o2, j_ho2, j_h2, j_no, j_no2, j_n2,
73     $           j_hdo_od, j_hdo_d
74
75      integer :: a_o2, a_co2, a_o3, a_h2o, a_h2o2, a_ho2, a_h2, a_no,
76     $           a_no2, a_n2
77      integer :: nlev, i, ilay, ilev, iw, ialt
78      real    :: deltaj
79
80!     absorbing gas numbering
81
82      a_o2      = 1      ! o2
83      a_co2     = 2      ! co2
84      a_o3      = 3      ! o3
85      a_h2o     = 4      ! h2o
86      a_h2o2    = 5      ! h2o2
87      a_ho2     = 6      ! ho2
88      a_h2      = 7      ! h2
89      a_no      = 8      ! no
90      a_no2     = 9      ! no2
91      a_n2      = 10     ! n2
92
93!     photodissociation rates numbering.
94!     photodissociations must be ordered the same way in subroutine "indice"
95
96      j_o2_o    = 1      ! o2 + hv     -> o + o
97      j_o2_o1d  = 2      ! o2 + hv     -> o + o(1d)
98      j_co2_o   = 3      ! co2 + hv    -> co + o
99      j_co2_o1d = 4      ! co2 + hv    -> co + o(1d)
100      j_o3_o1d  = 5      ! o3 + hv     -> o2 + o(1d)
101      j_o3_o    = 6      ! o3 + hv     -> o2 + o
102      j_h2o     = 7      ! h2o + hv    -> h + oh
103      j_h2o2    = 8      ! h2o2 + hv   -> oh + oh
104      j_ho2     = 9      ! ho2 + hv    -> oh + o
105      j_h2      = 10     ! h2 + hv     -> h + h
106      j_no      = 11     ! no + hv     -> n + o
107      j_no2     = 12     ! no2 + hv    -> no + o
108      j_n2      = 13     ! n2 + hv     -> n + n
109      j_hdo_od  = 14     ! hdo + hv    -> od + h
110      j_hdo_d   = 15     ! hdo + hv    -> d + oh
111
112!==== define working vertical grid for the uv radiative code
113
114      nlev = nlayer + 1
115
116      do ilev = 1,nlev-1
117         zpress(ilev) = press(ilev)
118         zalt(ilev)   = alt(ilev)
119         ztemp(ilev)  = temp(ilev)
120         zmmean(ilev) = mmean(ilev)
121      end do
122
123      zpress(nlev) = 0.                  ! top of the atmosphere
124      zalt(nlev)   = zalt(nlev-1) + (zalt(nlev-1) - zalt(nlev-2))
125      ztemp(nlev)  = ztemp(nlev-1)
126      zmmean(nlev) = zmmean(nlev-1)
127     
128!==== air column increments and rayleigh optical depth
129
130      call setair(nlev, nw, wl, wc, zpress, ztemp, zmmean, colinc, dtrl)
131
132!==== set temperature-dependent cross-sections and optical depths
133
134      dtgas(:,:,:) = 0.
135
136!     o2:
137
138      call seto2(nphot, nlayer, nw, wc, mopt, temp, xso2_150, xso2_200,
139     $           xso2_250, xso2_300, yieldo2, j_o2_o, j_o2_o1d,
140     $           colinc(1:nlayer), rm(:,i_o2),
141     $           dtgas(1:nlayer,:,a_o2), sj)
142
143!     co2:
144
145      call setco2(nphot, nlayer, nw, wc, temp, xsco2_195, xsco2_295,
146     $            xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
147     $            colinc(1:nlayer), rm(:,i_co2),
148     $            dtgas(1:nlayer,:,a_co2), sj)
149
150!     o3:
151
152      call seto3(nphot, nlayer, nw, wc, temp, xso3_218, xso3_298,
153     $           j_o3_o, j_o3_o1d,
154     $           colinc(1:nlayer), rm(:,i_o3),
155     $           dtgas(1:nlayer,:,a_o3), sj)
156
157!     h2o2:
158
159      call seth2o2(nphot, nlayer, nw, wc, temp, xsh2o2, j_h2o2,
160     $             colinc(1:nlayer), rm(:,i_h2o2),
161     $             dtgas(1:nlayer,:,a_h2o2), sj)
162
163!     no2:
164
165      call setno2(nphot, nlayer, nw, wc, temp, xsno2, xsno2_220,
166     $            xsno2_294, yldno2_248, yldno2_298, j_no2,
167     $            colinc(1:nlayer), rm(:,i_no2),
168     $            dtgas(1:nlayer,:,a_no2), sj)
169
170!==== temperature-independent optical depths and cross-sections
171
172!     optical depths
173
174      do ilay = 1,nlayer
175         do iw = 1,nw-1
176            dtgas(ilay,iw,a_h2o) = colinc(ilay)*rm(ilay,i_h2o)*xsh2o(iw)
177            dtgas(ilay,iw,a_ho2) = colinc(ilay)*rm(ilay,i_ho2)*xsho2(iw)
178            dtgas(ilay,iw,a_h2)  = colinc(ilay)*rm(ilay,i_h2)*xsh2(iw)
179            dtgas(ilay,iw,a_no)  = colinc(ilay)*rm(ilay,i_no)*xsno(iw)
180            dtgas(ilay,iw,a_n2)  = colinc(ilay)*rm(ilay,i_n2)*xsn2(iw)
181         end do
182      end do
183
184!     total gas optical depth
185
186      dagas(:,:) = 0.
187
188      do ilay = 1,nlayer
189         do iw = 1,nw-1
190            do i = 1,nabs
191               dagas(ilay,iw) = dagas(ilay,iw) + dtgas(ilay,iw,i)
192            end do
193         end do
194      end do
195
196!     cross-sections
197
198      do ilay = 1,nlayer
199         do iw = 1,nw-1
200            sj(ilay,iw,j_h2o) = xsh2o(iw)                ! h2o
201            sj(ilay,iw,j_ho2) = xsho2(iw)                ! ho2
202            sj(ilay,iw,j_h2)  = xsh2(iw)*yieldh2(iw)     ! h2
203            sj(ilay,iw,j_no)  = xsno(iw)*yieldno(iw)     ! no
204            sj(ilay,iw,j_n2)  = xsn2(iw)*yieldn2(iw)     ! n2
205         end do
206      end do
207
208!     if deuterium chemistry: hdo cross-section
209
210      if (deutchem) then
211         do ilay = 1,nlayer
212            do iw = 1,nw-1
213               ! two chanels for hdo, with same cross-section
214               sj(ilay,iw,j_hdo_od) = 0.5*xshdo(iw) ! hdo -> od + h
215               sj(ilay,iw,j_hdo_d)  = 0.5*xshdo(iw) ! hdo -> d + oh
216            end do
217         end do
218      end if
219
220!==== set aerosol properties and optical depth
221
222      call setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
223
224!==== set cloud properties and optical depth
225
226      call setcld(nlev,nw,dtcld,omcld,gcld)
227
228!==== slant path lengths in spherical geometry
229
230      call sphers(nlev,zalt,sza,dsdh,nid)
231
232!==== solar flux at mars
233
234      factor = (1./dist_sol)**2.
235      do iw = 1,nw-1
236         fmars(iw) = f(iw)*factor
237      end do
238
239!==== initialise photolysis rates
240
241      v_phot(:,1:nphot) = 0.
242
243!==== start of wavelength lopp
244
245      do iw = 1,nw-1
246
247!     monochromatic radiative transfer. outputs are:
248!     normalized irradiances     edir(nlev), edn(nlev), eup(nlev)
249!     normalized actinic fluxes  fdir(nlev), fdn(nlev), fup(nlev)
250!     where
251!     dir = direct beam, dn = down-welling diffuse, up = up-welling diffuse
252
253         call rtlink(nlev, nw, iw, albedo(iw), sza, dsdh, nid, dtrl,
254     $               dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
255     $               edir, edn, eup, fdir, fdn, fup)
256
257
258!     spherical actinic flux
259
260         do ilay = 1,nlayer
261           saflux(ilay) = fmars(iw)*(fdir(ilay) + fdn(ilay) + fup(ilay))
262         end do
263
264!     photolysis rate integration
265
266         do i = 1,nphot
267            do ilay = 1,nlayer
268               deltaj = saflux(ilay)*sj(ilay,iw,i)
269               v_phot(ilay,i) = v_phot(ilay,i) + deltaj*(wu(iw)-wl(iw))
270            end do
271         end do
272
273!     eliminate small values
274
275         where (v_phot(:,1:nphot) < 1.e-30)
276            v_phot(:,1:nphot) = 0.
277         end where
278
279      end do ! iw
280
281      contains
282
283!==============================================================================
284
285      subroutine setair(nlev, nw, wl, wc, press, temp, zmmean,
286     $                  colinc, dtrl)
287
288*-----------------------------------------------------------------------------*
289*=  PURPOSE:                                                                 =*
290*=  computes air column increments and rayleigh optical depth                =*
291*-----------------------------------------------------------------------------*
292
293      implicit none
294
295!     input:
296
297      integer :: nlev, nw
298
299      real, dimension(nw)   :: wl, wc              ! lower and central wavelength grid (nm)
300      real, dimension(nlev) :: press, temp, zmmean ! pressure (hpa), temperature (k), molecular mass (g.mol-1)
301
302!     output:
303
304      real, dimension(nlev) :: colinc    ! air column increments (molecule.cm-2)
305      real, dimension(nlev,nw) :: dtrl   ! rayleigh optical depth
306
307!     local:
308
309      real, parameter :: avo = 6.022e23
310      real, parameter :: g = 3.72
311      real :: dp, nu
312      real, dimension(nw) :: srayl
313      integer :: ilev, iw
314
315!     compute column increments
316
317      do ilev = 1, nlev-1
318         dp = (press(ilev) - press(ilev+1))*100.
319         colinc(ilev) = avo*0.1*dp/(zmmean(ilev)*g)
320      end do
321      colinc(nlev) = 0.
322
323      do iw = 1, nw - 1
324
325!        co2 rayleigh cross-section
326!        ityaksov et al., chem. phys. lett., 462, 31-34, 2008
327
328         nu = 1./(wc(iw)*1.e-7)
329         srayl(iw) = 1.78e-26*nu**(4. + 0.625)
330         srayl(iw) = srayl(iw)*1.e-20 ! cm2
331
332         do ilev = 1, nlev
333            dtrl(ilev,iw) = colinc(ilev)*srayl(iw)   ! cm2
334         end do
335      end do
336
337      end subroutine setair
338
339!==============================================================================
340
341      subroutine setco2(nd, nlayer, nw, wc, tlay, xsco2_195, xsco2_295,
342     $                  xsco2_370, yieldco2, j_co2_o, j_co2_o1d,
343     $                  colinc, rm, dt, sj)
344
345!-----------------------------------------------------------------------------*
346!=  PURPOSE:                                                                 =*
347!=  Set up the CO2 temperature-dependent cross-sections and optical depth    =*
348!-----------------------------------------------------------------------------*
349
350      implicit none
351
352!     input:
353
354      integer :: nd                                              ! number of photolysis rates
355      integer :: nlayer                                          ! number of vertical layers
356      integer :: nw                                              ! number of wavelength grid points
357      integer :: j_co2_o, j_co2_o1d                              ! photolysis indexes
358      real, dimension(nw)     :: wc                              ! central wavelength for each interval
359      real, dimension(nw)     :: xsco2_195, xsco2_295, xsco2_370 ! co2 cross-sections (cm2)
360      real, dimension(nw)     :: yieldco2                        ! co2 photodissociation yield
361      real, dimension(nlayer) :: tlay                            ! temperature (k)
362      real, dimension(nlayer) :: rm                              ! co2 mixing ratio
363      real, dimension(nlayer) :: colinc                          ! air column increment (molecule.cm-2)
364
365!     output:
366
367      real, dimension(nlayer,nw)    :: dt                        !  optical depth
368      real, dimension(nlayer,nw,nd) :: sj                        !  cross-section array (cm2)
369
370!     local:
371
372      integer :: extrapol
373      integer :: i, l
374      real    :: temp, sco2
375
376!     extrapol  = 0         no extrapolation below 195 k
377!     extrapol  = 1         extrapolation below 195 k
378
379      extrapol  = 0
380
381      do i = 1, nlayer
382         if (extrapol == 1) then
383            temp = tlay(i)             
384         else
385            temp = max(tlay(i), 195.)
386         end if
387         temp = min(temp, 370.)
388         do l = 1, nw-1
389            if (temp <= 295.) then
390               if (xsco2_195(l) /= 0. .and. xsco2_295(l) /= 0.) then
391                  sco2 =  alog(xsco2_195(l))
392     $                 + (alog(xsco2_295(l)) - alog(xsco2_195(l)))
393     $                   /(295. - 195.)*(temp - 195.)
394                  sco2 = exp(sco2)
395               else
396                  sco2 = 0.
397               end if
398            else
399               if (xsco2_295(l) /= 0. .and. xsco2_370(l) /= 0.) then
400                  sco2 =  alog(xsco2_295(l))
401     $                 + (alog(xsco2_370(l)) - alog(xsco2_295(l)))
402     $                   /(370. - 295.)*(temp - 295.)
403                  sco2 = exp(sco2)
404               else
405                  sco2 = 0.
406               end if
407            end if
408
409!           optical depth
410
411            dt(i,l) = colinc(i)*rm(i)*sco2
412   
413!           production of o(1d) for wavelengths shorter than 167 nm
414
415            if (wc(l) >= 167.) then
416               sj(i,l,j_co2_o) = sco2*yieldco2(l)
417               sj(i,l,j_co2_o1d) = 0.
418            else
419               sj(i,l,j_co2_o) = 0.
420               sj(i,l,j_co2_o1d) = sco2*yieldco2(l)
421            end if
422         end do
423      end do
424
425      end subroutine setco2
426
427!==============================================================================
428
429      subroutine seto2(nd, nlayer, nw, wc, mopt, tlay, xso2_150,
430     $                 xso2_200, xso2_250, xso2_300, yieldo2, j_o2_o,
431     $                 j_o2_o1d, colinc, rm, dt, sj)
432
433!-----------------------------------------------------------------------------*
434!=  PURPOSE:                                                                 =*
435!=  Set up the O2 temperature-dependent cross-sections and optical depth     =*
436!-----------------------------------------------------------------------------*
437
438      implicit none
439
440!     input:
441
442      integer :: nd                                            ! number of photolysis rates
443      integer :: nlayer                                        ! number of vertical layers
444      integer :: nw                                            ! number of wavelength grid points
445      integer :: mopt                                          ! high-res/low-res switch
446      integer :: j_o2_o, j_o2_o1d                              ! photolysis indexes
447      real, dimension(nw)     :: wc                            ! central wavelength for each interval
448      real, dimension(nw)     :: xso2_150, xso2_200, xso2_250, ! o2 cross-sections (cm2)
449     $                           xso2_300
450      real, dimension(nw)     :: yieldo2                       ! o2 photodissociation yield
451      real, dimension(nlayer) :: tlay                          ! temperature (k)
452      real, dimension(nlayer) :: rm                            ! o2 mixing ratio
453      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
454
455!     output:
456
457      real, dimension(nlayer,nw)    :: dt                      !  optical depth
458      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
459
460!     local:
461
462      integer :: ilev, iw
463      real    :: temp
464      real    :: xso2, factor
465
466!     correction by factor if low-resolution in schumann-runge bands
467
468      if (mopt == 1) then
469         factor = 1.
470      else if (mopt == 2) then
471         factor = 0.8
472      end if
473
474!     calculate temperature dependance
475
476      do ilev = 1,nlayer
477         temp = max(tlay(ilev),150.)
478         temp = min(temp, 300.)
479         do iw = 1, nw-1
480            if (tlay(ilev) > 250.) then
481               xso2 = xso2_250(iw) + (xso2_300(iw) - xso2_250(iw))
482     $                              /(300. - 250.)*(temp - 250.)
483            else if (tlay(ilev) > 200.) then
484               xso2 = xso2_200(iw) + (xso2_250(iw) - xso2_200(iw))
485     $                              /(250. - 200.)*(temp - 200.)
486            else
487               xso2 = xso2_150(iw) + (xso2_200(iw) - xso2_150(iw))
488     $                              /(200. - 150.)*(temp - 150.)
489            end if
490
491            if (wc(iw) > 180. .and. wc(iw) < 200.) then
492               xso2 = xso2*factor
493            end if
494
495!     optical depth
496
497            dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso2
498
499!     production of o(1d) for wavelengths shorter than 175 nm
500
501            if (wc(iw) >= 175.) then
502               sj(ilev,iw,j_o2_o)   = xso2*yieldo2(iw)
503               sj(ilev,iw,j_o2_o1d) = 0.
504            else
505               sj(ilev,iw,j_o2_o)   = 0.
506               sj(ilev,iw,j_o2_o1d) = xso2*yieldo2(iw)
507            end if
508
509         end do
510      end do
511
512      end subroutine seto2
513
514!==============================================================================
515
516      subroutine seto3(nd, nlayer, nw, wc, tlay, xso3_218, xso3_298,
517     $                 j_o3_o, j_o3_o1d,
518     $                 colinc, rm, dt, sj)
519
520!-----------------------------------------------------------------------------*
521!=  PURPOSE:                                                                 =*
522!=  Set up the O3 temperature dependent cross-sections and optical depth     =*
523!-----------------------------------------------------------------------------*
524
525      implicit none
526
527!     input:
528
529      integer :: nd                                            ! number of photolysis rates
530      integer :: nlayer                                        ! number of vertical layers
531      integer :: nw                                            ! number of wavelength grid points
532      integer :: j_o3_o, j_o3_o1d                              ! photolysis indexes
533      real, dimension(nw)     :: wc                            ! central wavelength for each interval
534      real, dimension(nw)     :: xso3_218, xso3_298            ! o3 cross-sections (cm2)
535      real, dimension(nlayer) :: tlay                          ! temperature (k)
536      real, dimension(nlayer) :: rm                            ! o3 mixing ratio
537      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
538
539!     output:
540
541      real, dimension(nlayer,nw)    :: dt                      !  optical depth
542      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
543
544!     local:
545!
546      integer :: ilev, iw
547      real :: temp
548      real, dimension(nw) :: xso3(nw)
549      real, dimension(nw) :: qy1d                 ! quantum yield for o(1d) production
550      real :: q1, q2, a1, a2, a3
551
552      do ilev = 1, nlayer
553         temp = max(tlay(ilev), 218.)
554         temp = min(temp,298.)
555         do iw = 1, nw-1
556            xso3(iw) = xso3_218(iw) + (xso3_298(iw) - xso3_218(iw))
557     $                                /(298. - 218.) *(temp - 218.)
558
559!     optical depth
560
561            dt(ilev,iw) = colinc(ilev)*rm(ilev)*xso3(iw)
562
563         end do
564
565!     calculate quantum yield for o(1d) production (jpl 2006)
566
567         temp = max(tlay(ilev),200.)
568         temp = min(temp,320.)
569         do iw = 1, nw-1
570            if (wc(iw) <= 306.) then
571               qy1d(iw) = 0.90
572            else if (wc(iw) > 306. .and. wc(iw) < 328.) then
573               q1 = 1.
574               q2 = exp(-825.518/(0.695*temp))
575               a1 = (304.225 - wc(iw))/5.576
576               a2 = (314.957 - wc(iw))/6.601
577               a3 = (310.737 - wc(iw))/2.187
578               qy1d(iw) = (q1/(q1 + q2))*0.8036*exp(-(a1*a1*a1*a1))
579     $                  + (q2/(q1 + q2))*8.9061*(temp/300.)**2.
580     $                   *exp(-(a2*a2))
581     $                  + 0.1192*(temp/300.)**1.5*exp(-(a3*a3))
582     $                  + 0.0765
583            else if (wc(iw) >= 328. .and. wc(iw) <= 340.) then
584               qy1d(iw) = 0.08
585            else
586               qy1d(iw) = 0.
587            endif
588         end do
589         do iw = 1, nw-1
590            sj(ilev,iw,j_o3_o)   = xso3(iw)*(1. - qy1d(iw))
591            sj(ilev,iw,j_o3_o1d) = xso3(iw)*qy1d(iw)
592         end do
593      end do
594
595      end subroutine seto3
596
597!==============================================================================
598
599      subroutine seth2o2(nd, nlayer, nw, wc, tlay, xsh2o2, j_h2o2,
600     $                   colinc, rm, dt, sj)
601
602!-----------------------------------------------------------------------------*
603!=  PURPOSE:                                                                 =*
604!=  Set up the h2o2 temperature dependent cross-sections and optical depth   =*
605!-----------------------------------------------------------------------------*
606
607      implicit none
608
609!     input:
610
611      integer :: nd                                            ! number of photolysis rates
612      integer :: nlayer                                        ! number of vertical layers
613      integer :: nw                                            ! number of wavelength grid points
614      integer :: j_h2o2                                        ! photolysis index
615      real, dimension(nw)     :: wc                            ! central wavelength for each interval
616      real, dimension(nw)     :: xsh2o2                        ! h2o2 cross-sections (cm2)
617      real, dimension(nlayer) :: tlay                          ! temperature (k)
618      real, dimension(nlayer) :: rm                            ! h2o2 mixing ratio
619      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
620
621!     output:
622
623      real, dimension(nlayer,nw)    :: dt                      !  optical depth
624      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
625
626!     local:
627
628      integer :: ilev, iw
629      real    :: a0, a1, a2, a3, a4, a5, a6, a7
630      real    :: b0, b1, b2, b3, b4
631      real    :: lambda, suma, sumb, chi, temp, xs
632
633      A0 = 6.4761E+04           
634      A1 = -9.2170972E+02       
635      A2 = 4.535649             
636      A3 = -4.4589016E-03       
637      A4 = -4.035101E-05         
638      A5 = 1.6878206E-07
639      A6 = -2.652014E-10
640      A7 = 1.5534675E-13
641
642      B0 = 6.8123E+03
643      B1 = -5.1351E+01
644      B2 = 1.1522E-01
645      B3 = -3.0493E-05
646      B4 = -1.0924E-07
647
648!     temperature dependance: jpl 2006
649
650      do ilev = 1,nlayer
651         temp = min(max(tlay(ilev),200.),400.)           
652         chi = 1./(1. + exp(-1265./temp))
653         do iw = 1, nw-1
654            if ((wc(iw) >= 260.) .and. (wc(iw) < 350.)) then
655               lambda = wc(iw)
656               sumA = ((((((A7*lambda + A6)*lambda + A5)*lambda +
657     $                      A4)*lambda +A3)*lambda + A2)*lambda +
658     $                      A1)*lambda + A0
659               sumB = (((B4*lambda + B3)*lambda + B2)*lambda +
660     $                   B1)*lambda + B0
661               xs = (chi*sumA + (1. - chi)*sumB)*1.e-21
662               sj(ilev,iw,j_h2o2) = xs
663            else
664               sj(ilev,iw,j_h2o2) = xsh2o2(iw)
665            end if
666
667!     optical depth
668
669            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_h2o2)
670         end do
671      end do
672
673      end subroutine seth2o2
674
675!==============================================================================
676
677      subroutine setno2(nd, nlayer, nw, wc, tlay, xsno2, xsno2_220,
678     $                  xsno2_294, yldno2_248, yldno2_298, j_no2,
679     $                  colinc, rm, dt, sj)
680
681!-----------------------------------------------------------------------------*
682!=  PURPOSE:                                                                 =*
683!=  Set up the no2 temperature-dependent cross-sections and optical depth    =*
684!-----------------------------------------------------------------------------*
685
686      implicit none
687
688!     input:
689
690      integer :: nd                                            ! number of photolysis rates
691      integer :: nlayer                                        ! number of vertical layers
692      integer :: nw                                            ! number of wavelength grid points
693      integer :: j_no2                                         ! photolysis index
694      real, dimension(nw)     :: wc                            ! central wavelength for each interval
695      real, dimension(nw) :: xsno2, xsno2_220, xsno2_294       ! no2 absorption cross-section at 220-294 k (cm2)
696      real, dimension(nw) :: yldno2_248, yldno2_298            ! no2 quantum yield at 248-298 k
697      real, dimension(nlayer) :: tlay                          ! temperature (k)
698      real, dimension(nlayer) :: rm                            ! no2 mixing ratio
699      real, dimension(nlayer) :: colinc                        ! air column increment (molecule.cm-2)
700
701!     output:
702
703      real, dimension(nlayer,nw)    :: dt                      !  optical depth
704      real, dimension(nlayer,nw,nd) :: sj                      !  cross-section array (cm2)
705
706!     local:
707
708      integer :: ilev, iw
709      real    :: temp, qy
710
711!     temperature dependance: jpl 2006
712
713      do ilev = 1,nlayer
714         temp = max(220.,min(tlay(ilev),294.))
715         do iw = 1, nw - 1
716            if (wc(iw) < 238.) then
717               sj(ilev,iw,j_no2) = xsno2(iw)
718            else
719               sj(ilev,iw,j_no2) = xsno2_220(iw)
720     $                           + (xsno2_294(iw) - xsno2_220(iw))
721     $                            /(294. - 220.)*(temp - 220.)
722            end if
723
724!     optical depth
725
726            dt(ilev,iw) = colinc(ilev)*rm(ilev)*sj(ilev,iw,j_no2)
727         end do
728      end do
729
730!     quantum yield: jpl 2006
731
732      do ilev = 1,nlayer
733         temp = max(248.,min(tlay(ilev),298.))
734         do iw = 1, nw - 1
735            qy = yldno2_248(iw) + (yldno2_298(iw) - yldno2_248(iw))
736     $                           /(298. - 248.)*(temp - 248.)
737            sj(ilev,iw,j_no2) = sj(ilev,iw,j_no2)*qy
738         end do
739      end do
740
741      end subroutine setno2
742
743!==============================================================================
744
745      subroutine setaer(nlev,zalt,tau,nw,dtaer,omaer,gaer)
746
747!-----------------------------------------------------------------------------*
748!=  PURPOSE:                                                                 =*
749!=  Set aerosol properties for each specified altitude layer.  Properties    =*
750!=  may be wavelength dependent.                                             =*
751!-----------------------------------------------------------------------------*
752
753      implicit none
754
755!     input
756
757      integer :: nlev                                ! number of vertical levels
758      integer :: nw                                  ! number of wavelength grid points
759      real, dimension(nlev) :: zalt                  ! altitude (km)
760      real :: tau                                    ! integrated aerosol optical depth at the surface
761
762!     output
763
764      real, dimension(nlev,nw) :: dtaer              ! aerosol optical depth
765      real, dimension(nlev,nw) :: omaer              ! aerosol single scattering albedo
766      real, dimension(nlev,nw) :: gaer               ! aerosol asymmetry parameter
767
768!     local
769
770      integer :: ilev, iw
771      real, dimension(nlev) :: aer                   ! dust extinction
772      real :: omega, g, scaleh, gamma
773      real :: dz, tautot, q0
774
775      omega  = 0.622 ! single scattering albedo : wolff et al.(2010) at 258 nm
776      g      = 0.88  ! asymmetry factor : mateshvili et al. (2007) at 210 nm
777      scaleh = 10.   ! scale height (km)
778      gamma  = 0.03  ! conrath parameter
779
780      dtaer(:,:) = 0.
781      omaer(:,:) = 0.
782      gaer(:,:)  = 0.
783
784!     optical depth profile:
785
786      tautot = 0.
787      do ilev = 1, nlev-1
788         dz = zalt(ilev+1) - zalt(ilev)
789         tautot = tautot + exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
790      end do
791     
792      q0 = tau/tautot
793      do ilev = 1, nlev-1
794         dz = zalt(ilev+1) - zalt(ilev)
795         dtaer(ilev,:) = q0*exp(gamma*(1. - exp(zalt(ilev)/scaleh)))*dz
796         omaer(ilev,:) = omega
797         gaer(ilev,:)  = g
798      end do
799
800      end subroutine setaer
801
802!==============================================================================
803
804      subroutine setcld(nlev,nw,dtcld,omcld,gcld)
805
806!-----------------------------------------------------------------------------*
807!=  PURPOSE:                                                                 =*
808!=  Set cloud properties for each specified altitude layer.  Properties      =*
809!=  may be wavelength dependent.                                             =*
810!-----------------------------------------------------------------------------*
811
812      implicit none
813
814!     input
815
816      integer :: nlev                                ! number of vertical levels
817      integer :: nw                                  ! number of wavelength grid points
818
819!     output
820
821      real, dimension(nlev,nw) :: dtcld            ! cloud optical depth
822      real, dimension(nlev,nw) :: omcld            ! cloud single scattering albedo
823      real, dimension(nlev,nw) :: gcld             ! cloud asymmetry parameter
824
825!     local
826
827      integer :: ilev, iw
828
829!     dtcld : optical depth
830!     omcld : single scattering albedo
831!     gcld  : asymmetry factor
832
833      dtcld(:,:) = 0.
834
835      do ilev = 1,nlev-1
836         do iw = 1,nw-1
837            dtcld(ilev,iw) = 0.       ! no clouds for the moment
838            omcld(ilev,iw) = 0.99
839            gcld(ilev,iw)  = 0.85
840         end do
841      end do
842
843      end subroutine setcld
844
845!==============================================================================
846
847      subroutine sphers(nlev, z, zen, dsdh, nid)
848
849!-----------------------------------------------------------------------------*
850!=  PURPOSE:                                                                 =*
851!=  Calculate slant path over vertical depth ds/dh in spherical geometry.    =*
852!=  Calculation is based on:  A.Dahlback, and K.Stamnes, A new spheric model =*
853!=  for computing the radiation field available for photolysis and heating   =*
854!=  at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B)  =*
855!-----------------------------------------------------------------------------*
856!=  PARAMETERS:                                                              =*
857!=  NZ      - INTEGER, number of specified altitude levels in the working (I)=*
858!=            grid                                                           =*
859!=  Z       - REAL, specified altitude working grid (km)                  (I)=*
860!=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
861!=  DSDH    - REAL, slant path of direct beam through each layer crossed  (O)=*
862!=            when travelling from the top of the atmosphere to layer i;     =*
863!=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
864!=  NID     - INTEGER, number of layers crossed by the direct beam when   (O)=*
865!=            travelling from the top of the atmosphere to layer i;          =*
866!=            NID(i), i = 0..NZ-1                                            =*
867!-----------------------------------------------------------------------------*
868
869      implicit none
870
871! input
872
873      integer, intent(in) :: nlev
874      real, dimension(nlev), intent(in) :: z
875      real, intent(in) :: zen
876
877! output
878
879      INTEGER nid(0:nlev)
880      REAL dsdh(0:nlev,nlev)
881
882! more program constants
883
884      REAL re, ze(nlev)
885      REAL  dr
886      real radius
887      parameter (radius = 3393.)
888
889! local
890
891      real :: pi, zenrad, rpsinz, rj, rjp1, dsj, dhj, ga, gb, sm
892      integer :: i, j, k, id, nlay
893
894      REAL zd(0:nlev-1)
895
896!-----------------------------------------------------------------------------
897
898      pi = acos(-1.0)
899      dr = pi/180.
900      zenrad = zen*dr
901
902! number of layers:
903
904      nlay = nlev - 1
905
906! include the elevation above sea level to the radius of Mars:
907
908      re = radius + z(1)
909
910! correspondingly z changed to the elevation above Mars surface:
911
912      DO k = 1, nlev
913         ze(k) = z(k) - z(1)
914      END DO
915
916! inverse coordinate of z
917
918      zd(0) = ze(nlev)
919      DO k = 1, nlay
920        zd(k) = ze(nlev - k)
921      END DO
922
923! initialise dsdh(i,j), nid(i)
924
925      nid(:) = 0.
926      dsdh(:,:) = 0.
927
928! calculate ds/dh of every layer
929
930      do i = 0,nlay
931         rpsinz = (re + zd(i))*sin(zenrad)
932 
933        IF ( (zen .GT. 90.0) .AND. (rpsinz .LT. re) ) THEN
934           nid(i) = -1
935        ELSE
936
937! Find index of layer in which the screening height lies
938
939           id = i
940           if (zen > 90.) then
941              do j = 1,nlay
942                 IF( (rpsinz .LT. ( zd(j-1) + re ) ) .AND.
943     $               (rpsinz .GE. ( zd(j) + re )) ) id = j
944              end do
945           end if
946 
947           do j = 1,id
948              sm = 1.0
949              IF (j .EQ. id .AND. id .EQ. i .AND. zen .GT. 90.0)
950     $            sm = -1.0
951 
952              rj = re + zd(j-1)
953              rjp1 = re + zd(j)
954 
955              dhj = zd(j-1) - zd(j)
956 
957              ga = rj*rj - rpsinz*rpsinz
958              gb = rjp1*rjp1 - rpsinz*rpsinz
959
960              ga = max(ga, 0.)
961              gb = max(gb, 0.)
962 
963              IF (id.GT.i .AND. j.EQ.id) THEN
964                 dsj = sqrt(ga)
965              ELSE
966                 dsj = sqrt(ga) - sm*sqrt(gb)
967              END IF
968              dsdh(i,j) = dsj/dhj
969            end do
970            nid(i) = id
971         end if
972      end do ! i = 0,nlay
973
974      end subroutine sphers
975
976!==============================================================================
977
978      SUBROUTINE rtlink(nlev, nw, iw, ag, zen, dsdh, nid, dtrl,
979     $                  dagas, dtcld, omcld, gcld, dtaer, omaer, gaer,
980     $                  edir, edn, eup, fdir, fdn, fup)
981
982      implicit none
983
984! input
985
986      integer, intent(in) :: nlev, nw, iw                       ! number of wavelength grid points
987      REAL ag
988      REAL zen
989      REAL dsdh(0:nlev,nlev)
990      INTEGER nid(0:nlev)
991
992      REAL dtrl(nlev,nw)
993      REAL dagas(nlev,nw)
994      REAL dtcld(nlev,nw), omcld(nlev,nw), gcld(nlev,nw)
995      REAL dtaer(nlev,nw), omaer(nlev,nw), gaer(nlev,nw)
996
997! output
998
999      REAL edir(nlev), edn(nlev), eup(nlev)
1000      REAL fdir(nlev), fdn(nlev), fup(nlev)
1001
1002! local:
1003
1004      REAL dt(nlev), om(nlev), g(nlev)
1005      REAL dtabs,dtsct,dscld,dsaer,dacld,daaer
1006      INTEGER i, ii
1007      real, parameter :: largest = 1.e+36
1008
1009! specific two ps2str
1010
1011      REAL ediri(nlev), edni(nlev), eupi(nlev)
1012      REAL fdiri(nlev), fdni(nlev), fupi(nlev)
1013
1014      logical, save :: delta = .true.
1015
1016!$OMP THREADPRIVATE(delta)
1017
1018!_______________________________________________________________________
1019
1020! initialize:
1021
1022      do i = 1, nlev
1023         fdir(i) = 0.
1024         fup(i)  = 0.
1025         fdn(i)  = 0.
1026         edir(i) = 0.
1027         eup(i)  = 0.
1028         edn(i)  = 0.
1029      end do
1030
1031      do i = 1, nlev - 1
1032         dscld = dtcld(i,iw)*omcld(i,iw)
1033         dacld = dtcld(i,iw)*(1.-omcld(i,iw))
1034
1035         dsaer = dtaer(i,iw)*omaer(i,iw)
1036         daaer = dtaer(i,iw)*(1.-omaer(i,iw))
1037
1038         dtsct = dtrl(i,iw) + dscld + dsaer
1039         dtabs = dagas(i,iw) + dacld + daaer
1040
1041         dtabs = amax1(dtabs,1./largest)
1042         dtsct = amax1(dtsct,1./largest)
1043
1044!     invert z-coordinate:
1045
1046         ii = nlev - i
1047         dt(ii) = dtsct + dtabs
1048         om(ii) = dtsct/(dtsct + dtabs)
1049         IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
1050         g(ii) = (gcld(i,iw)*dscld +
1051     $            gaer(i,iw)*dsaer)/dtsct
1052      end do
1053
1054!     call rt routine:
1055
1056      call ps2str(nlev, zen, ag, dt, om, g,
1057     $            dsdh, nid, delta,
1058     $            fdiri, fupi, fdni, ediri, eupi, edni)
1059
1060!     output (invert z-coordinate)
1061
1062      do i = 1, nlev
1063         ii = nlev - i + 1
1064         fdir(i) = fdiri(ii)
1065         fup(i) = fupi(ii)
1066         fdn(i) = fdni(ii)
1067         edir(i) = ediri(ii)
1068         eup(i) = eupi(ii)
1069         edn(i) = edni(ii)
1070      end do
1071
1072      end subroutine rtlink
1073
1074*=============================================================================*
1075
1076      subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
1077     $                  dsdh, nid, delta,
1078     $                  fdr, fup, fdn, edr, eup, edn)
1079
1080!-----------------------------------------------------------------------------*
1081!=  PURPOSE:                                                                 =*
1082!=  Solve two-stream equations for multiple layers.  The subroutine is based =*
1083!=  on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
1084!=  It contains 9 two-stream methods to choose from.  A pseudo-spherical     =*
1085!=  correction has also been added.                                          =*
1086!-----------------------------------------------------------------------------*
1087!=  PARAMETERS:                                                              =*
1088!=  NLEVEL  - INTEGER, number of specified altitude levels in the working (I)=*
1089!=            grid                                                           =*
1090!=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
1091!=  RSFC    - REAL, surface albedo at current wavelength                  (I)=*
1092!=  TAUU    - REAL, unscaled optical depth of each layer                  (I)=*
1093!=  OMU     - REAL, unscaled single scattering albedo of each layer       (I)=*
1094!=  GU      - REAL, unscaled asymmetry parameter of each layer            (I)=*
1095!=  DSDH    - REAL, slant path of direct beam through each layer crossed  (I)=*
1096!=            when travelling from the top of the atmosphere to layer i;     =*
1097!=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
1098!=  NID     - INTEGER, number of layers crossed by the direct beam when   (I)=*
1099!=            travelling from the top of the atmosphere to layer i;          =*
1100!=            NID(i), i = 0..NZ-1                                            =*
1101!=  DELTA   - LOGICAL, switch to use delta-scaling                        (I)=*
1102!=            .TRUE. -> apply delta-scaling                                  =*
1103!=            .FALSE.-> do not apply delta-scaling                           =*
1104!=  FDR     - REAL, contribution of the direct component to the total     (O)=*
1105!=            actinic flux at each altitude level                            =*
1106!=  FUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
1107!=            the total actinic flux at each altitude level                  =*
1108!=  FDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
1109!=            the total actinic flux at each altitude level                  =*
1110!=  EDR     - REAL, contribution of the direct component to the total     (O)=*
1111!=            spectral irradiance at each altitude level                     =*
1112!=  EUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
1113!=            the total spectral irradiance at each altitude level           =*
1114!=  EDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
1115*=            the total spectral irradiance at each altitude level           =*
1116!-----------------------------------------------------------------------------*
1117
1118      implicit none
1119
1120! input:
1121
1122      INTEGER nlev
1123      REAL zen, rsfc
1124      REAL tauu(nlev), omu(nlev), gu(nlev)
1125      REAL dsdh(0:nlev,nlev)
1126      INTEGER nid(0:nlev)
1127      LOGICAL delta
1128
1129! output:
1130
1131      REAL fup(nlev),fdn(nlev),fdr(nlev)
1132      REAL eup(nlev),edn(nlev),edr(nlev)
1133
1134! local:
1135
1136      REAL tausla(0:nlev), tauc(0:nlev)
1137      REAL mu2(0:nlev), mu, sum
1138
1139! internal coefficients and matrix
1140
1141      REAL lam(nlev),taun(nlev),bgam(nlev)
1142      REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
1143      REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
1144      REAL mu1(nlev)
1145      INTEGER row
1146      REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev)
1147
1148! other:
1149
1150      REAL pifs, fdn0
1151      REAL gi(nlev), omi(nlev), tempg
1152      REAL f, g, om
1153      REAL gam1, gam2, gam3, gam4
1154      real, parameter :: largest = 1.e+36
1155      real, parameter :: precis = 1.e-7
1156
1157! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1158! in delta-function, modified quadrature, hemispheric constant,
1159! Hybrid modified Eddington-delta function metods, p633,Table1.
1160! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1161! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
1162! uncomment the following two lines and the appropriate statements further
1163! down.
1164!     REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
1165!    >     BETA1, BETAn, amu1, subd
1166
1167      REAL expon, expon0, expon1, divisr, temp, up, dn
1168      REAL ssfc
1169      INTEGER nlayer, mrows, lev
1170
1171      INTEGER i, j
1172
1173! Some additional program constants:
1174
1175      real pi, dr
1176      REAL eps
1177      PARAMETER (eps = 1.E-3)
1178!_______________________________________________________________________
1179
1180! MU = cosine of solar zenith angle
1181! RSFC = surface albedo
1182! TAUU =  unscaled optical depth of each layer
1183! OMU  =  unscaled single scattering albedo
1184! GU   =  unscaled asymmetry factor
1185! KLEV = max dimension of number of layers in atmosphere
1186! NLAYER = number of layers in the atmosphere
1187! NLEVEL = nlayer + 1 = number of levels
1188
1189! initial conditions:  pi*solar flux = 1;  diffuse incidence = 0
1190
1191      pifs = 1.     
1192      fdn0 = 0.
1193
1194      nlayer = nlev - 1
1195
1196      pi = acos(-1.)
1197      dr = pi/180.
1198      mu = COS(zen*dr)
1199
1200!************* compute coefficients for each layer:
1201! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
1202! EXPON0 = calculation of e when TAU is zero
1203! EXPON1 = calculation of e when TAU is TAUN
1204! CUP and CDN = calculation when TAU is zero
1205! CUPTN and CDNTN = calc. when TAU is TAUN
1206! DIVISR = prevents division by zero
1207
1208      do j = 0, nlev
1209         tauc(j) = 0.
1210         tausla(j) = 0.
1211         mu2(j) = 1./SQRT(largest)
1212      end do
1213
1214      IF (.NOT. delta) THEN
1215         DO i = 1, nlayer
1216            gi(i) = gu(i)
1217            omi(i) = omu(i)
1218            taun(i) = tauu(i)
1219         END DO
1220      ELSE
1221
1222! delta-scaling. Have to be done for delta-Eddington approximation,
1223! delta discrete ordinate, Practical Improved Flux Method, delta function,
1224! and Hybrid modified Eddington-delta function methods approximations
1225
1226         DO i = 1, nlayer
1227            f = gu(i)*gu(i)
1228            gi(i) = (gu(i) - f)/(1 - f)
1229            omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)       
1230            taun(i) = (1 - omu(i)*f)*tauu(i)
1231         END DO
1232      END IF
1233
1234! calculate slant optical depth at the top of the atmosphere when zen>90.
1235! in this case, higher altitude of the top layer is recommended.
1236
1237      IF (zen .GT. 90.0) THEN
1238         IF (nid(0) .LT. 0) THEN
1239            tausla(0) = largest
1240         ELSE
1241            sum = 0.0
1242            DO j = 1, nid(0)
1243               sum = sum + 2.*taun(j)*dsdh(0,j)
1244            END DO
1245            tausla(0) = sum
1246         END IF
1247      END IF
1248
1249      DO 11, i = 1, nlayer
1250          g = gi(i)
1251          om = omi(i)
1252          tauc(i) = tauc(i-1) + taun(i)
1253
1254! stay away from 1 by precision.  For g, also stay away from -1
1255
1256          tempg = AMIN1(abs(g),1. - precis)
1257          g = SIGN(tempg,g)
1258          om = AMIN1(om,1.-precis)
1259
1260! calculate slant optical depth
1261             
1262          IF (nid(i) .LT. 0) THEN
1263             tausla(i) = largest
1264          ELSE
1265             sum = 0.0
1266             DO j = 1, MIN(nid(i),i)
1267                sum = sum + taun(j)*dsdh(i,j)
1268             END DO
1269             DO j = MIN(nid(i),i)+1,nid(i)
1270                sum = sum + 2.*taun(j)*dsdh(i,j)
1271             END DO
1272             tausla(i) = sum
1273             IF (tausla(i) .EQ. tausla(i-1)) THEN
1274                mu2(i) = SQRT(largest)
1275             ELSE
1276                mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
1277                mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
1278     $                     mu2(i) )
1279             END IF
1280          END IF
1281
1282!** the following gamma equations are from pg 16,289, Table 1
1283!** save mu1 for each approx. for use in converting irradiance to actinic flux
1284
1285! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
1286
1287c       gam1 =  (7. - om*(4. + 3.*g))/4.
1288c       gam2 = -(1. - om*(4. - 3.*g))/4.
1289c       gam3 = (2. - 3.*g*mu)/4.
1290c       gam4 = 1. - gam3
1291c       mu1(i) = 0.5
1292
1293* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
1294
1295c          gam1 = 1.7320508*(2. - om*(1. + g))/2.
1296c          gam2 = 1.7320508*om*(1. - g)/2.
1297c          gam3 = (1. - 1.7320508*g*mu)/2.
1298c          gam4 = 1. - gam3
1299c          mu1(i) = 1./sqrt(3.)
1300         
1301* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
1302
1303           gam1 = 2. - om*(1. + g)
1304           gam2 = om*(1. - g)
1305           gam3 = (2. - g*mu)/4.
1306           gam4 = 1. - gam3
1307           mu1(i) = 0.5
1308
1309* PIFM  (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
1310c         GAM1 = 0.25*(8. - OM*(5. + 3.*G))
1311c         GAM2 = 0.75*OM*(1.-G)
1312c         GAM3 = 0.25*(2.-3.*G*MU)
1313c         GAM4 = 1. - GAM3
1314c         mu1(i) = 0.5
1315
1316* delta discrete ordinates  (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
1317c         GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
1318c         GAM2 = 0.5*1.7320508*OM*(1.-G)
1319c         GAM3 = 0.5*(1.-1.7320508*G*MU)
1320c         GAM4 = 1. - GAM3
1321c         mu1(i) = 1./sqrt(3.)
1322
1323* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1324* in delta-function, modified quadrature, hemispheric constant,
1325* Hybrid modified Eddington-delta function metods, p633,Table1.
1326* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1327* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
1328c      YLM0 = 2.
1329c      YLM2 = -3.*G*MU
1330c      YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
1331c      YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
1332c     YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
1333c    *+429.*MU**6)
1334c     YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
1335c    *-25740.*MU**6+12155.*MU**8)
1336c     YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
1337c    *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
1338c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1339c      YLMS=0.25*YLMS
1340c      BETA0 = YLMS
1341c
1342c         amu1=1./1.7320508
1343c      YLM0 = 2.
1344c      YLM2 = -3.*G*amu1
1345c      YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
1346c      YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
1347c     YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
1348c    *+429.*amu1**6)
1349c     YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
1350c    *-25740.*amu1**6+12155.*amu1**8)
1351c     YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
1352c    *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
1353c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1354c      YLMS=0.25*YLMS
1355c      BETA1 = YLMS
1356c
1357c         BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
1358c    *-0.045776*G**7)
1359
1360
1361* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
1362c         subd=4.*(1.-G*G*(1.-MU))
1363c         GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
1364c         GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
1365c         GAM3 = BETA0
1366c         GAM4 = 1. - GAM3
1367c         mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)
1368
1369*****
1370* delta function  (Meador, and Weaver, 1980, JAS, 37, 630):
1371c         GAM1 = (1. - OM*(1. - beta0))/MU
1372c         GAM2 = OM*BETA0/MU
1373c         GAM3 = BETA0
1374c         GAM4 = 1. - GAM3
1375c         mu1(i) = mu
1376*****
1377* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
1378c         GAM1 = 1.7320508*(1. - OM*(1. - beta1))
1379c         GAM2 = 1.7320508*OM*beta1
1380c         GAM3 = BETA0
1381c         GAM4 = 1. - GAM3
1382c         mu1(i) = 1./sqrt(3.)
1383
1384* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
1385c         GAM1 = 2.*(1. - OM*(1. - betan))
1386c         GAM2 = 2.*OM*BETAn
1387c         GAM3 = BETA0
1388c         GAM4 = 1. - GAM3
1389c         mu1(i) = 0.5
1390
1391*****
1392
1393* lambda = pg 16,290 equation 21
1394* big gamma = pg 16,290 equation 22
1395* if gam2 = 0., then bgam = 0.
1396
1397         lam(i) = sqrt(gam1*gam1 - gam2*gam2)
1398
1399         IF (gam2 .NE. 0.) THEN
1400            bgam(i) = (gam1 - lam(i))/gam2
1401         ELSE
1402            bgam(i) = 0.
1403         END IF
1404
1405         expon = EXP(-lam(i)*taun(i))
1406
1407* e1 - e4 = pg 16,292 equation 44
1408         
1409         e1(i) = 1. + bgam(i)*expon
1410         e2(i) = 1. - bgam(i)*expon
1411         e3(i) = bgam(i) + expon
1412         e4(i) = bgam(i) - expon
1413
1414* the following sets up for the C equations 23, and 24
1415* found on page 16,290
1416* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
1417* which is approx equiv to shifting MU by 0.5*EPS* (MU)**3
1418
1419         expon0 = EXP(-tausla(i-1))
1420         expon1 = EXP(-tausla(i))
1421         
1422         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
1423         temp = AMAX1(eps,abs(divisr))
1424         divisr = SIGN(temp,divisr)
1425
1426         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
1427         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
1428         
1429* cup and cdn are when tau is equal to zero
1430* cuptn and cdntn are when tau is equal to taun
1431
1432         cup(i) = up*expon0
1433         cdn(i) = dn*expon0
1434         cuptn(i) = up*expon1
1435         cdntn(i) = dn*expon1
1436 
1437   11 CONTINUE
1438
1439***************** set up matrix ******
1440* ssfc = pg 16,292 equation 37  where pi Fs is one (unity).
1441
1442      ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs
1443
1444* MROWS = the number of rows in the matrix
1445
1446      mrows = 2*nlayer     
1447     
1448* the following are from pg 16,292  equations 39 - 43.
1449* set up first row of matrix:
1450
1451      i = 1
1452      a(1) = 0.
1453      b(1) = e1(i)
1454      d(1) = -e2(i)
1455      e(1) = fdn0 - cdn(i)
1456
1457      row=1
1458
1459* set up odd rows 3 thru (MROWS - 1):
1460
1461      i = 0
1462      DO 20, row = 3, mrows - 1, 2
1463         i = i + 1
1464         a(row) = e2(i)*e3(i) - e4(i)*e1(i)
1465         b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
1466         d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
1467         e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
1468     $        e1(i)*(cdntn(i) - cdn(i + 1))
1469   20 CONTINUE
1470
1471* set up even rows 2 thru (MROWS - 2):
1472
1473      i = 0
1474      DO 30, row = 2, mrows - 2, 2
1475         i = i + 1
1476         a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
1477         b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
1478         d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
1479         e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
1480     $        (cdn(i + 1) - cdntn(i))*e4(i + 1)
1481   30 CONTINUE
1482
1483* set up last row of matrix at MROWS:
1484
1485      row = mrows
1486      i = nlayer
1487     
1488      a(row) = e1(i) - rsfc*e3(i)
1489      b(row) = e2(i) - rsfc*e4(i)
1490      d(row) = 0.
1491      e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
1492
1493* solve tri-diagonal matrix:
1494
1495      CALL tridiag(a, b, d, e, y, mrows)
1496
1497**** unfold solution of matrix, compute output fluxes:
1498
1499      row = 1
1500      lev = 1
1501      j = 1
1502     
1503* the following equations are from pg 16,291  equations 31 & 32
1504
1505      fdr(lev) = EXP( -tausla(0) )
1506      edr(lev) = mu * fdr(lev)
1507      edn(lev) = fdn0
1508      eup(lev) =  y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
1509      fdn(lev) = edn(lev)/mu1(lev)
1510      fup(lev) = eup(lev)/mu1(lev)
1511
1512      DO 60, lev = 2, nlayer + 1
1513         fdr(lev) = EXP(-tausla(lev-1))
1514         edr(lev) =  mu *fdr(lev)
1515         edn(lev) =  y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
1516         eup(lev) =  y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
1517         fdn(lev) = edn(lev)/mu1(j)
1518         fup(lev) = eup(lev)/mu1(j)
1519
1520         row = row + 2
1521         j = j + 1
1522   60 CONTINUE
1523
1524      end subroutine ps2str
1525
1526*=============================================================================*
1527
1528      subroutine tridiag(a,b,c,r,u,n)
1529
1530!_______________________________________________________________________
1531! solves tridiagonal system.  From Numerical Recipies, p. 40
1532!_______________________________________________________________________
1533
1534      IMPLICIT NONE
1535
1536! input:
1537
1538      INTEGER n
1539      REAL a, b, c, r
1540      DIMENSION a(n),b(n),c(n),r(n)
1541
1542! output:
1543
1544      REAL u
1545      DIMENSION u(n)
1546
1547! local:
1548
1549      INTEGER j
1550
1551      REAL bet, gam
1552      DIMENSION gam(n)
1553!_______________________________________________________________________
1554
1555      IF (b(1) .EQ. 0.) STOP 1001
1556      bet   = b(1)
1557      u(1) = r(1)/bet
1558      DO 11, j = 2, n   
1559         gam(j) = c(j - 1)/bet
1560         bet = b(j) - a(j)*gam(j)
1561         IF (bet .EQ. 0.) STOP 2002
1562         u(j) = (r(j) - a(j)*u(j - 1))/bet
1563   11 CONTINUE
1564      DO 12, j = n - 1, 1, -1 
1565         u(j) = u(j) - gam(j + 1)*u(j + 1)
1566   12 CONTINUE
1567!_______________________________________________________________________
1568
1569      end subroutine tridiag
1570
1571      end subroutine photolysis_online
1572     
1573      end module photolysis_online_mod
Note: See TracBrowser for help on using the repository browser.