source: trunk/LMDZ.PLUTO/libf/phypluto/condense_n2.F90 @ 3347

subroutine condense_n2(klon,klev,nq,ptimestep, &
      pcapcal,pplay,pplev,ptsrf,pt, &
      pphi,pdt,pdu,pdv,pdtsrf,pu,pv,pq,pdq, &
      picen2,psolaralb,pemisurf, &
      pdtc,pdtsrfc,pdpsrf,pduc,pdvc, &
      pdqc,pdicen2)

  use radinc_h, only : naerkind
  use comgeomfi_h
  use comcstfi_mod, only: g, r, cpp, pi
  USE surfdat_h, only: phisfi,kp,p00
  USE tracer_h, only: noms, igcm_n2, lw_n2
  USE callkeys_mod, only: fast,ch4lag,latlag,nbsub,no_n2frost,tsurfmax,kmixmin,source_haze,vmrlag
  USE comvert_mod, ONLY: ap,bp
  use geometry_mod, only: latitude


  implicit none

!==================================================================
!     Purpose
!     -------
!     Condense and/or sublime N2 ice on the ground and in the
!     atmosphere, and sediment the ice.
!
!     Inputs
!     ------
!     klon                 Number of vertical columns
!     klev                Number of layers
!     pplay(klon,klev)   Pressure layers
!     pplev(klon,klev+1) Pressure levels
!     pt(klon,klev)      Temperature (in K)
!     ptsrf(klon)          Surface temperature
!
!     pdt(klon,klev)   Time derivative before condensation/sublimation of pt
!     pdtsrf(klon)         Time derivative before condensation/sublimation of ptsrf
!     picen2(klon)         n2 ice at the surface (kg/m2)
!
!     Outputs
!     -------
!     pdpsrf(klon)         \  Contribution of condensation/sublimation
!     pdtc(klon,klev)  /  to the time derivatives of Ps, pt, and ptsrf
!     pdtsrfc(klon)       /
!     pdicen2(klon)         Tendancy n2 ice at the surface (kg/m2)
!
!     Both
!     ----
!     psolaralb(klon)      Albedo at the surface
!     pemisurf(klon)       Emissivity of the surface
!
!     Authors
!     -------
!     Francois Forget (1996,2013)
!     Converted to Fortran 90 and slightly modified by R. Wordsworth (2009)
!
!
!==================================================================

!-----------------------------------------------------------------------
!     Arguments

  INTEGER klon, klev, nq

  REAL ptimestep
  REAL pplay(klon,klev),pplev(klon,klev+1)
  REAL pcapcal(klon)
  REAL pt(klon,klev)
  REAL ptsrf(klon),flu1(klon),flu2(klon),flu3(klon)
  REAL pphi(klon,klev)
  REAL pdt(klon,klev),pdtsrf(klon),pdtc(klon,klev)
  REAL pdtsrfc(klon),pdpsrf(klon)
  REAL picen2(klon),psolaralb(klon),pemisurf(klon)


  REAL pu(klon,klev) ,  pv(klon,klev)
  REAL pdu(klon,klev) , pdv(klon,klev)
  REAL pduc(klon,klev) , pdvc(klon,klev)
  REAL pq(klon,klev,nq),pdq(klon,klev,nq)
  REAL pdqc(klon,klev,nq)

!-----------------------------------------------------------------------
!     Local variables

  INTEGER l,ig,ilay,it,iq,ich4_gas

  REAL*8 zt(klon,klev)
  REAL tcond_n2
  REAL pcond_n2
  REAL glob_average2d         ! function
  REAL zqn2(klon,klev) ! N2 MMR used to compute Tcond/zqn2
  REAL ztcond (klon,klev)
  REAL ztcondsol(klon),zfallheat
  REAL pdicen2(klon)
  REAL zcondicea(klon,klev), zcondices(klon)
  REAL zfallice(klon,klev+1)
  REAL zmflux(klev+1)
  REAL zu(klev),zv(klev)
  REAL zq(klev,nq),zq1(klev)
  REAL ztsrf(klon)
  REAL ztc(klev), ztm(klev+1)
  REAL zum(klev+1) , zvm(klev+1)
  REAL zqm(klev+1,nq),zqm1(klev+1)
  LOGICAL condsub(klon)
  REAL subptimestep
  Integer Ntime
  real masse (klev),w(klev+1)
  real wq(klon,klev+1)
  real vstokes,reff
  real dWtotsn2
  real condnconsn2(klon)
  real nconsMAXn2
