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

Last change on this file since 3158 was 2968, checked in by flefevre, 3 years ago

on-line photolysis: bug fix in the vertical grid used by the UV code at the top of the atmosphere. This correction leads to significant changes in the chemistry with low top level PCM configurations. No signicant change is observed with PCM configurations including the thermosphere.

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