Changeset 2670 for trunk/LMDZ.GENERIC/libf/phystd/newsedim.F

Timestamp:

Apr 26, 2022, 3:37:00 PM (3 years ago)

Author:

bclmd

Message:

Adding calculation of viscosity/molrad, general stokes law and non-spherical water ice particles in newsedim.F

File:

• trunk/LMDZ.GENERIC/libf/phystd/newsedim.F (modified) (8 diffs)

trunk/LMDZ.GENERIC/libf/phystd/newsedim.F

@@ -1 +1 @@
       SUBROUTINE newsedim(ngrid,nlay,naersize,ptimestep,
-     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq)
+     &  pplev,masse,epaisseur,pt,rd,rho,pqi,wq,iq)
       
+      use ioipsl_getin_p_mod, only: getin_p
       use comcstfi_mod, only: r, g
+      use gases_h
+      use tracer_h, only : igcm_h2o_ice
+      use watercommon_h, only: T_h2O_ice_liq,T_h2O_ice_clouds     
+      use radii_mod, only: h2o_cloudrad
+
       IMPLICIT NONE
 
@@ -34 +40 @@
       real,intent(inout) :: pqi(ngrid,nlay)  ! tracer   (e.g. ?/kg)
       real,intent(out) :: wq(ngrid,nlay+1)  ! flux of tracer during timestep (?/m-2)
-      
+      integer,intent(in) :: iq ! tracer index
+
 c   local:
 c   ------
 
-      INTEGER l,ig, k, i
-      REAL rfall
+      INTEGER l,ig, k, i, igas
+      REAL rfall, rsurf, Reynolds, Cd, zfrac
+      REAL reffh2oliq(ngrid,nlay), reffh2oice(ngrid,nlay)
 
       LOGICAL,SAVE :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
+      LOGICAL,SAVE :: crystal_shape
+!$OMP THREADPRIVATE(crystal_shape)      
 
 c    Traceurs :
@@ -55 +65 @@
 c    ~~~~~~~~~~~~~~~~~
 c     Gas molecular viscosity (N.s.m-2)
-      real,parameter :: visc=1.e-5       ! CO2
+      real, allocatable, save :: visc(:,:)
+!$OMP THREADPRIVATE(visc)      
 c     Effective gas molecular radius (m)
-      real,parameter :: molrad=2.2e-10   ! CO2
+      real,save :: molrad
 
 c     local and saved variable
@@ -66 +77 @@
 c       --------------------------
 
-
       !print*,'temporary fixed particle rad in newsedim!!'
 
@@ -72 +82 @@
          firstcall=.false.
 
+c        Determination of the viscosity a(N.s.m-2) and the mean molecular radius (m)
+        write(*,*) "Calculation of the viscosity and the mean molecular"
+     &               ," radius from gases.def"
+         allocate(visc(ngrid,nlay))
+         visc(:,:)=0.0
+         molrad=0.
+         do igas=1, ngasmx
+           if(gfrac(igas).ge.0.0) then
+             if(igas.eq.igas_CO2) then
+               molrad = molrad + gfrac(igas)*2.2e-10                              ! CO2 
+               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! CO2
+             elseif(igas.eq.igas_N2) then
+               molrad = molrad + gfrac(igas)*1.8e-10                              ! N2 (Kunze et al. 2022)
+               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! N2
+             elseif(igas.eq.igas_H2) then
+               molrad = molrad + gfrac(igas)*1.41e-10                             ! H2 (Ackerman & Marley 2001)
+               visc(:,:) = visc(:,:) + gfrac(igas)*2.0d-07*pt(:,:)**0.66          ! H2 (from Rosner 2000)
+             elseif(igas.eq.igas_H2O) then 
+               molrad = molrad + gfrac(igas)*2.3e-10                              ! H2O (Crifo 1989 at 300K)
+               visc(:,:) = visc(:,:) + gfrac(igas)*8e-6                           ! H2O (Sengers & Kamgar-Parsi 1984)
+             elseif(igas.eq.igas_He) then
+               molrad = molrad + gfrac(igas)*1.1e-10                              ! He (Kunze et al. 2022)
+               visc(:,:) = visc(:,:) +                                            ! He
+     &                     gfrac(igas)*1.9e-5*(pt(:,:)/273.15)**0.7               ! He (Petersen 1970)
+             elseif(igas.eq.igas_CH4) then
+               molrad = molrad + gfrac(igas)*1.9e-10                              ! CH4 (Ismail et al. 2015)
+               visc(:,:) = visc(:,:) + gfrac(igas)*1.0e-5                         ! CH4 
+             else
+               molrad = molrad + gfrac(igas)*2.2e-10                              ! CO2 by default
+               visc(:,:) = visc(:,:) + gfrac(igas)*1.e-5                          ! CO2 by default
+               write(*,*) trim(gnom(igas))," is not included in"
+     &              ," newsedim, CO2 is used by default"
+             endif
+           endif
+         enddo
+         write(*,*) "visc(1,1)=",visc(1,1),"N.s.m-2; ",
+     &               "molrad=",molrad,"m"
+
+c Correction for non-spherical water ice particles 
+         write(*,*) "Use non-spherical water ice particles", 
+     &             " for the sedimentation ?"
+         crystal_shape=.false. !default
+         call getin_p("crystal_shape",crystal_shape)
+         write(*,*) " crystal_shape = ",crystal_shape
 
 
@@ -78 +132 @@
 
 !       - Constant to compute stokes speed simple formulae
