MODULE startvar
    !
    !
    !      There are three ways to access data from the database of atmospheric data which 
    !       can be used to initialize the model. This depends on the type of field which needs 
    !       to be extracted. In any case the call should come after a restget and should be of the type :
    !                CALL startget(...)
    !
    !       We will details the possible arguments to startget here :
    !
    !        - A 2D variable on the dynamical grid :
    !           CALL startget(varname, iml, jml, lon_in, lat_in, champ, val_ex)
    !            
    !
    !        - A 1D variable on the physical grid :
    !            CALL startget(varname, iml, jml, lon_in, lat_in, nbindex, champ, val_exp)
    !
    !
    !         - A 3D variable on the dynamical grid :
    !            CALL startget(varname, iml, jml, lon_in, lat_in, lml, pls, workvar, champ, val_exp)
    !
    !
    !         There is special constraint on the atmospheric data base except that the 
    !         the data needs to be in netCDF and the variables should have the the following 
    !        names in the file :
    !
    !      'RELIEF'  : High resolution orography 
    !       'ST'            : Surface temperature
    !       'CDSW'     : Soil moisture
    !       'Z'               : Surface geopotential
    !       'SP'            : Surface pressure
    !        'U'              : East ward wind
    !        'V'              : Northward wind
    !        'TEMP'             : Temperature
    !        'R'             : Relative humidity
    !      
      USE ioipsl
    !
    !
      IMPLICIT NONE
    !
    !
      PRIVATE
      PUBLIC startget
    !
    !
      INTERFACE startget
        MODULE PROCEDURE startget_phys2d, startget_phys1d, startget_dyn
      END INTERFACE
    !
      INTEGER, SAVE :: fid_phys, fid_dyn
      INTEGER, SAVE  :: iml_phys, iml_rel, iml_dyn
      INTEGER, SAVE :: jml_phys,  jml_rel, jml_dyn
      INTEGER, SAVE ::  llm_dyn, ttm_dyn
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lon_phys, lon_rug,
     . lon_alb, lon_rel, lon_dyn
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: lat_phys, lat_rug,
     . lat_alb, lat_rel, lat_dyn
      REAL, ALLOCATABLE, SAVE, DIMENSION (:)  :: lev_dyn
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: relief, zstd, zsig,
     . zgam, zthe, zpic, zval
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: rugo, masque, phis
    !
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:)  :: tsol, qsol, psol_dyn
      REAL, ALLOCATABLE, SAVE, DIMENSION (:,:,:)  ::   var_ana3d
    !
      CONTAINS
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in,
     . champ, val_exp)
    !
    !    There is a big mess with the size in logitude, should it be iml or iml+1.
    !    I have chosen to use the iml+1 as an argument to this routine and we declare
    !   internaly smaler fields when needed. This needs to be cleared once and for all in LMDZ. 
    !  A convention is required.
    !
    !
      CHARACTER*(*), INTENT(in) :: varname
      INTEGER, INTENT(in) :: iml, jml
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
      REAL, INTENT(inout) :: champ(iml,jml)
      REAL, INTENT(in) :: val_exp
    !
    !   This routine only works if the variable does not exist or is constant
    !
      IF ( MINVAL(champ(:,:)).EQ.MAXVAL(champ(:,:)) .AND. 
     .MINVAL(champ(:,:)).EQ.val_exp ) THEN
          !
          SELECTCASE(varname)
              !
              CASE ('relief')
                  !
                  !  If we do not have the orography we need to get it
                  !
                  IF ( .NOT.ALLOCATED(relief)) THEN
                      !
                      CALL start_init_orog( iml, jml, lon_in, lat_in)
                      !
                  ENDIF
                  !
                  IF (SIZE(relief) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*) 'STARTVAR module has been',
     .' initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  champ(:,:) = relief(:,:)
                  !
              CASE ('rugosite')
                  !
                  !  If we do not have the orography we need to get it
                  !
                  IF ( .NOT.ALLOCATED(rugo)) THEN
                      !
                      CALL start_init_orog( iml, jml, lon_in, lat_in)
                      !
                  ENDIF
                  !
                  IF (SIZE(rugo) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*) 
     .  'STARTVAR module has been initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  champ(:,:) = rugo(:,:)
                  !
              CASE ('masque')
                  !
                  !  If we do not have the orography we need to get it
                  !
                  IF ( .NOT.ALLOCATED(masque)) THEN
                      !
                      CALL start_init_orog( iml, jml, lon_in, lat_in)
                      !
                  ENDIF
                  !
                  IF (SIZE(masque) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*) 
     .   'STARTVAR module has been initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  champ(:,:) = masque(:,:)
                  !
              CASE ('surfgeo')
                  !
                  !  If we do not have the orography we need to get it
                  !
                  IF ( .NOT.ALLOCATED(phis)) THEN
                      !
                      CALL start_init_orog( iml, jml, lon_in, lat_in)
                      !
                  ENDIF
                  !
                  IF (SIZE(phis) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  champ(:,:) = phis(:,:)
                  !
              CASE ('psol')
                  !
                  !  If we do not have the orography we need to get it
                  !
                  IF ( .NOT.ALLOCATED(psol_dyn)) THEN
                      !
                      CALL start_init_dyn( iml, jml, lon_in, lat_in)
                      !
                  ENDIF
                  !
                  IF (SIZE(psol_dyn) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  champ(:,:) = psol_dyn(:,:)
                  !
              CASE DEFAULT
                  !
                  WRITE(*,*) 'startget_phys2d'
                  WRITE(*,*) 'No rule is present to extract variable', 
     .                 varname(:LEN_TRIM(varname)),' from any data set'
                  STOP
                  !
          END SELECT
          !
      ELSE
          !
          ! There are a few fields we might need if we need to interpolate 3D filed. Thus if they come through here we
          ! will catch them
          !
          SELECTCASE(varname)
              !
              CASE ('surfgeo')
                  !
                  IF ( .NOT.ALLOCATED(phis)) THEN
                      ALLOCATE(phis(iml,jml))
                  ENDIF
                  !
                  IF (SIZE(phis) .NE. SIZE(champ)) THEN
                      !
                      WRITE(*,*)
     .  'STARTVAR module has been initialized to the wrong size'
                      STOP
                      !
                  ENDIF
                  !
                  phis(:,:) = champ(:,:)
                  !
          END SELECT
          !
      ENDIF
    !
      END SUBROUTINE startget_phys2d
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE start_init_orog( iml, jml, lon_in, lat_in)
    !
      INTEGER, INTENT(in) :: iml, jml
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
    !
    !  LOCAL
    !
      REAL :: lev(1), date, dt
      INTEGER :: itau(1), fid
      INTEGER ::  llm_tmp, ttm_tmp
      INTEGER :: i, j
      INTEGER :: iret
      REAL, ALLOCATABLE :: relief_hi(:,:)
      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
      REAL, ALLOCATABLE :: tmp_var(:,:)
      INTEGER, ALLOCATABLE :: tmp_int(:,:)
    !
      CHARACTER*120 :: orogfname
      LOGICAL :: check=.TRUE.
    !
    !
      orogfname = 'Relief.nc'
    !
      IF ( check ) WRITE(*,*) 'Reading the high resolution orography'
    !
      CALL flininfo(orogfname,iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
    !
      ALLOCATE (lat_rel(iml_rel,jml_rel), stat=iret)
      ALLOCATE (lon_rel(iml_rel,jml_rel), stat=iret)
      ALLOCATE (relief_hi(iml_rel,jml_rel), stat=iret)
    !
      CALL flinopen(orogfname, .FALSE., iml_rel, jml_rel, 
     .llm_tmp, lon_rel, lat_rel, lev, ttm_tmp,
     .      itau, date, dt, fid)
    !
      CALL flinget(fid, 'RELIEF', iml_rel, jml_rel, llm_tmp, 
     . ttm_tmp, 1, 1, relief_hi)
    !
      CALL flinclo(fid)
    !
    !   In case we have a file which is in degrees we do the transformation
    !
      ALLOCATE(lon_rad(iml_rel))
      IF ( MAXVAL(lon_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lon_rad(:) = lon_rel(:,1) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lon_rad(:) = lon_rel(:,1) 
      ENDIF
      ALLOCATE(lat_rad(jml_rel))
      IF ( MAXVAL(lat_rel(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lat_rad(:) = lat_rel(1,:) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lat_rad(:) = lat_rel(1,:) 
      ENDIF
    !
    !
      IF ( check ) WRITE(*,*) 'Computes all the parameters needed',
     .' for the gravity wave drag code'
    !
    !    Allocate the data we need to put in the interpolated fields
    !
    !            RELIEF:  orographie moyenne
      ALLOCATE(relief(iml,jml))
    !            zphi :  orographie moyenne
      ALLOCATE(phis(iml,jml))
    !             zstd:  deviation standard de l'orographie sous-maille
      ALLOCATE(zstd(iml,jml))
    !             zsig:  pente de l'orographie sous-maille 
      ALLOCATE(zsig(iml,jml))
    !             zgam:  anisotropy de l'orographie sous maille
      ALLOCATE(zgam(iml,jml))
    !             zthe:  orientation de l'axe oriente dans la direction
    !                    de plus grande pente de l'orographie sous maille
      ALLOCATE(zthe(iml,jml))
    !             zpic:  hauteur pics de la SSO
      ALLOCATE(zpic(iml,jml))
    !             zval:  hauteur vallees de la SSO
      ALLOCATE(zval(iml,jml))
    !             masque : Masque terre ocean
      ALLOCATE(tmp_int(iml,jml))
      ALLOCATE(masque(iml,jml))
    !
      CALL grid_noro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
     $    iml-1, jml, lon_in, lat_in, 
    !     . phis, relief, zstd, zsig, zgam, zthe, zpic, zval, tmp_int)
    ! PB masque avec % terre mai 2000
     $    phis, relief, zstd, zsig, zgam, zthe, zpic, zval, masque)
      phis = phis * 9.81
    !
    !PB supression ligne suivant pour masque avec % terre
    !     masque(:,:) = FLOAT(tmp_int(:,:))
    !
    !  Compute surface roughness
    !
      IF ( check ) WRITE(*,*) 
     .'Compute surface roughness induced by the orography'
    !
      ALLOCATE(rugo(iml,jml))
      ALLOCATE(tmp_var(iml-1,jml))
    !
      CALL rugsoro(iml_rel, jml_rel, lon_rad, lat_rad, relief_hi,
     . iml-1, jml, lon_in, lat_in, tmp_var)
    !
      DO j = 1, jml
        DO i = 1, iml-1
          rugo(i,j) = tmp_var(i,j)
        ENDDO
        rugo(iml,j) = tmp_var(1,j)
      ENDDO
    !
    !   Build land-sea mask
    !
    !
      RETURN
    !
      END SUBROUTINE start_init_orog
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, 
     .lat_in, nbindex, champ, val_exp)
    !
      CHARACTER*(*), INTENT(in) :: varname
      INTEGER, INTENT(in) :: iml, jml, nbindex
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
      REAL, INTENT(inout) :: champ(nbindex)
      REAL, INTENT(in) :: val_exp
    !
    !
    !   This routine only works if the variable does not exist or is constant
    !
      IF ( MINVAL(champ(:)).EQ.MAXVAL(champ(:)) .AND. 
     .MINVAL(champ(:)).EQ.val_exp ) THEN
          SELECTCASE(varname)
            CASE ('tsol')
              IF ( .NOT.ALLOCATED(tsol)) THEN
                CALL start_init_phys( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(tsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*) 
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex, tsol, champ)
            CASE ('qsol')
              IF ( .NOT.ALLOCATED(qsol)) THEN
                CALL start_init_phys( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(qsol) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*) 
     . 'STARTVAR module has been initialized to the wrong size'
                STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex, qsol, champ)
            CASE ('psol')
              IF ( .NOT.ALLOCATED(psol_dyn)) THEN
                CALL start_init_dyn( iml, jml, lon_in, lat_in)
              ENDIF
              IF (SIZE(psol_dyn) .NE. SIZE(lon_in)*SIZE(lat_in)) THEN
                WRITE(*,*) 
     . 'STARTVAR module has been initialized to the wrong size'
                STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex, psol_dyn, champ)
    ! PB ajout pour masque terre mer fractionneiare    
          CASE ('zmasq')
              IF (.NOT. ALLOCATED(masque) ) THEN
                  CALL start_init_orog ( iml, jml,lon_in, lat_in)
              ENDIF
              IF ( SIZE(masque) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
             ENDIF
             CALL gr_dyn_fi(1, iml, jml, nbindex, masque, champ)
            CASE ('zmea')
              IF ( .NOT.ALLOCATED(relief)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(relief) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex, relief, champ)
            CASE ('zstd')
              IF ( .NOT.ALLOCATED(zstd)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zstd) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zstd, champ)
            CASE ('zsig')
              IF ( .NOT.ALLOCATED(zsig)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zsig) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zsig, champ)
            CASE ('zgam')
              IF ( .NOT.ALLOCATED(zgam)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zgam) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zgam, champ)
            CASE ('zthe')
              IF ( .NOT.ALLOCATED(zthe)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zthe) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zthe, champ)
            CASE ('zpic')
              IF ( .NOT.ALLOCATED(zpic)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zpic) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zpic, champ)
            CASE ('zval')
              IF ( .NOT.ALLOCATED(zval)) THEN
                CALL start_init_orog( iml, jml, lon_in, lat_in)
              ENDIF
              IF ( SIZE(zval) .NE. SIZE(lon_in)*SIZE(lat_in) ) THEN
                WRITE(*,*)
     . 'STARTVAR module has been initialized to the wrong size'
                 STOP
              ENDIF
              CALL gr_dyn_fi(1, iml, jml, nbindex,zval, champ)
             CASE ('rads')
                  champ(:) = 0.0
            CASE ('snow')
                  champ(:) = 0.0
            CASE ('deltat')
                  champ(:) = 0.0
            CASE ('rugmer')
                  champ(:) = 0.001
            CASE ('agsno')
                  champ(:) = 50.0
            CASE DEFAULT
              WRITE(*,*) 'startget_phys1d'
              WRITE(*,*) 'No rule is present to extract variable  ',
     . varname(:LEN_TRIM(varname)),' from any data set'
              STOP
          END SELECT
      ELSE
        !
        ! If we see tsol we catch it as we may need it for a 3D interpolation
        !
        SELECTCASE(varname)
          CASE ('tsol')
            IF ( .NOT.ALLOCATED(tsol)) THEN
              ALLOCATE(tsol(SIZE(lon_in),SIZE(lat_in) ))
            ENDIF
            CALL gr_fi_dyn(1, iml, jml, nbindex, champ, tsol)
        END SELECT
      ENDIF
      END SUBROUTINE startget_phys1d
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE start_init_phys( iml, jml, lon_in, lat_in)
    !
      INTEGER, INTENT(in) :: iml, jml
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
    !
    !  LOCAL
    !
      REAL :: lev(1), date, dt
      INTEGER :: itau(1)
      INTEGER ::  llm_tmp, ttm_tmp
      INTEGER :: i, j
    !
      CHARACTER*120 :: physfname
      LOGICAL :: check=.TRUE.
    !
      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:)
    !
      physfname = 'ECPHY.nc'
    !
      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
    !
      CALL flininfo(physfname, iml_phys, jml_phys, llm_tmp,
     . ttm_tmp, fid_phys)
    !
    !
      ALLOCATE (lat_phys(iml_phys,jml_phys))
      ALLOCATE (lon_phys(iml_phys,jml_phys))
    !
      CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, 
     . llm_tmp, lon_phys, lat_phys, lev, ttm_tmp, 
     . itau, date, dt, fid_phys)
    !
    ! Allocate the space we will need to get the data out of this file
    !
      ALLOCATE(var_ana(iml_phys, jml_phys))
    !
    !   In case we have a file which is in degrees we do the transformation
    !
      ALLOCATE(lon_rad(iml_phys))
      IF ( MAXVAL(lon_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lon_rad(:) = lon_phys(:,1) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lon_rad(:) = lon_phys(:,1) 
      ENDIF
      ALLOCATE(lat_rad(jml_phys))
      IF ( MAXVAL(lat_phys(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lat_rad(:) = lat_phys(1,:) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lat_rad(:) = lat_phys(1,:) 
      ENDIF
    !
    !   We get the two standard varibales
    !   Surface temperature
    !
      ALLOCATE(tsol(iml,jml))
      ALLOCATE(tmp_var(iml-1,jml))
    !
    !
      CALL flinget(fid_phys, 'ST', iml_phys, jml_phys, 
     .llm_tmp, ttm_tmp, 1, 1, var_ana)
      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
      CALL gr_int_dyn(tmp_var, tsol, iml-1, jml)
    !
    ! Soil moisture
    !
      ALLOCATE(qsol(iml,jml))
      CALL flinget(fid_phys, 'CDSW', iml_phys, jml_phys,
     . llm_tmp, ttm_tmp, 1, 1, var_ana)
      CALL grille_m(iml_phys, jml_phys, lon_rad, lat_rad,
     . var_ana, iml-1, jml, lon_in, lat_in, tmp_var)
      CALL gr_int_dyn(tmp_var, qsol, iml-1, jml)
    !
      CALL flinclo(fid_phys)
    !
      END SUBROUTINE start_init_phys
    !
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    !
      SUBROUTINE startget_dyn(varname, iml, jml, lon_in, lat_in,
     . lml, pls, workvar, champ, val_exp)
    !
    !   ARGUMENTS
    !
      CHARACTER*(*), INTENT(in) :: varname
      INTEGER, INTENT(in) :: iml, jml, lml
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
      REAL, INTENT(in) :: pls(iml, jml, lml)
      REAL, INTENT(in) :: workvar(iml, jml, lml)
      REAL, INTENT(inout) :: champ(iml, jml, lml)
      REAL, INTENT(in) :: val_exp
    !
    !    LOCAL
    !
      INTEGER :: il, ij, ii
      REAL :: xppn, xpps
    !
    !   C'est vraiment une galere de devoir rajouter tant de commons just pour avoir les aires.
    !   Il faudrait mettre une structure plus flexible et moins dangereuse.
    !
#include "dimensions.h"
#include "paramet.h"
#include "comgeom2.h"
#include "comconst.h"
    !
    !   This routine only works if the variable does not exist or is constant
    !
      IF ( MINVAL(champ(:,:,:)).EQ.MAXVAL(champ(:,:,:)) .AND.
     . MINVAL(champ(:,:,:)).EQ.val_exp ) THEN
        !
        SELECTCASE(varname)
          CASE ('u')
            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
              CALL start_init_dyn( iml, jml, lon_in, lat_in)
            ENDIF
            CALL start_inter_3d('U', iml, jml, lml, lon_in,
     .                           lat_in, pls, champ)
            DO il=1,lml
              DO ij=1,jml
                DO ii=1,iml-1
                  champ(ii,ij,il) = champ(ii,ij,il) * cu(ii,ij)
                ENDDO
                champ(iml,ij, il) = champ(1,ij, il)
              ENDDO
            ENDDO
          CASE ('v')
            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
              CALL start_init_dyn( iml, jml, lon_in, lat_in)
            ENDIF
            CALL start_inter_3d('V', iml, jml, lml, lon_in,
     .                          lat_in, pls, champ)
            DO il=1,lml
              DO ij=1,jml
                DO ii=1,iml-1
                  champ(ii,ij,il) = champ(ii,ij,il) * cv(ii,ij)
                ENDDO
                champ(iml,ij, il) = champ(1,ij, il)
              ENDDO
            ENDDO
          CASE ('t')
            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
              CALL start_init_dyn( iml, jml, lon_in, lat_in)
            ENDIF
            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
     .                           lat_in, pls, champ)
 
          CASE ('tpot')
            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
              CALL start_init_dyn( iml, jml, lon_in, lat_in)
            ENDIF
            CALL start_inter_3d('TEMP', iml, jml, lml, lon_in,
     .                           lat_in, pls, champ)
            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) )
     .                                    THEN
              DO il=1,lml
                DO ij=1,jml
                  DO ii=1,iml-1
                    champ(ii,ij,il) = champ(ii,ij,il) * cpp 
     .                                 / workvar(ii,ij,il)
                  ENDDO
                  champ(iml,ij,il) = champ(1,ij,il)
                ENDDO
              ENDDO
              DO il=1,lml
                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
                champ(:,1,il) = xppn
                champ(:,jml,il) = xpps
              ENDDO
            ELSE
              WRITE(*,*)'Could not compute potential temperature as the'
              WRITE(*,*)'Exner function is missing or constant.'
              STOP
            ENDIF
          CASE ('q')
            IF ( .NOT.ALLOCATED(psol_dyn)) THEN
              CALL start_init_dyn( iml, jml, lon_in, lat_in)
            ENDIF
            CALL start_inter_3d('R', iml, jml, lml, lon_in, lat_in,
     .                           pls, champ)
            IF ( MINVAL(workvar(:,:,:)) .NE. MAXVAL(workvar(:,:,:)) ) 
     .                                     THEN
              DO il=1,lml
                DO ij=1,jml
                  DO ii=1,iml-1
                    champ(ii,ij,il) = 0.01 * champ(ii,ij,il) *
     .                                       workvar(ii,ij,il)
                  ENDDO
                  champ(iml,ij,il) = champ(1,ij,il)
                ENDDO
              ENDDO
              WHERE ( champ .LT. 0.) champ = 1.0E-10
              DO il=1,lml
                xppn = SUM(aire(:,1)*champ(:,1,il))/apoln
                xpps = SUM(aire(:,jml)*champ(:,jml,il))/apols
                champ(:,1,il) = xppn
                champ(:,jml,il) = xpps
              ENDDO
            ELSE
              WRITE(*,*)'Could not compute specific humidity as the'
              WRITE(*,*)'saturated humidity is missing or constant.'
              STOP
            ENDIF
          CASE DEFAULT
            WRITE(*,*) 'startget_dyn'
            WRITE(*,*) 'No rule is present to extract variable  ',
     . varname(:LEN_TRIM(varname)),' from any data set'
            STOP
          END SELECT
      ENDIF
      END SUBROUTINE startget_dyn
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE start_init_dyn( iml, jml, lon_in, lat_in)
    !
      INTEGER, INTENT(in) :: iml, jml
      REAL, INTENT(in) :: lon_in(iml), lat_in(jml)
    !
    !  LOCAL
    !
      REAL :: lev(1), date, dt
      INTEGER :: itau(1)
      INTEGER :: i, j
      integer :: iret
    !
      CHARACTER*120 :: physfname
      LOGICAL :: check=.TRUE.
    !
      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
      REAL, ALLOCATABLE :: var_ana(:,:), tmp_var(:,:), z(:,:)
      REAL, ALLOCATABLE :: xppn(:), xpps(:)
      LOGICAL :: allo
    !
    !   Ce n'est pas tres pratique d'avoir a charger 3 include pour avoir la grille du modele
    !
#include "dimensions.h"
#include "paramet.h"
#include "comgeom2.h"
    !
      physfname = 'ECDYN.nc'
    !
      IF ( check ) WRITE(*,*) 'Opening the surface analysis'
    !
      CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn,
     .                            ttm_dyn, fid_dyn)
      IF ( check ) WRITE(*,*) 'Values read: ', iml_dyn, jml_dyn, 
     .                                         llm_dyn, ttm_dyn
    !
      ALLOCATE (lat_dyn(iml_dyn,jml_dyn), stat=iret)
      ALLOCATE (lon_dyn(iml_dyn,jml_dyn), stat=iret)
      ALLOCATE (lev_dyn(llm_dyn), stat=iret)
    !
      CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn,
     . lon_dyn, lat_dyn, lev_dyn, ttm_dyn, 
     . itau, date, dt, fid_dyn)
    !

      allo = allocated (var_ana)
      if (allo) then
        DEALLOCATE(var_ana, stat=iret)
      endif
      ALLOCATE(var_ana(iml_dyn, jml_dyn), stat=iret)

      allo = allocated (lon_rad)
      if (allo) then
        DEALLOCATE(lon_rad, stat=iret)
      endif
        ALLOCATE(lon_rad(iml_dyn), stat=iret)
       
      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lon_rad(:) = lon_dyn(:,1) 
      ENDIF
      ALLOCATE(lat_rad(jml_dyn))
      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lat_rad(:) = lat_dyn(1,:) 
      ENDIF
    !
      ALLOCATE(z(iml, jml))
      ALLOCATE(tmp_var(iml-1,jml))
    !
      CALL flinget(fid_dyn, 'Z', iml_dyn, jml_dyn, 0, ttm_dyn,
     .              1, 1, var_ana)
      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
     .               iml-1, jml, lon_in, lat_in, tmp_var)
      CALL gr_int_dyn(tmp_var, z, iml-1, jml)
    !
      ALLOCATE(psol_dyn(iml, jml))
    !
      CALL flinget(fid_dyn, 'SP', iml_dyn, jml_dyn, 0, ttm_dyn,
     .              1, 1, var_ana)
      CALL grille_m(iml_dyn, jml_dyn , lon_rad, lat_rad, var_ana,
     .               iml-1, jml, lon_in, lat_in, tmp_var)
      CALL gr_int_dyn(tmp_var, psol_dyn, iml-1, jml)
    !
      IF ( .NOT.ALLOCATED(tsol)) THEN
    !   These variables may have been allocated by the need to 
    !   create a start field for them or by the varibale
    !   coming out of the restart file. In case we dor have it we will initialize it.
    !
        CALL start_init_phys( iml, jml, lon_in, lat_in)
      ELSE
        IF ( SIZE(tsol) .NE. SIZE(psol_dyn) ) THEN
        WRITE(*,*) 'start_init_dyn :'
        WRITE(*,*) 'The temperature field we have does not ',
     .             'have the right size'
        STOP
      ENDIF
      ENDIF
      IF ( .NOT.ALLOCATED(phis)) THEN
            !
            !    These variables may have been allocated by the need to create a start field for them or by the varibale
            !     coming out of the restart file. In case we dor have it we will initialize it.
            !
          CALL start_init_orog( iml, jml, lon_in, lat_in)
            !
      ELSE
            !
          IF (SIZE(phis) .NE. SIZE(psol_dyn)) THEN
                !
              WRITE(*,*) 'start_init_dyn :'
              WRITE(*,*) 'The orography field we have does not ',
     .                   ' have the right size'
              STOP
          ENDIF
            !
      ENDIF
    !
    !     PSOL is computed in Pascals
    !
    !
      DO j = 1, jml
        DO i = 1, iml-1
          psol_dyn(i,j) = psol_dyn(i,j)*(1.0+(z(i,j)-phis(i,j)) 
     .                    /287.0/tsol(i,j))
        ENDDO
        psol_dyn(iml,j) = psol_dyn(1,j)
      ENDDO
    !
    !
      ALLOCATE(xppn(iml-1))
      ALLOCATE(xpps(iml-1)) 
    !
      DO  i   = 1, iml-1
        xppn(i) = aire( i,1) * psol_dyn( i,1)
        xpps(i) = aire( i,jml) * psol_dyn( i,jml)
      ENDDO
    !
      DO i   = 1, iml
        psol_dyn(i,1    )  = SUM(xppn)/apoln
        psol_dyn(i,jml)  = SUM(xpps)/apols
      ENDDO
    !
      RETURN
    !
      END SUBROUTINE start_init_dyn
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
      SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in,
     .                           lat_in, pls_in, var3d)
    !
    !    This subroutine gets a variables from a 3D file and does the interpolations needed
    !
    !
    !    ARGUMENTS
    !
      CHARACTER*(*) :: varname
      INTEGER :: iml, jml, lml
      REAL :: lon_in(iml), lat_in(jml), pls_in(iml, jml, lml)
      REAL :: var3d(iml, jml, lml)
    !
    !  LOCAL
    !
      INTEGER :: ii, ij, il
      REAL :: bx, by
      REAL, ALLOCATABLE :: lon_rad(:), lat_rad(:)
      REAL, ALLOCATABLE :: var_tmp2d(:,:), var_tmp3d(:,:,:)
      REAL, ALLOCATABLE :: ax(:), ay(:), yder(:)
      INTEGER, ALLOCATABLE :: lind(:)
    !
      LOGICAL :: check = .TRUE.
    !
      IF ( .NOT. ALLOCATED(var_ana3d)) THEN
          ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
      ENDIF
    !
    !
      IF ( check) WRITE(*,*) 'Going into flinget to extract the 3D ',
     .  ' field.', fid_dyn
      IF ( check) WRITE(*,*) fid_dyn, varname, iml_dyn, jml_dyn,
     .                        llm_dyn,ttm_dyn
    !
      CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, 
     . ttm_dyn, 1, 1, var_ana3d)
    !
      IF ( check) WRITE(*,*) 'Allocating space for the interpolation',
     . iml, jml, llm_dyn
    !
      ALLOCATE(lon_rad(iml_dyn))
      IF ( MAXVAL(lon_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lon_rad(:) = lon_dyn(:,1) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lon_rad(:) = lon_dyn(:,1) 
      ENDIF
      ALLOCATE(lat_rad(jml_dyn))
      IF ( MAXVAL(lat_dyn(:,:)) .GT. 2.0 * ASIN(1.0) ) THEN
          lat_rad(:) = lat_dyn(1,:) * 2.0 * ASIN(1.0) / 180.0
      ELSE
          lat_rad(:) = lat_dyn(1,:) 
      ENDIF
    !
      ALLOCATE(var_tmp2d(iml-1, jml))
      ALLOCATE(var_tmp3d(iml, jml, llm_dyn))
      ALLOCATE(ax(llm_dyn))
      ALLOCATE(ay(llm_dyn))
      ALLOCATE(yder(llm_dyn))
      ALLOCATE(lind(llm_dyn))
    !
      DO il=1,llm_dyn
        !
        CALL grille_m(iml_dyn, jml_dyn, lon_rad, lat_rad, 
     .var_ana3d(:,:,il), iml-1, jml, lon_in, lat_in, var_tmp2d)
        !
        CALL gr_int_dyn(var_tmp2d, var_tmp3d(:,:,il), iml-1, jml)
        !
       ENDDO
       !
       !  IF needed we return the vertical axis. The spline interpolation
       !  Requires the coordinate to be in increasing order.
    !
      IF ( lev_dyn(1) .LT. lev_dyn(llm_dyn)) THEN
          DO il=1,llm_dyn
            lind(il) = il
          ENDDO        
      ELSE
          DO il=1,llm_dyn
            lind(il) = llm_dyn-il+1
          ENDDO
      ENDIF
    !
      DO ij=1,jml
        DO ii=1,iml-1
          !
          ax(:) = lev_dyn(lind(:)) * 100
          ay(:) = var_tmp3d(ii, ij, lind(:))
          !
          CALL SPLINE(ax, ay, llm_dyn, 1.e30, 1.e30, yder)
          !
          DO il=1,lml
            bx = pls_in(ii, ij, il)
            CALL SPLINT(ax, ay, yder, llm_dyn, bx, by)
            var3d(ii, ij, il) = by
          ENDDO
          !
        ENDDO
        var3d(iml, ij, :) = var3d(1, ij, :) 
      ENDDO

      DEALLOCATE(lon_rad)
      DEALLOCATE(lat_rad)
      DEALLOCATE(var_tmp2d)
      DEALLOCATE(var_tmp3d)
      DEALLOCATE(ax)
      DEALLOCATE(ay)
      DEALLOCATE(yder)
      DEALLOCATE(lind)

    !
      RETURN
    !
      END SUBROUTINE start_inter_3d
    !
      END MODULE startvar