!     Special diagnostic variables
  real tconda1(klon,klev)
  real tconda2(klon,klev)
  real zdtsig (klon,klev)
  real zdtlatent (klon,klev)
  real zdt (klon,klev)
  REAL albediceF(klon)
!   SAVE albediceF
  INTEGER nsubtimestep,itsub    !number of subtimestep when calling vl1d
  REAL subtimestep        !ptimestep/nsubtimestep
  REAL dtmax

  REAL zplevhist(klon)
  REAL zplevnew(klon)
  REAL zplev(klon)
  REAL zpicen2(klon)
  REAL ztsrfhist(klon)
  REAL zdtsrf(klon)
  REAL globzplevnew

!   REAL vmrn2(klon)
!   SAVE vmrn2
  REAL stephan
  DATA stephan/5.67e-08/  ! Stephan Boltzman constant
  SAVE stephan
!-----------------------------------------------------------------------
!     Saved local variables

!      REAL latcond
  REAL ccond
  REAL cpice  ! for atm condensation
  SAVE cpice
!      SAVE latcond,ccond
  SAVE ccond

  LOGICAL firstcall
  SAVE firstcall
  REAL SSUM
  EXTERNAL SSUM

!      DATA latcond /2.5e5/
!      DATA latcond /1.98e5/
  DATA cpice /1300./
  DATA firstcall/.true./

  INTEGER,SAVE :: i_n2ice=0       ! n2 ice
  CHARACTER(LEN=20) :: tracername ! to temporarily store text
  logical olkin  ! option to prevent N2 ice effect in the south
  DATA olkin/.false./
  save olkin

  CHARACTER(len=10) :: tname

!-----------------------------------------------------------------------

!     Initialisation
  IF (firstcall) THEN
     ccond=cpp/(g*lw_n2)
     print*,'In condense_n2cloud: ccond=',ccond,' latcond=',lw_n2

     ! calculate global mean surface pressure for the fast mode
     IF (fast) THEN
        IF (.not. ALLOCATED(kp)) ALLOCATE(kp(klon))
        DO ig=1,klon
           kp(ig) = exp(-phisfi(ig)/(r*38.))
        ENDDO
        p00=glob_average2d(kp) ! mean pres at ref level
     ENDIF

     !vmrn2(:) = 1.
     !IF (ch4lag) then
     !   DO ig=1,klon
     !      if (latitude(ig)*180./pi.ge.latlag) then
     !         vmrn2(ig) = vmrlag
     !      endif
     !   ENDDO
     !ENDIF
     !IF (no_n2frost) then
     !   DO ig=1,klon
     !      if (picen2(ig).eq.0.) then
     !         vmrn2(ig) = 1.e-15
     !      endif
     !   ENDDO
     !ENDIF
     firstcall=.false.
  ENDIF

!-----------------------------------------------------------------------
!     Calculation of n2 condensation / sublimation

!     Variables used:
!     picen2(klon)            : amount of n2 ice on the ground     (kg/m2)
!     zcondicea(klon,klev): condensation rate in layer l       (kg/m2/s)
!     zcondices(klon)         : condensation rate on the ground    (kg/m2/s)
!     zfallice(klon,klev) : amount of ice falling from layer l (kg/m2/s)
!     zdtlatent(klon,klev): dT/dt due to phase changes         (K/s)

!     Tendencies initially set to 0
  zcondices(1:klon) = 0.
  pdtsrfc(1:klon) = 0.
  pdpsrf(1:klon) = 0.
  ztsrfhist(1:klon) = 0.
  condsub(1:klon) = .false.
  pdicen2(1:klon) = 0.
  zfallheat=0
  pdqc(1:klon,1:klev,1:nq)=0
  pdtc(1:klon,1:klev)=0
  pduc(1:klon,1:klev)=0
  pdvc(1:klon,1:klev)=0
  zfallice(1:klon,1:klev+1)=0
  zcondicea(1:klon,1:klev)=0
  zdtlatent(1:klon,1:klev)=0
  zt(1:klon,1:klev)=0.

!-----------------------------------------------------------------------
!     Atmospheric condensation

!     Condensation / sublimation in the atmosphere
!     --------------------------------------------
!      (calcul of zcondicea , zfallice and pdtc)

  zt(1:klon,1:klev)=pt(1:klon,1:klev)+ pdt(1:klon,1:klev)*ptimestep
  if (igcm_n2.ne.0) then
     zqn2(1:klon,1:klev) = 1. ! & temporaire
!        zqn2(1:klon,1:klev)= pq(1:klon,1:klev,igcm_n2) + pdq(1:klon,1:klev,igcm_n2)*ptimestep
  else
     zqn2(1:klon,1:klev) = 1.
  end if

  if (.not.fast) then