-!        (Vstokes =  b * rho* r**2   avec   b= (2/9) * rho * g / visc
-         b = 2./9. * g / visc
+!        (Vstokes =  b / visc * rho* r**2   avec   b= (2/9) * rho * g
+         b = 2./9. * g
 
 !       - Constant  to compute gas mean free path
@@ -100 +154 @@
 c     (stokes law corrected for low pressure by the Cunningham
 c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
-       
+
+c Compute liquid and ice particle radii
+        if((iq.eq.igcm_h2o_ice).and.crystal_shape) then
+            call h2o_cloudrad(ngrid,nlay,pqi,reffh2oliq,reffh2oice)
+        endif
+
+
         do  l=1,nlay
           do ig=1, ngrid
@@ -110 +170 @@
             endif  
 
-            vstokes(ig,l) = b * rho * rfall**2 *
-     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
+c Correction for non-spherical water ice particles            
+            if((iq.eq.igcm_h2o_ice).and.crystal_shape) then
+              zfrac= (pt(ig,l)-T_h2O_ice_clouds) / 
+     &               (T_h2O_ice_liq-T_h2O_ice_clouds)
+              zfrac= MAX(zfrac, 0.0)
+              zfrac= MIN(zfrac, 1.0)
+              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)
+              rsurf=1/(zfrac/reffh2oice(ig,l)+(1-zfrac)/rsurf)              ! radius giving the mean velocity between liquid and ice particles
+            else
+              rsurf=rfall
+            endif
+
+            vstokes(ig,l) = b / visc(ig,l) * rho * rfall**3 / rsurf *
+     &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rsurf)
+
+c Correction for high Reynolds number
+            Reynolds=2. * pplev(ig,l) / r / pt(ig,l) *
+     &        rsurf * vstokes(ig,l) / visc(ig,l)
+            if(Reynolds.ge.1.0) then
+              do i=1,5
+              Cd=24. / Reynolds * (1. + 0.15 * Reynolds**0.687) + 
+     &           0.42 / (1. + 42500 / Reynolds**1.16)                       ! (Formula from Bagheri 2018)
+              vstokes(ig,l) =(8./3.*pplev(ig,l)/r/pt(ig,l)*g*rfall**3 / 
+     &           rsurf**2/rho/Cd *
+     &           (1.+1.333*(a*pt(ig,l)/pplev(ig,l))/rsurf))**0.5
+              Reynolds=2. * pplev(ig,l) / r / pt(ig,l) * 
+     &           rsurf * vstokes(ig,l) / visc(ig,l)
+              enddo
+            endif    
 
 c           Layer crossing time (s) :