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

Last change on this file since 3696 was 3663, checked in by gmilcareck, 4 months ago

Generic PCM:
Cleaning up the code. The gas constant and Avogadro number values
was not the same between the files in the model.
We choose the value recommended by the 2019 revision of the SI.
R=8.314463 J.K-1.mol-1 and NA=6.022141e23 mol-1.
GM
File size: 11.7 KB
Line 
1      MODULE newsedim_mod
2     
3      IMPLICIT NONE
4     
5      CONTAINS
6
7      SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep,
8     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,iq)
9     
10      use ioipsl_getin_p_mod, only: getin_p
11      use comcstfi_mod, only: r, g
12      use gases_h, only: gfrac, gnom, ngasmx
13      use gases_h, only: igas_CH4, igas_CO2, igas_H2, igas_H2O
14      use gases_h, only: igas_He, igas_N2
15      use tracer_h, only : igcm_h2o_ice
16      use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds     
17      use radii_mod, only: h2o_cloudrad
18
19      IMPLICIT NONE
20
21!==================================================================
22!     
23!     Purpose
24!     -------
25!      Calculates sedimentation of 1 tracer
26!      of radius rd (m) and density rho (kg.m-3)
27!     
28!==================================================================
29
30!-----------------------------------------------------------------------
31!   declarations
32!   ------------
33
34!   arguments
35!   ---------
36
37      integer,intent(in) :: ngrid  ! number of atmospheric columns
38      integer,intent(in) :: nlay  ! number of atmospheric layers
39      integer,intent(in) :: naersize   ! number of particle sizes (1 or number
40                                       ! of grid boxes)
41      real,intent(in) :: ptimestep ! physics time step (s)
42      real,intent(in) :: pplev(ngrid,nlay+1)   ! inter-layer pressures (Pa)
43      real,intent(in) :: pt(ngrid,nlay)    ! mid-layer temperatures (K)
44      real,intent(in) :: masse (ngrid,nlay) ! atmospheric mass (kg)
45      real,intent(in) :: epaisseur (ngrid,nlay)  ! thickness of atm. layers (m)
46      real,intent(in) :: rd(naersize) ! particle radius (m)
47      real,intent(in) :: rho ! particle density (kg.m-3)
48      real,intent(inout) :: pqi(ngrid,nlay)  ! tracer   (e.g. ?/kg)
49      real,intent(out) :: wq(ngrid,nlay+1)  ! flux of tracer during timestep (?/m-2)
50      integer,intent(in) :: iq ! tracer index
51
52c   local:
53c   ------
54
55      INTEGER l,ig, k, i, igas
56      REAL rfall, rsurf, Reynolds, Cd, zfrac
57      REAL reffh2oliq(ngrid,nlay), reffh2oice(ngrid,nlay)
58
59      LOGICAL,SAVE :: firstcall=.true.
60!$OMP THREADPRIVATE(firstcall)
61      LOGICAL,SAVE :: crystal_shape
62!$OMP THREADPRIVATE(crystal_shape)     
63
64c    Traceurs :
65c    ~~~~~~~~
66      real traversee (ngrid,nlay)
67      real vstokes(ngrid,nlay)
68      real w(ngrid,nlay)
69      real ptop, dztop, Ep, Stra
70
71
72c    Physical constant
73c    ~~~~~~~~~~~~~~~~~
74c     Gas molecular viscosity (N.s.m-2)
75      real, allocatable, save :: visc(:,:)
76!$OMP THREADPRIVATE(visc)     
77c     Effective gas molecular radius (m)
78      real,save :: molrad
79!$OMP THREADPRIVATE(molrad)
80
81c     local and saved variable
82      real,save :: a,b
83!$OMP THREADPRIVATE(a,b)
84
85c    ** un petit test de coherence
86c       --------------------------
87
88      !print*,'temporary fixed particle rad in newsedim!!'
89
90      IF (firstcall) THEN
91         firstcall=.false.
92
93c        Determination of the viscosity a(N.s.m-2) and the mean molecular radius (m)
94        write(*,*) "Calculation of the viscosity and the mean molecular"
95     &               ," radius from gases.def"
96         allocate(visc(ngrid,nlay))
97         visc(:,:)=0.0
98         molrad=0.
99         do igas=1, ngasmx
100           if(gfrac(igas).ge.0.0) then
101             if(igas.eq.igas_CO2) then
102               molrad = molrad + gfrac(igas)*2.2e-10                              ! CO2
103               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! CO2
104             elseif(igas.eq.igas_N2) then
105               molrad = molrad + gfrac(igas)*1.8e-10                              ! N2 (Kunze et al. 2022)
106               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! N2
107             elseif(igas.eq.igas_H2) then
108               molrad = molrad + gfrac(igas)*1.41e-10                             ! H2 (Ackerman & Marley 2001)
109               visc(:,:) = visc(:,:) + gfrac(igas)*2.0d-07*pt(:,:)**0.66          ! H2 (from Rosner 2000)
110             elseif(igas.eq.igas_H2O) then
111               molrad = molrad + gfrac(igas)*2.3e-10                              ! H2O (Crifo 1989 at 300K)
112               visc(:,:) = visc(:,:) + gfrac(igas)*8e-6                           ! H2O (Sengers & Kamgar-Parsi 1984)
113             elseif(igas.eq.igas_He) then
114               molrad = molrad + gfrac(igas)*1.1e-10                              ! He (Kunze et al. 2022)
115               visc(:,:) = visc(:,:) +                                            ! He
116     &                     gfrac(igas)*1.9e-5*(pt(:,:)/273.15)**0.7               ! He (Petersen 1970)
117             elseif(igas.eq.igas_CH4) then
118               molrad = molrad + gfrac(igas)*1.9e-10                              ! CH4 (Ismail et al. 2015)
119               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! CH4
120             else
121               molrad = molrad + gfrac(igas)*2.2e-10                              ! CO2 by default
122               visc(:,:) = visc(:,:) + gfrac(igas)*1.e-5                          ! CO2 by default
123               write(*,*) trim(gnom(igas))," is not included in"
124     &              ," newsedim, CO2 is used by default"
125             endif
126           endif
127         enddo
128         write(*,*) "visc(1,1)=",visc(1,1),"N.s.m-2; ",
129     &               "molrad=",molrad,"m"
130
131c Correction for non-spherical water ice particles
132         write(*,*) "Use non-spherical water ice particles",
133     &             " for the sedimentation ?"
134         crystal_shape=.false. !default
135         call getin_p("crystal_shape",crystal_shape)
136         write(*,*) " crystal_shape = ",crystal_shape
137
138
139!=======================================================================
140!     Preliminary calculations for sedimenation velocity
141
142!       - Constant to compute stokes speed simple formulae
143!        (Vstokes =  b / visc * rho* r**2   avec   b= (2/9) * rho * g
144         b = 2./9. * g
145
146!       - Constant  to compute gas mean free path
147!        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
148         a = 0.707*8.314463/(4*3.1416* molrad**2  * 6.022141e23)
149
150c       - Correction to account for non-spherical shape (Murphy et al.  1990)
151c   (correction = 0.85 for irregular particles, 0.5 for disk shaped particles)
152c        a = a    * 0.85
153
154
155      ENDIF
156
157c-----------------------------------------------------------------------
158c    1. initialisation
159c    -----------------
160
161c     Sedimentation velocity (m/s)
162c     ~~~~~~~~~~~~~~~~~~~~~~
163c     (stokes law corrected for low pressure by the Cunningham
164c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
165
166c Compute liquid and ice particle radii
167        if((iq.eq.igcm_h2o_ice).and.crystal_shape) then
168            call h2o_cloudrad(ngrid,nlay,pqi,reffh2oliq,reffh2oice)
169        endif
170
171
172        do  l=1,nlay
173          do ig=1, ngrid
174            if (naersize.eq.1) then
175              rfall=rd(1)
176            else
177              i=ngrid*(l-1)+ig
178              rfall=rd(i) ! how can this be correct!!?
179            endif 
180
181c Correction for non-spherical water ice particles           
182            if((iq.eq.igcm_h2o_ice).and.crystal_shape) then
183              zfrac= (pt(ig,l)-T_h2O_ice_clouds) /
184     &               (T_h2O_ice_liq-T_h2O_ice_clouds)
185              zfrac= MAX(zfrac, 0.0)
186              zfrac= MIN(zfrac, 1.0)
187              rsurf=max(reffh2oice(ig,l),45.6*reffh2oice(ig,l)**1.3)        ! surface radius (formula for rimed dendrites from Hemsfield 1977, transition at around 30 microns)
188              rsurf=1/(zfrac/reffh2oice(ig,l)+(1-zfrac)/rsurf)              ! radius giving the mean velocity between liquid and ice particles
189            else
190              rsurf=rfall
191            endif
192
193            vstokes(ig,l) = b / visc(ig,l) * rho * rfall**3 / rsurf *
194     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rsurf)
195
196c Correction for high Reynolds number
197            Reynolds=2. * pplev(ig,l) / r / pt(ig,l) *
198     &        rsurf * vstokes(ig,l) / visc(ig,l)
199            if(Reynolds.ge.1.0) then
200              do i=1,5
201              Cd=24. / Reynolds * (1. + 0.15 * Reynolds**0.687) +
202     &           0.42 / (1. + 42500 / Reynolds**1.16)                       ! (Formula from Bagheri 2018)
203              vstokes(ig,l) =(8./3.*pplev(ig,l)/r/pt(ig,l)*g*rfall**3 /
204     &           rsurf**2/rho/Cd *
205     &           (1.+1.333*(a*pt(ig,l)/pplev(ig,l))/rsurf))**0.5
206              Reynolds=2. * pplev(ig,l) / r / pt(ig,l) *
207     &           rsurf * vstokes(ig,l) / visc(ig,l)
208              enddo
209            endif   
210
211c           Layer crossing time (s) :
212            traversee(ig,l)= epaisseur(ig,l)/vstokes(ig,l)
213          end do
214        end do
215
216
217c     Calcul de la masse d'atmosphere correspondant a q transferee
218c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
219c     (e.g. on recherche le niveau  en dessous de laquelle le traceur
220c      va traverser le niveau intercouche l : "dztop" est sa hauteur
221c      au dessus de l (m), "ptop" est sa pression (Pa))
222      do  l=1,nlay
223        do ig=1, ngrid
224             
225             dztop = vstokes(ig,l)*  ptimestep
226             Ep=0
227             k=0
228           w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS)
229c **************************************************************
230c            Simple Method
231cc             w(ig,l) =
232cc     &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
233cc             write(*,*) 'OK simple method l,w =', l, w(ig,l)
234cc            write(*,*) 'OK simple method dztop =', dztop
235           w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
236             !!! Diagnostic: JF. Fix: AS. Date: 05/11
237             !!! Probleme arrondi avec la quantite ci-dessus
238             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
239             !!! ---> dans ce cas on utilise le developpement limite !
240             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2 
241
242             IF ( w(ig,l) .eq. 0. ) THEN
243                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
244             ELSE
245                w(ig,l) = w(ig,l) * pplev(ig,l) / g
246             ENDIF
247! LK borrowed simple method from Mars model (AS/JF)
248
249!**************************************************************
250cccc         Complex method :
251           if (dztop.gt.epaisseur(ig,l)) then
252cccc            Cas ou on "epuise" la couche l : On calcule le flux
253cccc            Venant de dessus en tenant compte de la variation de Vstokes
254
255               Ep= epaisseur(ig,l)
256               Stra= traversee(ig,l)
257               do while(dztop.gt.Ep.and.l+k+1.le.nlay)
258                 k=k+1
259                 dztop= Ep + vstokes(ig,l+k)*(ptimestep -Stra)
260                 Ep = Ep + epaisseur(ig,l+k)
261                 Stra = Stra + traversee(ig,l+k)
262               enddo
263               Ep = Ep - epaisseur(ig,l+k)
264!             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
265             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
266             IF ( ptop .eq. 1. ) THEN
267                !PRINT*, 'newsedim: exposant trop petit ', ig, l
268                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
269             ELSE
270                ptop=pplev(ig,l+k) * ptop
271             ENDIF
272
273             w(ig,l) = (pplev(ig,l) - ptop)/g
274
275            endif   !!! complex method
276c
277cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
278cc           write(*,*) 'OK new    method dztop =', dztop
279cc       if(l.eq.7)write(*,*)'l=7,k,pplev,Ptop',pplev(ig,l),Ptop
280cc       if(l.eq.7)write(*,*)'l=7,dztop,Ep',dztop,Ep
281cc            if(l.eq.6)write(*,*)'l=6,k, w',k, w(1,l)
282cc            if(l.eq.7)write(*,*)'l=7,k, w',k, w(1,l)
283cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
284c **************************************************************
285
286        end do
287      end do
288
289      call vlz_fi(ngrid,nlay,pqi,2.,masse,w,wq)
290c         write(*,*) ' newsed: wq(6), wq(7), q(6)',
291c    &                wq(1,6),wq(1,7),pqi(1,6)
292
293      END SUBROUTINE newsedim
294     
295      END MODULE newsedim_mod
Note: See TracBrowser for help on using the repository browser.