!     Forecast the atmospheric frost temperature 'ztcond' with function tcond_n2
    DO l=1,klev
      DO ig=1,klon
          ztcond (ig,l) = tcond_n2(pplay(ig,l),zqn2(ig,l))
      ENDDO
    ENDDO

    DO l=klev,1,-1
      DO ig=1,klon
       pdtc(ig,l)=0.  ! final tendancy temperature set to 0

       IF((zt(ig,l).LT.ztcond(ig,l)).or.(zfallice(ig,l+1).gt.0))THEN
           condsub(ig)=.true.  !condensation at level l
           IF (zfallice(ig,l+1).gt.0) then
             zfallheat=zfallice(ig,l+1)*&
            (pphi(ig,l+1)-pphi(ig,l) +&
           cpice*(ztcond(ig,l+1)-ztcond(ig,l)))/lw_n2
           ELSE
               zfallheat=0.
           ENDIF
           zdtlatent(ig,l)=(ztcond(ig,l) - zt(ig,l))/ptimestep
           zcondicea(ig,l)=(pplev(ig,l)-pplev(ig,l+1))&
                         *ccond*zdtlatent(ig,l)- zfallheat
!              Case when the ice from above sublimes entirely
!              """""""""""""""""""""""""""""""""""""""""""""""
           IF ((zfallice(ig,l+1).lt.-zcondicea(ig,l)) &
                .AND. (zfallice(ig,l+1).gt.0)) THEN

              zdtlatent(ig,l)=(-zfallice(ig,l+1)+zfallheat)/&
               (ccond*(pplev(ig,l)-pplev(ig,l+1)))
              zcondicea(ig,l)= -zfallice(ig,l+1)
           END IF

           zfallice(ig,l) = zcondicea(ig,l)+zfallice(ig,l+1)

        END IF

      ENDDO
    ENDDO
  endif  ! not fast

!-----------------------------------------------------------------------
!     Condensation/sublimation on the ground
!     (calculation of zcondices and pdtsrfc)

!     Adding subtimesteps : in fast version, pressures too low lead to
!     instabilities.
  IF (fast) THEN
     IF (pplev(1,1).gt.0.3) THEN
         nsubtimestep= 1
     ELSE
         nsubtimestep= nbsub !max(nint(ptimestep/dtmax),1)
     ENDIF
  ELSE
     nsubtimestep= 1 ! if more then kp and p00 have to be calculated
                     ! for nofast mode
  ENDIF
  subtimestep=ptimestep/float(nsubtimestep)

  do itsub=1,nsubtimestep
     ! first loop : getting zplev, ztsurf
     IF (itsub.eq.1) then
       DO ig=1,klon
          zplev(ig)=pplev(ig,1)
          ztsrfhist(ig)=ptsrf(ig) + pdtsrf(ig)*ptimestep
          ztsrf(ig)=ptsrf(ig) + pdtsrf(ig)*subtimestep    !!
          zpicen2(ig)=picen2(ig)
       ENDDO
     ELSE
     ! additional loop :
     ! 1)  pressure updated
     ! 2)  direct redistribution of pressure over the globe
     ! 3)  modification pressure for unstable cases
     ! 4)  pressure update to conserve mass
     ! 5)  temperature updated with radiative tendancies
       DO ig=1,klon
          zplevhist(ig)=zplev(ig)
          zplevnew(ig)=zplev(ig)+pdpsrf(ig)*subtimestep   ! 1)
          !IF (zplevnew(ig).le.0.0001) then
          !   zplevnew(ig)=0.0001*kp(ig)/p00
          !ENDIF
       ENDDO
       ! intermediaire de calcul: valeur moyenne de zplevnew (called twice in the code)
       globzplevnew=glob_average2d(zplevnew)
       DO ig=1,klon
         zplev(ig)=kp(ig)*globzplevnew/p00  ! 2)
       ENDDO
       DO ig=1,klon     ! 3) unstable case condensation and sublimation: zplev=zplevhist
          IF (((pdpsrf(ig).lt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).le.ztsrf(ig))).or.  &
          ((pdpsrf(ig).gt.0.).and.(tcond_n2(zplev(ig),zqn2(ig,1)).ge.ztsrf(ig)))) then
             zplev(ig)=zplevhist(ig)
          ENDIF
          zplevhist(ig)=zplev(ig)
       ENDDO
       zplev=zplev*globzplevnew/glob_average2d(zplevhist)   ! 4)
       DO ig=1,klon    ! 5)
         zdtsrf(ig)=pdtsrf(ig) + (stephan/pcapcal(ig))*(ptsrf(ig)**4-ztsrf(ig)**4)
         ztsrf(ig)=ztsrf(ig)+pdtsrfc(ig)*subtimestep+zdtsrf(ig)*subtimestep
         zpicen2(ig)=zpicen2(ig)+pdicen2(ig)*subtimestep
       ENDDO
     ENDIF  ! (itsub=1)

   DO ig=1,klon
     ! forecast of frost temperature ztcondsol
     ztcondsol(ig) = tcond_n2(zplev(ig),zqn2(ig,1))
     !ztcondsol(ig) = tcond_n2(zplev(ig),vmrn2(ig))

!     Loop over where we have condensation / sublimation
     IF ((ztsrf(ig) .LT. ztcondsol(ig)) .OR.           &    ! ground cond
                ((ztsrf(ig) .GT. ztcondsol(ig)) .AND.  &    ! ground sublim
               (zpicen2(ig) .GT. 0.))) THEN
        condsub(ig) = .true.    ! condensation or sublimation

!     Condensation or partial sublimation of N2 ice
        if (ztsrf(ig) .LT. ztcondsol(ig)) then  ! condensation
!             Include a correction to account for the cooling of air near
!             the surface before condensing:
         zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
         /((lw_n2+cpp*(zt(ig,1)-ztcondsol(ig)))*subtimestep)
        else                                    ! sublimation
          zcondices(ig)=pcapcal(ig)*(ztcondsol(ig)-ztsrf(ig)) &
          /(lw_n2*subtimestep)
        end if

        pdtsrfc(ig) = (ztcondsol(ig) - ztsrf(ig))/subtimestep

!     partial sublimation of N2 ice
!     If the entire N_2 ice layer sublimes
!     (including what has just condensed in the atmosphere)
        IF((zpicen2(ig)/subtimestep).LE. &
            -zcondices(ig))THEN
           zcondices(ig) = -zpicen2(ig)/subtimestep
           pdtsrfc(ig)=(lw_n2/pcapcal(ig))*       &
               (zcondices(ig))
        END IF

!     Changing N2 ice amount and pressure

        pdicen2(ig) = zcondices(ig)
        pdpsrf(ig)   = -pdicen2(ig)*g
    !    pdpsrf(ig)   = 0. ! OPTION to check impact N2 sub/cond
        IF (zplev(ig)+pdpsrf(ig)*subtimestep.le.0.0000001) then
            pdpsrf(ig)=(0.0000001*kp(ig)/p00-zplev(ig))/subtimestep
            pdicen2(ig)=-pdpsrf(ig)/g
        ENDIF

     ELSE  ! no condsub
        pdpsrf(ig)=0.
        pdicen2(ig)=0.
        pdtsrfc(ig)=0.
     ENDIF
   ENDDO ! ig
  enddo ! subtimestep

!     Updating pressure, temperature and ice reservoir
  DO ig=1,klon
     pdpsrf(ig)=(zplev(ig)+pdpsrf(ig)*subtimestep-pplev(ig,1))/ptimestep
     ! Two options here : 1 ok, 2 is wrong
     pdicen2(ig)=(zpicen2(ig)+pdicen2(ig)*subtimestep-picen2(ig))/ptimestep
     !pdicen2(ig)=-pdpsrf(ig)/g

     pdtsrfc(ig)=((ztsrf(ig)+pdtsrfc(ig)*subtimestep)-(ztsrfhist(ig)))/ptimestep

! security
     if (picen2(ig) + pdicen2(ig)*ptimestep.lt.0.) then
       write(*,*) 'WARNING in condense_n2:'
       write(*,*) picen2(ig),pdicen2(ig)*ptimestep
       pdicen2(ig)= -picen2(ig)/ptimestep
       pdpsrf(ig)=-pdicen2(ig)*g
     endif

     if(.not.picen2(ig).ge.0.) THEN
!           if(picen2(ig) + pdicen2(ig)*ptimestep.le.-1.e-8) then
          print*, 'WARNING NEG RESERVOIR in condense_n2: picen2(',ig,')=', picen2(ig) + pdicen2(ig)*ptimestep
!              pdicen2(ig)= -picen2(ig)/ptimestep
!           else
          picen2(ig)=0.0
!           endif
     endif
  ENDDO

! ***************************************************************
!  Correction to account for redistribution between sigma or hybrid
!  layers when changing surface pressure (and warming/cooling
!  of the n2 currently changing phase).
! *************************************************************
  if (.not.fast) then
   DO ig=1,klon
    if (condsub(ig)) then

!  Mass fluxes through the sigma levels (kg.m-2.s-1)  (>0 when up)
!  """"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
        zmflux(1) = -zcondices(ig)
        DO l=1,klev
         zmflux(l+1) = zmflux(l) -zcondicea(ig,l)  &
         + (bp(l)-bp(l+1))*(zfallice(ig,1)-zmflux(1))
