Changeset 2091 for trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

Timestamp:

Feb 6, 2019, 2:06:01 PM (6 years ago)

Author:

mvals

Message:

Mars GCM:

• follow-up of the last change in rocketduststorm_mod.F90: corrections in vector definitions

MV

File:

• trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90 (modified) (4 diffs)

trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90

@@ -92 +92 @@
       REAL zdtsw1(ngrid,nlayer)    ! (K/s) storm
       REAL zt(ngrid,nlayer)        ! actual temperature at mid-layer (K)
-      REAL zdtvert(nlayer)   ! dT/dz , lapse rate
-      REAL ztlev(nlayer)     ! temperature at intermediate levels l+1/2 without last level
+      REAL zdtvert(ngrid,nlayer)   ! dT/dz , lapse rate
+      REAL ztlev(ngrid,nlayer)     ! temperature at intermediate levels l+1/2 without last level
 
       REAL zdtlw1_lev(nlayer),zdtsw1_lev(nlayer) ! rad. heating rate at intermediate levels l+1/2 for stormdust
@@ -193 +193 @@
       wqmass(:,:)=0.
       wqnumber(:,:)=0.
+      zdtvert(:,:)=0.
       lapserate(:,:)=0.
       deltahr(:,:)=0.
@@ -251 +252 @@
         !! gradient at boundaries of each layer, start from surface
         DO ig=1,ngrid
-          zdtvert(1)=0. !This is the env. lapse rate 
-          DO l=1,nlayer-1
-            zdtvert(l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
-            lapserate(ig,l+1)=zdtvert(l+1) !for diagnosing
-          ENDDO
-      
-          ! computing heating rates gradient at boundraies of each layer
-          ! start from surface
-          zdtlw1_lev(1)=0.
-          zdtsw1_lev(1)=0.
-          zdtlw_lev(1)=0.
-          zdtsw_lev(1)=0.
-          ztlev(1)=zt(ig,1)
-
-          DO l=1,nlayer-1
-            ! Calculation for the dust storm fraction
-            zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
-                         zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
-                            (pzlay(ig,l+1)-pzlay(ig,l))
+          IF (storm(ig)) THEN
+
+            scheme(ig)=2
+            !! This is the env. lapse rate 
+            zdtvert(ig,1)=0.
+            DO l=1,nlayer-1
+              zdtvert(ig,l+1)=(zt(ig,l+1)-zt(ig,l))/(pzlay(ig,l+1)-pzlay(ig,l))
+              lapserate(ig,l+1)=zdtvert(ig,l+1)
+            ENDDO
+      
+            !! computing heating rates gradient at boundraies of each layer
+            !! start from surface
+            zdtlw1_lev(1)=0.
+            zdtsw1_lev(1)=0.
+            zdtlw_lev(1)=0.
+            zdtsw_lev(1)=0.
+            ztlev(ig,1)=zt(ig,1)
+
+            DO l=1,nlayer-1
+              !! Calculation for the dust storm fraction
+              zdtlw1_lev(l+1)=(zdtlw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
+                           zdtlw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
+                              (pzlay(ig,l+1)-pzlay(ig,l))
+            
+              zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
+                           zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
+                              (pzlay(ig,l+1)-pzlay(ig,l))
+              !! Calculation for the background dust fraction
+              zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
+                           pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
+                              (pzlay(ig,l+1)-pzlay(ig,l))
            
-            zdtsw1_lev(l+1)=(zdtsw1(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+ &
-                         zdtsw1(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /  &
-                            (pzlay(ig,l+1)-pzlay(ig,l))
-            !! Calculation for the background dust fraction
-            zdtlw_lev(l+1)=(pdtlw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
-                         pdtlw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
-                            (pzlay(ig,l+1)-pzlay(ig,l))
+              zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
+                           pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
+                              (pzlay(ig,l+1)-pzlay(ig,l))
            
-            zdtsw_lev(l+1)=(pdtsw(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+   &
-                         pdtsw(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /   &
-                            (pzlay(ig,l+1)-pzlay(ig,l))
-           
-            ztlev(l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
-                         zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
-                            (pzlay(ig,l+1)-pzlay(ig,l))
-          ENDDO ! DO l=1,nlayer-1
-        ENDDO ! DO ig=1,ngrid     
-
-        !! **********************************************************************
-        !! 2.3 Calculation of the vertical velocity : extra heating
-        !! balanced by adiabatic cooling 
-        DO ig=1,ngrid
-          IF (storm(ig)) THEN       
-            scheme(ig)=2         
+              ztlev(ig,l+1)=(zt(ig,l)*(pzlay(ig,l+1)-pzlev(ig,l+1))+          &
+                           zt(ig,l+1)*(pzlev(ig,l+1)-pzlay(ig,l)))  /      &
+                              (pzlay(ig,l+1)-pzlay(ig,l))
+            ENDDO ! DO l=1,nlayer-1    
+
+            !! **********************************************************************
+            !! 2.3 Calculation of the vertical velocity : extra heating
+            !! balanced by adiabatic cooling 
+            
             DO l=1,nlayer
               deltahr(ig,l)=(zdtlw1_lev(l)+zdtsw1_lev(l))  &
                                           -(zdtlw_lev(l)+zdtsw_lev(l))
               wrad(ig,l)=-deltahr(ig,l)/(g/cpp+   &
-                                         max(zdtvert(l),-0.99*g/cpp))       
+                                         max(zdtvert(ig,l),-0.99*g/cpp))       
               !! Limit vertical wind in case of lapse rate close to adiabatic
               wrad(ig,l)=max(wrad(ig,l),-wmax) 
               wrad(ig,l)=min(wrad(ig,l),wmax)
             ENDDO
+
           ENDIF ! IF (storm(ig))
         ENDDO ! DO ig=1,ngrid
@@ -337 +340 @@
             !! Mass flux in kg/m2
             DO l=1,nlayer
-               w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(l)))*ptimestep
+               w(ig,l)=wrad(ig,l)*(pplev(ig,l)/(r*ztlev(ig,l)))*ptimestep
             ENDDO