Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2598)
+++ trunk/LMDZ.MARS/README	(revision 2599)
@@ -3556,2 +3556,5 @@
 Also cleaned up and commented comsaison_h in the process.
 
+== 04/01/2021 == CM
+- Clean co2condens_mod.F
+- remove dqsurf duplication after call co2condens
Index: trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2598)
+++ trunk/LMDZ.MARS/libf/phymars/co2condens_mod.F	(revision 2599)
@@ -323,4 +323,7 @@
       ENDDO
       
+      condens_layer(:,:) = zcondicea(:,:)
+      condens_column(:) = 0.
+
       if (scavco2cond) then
         call scavenging_by_co2(ngrid,nlayer,nq,ptimestep,pplev,zq,
@@ -339,4 +342,5 @@
         DO ig=1,ngrid
           condens_column(ig) = sum(condens_layer(ig,:))
+          zfallice(ig,1) = zdqssed_co2(ig)
           piceco2(ig) = piceco2(ig) + zdqssed_co2(ig)*ptimestep
         ENDDO
@@ -347,7 +351,7 @@
      &         " ",3,pdtc) 
      
-       call WRITEdiagfi(ngrid,"zcondicea",
+       call WRITEdiagfi(ngrid,"condens_layer",
      &         "",
-     &         " ",3,zcondicea)         
+     &         " ",3,condens_layer)
       
        call WRITEdiagfi(ngrid,"zfallice",
@@ -374,5 +378,4 @@
 c
 c     No microphysics of CO2 clouds
-      if (.not.co2clouds)then
       DO ig=1,ngrid
          IF(latitude(ig).lt.0) THEN
@@ -383,5 +386,4 @@
             icap=1
          ENDIF
-       
 c
 c        Loop on where we have  condensation/ sublimation 
@@ -390,19 +392,17 @@
      $      ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  ! ground sublim.
      $      ((piceco2(ig)+zfallice(ig,1)*ptimestep) .NE. 0.))) THEN
