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

Last change on this file since 2833 was 2615, checked in by romain.vande, 3 years ago

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