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

Last change on this file since 2599 was 2489, checked in by flefevre, 4 years ago

Photolysis : minor cosmetic changes

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!_______________________________________________________________________
983
984! initialize:
985
986      do i = 1, nlev
987         fdir(i) = 0.
988         fup(i)  = 0.
989         fdn(i)  = 0.
990         edir(i) = 0.
991         eup(i)  = 0.
992         edn(i)  = 0.
993      end do
994
995      do i = 1, nlev - 1
996         dscld = dtcld(i,iw)*omcld(i,iw)
997         dacld = dtcld(i,iw)*(1.-omcld(i,iw))
998
999         dsaer = dtaer(i,iw)*omaer(i,iw)
1000         daaer = dtaer(i,iw)*(1.-omaer(i,iw))
1001
1002         dtsct = dtrl(i,iw) + dscld + dsaer
1003         dtabs = dagas(i,iw) + dacld + daaer
1004
1005         dtabs = amax1(dtabs,1./largest)
1006         dtsct = amax1(dtsct,1./largest)
1007
1008!     invert z-coordinate:
1009
1010         ii = nlev - i
1011         dt(ii) = dtsct + dtabs
1012         om(ii) = dtsct/(dtsct + dtabs)
1013         IF(dtsct .EQ. 1./largest) om(ii) = 1./largest
1014         g(ii) = (gcld(i,iw)*dscld +
1015     $            gaer(i,iw)*dsaer)/dtsct
1016      end do
1017
1018!     call rt routine:
1019
1020      call ps2str(nlev, zen, ag, dt, om, g,
1021     $            dsdh, nid, delta,
1022     $            fdiri, fupi, fdni, ediri, eupi, edni)
1023
1024!     output (invert z-coordinate)
1025
1026      do i = 1, nlev
1027         ii = nlev - i + 1
1028         fdir(i) = fdiri(ii)
1029         fup(i) = fupi(ii)
1030         fdn(i) = fdni(ii)
1031         edir(i) = ediri(ii)
1032         eup(i) = eupi(ii)
1033         edn(i) = edni(ii)
1034      end do
1035
1036      end subroutine rtlink
1037
1038*=============================================================================*
1039
1040      subroutine ps2str(nlev,zen,rsfc,tauu,omu,gu,
1041     $                  dsdh, nid, delta,
1042     $                  fdr, fup, fdn, edr, eup, edn)
1043
1044!-----------------------------------------------------------------------------*
1045!=  PURPOSE:                                                                 =*
1046!=  Solve two-stream equations for multiple layers.  The subroutine is based =*
1047!=  on equations from:  Toon et al., J.Geophys.Res., v94 (D13), Nov 20, 1989.=*
1048!=  It contains 9 two-stream methods to choose from.  A pseudo-spherical     =*
1049!=  correction has also been added.                                          =*
1050!-----------------------------------------------------------------------------*
1051!=  PARAMETERS:                                                              =*
1052!=  NLEVEL  - INTEGER, number of specified altitude levels in the working (I)=*
1053!=            grid                                                           =*
1054!=  ZEN     - REAL, solar zenith angle (degrees)                          (I)=*
1055!=  RSFC    - REAL, surface albedo at current wavelength                  (I)=*
1056!=  TAUU    - REAL, unscaled optical depth of each layer                  (I)=*
1057!=  OMU     - REAL, unscaled single scattering albedo of each layer       (I)=*
1058!=  GU      - REAL, unscaled asymmetry parameter of each layer            (I)=*
1059!=  DSDH    - REAL, slant path of direct beam through each layer crossed  (I)=*
1060!=            when travelling from the top of the atmosphere to layer i;     =*
1061!=            DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1                            =*
1062!=  NID     - INTEGER, number of layers crossed by the direct beam when   (I)=*
1063!=            travelling from the top of the atmosphere to layer i;          =*
1064!=            NID(i), i = 0..NZ-1                                            =*
1065!=  DELTA   - LOGICAL, switch to use delta-scaling                        (I)=*
1066!=            .TRUE. -> apply delta-scaling                                  =*
1067!=            .FALSE.-> do not apply delta-scaling                           =*
1068!=  FDR     - REAL, contribution of the direct component to the total     (O)=*
1069!=            actinic flux at each altitude level                            =*
1070!=  FUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
1071!=            the total actinic flux at each altitude level                  =*
1072!=  FDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
1073!=            the total actinic flux at each altitude level                  =*
1074!=  EDR     - REAL, contribution of the direct component to the total     (O)=*
1075!=            spectral irradiance at each altitude level                     =*
1076!=  EUP     - REAL, contribution of the diffuse upwelling component to    (O)=*
1077!=            the total spectral irradiance at each altitude level           =*
1078!=  EDN     - REAL, contribution of the diffuse downwelling component to  (O)=*
1079*=            the total spectral irradiance at each altitude level           =*
1080!-----------------------------------------------------------------------------*
1081
1082      implicit none
1083
1084! input:
1085
1086      INTEGER nlev
1087      REAL zen, rsfc
1088      REAL tauu(nlev), omu(nlev), gu(nlev)
1089      REAL dsdh(0:nlev,nlev)
1090      INTEGER nid(0:nlev)
1091      LOGICAL delta
1092
1093! output:
1094
1095      REAL fup(nlev),fdn(nlev),fdr(nlev)
1096      REAL eup(nlev),edn(nlev),edr(nlev)
1097
1098! local:
1099
1100      REAL tausla(0:nlev), tauc(0:nlev)
1101      REAL mu2(0:nlev), mu, sum
1102
1103! internal coefficients and matrix
1104
1105      REAL lam(nlev),taun(nlev),bgam(nlev)
1106      REAL e1(nlev),e2(nlev),e3(nlev),e4(nlev)
1107      REAL cup(nlev),cdn(nlev),cuptn(nlev),cdntn(nlev)
1108      REAL mu1(nlev)
1109      INTEGER row
1110      REAL a(2*nlev),b(2*nlev),d(2*nlev),e(2*nlev),y(2*nlev)
1111
1112! other:
1113
1114      REAL pifs, fdn0
1115      REAL gi(nlev), omi(nlev), tempg
1116      REAL f, g, om
1117      REAL gam1, gam2, gam3, gam4
1118      real, parameter :: largest = 1.e+36
1119      real, parameter :: precis = 1.e-7
1120
1121! For calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1122! in delta-function, modified quadrature, hemispheric constant,
1123! Hybrid modified Eddington-delta function metods, p633,Table1.
1124! W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1125! W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440,
1126! uncomment the following two lines and the appropriate statements further
1127! down.
1128!     REAL YLM0, YLM2, YLM4, YLM6, YLM8, YLM10, YLM12, YLMS, BETA0,
1129!    >     BETA1, BETAn, amu1, subd
1130
1131      REAL expon, expon0, expon1, divisr, temp, up, dn
1132      REAL ssfc
1133      INTEGER nlayer, mrows, lev
1134
1135      INTEGER i, j
1136
1137! Some additional program constants:
1138
1139      real pi, dr
1140      REAL eps
1141      PARAMETER (eps = 1.E-3)
1142!_______________________________________________________________________
1143
1144! MU = cosine of solar zenith angle
1145! RSFC = surface albedo
1146! TAUU =  unscaled optical depth of each layer
1147! OMU  =  unscaled single scattering albedo
1148! GU   =  unscaled asymmetry factor
1149! KLEV = max dimension of number of layers in atmosphere
1150! NLAYER = number of layers in the atmosphere
1151! NLEVEL = nlayer + 1 = number of levels
1152
1153! initial conditions:  pi*solar flux = 1;  diffuse incidence = 0
1154
1155      pifs = 1.     
1156      fdn0 = 0.
1157
1158      nlayer = nlev - 1
1159
1160      pi = acos(-1.)
1161      dr = pi/180.
1162      mu = COS(zen*dr)
1163
1164!************* compute coefficients for each layer:
1165! GAM1 - GAM4 = 2-stream coefficients, different for different approximations
1166! EXPON0 = calculation of e when TAU is zero
1167! EXPON1 = calculation of e when TAU is TAUN
1168! CUP and CDN = calculation when TAU is zero
1169! CUPTN and CDNTN = calc. when TAU is TAUN
1170! DIVISR = prevents division by zero
1171
1172      do j = 0, nlev
1173         tauc(j) = 0.
1174         tausla(j) = 0.
1175         mu2(j) = 1./SQRT(largest)
1176      end do
1177
1178      IF (.NOT. delta) THEN
1179         DO i = 1, nlayer
1180            gi(i) = gu(i)
1181            omi(i) = omu(i)
1182            taun(i) = tauu(i)
1183         END DO
1184      ELSE
1185
1186! delta-scaling. Have to be done for delta-Eddington approximation,
1187! delta discrete ordinate, Practical Improved Flux Method, delta function,
1188! and Hybrid modified Eddington-delta function methods approximations
1189
1190         DO i = 1, nlayer
1191            f = gu(i)*gu(i)
1192            gi(i) = (gu(i) - f)/(1 - f)
1193            omi(i) = (1 - f)*omu(i)/(1 - omu(i)*f)       
1194            taun(i) = (1 - omu(i)*f)*tauu(i)
1195         END DO
1196      END IF
1197
1198! calculate slant optical depth at the top of the atmosphere when zen>90.
1199! in this case, higher altitude of the top layer is recommended.
1200
1201      IF (zen .GT. 90.0) THEN
1202         IF (nid(0) .LT. 0) THEN
1203            tausla(0) = largest
1204         ELSE
1205            sum = 0.0
1206            DO j = 1, nid(0)
1207               sum = sum + 2.*taun(j)*dsdh(0,j)
1208            END DO
1209            tausla(0) = sum
1210         END IF
1211      END IF
1212
1213      DO 11, i = 1, nlayer
1214          g = gi(i)
1215          om = omi(i)
1216          tauc(i) = tauc(i-1) + taun(i)
1217
1218! stay away from 1 by precision.  For g, also stay away from -1
1219
1220          tempg = AMIN1(abs(g),1. - precis)
1221          g = SIGN(tempg,g)
1222          om = AMIN1(om,1.-precis)
1223
1224! calculate slant optical depth
1225             
1226          IF (nid(i) .LT. 0) THEN
1227             tausla(i) = largest
1228          ELSE
1229             sum = 0.0
1230             DO j = 1, MIN(nid(i),i)
1231                sum = sum + taun(j)*dsdh(i,j)
1232             END DO
1233             DO j = MIN(nid(i),i)+1,nid(i)
1234                sum = sum + 2.*taun(j)*dsdh(i,j)
1235             END DO
1236             tausla(i) = sum
1237             IF (tausla(i) .EQ. tausla(i-1)) THEN
1238                mu2(i) = SQRT(largest)
1239             ELSE
1240                mu2(i) = (tauc(i)-tauc(i-1))/(tausla(i)-tausla(i-1))
1241                mu2(i) = SIGN( AMAX1(ABS(mu2(i)),1./SQRT(largest)),
1242     $                     mu2(i) )
1243             END IF
1244          END IF
1245
1246!** the following gamma equations are from pg 16,289, Table 1
1247!** save mu1 for each approx. for use in converting irradiance to actinic flux
1248
1249! Eddington approximation(Joseph et al., 1976, JAS, 33, 2452):
1250
1251c       gam1 =  (7. - om*(4. + 3.*g))/4.
1252c       gam2 = -(1. - om*(4. - 3.*g))/4.
1253c       gam3 = (2. - 3.*g*mu)/4.
1254c       gam4 = 1. - gam3
1255c       mu1(i) = 0.5
1256
1257* quadrature (Liou, 1973, JAS, 30, 1303-1326; 1974, JAS, 31, 1473-1475):
1258
1259c          gam1 = 1.7320508*(2. - om*(1. + g))/2.
1260c          gam2 = 1.7320508*om*(1. - g)/2.
1261c          gam3 = (1. - 1.7320508*g*mu)/2.
1262c          gam4 = 1. - gam3
1263c          mu1(i) = 1./sqrt(3.)
1264         
1265* hemispheric mean (Toon et al., 1089, JGR, 94, 16287):
1266
1267           gam1 = 2. - om*(1. + g)
1268           gam2 = om*(1. - g)
1269           gam3 = (2. - g*mu)/4.
1270           gam4 = 1. - gam3
1271           mu1(i) = 0.5
1272
1273* PIFM  (Zdunkovski et al.,1980, Conrib.Atmos.Phys., 53, 147-166):
1274c         GAM1 = 0.25*(8. - OM*(5. + 3.*G))
1275c         GAM2 = 0.75*OM*(1.-G)
1276c         GAM3 = 0.25*(2.-3.*G*MU)
1277c         GAM4 = 1. - GAM3
1278c         mu1(i) = 0.5
1279
1280* delta discrete ordinates  (Schaller, 1979, Contrib.Atmos.Phys, 52, 17-26):
1281c         GAM1 = 0.5*1.7320508*(2. - OM*(1. + G))
1282c         GAM2 = 0.5*1.7320508*OM*(1.-G)
1283c         GAM3 = 0.5*(1.-1.7320508*G*MU)
1284c         GAM4 = 1. - GAM3
1285c         mu1(i) = 1./sqrt(3.)
1286
1287* Calculations of Associated Legendre Polynomials for GAMA1,2,3,4
1288* in delta-function, modified quadrature, hemispheric constant,
1289* Hybrid modified Eddington-delta function metods, p633,Table1.
1290* W.E.Meador and W.R.Weaver, GAS,1980,v37,p.630
1291* W.J.Wiscombe and G.W. Grams, GAS,1976,v33,p2440
1292c      YLM0 = 2.
1293c      YLM2 = -3.*G*MU
1294c      YLM4 = 0.875*G**3*MU*(5.*MU**2-3.)
1295c      YLM6=-0.171875*G**5*MU*(15.-70.*MU**2+63.*MU**4)
1296c     YLM8=+0.073242*G**7*MU*(-35.+315.*MU**2-693.*MU**4
1297c    *+429.*MU**6)
1298c     YLM10=-0.008118*G**9*MU*(315.-4620.*MU**2+18018.*MU**4
1299c    *-25740.*MU**6+12155.*MU**8)
1300c     YLM12=0.003685*G**11*MU*(-693.+15015.*MU**2-90090.*MU**4
1301c    *+218790.*MU**6-230945.*MU**8+88179.*MU**10)
1302c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1303c      YLMS=0.25*YLMS
1304c      BETA0 = YLMS
1305c
1306c         amu1=1./1.7320508
1307c      YLM0 = 2.
1308c      YLM2 = -3.*G*amu1
1309c      YLM4 = 0.875*G**3*amu1*(5.*amu1**2-3.)
1310c      YLM6=-0.171875*G**5*amu1*(15.-70.*amu1**2+63.*amu1**4)
1311c     YLM8=+0.073242*G**7*amu1*(-35.+315.*amu1**2-693.*amu1**4
1312c    *+429.*amu1**6)
1313c     YLM10=-0.008118*G**9*amu1*(315.-4620.*amu1**2+18018.*amu1**4
1314c    *-25740.*amu1**6+12155.*amu1**8)
1315c     YLM12=0.003685*G**11*amu1*(-693.+15015.*amu1**2-90090.*amu1**4
1316c    *+218790.*amu1**6-230945.*amu1**8+88179.*amu1**10)
1317c      YLMS=YLM0+YLM2+YLM4+YLM6+YLM8+YLM10+YLM12
1318c      YLMS=0.25*YLMS
1319c      BETA1 = YLMS
1320c
1321c         BETAn = 0.25*(2. - 1.5*G-0.21875*G**3-0.085938*G**5
1322c    *-0.045776*G**7)
1323
1324
1325* Hybrid modified Eddington-delta function(Meador and Weaver,1980,JAS,37,630):
1326c         subd=4.*(1.-G*G*(1.-MU))
1327c         GAM1 = (7.-3.*G*G-OM*(4.+3.*G)+OM*G*G*(4.*BETA0+3.*G))/subd
1328c         GAM2 =-(1.-G*G-OM*(4.-3.*G)-OM*G*G*(4.*BETA0+3.*G-4.))/subd
1329c         GAM3 = BETA0
1330c         GAM4 = 1. - GAM3
1331c         mu1(i) = (1. - g*g*(1.- mu) )/(2. - g*g)
1332
1333*****
1334* delta function  (Meador, and Weaver, 1980, JAS, 37, 630):
1335c         GAM1 = (1. - OM*(1. - beta0))/MU
1336c         GAM2 = OM*BETA0/MU
1337c         GAM3 = BETA0
1338c         GAM4 = 1. - GAM3
1339c         mu1(i) = mu
1340*****
1341* modified quadrature (Meador, and Weaver, 1980, JAS, 37, 630):
1342c         GAM1 = 1.7320508*(1. - OM*(1. - beta1))
1343c         GAM2 = 1.7320508*OM*beta1
1344c         GAM3 = BETA0
1345c         GAM4 = 1. - GAM3
1346c         mu1(i) = 1./sqrt(3.)
1347
1348* hemispheric constant (Toon et al., 1989, JGR, 94, 16287):
1349c         GAM1 = 2.*(1. - OM*(1. - betan))
1350c         GAM2 = 2.*OM*BETAn
1351c         GAM3 = BETA0
1352c         GAM4 = 1. - GAM3
1353c         mu1(i) = 0.5
1354
1355*****
1356
1357* lambda = pg 16,290 equation 21
1358* big gamma = pg 16,290 equation 22
1359* if gam2 = 0., then bgam = 0.
1360
1361         lam(i) = sqrt(gam1*gam1 - gam2*gam2)
1362
1363         IF (gam2 .NE. 0.) THEN
1364            bgam(i) = (gam1 - lam(i))/gam2
1365         ELSE
1366            bgam(i) = 0.
1367         END IF
1368
1369         expon = EXP(-lam(i)*taun(i))
1370
1371* e1 - e4 = pg 16,292 equation 44
1372         
1373         e1(i) = 1. + bgam(i)*expon
1374         e2(i) = 1. - bgam(i)*expon
1375         e3(i) = bgam(i) + expon
1376         e4(i) = bgam(i) - expon
1377
1378* the following sets up for the C equations 23, and 24
1379* found on page 16,290
1380* prevent division by zero (if LAMBDA=1/MU, shift 1/MU^2 by EPS = 1.E-3
1381* which is approx equiv to shifting MU by 0.5*EPS* (MU)**3
1382
1383         expon0 = EXP(-tausla(i-1))
1384         expon1 = EXP(-tausla(i))
1385         
1386         divisr = lam(i)*lam(i) - 1./(mu2(i)*mu2(i))
1387         temp = AMAX1(eps,abs(divisr))
1388         divisr = SIGN(temp,divisr)
1389
1390         up = om*pifs*((gam1 - 1./mu2(i))*gam3 + gam4*gam2)/divisr
1391         dn = om*pifs*((gam1 + 1./mu2(i))*gam4 + gam2*gam3)/divisr
1392         
1393* cup and cdn are when tau is equal to zero
1394* cuptn and cdntn are when tau is equal to taun
1395
1396         cup(i) = up*expon0
1397         cdn(i) = dn*expon0
1398         cuptn(i) = up*expon1
1399         cdntn(i) = dn*expon1
1400 
1401   11 CONTINUE
1402
1403***************** set up matrix ******
1404* ssfc = pg 16,292 equation 37  where pi Fs is one (unity).
1405
1406      ssfc = rsfc*mu*EXP(-tausla(nlayer))*pifs
1407
1408* MROWS = the number of rows in the matrix
1409
1410      mrows = 2*nlayer     
1411     
1412* the following are from pg 16,292  equations 39 - 43.
1413* set up first row of matrix:
1414
1415      i = 1
1416      a(1) = 0.
1417      b(1) = e1(i)
1418      d(1) = -e2(i)
1419      e(1) = fdn0 - cdn(i)
1420
1421      row=1
1422
1423* set up odd rows 3 thru (MROWS - 1):
1424
1425      i = 0
1426      DO 20, row = 3, mrows - 1, 2
1427         i = i + 1
1428         a(row) = e2(i)*e3(i) - e4(i)*e1(i)
1429         b(row) = e1(i)*e1(i + 1) - e3(i)*e3(i + 1)
1430         d(row) = e3(i)*e4(i + 1) - e1(i)*e2(i + 1)
1431         e(row) = e3(i)*(cup(i + 1) - cuptn(i)) +
1432     $        e1(i)*(cdntn(i) - cdn(i + 1))
1433   20 CONTINUE
1434
1435* set up even rows 2 thru (MROWS - 2):
1436
1437      i = 0
1438      DO 30, row = 2, mrows - 2, 2
1439         i = i + 1
1440         a(row) = e2(i + 1)*e1(i) - e3(i)*e4(i + 1)
1441         b(row) = e2(i)*e2(i + 1) - e4(i)*e4(i + 1)
1442         d(row) = e1(i + 1)*e4(i + 1) - e2(i + 1)*e3(i + 1)
1443         e(row) = (cup(i + 1) - cuptn(i))*e2(i + 1) -
1444     $        (cdn(i + 1) - cdntn(i))*e4(i + 1)
1445   30 CONTINUE
1446
1447* set up last row of matrix at MROWS:
1448
1449      row = mrows
1450      i = nlayer
1451     
1452      a(row) = e1(i) - rsfc*e3(i)
1453      b(row) = e2(i) - rsfc*e4(i)
1454      d(row) = 0.
1455      e(row) = ssfc - cuptn(i) + rsfc*cdntn(i)
1456
1457* solve tri-diagonal matrix:
1458
1459      CALL tridiag(a, b, d, e, y, mrows)
1460
1461**** unfold solution of matrix, compute output fluxes:
1462
1463      row = 1
1464      lev = 1
1465      j = 1
1466     
1467* the following equations are from pg 16,291  equations 31 & 32
1468
1469      fdr(lev) = EXP( -tausla(0) )
1470      edr(lev) = mu * fdr(lev)
1471      edn(lev) = fdn0
1472      eup(lev) =  y(row)*e3(j) - y(row + 1)*e4(j) + cup(j)
1473      fdn(lev) = edn(lev)/mu1(lev)
1474      fup(lev) = eup(lev)/mu1(lev)
1475
1476      DO 60, lev = 2, nlayer + 1
1477         fdr(lev) = EXP(-tausla(lev-1))
1478         edr(lev) =  mu *fdr(lev)
1479         edn(lev) =  y(row)*e3(j) + y(row + 1)*e4(j) + cdntn(j)
1480         eup(lev) =  y(row)*e1(j) + y(row + 1)*e2(j) + cuptn(j)
1481         fdn(lev) = edn(lev)/mu1(j)
1482         fup(lev) = eup(lev)/mu1(j)
1483
1484         row = row + 2
1485         j = j + 1
1486   60 CONTINUE
1487
1488      end subroutine ps2str
1489
1490*=============================================================================*
1491
1492      subroutine tridiag(a,b,c,r,u,n)
1493
1494!_______________________________________________________________________
1495! solves tridiagonal system.  From Numerical Recipies, p. 40
1496!_______________________________________________________________________
1497
1498      IMPLICIT NONE
1499
1500! input:
1501
1502      INTEGER n
1503      REAL a, b, c, r
1504      DIMENSION a(n),b(n),c(n),r(n)
1505
1506! output:
1507
1508      REAL u
1509      DIMENSION u(n)
1510
1511! local:
1512
1513      INTEGER j
1514
1515      REAL bet, gam
1516      DIMENSION gam(n)
1517!_______________________________________________________________________
1518
1519      IF (b(1) .EQ. 0.) STOP 1001
1520      bet   = b(1)
1521      u(1) = r(1)/bet
1522      DO 11, j = 2, n   
1523         gam(j) = c(j - 1)/bet
1524         bet = b(j) - a(j)*gam(j)
1525         IF (bet .EQ. 0.) STOP 2002
1526         u(j) = (r(j) - a(j)*u(j - 1))/bet
1527   11 CONTINUE
1528      DO 12, j = n - 1, 1, -1 
1529         u(j) = u(j) - gam(j + 1)*u(j + 1)
1530   12 CONTINUE
1531!_______________________________________________________________________
1532
1533      end subroutine tridiag
1534
1535      end subroutine photolysis_online
Note: See TracBrowser for help on using the repository browser.