source: trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff.F @ 3026

Last change on this file since 3026 was 2615, checked in by romain.vande, 3 years ago

LMDZ_MARS RV : Open_MP;
Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars
File size: 6.5 KB
Line 
1      subroutine moldiffcoeff(dij)
2
3      use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
4     &                      igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
5     &                      igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar,
6     &                      igcm_h2o_vap, mmol
7
8       IMPLICIT NONE
9c=======================================================================
10c   subject:
11c   --------
12c   Computing molecular diffusion coefficients
13c   following Nair 94 (pg 131)
14c   author:  MAC 2002
15c   ------
16c
17c=======================================================================
18#include "callkeys.h"
19
20c-----------------------------------------------------------------------
21c    Input/Output
22c    ------------
23       integer,parameter :: ncompmoldiff = 14
24      real dij(ncompmoldiff,ncompmoldiff)
25
26c    Local variables:
27c    ---------------
28      INTEGER nq, n, nn, i,iq
29cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
30c     tracer numbering in the molecular diffusion
31cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
32c  Atomic oxygen must always be the LAST species of the list as
33c it is the dominant species at high altitudes. 
34      integer,parameter :: i_co   = 1
35      integer,parameter :: i_n2   = 2
36      integer,parameter :: i_o2   = 3
37      integer,parameter :: i_co2  = 4
38      integer,parameter :: i_h2   = 5
39      integer,parameter :: i_h    = 6
40      integer,parameter :: i_oh   = 7
41      integer,parameter :: i_ho2  = 8
42      integer,parameter :: i_h2o  = 9
43      integer,parameter :: i_h2o2 = 10
44      integer,parameter :: i_o1d  = 11
45      integer,parameter :: i_o3   = 12
46      integer,parameter :: i_ar   = 13
47      integer,parameter :: i_o    = 14
48
49! Tracer indexes in the GCM:
50      integer,save :: g_co2=0
51      integer,save :: g_co=0
52      integer,save :: g_o=0
53      integer,save :: g_o1d=0
54      integer,save :: g_o2=0
55      integer,save :: g_o3=0
56      integer,save :: g_h=0
57      integer,save :: g_h2=0
58      integer,save :: g_oh=0
59      integer,save :: g_ho2=0
60      integer,save :: g_h2o2=0
61      integer,save :: g_n2=0
62      integer,save :: g_ar=0
63      integer,save :: g_h2o=0
64
65!$OMP THREADPRIVATE(g_co2,g_co,g_o,g_o1d,g_o2,g_o3,g_h,g_h2,g_oh)
66!$OMP THREADPRIVATE(g_ho2,g_h2o2,g_n2,g_ar,g_h2o)
67
68      integer,save :: gcmind(ncompmoldiff)
69
70      real dnh
71      logical,save :: firstcall=.true.
72
73!$OMP THREADPRIVATE(gcmind,firstcall)
74
75! Initializations at first call (and some sanity checks)
76      if (firstcall) then
77        ! identify the indexes of the tracers we'll need
78        g_co2=igcm_co2
79        if (g_co2.eq.0) then
80          write(*,*) "moldiffcoeff: Error; no CO2 tracer !!!"
81          stop
82        endif
83        g_co=igcm_co
84        if (g_co.eq.0) then
85          write(*,*) "moldiffcoeff: Error; no CO tracer !!!"
86          stop
87        endif
88        g_o=igcm_o
89        if (g_o.eq.0) then
90          write(*,*) "moldiffcoeff: Error; no O tracer !!!"
91          stop
92        endif
93        g_o1d=igcm_o1d
94        if (g_o1d.eq.0) then
95          write(*,*) "moldiffcoeff: Error; no O1D tracer !!!"
96          stop
97        endif
98        g_o2=igcm_o2
99        if (g_o2.eq.0) then
100          write(*,*) "moldiffcoeff: Error; no O2 tracer !!!"
101          stop
102        endif
103        g_o3=igcm_o3
104        if (g_o3.eq.0) then
105          write(*,*) "moldiffcoeff: Error; no O3 tracer !!!"
106          stop
107        endif
108        g_h=igcm_h
109        if (g_h.eq.0) then
110          write(*,*) "moldiffcoeff: Error; no H tracer !!!"
111          stop
112        endif
113        g_h2=igcm_h2
114        if (g_h2.eq.0) then
115          write(*,*) "moldiffcoeff: Error; no H2 tracer !!!"
116          stop
117        endif
118        g_oh=igcm_oh
119        if (g_oh.eq.0) then
120          write(*,*) "moldiffcoeff: Error; no OH tracer !!!"
121          stop
122        endif
123        g_ho2=igcm_ho2
124        if (g_ho2.eq.0) then
125          write(*,*) "moldiffcoeff: Error; no HO2 tracer !!!"
126          stop
127        endif
128        g_h2o2=igcm_h2o2
129        if (g_h2o2.eq.0) then
130          write(*,*) "moldiffcoeff: Error; no H2O2 tracer !!!"
131          stop
132        endif
133        g_n2=igcm_n2
134        if (g_n2.eq.0) then
135          write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
136          stop
137        endif
138        g_ar=igcm_ar
139        if (g_ar.eq.0) then
140          write(*,*) "moldiffcoeff: Error; no AR tracer !!!"
141          stop
142        endif
143        g_h2o=igcm_h2o_vap
144        if (g_h2o.eq.0) then
145          write(*,*) "moldiffcoeff: Error; no water vapor tracer !!!"
146          stop
147        endif
148
149c
150cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
151c    fill array to relate local indexes to gcm indexes
152cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
153
154        gcmind(i_co)  =   g_co
155        gcmind(i_n2)  =   g_n2
156        gcmind(i_o2)  =   g_o2
157        gcmind(i_co2) =   g_co2
158        gcmind(i_h2)  =   g_h2
159        gcmind(i_h)   =   g_h
160        gcmind(i_oh)  =   g_oh
161        gcmind(i_ho2) =   g_ho2
162        gcmind(i_h2o) =   g_h2o
163        gcmind(i_h2o2)=   g_h2o2
164        gcmind(i_o1d) =   g_o1d
165        gcmind(i_o3)  =   g_o3
166        gcmind(i_o)   =   g_o
167        gcmind(i_ar)   =  g_ar
168c
169cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
170        firstcall= .false.
171      endif ! of if (firstcall)
172
173
174      dij(i_h2,i_co)   = 0.0000651
175      dij(i_h2,i_n2)   = 0.0000674
176      dij(i_h2,i_o2)   = 0.0000697
177      dij(i_h2,i_co2)  = 0.0000550
178      dij(i_h2,i_h2)   = 0.0
179      dij(i_h2,i_h)    = 0.0
180      dij(i_h2,i_oh)   = 0.0    !0003
181      dij(i_h2,i_ho2)  = 0.0    !0003
182      dij(i_h2,i_h2o)  = 0.0    !0003
183      dij(i_h2,i_h2o2) = 0.0    !0003
184      dij(i_h2,i_o1d)  = 0.0
185      dij(i_h2,i_o3)   = 0.0    !0003
186      dij(i_h2,i_o)    = 0.0
187      dij(i_h2,i_ar)   = 0.0
188
189c      dij(i_h,i_o)     = 0.0000144
190      dij(i_h,i_o)     = 0.000114
191
192       print*,'moldiffcoef: COEFF CALC'
193       open(56,file='coeffs.dat',status='unknown')
194      do n=1,ncompmoldiff
195        if (dij(i_h2,n).gt.0.0) then
196          do nn=n,ncompmoldiff
197            dij(nn,n)=dij(i_h2,n)
198     &                  *sqrt(mmol(g_h2)/mmol(gcmind(nn)))
199            if(n.eq.nn) dij(nn,n)=1.0
200            dij(n,nn)=dij(nn,n)
201          enddo
202        endif
203        if (dij(i_h2,n).eq.0.0) then
204          dnh=dij(i_h,i_o)*sqrt(mmol(g_o)/mmol(gcmind(n)))
205          do nn=n,ncompmoldiff
206            dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn)))
207            if(n.eq.nn) dij(nn,n)=1.0
208            dij(n,nn)=dij(nn,n)
209          enddo
210        endif
211      enddo
212
213      do n=1,ncompmoldiff
214        do nn=n,ncompmoldiff
215          write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
216        enddo
217      enddo
218      close(56)
219
220
221      return   
222      end
Note: See TracBrowser for help on using the repository browser.