source: trunk/LMDZ.MARS/libf/phymars/newsedim.F @ 259

Last change on this file since 259 was 147, checked in by aslmd, 13 years ago

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