source: trunk/LMDZ.MARS/libf/aeronomars/photolysis_asis.F90 @ 1495

Last change on this file since 1495 was 1495, checked in by flefevre, 9 years ago

Implementation du coeur chimique Adaptative Semi-Implicit Symmetric (ASIS)
Par defaut non actif pour l'instant (variable logique asis = .false.)

File size: 16.1 KB
Line 
1!==========================================================================
2
3      subroutine photolysis_asis(nlayer,                              &
4                                 lswitch, press, temp, sza, tauref,   &
5                                 zmmean, dist_sol, rmco2, rmo3, v_phot)
6
7!==========================================================================
8
9      use comcstfi_h
10
11      implicit none
12
13#include "chimiedata.h"
14
15!==========================================================================
16!     input:
17!==========================================================================
18       
19      integer, intent(in) :: nlayer ! number of atmospheric layers
20      integer :: lswitch            ! interface level between chemistries
21      real :: press(nlayer)         ! pressure (hPa)
22      real :: temp(nlayer)          ! temperature (K)
23      real :: sza                   ! solar zenith angle (deg)
24      real :: tauref                ! optical depth at 7 hpa
25      real :: zmmean(nlayer)        ! mean molecular mass (g)
26      real :: dist_sol              ! sun distance (AU)
27      real :: rmco2(nlayer)         ! co2 mixing ratio
28      real :: rmo3(nlayer)          ! ozone mixing ratio
29
30!==========================================================================
31!     output: interpolated photodissociation rates (s-1)
32!==========================================================================
33
34      real (kind = 8), dimension(nlayer,nb_phot_max) :: v_phot
35
36!==========================================================================
37!     local:
38!==========================================================================
39
40      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
41                 iozo, isza, itau, it, l
42
43      integer :: j_o2_o, j_o2_o1d, j_co2_o, j_co2_o1d, j_o3_o1d,        &
44                 j_o3_o, j_h2o, j_hdo, j_h2o2, j_ho2, j_no, j_no2,      &
45                 j_hno3, j_hno4,                                        &
46                 j_ch4_ch3_h, j_ch4_1ch2_h2, j_ch4_3ch2_h_h,            &
47                 j_ch4_ch_h2_h, j_ch3o2h, j_ch2o_hco, j_ch2o_co,        &
48                 j_ch3oh, j_c2h6, j_hcl, j_hocl, j_clo, j_so2, j_so,    &
49                 j_h2s, j_so3
50
51      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
52      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
53      real :: poids(2,2,2,2,2)            ! 5D interpolation weights
54      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
55      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
56      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau
57      real :: colo3min, dp, coef
58      real :: ratio_o3(nlayer)
59      real :: tau
60      real :: j(nlayer,nd)
61
62!==========================================================================
63!     day/night criterion
64!==========================================================================
65
66      if (sza <= 95.) then
67
68!==========================================================================
69!     temperatures at 1.9 hPa in lookup table
70!==========================================================================
71     
72      table_temp(1) = 226.2
73      table_temp(2) = 206.2
74      table_temp(3) = 186.2
75      table_temp(4) = 169.8
76
77!==========================================================================
78!     interpolation in solar zenith angle
79!==========================================================================
80 
81      indsza = nsza - 1
82      do isza = 1,nsza
83         if (szatab(isza) >= sza) then
84            indsza = isza - 1
85            indsza = max(indsza, 1)
86            exit
87         end if
88      end do
89      cisza = (sza - szatab(indsza))  &
90             /(szatab(indsza + 1) - szatab(indsza))
91
92!==========================================================================
93!     interpolation in dust (tau)
94!==========================================================================
95
96      tau = min(tauref, tautab(ntau))
97      tau = max(tau, tautab(1))
98
99      indtau = ntau - 1
100      do itau = 1,ntau
101         if (tautab(itau) >= tau) then
102            indtau = itau - 1
103            indtau = max(indtau, 1)
104            exit
105         end if
106      end do
107      citau = (tau - tautab(indtau))     &
108             /(tautab(indtau + 1) - tautab(indtau))
109
110!==========================================================================
111!     co2 and ozone columns
112!==========================================================================
113
114!     co2 column at model top (molecule.cm-2)
115
116      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
117                       /(zmmean(lswitch-1)*g)
118
119!     ozone column at model top
120
121      colo3(lswitch-1) = 0.
122
123!     co2 and ozone columns for other levels (molecule.cm-2)
124
125      do l = lswitch-2,1,-1
126         dp = (press(l) - press(l+1))*100.
127         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
128                             *6.022e22*dp/(zmmean(l)*g)
129         col(l) = min(col(l), colairtab(1))
130         colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
131                                 *6.022e22*dp/(zmmean(l)*g)
132      end do
133
134!     ratio ozone column/minimal theoretical column (0.1 micron-atm)
135
136!     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
137!     profile giving a column 0.1 micron-atmosphere at
138!     a surface pressure of 10 hpa.
139
140      do l = 1,lswitch-1
141         colo3min    = col(l)*7.171e-10
142         ratio_o3(l) = colo3(l)/colo3min
143         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
144         ratio_o3(l) = max(ratio_o3(l), 1.)
145      end do
146
147!==========================================================================
148!     temperature dependence
149!==========================================================================
150
151!     1) search for temperature at 1.9 hPa (tref): vertical interpolation
152
153      tref = temp(1)
154      do l = (lswitch-1)-1,1,-1
155         if (press(l) > 1.9) then
156            cinf = (press(l) - 1.9)/(press(l) - press(l+1))
157            csup = 1. - cinf
158            tref = cinf*temp(l+1) + csup*temp(l)
159            exit
160         end if
161      end do
162
163!     2) interpolation in temperature
164
165      tref = min(tref,table_temp(1))
166      tref = max(tref,table_temp(ntemp))
167
168      do it = 2, ntemp
169         if (table_temp(it) <= tref) then
170            citemp = (log(tref) - log(table_temp(it)))              &
171                    /(log(table_temp(it-1)) - log(table_temp(it)))
172            indtemp = it - 1
173            exit
174         end if
175      end do
176
177!==========================================================================
178!     loop over vertical levels
179!==========================================================================
180
181      do l = 1,lswitch-1
182
183!     interpolation in air column
184
185         indcol = nz - 1
186         do icol = 1,nz
187            if (colairtab(icol) < col(l)) then
188               indcol = icol - 1
189               exit
190            end if
191         end do
192         cicol = (log(col(l)) - log(colairtab(indcol + 1)))              &
193                /(log(colairtab(indcol)) - log(colairtab(indcol + 1)))
194
195!     interpolation in ozone column
196
197         indozo = nozo - 1
198         do iozo = 1,nozo
199            if (table_ozo(iozo)*10. >= ratio_o3(l)) then
200               indozo = iozo - 1
201               indozo = max(indozo, 1)
202               exit
203            end if
204         end do
205         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)             &
206                /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
207
208!     4-dimensional interpolation weights
209
210!     poids(temp,sza,co2,o3,tau)
211
212         poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
213         poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)
214         poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
215         poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
216         poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)
217         poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau)
218         poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
219         poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)
220         poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
221         poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)
222         poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
223         poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
224         poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)
225         poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)
226         poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
227         poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)
228!
229         poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau
230         poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau
231         poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
232         poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau
233         poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau
234         poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau
235         poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau
236         poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau
237         poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau
238         poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau
239         poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
240         poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau
241         poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau
242         poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau
243         poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau
244         poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau
245
246!     4-dimensional interpolation in the lookup table
247
248         do ij = 1,nd
249            j(l,ij) =                                                                &
250            poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij)           &
251          + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)         &
252          + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)         &
253          + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)       &
254          + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)         &
255          + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)       &
256          + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)       &
257          + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)     &
258          + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)         &
259          + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)       &
260          + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)       &
261          + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)     &
262          + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)       &
263          + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)     &
264          + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)     &
265          + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)   &
266!
267          + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)         &
268          + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)       &
269          + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)       &
270          + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)     &
271          + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)       &
272          + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)     &
273          + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)     &
274          + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)   &
275          + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)       &
276          + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)     &
277          + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)     &
278          + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)   &
279          + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)     &
280          + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)   &
281          + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)   &
282          + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
283         end do
284
285!     correction for sun distance
286
287         do ij = 1,nd
288            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
289         end do
290
291!==========================================================================
292!     end of loop over vertical levels
293!==========================================================================
294
295      end do
296
297      else
298
299!==========================================================================
300!     night
301!==========================================================================
302
303         j(:,:) = 0.
304
305      end if
306
307! photodissociation rates numbering in the lookup table
308
309! jmars.20140930
310
311      j_o2_o         =  1      ! o2 + hv     -> o + o
312      j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
313      j_co2_o        =  3      ! co2 + hv    -> co + o
314      j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
315      j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
316      j_o3_o         =  6      ! o3 + hv     -> o2 + o
317      j_h2o          =  7      ! h2o + hv    -> h + oh
318      j_h2o2         =  8      ! h2o2 + hv   -> oh + oh
319      j_ho2          =  9      ! ho2 + hv    -> oh + o
320      j_no           =  10     ! no + hv     -> n + o
321      j_no2          =  11     ! no2 + hv    -> no + o
322      j_hno3         =  12     ! hno3 + hv   -> no2 + oh
323      j_hno4         =  13     ! hno4 + hv   -> no2 + ho2
324
325! jmars.20111014
326
327!     j_o2_o         =  1      ! o2 + hv     -> o + o
328!     j_o2_o1d       =  2      ! o2 + hv     -> o + o(1d)
329!     j_co2_o        =  3      ! co2 + hv    -> co + o
330!     j_co2_o1d      =  4      ! co2 + hv    -> co + o(1d)
331!     j_o3_o1d       =  5      ! o3 + hv     -> o2 + o(1d)
332!     j_o3_o         =  6      ! o3 + hv     -> o2 + o
333!     j_h2o          =  7      ! h2o + hv    -> h + oh
334!     j_hdo          =  8      ! hdo + hv    -> d + oh
335!     j_h2o2         =  9      ! h2o2 + hv   -> oh + oh
336!     j_ho2          =  10     ! ho2 + hv    -> oh + o
337!     j_no2          =  11     ! no2 + hv    -> no + o
338!     j_ch4_ch3_h    =  12     ! ch4 + hv    -> ch3 + h
339!     j_ch4_1ch2_h2  =  13     ! ch4 + hv    -> 1ch2 + h2
340!     j_ch4_3ch2_h_h =  14     ! ch4 + hv    -> 3ch2 + h + h
341!     j_ch4_ch_h2_h  =  15     ! ch4 + hv    -> ch + h2 + h
342!     j_ch3o2h       =  16     ! ch3o2h + hv -> ch3o + oh
343!     j_ch2o_hco     =  17     ! ch2o + hv   -> h + hco
344!     j_ch2o_co      =  18     ! ch2o + hv   -> h2 + co
345!     j_ch3oh        =  19     ! ch3oh + hv  -> ch3o + h
346!     j_c2h6         =  20     ! c2h6 + hv   -> products
347!     j_hcl          =  21     ! hcl + hv    -> h + cl
348!     j_hocl         =  22     ! hocl + hv   -> oh + cl
349!     j_clo          =  23     ! clo + hv    -> o + cl
350!     j_so2          =  24     ! so2 + hv    -> so + o
351!     j_so           =  25     ! so + hv     -> s + o
352!     j_h2s          =  26     ! h2s + hv    -> hs + s
353!     j_so3          =  27     ! so2 + hv    -> so2 + o
354!     j_hno3         =  28     ! hno3 + hv   -> oh + no2
355!     j_hno4         =  29     ! hno4 + hv   -> ho2 + no2
356
357! fill v_phot array
358
359      v_phot(:,:) = 0.
360
361      do l = 1,lswitch-1
362         v_phot(l, 1) = j(l,j_o2_o)
363         v_phot(l, 2) = j(l,j_o2_o1d)
364         v_phot(l, 3) = j(l,j_co2_o)
365         v_phot(l, 4) = j(l,j_co2_o1d)
366         v_phot(l, 5) = j(l,j_o3_o1d)
367         v_phot(l, 6) = j(l,j_o3_o)
368         v_phot(l, 7) = j(l,j_h2o)
369         v_phot(l, 8) = j(l,j_h2o2)
370         v_phot(l, 9) = j(l,j_ho2)
371         v_phot(l,10) = j(l,j_no)
372         v_phot(l,11) = j(l,j_no2)
373      end do
374
375      return
376      end
377
378!*****************************************************************
Note: See TracBrowser for help on using the repository browser.