source: LMDZ4/trunk/libf/dyn3d/iniacademic.F @ 711

!
! $Header$
!
c
c
      SUBROUTINE iniacademic(nq,vcov,ucov,teta,q,masse,ps,phis,time_0)

c%W%    %G%
c=======================================================================
c
c   Author:    Frederic Hourdin      original: 15/01/93
c   -------
c
c   Subject:
c   ------
c
c   Method:
c   --------
c
c   Interface:
c   ----------
c
c      Input:
c      ------
c
c      Output:
c      -------
c
c=======================================================================
      IMPLICIT NONE
c-----------------------------------------------------------------------
c   Declararations:
c   ---------------

#include "dimensions.h"
#include "paramet.h"
#include "comvert.h"
#include "comconst.h"
#include "comgeom.h"
#include "academic.h"
#include "ener.h"
#include "temps.h"
#include "control.h"

c   Arguments:
c   ----------

      integer nq
      real time_0

c   variables dynamiques
      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
      REAL q(ip1jmp1,llm,nq)               ! champs advectes
      REAL ps(ip1jmp1)                       ! pression  au sol
      REAL masse(ip1jmp1,llm)                ! masse d'air
      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
      REAL pks(ip1jmp1)                      ! exner au  sol
      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
      REAL phis(ip1jmp1)                     ! geopotentiel au sol
      REAL phi(ip1jmp1,llm)                  ! geopotentiel





c   Local:
c   ------

      REAL ddsin,tetarappelj,tetarappell,zsig
      real tetajl(jjp1,llm)
      INTEGER i,j,l,lsup,ij

      real zz,ran1
      integer idum

      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm),zdtvr

c-----------------------------------------------------------------------

c
      time_0=0.

      im         = iim
      jm         = jjm
      day_ini    = 0
      omeg       = 4.*asin(1.)/86400.
      rad    = 6371229.
      g      = 9.8
      daysec = 86400.
      dtvr    = daysec/FLOAT(day_step)
      zdtvr=dtvr
      kappa  = 0.2857143
      cpp    = 1004.70885
      preff     = 101325.
      pa        =  50 000.
      etot0      = 0.
      ptot0      = 0.
      ztot0      = 0.
      stot0      = 0.
      ang0       = 0.
      pa         = 0.

      CALL inicons0
      CALL inigeom
      CALL inifilr

      ps=0.
      phis=0.
c---------------------------------------------------------------------

      taurappel=10.*daysec

c---------------------------------------------------------------------
c   Calcul de la temperature potentielle :
c   --------------------------------------

      DO l=1,llm
       zsig=ap(l)/preff+bp(l)
       if (zsig.gt.0.3) then
         lsup=l
         tetarappell=1./8.*(-log(zsig)-.5)
         DO j=1,jjp1
            ddsin=sin(rlatu(j))-sin(pi/20.)
            tetajl(j,l)=300.*(1+1./18.*(1.-3.*ddsin*ddsin)+tetarappell)
         ENDDO
        else
c   Choix isotherme au-dessus de 300 mbar
         do j=1,jjp1
            tetajl(j,l)=tetajl(j,lsup)*(0.3/zsig)**kappa
         enddo
        endif
      ENDDO

      do l=1,llm
         do j=1,jjp1
            do i=1,iip1
               ij=(j-1)*iip1+i
               tetarappel(ij,l)=tetajl(j,l)
            enddo
         enddo
      enddo

c     call dump2d(jjp1,llm,tetajl,'TEQ   ')

      ps=1.e5
      phis=0.
      CALL pression ( ip1jmp1, ap, bp, ps, p       )
      CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
      CALL massdair(p,masse)

c  intialisation du vent et de la temperature
      teta(:,:)=tetarappel(:,:)
      CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
      call ugeostr(phi,ucov)
      vcov=0.
      q(:,:,1   )=1.e-10
      q(:,:,2   )=1.e-15
      q(:,:,3:nq)=0.


c   perturbation al\351atoire sur la temp\351rature
      idum  = -1
      zz = ran1(idum)
      idum  = 0
      do l=1,llm
         do ij=iip2,ip1jm
            teta(ij,l)=teta(ij,l)*(1.+0.005*ran1(idum))
         enddo
      enddo

      do l=1,llm
         do ij=1,ip1jmp1,iip1
            teta(ij+iim,l)=teta(ij,l)
         enddo
      enddo



c     PRINT *,' Appel test_period avec tetarappel '
c     CALL  test_period ( ucov,vcov,tetarappel,q,p,phis )
c     PRINT *,' Appel test_period avec teta '
c     CALL  test_period ( ucov,vcov,teta,q,p,phis )

c   initialisation d'un traceur sur une colonne
      j=jjp1*3/4
      i=iip1/2
      ij=(j-1)*iip1+i
      q(ij,:,3)=1.

      return
      END
c-----------------------------------------------------------------------