! 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.) zmflux(l+1)=0.
        END DO

! Mass of each layer
! ------------------
       DO l=1,klev
         masse(l)=(pplev(ig,l) - pplev(ig,l+1))/g
       END DO


!  Corresponding fluxes for T,U,V,Q
!  """"""""""""""""""""""""""""""""
!           averaging operator for TRANSPORT
!           """"""""""""""""""""""""""""""""

!           Subtimestep loop to perform the redistribution gently and simultaneously with
!           the other tendencies
!           Estimation of subtimestep (using only  the first layer, the most critical)
        dtmax=ptimestep
        if (zmflux(1).gt.1.e-20) then
            dtmax=min(dtmax,masse(1)*zqn2(ig,1)/abs(zmflux(1)))
        endif
        nsubtimestep= max(nint(ptimestep/dtmax),nint(2.))
        subtimestep=ptimestep/float(nsubtimestep)

!          New flux for each subtimestep
       do l=1,klev+1
            w(l)=-zmflux(l)*subtimestep
       enddo
!          initializing variables that will vary during subtimestep:
       do l=1,klev
           ztc(l)  =pt(ig,l)
           zu(l)   =pu(ig,l)
           zv(l)   =pv(ig,l)
           do iq=1,nq
              zq(l,iq) = pq(ig,l,iq)
           enddo
       end do

