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

!
! $Id: iniacademic.F 1403 2010-07-01 09:02:53Z nfevrier $
!
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"
#include "logic.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.

        if (llm.eq.1) then
          ! specific initializations for the shallow water case
          kappa=1
        endif
        
        CALL iniconst
        CALL inigeom
        CALL inifilr

        if (llm.eq.1) then
          ! initialize fields for the shallow water case, if required
          if (.not.read_start) then
            phis(:)=0.
            q(:,:,1)=1.e-10
            q(:,:,2)=1.e-15
            q(:,:,3:nqtot)=0.
            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
          endif
        endif

        if (iflag_phys.eq.2) then
          ! initializations for the academic case
          ps(:)=1.e5
          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   ')

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