Changeset 2124 for trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F

Timestamp:

Apr 24, 2019, 12:11:24 PM (6 years ago)

Author:

emillour

Message:

Mars GCM:
Updated co2condens to correctly conserve tracer mass
EM+FF

File:

• trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F (modified) (14 diffs)

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

@@ -16 +16 @@
        use tracer_mod, only: noms
        use surfdat_h, only: emissiv, phisfi
-       use geometry_mod, only: latitude ! grid point latitudes (rad)
+       use geometry_mod, only: latitude, ! grid point latitudes (rad)
+     &                         longitude_deg, latitude_deg
        use planete_h, only: obliquit
        use comcstfi_h, only: cpp, g, r, pi 
 #ifndef MESOSCALE
-       USE vertical_layers_mod, ONLY: bp
+       USE vertical_layers_mod, ONLY: ap, bp
 #endif
        IMPLICIT NONE
@@ -113 +114 @@
       REAL masse(nlayer),w(nlayer+1)
       REAL Sm(nlayer),Smq(nlayer,nq),mixmas,qmix
+      REAL availco2
       LOGICAL condsub(ngrid)
 
@@ -130 +132 @@
 c   local saved variables
       integer,save :: ico2 ! index of CO2 tracer
-      real,save :: qco2min,qco2,mmean
+      real,save :: qco2,mmean
       real,parameter :: latcond=5.9e5 ! (J/kg) Latent heat of solid CO2 ice
       real,parameter :: tcond1mb=136.27 ! condensation temperature (K) at 1 mbar
@@ -139 +141 @@
       LOGICAL,SAVE :: firstcall = .true. !,firstcall2=.true.
 
-c D.BARDET: to debug 
+c D.BARDET: for debug 
       real ztc3D(ngrid,nlayer)
       REAL ztm3D(ngrid,nlayer)
@@ -156 +158 @@
 
          firstcall=.false.
-         write(*,*) 'Newcondens: improved_ztcond=',improved_ztcond
-         PRINT*,'In newcondens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
+         write(*,*) 'CO2condens: improved_ztcond=',improved_ztcond
+         PRINT*,'In co2condens:Tcond(P=1mb)=',tcond1mb,' Lcond=',latcond
          PRINT*,'acond,bcond,ccond',acond,bcond,ccond
 
@@ -177 +179 @@
              endif
            enddo
-c          minimum CO2 mix. ratio below which mixing occur with layer above: 
-           qco2min =0.75  
          end if
       ENDIF ! of IF (firstcall)
@@ -378 +378 @@
      &      /(latcond*ptimestep)         - zfallheat
             pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-       
+            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
+
+c           If there is not enough CO2 tracer in 1st layer to condense
+c           """"""""""""""""""""""""""""""""""""""""""""""""""""""
+            IF(ico2.ne.0) then
+c              Available CO2 tracer in layer 1 at end of timestep (kg/m2)
+                availco2= pq(ig,1,ico2)*((ap(1)-ap(2))+
+     &          (bp(1)-bp(2))*(pplev(ig,1)/g-zdiceco2(ig)*ptimestep))
+
+               IF ((zcondices(ig) + zcondicea(ig,1))*ptimestep
+     &           .gt.availco2) then
+                   zcondices(ig) = availco2/ptimestep -zcondicea(ig,1)
+                   zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
+                   pdtsrfc(ig)=(latcond/pcapcal(ig))*
+     &                          (zcondices(ig)+zfallheat)
+               ENDIF 
+            ENDIF 
+
+
 c           If the entire CO_2 ice layer sublimes
 c           """"""""""""""""""""""""""""""""""""""""""""""""""""
@@ -388 +406 @@
                pdtsrfc(ig)=(latcond/pcapcal(ig))*
      &         (zcondices(ig)+zfallheat)
+               zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
             END IF
 
@@ -393 +412 @@
 c           """"""""""""""""""""""""""""""""""""
 
-            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
             piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
             pdpsrf(ig) = -zdiceco2(ig)*g
@@ -401 +419 @@
                PRINT*,'condensing more than total mass'
                PRINT*,'Grid point ',ig
+               PRINT*,'Longitude(degrees): ',longitude_deg(ig)
+               PRINT*,'Latitude(degrees): ',latitude_deg(ig)
                PRINT*,'Ps = ',pplev(ig,1)
                PRINT*,'d Ps = ',pdpsrf(ig)
@@ -416 +436 @@
            if(.not.piceco2(ig).ge.0.) THEN
               if(piceco2(ig).le.-5.e-8) print*,
-     $         'WARNING newcondens piceco2(',ig,')=', piceco2(ig)
+     $         'WARNING co2condens piceco2(',ig,')=', piceco2(ig)
                piceco2(ig)=0.
            endif
@@ -465 +485 @@
             END DO
 
-c Mass of each layer
-c ------------------ 
+c Mass of each layer at the end of timestep
+c ----------------------------------------- 
             DO l=1,nlayer
-              masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
+              masse(l)=( pplev(ig,l) - pplev(ig,l+1) +  
+     &                 (bp(l)-bp(l+1))*pdpsrf(ig)*ptimestep)/g
             END DO
 
@@ -628 +649 @@
            ! NB: Updated surface pressure, at grid point 'ngrid', is
            !     ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep
-!             write(*,*) "newcondens: South pole: latitude(ngrid)=",
+!             write(*,*) "co2condens: South pole: latitude(ngrid)=",
 !     &                                           latitude(ngrid)
              ztcondsol(ngrid)=
@@ -754 +775 @@
        enddo
 
-c     boundary conditions (not used in newcondens !!)
+c     boundary conditions (not used in co2condens !!)
 c         qm(nlayer+1)=0.
 c         if(w(1).gt.0.) then