!          loop over nsubtimestep
!          """"""""""""""""""""""
       do itsub=1,nsubtimestep
!             Progressively adding tendancies from other processes.
          do l=1,klev
             ztc(l)  =ztc(l)   +(pdt(ig,l) + zdtlatent(ig,l))*subtimestep
             zu(l)   =zu(l)   +pdu( ig,l) * subtimestep
             zv(l)   =zv(l)   +pdv( ig,l) * subtimestep
             do iq=1,nq
                zq(l,iq) = zq(l,iq) + pdq(ig,l,iq)* subtimestep
             enddo
           end do

!             Change of mass in each layer
          do l=1,klev
             masse(l)=masse(l)+pdpsrf(ig)*subtimestep*(pplev(ig,l) - pplev(ig,l+1))&
                                             /(g*pplev(ig,1))
          end do

!             Value transfert at the surface interface when condensation sublimation:

          if (zmflux(1).lt.0) then
!               Surface condensation
            zum(1)= zu(1)
            zvm(1)= zv(1)
            ztm(1) = ztc(1)
          else
!               Surface sublimation:
            ztm(1) = ztsrf(ig) + pdtsrfc(ig)*ptimestep
            zum(1) = 0
            zvm(1) = 0
          end if
          do iq=1,nq
              zqm(1,iq)=0. ! most tracer do not condense !
          enddo
!             Special case if the tracer is n2 gas
          if (igcm_n2.ne.0) zqm(1,igcm_n2)=1.

           ztm(2:klev+1)=0.
           zum(2:klev+1)=0.
           zvm(2:klev+1)=0.
           zqm1(1:klev+1)=0.

!             Van Leer scheme:
          call vl1d(klev,ztc,2.,masse,w,ztm)
          call vl1d(klev,zu ,2.,masse,w,zum)
          call vl1d(klev,zv ,2.,masse,w,zvm)
         do iq=1,nq
           do l=1,klev
            zq1(l)=zq(l,iq)
           enddo
           zqm1(1)=zqm(1,iq)
           call vl1d(klev,zq1,2.,masse,w,zqm1)
           do l=2,klev
            zqm(l,iq)=zqm1(l)
           enddo
          enddo

!             Correction to prevent negative value for qn2
          if (igcm_n2.ne.0) then
            zqm(1,igcm_n2)=1.
            do l=1,klev-1
              if (w(l)*zqm(l,igcm_n2).gt.zq(l,igcm_n2)*masse(l)) then
                 zqm(l+1,igcm_n2)=max(zqm(l+1,igcm_n2),    &
                 (zqm(l,igcm_n2)*w(l) -zq(l,igcm_n2)*masse(l))/w(l+1) )
              else
                exit
              endif
            end do
           end if


          !!! Source haze: 0.02 pourcent when n2 sublimes
          IF (source_haze) THEN
           IF (pdicen2(ig).lt.0) THEN
            DO iq=1,nq
             tname=noms(iq)
             if (tname(1:4).eq."haze") then
                   !zqm(1,iq)=0.02
                   !zqm(1,iq)=-pdicen2(ig)*0.02
                   zqm(1,iq)=-pdicen2(ig)*ptimestep*0.02
                   !zqm(10,iq)=-pdicen2(ig)*ptimestep*100.
                   !zqm(1,iq)=-pdicen2(ig)*1000000.

             endif
            ENDDO
           ENDIF
          ENDIF
          ztm(klev+1)= ztc(klev) ! should not be used, but...
          zum(klev+1)= zu(klev)  ! should not be used, but...
          zvm(klev+1)= zv(klev)  ! should not be used, but...
          do iq=1,nq
           zqm(klev+1,iq)= zq(klev,iq)
          enddo

!             Tendencies on T, U, V, Q
!             """""""""""""""""""""""
          DO l=1,klev

