source: trunk/LMDZ.MARS/libf/aeronomars/photolysis.F90 @ 3026

Last change on this file since 3026 was 3012, checked in by emillour, 17 months ago

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