Changeset 3157

Timestamp:

Dec 12, 2023, 10:45:40 AM (14 months ago)

Author:

csegonne

Message:

MARS PCM

• zzlay and zzlev are calculated at the begining of physics taking into account the evolution of gravitational acceleration with altitude (g becomes gz) and varying reduced gas constant with composition of the atmosphere (r becomes rnew).
• zzlay and zzlev are updated at the end of physics after call co2condens with updated pressure and temperature.
• call concentrations, when photochem or callthermos is true, has been moved before the first calculation of zzlay and zzlev to be able to use varying reduced gas constant rnew.

Location:

trunk/LMDZ.MARS

Files:

• changelog.txt (modified) (1 diff)
• libf/phymars/physiq_mod.F (modified) (9 diffs)

trunk/LMDZ.MARS/changelog.txt

@@ -4399 +4399 @@
 == 12/12/2023 == CS
 rnew and cpnew are already initialized in physiq_mod.F with call concentrations.F if callthermos is true.
+
+== 12/12/2023 == CS
+* zzlay and zzlev are calculated at the begining of physics taking into account the evolution of gravitational acceleration with altitude (g becomes gz) and varying reduced gas constant with composition of the atmosphere (r becomes rnew).
+* zzlay and zzlev are updated at the end of physics after call co2condens with updated pressure and temperature.
+* call concentrations, when photochem or callthermos is true, has been moved before the first calculation of zzlay and zzlev to be able to use varying reduced gas constant rnew.

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -343 +343 @@
       REAL zzlay(ngrid,nlayer)     ! altitude at the middle of the layers
       REAL zzlev(ngrid,nlayer+1)   ! altitude at layer boundaries
+      REAL gz(ngrid,nlayer)        ! variation of g with altitude from aeroid surface
 !      REAL latvl1,lonvl1             ! Viking Lander 1 point (for diagnostic)
 
@@ -393 +394 @@
       INTEGER length
       PARAMETER (length=100)
+      REAL tlaymean ! temporary value of mean layer temperature for zzlay
+
 
 c     Variables for the total dust for diagnostics
@@ -836 +839 @@
       dwatercap(:,:)=0
 
+
+
       call compute_meshgridavg(ngrid,nq,albedo,emis,tsurf,qsurf,
      &       albedo_meshavg,emis_meshavg,tsurf_meshavg,qsurf_meshavg)
@@ -869 +874 @@
       ps(:) = pplev(:,1)
 
+
+#ifndef MESOSCALE
+c-----------------------------------------------------------------------
+c    1.2.1 Compute mean mass, cp, and R
+c    concentrations outputs rnew(ngrid,nlayer), cpnew(ngrid,nlayer)
+c       , mmean(ngrid,nlayer) and Akknew(ngrid,nlayer)
+c    --------------------------------
+
+      if(photochem.or.callthermos) then
+         call concentrations(ngrid,nlayer,nq,
+     &                       zplay,pt,pdt,pq,pdq,ptimestep)
+      endif
+#endif
+
 c     Compute geopotential at interlayers
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 c     ponderation des altitudes au niveau des couches en dp/p
-
-      DO l=1,nlayer
-         DO ig=1,ngrid
-            zzlay(ig,l)=pphi(ig,l)/g
-         ENDDO
+cc ------------------------------------------
+        !Calculation zzlev & zzlay for first layer      
+      DO ig=1,ngrid
+         zzlay(ig,1)=-(log(pplay(ig,1)/ps(ig)))*rnew(ig,1)*pt(ig,1)/g                                      
+         zzlev(ig,1)=0
+         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
+         gz(ig,1)=g
+        
+        DO l=2,nlayer     
+         ! compute "mean" temperature of the layer
+          if(pt(ig,l) .eq. pt(ig,l-1)) then
+             tlaymean=pt(ig,l)
+          else
+             tlaymean=(pt(ig,l)- pt(ig,l-1))/log(pt(ig,l)/pt(ig,l-1))
+          endif
+                     
+          ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
+          gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
+          
+          zzlay(ig,l)=zzlay(ig,l-1)- 
+     &   (log(pplay(ig,l)/pplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
+                 
+          z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
+          z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
+          zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+        ENDDO
       ENDDO
-      DO ig=1,ngrid
-         zzlev(ig,1)=0.
-         zzlev(ig,nlayer+1)=1.e7    ! dummy top of last layer above 10000 km...
-      ENDDO
-      DO l=2,nlayer
-         DO ig=1,ngrid
-            z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
-            z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
-            zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
-         ENDDO
-      ENDDO
-
 
 !     Potential temperature calculation not the same in physiq and dynamic
@@ -902 +930 @@
       ENDDO
 
-#ifndef MESOSCALE
-c-----------------------------------------------------------------------
-c    1.2.5 Compute mean mass, cp, and R
-c    --------------------------------
-
-      if(photochem.or.callthermos) then
-         call concentrations(ngrid,nlayer,nq,
-     &                       zplay,pt,pdt,pq,pdq,ptimestep)
-      endif
-#endif
 
       ! Compute vertical velocity (m/s) from vertical mass flux
@@ -1349 +1367 @@
      &                      pdqrds,wspeed,dsodust,dsords,dsotop,
      &                      tau_pref_scenario,tau_pref_gcm)
-
+     
 c      update the tendencies of both dust after vertical transport
          DO l=1,nlayer
@@ -1609 +1627 @@
      $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th,
      $ dtke_th,zdhdif,hfmax_th,wstar,sensibFlux)
