source: trunk/LMDZ.MARS/libf/phymars/lwi.F @ 1112

      subroutine lwi (ig0,kdlon,kflev
     .                ,psi,zdblay,pdp
     .                ,newpcolc )


      use dimradmars_mod, only: ndlo2, ndlon, nflev, nir
      use yomlw_h, only: gcp, nlaylte, xi
      implicit none


!#include "dimensions.h"
!#include "dimphys.h"
!#include "dimradmars.h"
#include "comg1d.h"
#include "callkeys.h"
#include "comcstfi.h"
!#include "yomlw.h"
 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C                             -   lwi    -    
C
C     PURPOSE:       Shema semi - implicite 
C
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


c************************************************************************
c
c        0.    Declarations
c              ------------
c
c-------------------------------------------------------------------------
c        0.1   Arguments
c              ---------
c
 
      integer ig0,kdlon,kflev

      real    psi(ndlo2,kflev)
     .     ,  zdblay(ndlo2,nir,kflev)
     .     ,  pdp(ndlo2,kflev)


      real    newpcolc(ndlo2,kflev)

c-------------------------------------------------------------------------
c        0.2   local arrays
c              ------------
c
      real    di(ndlon,nflev) 
     .      , hi(ndlon,nflev) 
     .      , bi(ndlon,nflev) 

      real    ci(ndlon,nflev) 
     .      , ai(ndlon,nflev) 
      real    deltat

      real   semit, denom

      integer i, jl

c************************************************************************
c
c        1.    Initialisations
c              ---------------
c
c-----------------------------------------------------------------------
 
        deltat = dtphys * iradia
c       print*,'SEMI = ',semi, '(expl:0  semi-implicite:0.5  impl:1)'
        semit = semi * deltat
c       semi = 0.

c       print*,'dtphys,iradia,deltat,semit:',dtphys,iradia,deltat,semit
c       print*,'g,cpp',g,cpp


c************************************************************************
c
c        2.    
c              ---------------
c
c-------------------------------------------------------------------------
c        2.1   Calcul des di
c              -------------
c


      do i = 1 , nlaylte-1
        do jl = 1 , kdlon
c     -------------------
      di(jl,i) =  1 + semit * (g / pdp(jl,i) / cpp) * (
     .    ( xi(ig0+jl,1,i,nlaylte+1)
     .    + xi(ig0+jl,1,i,i+1)
     .    + xi(ig0+jl,1,i,i-1) )
     .    *    zdblay(jl,1,i)
     .  + ( xi(ig0+jl,2,i,nlaylte+1)
     .    + xi(ig0+jl,2,i,i+1)
     .    + xi(ig0+jl,2,i,i-1) )
     .    *    zdblay(jl,2,i)
     .     )
c     -------------------
        enddo
      enddo

c couche nlaylte
c ------------
c      , on enleve i,i+1 sinon on a 2 fois le cooling2space

      do jl = 1 , kdlon
c     -------------------
      di(jl,nlaylte) =  1 + semit * (g / pdp(jl,nlaylte) / cpp) * (
     .    ( xi(ig0+jl,1,nlaylte,nlaylte+1)
     .    + xi(ig0+jl,1,nlaylte,nlaylte-1) )
     .    *    zdblay(jl,1,nlaylte)
     .  + ( xi(ig0+jl,2,nlaylte,nlaylte+1)
     .    + xi(ig0+jl,2,nlaylte,nlaylte-1) )
     .    *    zdblay(jl,2,nlaylte)
     .     )
c     -------------------
      enddo

c-------------------------------------------------------------------------
c        2.2   Calcul des hi
c              -------------
c

      do i = 1 , nlaylte-1
        do jl = 1 , kdlon
c     -------------------
      hi(jl,i) =    - semit * (g / pdp(jl,i) / cpp) *
     .            (    xi(ig0+jl,1,i,i+1) * zdblay(jl,1,i+1)   
     .               + xi(ig0+jl,2,i,i+1) * zdblay(jl,2,i+1)   )
c     -------------------
        enddo
      enddo

c-------------------------------------------------------------------------
c        2.3   Calcul des bi
c              -------------
c


      do i = 2 , nlaylte
        do jl = 1 , kdlon
c     -------------------
      bi(jl,i) =   - semit * (g / pdp(jl,i) / cpp) * 
     .           (     xi(ig0+jl,1,i,i-1) * zdblay(jl,1,i-1)   
     .               + xi(ig0+jl,2,i,i-1) * zdblay(jl,2,i-1)   )
c     -------------------
        enddo
      enddo


c couche 1
c --------
c  tant qu'on a pas un calcul propre de zdblay(0) qui tienne compte de 
c    la discontinuite de temperature au sol , on met  b1 = 0


      do jl = 1 , kdlon
        bi(jl,1) = 0 
      enddo

c-------------------------------------------------------------------------
c        2.4   
c              -------------
c

c couche nlaylte
c ------------

      do jl = 1 , kdlon
c     -------------------
      ci(jl,nlaylte) = (gcp * psi(jl,nlaylte) / pdp(jl,nlaylte))
     .                   / di(jl,nlaylte)

      ai(jl,nlaylte) = - bi(jl,nlaylte) / di(jl,nlaylte)
c     -------------------
      enddo



      do i = nlaylte-1 , 1 , -1
        do jl = 1 , kdlon
c     -------------------
      denom = di(jl,i) + hi(jl,i) * ai(jl,i+1)

      ci(jl,i) = (  gcp * psi(jl,i) / pdp(jl,i)
     .             - hi(jl,i) * ci(jl,i+1)  )  / denom
 
      ai(jl,i) = -bi(jl,i) / denom
c     -------------------
        enddo
      enddo


c-------------------------------------------------------------------------
c        2.5   
c              -------------
c

c couche 1
c -------
      do jl = 1 , kdlon
        newpcolc(jl,1) = ci(jl,1)
      enddo


      do i = 2 , nlaylte
        do jl = 1 , kdlon
           newpcolc(jl,i) = ci(jl,i) + ai(jl,i) * newpcolc(jl,i-1)
        enddo
      enddo



c-------------------------------------------------------------------------
      RETURN
      END