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

Timestamp:

Dec 21, 2011, 10:09:07 AM (13 years ago)

Author:

rwordsworth

Message:

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:

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

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

@@ -79 +79 @@
             PRINT*,'ngrid  =',ngrid
             PRINT*,'ngridmx  =',ngridmx
-            STOP
+            STOP 
          ENDIF
          firstcall=.false.
@@ -111 +111 @@
 c     (stokes law corrected for low pressure by the Cunningham
 c     slip-flow correction  according to Rossow (Icarus 36, 1-50, 1978)
-
+       
         do  l=1,nlay
           do ig=1, ngrid
@@ -121 +121 @@
             endif  
 
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! TEMPORARY MODIF BY RDW !!!!
-            !rfall=5.e-6
-            if(rfall.lt.1.e-7)then
-               rfall=1.e-7
-            endif
-            !if(rfall.gt.5.e-5)then
-            !   rfall=5.e-5
-            !endif
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
             vstokes(ig,l) = b * rho * rfall**2 *
      &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
@@ -147 +135 @@
 c      va traverser le niveau intercouche l : "dztop" est sa hauteur
 c      au dessus de l (m), "ptop" est sa pression (Pa))
-
       do  l=1,nlay
         do ig=1, ngrid
@@ -154 +141 @@
              Ep=0
              k=0
-
+           w(ig,l) = 0. !! JF+AS ajout initialisation (LK MARS)
 c **************************************************************
 c            Simple Method
-             w(ig,l) =
-     &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
-cc           write(*,*) 'OK simple method l,w =', l, w(ig,l)
-cc           write(*,*) 'OK simple method dztop =', dztop
-c **************************************************************
+cc             w(ig,l) =
+cc     &       (1- exp(-dztop*g/(r*pt(ig,l))))*pplev(ig,l) / g
+cc             write(*,*) 'OK simple method l,w =', l, w(ig,l)
+cc            write(*,*) 'OK simple method dztop =', dztop
+           w(ig,l) = 1. - exp(-dztop*g/(r*pt(ig,l)))
+             !!! Diagnostic: JF. Fix: AS. Date: 05/11
+             !!! Probleme arrondi avec la quantite ci-dessus
+             !!! ---> vaut 0 pour -dztop*g/(r*pt(ig,l)) trop petit
+             !!! ---> dans ce cas on utilise le developpement limite !
+             !!! ---> exp(-x) = 1 - x lorsque x --> 0 avec une erreur de x^2 / 2  
+
+             IF ( w(ig,l) .eq. 0. ) THEN
+                w(ig,l) = ( dztop*g/(r*pt(ig,l)) ) * pplev(ig,l) / g
+             ELSE
+                w(ig,l) = w(ig,l) * pplev(ig,l) / g
+             ENDIF
+! LK borrowed simple method from Mars model (AS/JF)
+
+!**************************************************************
 cccc         Complex method :
-            if (dztop.gt.epaisseur(ig,l)) then
+           if (dztop.gt.epaisseur(ig,l)) then
 cccc            Cas ou on "epuise" la couche l : On calcule le flux
 cccc            Venant de dessus en tenant compte de la variation de Vstokes
@@ -176 +177 @@
                enddo 
                Ep = Ep - epaisseur(ig,l+k)
-             end if
-             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
-             w(ig,l) = (pplev(ig,l) -Ptop)/g
+!             ptop=pplev(ig,l+k)*exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
+             ptop=exp(-(dztop-Ep)*g/(r*pt(ig,l+k)))
+             IF ( ptop .eq. 1. ) THEN
+                !PRINT*, 'newsedim: exposant trop petit ', ig, l
+                ptop=pplev(ig,l+k) * ( 1. - (dztop-Ep)*g/(r*pt(ig,l+k)))
+             ELSE
+                ptop=pplev(ig,l+k) * ptop
+             ENDIF
+
+             w(ig,l) = (pplev(ig,l) - ptop)/g
+
+            endif   !!! complex method
 c
 cc           write(*,*) 'OK new    method l,w =', l, w(ig,l)
@@ -188 +198 @@
 cc            if(l.eq.8)write(*,*)'l=8,k, w',k, w(1,l)
 c **************************************************************
+
         end do
       end do
@@ -195 +206 @@
 c    &                wq(1,6),wq(1,7),pqi(1,6)
 
-
       RETURN
       END
-