!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
*** S/P LIN_PHYEXE1

      subroutine LIN_PHYEXE1 (e,   d,   f,   v, 1,11
     $                        esiz,dsiz,fsiz,vsiz,
     $                        dt,trnch,kount,task,ni,nk)
*
#include "impnone.cdk"
*
      integer esiz,dsiz,fsiz,vsiz,trnch,kount,task,ni,nk
      real e(esiz), d(dsiz), f(fsiz), v(vsiz)
      real dt
*
*Author
*     Stephane Laroche - Janvier 2001
*
*Revisions
* 001   S. Laroche   - phyexe1 for the simplified physics
*                      from version 3.6.8
* 002   A. Zadra (May 2002) Add subgrid orography component
* 003   J.-F. Mahfouf (September 2002) 
*       Driver for moist processes (lin_precip1) 
*       [deep convection + large scale condensation]              
*
*Object
*          this is the main interface subroutine for the
*          CMC/RPN unified physics
*
*Arguments
*
*          - Input -
* E        entry    input field
* D        dynamics input field
*
*          - Input/Output -
* F        historic variables for the physics
*
*          - Output -
* V        physics tendencies and other output fields from the physics
*
*          - Input -
* ESIZ     dimension of e
* DSIZ     dimension of d
* FSIZ     dimension of f
* VSIZ     dimension of v
* DT       timestep (sec.)
* TRNCH    slice number
* KOUNT    timestep number
* ICPU     cpu number executing slice "trnch"
* N        horizontal running length
* NK       vertical dimension
*
*Notes
*          LIN_PHYEXE1 is called by the GEM model for the simplified physics.
*          It returns tendencies to the dynamics.
*
*IMPLICITES
*
#include "phy_macros_f.h"
#include "phybus.cdk"
#include "workspc.cdk"
#include "options.cdk"
#include "consphy.cdk"
*
*MODULES
**
*
      integer i,j,k,icpu
      real cdt1, rcdt1
      logical tvecst
*
*-----------------------------------------------------------------------
*
************************************************************************
*     AUTOMATIC ARRAYS
************************************************************************
*
      AUTOMATIC ( HUPLUS0 , REAL , (NI,NK  ) )
      AUTOMATIC (  TPLUS0 , REAL , (NI,NK  ) )
      AUTOMATIC (  UPLUS0 , REAL , (NI,NK  ) )
      AUTOMATIC (  VPLUS0 , REAL , (NI,NK  ) )
*
      AUTOMATIC ( GZMOINS , REAL , (NI,NK  ) )
      AUTOMATIC ( QE      , REAL , (NI,NK  ) )
      AUTOMATIC ( SELOC   , REAL , (NI,NK  ) )
      AUTOMATIC ( TRAV1D  , REAL , (NI     ) )
      AUTOMATIC ( TRAV2D  , REAL , (NI,NK,4) )
      AUTOMATIC ( TVIRT   , REAL , (NI,NK  ) )
      AUTOMATIC ( WORK    , REAL , (ESPWORK) )
*
************************************************************************
*
*
      external difuvd9,inichamp1,mfotvt,integ2
      external lin_kdif_simp1,lin_difver1,lin_conds1
      external lin_sgo1
*
*
      integer ik
      ik(i,k) = (k-1)*ni + i - 1
*
************************************************************************
*     preliminary calculations                                         *
*     ------------------------                                         *
************************************************************************
*
*
      icpu = task
*
*
************************************************************************
*     save fieds at time t*                                            *
*     ---------------------                                            *
************************************************************************
*
      do k=1,nk
         do i = 1,ni
           huplus0(i,k) =  d(huplus+ik(i,k))
            tplus0(i,k) =  d( tplus+ik(i,k))
            uplus0(i,k) =  d( uplus+ik(i,k))
            vplus0(i,k) =  d( vplus+ik(i,k))
         enddo
      enddo
*
*
*
************************************************************************
*     constant related to time step                                    *
*     -----------------------------                                    *
************************************************************************
*
      cdt1  = factdt * dt
      rcdt1 = 1./cdt1
*
*
************************************************************************
*     initialisations                                                  *
*     ---------------                                                  *
************************************************************************
*
*     calcul des niveaux intermediaires
      call difuvd9(seloc,.true.,d(sigm),ni,nk,nk-1)
