source: trunk/LMDZ.MARS/libf/phymars/newsedim_mod.F @ 3026

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

LMDZ_MARS RV : Open_MP;
Put all the "save" variables as "!$OMP THREADPRIVATE" in phymars.
The code can now be tested, see README for more info
File size: 7.4 KB
Line 
1      MODULE newsedim_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6
7      SUBROUTINE newsedim(ngrid,nlay,naersize,nrhosize,ptimestep,
8     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,beta)
9      USE comcstfi_h, ONLY: r,g
10      IMPLICIT NONE
11
12c=======================================================================
13c
14c      Compute sedimentation of 1 tracer
15c      of radius rd (m) and density rho (kg.m-3)
16c
17c=======================================================================
18
19c-----------------------------------------------------------------------
20c   declarations:
21c   -------------
22
23c
24c   arguments:
25c   ----------
26
27      INTEGER,INTENT(IN) :: ngrid,nlay,naersize,nrhosize
28      REAL,INTENT(IN) :: ptimestep            ! pas de temps physique (s)
29      REAL,INTENT(IN) :: pplev(ngrid,nlay+1) ! pression aux inter-couches (Pa)
30      REAL,INTENT(IN) :: pt(ngrid,nlay) ! temperature au centre des couches (K)
31      real,intent(in) :: masse (ngrid,nlay) ! masse d'une couche (kg)
32      real,intent(in) :: epaisseur (ngrid,nlay)  ! epaisseur d'une couche (m)
33      real,intent(in) :: rd(naersize)             ! particle radius (m)
34      real,intent(in) :: rho(nrhosize)             ! particle density (kg.m-3)
35
36
37c    Traceurs :
38      real,intent(inout) :: pqi(ngrid,nlay)  ! traceur   (e.g. ?/kg)
39      real,intent(out) :: wq(ngrid,nlay+1)  ! flux de traceur durant timestep (?/m-2)
40      real,intent(in) :: beta ! correction for the shape of the particles
41                !   (see Murphy et al. JGR 1990 vol.95)
42                !   beta=1 for spheres
43                !   beta=0.85 for irregular particles
44                !   beta=0.5 for disk shaped particles
45     
46c   local:
47c   ------
48
49      INTEGER l,ig, k, i
50      REAL rfall,rhofall
51
52      LOGICAL,SAVE :: firstcall=.true.
53
54!$OMP THREADPRIVATE(firstcall)
55
56c    Traceurs :
57c    ~~~~~~~~
58      real traversee (ngrid,nlay)
59      real vstokes(ngrid,nlay)
60      real w(ngrid,nlay)
61      real ptop, dztop, Ep, Stra
62
63
64c    Physical constant
65c    ~~~~~~~~~~~~~~~~~
66c     Gas molecular viscosity (N.s.m-2)
67      real,parameter :: visc=1.e-5       ! CO2
68c     Effective gas molecular radius (m)
69      real,parameter :: molrad=2.2e-10   ! CO2
70
71c     local and saved variable
72      real,save :: a,b
73
74!$OMP THREADPRIVATE(a,b)
75
76
77c    ** un petit test de coherence
78c       --------------------------
79
80      ! AS: OK firstcall absolute
81      IF (firstcall) THEN
82
83         firstcall=.false.
84
85
86c       Preliminary calculations for sedimenation velocity :
87c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
88
89c       - Constant to compute stokes speed simple formulae
90c        (Vstokes =  b * rho* r**2   avec   b= (2/9) * rho * g / visc
91         b = 2./9. * g / visc
92      ENDIF ! of IF(firstcall)
93     
94c       - Constant  to compute gas mean free path
95c        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
96         a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23)
97
98c       - Correction to account for non-spherical shape (Murphy et al.  1990)
99         a = a * beta
100
101
102
103c-----------------------------------------------------------------------
104c    1. initialisation
105c    -----------------
106
107c     Sedimentation velocity (m/s)
108c     ~~~~~~~~~~~~~~~~~~~~~~
109c     (stokes law corrected for low pressure by the Cunningham
110c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
111
112        do  l=1,nlay
113          do ig=1, ngrid
114c           radius
115            if (naersize.eq.1) then
116              rfall=rd(1)
117            else
118              i=ngrid*(l-1)+ig
119              rfall=rd(i)
120            endif 
121c           density
122            if (nrhosize.eq.1) then
123              rhofall=rho(1)
124            else
125              i=ngrid*(l-1)+ig
126              rhofall=rho(i)
127            endif 
128c           vstokes
129            vstokes(ig,l) = b * rhofall * rfall*rfall *
130     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
131
132c           Layer crossing time (s) :
133            traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l)
134          end do
135        end do
136
137
138c     Calcul de la masse d'atmosphere correspondant a q transferee
139c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
140c     (e.g. on recherche le niveau  en dessous de laquelle le traceur
141c      va traverser le niveau intercouche l : "dztop" est sa hauteur
142c      au dessus de l (m), "ptop" est sa pression (Pa))
143
144      do  l=1,nlay
145        do ig=1, ngrid
146             
147             dztop = vstokes(ig,l)*  ptimestep
148             Ep=0
149             k=0
150
151            w(ig,l) = 0. !! JF+AS ajout initialisation
152c **************************************************************
153c            Simple Method
154
155cc             w(ig,l) =
156cc     &       (1.- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
157cccc           write(*,*) 'OK simple method l,w =', l, w(ig,l)
158cccc           write(*,*) 'OK simple method dztop =', dztop
159
160             w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
161             !!! Diagnostic: JF. Fix: AS. Date: 05/11
162             !!! Probleme arrondi avec la quantite ci-dessus
163             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
164             !!! ---> dans ce cas on utilise le developpement limite !
165             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2           
166             IF ( w(ig,l) .eq. 0. ) THEN
167                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
168             ELSE
169                w(ig,l) = w(ig,l) * pplev(ig,l) / g 
170             ENDIF
171
172
173c **************************************************************
174cccc         Complex method :
175            if (dztop.gt.epaisseur(ig,l)) then                !!!if on traverse plus d'une couche
176cccc            Cas ou on "epuise" la couche l : On calcule le flux
177cccc            Venant de dessus en tenant compte de la variation de Vstokes
178c **************************************************************
179               Ep= epaisseur(ig,l)
180               Stra= traversee(ig,l)
181               do while(dztop.gt.Ep.and.l+k+1.le.nlay)
182                 k=k+1
183                 dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra)
184                 Ep = Ep + epaisseur(ig,l+k)
185                 Stra = Stra + traversee(ig,l+k)
186               enddo
187               Ep = Ep - epaisseur(ig,l+k)
188             !ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
189
190             !!! JF+AS 05/11 Probleme arrondi potentiel, meme solution que ci-dessus
191             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
192             IF ( ptop .eq. 1. ) THEN
193                !PRINT*, 'newsedim: exposant trop petit ', ig, l
194                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
195             ELSE
196                ptop=pplev(ig,l+k) * ptop
197             ENDIF
198             w(ig,l) = (pplev(ig,l) - Ptop)/g
199
200            endif                !!!!!if complex method
201
202
203cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
204cc           write(*,*) 'OK new    method dztop =', dztop
205cc       if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop
206cc       if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep
207cc            if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l)
208cc            if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l)
209cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
210c **************************************************************
211
212
213        end do
214      end do
215      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
216c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
217c    &                wq(1,6),wq(1,7),pqi(1,6)
218
219
220      END SUBROUTINE newsedim
221     
222      END MODULE newsedim_mod
223
Note: See TracBrowser for help on using the repository browser.