source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/gcm.F @ 1140

!
! $Header$
!
c
c
      PROGRAM gcm

#ifdef CPP_IOIPSL
      USE IOIPSL
#endif

      USE mod_const_mpi, ONLY: init_const_mpi
      USE parallel
      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
      USE infotrac
      USE mod_interface_dyn_phys
      USE mod_hallo
      USE Bands

      USE filtreg_mod

! Ehouarn: for now these only apply to Earth:
#ifdef CPP_EARTH
      USE mod_grid_phy_lmdz
      USE mod_phys_lmdz_omp_data, ONLY: klon_omp 
      USE dimphy
      USE comgeomphy
#endif
      IMPLICIT NONE

c      ......   Version  du 10/01/98    ..........

c             avec  coordonnees  verticales hybrides 
c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )

c=======================================================================
c
c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
c   -------
c
c   Objet:
c   ------
c
c   GCM LMD nouvelle grille
c
c=======================================================================
c
c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
c      et possibilite d'appeler une fonction f(y)  a derivee tangente
c      hyperbolique a la  place de la fonction a derivee sinusoidale.
c  ... Possibilite de choisir le schema pour l'advection de
c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
c
c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
c      Pour Van-Leer iadv=10
c
c-----------------------------------------------------------------------
c   Declarations:
c   -------------
#include "dimensions.h"
#include "paramet.h"
#include "comconst.h"
#include "comdissnew.h"
#include "comvert.h"
#include "comgeom.h"
#include "logic.h"
#include "temps.h"
#include "control.h"
#include "ener.h"
#include "description.h"
#include "serre.h"
#include "com_io_dyn.h"
#include "iniprint.h"
#include "tracstoke.h"

      INTEGER         longcles
      PARAMETER     ( longcles = 20 )
      REAL  clesphy0( longcles )
      SAVE  clesphy0



      REAL zdtvr
c      INTEGER nbetatmoy, nbetatdem,nbetat
      INTEGER nbetatmoy, nbetatdem

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

c variables dynamiques intermediaire pour le transport

c   variables pour le fichier histoire
      REAL dtav      ! intervalle de temps elementaire

      REAL time_0

      LOGICAL lafin
c      INTEGER ij,iq,l,i,j
      INTEGER i,j


      real time_step, t_wrt, t_ops

c      REAL rdayvrai,rdaym_ini,rday_ecri
c      LOGICAL first

      LOGICAL call_iniphys
      data call_iniphys/.true./

c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
c+jld variables test conservation energie
c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
C     Tendance de la temp. potentiel d (theta)/ d t due a la 
C     tansformation d'energie cinetique en energie thermique
C     cree par la dissipation
c      REAL dhecdt(ip1jmp1,llm)
c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
c      CHARACTER (len=15) :: ztit
c-jld 


      character (len=80) :: dynhist_file, dynhistave_file
      character (len=20) ::modname
      character (len=80) ::abort_message

C Calendrier
      LOGICAL true_calendar
      PARAMETER (true_calendar = .false.)

c-----------------------------------------------------------------------
c    variables pour l'initialisation de la physique :
c    ------------------------------------------------
      INTEGER ngridmx
      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
      REAL zcufi(ngridmx),zcvfi(ngridmx)
      REAL latfi(ngridmx),lonfi(ngridmx)
      REAL airefi(ngridmx)
      SAVE latfi, lonfi, airefi
      
      INTEGER :: ierr


c-----------------------------------------------------------------------
c   Initialisations:
c   ----------------

      abort_message = 'last timestep reached'
      modname = 'gcm'
      descript = 'Run GCM LMDZ'
      lafin    = .FALSE.
      dynhist_file = 'dyn_hist'
      dynhistave_file = 'dyn_hist_ave'


c-----------------------------------------------------------------------
c   Choix du calendrier
c   -------------------

#ifdef CPP_IOIPSL
      if (true_calendar) then
        call ioconf_calendar('gregorian')
      else
        call ioconf_calendar('360d')
      endif