*
       call lin_inichamp1 (e, esiz, f, fsiz, 
     $                     v, vsiz, d, dsiz,
     $                     trav2d, seloc, kount, trnch,
     $                     dt, cdt1, ni, nk)
*
*
************************************************************************
*     extraction des tendances de la dynamique                         *
*     ----------------------------------------                         *
************************************************************************
*
*
CSTL  Missing section
*
*
************************************************************************
*     calculs radiatifs                                                *
*     -----------------                                                *
************************************************************************
*
*
CSTL  Missing section
*
*
************************************************************************
*     certains preparatifs                                             *
*     --------------------                                             *
************************************************************************
*     calcul de la temperature virtuelle (tve),
*     de l'humidite specifique (qe) et des hauteurs
*     geopotentielles (ze) aux niveaux decales
*     [plus facteur de coriolis et za (Stephane Laroche)]
*
      do k=1,nk-2
         do i=1,ni
           v(tve+ik(i,k)) = ( d(tmoins+ik(i,k))       +
     $                        d(tmoins+ik(i,k+1)))/2.
           qe(i,k       ) = ( d(humoins+ik(i,k  ))    +
     $                        d(humoins+ik(i,k+1)))/2.
         end do
      end do
*
      do i=1,ni
        v(tve+ik(i,nk-1)) = d( tmoins+ik(i,nk-1))
       qe(       i,nk-1)  = d(humoins+ik(i,nk-1))
        v(fcor+i-1)       = 1.45441e-4*sin(f(dlat+i-1))
        v(za+i-1)         = -rgasd/grav* v(tve+ik(i,nk-1)) *
     $                      alog(d(sigm+ik(i,nk-1)))
      enddo
*
      call mfotvt(v(tve),v(tve),qe,ni,nk-1,ni)
*
************************************************************************
*    champ de temparature invariante si tvecst.eq.true pour
*    definir les chamgements de variable de z a sigma
*
      tvecst=.false.
      if(tvecst) then

        do k=1,nk-1
           do i=1,ni
             v(tve+ik(i,k)) = 288.15*d(sigm+ik(i,nk-1))**0.19
           end do
        end do
        do i=1,ni
          v(za+i-1)         = -rgasd/grav* v(tve+ik(i,nk-1)) *
     $                        alog(d(sigm+ik(i,nk-1)))
        enddo

      endif
*
*
************************************************************************
*
      do i=1,ni
         v(ze+ik(i,nk-1)) = 0.
      end do
*
      call integ2 ( v(ze), v(tve), -rgasd/grav, seloc,
     $              trav2d(1,1,1),trav2d(1,1,2),trav2d(1,1,3),
     $              ni, ni, ni, nk-1, .true. )
*
*
*
************************************************************************
*     processus de surface                                             *
*     --------------------                                             *
************************************************************************
*
*
CSTL  Missing section
*
*
***********************************************************************
*     flux de surface et coefficients de diffusion                    *
*     simplifies et lineaires                                         *
*     --------------------------------------------                    *
***********************************************************************
*
      if(lin_pbl.eq.1.or.lin_pbl.eq.2) then

        call lin_kdif_simp1 ( d, dsiz, f, fsiz, v, vsiz, ni, nk-1 )

      endif
* 
*
*
*
************************************************************************
*     sub-grid scale orography                                         *
*     ------------------------                                         *
************************************************************************
*
      if (lin_sgo.ne.0) then
*
*     calcul de la temperature virtuelle au temps moins
        call mfotvt(tvirt,d(tplus),d(huplus),ni,nk,ni)        
*
        call lin_sgo1 ( f, fsiz, d(uplus), d(vplus), tvirt, 
     +                  d(pplus), d(sigm), v(ugwd),  v(vgwd), 
     +                  cdt1,     kount,   trnch,    ni, 
     +                  ni,       nk-1,    icpu,     lin_sgo )
*
        do i=1,ni
           v(ugwd+ik(i,nk))  = 0.
           v(vgwd+ik(i,nk))  = 0.
        end do
*
      endif
