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

Last change on this file was 3015, checked in by emillour, 17 months ago

Mars PCM:
Code cleanup in diffusion, turn variables from diffusion.h into module
variables in moldiff_red.F90 (turn it into a module in the process).
Also turn moldiffcoeff_red.F and thermosphere.F into modules.
EM
File size: 8.8 KB
RevLine 
[3015]1      MODULE moldiffcoeff_red_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6
[645]7      subroutine moldiffcoeff_red(dij,indic,gcmind,ncompdiff2)
[517]8
[1036]9       use tracer_mod, only: nqmx, noms, mmol
[517]10       IMPLICIT NONE
11c=======================================================================
12c   subject:
13c   --------
14c   Computing molecular diffusion coefficients
15c   following Nair 94 (pg 131)
16c   author:  MAC 2002
17c   ------
18c
19c=======================================================================
[3015]20      include "callkeys.h"
[517]21
22c-----------------------------------------------------------------------
23c    Input/Output
24c    ------------
25c       integer,parameter :: ncompmoldiff = 12
[645]26        integer ncompdiff2
27        integer gcmind(ncompdiff2)
28      real dij(ncompdiff2,ncompdiff2)
29      integer indic(nqmx)
[517]30
31c    Local variables:
32c    ---------------
33      INTEGER nq, n, nn, i,iq
34cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
35c     tracer numbering in the molecular diffusion
36cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
37
[645]38      real :: dijh2co,dijh2n2,dijh2co2,dijh2o2,dijho,dijref
39!       integer :: i_h2,i_h,i_o
40        integer :: g_h2,g_h,g_o
41!      integer,parameter :: i_o   = 1
42!      integer,parameter :: i_n2   = 2
43!      integer,parameter :: i_co   = 3
44!      integer,parameter :: i_ar  = 4
45!      integer,parameter :: i_h2   = 5
46!      integer,parameter :: i_h    = 6
47!      integer,parameter :: i_o2   = 7
48!      integer,parameter :: i_oh  = 8
49!      integer,parameter :: i_ho2  = 9
50!      integer,parameter :: i_h2o = 10
51!      integer,parameter :: i_h2o2  = 11
52!      integer,parameter :: i_o1d   = 12
[517]53!      integer,parameter :: i_o3   = 13
[645]54!      integer,parameter :: i_n    = 13
55!      integer,parameter :: i_no   = 14
56!      integer,parameter :: i_no2  = 15
[517]57!      integer,parameter :: i_n2d  = 17
58!      integer,parameter :: i_oplus = 18
[645]59!      integer,parameter :: i_co2    = 16
[517]60!      integer,parameter :: i_oplus = 17
61!      integer,parameter :: i_hplus = 18
62
63! Tracer indexes in the GCM:
[645]64!      integer,save :: g_co2=0
65!      integer,save :: g_co=0
66!      integer,save :: g_o=0
67!      integer,save :: g_o1d=0
68!      integer,save :: g_o2=0
69!      integer,save :: g_o3=0
70!      integer,save :: g_h=0
71!      integer,save :: g_h2=0
72!      integer,save :: g_oh=0
73!      integer,save :: g_ho2=0
74!      integer,save :: g_h2o2=0
75!      integer,save :: g_n2=0
76!      integer,save :: g_ar=0
77!      integer,save :: g_h2o=0
78!      integer,save :: g_n=0
79!      integer,save :: g_no=0
80!      integer,save :: g_no2=0
81!      integer,save :: g_n2d=0
[517]82!      integer,save :: g_oplus=0
83!      integer,save :: g_hplus=0
84
[645]85!      integer,save :: gcmind(ncompdiff)
[517]86
87      real dnh
88      logical,save :: firstcall=.true.
[2615]89
90!$OMP THREADPRIVATE(firstcall)
91
[710]92      logical,parameter :: outputcoeffs=.false. ! to output 'coeffs.dat' file,
93                                                ! set outputcoeffs=.true.
[517]94
95! Initializations at first call (and some sanity checks)
96      if (firstcall) then
97        ! identify the indexes of the tracers we'll need
[645]98!        g_co2=igcm_co2
99!        if (g_co2.eq.0) then
100!          write(*,*) "moldiffcoeff: Error; no CO2 tracer !!!"
101!          stop
102!        endif
103!        g_n2=igcm_n2
104!        if (g_n2.eq.0) then
105!          write(*,*) "moldiffcoeff: Error; no N2 tracer !!!"
106!          stop
107!        endif
108!        g_ar=igcm_ar
109!        if (g_ar.eq.0) then
110!          write(*,*) "moldiffcoeff: Error; no Ar 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_h=igcm_h
119!        if (g_h.eq.0) then
120!          write(*,*) "moldiffcoeff: Error; no H tracer !!!"
121!          stop
122!        endif
123!        g_co=igcm_co
124!        if (g_co.eq.0) then
125!          write(*,*) "moldiffcoeff: Error; no CO tracer !!!"
126!          stop
127!        endif
128!        g_o2=igcm_o2
129!        if (g_o2.eq.0) then
130!          write(*,*) "moldiffcoeff: Error; no O2 tracer !!!"
131!          stop
132!        endif
133!        g_oh=igcm_oh
134!        if (g_oh.eq.0) then
135!          write(*,*) "moldiffcoeff: Error; no OH tracer !!!"
136!          stop
137!        endif
138!        g_ho2=igcm_ho2
139!        if (g_ho2.eq.0) then
140!          write(*,*) "moldiffcoeff: Error; no HO2 tracer !!!"
141!          stop
142!        endif
143!        g_h2o=igcm_h2o_vap
144!        if (g_h2o.eq.0) then
145!          write(*,*) "moldiffcoeff: Error; no H2O tracer !!!"
146!          stop
147!        endif
148!        g_h2o2=igcm_h2o2
149!        if (g_h2o2.eq.0) then
150!          write(*,*) "moldiffcoeff: Error; no H2O2 tracer !!!"
151!          stop
152!        endif
153!        g_o1d=igcm_o1d
154!        if (g_h.eq.0) then
155!          write(*,*) "moldiffcoeff: Error; no O1d tracer !!!"
156!          stop
157!        endif
158!        g_o3=igcm_o3
159!        if (g_o3.eq.0) then
160!          write(*,*) "moldiffcoeff: Error; no O3 tracer !!!"
161!          stop
162!        endif
163!        g_n=igcm_n
164!        if (g_n.eq.0) then
165!          write(*,*) "moldiffcoeff: Error; no N tracer !!!"
166!          stop
167!        endif
168!        g_no=igcm_no
169!        if (g_no.eq.0) then
170!          write(*,*) "moldiffcoeff: Error; no NO tracer !!!"
171!          stop
172!        endif
173!        g_no2=igcm_no2
174!        if (g_no2.eq.0) then
175!          write(*,*) "moldiffcoeff: Error; no NO2 tracer !!!"
176!          stop
177!        endif
178!        g_n2d=igcm_n2d
179!        if (g_n2d.eq.0) then
180!          write(*,*) "moldiffcoeff: Error; no N2(D) tracer !!!"
181!          stop
182!        endif
[517]183!        g_oplus=igcm_oplus
184!        if (g_oplus .eq. 0) then
185!        write(*,*) "moldiffcoeff: Error; no Oplus tracer !!!"
186!        stop
187!        endif
188!       g_hplus=igcm_hplus
189!        if (g_hplus .eq. 0) then
190!        write(*,*) "moldiffcoeff: Error; no Hplus tracer !!!"
191!        stop
192!        endif
[645]193!        g_o=igcm_o
194!        if (g_o.eq.0) then
195!          write(*,*) "moldiffcoeff: Error; no O tracer !!!"
196!          stop
197!        endif
[517]198       
199c
200cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
201c    fill array to relate local indexes to gcm indexes
202cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
203
[645]204!        gcmind(i_co2)  =   g_co2
205!        gcmind(i_n2)  =   g_n2
206!        gcmind(i_ar)  =   g_ar
207!        gcmind(i_h2) =   g_h2
208!        gcmind(i_h)  =   g_h
209!        gcmind(i_co)   =   g_co
210!        gcmind(i_o2) =   g_o2
211!        gcmind(i_oh)=   g_oh
212!        gcmind(i_ho2)  =   g_ho2
213!        gcmind(i_h2o) = g_h2o
214!        gcmind(i_h2o2)= g_h2o2
215!        gcmind(i_o1d) = g_o1d
[517]216!        gcmind(i_o3) = g_o3
[645]217!        gcmind(i_n)= g_n
218!        gcmind(i_no) = g_no
219!        gcmind(i_no2) = g_no2
[517]220!        gcmind(i_n2d) = g_n2d
221!        gcmind(i_oplus) =  g_oplus
222!        gcmind(i_hplus) = g_hplus
[645]223!        gcmind(i_o)   =   g_o
[517]224
225c
226cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
227        firstcall= .false.
228      endif ! of if (firstcall)
229
[645]230        dijh2co = 0.0000651
231        dijh2n2 = 0.0000674
232        dijh2o2 = 0.0000697
233        dijh2co2 = 0.0000550
234        dijho = 0.000114
[517]235
[645]236!      dij(i_h2,i_co)   = 0.0000651
237!      dij(i_h2,i_n2)   = 0.0000674
238!      dij(i_h2,i_o2)   = 0.0000697
239!      dij(i_h2,i_co2)  = 0.0000550
240!      dij(i_h2,i_h2)   = 0.0
241!      dij(i_h2,i_h)    = 0.0
242!      dij(i_h2,i_h2o)  = 0.0   !0003
243!      dij(i_h2,i_h2o2) = 0.0   !0003
[517]244!      dij(i_h2,i_o3)   = 0.0   !0003
[645]245!      dij(i_h2,i_o)    = 0.0
246!      dij(i_h2,i_ar)   = 0.0
247!      dij(i_h2,i_n)    = 0.0
[517]248
[645]249!c      dij(i_h,i_o)     = 0.0000144
250!      dij(i_h,i_o)     = 0.000114
[517]251
[645]252! find h2, h and o index in gcm
253! these species are used to define the diffusion coefficients
254
255        do n=1,nqmx
256        if (noms(n) .eq. 'h2') g_h2=n
257        if (noms(n) .eq. 'h') g_h=n
258        if (noms(n) .eq. 'o') g_o=n
259        enddo
[710]260        print*,'moldiffcoeff_red: gh2',g_h2,g_h,g_o
[645]261
[710]262       print*,'moldiffcoeff_red: COEFF CALC'
263
[645]264      do n=1,ncompdiff2
265        dijref=0.
266        if (noms(gcmind(n)) .eq. 'co') dijref=dijh2co
267        if (noms(gcmind(n)) .eq. 'n2') dijref=dijh2n2
268        if (noms(gcmind(n)) .eq. 'o2') dijref=dijh2o2
269        if (noms(gcmind(n)) .eq. 'co2') dijref=dijh2co2
270!       print*,'test',n,dijref
271        if (dijref .gt. 0.0) then
272          do nn=n,ncompdiff2
273            dij(nn,n)=dijref
[517]274     &                  *sqrt(mmol(g_h2)/mmol(gcmind(nn)))
275            if(n.eq.nn) dij(nn,n)=1.0
276            dij(n,nn)=dij(nn,n)
277          enddo
278        endif
[645]279        if (dijref .eq. 0.0) then
280        dijref=dijho
281          dnh=dijref*sqrt(mmol(g_o)/mmol(gcmind(n)))
282          do nn=n,ncompdiff2
[517]283            dij(nn,n)=dnh*sqrt(mmol(g_h)/mmol(gcmind(nn)))
284            if(n.eq.nn) dij(nn,n)=1.0
285            dij(n,nn)=dij(nn,n)
286          enddo
287        endif
288      enddo
289
[710]290      if (outputcoeffs) then
291       ! output coefficients in 'coeffs.dat' file
292       open(56,file='coeffs.dat',status='unknown')
293       do n=1,ncompdiff2
[645]294        do nn=n,ncompdiff2
[517]295          write(56,*) n,nn,dij(n,nn)    !*1.e5/1.381e-23/(273**1.75)
296        enddo
[710]297       enddo
298       close(56)
299      endif
[517]300
[3015]301      end subroutine moldiffcoeff_red
[645]302     
[3015]303      END MODULE moldiffcoeff_red_mod
Note: See TracBrowser for help on using the repository browser.