source: trunk/LMDZ.MARS/libf/aeronomars/jthermcalc.F @ 226

Last change on this file since 226 was 38, checked in by emillour, 14 years ago

Ajout du modè Martien (mon LMDZ.MARS.BETA, du 28/01/2011) dans le rértoire mars, pour pouvoir suivre plus facilement les modifs.
EM

File size: 16.1 KB
Line 
1c**********************************************************************
2
3      subroutine jthermcalc
4     $ (co2x,o2x,o3px,h2x,h2ox,h2o2x,aux1,aux2,tx,nz,iz,date,zenit)
5
6
7c     feb 2002        fgg           first version
8c     nov 2002        fgg           second version
9c
10c modified from paramhr.F
11c MAC July 2003
12c**********************************************************************
13
14      implicit none
15
16c     common variables and constants
17
18      include 'param.h'
19      include 'param_v3.h'
20
21c    local parameters and variables
22
23      real       xabsi(nabs,nzmax) !densities
24      real       nada
25      real       co2colx2(nzmax,ninter)
26      real       o2colx2(nzmax,ninter)
27      real       o3pcolx2(nzmax,ninter)
28      real       co2colx(nzmax)              !column density of CO2 (cm^-2)
29      real       o2colx(nzmax)               !column density of O2(cm^-2)
30      real       o3pcolx(nzmax)              !column density of O(3P)(cm^-2)
31      real       h2o2colx(nzmax)             !column density of H2O2(cm^-2)
32      real       t2(nzmax)
33      real       coltemp(nzmax)
34      real       dt(nzmax)
35      real       sigma(ninter,nzmax)
36      real       alfa(ninter,nzmax)
37      real       xx
38      real       grav(nzmax)
39     
40      integer    i,j,k,icol,indexint          !indexes
41      character  dn
42
43c     input and output variables
44
45      integer    nz                          !number of layers
46      real       co2x(nz)                    !density of CO2(cm^-3)
47      real       o2x(nz)                     !density of O2(cm^-3)
48      real       o3px(nz)                    !density of O(3P)(cm^-3)
49      real       h2x(nz)                     !density of H2(cm^-3)
50      real       h2ox(nz)                    !density of H2O(cm^-3)
51      real       h2o2x(nz)                   !density of H2O2(cm^-3)
52      real       aux1(nz)                    !auxiliar variable
53      real       aux2(nz)                    !auxiliar variable
54      real       tx(nz)                      !temperature
55      real       date
56      real       zenit
57      real       iz(nz)
58
59
60c     variables used in interpolation
61
62      real aux3(nz2), aux4(nz2)
63                                              !auxiliar variables
64      real       limdown                      !limits for interpolation
65      real       limup                        !        ""
66
67c*************************PROGRAM STARTS*******************************
68
69      if(zenit.gt.100.) then
70         dn='n'
71         else
72         dn='d'
73      end if
74      if(dn.eq.'n') then
75        return
76      endif
77     
78      do i=1,nz
79         xabsi(1,i) = co2x(i)
80         xabsi(2,i) = o2x(i)
81         xabsi(3,i) = o3px(i)
82         xabsi(4,i) = h2ox(i)
83         xabsi(5,i) = h2x(i)
84         xabsi(6,i) = h2o2x(i)
85      end do
86
87      do i=1,nz
88         t2(i)=tx(i)
89         if(t0(i).lt.195.0) t0(i)=195.0
90         if(t0(i).gt.295.0) t0(i)=295.0
91         if(t2(i).lt.195.0) t2(i)=195.0
92         if(t2(i).gt.295.0) t2(i)=295.0
93      end do
94
95
96c    COLUMN CALCULATION
97
98      call column(co2x,o2x,o3px,h2x,h2ox,h2o2x,tx,nz,iz,zenit,co2colx
99     $,o2colx,o3pcolx,h2o2colx)
100
101      coltemp(nz)=co2colx(nz)*abs(t2(nz)-t0(nz))
102      dt(nz)=coltemp(nz)/co2colx(nz)
103      do i=nz-1,1,-1
104        coltemp(i)=coltemp(i+1)+
105     $         ( co2x(i) + co2x(i+1) ) * 0.5
106     $         * 1e5 * (iz(i+1)-iz(i)) * abs(t2(i)-t0(i))
107        dt(i)=coltemp(i)/co2colx(i)
108      end do
109
110
111c     Calculo la seccion eficaz de CO2 a la temperatura t0 a cada altura
112c CALCULATION OF CO2 CROSS-SECTION AT TEMPERATURE T0 AT EACH HEIGHT
113
114      do i=1,nz
115         do indexint=24,32
116           sigma(indexint,i)=co2crsc195(indexint-23)*
117     $         (1+((co2crsc295(indexint-23)/co2crsc195(indexint-23))-1.)
118     $         *(t0(i)-195.)/100.)
119           if(t0(i).ne.295.) then
120             alfa(indexint,i)=((co2crsc295(indexint-23)
121     $                        /sigma(indexint,i))-1.)/(295.-t0(i))
122           else
123             alfa(indexint,i)=0.
124           end if
125         end do
126      end do
127
128
129c     Comienza la interpolacion
130c INTERPOLATION BEGINS
131c************************************************
132c     o2,0.1,5.0
133c************************************************
134     
135      indexint=1
136      limdown=1e-20
137      limup=1e26
138       do i=1,nz
139         aux2(nz-i+1) = o2colx(i) + o3pcolx(i)
140       end do
141       do i=1,nz2
142         aux3(i) = jabsifotsint(indexint,2,nz2-i+1)
143         aux4(i) = c23(nz2-i+1)
144       enddo
145       call interpfast
146     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup)
147       do i=1,nz
148         jfotsout(indexint,2,i) = aux1(nz-i+1)
149       enddo
150
151c************************************************
152c     o3p,0.1,5.0
153c************************************************
154
155      indexint=1
156      limdown=1e-20
157      limup=1e26
158         do i=1,nz
159            aux2(nz-i+1) = o2colx(i) + o3pcolx(i)
160         end do
161         do i=1,nz2
162            aux3(i) = jabsifotsint(indexint,3,nz2-i+1)
163            aux4(i) = c23(nz2-i+1)
164         enddo
165         call interpfast
166     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup)
167         do i=1,nz
168            jfotsout(indexint,3,i) = aux1(nz-i+1)
169         enddo
170
171
172c************************************************
173c     h2,0.1,5.0
174c************************************************
175
176      indexint=1
177      limdown=1e-20
178      limup=1e26
179         do i=1,nz
180            aux2(nz-i+1) = o2colx(i) + o3pcolx(i)
181         end do
182         do i=1,nz2
183            aux3(i) = jabsifotsint(indexint,5,nz2-i+1)
184            aux4(i) = c23(nz2-i+1)
185         enddo
186         call interpfast
187     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
188         do i=1,nz
189            jfotsout(indexint,5,i) = aux1(nz-i+1)
190         enddo
191
192
193c************************************************
194c     interpola para el o2 entre 5 y 90.8nm
195c************************************************
196
197      limdown=1e-20
198      limup=1e12
199      do indexint=2,16
200         co2colx2(nz,indexint)= co2colx(nz)
201     $              * crscabsi2(1,indexint)
202         o2colx2(nz,indexint)= o2colx(nz)
203     $              * crscabsi2(2,indexint)
204         o3pcolx2(nz,indexint)= o3pcolx(nz)
205     $              * crscabsi2(3,indexint)
206         do i=nz-1,1,-1
207            co2colx2(i,indexint) = co2colx2(i+1,indexint) +
208     $         ( xabsi(1,i) + xabsi(1,i+1) ) * 0.5
209     $           * 1e5 * (iz(i+1)-iz(i)) * crscabsi2(1,indexint)
210            o2colx2(i,indexint) = o2colx2(i+1,indexint) +
211     $         ( xabsi(2,i) + xabsi(2,i+1) ) * 0.5
212     $         * 1e5 * (iz(i+1)-iz(i)) * crscabsi2(2,indexint)
213            o3pcolx2(i,indexint) = o3pcolx2(i+1,indexint) +
214     $         ( xabsi(3,i) + xabsi(3,i+1) ) * 0.5
215     $         * 1e5 * (iz(i+1)-iz(i)) * crscabsi2(3,indexint)
216         end do
217         do i=1,nz
218            aux2(nz-i+1) = co2colx2(i,indexint)+
219     $      o2colx2(i,indexint)+
220     $      o3pcolx2(i,indexint)
221         end do
222         do i=1,nz2
223            aux3(i) = jabsifotsint(indexint,2,nz2-i+1)
224            aux4(i) = c123(indexint,nz2-i+1)
225         enddo
226          call interpfast
227     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
228          do i=1,nz
229            jfotsout(indexint,2,i) = aux1(nz-i+1)
230          enddo
231      end do
232     
233
234c************************************************
235c     idem para o3p entre 5 y 90.8nm
236c************************************************
237
238      do indexint=2,16
239         do i=1,nz
240            aux2(nz-i+1) = co2colx2(i,indexint) +
241     $      o2colx2(i,indexint) +
242     $      o3pcolx2(i,indexint)
243         end do
244         do i=1,nz2
245            aux3(i) = jabsifotsint(indexint,3,nz2-i+1)
246            aux4(i) = c123(indexint,nz2-i+1)
247         enddo
248         call interpfast
249     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
250         do i=1,nz
251            jfotsout(indexint,3,i) = aux1(nz-i+1)
252         enddo
253      end do
254
255
256c************************************************
257c     idem para co2 entre 5 y 90.8nm
258c************************************************
259
260      do indexint=2,16
261         do i=1,nz
262            aux2(nz-i+1) = co2colx2(i,indexint) +
263     $      o2colx2(i,indexint) +
264     $      o3pcolx2(i,indexint)
265         end do
266         do i=1,nz2
267            aux3(i) = jabsifotsint(indexint,1,nz2-i+1)
268            aux4(i) = c123(indexint,nz2-i+1)
269         enddo
270         call interpfast
271     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
272         do i=1,nz
273            jfotsout(indexint,1,i) = aux1(nz-i+1)
274         enddo
275      end do
276
277
278c************************************************
279c     idem para h2 entre 5 y 80.6nm
280c************************************************
281
282      do indexint=2,15
283         do i=1,nz
284            aux2(nz-i+1) = co2colx2(i,indexint) +
285     $      o2colx2(i,indexint) +
286     $      o3pcolx2(i,indexint)
287         end do
288         do i=1,nabs
289         end do
290         do i=1,nz2
291            aux3(i) = jabsifotsint(indexint,5,nz2-i+1)
292            aux4(i) = c123(indexint,nz2-i+1)
293         enddo
294         call interpfast
295     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
296         do i=1,nz
297            jfotsout(indexint,5,i) = aux1(nz-i+1)
298         enddo
299      end do
300
301
302c************************************************
303c      co2,90.9,119.5
304c************************************************
305
306      limdown=1e-20
307      limup=1e26
308      do indexint=17,24
309         do i=1,nz
310            aux2(nz-i+1) = co2colx(i) + o2colx(i)
311         end do
312         do i=1,nz2
313            aux3(i) = jabsifotsint(indexint,1,nz2-i+1)
314            aux4(i) = c12(nz2-i+1)
315         enddo
316         call interpfast
317     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
318         do i=1,nz
319c           print*,indexint,i,aux2(i),aux1(i),aux3(i),aux4(i)
320           jfotsout(indexint,1,i) = aux1(nz-i+1)
321           if(indexint.eq.24) then
322           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
323             jfotsout(indexint,1,i) = aux1(nz-i+1)
324     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))*
325     $           (1+alfa(indexint,i)*(t2(i)-t0(i)))
326           else
327             jfotsout(indexint,1,i)=0.
328           end if
329           end if
330c           print*,indexint,i,(jfotsout(indexint,j,i),j=1,nabs)
331         enddo
332      end do
333
334
335c************************************************
336c     o2,90.9,119.5
337c************************************************
338     
339      do indexint=17,24
340         do i=1,nz
341            aux2(nz-i+1) = o2colx(i) + co2colx(i)
342         end do
343         do i=1,nz2
344            aux3(i) = jabsifotsint(indexint,2,nz2-i+1)
345            aux4(i) = c12(nz2-i+1)
346         enddo
347         call interpfast
348     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
349         do i=1,nz
350            jfotsout(indexint,2,i) = aux1(nz-i+1)
351            if(indexint.eq.24) then
352           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
353               jfotsout(indexint,2,i) = aux1(nz-i+1)
354     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
355               else
356                  jfotsout(indexint,2,i)=0.
357               end if
358            end if
359         enddo
360      end do
361
362
363c************************************************
364c      co2,119.6,202.5
365c************************************************
366
367      do indexint=25,31
368         do i=1,nz
369            aux2(nz-i+1) = co2colx(i) + o2colx(i)
370         end do
371         do i=1,nz2
372            aux3(i) = jabsifotsint(indexint,1,nz2-i+1)
373            aux4(i) = c12(nz2-i+1)
374         enddo
375         call interpfast
376     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
377         do i=1,nz
378           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
379               jfotsout(indexint,1,i) = aux1(nz-i+1)
380     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))*
381     $           (1+alfa(indexint,i)*(t2(i)-t0(i)))
382            else
383               jfotsout(indexint,1,i) = 0.
384            end if
385         enddo
386      end do
387
388
389c************************************************
390c      o2,119.6,202.5
391c************************************************
392
393      do indexint=25,31
394         do i=1,nz
395            aux2(nz-i+1) = co2colx(i) + o2colx(i)
396         end do
397         do i=1,nz2
398            aux3(i) = jabsifotsint(indexint,2,nz2-i+1)
399            aux4(i) = c12(nz2-i+1)
400         enddo
401         call interpfast
402     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
403         do i=1,nz
404           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
405            jfotsout(indexint,2,i) = aux1(nz-i+1)
406     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
407            else
408               jfotsout(indexint,2,i)=0.
409            end if
410         enddo
411      end do
412
413
414c************************************************
415c      h2o,119.6,202.5
416c************************************************
417
418      do indexint=25,31
419         do i=1,nz
420            aux2(nz-i+1) = co2colx(i) + o2colx(i)
421         end do
422         do i=1,nz2
423            aux3(i) = jabsifotsint(indexint,4,nz2-i+1)
424            aux4(i) = c12(nz2-i+1)
425         enddo
426         call interpfast
427     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
428         do i=1,nz
429           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
430            jfotsout(indexint,4,i) = aux1(nz-i+1)
431     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
432            else
433               jfotsout(indexint,4,i) = 0.
434            end if
435         enddo
436      end do
437
438
439c************************************************
440c      h2o2,119.6,202.5
441c************************************************
442
443      do indexint=25,31
444         do i=1,nz
445            aux2(nz-i+1) = co2colx(i) + o2colx(i)
446         end do
447         do i=1,nz2
448            aux3(i) = jabsifotsint(indexint,6,nz2-i+1)
449            aux4(i) = c12(nz2-i+1)
450         enddo
451         call interpfast
452     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
453         do i=1,nz
454           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
455            jfotsout(indexint,6,i) = aux1(nz-i+1)
456     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
457            else
458               jfotsout(indexint,6,i) = 0.
459            end if
460         enddo
461      end do
462
463
464c************************************************
465c     co2,202.6,210.0
466c************************************************
467
468      indexint=32
469         do i=1,nz
470            aux2(nz-i+1) = co2colx(i)
471         end do
472         do i=1,nz2
473            aux3(i) = jabsifotsint(indexint,1,nz2-i+1)
474            aux4(i) = c1(nz2-i+1)
475         enddo
476         call interpfast
477     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
478         do i=1,nz
479           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
480               jfotsout(indexint,1,i) = aux1(nz-i+1)
481     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))*
482     $           (1+alfa(indexint,i)*(t2(i)-t0(i)))
483            else
484               jfotsout(indexint,1,i)=0.
485            end if
486         enddo
487
488
489c************************************************
490c     h2o2,202.6,210.0
491c************************************************
492
493      indexint=32
494         do i=1,nz
495            aux2(nz-i+1) = co2colx(i)
496         end do
497         do i=1,nz2
498            aux3(i) = jabsifotsint(indexint,6,nz2-i+1)
499            aux4(i) = c1(nz2-i+1)
500         enddo
501         call interpfast
502     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
503         do i=1,nz
504           if(sigma(indexint,i)*alfa(indexint,i)*coltemp(i).lt.60.) then
505            jfotsout(indexint,6,i) = aux1(nz-i+1)
506     $            *exp(-sigma(indexint,i)*alfa(indexint,i)*coltemp(i))
507            else
508               jfotsout(indexint,6,i)=0.
509            end if
510         enddo
511
512
513c************************************************
514c     h2o2,210.1,337.5
515c************************************************
516
517      indexint=33
518      limdown=1.e-20
519      limup=1e25
520         do i=1,nz
521            aux2(nz-i+1) = h2o2colx(i)
522         end do
523         do i=1,nz2
524            aux3(i) = jabsifotsint(indexint,6,nz2-i+1)
525            aux4(i) = ch2o2(nz2-i+1)
526         enddo
527         call interpfast
528     $        ( aux1, aux2, nz,  aux3, aux4, nz2 , limdown , limup )
529         do i=1,nz
530            jfotsout(indexint,6,i) = aux1(nz-i+1)
531         enddo
532
533      return
534
535      end
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
Note: See TracBrowser for help on using the repository browser.