#endif
c----------------------------------------------------------------------
c  lecture des fichiers gcm.def ou run.def
c  ---------------------------------------
c
! Ehouarn: dump possibility of using defrun
!#ifdef CPP_IOIPSL
      CALL conf_gcm( 99, .TRUE. , clesphy0 )
!#else
!      CALL defrun( 99, .TRUE. , clesphy0 )
!#endif
c
c
c------------------------------------
c   Initialisation partie parallele
c------------------------------------
      CALL init_const_mpi

      call init_parallel
      call Read_Distrib
! Ehouarn : temporarily (?) keep this only for Earth
      if (planet_type.eq."earth") then
#ifdef CPP_EARTH
        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
#endif
      endif ! of if (planet_type.eq."earth")
      CALL set_bands
      CALL Init_interface_dyn_phys
      CALL barrier

      if (mpi_rank==0) call WriteBands
      call SetDistrib(jj_Nb_Caldyn)

c$OMP PARALLEL
      call Init_Mod_hallo
c$OMP END PARALLEL

! Ehouarn : temporarily (?) keep this only for Earth
      if (planet_type.eq."earth") then
#ifdef CPP_EARTH
c$OMP PARALLEL
      call InitComgeomphy
c$OMP END PARALLEL 
#endif
      endif ! of if (planet_type.eq."earth")

      IF (config_inca /= 'none') THEN
#ifdef INCA
         call init_const_lmdz(
     $        nbtr,anneeref,dayref,
     $        iphysiq,day_step,nday)

         call init_inca_para(
     $        iim,jjm+1,llm,klon_glo,mpi_size,
     $        distrib_phys,COMM_LMDZ)
#endif
      END IF

c-----------------------------------------------------------------------
c   Initialisation des traceurs
c   ---------------------------
c  Choix du nombre de traceurs et du schema pour l'advection
c  dans fichier traceur.def, par default ou via INCA
      call infotrac_init

c Allocation de la tableau q : champs advectes   
      ALLOCATE(q(ip1jmp1,llm,nqtot))

c-----------------------------------------------------------------------
c   Lecture de l'etat initial :
c   ---------------------------

c  lecture du fichier start.nc
      if (read_start) then
      ! we still need to run iniacademic to initialize some
      ! constants & fields, if we run the 'newtonian' case:
        if (iflag_phys.eq.2) then
          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
        endif
!#ifdef CPP_IOIPSL
        if (planet_type.eq."earth") then
#ifdef CPP_EARTH
! Load an Earth-format start file
         CALL dynetat0("start.nc",vcov,ucov,
     .              teta,q,masse,ps,phis, time_0)
#endif
        endif ! of if (planet_type.eq."earth")
c       write(73,*) 'ucov',ucov
c       write(74,*) 'vcov',vcov
c       write(75,*) 'teta',teta
c       write(76,*) 'ps',ps
c       write(77,*) 'q',q

      endif ! of if (read_start)

c le cas echeant, creation d un etat initial
      IF (prt_level > 9) WRITE(lunout,*)
     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
      if (.not.read_start) then
         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
      endif

c-----------------------------------------------------------------------
c   Lecture des parametres de controle pour la simulation :
c   -------------------------------------------------------
c  on recalcule eventuellement le pas de temps

      IF(MOD(day_step,iperiod).NE.0) THEN
        abort_message = 
     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
        call abort_gcm(modname,abort_message,1)
      ENDIF

      IF(MOD(day_step,iphysiq).NE.0) THEN
        abort_message = 
     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
        call abort_gcm(modname,abort_message,1)
      ENDIF

      zdtvr    = daysec/FLOAT(day_step)
        IF(dtvr.NE.zdtvr) THEN
         WRITE(lunout,*)
     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
        ENDIF