-            condsub(ig) = .true. 
-
-            IF (zfallice(ig,1).gt.0) then  
-                 zfallheat=zfallice(ig,1)*
-     &           (pphi(ig,1)- phisfi(ig) +
-     &           cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
+            condsub(ig) = .true.
+
+            IF (zfallice(ig,1).gt.0 .and. .not. co2clouds) then
+              zfallheat = zfallice(ig,1) * (pphi(ig,1)- phisfi(ig) +
+     &                    cpice*(ztcond(ig,1)-ztcondsol(ig)))/latcond
             ELSE
-                 zfallheat=0.
+              zfallheat = 0.
             ENDIF
 c           condensation or partial sublimation of CO2 ice
 c           """""""""""""""""""""""""""""""""""""""""""""""
-            zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) 
-     &      /(latcond*ptimestep)         - zfallheat
+            zcondices(ig) = pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig))
+     &                      /(latcond*ptimestep)  - zfallheat
             pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/ptimestep
-            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
 #ifdef MESOSCALE
       print*, "not enough CO2 tracer in 1st layer to condense"
@@ -413,8 +413,16 @@
 c           """"""""""""""""""""""""""""""""""""""""""""""""""""""
             IF(ico2.ne.0) then
-c              Available CO2 tracer in layer 1 at end of timestep (kg/m2)
+c             Available CO2 tracer in layer 1 at end of timestep (kg/m2)
 #ifndef MESOSCALE
-                availco2= pq(ig,1,ico2)*((ap(1)-ap(2))+
-     &          (bp(1)-bp(2))*(pplev(ig,1)/g-zdiceco2(ig)*ptimestep))
+                availco2 = pq(ig,1,ico2)*((ap(1)-ap(2))+
+     &                     (bp(1)-bp(2))*(pplev(ig,1)/g -
+     &                     (zcondices(ig) + zfallice(ig,1))*ptimestep))
+                if ((zcondices(ig) + condens_layer(ig,1))*ptimestep
+     &           .gt.availco2) then
+                  zcondices(ig) = availco2/ptimestep -
+     &                            condens_layer(ig,1)
+                  pdtsrfc(ig) = (latcond/pcapcal(ig))*
+     &                          (zcondices(ig)+zfallheat)
+                end if
 #else
                 availco2 = pq(ig,1,igcm_co2)
@@ -422,11 +430,4 @@
      &                  ' CORRECTED FROM SIGMA LEVELS"
 #endif
-               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
 #endif
@@ -437,13 +438,20 @@
             IF((piceco2(ig)/ptimestep+zfallice(ig,1)).LE.
      &          -zcondices(ig))THEN
-               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1) 
+               zcondices(ig) = -piceco2(ig)/ptimestep - zfallice(ig,1)
                pdtsrfc(ig)=(latcond/pcapcal(ig))*
      &         (zcondices(ig)+zfallheat)
-               zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
             END IF
 
 c           Changing CO2 ice amount and pressure :
 c           """"""""""""""""""""""""""""""""""""
+            zdiceco2(ig) = zcondices(ig) + zfallice(ig,1)
+     &        + condens_column(ig)
+            if (co2clouds) then
+             ! add here only direct condensation/sublimation
+            piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
+            else
+             ! add here also CO2 ice in the atmosphere
             piceco2(ig) = piceco2(ig) + zdiceco2(ig)*ptimestep
+            end if
             pdpsrf(ig) = -zdiceco2(ig)*g
 
@@ -461,64 +469,5 @@
          END IF ! if there is condensation/sublimation
       ENDDO ! of DO ig=1,ngrid
-      else
-        do ig=1,ngrid
-          if ((ztsrf(ig)<ztcondsol(ig)) .OR. (zdqssed_co2(ig)/=0.) .OR.
-     &       ((ztsrf(ig)>ztcondsol(ig)) .AND. (piceco2(ig)/=0.))) then
-             condsub(ig) = .true.
-             zfallheat = 0.
-             !----------------------------------------------------------------------------------------------------------------------!
-             ! Compute direct condensation/sublimation of CO2 ice
-             !----------------------------------------------------------------------------------------------------------------------!
-             zcondices(ig) = pcapcal(ig) * (ztcondsol(ig)-ztsrf(ig)) /
-     &                       (latcond*ptimestep) - zfallheat
-             pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig)) / ptimestep
-             !----------------------------------------------------------------------------------------------------------------------!
-             ! If there is not enough CO2 tracer in 1st layer to condense
-             !----------------------------------------------------------------------------------------------------------------------!
-             !     Available CO2 tracer in layer 1 at end of timestep (kg/m2)
-#ifndef MESOSCALE
-             availco2 = pq(ig,1,igcm_co2) * ( (ap(1)-ap(2)) +
-     &                  (bp(1)-bp(2)) * (pplev(ig,1)/g - zcondices(ig)
-     &                  *ptimestep) )
-#else
-             availco2 = pq(ig,1,igcm_co2)
-             PRINT*, "MESOSCALE: CO2 tracer AT FIRST LEVEL IS NOT CORRECTED FROM SIGMA LEVELS"
-#endif
-             if ( zcondices(ig) * ptimestep>availco2 ) then
-               zcondices(ig) = availco2/ptimestep
-               zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig)
-               pdtsrfc(ig) = (latcond/pcapcal(ig)) *
-     &                       (zcondices(ig)+zfallheat)
-             end if
-             !----------------------------------------------------------------------------------------------------------------------!
-             ! If the entire CO2 ice layer sublimes (including what has just condensed in the atmosphere)
-             !----------------------------------------------------------------------------------------------------------------------!
-             if ( (piceco2(ig)/ptimestep) <= -zcondices(ig) ) then
-                zcondices(ig) = -piceco2(ig)/ptimestep
-                pdtsrfc(ig) = (latcond/pcapcal(ig)) *
-     &                        (zcondices(ig)+zfallheat)
-             end if
-             !----------------------------------------------------------------------------------------------------------------------!
-             ! Changing CO2 ice amount and pressure
-             !----------------------------------------------------------------------------------------------------------------------!
-             zdiceco2(ig) = zcondices(ig) + zdqssed_co2(ig) +
-     &                      condens_column(ig)
-             piceco2(ig) = piceco2(ig) + zcondices(ig)*ptimestep
-             pdpsrf(ig) = -zdiceco2(ig) * g
-
-             if (abs(pdpsrf(ig)*ptimestep)>pplev(ig,1)) then
-             print *, 'STOP in condens'
-             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)
-             call abort_physic('co2condens',
-     &                         'condensing more than total mass', 1)
-             end if
-          end if
-        end do
-      end if
+
 c  ********************************************************************
 c  Surface albedo and emissivity of the surface below the snow (emisref)
@@ -565,10 +514,9 @@
 c  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
 c  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
-            if (.not. co2clouds) then
-              zmflux(1) = -zcondices(ig)
+              zmflux(1) = -zcondices(ig) - zdqssed_co2(ig)
               DO l=1,nlayer
-                zmflux(l+1) = zmflux(l) -zcondicea(ig,l)
+                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
 #ifndef MESOSCALE
-     &          + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
+     &          + (bp(l)-bp(l+1))*(-pdpsrf(ig)/g)
 c zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
                 if (abs(zmflux(l+1)).lt.1E-13.OR.bp(l+1).eq.0.) then
@@ -576,5 +524,5 @@
                 end if
 #else
-                zmflux(l+1) = zmflux(l) - zcondicea(ig,l)
+                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
                 if (abs(zmflux(l+1)).lt.1E-13) zmflux(l+1)=0.
                 PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
@@ -582,24 +530,4 @@
 #endif
               END DO
-            else
-              zmflux(1) = - zcondices(ig) - zdqssed_co2(ig)
-              do l = 1, nlayer
-#ifndef MESOSCALE
-                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
-     &                        + (bp(l)-bp(l+1)) * (-pdpsrf(ig)/g)
-! zmflux set to 0 if very low to avoid: top layer is disappearing in v1ld
-                if (abs(zmflux(l+1))<1E-13.OR.bp(l+1)==0.) then
-                  zmflux(l+1) = 0.
-                end if
-#else
-                zmflux(l+1) = zmflux(l) - condens_layer(ig,l)
-                if (abs(zmflux(l+1))<1E-13) then
-                  zmflux(l+1) = 0.
-                end if
-                PRINT*, "MESOSCALE: FLUX THROUGH SIGMA LEVELS from"//
-     &            "dPS HAVE TO BE IMPLEMENTED"
-#endif
-              end do
-            end if
 #ifdef MESOSCALE
          print*,"absurd mass set because bp not available"
@@ -705,5 +633,5 @@
      &        ( zmflux(l)*(ztm(l) - ztc(l)) 
      &        - zmflux(l+1)*(ztm(l+1) - ztc(l))
-     &        + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )
+     &        + condens_layer(ig,l)*(ztcond(ig,l)-ztc(l))  )
                ELSE
                 zdtsig(ig,l) = (1/masse(l)) *
@@ -733,30 +661,6 @@
 #endif
 
-            if (.not. co2clouds) then
               do iq=1,nq
 !                if (noms(iq).eq.'co2') then
-                if (iq.eq.ico2) then
-c                 SPECIAL Case when the tracer IS CO2 :
-                  DO l=1,nlayer
-                    pdqc(ig,l,iq)= (1/masse(l)) *
-     &              ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
-     &              - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
-     &              + zcondicea(ig,l)*(zqc(l,iq)-1.) )
-                  END DO
-                else
-                  DO l=1,nlayer
-                    pdqc(ig,l,iq)= (1/masse(l)) *
-     &             ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
-     &             - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
-     &             + zcondicea(ig,l)*zqc(l,iq) )
-
-                    pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if watercloud=false
-                  END DO
-                end if
-              enddo
-            else
-c           Tendencies on Q
-              do iq=1,nq
-!               if (noms(iq).eq.'co2') then
                 if (iq.eq.ico2) then
 c                 SPECIAL Case when the tracer IS CO2 :
@@ -770,12 +674,12 @@
                   DO l=1,nlayer
                     pdqc(ig,l,iq)= (1/masse(l)) *
-     &              ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
-     &              - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
-     &              + condens_layer(ig,l)*zqc(l,iq) )
+     &             ( zmflux(l)*(zqm(l,iq) - zqc(l,iq))
+     &             - zmflux(l+1)*(zqm(l+1,iq) - zqc(l,iq))
+     &             + condens_layer(ig,l)*zqc(l,iq) )
+
+                    pdqc(ig,l,iq)=pdqc(ig,l,iq)+zdq_scav(ig,l,iq) ! zdq_scav=0 if watercloud=false
                   END DO
                 end if
               enddo
-            end if
-
 
        end if ! if (condsub)
@@ -785,11 +689,6 @@
 c   CO2 snow / clouds scheme
 c ***************************************************************
-      if (.not. co2clouds) then
         call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
-     &               zcondicea,zcondices,zfallice,pemisurf)
-      else
-        call co2snow(ngrid,nlayer,ptimestep,emisref,condsub,pplev,
-     &               zcondicea,zcondices,zdqssed_co2,pemisurf)
-      end if
+     &               condens_layer,zcondices,zfallice,pemisurf)
 c ***************************************************************
 c Ecriture des diagnostiques
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2598)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2599)
@@ -874,4 +874,14 @@
           co2totA = co2totA + co2ice(ig)
         end do
+      else
+        co2totA = 0.
+        do ig=1,ngrid
+          do l=1,nlayer
+             co2totA = co2totA + (zplev(ig,l)-zplev(ig,l+1))/g*
+     &             (pq(ig,l,igcm_co2)
+     &        +pdq(ig,l,igcm_co2)*ptimestep)
+          end do
+          co2totA = co2totA + co2ice(ig)
+        end do
       endif ! of if (igcm_co2_ice.ne.0)
 c-----------------------------------------------------------------------
@@ -2059,11 +2069,4 @@
             ENDDO
            ENDDO
-
-          DO iq=1, nq
-            DO ig=1,ngrid
-              dqsurf(ig,iq)=dqsurf(ig,iq)+zdqsc(ig,iq)
-            ENDDO  ! (ig)
-           ENDDO    ! (iq)
-
          ENDIF ! of IF (tracer)
 
@@ -3953,9 +3956,17 @@
           co2totB = co2totB + co2ice(ig)
         enddo
-
+      else
+        co2totB = 0. ! added by C.M.
+        do ig=1,ngrid
+          do l=1,nlayer
+            co2totB = co2totB + (zplev(ig,l)-zplev(ig,l+1))/g*
+     &             (pq(ig,l,igcm_co2)+pdq(ig,l,igcm_co2)*ptimestep)
+          enddo
+          co2totB = co2totB + co2ice(ig)
+        enddo
+      endif ! of if (igcm_co2_ice.ne.0)
         call WRITEDIAGFI(ngrid,'co2conservation',
      &                     'Total CO2 mass conservation in physic',
      &                     '%',0,[(co2totA-co2totB)/co2totA])
-      endif ! of if (igcm_co2_ice.ne.0)
 ! XIOS outputs
 #ifdef CPP_XIOS      
