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

Last change on this file since 1621 was 1430, checked in by flefevre, 10 years ago

aeronomars : preparation de l'implementation du solveur chimique ASIS

1) extension de la table de photodissociation "jmars" a la thermosphere

et restriction aux seules especes necessaires.

2) modernisation de la subroutine d'interpolation dans cette

table de photodissociation.

3) chimiedata.h compatible avec solveur ASIS

ATTENTION!! la nouvelle table de photodissociation

jmars.20140930

doit etre presente dans /datadir sinon crash modele...
File size: 12.6 KB
Line 
1!==========================================================================
2
3      subroutine photolysis(nlayer,                              &
4                            lswitch, press, temp, sza, tauref,   &
5                            dist_sol, rmco2, rmo3, j)
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 :: dist_sol              ! sun distance (AU)
26      real :: rmco2(nlayer)         ! co2 mixing ratio
27      real :: rmo3(nlayer)          ! ozone mixing ratio
28
29!==========================================================================
30!     output: interpolated photodissociation rates (s-1)
31!==========================================================================
32
33      real :: j(nlayer,nd)          ! interpolated photolysis rates (s-1)
34
35!==========================================================================
36!     local:
37!==========================================================================
38
39      integer :: icol, ij, indsza, indtau, indcol, indozo, indtemp,     &
40                 iozo, isza, itau, it, l
41
42      real :: col(nlayer)                 ! overhead air column   (molecule cm-2)
43      real :: colo3(nlayer)               ! overhead ozone column (molecule cm-2)
44      real :: poids(2,2,2,2,2)            ! 5D interpolation weights
45      real :: tref                        ! temperature  at 1.9 hPa in the gcm (K)
46      real :: table_temp(ntemp)           ! temperatures at 1.9 hPa in jmars   (K)
47      real :: cinf, csup, cicol, ciozo, cisza, citemp, citau
48      real :: colo3min, dp, coef
49      real :: ratio_o3(nlayer)
50      real :: tau
51
52!==========================================================================
53!     day/night criterion
54!==========================================================================
55
56      if (sza <= 95.) then
57
58!==========================================================================
59!     temperatures at 1.9 hPa in lookup table
60!==========================================================================
61     
62      table_temp(1) = 226.2
63      table_temp(2) = 206.2
64      table_temp(3) = 186.2
65      table_temp(4) = 169.8
66
67!==========================================================================
68!     interpolation in solar zenith angle
69!==========================================================================
70 
71      indsza = nsza - 1
72      do isza = 1,nsza
73         if (szatab(isza) >= sza) then
74            indsza = isza - 1
75            indsza = max(indsza, 1)
76            exit
77         end if
78      end do
79      cisza = (sza - szatab(indsza))  &
80             /(szatab(indsza + 1) - szatab(indsza))
81
82!==========================================================================
83!     interpolation in dust (tau)
84!==========================================================================
85
86      tau = min(tauref, tautab(ntau))
87      tau = max(tau, tautab(1))
88
89      indtau = ntau - 1
90      do itau = 1,ntau
91         if (tautab(itau) >= tau) then
92            indtau = itau - 1
93            indtau = max(indtau, 1)
94            exit
95         end if
96      end do
97      citau = (tau - tautab(indtau))     &
98             /(tautab(indtau + 1) - tautab(indtau))
99
100!==========================================================================
101!     co2 and ozone columns
102!==========================================================================
103
104!     co2 column at model top (molecule.cm-2)
105
106      col(lswitch-1) = 6.022e22*rmco2(lswitch-1)*press(lswitch-1)*100.  &
107                       /(mugaz*g)
108
109!     ozone column at model top
110
111      colo3(lswitch-1) = 0.
112
113!     co2 and ozone columns for other levels (molecule.cm-2)
114
115      do l = lswitch-2,1,-1
116         dp = (press(l) - press(l+1))*100.
117         col(l) = col(l+1) + (rmco2(l+1) + rmco2(l))*0.5   &
118                             *6.022e22*dp/(mugaz*g)
119         col(l) = min(col(l), colairtab(1))
120         colo3(l) = colo3(l+1) + (rmo3(l+1) + rmo3(l))*0.5 &
121                                 *6.022e22*dp/(mugaz*g)
122      end do
123
124!     ratio ozone column/minimal theoretical column (0.1 micron-atm)
125
126!     ro3 = 7.171e-10 is the o3 mixing ratio for a uniform
127!     profile giving a column 0.1 micron-atmosphere at
128!     a surface pressure of 10 hpa.
129
130      do l = 1,lswitch-1
131         colo3min    = col(l)*7.171e-10
132         ratio_o3(l) = colo3(l)/colo3min
133         ratio_o3(l) = min(ratio_o3(l), table_ozo(nozo)*10.)
134         ratio_o3(l) = max(ratio_o3(l), 1.)
135      end do
136
137!==========================================================================
138!     temperature dependence
139!==========================================================================
140
141!     1) search for temperature at 1.9 hPa (tref): vertical interpolation
142
143      tref = temp(1)
144      do l = (lswitch-1)-1,1,-1
145         if (press(l) > 1.9) then
146            cinf = (press(l) - 1.9)/(press(l) - press(l+1))
147            csup = 1. - cinf
148            tref = cinf*temp(l+1) + csup*temp(l)
149            exit
150         end if
151      end do
152
153!     2) interpolation in temperature
154
155      tref = min(tref,table_temp(1))
156      tref = max(tref,table_temp(ntemp))
157
158      do it = 2, ntemp
159         if (table_temp(it) <= tref) then
160            citemp = (log(tref) - log(table_temp(it)))              &
161                    /(log(table_temp(it-1)) - log(table_temp(it)))
162            indtemp = it - 1
163            exit
164         end if
165      end do
166
167!==========================================================================
168!     loop over vertical levels
169!==========================================================================
170
171      do l = 1,lswitch-1
172
173!     interpolation in air column
174
175         indcol = nz - 1
176         do icol = 1,nz
177            if (colairtab(icol) < col(l)) then
178               indcol = icol - 1
179               exit
180            end if
181         end do
182         cicol = (log(col(l)) - log(colairtab(indcol + 1)))              &
183                /(log(colairtab(indcol)) - log(colairtab(indcol + 1)))
184
185!     interpolation in ozone column
186
187         indozo = nozo - 1
188         do iozo = 1,nozo
189            if (table_ozo(iozo)*10. >= ratio_o3(l)) then
190               indozo = iozo - 1
191               indozo = max(indozo, 1)
192               exit
193            end if
194         end do
195         ciozo = (ratio_o3(l) - table_ozo(indozo)*10.)             &
196                /(table_ozo(indozo + 1)*10. - table_ozo(indozo)*10.)
197
198!     4-dimensional interpolation weights
199
200!     poids(temp,sza,co2,o3,tau)
201
202         poids(1,1,1,1,1) = citemp*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
203         poids(1,1,1,2,1) = citemp*(1.-cisza)*cicol*ciozo*(1.-citau)
204         poids(1,1,2,1,1) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
205         poids(1,1,2,2,1) = citemp*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
206         poids(1,2,1,1,1) = citemp*cisza*cicol*(1.-ciozo)*(1.-citau)
207         poids(1,2,1,2,1) = citemp*cisza*cicol*ciozo*(1.-citau)
208         poids(1,2,2,1,1) = citemp*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
209         poids(1,2,2,2,1) = citemp*cisza*(1.-cicol)*ciozo*(1.-citau)
210         poids(2,1,1,1,1) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*(1.-citau)
211         poids(2,1,1,2,1) = (1.-citemp)*(1.-cisza)*cicol*ciozo*(1.-citau)
212         poids(2,1,2,1,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*(1.-citau)
213         poids(2,1,2,2,1) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*(1.-citau)
214         poids(2,2,1,1,1) = (1.-citemp)*cisza*cicol*(1.-ciozo)*(1.-citau)
215         poids(2,2,1,2,1) = (1.-citemp)*cisza*cicol*ciozo*(1.-citau)
216         poids(2,2,2,1,1) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*(1.-citau)
217         poids(2,2,2,2,1) = (1.-citemp)*cisza*(1.-cicol)*ciozo*(1.-citau)
218!
219         poids(1,1,1,1,2) = citemp*(1.-cisza)*cicol*(1.-ciozo)*citau
220         poids(1,1,1,2,2) = citemp*(1.-cisza)*cicol*ciozo*citau
221         poids(1,1,2,1,2) = citemp*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
222         poids(1,1,2,2,2) = citemp*(1.-cisza)*(1.-cicol)*ciozo*citau
223         poids(1,2,1,1,2) = citemp*cisza*cicol*(1.-ciozo)*citau
224         poids(1,2,1,2,2) = citemp*cisza*cicol*ciozo*citau
225         poids(1,2,2,1,2) = citemp*cisza*(1.-cicol)*(1.-ciozo)*citau
226         poids(1,2,2,2,2) = citemp*cisza*(1.-cicol)*ciozo*citau
227         poids(2,1,1,1,2) = (1.-citemp)*(1.-cisza)*cicol*(1.-ciozo)*citau
228         poids(2,1,1,2,2) = (1.-citemp)*(1.-cisza)*cicol*ciozo*citau
229         poids(2,1,2,1,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*(1.-ciozo)*citau
230         poids(2,1,2,2,2) = (1.-citemp)*(1.-cisza)*(1.-cicol)*ciozo*citau
231         poids(2,2,1,1,2) = (1.-citemp)*cisza*cicol*(1.-ciozo)*citau
232         poids(2,2,1,2,2) = (1.-citemp)*cisza*cicol*ciozo*citau
233         poids(2,2,2,1,2) = (1.-citemp)*cisza*(1.-cicol)*(1.-ciozo)*citau
234         poids(2,2,2,2,2) = (1.-citemp)*cisza*(1.-cicol)*ciozo*citau
235
236!     4-dimensional interpolation in the lookup table
237
238         do ij = 1,nd
239            j(l,ij) =                                                                &
240            poids(1,1,1,1,1)*jphot(indtemp,indsza,indcol,indozo,indtau,ij)           &
241          + poids(1,1,1,2,1)*jphot(indtemp,indsza,indcol,indozo+1,indtau,ij)         &
242          + poids(1,1,2,1,1)*jphot(indtemp,indsza,indcol+1,indozo,indtau,ij)         &
243          + poids(1,1,2,2,1)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau,ij)       &
244          + poids(1,2,1,1,1)*jphot(indtemp,indsza+1,indcol,indozo,indtau,ij)         &
245          + poids(1,2,1,2,1)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau,ij)       &
246          + poids(1,2,2,1,1)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau,ij)       &
247          + poids(1,2,2,2,1)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau,ij)     &
248          + poids(2,1,1,1,1)*jphot(indtemp+1,indsza,indcol,indozo,indtau,ij)         &
249          + poids(2,1,1,2,1)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau,ij)       &
250          + poids(2,1,2,1,1)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau,ij)       &
251          + poids(2,1,2,2,1)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau,ij)     &
252          + poids(2,2,1,1,1)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau,ij)       &
253          + poids(2,2,1,2,1)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau,ij)     &
254          + poids(2,2,2,1,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau,ij)     &
255          + poids(2,2,2,2,1)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau,ij)   &
256!
257          + poids(1,1,1,1,2)*jphot(indtemp,indsza,indcol,indozo,indtau+1,ij)         &
258          + poids(1,1,1,2,2)*jphot(indtemp,indsza,indcol,indozo+1,indtau+1,ij)       &
259          + poids(1,1,2,1,2)*jphot(indtemp,indsza,indcol+1,indozo,indtau+1,ij)       &
260          + poids(1,1,2,2,2)*jphot(indtemp,indsza,indcol+1,indozo+1,indtau+1,ij)     &
261          + poids(1,2,1,1,2)*jphot(indtemp,indsza+1,indcol,indozo,indtau+1,ij)       &
262          + poids(1,2,1,2,2)*jphot(indtemp,indsza+1,indcol,indozo+1,indtau+1,ij)     &
263          + poids(1,2,2,1,2)*jphot(indtemp,indsza+1,indcol+1,indozo,indtau+1,ij)     &
264          + poids(1,2,2,2,2)*jphot(indtemp,indsza+1,indcol+1,indozo+1,indtau+1,ij)   &
265          + poids(2,1,1,1,2)*jphot(indtemp+1,indsza,indcol,indozo,indtau+1,ij)       &
266          + poids(2,1,1,2,2)*jphot(indtemp+1,indsza,indcol,indozo+1,indtau+1,ij)     &
267          + poids(2,1,2,1,2)*jphot(indtemp+1,indsza,indcol+1,indozo,indtau+1,ij)     &
268          + poids(2,1,2,2,2)*jphot(indtemp+1,indsza,indcol+1,indozo+1,indtau+1,ij)   &
269          + poids(2,2,1,1,2)*jphot(indtemp+1,indsza+1,indcol,indozo,indtau+1,ij)     &
270          + poids(2,2,1,2,2)*jphot(indtemp+1,indsza+1,indcol,indozo+1,indtau+1,ij)   &
271          + poids(2,2,2,1,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo,indtau+1,ij)   &
272          + poids(2,2,2,2,2)*jphot(indtemp+1,indsza+1,indcol+1,indozo+1,indtau+1,ij)
273         end do
274
275!     correction for sun distance
276
277         do ij = 1,nd
278            j(l,ij) = j(l,ij)*(1.52/dist_sol)**2.
279         end do
280
281!==========================================================================
282!     end of loop over vertical levels
283!==========================================================================
284
285      end do
286
287      else
288
289!==========================================================================
290!     night
291!==========================================================================
292
293         j(:,:) = 0.
294
295      end if
296
297      return
298      end
Note: See TracBrowser for help on using the repository browser.