!               Tendencies on T
            zdtsig(ig,l) = (1/masse(l)) *   &
           ( zmflux(l)*(ztm(l) - ztc(l))    &
           - zmflux(l+1)*(ztm(l+1) - ztc(l)) &
           + zcondicea(ig,l)*(ztcond(ig,l)-ztc(l))  )

!               Tendencies on U
              pduc(ig,l)   = (1/masse(l)) *   &
           ( zmflux(l)*(zum(l) - zu(l))&
           - zmflux(l+1)*(zum(l+1) - zu(l)) )

!               Tendencies on V
              pdvc(ig,l)   = (1/masse(l)) *   &
           ( zmflux(l)*(zvm(l) - zv(l))   &
           - zmflux(l+1)*(zvm(l+1) - zv(l)) )

          END DO

!             Tendencies on Q
          do iq=1,nq
            if (iq.eq.igcm_n2) then
!                 SPECIAL Case when the tracer IS N2 :
              DO l=1,klev
              pdqc(ig,l,iq)= (1/masse(l)) *        &
               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))    &
               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq))&
               + zcondicea(ig,l)*(zq(l,iq)-1.) )
              END DO
            else
              DO l=1,klev
                pdqc(ig,l,iq)= (1/masse(l)) *         &
               ( zmflux(l)*(zqm(l,iq) - zq(l,iq))     &
               - zmflux(l+1)*(zqm(l+1,iq) - zq(l,iq)) &
               + zcondicea(ig,l)*zq(l,iq) )
              END DO
            end if
          enddo
!             Update variables at the end of each subtimestep.
          do l=1,klev
            ztc(l)  =ztc(l)  + zdtsig(ig,l)  *subtimestep
            zu(l)   =zu(l)   + pduc(ig,l)  *subtimestep
            zv(l)   =zv(l)   + pdvc(ig,l)  *subtimestep
            do iq=1,nq
              zq(l,iq) = zq(l,iq) + pdqc(ig,l,iq)*subtimestep
            enddo
          end do
       enddo   ! loop on nsubtimestep
!          Recomputing Total tendencies
       do l=1,klev
         pdtc(ig,l)   = (ztc(l) - zt(ig,l) )/ptimestep
         pduc(ig,l)   = (zu(l) - (pu(ig,l) + pdu(ig,l)*ptimestep))/ptimestep
         pdvc(ig,l)   = (zv(l) - (pv(ig,l) + pdv(ig,l)*ptimestep))/ptimestep
         do iq=1,nq
           pdqc(ig,l,iq) = (zq(l,iq) - (pq(ig,l,iq) + pdq(ig,l,iq)*ptimestep))/ptimestep


!           Correction temporaire
            if (iq.eq.igcm_n2) then
             if((pq(ig,l,iq) +(pdqc(ig,l,iq)+ pdq(ig,l,iq))*ptimestep) &
                     .lt.0.01) then ! if n2 < 1 %  !
                     write(*,*) 'Warning: n2 < 1%'
                     pdqc(ig,l,iq)=(0.01-pq(ig,l,iq))/ptimestep-pdq(ig,l,iq)
             end if
            end if

         enddo
       end do
