source: LMDZ5/trunk/libf/dyn3dpar/iniacademic.F @ 1454

!
! $Id: iniacademic.F 1454 2010-11-18 12:01:24Z fairhead $
!
c
c
      SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)

      USE filtreg_mod
      USE infotrac, ONLY : nqtot
      USE control_mod, ONLY: day_step,planet_type
#ifdef CPP_IOIPSL
      USE IOIPSL
#else
! if not using IOIPSL, we still need to use (a local version of) getin
      USE ioipsl_getincom
#endif
      USE Write_Field
 

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,tetastrat,zsig,tetapv,w_pv  ! variables auxiliaires
      real tetajl(jjp1,llm)
      INTEGER i,j,l,lsup,ij

      REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
      REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
      LOGICAL ok_geost             ! Initialisation vent geost. ou nul
      LOGICAL ok_pv                ! Polar Vortex
      REAL phi_pv,dphi_pv,gam_pv   ! Constantes pour polar vortex 
      
      real zz,ran1
      integer idum

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

c-----------------------------------------------------------------------
! 1. Initializations for Earth-like case
! --------------------------------------
c
        ! initialize planet radius, rotation rate,...
        call conf_planete

        time_0=0.
        day_ref=1
        annee_ref=0

        im         = iim
        jm         = jjm
        day_ini    = 1
        dtvr    = daysec/REAL(day_step)
        zdtvr=dtvr
        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(:,:,:)=0
            CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
          endif
        endif

        if (iflag_phys.eq.2) then
          ! initializations for the academic case
          
!         if (planet_type=="earth") then

          ! 1. local parameters
          ! by convention, winter is in the southern hemisphere
          ! Geostrophic wind or no wind?
          ok_geost=.TRUE.
          CALL getin('ok_geost',ok_geost)
          ! Constants for Newtonian relaxation and friction
          k_f=1.                !friction 
          CALL getin('k_j',k_f)
          k_f=1./(daysec*k_f)
          k_c_s=4.  !cooling surface
          CALL getin('k_c_s',k_c_s)
          k_c_s=1./(daysec*k_c_s)
          k_c_a=40. !cooling free atm
          CALL getin('k_c_a',k_c_a)
          k_c_a=1./(daysec*k_c_a)
          ! Constants for Teta equilibrium profile
          teta0=315.     ! mean Teta (S.H. 315K)
          CALL getin('teta0',teta0)
          ttp=200.       ! Tropopause temperature (S.H. 200K)
          CALL getin('ttp',ttp)
          eps=0.         ! Deviation to N-S symmetry(~0-20K)
          CALL getin('eps',eps)
          delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
          CALL getin('delt_y',delt_y)
          delt_z=10.     ! Vertical Gradient (S.H. 10K)
          CALL getin('delt_z',delt_z)
          ! Polar vortex
          ok_pv=.false.
          CALL getin('ok_pv',ok_pv)
          phi_pv=-50.            ! Latitude of edge of vortex
          CALL getin('phi_pv',phi_pv)
          phi_pv=phi_pv*pi/180.
          dphi_pv=5.             ! Width of the edge
          CALL getin('dphi_pv',dphi_pv)
          dphi_pv=dphi_pv*pi/180.
          gam_pv=4.              ! -dT/dz vortex (in K/km)
          CALL getin('gam_pv',gam_pv)
          
          ! 2. Initialize fields towards which to relax
          ! Friction
          knewt_g=k_c_a
          DO l=1,llm
           zsig=presnivs(l)/preff
           knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
           kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
          ENDDO
          DO j=1,jjp1
            clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
          ENDDO
          
          ! Potential temperature 
          DO l=1,llm
           zsig=presnivs(l)/preff
           tetastrat=ttp*zsig**(-kappa)
           tetapv=tetastrat
           IF ((ok_pv).AND.(zsig.LT.0.1)) THEN
               tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
           ENDIF
           DO j=1,jjp1
             ! Troposphere
             ddsin=sin(rlatu(j))
             tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin
     &           -delt_z*(1.-ddsin*ddsin)*log(zsig)
             ! Profil stratospherique isotherme (+vortex)
             w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
             tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
             tetajl(j,l)=MAX(tetajl(j,l),tetastrat)  
           ENDDO
          ENDDO ! of DO l=1,llm
 
!          CALL writefield('theta_eq',tetajl)

          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


!         else
!          write(lunout,*)"iniacademic: planet types other than earth",
!     &                   " not implemented (yet)."
!          stop
!         endif ! of if (planet_type=="earth")

          ! 3. Initialize fields (if necessary)
          IF (.NOT. read_start) THEN
            ! surface pressure
            ps(:)=preff
            ! ground geopotential
            phis(:)=0.
            
            CALL pression ( ip1jmp1, ap, bp, ps, p       )
            CALL exner_hyb( ip1jmp1, ps, p,alpha,beta, pks, pk, pkf )
            CALL massdair(p,masse)

            ! bulk initialization of temperature
            teta(:,:)=tetarappel(:,:)
            
            ! geopotential
            CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
            
            ! winds
            if (ok_geost) then
              call ugeostr(phi,ucov)
            else
              ucov(:,:)=0.
            endif
            vcov(:,:)=0.
            
            ! bulk initialization of tracers
            if (planet_type=="earth") then
              ! Earth: first two tracers will be water
              do i=1,nqtot
                if (i.eq.1) q(:,:,i)=1.e-10
                if (i.eq.2) q(:,:,i)=1.e-15
                if (i.gt.2) q(:,:,i)=0.
              enddo
            else
              q(:,:,:)=0
            endif ! of if (planet_type=="earth")

            ! add random perturbation to 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

            ! maintain periodicity in longitude
            do l=1,llm
              do ij=1,ip1jmp1,iip1
                teta(ij+iim,l)=teta(ij,l)
              enddo
            enddo

          ENDIF ! of IF (.NOT. read_start)
        endif ! of if (iflag_phys.eq.2)
        
      END
c-----------------------------------------------------------------------