*
************************************************************************
*     diffusion verticale                                              *
*     -------------------                                              *
*                                                                      *
*     note: la diffusion n'est pas appliquee sur les champs            *
*           dynamiques comme cela est fait dans la physique            *
*           complete                                                   *
*                                                                      *
************************************************************************
*
      if(lin_pbl.ge.1) then
*
        call lin_difver1 (d, dsiz, f, fsiz, v, vsiz, 
     $                    work, espwork, seloc,
     $                    cdt1, kount, trnch, ni, nk-1, icpu )
*
*
*      calcul des tendances de la diffusion au niveau diagnostique
*
         do i=1,ni
              v( qdifv+ik(i,nk)) = (f(qdiag+i-1)-d(huplus+ik(i,nk)))*rcdt1
              v( tdifv+ik(i,nk)) = (f(tdiag+i-1)-d( tplus+ik(i,nk)))*rcdt1
              v( udifv+ik(i,nk)) = (f(udiag+i-1)-d( uplus+ik(i,nk)))*rcdt1
              v( vdifv+ik(i,nk)) = (f(vdiag+i-1)-d( vplus+ik(i,nk)))*rcdt1
         end do

       endif
*
*
************************************************************************
*     cconvection/condensation processes                               *
*     ----------------------------------                               *
************************************************************************
*
*
      if(lin_lsc.ge.1) then
*
*     calcul de la temperature virtuelle au temps moins
        call mfotvt(tvirt,d(tmoins),d(humoins),ni,nk,ni)
*
*     calcul du geopotentiel au temps moins
        do  i=1,ni
           gzmoins(i,nk) = 0.
        end do
*
        call integ2 (gzmoins, tvirt, -rgasd, d(sigm),
     $               trav2d(1,1,1), trav2d(1,1,2), trav2d(1,1,3),
     $               ni, ni, ni, nk, .true. )
*
        call lin_precip1(d, dsiz, f, fsiz, v, vsiz, 
     $                   work, espwork, gzmoins,
     $                   dt, ni, ni, nk-1, 
     $                   kount, trnch, icpu)
*
      endif
**
************************************************************************
*     Accumulation des flux                                            *
*     ---------------------                                            *
************************************************************************
*
      if (kount.gt.0) then
*
*VDIR NODEP
         do i=0,ni-1
*
*           taux des precipitations convectives
            v(rc+i) = (f(tsc+i) + f(tlc+i)) * 1.0E-3
*
*           taux des precipitations stratiformes
            v(rr+i) = (f(tss+i) + f(tls+i)) * 1.0E-3
*
*           taux total
            v(rt+i) = v(rc+i) + v(rr+i)
*
*           asc : accumulation des precipitations solides convectives
            f(asc+i) = f(asc+i) + f(tsc+i) * dt * 1.0E-3
*
*           ass : accumulation des precipitations solides stratiformes
            f(ass+i) = f(ass+i) + f(tss+i) * dt * 1.0E-3
*
*           alc : accumulation des precipitations liquides convectives
            f(alc+i) = f(alc+i) + f(tlc+i) * dt * 1.0E-3
*
*           als : accumulation des precipitations liquides stratiformes
            f(als+i) = f(als+i) + f(tls+i) * dt * 1.0E-3
*
*           pc : accumulation des precipitations convectives
            v(pc+i) = f(alc+i) + f(asc+i)
*
*           ae : accumulation des precipitations stratiformes
            v(ae+i) = f(als+i) + f(ass+i)
*
*           pr : accumulation des precipitations totales
            v(pr+i) = v(pc+i) + v(ae+i)
*
         end do
*
      endif
*      
*
************************************************************************
*     Tendances moyennees                                              *
*     -------------------                                              *
************************************************************************
*
*
CSTL  Missing section
*
*
************************************************************************
*     extraction des variables dynamiques du temps plus                *
*     -------------------------------------------------                *
************************************************************************
*
*
CSTL  Missing section
*
*
************************************************************************
*     Bring back the fields at time t*                                 *
*     --------------------------------                                 *
************************************************************************
*
      do k=1,nk
         do i = 1,ni
           d(huplus+ik(i,k)) = huplus0(i,k)
           d( tplus+ik(i,k)) =  tplus0(i,k)
           d( uplus+ik(i,k)) =  uplus0(i,k)
           d( vplus+ik(i,k)) =  vplus0(i,k)
         enddo
      enddo
*
      return
      end