! *******************************TEMPORAIRE ******************
       if (klon.eq.1) then
              write(*,*) 'nsubtimestep=' ,nsubtimestep
           write(*,*) 'masse avant' , (pplev(ig,1) - pplev(ig,2))/g
              write(*,*) 'masse apres' , masse(1)
              write(*,*) 'zmflux*DT, l=1' ,  zmflux(1)*ptimestep
              write(*,*) 'zmflux*DT, l=2' ,  zmflux(2)*ptimestep
          write(*,*) 'pq, l=1,2,3' ,  pq(1,1,1), pq(1,2,1),pq(1,3,1)
          write(*,*) 'zq, l=1,2,3' ,  zq(1,1), zq(2,1),zq(3,1)
              write(*,*) 'dq*Dt l=1' ,  pdq(1,1,1)*ptimestep
              write(*,*) 'dqcond*Dt l=1' ,  pdqc(1,1,1)*ptimestep
       end if

! ***********************************************************
    end if ! if (condsub)
   END DO  ! loop on ig
  endif ! not fast

! ************ Option Olkin to prevent N2 effect in the south********
112   continue
  if (olkin) then
  DO ig=1,klon
     if (latitude(ig).lt.0) then
       pdtsrfc(ig) = max(0.,pdtsrfc(ig))
       pdpsrf(ig) = 0.
       pdicen2(ig)  = 0.
       do l=1,klev
         pdtc(ig,l)  =  max(0.,zdtlatent(ig,l))
         pduc(ig,l)  = 0.
         pdvc(ig,l)  = 0.
         do iq=1,nq
           pdqc(ig,l,iq) = 0
         enddo
       end do
     end if
  END DO
  end if
! *******************************************************************

! ***************************************************************
! Ecriture des diagnostiques
! ***************************************************************

!     DO l=1,klev
!        DO ig=1,klon
!         Taux de cond en kg.m-2.pa-1.s-1
!          tconda1(ig,l)=zcondicea(ig,l)/(pplev(ig,l)-pplev(ig,l+1))
!          Taux de cond en kg.m-3.s-1
!          tconda2(ig,l)=tconda1(ig,l)*pplay(ig,l)*g/(r*pt(ig,l))
!        END DO
!     END DO
!     call WRITEDIAGFI(klon,'tconda1',              &
!     'Taux de condensation N2 atmospherique /Pa',    &
!      'kg.m-2.Pa-1.s-1',3,tconda1)
!     call WRITEDIAGFI(klon,'tconda2',              &
!     'Taux de condensation N2 atmospherique /m',     &
!      'kg.m-3.s-1',3,tconda2)


  return
  end subroutine condense_n2

!-------------------------------------------------------------------------

  real function tcond_n2(p,vmr)
!   Calculates the condensation temperature for N2 at pressure P and vmr
  implicit none
  real, intent(in):: p,vmr

!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
  IF (p.ge.0.529995) then
!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
  ! tcond_n2 = (1.)/(1./63.1470-296.925/(2.5e5*0.98)*log(1./(0.125570*1.e5)*p*vmr))
    tcond_n2 = (1.)/(0.01583606505-1.211938776e-3*log(7.963685594e-5*p*vmr))
  ELSE
!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT BT
  ! tcond_n2 = (1.)/(1./35.6-296.925/(2.5e5*1.09)*log(1./(0.508059)*p*vmr))
    tcond_n2 = (1.)/(1./35.6-1.089633028e-3*log(1.968275338*p*vmr))
  ENDIF
  return
  end function tcond_n2

!-------------------------------------------------------------------------

  real function pcond_n2(t,vmr)
!   Calculates the condensation pressure for N2 at temperature T and vmr
  implicit none
  real, intent(in):: t,vmr

!       tcond_n2 = (1.)/(0.026315-0.0011877*log(.7143*p*vmr))
  IF (t.ge.35.6) then
!     tcond Fray and Schmitt for N2 phase beta (T>35.6 K) FIT TB
  ! pcond_n2 = 0.125570*1.e5/vmr*exp((2.5e5*0.98)/296.925*(1./63.1470-1./t))
    pcond_n2 = 0.125570e5/vmr*exp(825.1241896*(1./63.147-1./t))
  ELSE
!     tcond Fray and Schmitt for N2 phase alpha(T<35.6 K) FIT TB
  ! pcond_n2 = 0.508059/vmr*exp((2.5e5*1.09)/296.925*(1./35.6-1./t))
    pcond_n2 = 0.508059/vmr*exp(917.7401701*(1./35.6-1./t))
  ENDIF
  return
  end function pcond_n2

!-------------------------------------------------------------------------

  real function glob_average2d(var)
!   Calculates the global average of variable var
  use comgeomfi_h
  use dimphy, only: klon
  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
  use geometry_mod, only: cell_area, latitude

  implicit none