-
+       
          DO l=1,nlayer
            DO ig=1,ngrid
@@ -2061 +2079 @@
      &                tau,tauscaling)
 
+
 c Flux at the surface of co2 ice computed in co2cloud microtimestep
            IF (rdstorm) THEN
@@ -2282 +2301 @@
         ENDDO
         zplev(:,nlayer+1) = 0.
-        ! update layers altitude
-        DO l=2,nlayer
-          DO ig=1,ngrid
-           z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
-           z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
-           zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
+        
+        ! Calculation of zzlay and zzlay with udpated pressure and temperature
+        DO ig=1,ngrid
+         zzlay(ig,1)=-(log(zplay(ig,1)/ps(ig)))*rnew(ig,1)*
+     &          (pt(ig,1)+pdt(ig,1)*ptimestep) /g                                      
+
+          DO l=2,nlayer
+                 
+           ! compute "mean" temperature of the layer
+            if((pt(ig,l)+pdt(ig,l)*ptimestep) .eq.
+     &               (pt(ig,l-1)+pdt(ig,l-1)*ptimestep)) then
+               tlaymean= pt(ig,l)+pdt(ig,l)*ptimestep
+            else
+               tlaymean=((pt(ig,l)+pdt(ig,l)*ptimestep)- 
+     &            (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))/
+     &            log((pt(ig,l)+pdt(ig,l)*ptimestep)/ 
+     &            (pt(ig,l-1)+pdt(ig,l-1)*ptimestep))
+            endif
+                     
+            ! compute gravitational acceleration (at altitude zaeroid(nlayer-1))
+            gz(ig,l)=g*(rad**2)/(rad+zzlay(ig,l-1)+(phisfi(ig)/g))**2
+          
+
+            zzlay(ig,l)=zzlay(ig,l-1)- 
+     &     (log(zplay(ig,l)/zplay(ig,l-1))*rnew(ig,l)*tlaymean/gz(ig,l))
+              
+        
+        !   update layers altitude
+             z1=(zplay(ig,l-1)+zplev(ig,l))/(zplay(ig,l-1)-zplev(ig,l))
+             z2=(zplev(ig,l)+zplay(ig,l))/(zplev(ig,l)-zplay(ig,l))
+             zzlev(ig,l)=(z1*zzlay(ig,l-1)+z2*zzlay(ig,l))/(z1+z2)
           ENDDO
-        ENDDO
+        ENDDO 
 #endif
       ENDIF  ! of IF (callcond)