C
C on remet le calendrier à zero si demande
c
      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
        write(lunout,*)
     .  ' Attention les dates initiales lues dans le fichier'
        write(lunout,*)
     .  ' restart ne correspondent pas a celles lues dans '
        write(lunout,*)' gcm.def'
        if (raz_date .ne. 1) then
          write(lunout,*)
     .    ' On garde les dates du fichier restart'
        else
          annee_ref = anneeref
          day_ref = dayref
          day_ini = dayref
          itau_dyn = 0
          itau_phy = 0
          time_0 = 0.
          write(lunout,*)
     .   ' On reinitialise a la date lue dans gcm.def'
        endif
      ELSE
        raz_date = 0
      endif


c  nombre d'etats dans les fichiers demarrage et histoire
      nbetatdem = nday / iecri
      nbetatmoy = nday / periodav + 1

c-----------------------------------------------------------------------
c   Initialisation des constantes dynamiques :
c   ------------------------------------------
      dtvr = zdtvr
      CALL iniconst

c-----------------------------------------------------------------------
c   Initialisation de la geometrie :
c   --------------------------------
      CALL inigeom

c-----------------------------------------------------------------------
c   Initialisation du filtre :
c   --------------------------
      CALL inifilr
c
c-----------------------------------------------------------------------
c   Initialisation de la dissipation :
c   ----------------------------------

      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
     *                tetagdiv, tetagrot , tetatemp              )

c-----------------------------------------------------------------------
c   Initialisation de la physique :
c   -------------------------------
      IF (call_iniphys.and.iflag_phys.eq.1) THEN
         latfi(1)=rlatu(1)
         lonfi(1)=0.
         zcufi(1) = cu(1)
         zcvfi(1) = cv(1)
         DO j=2,jjm
            DO i=1,iim
               latfi((j-2)*iim+1+i)= rlatu(j)
               lonfi((j-2)*iim+1+i)= rlonv(i)
               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
            ENDDO
         ENDDO
         latfi(ngridmx)= rlatu(jjp1)
         lonfi(ngridmx)= 0.
         zcufi(ngridmx) = cu(ip1jm+1)
         zcvfi(ngridmx) = cv(ip1jm-iim)
         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)

         WRITE(lunout,*)
     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
! Earth:
         if (planet_type.eq."earth") then
#ifdef CPP_EARTH
         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
#endif
         endif ! of if (planet_type.eq."earth")
         call_iniphys=.false.
      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))


c-----------------------------------------------------------------------
c   Initialisation des dimensions d'INCA :
c   --------------------------------------
      IF (config_inca /= 'none') THEN
!$OMP PARALLEL
#ifdef INCA
         CALL init_inca_dim(klon_omp,llm,iim,jjm,
     $        rlonu,rlatu,rlonv,rlatv)
#endif
!$OMP END PARALLEL
      END IF

c-----------------------------------------------------------------------
c   Initialisation des I/O :
c   ------------------------


      day_end = day_ini + nday
      WRITE(lunout,300)day_ini,day_end
 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)

!#ifdef CPP_IOIPSL
      if (planet_type.eq."earth") then
#ifdef CPP_EARTH
      CALL dynredem0_p("restart.nc", day_end, phis)
#endif
      endif

      ecripar = .TRUE.

#ifdef CPP_IOIPSL
      if ( 1.eq.1) then
      time_step = zdtvr
      t_ops = iecri * daysec
      t_wrt = iecri * daysec
      CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
     .              t_ops, t_wrt, histid, histvid)

      t_ops = iperiod * time_step
      t_wrt = periodav * daysec
      CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
     .              t_ops, t_wrt, histaveid)

      dtav = iperiod*dtvr/daysec
      endif


#endif
! #endif of #ifdef CPP_IOIPSL

c  Choix des frequences de stokage pour le offline
c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
      istphy=istdyn/iphysiq     


c
c-----------------------------------------------------------------------
c   Integration temporelle du modele :
c   ----------------------------------

c       write(78,*) 'ucov',ucov
c       write(78,*) 'vcov',vcov
c       write(78,*) 'teta',teta
c       write(78,*) 'ps',ps
c       write(78,*) 'q',q

c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
     .              time_0)
c$OMP END PARALLEL


      END