!      INTEGER klon
  REAL var(klon)
  INTEGER ig

  glob_average2d = 0.
  DO ig=2,klon-1
       glob_average2d = glob_average2d + var(ig)*cell_area(ig)
  END DO
  glob_average2d = glob_average2d + var(1)*cell_area(1)*nbp_lon
  glob_average2d = glob_average2d + var(klon)*cell_area(klon)*nbp_lon
  glob_average2d = glob_average2d/(totarea+(cell_area(1)+cell_area(klon))*(nbp_lon-1))

  end function glob_average2d

! *****************************************************************

  subroutine vl1d(klev,q,pente_max,masse,w,qm)
!
!     Operateur de moyenne inter-couche pour calcul de transport type
!     Van-Leer " pseudo amont " dans la verticale
!    q,w sont des arguments d'entree  pour le s-pg ....
!    masse : masse de la couche Dp/g
!    w : masse d'atm ``transferee'' a chaque pas de temps (kg.m-2)
!    pente_max = 2 conseillee
!   --------------------------------------------------------------------
  IMPLICIT NONE

!   Arguments:
!   ----------
  integer klev
  real masse(klev),pente_max
  REAL q(klev),qm(klev+1)
  REAL w(klev+1)
!
!      Local
!   ---------
!
  INTEGER l
!
  real dzq(klev),dzqw(klev),adzqw(klev),dzqmax
  real sigw, Mtot, MQtot
  integer m


!    On oriente tout dans le sens de la pression
!     W > 0 WHEN DOWN !!!!!!!!!!!!!

  do l=2,klev
        dzqw(l)=q(l-1)-q(l)
        adzqw(l)=abs(dzqw(l))
  enddo

  do l=2,klev-1
        if(dzqw(l)*dzqw(l+1).gt.0.) then
            dzq(l)=0.5*(dzqw(l)+dzqw(l+1))
        else
            dzq(l)=0.
        endif
        dzqmax=pente_max*min(adzqw(l),adzqw(l+1))
        dzq(l)=sign(min(abs(dzq(l)),dzqmax),dzq(l))
  enddo

     dzq(1)=0.
     dzq(klev)=0.

   do l = 1,klev-1

!         Regular scheme (transfered mass < layer mass)
!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      if(w(l+1).gt.0. .and. w(l+1).le.masse(l+1)) then
         sigw=w(l+1)/masse(l+1)
         qm(l+1)=(q(l+1)+0.5*(1.-sigw)*dzq(l+1))
      else if(w(l+1).le.0. .and. -w(l+1).le.masse(l)) then
         sigw=w(l+1)/masse(l)
         qm(l+1)=(q(l)-0.5*(1.+sigw)*dzq(l))

!         Extended scheme (transfered mass > layer mass)
!         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      else if(w(l+1).gt.0.) then
         m=l+1
         Mtot = masse(m)
         MQtot = masse(m)*q(m)
         !if (m.lt.klev) then ! because some compilers will have problems
         !                      ! evaluating masse(klev+1)
         do while ((m.lt.klev).and.(w(l+1).gt.(Mtot+masse(m+1))))
            m=m+1
            Mtot = Mtot + masse(m)
            MQtot = MQtot + masse(m)*q(m)
         !   if (m.eq.klev) exit
         end do
         !endif
         if (m.lt.klev) then
            sigw=(w(l+1)-Mtot)/masse(m+1)
            qm(l+1)= (1/w(l+1))*(MQtot + (w(l+1)-Mtot)* &
           (q(m+1)+0.5*(1.-sigw)*dzq(m+1)) )
         else
            w(l+1) = Mtot
            qm(l+1) = Mqtot / Mtot
            write(*,*) 'top layer is disapearing !'
            stop
         end if
      else      ! if(w(l+1).lt.0)
         m = l-1
         Mtot = masse(m+1)
         MQtot = masse(m+1)*q(m+1)
         if (m.gt.0) then ! because some compilers will have problems
                          ! evaluating masse(0)
          do while ((m.gt.0).and.(-w(l+1).gt.(Mtot+masse(m))))
            m=m-1
            Mtot = Mtot + masse(m+1)
            MQtot = MQtot + masse(m+1)*q(m+1)
            if (m.eq.0) exit
          end do
         endif
         if (m.gt.0) then
            sigw=(w(l+1)+Mtot)/masse(m)
            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)* &
           (q(m)-0.5*(1.+sigw)*dzq(m))  )
         else
            qm(l+1)= (-1/w(l+1))*(MQtot + (-w(l+1)-Mtot)*qm(1))
         end if
      end if
   enddo

   return
   end  subroutine vl1d