source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/config/ppsrc/phys/forceWCfn.f @ 224

      subroutine forceWCfn(ngrid,nq,pplev,pt,dq,dqs)

      USE tracer_h

      implicit none

!==================================================================
!     
!     Purpose
!     -------
!     Force tracer conservation in a column for a given pair of 
!     delta q, delta q_s
!
!     Authors
!     ------- 
!     R. Wordsworth
!     
!==================================================================

!-----------------------------------------------------------------------
!   INCLUDE 'dimensions.h'
!
!   dimensions.h contient les dimensions du modele
!   ndm est tel que iim=2**ndm
!-----------------------------------------------------------------------

      INTEGER iim,jjm,llm,ndm

      PARAMETER (iim= 128,jjm=96,llm=64,ndm=1)

!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!   INCLUDE 'dimphys.h'

! ngridmx : number of horizontal grid points
! note: the -1/jjm term will be 0; unless jj=1
      integer, parameter :: ngridmx = (2+(jjm-1)*iim - 1/jjm)   
! nlayermx : number of atmospheric layers
      integer, parameter :: nlayermx = llm 
! nsoilmx : number of subterranean layers ! nsoilmx is now in comsoil_h
      !integer, parameter :: nsoilmx = 4 ! for a test 
      !integer, parameter :: nsoilmx = 18 ! for z1=0.0002 m, depth = 18 m => mars case 
      !integer, parameter :: nsoilmx = 13 ! for z1=0.03 m, depth = 104.8 m => earth case
!-----------------------------------------------------------------------

!-----------------------------------------------------------------------
! INCLUDE "comcstfi.h"

      common/comcstfi/pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
      common/comcstfi/avocado!,molrad,visc
      
      real pi,rad,g,r,cpp,rcp,dtphys,daysec,mugaz,omeg
      real avocado!,molrad,visc


      INTEGER ngrid,nq

      real masse, Wtot, Wdiff

      real pplev(ngrid,nlayermx+1)
      real pt(ngrid)

      real dqs(ngrid,nq) 
      real dq(ngrid,nlayermx,nq)

      integer iq, ig, ilay

      do iq=1,nq 
        do ig=1,ngrid
           Wtot = 0.0
           do ilay=1,nlayermx
              masse = (pplev(ig,ilay) - pplev(ig,ilay+1))/g
              Wtot  = Wtot + masse*dq(ig,ilay,iq)
           enddo
           Wdiff = Wtot + dqs(ig,iq)
         
           dqs(ig,iq) = dqs(ig,iq) - Wdiff
        enddo
      enddo

      end