source: trunk/LMDZ.GENERIC/libf/phystd/newsedim.F @ 486

Last change on this file since 486 was 486, checked in by rwordsworth, 13 years ago

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