source: LMDZ4/branches/LMDZ4V5.0-dev/libf/dyn3d/iniacademic.F @ 1299

!
! $Id: iniacademic.F 1299 2010-01-20 14:27:21Z fairhead $
!
c
c
      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)

      USE filtreg_mod
      USE infotrac, ONLY : nqtot
      USE control_mod
 

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 "iniprint.h"

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

      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,nqtot)               ! champs advectes
      REAL ps(ip1jmp1)                       ! pression  au sol
      REAL masse(ip1jmp1,llm)                ! masse d'air
      REAL phis(ip1jmp1)                     ! geopotentiel au sol

c   Local:
c   ------

      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 phi(ip1jmp1,llm)                  ! geopotentiel
      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-----------------------------------------------------------------------
! 1. Initializations for Earth-like case
! --------------------------------------
      if (planet_type=="earth") then
c
        time_0=0.
        day_ref=0
        annee_ref=0

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

        CALL iniconst
        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 ! of if (zsig.gt.0.3)
        ENDDO ! of DO l=1,llm

        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:nqtot)=0.


c   perturbation aleatoire sur la temperature
        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.
      
      else
        write(lunout,*)"iniacademic: planet types other than earth",
     &                 " not implemented (yet)."
        stop
      endif ! of if (planet_type=="earth")
      return
      END
c-----------------------------------------------------------------------