source: trunk/LMDZ.MARS/libf/aeronomars/moldiffcoeff_red.F @ 545

Last change on this file since 545 was 517, checked in by emillour, 13 years ago

Mars GCM: Adding the new molecular diffusion routines (old "moldiff" and "moldiffcoeff" are kept for now and if further tests are needed) by JYC:
moldiffcoeff_red.F, moldiff_red.F90, diffusion.h
Also updated "phymars/comdiurn.h", "aeronomars/conc.h" and "aeronomars/chimiedata.h" to be fixed/free form compatible.
File size: 7.5 KB
RevLine 
[517]1      subroutine moldiffcoeff_red(dij)
2
3       IMPLICIT NONE
4c=======================================================================
5c   subject:
6c   --------
7c   Computing molecular diffusion coefficients
8c   following Nair 94 (pg 131)
9c   author:  MAC 2002
10c   ------
11c
12c=======================================================================
13#include "dimensions.h"
14#include "dimphys.h"
15#include "callkeys.h"
16#include "comdiurn.h"
17#include "chimiedata.h"
18#include "tracer.h"
19#include "conc.h"
20#include "diffusion.h"
21
22c-----------------------------------------------------------------------
23c    Input/Output
24c    ------------
25c       integer,parameter :: ncompmoldiff = 12
26      real dij(ncompdiff,ncompdiff)
27
28c    Local variables:
29c    ---------------
30      INTEGER nq, n, nn, i,iq
31cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
32c     tracer numbering in the molecular diffusion
33cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
34
35      integer,parameter :: i_o   = 1
36      integer,parameter :: i_n2   = 2
37      integer,parameter :: i_co   = 3
38      integer,parameter :: i_ar  = 4
39      integer,parameter :: i_h2   = 5
40      integer,parameter :: i_h    = 6
41      integer,parameter :: i_o2   = 7
42      integer,parameter :: i_oh  = 8
43      integer,parameter :: i_ho2  = 9
44      integer,parameter :: i_h2o = 10
45      integer,parameter :: i_h2o2  = 11
46      integer,parameter :: i_o1d   = 12
47!      integer,parameter :: i_o3   = 13
48      integer,parameter :: i_n    = 13
49      integer,parameter :: i_no   = 14
50      integer,parameter :: i_no2  = 15
51!      integer,parameter :: i_n2d  = 17
52!      integer,parameter :: i_oplus = 18
53      integer,parameter :: i_co2    = 16
54!      integer,parameter :: i_oplus = 17
55!      integer,parameter :: i_hplus = 18
56
57! Tracer indexes in the GCM:
58      integer,save :: g_co2=0
59      integer,save :: g_co=0
60      integer,save :: g_o=0
61      integer,save :: g_o1d=0
62      integer,save :: g_o2=0
63      integer,save :: g_o3=0
64      integer,save :: g_h=0
65      integer,save :: g_h2=0
66      integer,save :: g_oh=0
67      integer,save :: g_ho2=0
68      integer,save :: g_h2o2=0
69      integer,save :: g_n2=0
70      integer,save :: g_ar=0
71      integer,save :: g_h2o=0
72      integer,save :: g_n=0
73      integer,save :: g_no=0
74      integer,save :: g_no2=0
75      integer,save :: g_n2d=0
76!      integer,save :: g_oplus=0
77!      integer,save :: g_hplus=0
78
79      integer,save :: gcmind(ncompdiff)
80
81      real dnh
82      logical,save :: firstcall=.true.
83
84! Initializations at first call (and some sanity checks)
85      if (firstcall) then
86        ! identify the indexes of the tracers we'll need
87        g_co2=igcm_co2
88        if (g_co2.eq.0) then
89          write(*,*) "moldiffcoeff: Error; no CO2 tracer !!!"
90          stop
91        endif
92        g_n2=igcm_n2
93        if (g_n2.eq.0) then
94          write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
95          stop
96        endif
97        g_ar=igcm_ar
98        if (g_ar.eq.0) then
99          write(*,*) "moldiffcoeff: Error; no Ar tracer !!!"
100          stop
101        endif       
102        g_h2=igcm_h2
103        if (g_h2.eq.0) then
104          write(*,*) "moldiffcoeff: Error; no H2 tracer !!!"
105          stop
106        endif
107        g_h=igcm_h
108        if (g_h.eq.0) then
109          write(*,*) "moldiffcoeff: Error; no H tracer !!!"
110          stop
111        endif
112        g_co=igcm_co
113        if (g_co.eq.0) then
114          write(*,*) "moldiffcoeff: Error; no CO tracer !!!"
115          stop
116        endif
117        g_o2=igcm_o2
118        if (g_o2.eq.0) then
119          write(*,*) "moldiffcoeff: Error; no O2 tracer !!!"
120          stop
121        endif
122        g_oh=igcm_oh
123        if (g_oh.eq.0) then
124          write(*,*) "moldiffcoeff: Error; no OH tracer !!!"
125          stop
126        endif
127        g_ho2=igcm_ho2
128        if (g_ho2.eq.0) then
129          write(*,*) "moldiffcoeff: Error; no HO2 tracer !!!"
130          stop
131        endif
132        g_h2o=igcm_h2o_vap
133        if (g_h2o.eq.0) then
134          write(*,*) "moldiffcoeff: Error; no H2O tracer !!!"
135          stop
136        endif
137        g_h2o2=igcm_h2o2
138        if (g_h2o2.eq.0) then
139          write(*,*) "moldiffcoeff: Error; no H2O2 tracer !!!"
140          stop
141        endif
142        g_o1d=igcm_o1d
143        if (g_h.eq.0) then
144          write(*,*) "moldiffcoeff: Error; no O1d tracer !!!"
145          stop
146        endif
147        g_o3=igcm_o3
148        if (g_o3.eq.0) then
149          write(*,*) "moldiffcoeff: Error; no O3 tracer !!!"
150          stop
151        endif
152        g_n=igcm_n
153        if (g_n.eq.0) then
154          write(*,*) "moldiffcoeff: Error; no N tracer !!!"
155          stop
156        endif
157        g_no=igcm_no
158        if (g_no.eq.0) then
159          write(*,*) "moldiffcoeff: Error; no NO tracer !!!"
160          stop
161        endif
162        g_no2=igcm_no2
163        if (g_no2.eq.0) then
164          write(*,*) "moldiffcoeff: Error; no NO2 tracer !!!"
165          stop
166        endif
167        g_n2d=igcm_n2d
168        if (g_n2d.eq.0) then
169          write(*,*) "moldiffcoeff: Error; no N2(D) tracer !!!"
170          stop
171        endif
172!        g_oplus=igcm_oplus
173!        if (g_oplus .eq. 0) then
174!        write(*,*) "moldiffcoeff: Error; no Oplus tracer !!!"
175!        stop
176!        endif
177!       g_hplus=igcm_hplus
178!        if (g_hplus .eq. 0) then
179!        write(*,*) "moldiffcoeff: Error; no Hplus tracer !!!"
180!        stop
181!        endif
182        g_o=igcm_o
183        if (g_o.eq.0) then
184          write(*,*) "moldiffcoeff: Error; no O tracer !!!"
185          stop
186        endif
187       
188c
189cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
190c    fill array to relate local indexes to gcm indexes
191cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
192
193        gcmind(i_co2)  =   g_co2
194        gcmind(i_n2)  =   g_n2
195        gcmind(i_ar)  =   g_ar
196        gcmind(i_h2) =   g_h2
197        gcmind(i_h)  =   g_h
198        gcmind(i_co)   =   g_co
199        gcmind(i_o2) =   g_o2
200        gcmind(i_oh)=   g_oh
201        gcmind(i_ho2)  =   g_ho2
202        gcmind(i_h2o) = g_h2o
203        gcmind(i_h2o2)= g_h2o2
204        gcmind(i_o1d) = g_o1d
205!        gcmind(i_o3) = g_o3
206        gcmind(i_n)= g_n
207        gcmind(i_no) = g_no
208        gcmind(i_no2) = g_no2
209!        gcmind(i_n2d) = g_n2d
210!        gcmind(i_oplus) =  g_oplus
211!        gcmind(i_hplus) = g_hplus
212        gcmind(i_o)   =   g_o
213
214c
215cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
216        firstcall= .false.
217      endif ! of if (firstcall)
218
219
220      dij(i_h2,i_co)   = 0.0000651
221      dij(i_h2,i_n2)   = 0.0000674
222      dij(i_h2,i_o2)   = 0.0000697
223      dij(i_h2,i_co2)  = 0.0000550
224      dij(i_h2,i_h2)   = 0.0
225      dij(i_h2,i_h)    = 0.0
226      dij(i_h2,i_h2o)  = 0.0    !0003
227      dij(i_h2,i_h2o2) = 0.0    !0003
228!      dij(i_h2,i_o3)   = 0.0   !0003
229      dij(i_h2,i_o)    = 0.0
230      dij(i_h2,i_ar)   = 0.0
231      dij(i_h2,i_n)    = 0.0
232
233c      dij(i_h,i_o)     = 0.0000144
234      dij(i_h,i_o)     = 0.000114
235
236       print*,'moldiffcoef: COEFF CALC'
237       open(56,file='coeffs.dat',status='unknown')
238      do n=1,ncompdiff
239        if (dij(i_h2,n).gt.0.0) then
240          do nn=n,ncompdiff
241            dij(nn,n)=dij(i_h2,n)
242     &                  *sqrt(mmol(g_h2)/mmol(gcmind(nn)))
243            if(n.eq.nn) dij(nn,n)=1.0
244            dij(n,nn)=dij(nn,n)
245          enddo
246        endif
247        if (dij(i_h2,n).eq.0.0) then
248          dnh=dij(i_h,i_o)*sqrt(mmol(g_o)/mmol(gcmind(n)))
249          do nn=n,ncompdiff
250            dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn)))
251            if(n.eq.nn) dij(nn,n)=1.0
252            dij(n,nn)=dij(nn,n)
253          enddo
254        endif
255      enddo
256
257      do n=1,ncompdiff
258        do nn=n,ncompdiff
259          write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
260        enddo
261      enddo
262      close(56)
263
264
265      return   
266      end
Note: See TracBrowser for help on using the repository browser.