Index: trunk/LMDZ.COMMON/libf/dyn3d/grid_noro.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/grid_noro.F	(revision 780)
+++ 	(revision )
@@ -1,521 +1,0 @@
-!
-! $Id: grid_noro.F 1442 2010-10-18 08:31:31Z jghattas $
-!
-c
-c
-      SUBROUTINE grid_noro(imdep, jmdep, xdata, ydata, zdata,
-     .             imar, jmar, x, y,
-     .             zphi,zmea,zstd,zsig,zgam,zthe,
-     .             zpic,zval,mask)
-c=======================================================================
-c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
-c
-c      Compute the Parameters of the SSO scheme as described in
-c      LOTT & MILLER (1997) and LOTT(1999).
-c      Target points are on a rectangular grid:
-c      iim+1 latitudes including North and South Poles;
-c      jjm+1 longitudes, with periodicity jjm+1=1.
-c      aux poles.  At the poles the fields value is repeated
-c      jjm+1 time.
-c      The parameters a,b,c,d represent the limite of the target
-c      gridpoint region. The means over this region are calculated
-c      from USN data, ponderated by a weight proportional to the 
-c      surface occupated by the data inside the model gridpoint area.
-c      In most circumstances, this weight is the ratio between the
-c      surface of the USN gridpoint area and the surface of the
-c      model gridpoint area. 
-c
-c           (c)
-c        ----d-----
-c        | . . . .|
-c        |        |
-c     (b)a . * . .b(a)
-c        |        |
-c        | . . . .|
-c        ----c-----
-c           (d)
-C=======================================================================
-c INPUT:
-c        imdep, jmdep: dimensions X and Y input field
-c        xdata, ydata: coordinates X and Y input field
-c        zdata: Input field
-c        In this version it is assumed that the entry data come from
-c        the USNavy dataset: imdep=iusn=2160, jmdep=jusn=1080.
-c OUTPUT:
-c        imar, jmar: dimensions X and Y Output field
-c        x, y: ccordinates  X and Y Output field.
-c             zmea:  Mean orographie   
-c             zstd:  Standard deviation
-c             zsig:  Slope
-c             zgam:  Anisotropy
-c             zthe:  Orientation of the small axis
-c             zpic:  Maximum altitude
-c             zval:  Minimum altitude
-C=======================================================================
-
-      IMPLICIT INTEGER (I,J)
-      IMPLICIT REAL(X,Z) 
-      
-	  parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
-#include "dimensions.h"
-	  REAL xusn(iusn+2*iext),yusn(jusn+2)	
-      REAL zusn(iusn+2*iext,jusn+2)
-
-      INTEGER imdep, jmdep
-      REAL xdata(imdep),ydata(jmdep) 
-      REAL zdata(imdep,jmdep)
-c
-      INTEGER imar, jmar
-  
-C INTERMEDIATE FIELDS  (CORRELATIONS OF OROGRAPHY GRADIENT)
-
-      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
-      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
-      REAL weight(iim+1,jjm+1)
-
-C CORRELATIONS OF USN OROGRAPHY GRADIENTS
-
-      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
-      REAL zxtzyusn(iusn+2*iext,jusn+2)
-      REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
-      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
-      REAL zmea0(imar+1,jmar) ! GK211005 (CG)
-      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
-      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
-cx$$ PB     integer mask(imar+1,jmar)
-      real mask(imar+1,jmar), mask_tmp(imar+1,jmar)
-      real num_tot(2200,1100),num_lan(2200,1100)
-c
-      REAL a(2200),b(2200),c(1100),d(1100)
-      logical masque_lu
-c
-      print *,' parametres de l orographie a l echelle sous maille' 
-      xpi=acos(-1.)
-      rad    = 6 371 229.
-      zdeltay=2.*xpi/REAL(jusn)*rad
-c
-c utilise-t'on un masque lu?
-c
-      masque_lu = .true.
-      if (maxval(mask) == -99999 .and. minval(mask) == -99999) then
-        masque_lu= .false.
-        masque = 0.0
-      endif
-      write(*,*)'Masque lu', masque_lu
-c
-c  quelques tests de dimensions:
-c    
-c
-      if(iim.ne.imar) STOP 'Problem dim. x'
-      if(jjm.ne.jmar-1) STOP 'Problem dim. y'
-      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
-         PRINT*, 'imar or jmar too big', imar, jmar
-         CALL ABORT
-      ENDIF
-
-      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
-         print *,' imdep or jmdep bad dimensions:',imdep,jmdep
-         call abort
-      ENDIF
-
-      IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
-        print *,' imar or jmar bad dimensions:',imar,jmar
-        call abort
-      ENDIF
-
-
-c      print *,'xdata:',xdata
-c      print *,'ydata:',ydata
-c      print *,'x:',x
-c      print *,'y:',y
-c
-C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
-C  BOUNDARIES:
-c
-      DO j=1,jusn
-        yusn(j+1)=ydata(j)
-      DO i=1,iusn
-        zusn(i+iext,j+1)=zdata(i,j)
-        xusn(i+iext)=xdata(i)
-      ENDDO
-      DO i=1,iext
-        zusn(i,j+1)=zdata(iusn-iext+i,j)
-        xusn(i)=xdata(iusn-iext+i)-2.*xpi
-        zusn(iusn+iext+i,j+1)=zdata(i,j)
-        xusn(iusn+iext+i)=xdata(i)+2.*xpi
-      ENDDO
-      ENDDO
-
-        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
-        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
-       DO i=1,iusn/2+iext
-        zusn(i,1)=zusn(i+iusn/2,2)
-        zusn(i+iusn/2+iext,1)=zusn(i,2)
-        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
-        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
-       ENDDO
-c  
-c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
-C     ( REGULAR GRID)
-c
-      a(1) = x(1) - (x(2)-x(1))/2.0
-      b(1) = (x(1)+x(2))/2.0
-      DO i = 2, imar
-         a(i) = b(i-1)
-         b(i) = (x(i)+x(i+1))/2.0
-      ENDDO
-      a(imar+1) = b(imar)
-      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0
-
-      c(1) = y(1) - (y(2)-y(1))/2.0
-      d(1) = (y(1)+y(2))/2.0
-      DO j = 2, jmar-1
-         c(j) = d(j-1)
-         d(j) = (y(j)+y(j+1))/2.0
-      ENDDO
-      c(jmar) = d(jmar-1)
-      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
-c
-c  initialisations:
-c
-      DO i = 1, imar+1
-      DO j = 1, jmar
-         weight(i,j) = 0.0
-         zxtzx(i,j)  = 0.0
-         zytzy(i,j)  = 0.0
-         zxtzy(i,j)  = 0.0
-         ztz(i,j)    = 0.0
-         zmea(i,j)   = 0.0
-         zpic(i,j)  =-1.E+10
-         zval(i,j)  = 1.E+10
-      ENDDO
-      ENDDO
-c
-c  COMPUTE SLOPES CORRELATIONS ON USN GRID
-c
-         DO j = 1,jusn+2 
-         DO i = 1, iusn+2*iext
-            zytzyusn(i,j)=0.0
-            zxtzxusn(i,j)=0.0
-            zxtzyusn(i,j)=0.0
-         ENDDO
-         ENDDO
-
-
-         DO j = 2,jusn+1 
-            zdeltax=zdeltay*cos(yusn(j))
-         DO i = 2, iusn+2*iext-1
-            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
-            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
-            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
-     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
-         ENDDO
-         ENDDO
-c
-c  SUMMATION OVER GRIDPOINT AREA
-c 
-      zleny=xpi/REAL(jusn)*rad
-      xincr=xpi/2./REAL(jusn)
-       DO ii = 1, imar+1
-       DO jj = 1, jmar
-       num_tot(ii,jj)=0.
-       num_lan(ii,jj)=0.
-c        PRINT *,' iteration ii jj:',ii,jj
-         DO j = 2,jusn+1 
-c         DO j = 3,jusn 
-            zlenx=zleny*cos(yusn(j))
-            zdeltax=zdeltay*cos(yusn(j))
-            zbordnor=(c(jj)-yusn(j)+xincr)*rad
-            zbordsud=(yusn(j)-d(jj)+xincr)*rad
-            weighy=AMAX1(0.,
-     *             amin1(zbordnor,zbordsud,zleny))
-         IF(weighy.ne.0)THEN
-         DO i = 2, iusn+2*iext-1
-            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
-            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
-            weighx=AMAX1(0.,
-     *             amin1(zbordest,zbordoue,zlenx))
-            IF(weighx.ne.0)THEN
-            num_tot(ii,jj)=num_tot(ii,jj)+1.0
-            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
-            weight(ii,jj)=weight(ii,jj)+weighx*weighy
-            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
-            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
-            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
-            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
-c mean
-            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
-c peacks
-            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
-c valleys
-            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
-            ENDIF
-         ENDDO
-         ENDIF
-         ENDDO
-       ENDDO
-       ENDDO
-c
-c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
-C  LOTT (1999) SSO SCHEME.
-c
-      zllmmea=0.
-      zllmstd=0.
-      zllmsig=0.
-      zllmgam=0.
-      zllmpic=0.
-      zllmval=0.
-      zllmthe=0.
-      zminthe=0.
-c     print 100,' '
-c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH') 
-       DO ii = 1, imar+1
-       DO jj = 1, jmar
-         IF (weight(ii,jj) .NE. 0.0) THEN
-c  Mask
-cx$$           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
-cx$$             mask(ii,jj)=1
-cx$$           else
-cx$$             mask(ii,jj)=0
-cx$$           ENDIF
-             if (.not. masque_lu) then
-               mask(ii,jj) = num_lan(ii,jj)/num_tot(ii,jj)
-             endif
-c  Mean Orography:
-           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
-           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
-           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
-           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
-           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
-c  Standard deviation:
-           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
-         ELSE
-            PRINT*, 'probleme,ii,jj=', ii,jj
-         ENDIF
-       ENDDO
-       ENDDO
-
-C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
-
-       DO ii = 1, imar+1
-         zxtzx(ii,1)=zxtzx(ii,2)
-         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
-         zxtzy(ii,1)=zxtzy(ii,2)
-         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
-         zytzy(ii,1)=zytzy(ii,2)
-         zytzy(ii,jmar)=zytzy(ii,jmar-1)
-       ENDDO
-
-C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
-
-C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
-
-       zmea0(:,:) = zmea(:,:) ! GK211005 (CG) on sauvegarde la topo non lissee
-       CALL MVA9(zmea,iim+1,jjm+1)
-       CALL MVA9(zstd,iim+1,jjm+1)
-       CALL MVA9(zpic,iim+1,jjm+1)
-       CALL MVA9(zval,iim+1,jjm+1)
-       CALL MVA9(zxtzx,iim+1,jjm+1)
-       CALL MVA9(zxtzy,iim+1,jjm+1) 
-       CALL MVA9(zytzy,iim+1,jjm+1)
-Cx$$   Masque prenant en compte maximum de terre
-Cx$$  On seuil a 10% de terre de terre car en dessous les parametres de surface n'on
-Cx$$ pas de sens (PB)
-       mask_tmp= 0.0
-       WHERE(mask .GE. 0.1) mask_tmp = 1.
-
-       DO ii = 1, imar
-       DO jj = 1, jmar
-         IF (weight(ii,jj) .NE. 0.0) THEN
-c  Coefficients K, L et M:
-           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
-           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
-           xm=zxtzy(ii,jj)
-           xp=xk-sqrt(xl**2+xm**2)
-           xq=xk+sqrt(xl**2+xm**2)
-           xw=1.e-8
-           if(xp.le.xw) xp=0.
-           if(xq.le.xw) xq=xw
-           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
-c slope: 
-cx$$           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
-cx$$c isotropy:
-cx$$           zgam(ii,jj)=xp/xq*mask(ii,jj)
-cx$$c angle theta:
-cx$$           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
-cx$$           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
-cx$$           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
-cx$$           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
-cx$$           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
-cx$$           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
-Cx$* PB modif pour maque de terre fractionnaire
-c slope: 
-           zsig(ii,jj)=sqrt(xq)*mask_tmp(ii,jj)
-c isotropy:
-           zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
-c angle theta:
-           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask_tmp(ii,jj)
-           ! GK211005 (CG) ne pas forcement lisser la topo
-           ! zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
-           zphi(ii,jj)=zmea0(ii,jj)*mask_tmp(ii,jj)
-           !
-           zmea(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
-           zpic(ii,jj)=zpic(ii,jj)*mask_tmp(ii,jj)
-           zval(ii,jj)=zval(ii,jj)*mask_tmp(ii,jj)
-           zstd(ii,jj)=zstd(ii,jj)*mask_tmp(ii,jj)
-c          print 101,ii,jj,
-c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
-c    *           zthe(ii,jj)
-c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
-         ELSE
-c           PRINT*, 'probleme,ii,jj=', ii,jj
-         ENDIF
-      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
-      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
-      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
-      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
-      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
-      zminthe=amin1(zthe(ii,jj),zminthe)
-      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
-      zllmval=AMAX1(zval(ii,jj),zllmval)
-       ENDDO
-       ENDDO
-      print *,'  MEAN ORO:',zllmmea
-      print *,'  ST. DEV.:',zllmstd
-      print *,'  PENTE:',zllmsig
-      print *,' ANISOTROP:',zllmgam
-      print *,'  ANGLE:',zminthe,zllmthe	
-      print *,'  pic:',zllmpic
-      print *,'  val:',zllmval
-      
-C
-c gamma and theta a 1. and 0. at poles
-c
-      DO jj=1,jmar
-      zmea(imar+1,jj)=zmea(1,jj)
-      zphi(imar+1,jj)=zphi(1,jj)
-      zpic(imar+1,jj)=zpic(1,jj)
-      zval(imar+1,jj)=zval(1,jj)
-      zstd(imar+1,jj)=zstd(1,jj)
-      zsig(imar+1,jj)=zsig(1,jj)
-      zgam(imar+1,jj)=zgam(1,jj)
-      zthe(imar+1,jj)=zthe(1,jj)
-      ENDDO
-
-
-      zmeanor=0.0
-      zmeasud=0.0
-      zstdnor=0.0
-      zstdsud=0.0
-      zsignor=0.0
-      zsigsud=0.0
-      zweinor=0.0
-      zweisud=0.0
-      zpicnor=0.0
-      zpicsud=0.0                                   
-      zvalnor=0.0
-      zvalsud=0.0 
-
-      DO ii=1,imar
-      zweinor=zweinor+              weight(ii,   1)
-      zweisud=zweisud+              weight(ii,jmar)
-      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
-      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
-      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
-      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
-      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
-      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
-      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
-      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
-      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
-      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
-      ENDDO
-
-      DO ii=1,imar+1
-      zmea(ii,   1)=zmeanor/zweinor
-      zmea(ii,jmar)=zmeasud/zweisud
-      zphi(ii,   1)=zmeanor/zweinor
-      zphi(ii,jmar)=zmeasud/zweisud
-      zpic(ii,   1)=zpicnor/zweinor
-      zpic(ii,jmar)=zpicsud/zweisud
-      zval(ii,   1)=zvalnor/zweinor
-      zval(ii,jmar)=zvalsud/zweisud
-      zstd(ii,   1)=zstdnor/zweinor
-      zstd(ii,jmar)=zstdsud/zweisud
-      zsig(ii,   1)=zsignor/zweinor
-      zsig(ii,jmar)=zsigsud/zweisud
-      zgam(ii,   1)=1.
-      zgam(ii,jmar)=1.
-      zthe(ii,   1)=0.
-      zthe(ii,jmar)=0.
-      ENDDO
-
-      RETURN
-      END
-
-      SUBROUTINE MVA9(X,IMAR,JMAR)
-
-C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
-
-      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
-      real WEIGHTpb(-1:1,-1:1)
-
-
-      SUM=0.
-      DO IS=-1,1
-        DO JS=-1,1
-          WEIGHTpb(IS,JS)=1./REAL((1+IS**2)*(1+JS**2))
-          SUM=SUM+WEIGHTpb(IS,JS)
-        ENDDO
-      ENDDO
-      
-c     WRITE(*,*) 'MVA9 ', IMAR, JMAR
-c     WRITE(*,*) 'MVA9 ', WEIGHTpb
-c     WRITE(*,*) 'MVA9 SUM ', SUM
-      DO IS=-1,1
-        DO JS=-1,1
-          WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
-        ENDDO
-      ENDDO
-
-      DO J=2,JMAR-1
-        DO I=2,IMAR-1
-          XF(I,J)=0.
-          DO IS=-1,1
-            DO JS=-1,1
-              XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
-            ENDDO
-          ENDDO
-        ENDDO
-      ENDDO
-
-      DO J=2,JMAR-1
-        XF(1,J)=0.
-        IS=IMAR-1
-        DO JS=-1,1 
-          XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
-        ENDDO
-        DO IS=0,1 
-          DO JS=-1,1 
-            XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
-          ENDDO
-        ENDDO
-        XF(IMAR,J)=XF(1,J)
-      ENDDO
-
-      DO I=1,IMAR
-        XF(I,1)=XF(I,2)
-        XF(I,JMAR)=XF(I,JMAR-1)
-      ENDDO
-
-      DO I=1,IMAR
-        DO J=1,JMAR
-          X(I,J)=XF(I,J)
-        ENDDO
-      ENDDO
-
-      RETURN
-      END
-
-
-
Index: trunk/LMDZ.COMMON/libf/dyn3d/startvar.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3d/startvar.F90	(revision 780)
+++ 	(revision )
@@ -1,783 +1,0 @@
-!
-! $Id: startvar.F90 1425 2010-09-02 13:45:23Z lguez $
-!
-!*******************************************************************************
-!
-MODULE startvar
-!
-!*******************************************************************************
-!
-!-------------------------------------------------------------------------------
-! Purpose: Access data from the database of atmospheric to initialize the model.
-!-------------------------------------------------------------------------------
-! Comments:
-!
-!    *  This module is designed to work for Earth (and with ioipsl)
-!
-!    *  There are three ways to acces data, depending 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(...)
-!
-!  - A 1D variable on the physical grid :
-!    CALL startget_phys1d((varname, iml, jml,  lon_in,  lat_in,  nbindex,              &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
-!
-!  - A 2D variable on the dynamical grid :
-!    CALL startget_phys2d(varname, iml, jml,  lon_in,  lat_in,                        &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )             
-!
-!  - A 3D variable on the dynamical grid :
-!    CALL startget_dyn((varname, iml, jml,  lon_in,  lat_in,  lml, pls, workvar,    &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
-!
-!    *  Data needs to be in NetCDF format
-!
-!    *  Variables should have the following names in the files:
-!            '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
-!
-!   *   There is a big mess with the longitude size. 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 smaller fields when needed. This needs to be cleared once and for
-!  all in LMDZ. A convention is required.
-!-------------------------------------------------------------------------------
-#ifdef CPP_EARTH
-  USE ioipsl
-  IMPLICIT NONE
-
-  PRIVATE
-  PUBLIC startget_phys2d, startget_phys1d, startget_dyn
-!  INTERFACE startget
-!    MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn
-!  END INTERFACE
-
-  REAL,    SAVE :: deg2rad,  pi
-  INTEGER, SAVE ::           iml_rel,  jml_rel
-  INTEGER, SAVE :: fid_phys, iml_phys, jml_phys
-  INTEGER, SAVE :: fid_dyn,  iml_dyn,  jml_dyn,  llm_dyn,  ttm_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel
-  REAL, DIMENSION(:),     ALLOCATABLE, TARGET, SAVE :: levdyn_ini
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: psol_dyn
-  REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d
-
-   CONTAINS
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ,  &
-                           val_exp ,jml2, lon_in2, lat_in2, ibar)
-!
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),         INTENT(IN)    :: varname
-  INTEGER,                  INTENT(IN)    :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)    :: lat_in
-  INTEGER,                  INTENT(IN)    :: nbindex
-  REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ
-  REAL,                     INTENT(IN)    :: val_exp
-  INTEGER,                  INTENT(IN)    :: jml2
-  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in2
-  REAL, DIMENSION(jml2),    INTENT(IN)    :: lat_in2
-  LOGICAL,                  INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(:,:), POINTER :: v2d
-!-------------------------------------------------------------------------------
-  v2d=>NULL()
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
-    SELECT CASE(varname)
-      CASE('tsol')
-        IF(.NOT.ALLOCATED(tsol))                                               &
-         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('qsol')
-        IF(.NOT.ALLOCATED(qsol))                                               &
-         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('psol')
-        IF(.NOT.ALLOCATED(psol_dyn))                                           &
-         CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval')
-        IF(.NOT.ALLOCATED(relief))                                             &
-         CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('rads','snow','tslab','seaice','rugmer','agsno')
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_phys1d'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                     //' from any data set'; STOP
-    END SELECT
-
-!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
-    SELECT CASE(varname)
-      CASE('rads','snow','tslab','seaice');  champ=0.0
-      CASE('rugmer');                        champ(:)=0.001
-      CASE('agsno');                         champ(:)=50.0
-      CASE DEFAULT
-        SELECT CASE(varname)
-          CASE('tsol'); v2d=>tsol
-          CASE('qsol'); v2d=>qsol
-          CASE('psol'); v2d=>psol_dyn
-          CASE('zmea'); v2d=>relief
-          CASE('zstd'); v2d=>zstd
-          CASE('zsig'); v2d=>zsig
-          CASE('zgam'); v2d=>zgam
-          CASE('zthe'); v2d=>zthe
-          CASE('zpic'); v2d=>zpic
-          CASE('zval'); v2d=>zval
-        END SELECT
-        IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN
-         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
-         STOP
-        END IF
-        CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ)
-    END SELECT
-
-  ELSE
-
-!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
-    SELECT CASE(varname)
-      CASE('tsol')
-        IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml))
-        CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol)
-    END SELECT
-
-  END IF
-
-END SUBROUTINE  startget_phys1d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp,  &
-                           jml2, lon_in2, lat_in2 , ibar, msk)
-!
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),         INTENT(IN)           :: varname
-  INTEGER,                  INTENT(IN)           :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
-  REAL, DIMENSION(iml,jml), INTENT(INOUT)        :: champ
-  REAL,                     INTENT(IN)           :: val_exp
-  INTEGER,                  INTENT(IN)           :: jml2
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in2
-  REAL, DIMENSION(jml2),    INTENT(IN)           :: lat_in2
-  LOGICAL,                  INTENT(IN)           :: ibar
-  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(:,:), POINTER :: v2d=>NULL()
-  LOGICAL                       :: lrelief1, lrelief2
-!-------------------------------------------------------------------------------
-  v2d=>NULL()
-  lrelief1=(.NOT.ALLOCATED(relief).AND.     PRESENT(msk))
-  lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk))
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
-    SELECT CASE(varname)
-      CASE('psol')
-        IF(.NOT.ALLOCATED(psol_dyn))                                           &
-          CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('relief')
-        IF(lrelief1)             CALL start_init_orog(iml,jml,lon_in,lat_in,msk)
-        IF(lrelief2)             CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('surfgeo')
-        IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('rugosite','masque')
-        IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_phys2d'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                     //' from any data set'; STOP
-    END SELECT
-
-!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
-    SELECT CASE(varname)
-      CASE('psol');     v2d=>psol_dyn
-      CASE('relief');   v2d=>relief
-      CASE('rugosite'); v2d=>rugo
-      CASE('masque');   v2d=>masque
-      CASE('surfgeo');  v2d=>phis
-    END SELECT
-    IF(SIZE(champ)/=SIZE(v2d)) THEN
-      WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size'
-      STOP
-    END IF
-
-    champ(:,:)=v2d(:,:)
-
-  ELSE
-
-!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
-    SELECT CASE(varname)
-      CASE ('surfgeo')
-        IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml))
-        IF(SIZE(phis)/=SIZE(champ)) THEN
-         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
-         STOP
-        END IF
-        phis(:,:)=champ(:,:)
-    END SELECT
-
-  END IF
-
-END SUBROUTINE startget_phys2d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_dyn(varname,  lon_in,  lat_in, pls,workvar,&
-                     champ, val_exp, lon_in2, lat_in2, ibar)
-
-      use assert_eq_m, only: assert_eq
-
-
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN)    :: varname
-  REAL, INTENT(IN)    :: lon_in(:) ! dim(iml)
-  REAL, INTENT(IN)    :: lat_in(:) ! dim(jml)
-  REAL, INTENT(IN)    :: pls(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(IN)    :: workvar(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(INOUT) :: champ(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(IN)    :: val_exp
-  REAL, INTENT(IN)    :: lon_in2(:) ! dim(iml)
-  REAL, INTENT(IN)    :: lat_in2(:) ! dim(jml2)
-  LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-#include "dimensions.h"
-#include "comconst.h"
-#include "paramet.h"
-#include "comgeom2.h"
-  INTEGER    :: iml, jml
-  INTEGER    :: lml
-  INTEGER    :: jml2
-  REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
-  CHARACTER(LEN=10) :: vname
-  INTEGER :: il
-  REAL    :: xppn, xpps
-!-------------------------------------------------------------------------------
-  NULLIFY(v3d)
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-      iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), &
-     &     size(champ, 1), size(lon_in2)/), "startget_dyn iml")
-      jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2),   &
-     &     size(champ, 2), "startget_dyn jml")
-      lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), &
-     &     "startget_dyn lml")
-      jml2 = size(lat_in2)
-
-!--- READING UNALLOCATED FILES
-    IF(.NOT.ALLOCATED(psol_dyn))                                               &
-      CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-
-!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
-    SELECT CASE(varname)
-      CASE('u');        vname='U'
-      CASE('v');        vname='V'
-      CASE('t','tpot'); vname='TEMP'
-      CASE('q');        vname='R'
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_dyn'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                //' from any data set'; STOP
-    END SELECT
-    CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2,      &
-                        lon_in2, lat_in2,  pls, champ,ibar )
-
-!--- COMPUTING THE REQUIRED FILED
-    SELECT CASE(varname)
-      CASE('u')                                            !--- Eastward wind
-        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
-        champ(iml,:,:)=champ(1,:,:)
-
-      CASE('v')                                            !--- Northward wind
-        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
-        champ(iml,:,:)=champ(1,:,:)
-
-      CASE('tpot')                                         !--- Temperature
-        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
-          champ=champ*cpp/workvar
-          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
-          END DO
-        ELSE
-          WRITE(lunout,*)'Could not compute potential temperature as the'
-          WRITE(lunout,*)'Exner function is missing or constant.'; STOP
-        END IF
-
-      CASE('q')                                            !--- Relat. humidity
-        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
-          champ=0.01*champ*workvar
-          WHERE(champ<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
-          END DO
-        ELSE
-          WRITE(lunout,*)'Could not compute specific humidity as the'
-          WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
-        END IF
-
-    END SELECT
-
-  END IF
-
-END SUBROUTINE startget_dyn
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,                  INTENT(IN)           :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
-  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: orofname
-  LOGICAL               :: check=.TRUE.
-  REAL,    DIMENSION(1) :: lev
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: fid, llm_tmp, ttm_tmp
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-!-------------------------------------------------------------------------------
-  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-
-  orofname = 'Relief.nc'; title='RELIEF'
-  IF(check) WRITE(lunout,*)'Reading the high resolution orography'
-  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
-
-  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
-  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
-                lev, ttm_tmp, itau, date, dt, fid)
-  ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
-  CALL flinclo(fid)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
-  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  relief_hi, .FALSE.)
-  DEALLOCATE(lon_ini,lat_ini)
-
-!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
-     &g code'
-
-  ALLOCATE(phis(iml,jml))      ! Geopotentiel au sol
-  ALLOCATE(zstd(iml,jml))      ! Deviation standard de l'orographie sous-maille
-  ALLOCATE(zsig(iml,jml))      ! Pente de l'orographie sous-maille 
-  ALLOCATE(zgam(iml,jml))      ! Anisotropie de l'orographie sous maille
-  ALLOCATE(zthe(iml,jml))      ! Orientation axe +grande pente d'oro sous maille
-  ALLOCATE(zpic(iml,jml))      ! Hauteur pics de la SSO
-  ALLOCATE(zval(iml,jml))      ! Hauteur vallees de la SSO
-  ALLOCATE(relief(iml,jml))    ! Orographie moyenne
-  ALLOCATE(masque(iml,jml))    ! Masque terre ocean
-  masque = -99999.
-  IF(PRESENT(masque_lu)) masque=masque_lu
-
-  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, masque)
-  phis = phis * 9.81
-
-!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
-  IF(check) WRITE(lunout,*)'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)
-  rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
-  DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
-  RETURN
-
-END SUBROUTINE start_init_orog
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,               INTENT(IN) :: iml, jml
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
-  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
-  INTEGER,               INTENT(IN) :: jml2
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
-  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
-  LOGICAL,               INTENT(IN) :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: physfname
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: llm_tmp, ttm_tmp
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),   ALLOCATABLE :: levphys_ini
-!-------------------------------------------------------------------------------
-  physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(check) WRITE(lunout,*)'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))
-  ALLOCATE(levphys_ini(llm_tmp))
-  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys,     &
-                lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
-  DEALLOCATE(levphys_ini)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
-  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
-
-  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
-
-!--- SURFACE TEMPERATURE
-  title='ST'
-  ALLOCATE(tsol(iml,jml))
-  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
-  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
-                  var_ana , ibar  )
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
-      lon_in,   lat_in,   lon_in2, lat_in2, tsol)
-
-!--- SOIL MOISTURE
-  title='CDSW'
-  ALLOCATE(qsol(iml,jml))
-  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
-  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
-                  var_ana, ibar  )
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
-      lon_in,   lat_in,   lon_in2, lat_in2, qsol)
-
-  CALL flinclo(fid_phys)
-
-  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
-
-END SUBROUTINE start_init_phys
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,               INTENT(IN) :: iml, jml
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
-  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
-  INTEGER,               INTENT(IN) :: jml2
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
-  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
-  LOGICAL,               INTENT(IN) :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: physfname
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: i, j
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana, z
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),   ALLOCATABLE :: xppn, xpps
-!-------------------------------------------------------------------------------
-
-!--- KINETIC ENERGY
-  physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(check) WRITE(lunout,*) 'Opening the surface analysis'
-  CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
-  IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-
-  ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
-  ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
-  ALLOCATE(levdyn_ini(llm_dyn))
-  CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
-                levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
-
-  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
-
-!--- SURFACE GEOPOTENTIAL
-  title='Z'
-  ALLOCATE(z(iml,jml))
-  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
-  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
-      lon_in,  lat_in,  lon_in2, lat_in2, z)
-
-!--- SURFACE PRESSURE
-  title='SP'
-  ALLOCATE(psol_dyn(iml,jml))
-  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
-  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
-      lon_in,  lat_in,  lon_in2, lat_in2, psol_dyn)
-
-  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
-
-!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
-  IF(.NOT.ALLOCATED(tsol)) THEN
-    CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-  ELSE
-    IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
-      WRITE(lunout,*)'start_init_dyn :'
-      WRITE(lunout,*)'The temperature field we have does not have the right size'
-      STOP
-    END IF
-  END IF
-
-  IF(.NOT.ALLOCATED(phis)) THEN
-    CALL start_init_orog(iml,jml,lon_in,lat_in)
-  ELSE
-    IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
-      WRITE(lunout,*)'start_init_dyn :'
-      WRITE(lunout,*)'The orography field we have does not have the right size'
-      STOP
-    END IF
-  END IF
-
-!--- 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))
-    END DO
-    psol_dyn(iml,j) = psol_dyn(1,j)
-  END DO
-  DEALLOCATE(z)
-
-  ALLOCATE(xppn(iml-1),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)
-  END DO
-  DO i = 1, iml
-    psol_dyn(i,1  ) = SUM(xppn)/apoln
-    psol_dyn(i,jml) = SUM(xpps)/apols
-  END DO
-  DEALLOCATE(xppn,xpps) 
-
-  RETURN
-
-END SUBROUTINE start_init_dyn
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
-     lon_in2, lat_in2, pls_in, var3d, ibar)
-
-  use pchsp_95_m, only: pchsp_95
-  use pchfe_95_m, only: pchfe_95
-
-! Arguments:
-  CHARACTER(LEN=*),             INTENT(IN)    :: varname
-  INTEGER,                      INTENT(IN)    :: iml, jml, lml
-  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
-  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
-  INTEGER,                      INTENT(IN)    :: jml2
-  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
-  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
-  LOGICAL,                      INTENT(IN)    :: ibar
-!----------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  LOGICAL:: check=.TRUE., skip
-  REAL                  chmin, chmax
-  INTEGER ii, ij, il, ierr
-  integer n_extrap ! number of extrapolated points
-  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
-
-!---------------------------------------------------------------------------
-  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
-       ttm_dyn
-  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
-       llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
-  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
-       var_ana3d)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
-  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
-                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini, lat_ini)
-
-!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  DO il=1, llm_dyn
-    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
-         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
-         lon_in2, lat_in2, var_tmp3d(:, :, il))
-  END DO
-  DEALLOCATE(lon_rad, lat_rad)
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ax = lev_dyn(llm_dyn:1:-1) 
-  skip = .false.
-  n_extrap = 0
-  DO ij=1, jml
-    DO ii=1, iml-1
-      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
-      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
-      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
-           var3d(ii, ij, lml:1:-1), ierr)
-      if (ierr < 0) stop 1
-      n_extrap = n_extrap + ierr
-    END DO
-  END DO
-  if (n_extrap /= 0) then
-     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
-  end if
-  var3d(iml, :, :) = var3d(1, :, :) 
-
-  DO il=1, lml
-    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
-    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
-  END DO
-
-END SUBROUTINE start_inter_3d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj,    lon,  lat,  vari,     &
-                                 i1, j1, j2, lon1, lat1, lon2, lat2, varo)
-!
-!-------------------------------------------------------------------------------
-
-  USE inter_barxy_m, only: inter_barxy
-
-! Arguments:
-  CHARACTER(LEN=*),       INTENT(IN)  :: vname
-  LOGICAL,                INTENT(IN)  :: ibar, ibeg
-  INTEGER,                INTENT(IN)  :: ii, jj
-  REAL, DIMENSION(ii),    INTENT(IN)  :: lon
-  REAL, DIMENSION(jj),    INTENT(IN)  :: lat
-  REAL, DIMENSION(ii,jj), INTENT(IN)  :: vari
-  INTEGER,                INTENT(IN)  :: i1, j1, j2
-  REAL, DIMENSION(i1),    INTENT(IN)  :: lon1
-  REAL, DIMENSION(j1),    INTENT(IN)  :: lat1
-  REAL, DIMENSION(i1),    INTENT(IN)  :: lon2
-  REAL, DIMENSION(j2),    INTENT(IN)  :: lat2
-  REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(i1-1,j1) :: vtmp
-!-------------------------------------------------------------------------------
-  IF(ibar) THEN
-    IF(ibeg) THEN
-      WRITE(lunout,*)                                                          &
-               '---------------------------------------------------------------'
-      WRITE(lunout,*)                                                          &
- '$$$ Utilisation de l interpolation barycentrique  pour  '//TRIM(vname)//' $$$'
-      WRITE(lunout,*)                                                          &
-               '---------------------------------------------------------------'
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
-  ELSE
-    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
-  END IF
-  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
-
-END SUBROUTINE interp_startvar
-!
-!-------------------------------------------------------------------------------
-
-#endif
-! of #ifdef CPP_EARTH
-
-END MODULE startvar
-!
-!*******************************************************************************
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/grid_noro.F
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/grid_noro.F	(revision 780)
+++ 	(revision )
@@ -1,521 +1,0 @@
-!
-! $Id: grid_noro.F 1442 2010-10-18 08:31:31Z jghattas $
-!
-c
-c
-      SUBROUTINE grid_noro(imdep, jmdep, xdata, ydata, zdata,
-     .             imar, jmar, x, y,
-     .             zphi,zmea,zstd,zsig,zgam,zthe,
-     .             zpic,zval,mask)
-c=======================================================================
-c (F. Lott) (voir aussi z.x. Li, A. Harzallah et L. Fairhead)
-c
-c      Compute the Parameters of the SSO scheme as described in
-c      LOTT & MILLER (1997) and LOTT(1999).
-c      Target points are on a rectangular grid:
-c      iim+1 latitudes including North and South Poles;
-c      jjm+1 longitudes, with periodicity jjm+1=1.
-c      aux poles.  At the poles the fields value is repeated
-c      jjm+1 time.
-c      The parameters a,b,c,d represent the limite of the target
-c      gridpoint region. The means over this region are calculated
-c      from USN data, ponderated by a weight proportional to the 
-c      surface occupated by the data inside the model gridpoint area.
-c      In most circumstances, this weight is the ratio between the
-c      surface of the USN gridpoint area and the surface of the
-c      model gridpoint area. 
-c
-c           (c)
-c        ----d-----
-c        | . . . .|
-c        |        |
-c     (b)a . * . .b(a)
-c        |        |
-c        | . . . .|
-c        ----c-----
-c           (d)
-C=======================================================================
-c INPUT:
-c        imdep, jmdep: dimensions X and Y input field
-c        xdata, ydata: coordinates X and Y input field
-c        zdata: Input field
-c        In this version it is assumed that the entry data come from
-c        the USNavy dataset: imdep=iusn=2160, jmdep=jusn=1080.
-c OUTPUT:
-c        imar, jmar: dimensions X and Y Output field
-c        x, y: ccordinates  X and Y Output field.
-c             zmea:  Mean orographie   
-c             zstd:  Standard deviation
-c             zsig:  Slope
-c             zgam:  Anisotropy
-c             zthe:  Orientation of the small axis
-c             zpic:  Maximum altitude
-c             zval:  Minimum altitude
-C=======================================================================
-
-      IMPLICIT INTEGER (I,J)
-      IMPLICIT REAL(X,Z) 
-      
-	  parameter(iusn=2160,jusn=1080,iext=216, epsfra = 1.e-5)
-#include "dimensions.h"
-	  REAL xusn(iusn+2*iext),yusn(jusn+2)	
-      REAL zusn(iusn+2*iext,jusn+2)
-
-      INTEGER imdep, jmdep
-      REAL xdata(imdep),ydata(jmdep) 
-      REAL zdata(imdep,jmdep)
-c
-      INTEGER imar, jmar
-  
-C INTERMEDIATE FIELDS  (CORRELATIONS OF OROGRAPHY GRADIENT)
-
-      REAL ztz(iim+1,jjm+1),zxtzx(iim+1,jjm+1)
-      REAL zytzy(iim+1,jjm+1),zxtzy(iim+1,jjm+1)
-      REAL weight(iim+1,jjm+1)
-
-C CORRELATIONS OF USN OROGRAPHY GRADIENTS
-
-      REAL zxtzxusn(iusn+2*iext,jusn+2),zytzyusn(iusn+2*iext,jusn+2)
-      REAL zxtzyusn(iusn+2*iext,jusn+2)
-      REAL x(imar+1),y(jmar),zphi(imar+1,jmar)
-      REAL zmea(imar+1,jmar),zstd(imar+1,jmar)
-      REAL zmea0(imar+1,jmar) ! GK211005 (CG)
-      REAL zsig(imar+1,jmar),zgam(imar+1,jmar),zthe(imar+1,jmar)
-      REAL zpic(imar+1,jmar),zval(imar+1,jmar)
-cxxx PB     integer mask(imar+1,jmar)
-      real mask(imar+1,jmar), mask_tmp(imar+1,jmar)
-      real num_tot(2200,1100),num_lan(2200,1100)
-c
-      REAL a(2200),b(2200),c(1100),d(1100)
-      logical masque_lu
-c
-      print *,' parametres de l orographie a l echelle sous maille' 
-      xpi=acos(-1.)
-      rad    = 6 371 229.
-      zdeltay=2.*xpi/REAL(jusn)*rad
-c
-c utilise-t'on un masque lu?
-c
-      masque_lu = .true.
-      if (maxval(mask) == -99999 .and. minval(mask) == -99999) then
-        masque_lu= .false.
-        masque = 0.0
-      endif
-      write(*,*)'Masque lu', masque_lu
-c
-c  quelques tests de dimensions:
-c    
-c
-      if(iim.ne.imar) STOP 'Problem dim. x'
-      if(jjm.ne.jmar-1) STOP 'Problem dim. y'
-      IF (imar.GT.2200 .OR. jmar.GT.1100) THEN
-         PRINT*, 'imar or jmar too big', imar, jmar
-         CALL ABORT
-      ENDIF
-
-      IF(imdep.ne.iusn.or.jmdep.ne.jusn)then
-         print *,' imdep or jmdep bad dimensions:',imdep,jmdep
-         call abort
-      ENDIF
-
-      IF(imar+1.ne.iim+1.or.jmar.ne.jjm+1)THEN
-        print *,' imar or jmar bad dimensions:',imar,jmar
-        call abort
-      ENDIF
-
-
-c      print *,'xdata:',xdata
-c      print *,'ydata:',ydata
-c      print *,'x:',x
-c      print *,'y:',y
-c
-C  EXTENSION OF THE USN DATABASE TO POCEED COMPUTATIONS AT
-C  BOUNDARIES:
-c
-      DO j=1,jusn
-        yusn(j+1)=ydata(j)
-      DO i=1,iusn
-        zusn(i+iext,j+1)=zdata(i,j)
-        xusn(i+iext)=xdata(i)
-      ENDDO
-      DO i=1,iext
-        zusn(i,j+1)=zdata(iusn-iext+i,j)
-        xusn(i)=xdata(iusn-iext+i)-2.*xpi
-        zusn(iusn+iext+i,j+1)=zdata(i,j)
-        xusn(iusn+iext+i)=xdata(i)+2.*xpi
-      ENDDO
-      ENDDO
-
-        yusn(1)=ydata(1)+(ydata(1)-ydata(2))
-        yusn(jusn+2)=ydata(jusn)+(ydata(jusn)-ydata(jusn-1))
-       DO i=1,iusn/2+iext
-        zusn(i,1)=zusn(i+iusn/2,2)
-        zusn(i+iusn/2+iext,1)=zusn(i,2)
-        zusn(i,jusn+2)=zusn(i+iusn/2,jusn+1)
-        zusn(i+iusn/2+iext,jusn+2)=zusn(i,jusn+1)
-       ENDDO
-c  
-c COMPUTE LIMITS OF MODEL GRIDPOINT AREA
-C     ( REGULAR GRID)
-c
-      a(1) = x(1) - (x(2)-x(1))/2.0
-      b(1) = (x(1)+x(2))/2.0
-      DO i = 2, imar
-         a(i) = b(i-1)
-         b(i) = (x(i)+x(i+1))/2.0
-      ENDDO
-      a(imar+1) = b(imar)
-      b(imar+1) = x(imar+1) + (x(imar+1)-x(imar))/2.0
-
-      c(1) = y(1) - (y(2)-y(1))/2.0
-      d(1) = (y(1)+y(2))/2.0
-      DO j = 2, jmar-1
-         c(j) = d(j-1)
-         d(j) = (y(j)+y(j+1))/2.0
-      ENDDO
-      c(jmar) = d(jmar-1)
-      d(jmar) = y(jmar) + (y(jmar)-y(jmar-1))/2.0
-c
-c  initialisations:
-c
-      DO i = 1, imar+1
-      DO j = 1, jmar
-         weight(i,j) = 0.0
-         zxtzx(i,j)  = 0.0
-         zytzy(i,j)  = 0.0
-         zxtzy(i,j)  = 0.0
-         ztz(i,j)    = 0.0
-         zmea(i,j)   = 0.0
-         zpic(i,j)  =-1.E+10
-         zval(i,j)  = 1.E+10
-      ENDDO
-      ENDDO
-c
-c  COMPUTE SLOPES CORRELATIONS ON USN GRID
-c
-         DO j = 1,jusn+2 
-         DO i = 1, iusn+2*iext
-            zytzyusn(i,j)=0.0
-            zxtzxusn(i,j)=0.0
-            zxtzyusn(i,j)=0.0
-         ENDDO
-         ENDDO
-
-
-         DO j = 2,jusn+1 
-            zdeltax=zdeltay*cos(yusn(j))
-         DO i = 2, iusn+2*iext-1
-            zytzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))**2/zdeltay**2
-            zxtzxusn(i,j)=(zusn(i+1,j)-zusn(i-1,j))**2/zdeltax**2
-            zxtzyusn(i,j)=(zusn(i,j+1)-zusn(i,j-1))/zdeltay
-     *                   *(zusn(i+1,j)-zusn(i-1,j))/zdeltax
-         ENDDO
-         ENDDO
-c
-c  SUMMATION OVER GRIDPOINT AREA
-c 
-      zleny=xpi/REAL(jusn)*rad
-      xincr=xpi/2./REAL(jusn)
-       DO ii = 1, imar+1
-       DO jj = 1, jmar
-       num_tot(ii,jj)=0.
-       num_lan(ii,jj)=0.
-c        PRINT *,' iteration ii jj:',ii,jj
-         DO j = 2,jusn+1 
-c         DO j = 3,jusn 
-            zlenx=zleny*cos(yusn(j))
-            zdeltax=zdeltay*cos(yusn(j))
-            zbordnor=(c(jj)-yusn(j)+xincr)*rad
-            zbordsud=(yusn(j)-d(jj)+xincr)*rad
-            weighy=AMAX1(0.,
-     *             amin1(zbordnor,zbordsud,zleny))
-         IF(weighy.ne.0)THEN
-         DO i = 2, iusn+2*iext-1
-            zbordest=(xusn(i)-a(ii)+xincr)*rad*cos(yusn(j))
-            zbordoue=(b(ii)+xincr-xusn(i))*rad*cos(yusn(j))
-            weighx=AMAX1(0.,
-     *             amin1(zbordest,zbordoue,zlenx))
-            IF(weighx.ne.0)THEN
-            num_tot(ii,jj)=num_tot(ii,jj)+1.0
-            if(zusn(i,j).ge.1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
-            weight(ii,jj)=weight(ii,jj)+weighx*weighy
-            zxtzx(ii,jj)=zxtzx(ii,jj)+zxtzxusn(i,j)*weighx*weighy
-            zytzy(ii,jj)=zytzy(ii,jj)+zytzyusn(i,j)*weighx*weighy
-            zxtzy(ii,jj)=zxtzy(ii,jj)+zxtzyusn(i,j)*weighx*weighy
-            ztz(ii,jj)  =ztz(ii,jj)  +zusn(i,j)*zusn(i,j)*weighx*weighy
-c mean
-            zmea(ii,jj) =zmea(ii,jj)+zusn(i,j)*weighx*weighy
-c peacks
-            zpic(ii,jj)=amax1(zpic(ii,jj),zusn(i,j))
-c valleys
-            zval(ii,jj)=amin1(zval(ii,jj),zusn(i,j))
-            ENDIF
-         ENDDO
-         ENDIF
-         ENDDO
-       ENDDO
-       ENDDO
-c
-c  COMPUTE PARAMETERS NEEDED BY THE LOTT & MILLER (1997) AND
-C  LOTT (1999) SSO SCHEME.
-c
-      zllmmea=0.
-      zllmstd=0.
-      zllmsig=0.
-      zllmgam=0.
-      zllmpic=0.
-      zllmval=0.
-      zllmthe=0.
-      zminthe=0.
-c     print 100,' '
-c100  format(1X,A1,'II JJ',4X,'H',8X,'SD',8X,'SI',3X,'GA',3X,'TH') 
-       DO ii = 1, imar+1
-       DO jj = 1, jmar
-         IF (weight(ii,jj) .NE. 0.0) THEN
-c  Mask
-cXXX           if(num_lan(ii,jj)/num_tot(ii,jj).ge.0.5)then
-cXXX             mask(ii,jj)=1
-cXXX           else
-cXXX             mask(ii,jj)=0
-cXXX           ENDIF
-             if (.not. masque_lu) then
-               mask(ii,jj) = num_lan(ii,jj)/num_tot(ii,jj)
-             endif
-c  Mean Orography:
-           zmea (ii,jj)=zmea (ii,jj)/weight(ii,jj)
-           zxtzx(ii,jj)=zxtzx(ii,jj)/weight(ii,jj)
-           zytzy(ii,jj)=zytzy(ii,jj)/weight(ii,jj)
-           zxtzy(ii,jj)=zxtzy(ii,jj)/weight(ii,jj)
-           ztz(ii,jj)  =ztz(ii,jj)/weight(ii,jj)
-c  Standard deviation:
-           zstd(ii,jj)=sqrt(AMAX1(0.,ztz(ii,jj)-zmea(ii,jj)**2))
-         ELSE
-            PRINT*, 'probleme,ii,jj=', ii,jj
-         ENDIF
-       ENDDO
-       ENDDO
-
-C CORRECT VALUES OF HORIZONTAL SLOPE NEAR THE POLES:
-
-       DO ii = 1, imar+1
-         zxtzx(ii,1)=zxtzx(ii,2)
-         zxtzx(ii,jmar)=zxtzx(ii,jmar-1)
-         zxtzy(ii,1)=zxtzy(ii,2)
-         zxtzy(ii,jmar)=zxtzy(ii,jmar-1)
-         zytzy(ii,1)=zytzy(ii,2)
-         zytzy(ii,jmar)=zytzy(ii,jmar-1)
-       ENDDO
-
-C  FILTERS TO SMOOTH OUT FIELDS FOR INPUT INTO SSO SCHEME.
-
-C  FIRST FILTER, MOVING AVERAGE OVER 9 POINTS.
-
-       zmea0(:,:) = zmea(:,:) ! GK211005 (CG) on sauvegarde la topo non lissee
-       CALL MVA9(zmea,iim+1,jjm+1)
-       CALL MVA9(zstd,iim+1,jjm+1)
-       CALL MVA9(zpic,iim+1,jjm+1)
-       CALL MVA9(zval,iim+1,jjm+1)
-       CALL MVA9(zxtzx,iim+1,jjm+1)
-       CALL MVA9(zxtzy,iim+1,jjm+1) 
-       CALL MVA9(zytzy,iim+1,jjm+1)
-CXXX   Masque prenant en compte maximum de terre
-CXXX  On seuil a 10% de terre de terre car en dessous les parametres de surface n'on
-CXXX pas de sens (PB)
-       mask_tmp= 0.0
-       WHERE(mask .GE. 0.1) mask_tmp = 1.
-
-       DO ii = 1, imar
-       DO jj = 1, jmar
-         IF (weight(ii,jj) .NE. 0.0) THEN
-c  Coefficients K, L et M:
-           xk=(zxtzx(ii,jj)+zytzy(ii,jj))/2.
-           xl=(zxtzx(ii,jj)-zytzy(ii,jj))/2.
-           xm=zxtzy(ii,jj)
-           xp=xk-sqrt(xl**2+xm**2)
-           xq=xk+sqrt(xl**2+xm**2)
-           xw=1.e-8
-           if(xp.le.xw) xp=0.
-           if(xq.le.xw) xq=xw
-           if(abs(xm).le.xw) xm=xw*sign(1.,xm)
-c slope: 
-cXXX           zsig(ii,jj)=sqrt(xq)*mask(ii,jj)
-cXXXc isotropy:
-cXXX           zgam(ii,jj)=xp/xq*mask(ii,jj)
-cXXXc angle theta:
-cXXX           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask(ii,jj)
-cXXX           zphi(ii,jj)=zmea(ii,jj)*mask(ii,jj)
-cXXX           zmea(ii,jj)=zmea(ii,jj)*mask(ii,jj)
-cXXX           zpic(ii,jj)=zpic(ii,jj)*mask(ii,jj)
-cXXX           zval(ii,jj)=zval(ii,jj)*mask(ii,jj)
-cXXX           zstd(ii,jj)=zstd(ii,jj)*mask(ii,jj)
-CXX* PB modif pour maque de terre fractionnaire
-c slope: 
-           zsig(ii,jj)=sqrt(xq)*mask_tmp(ii,jj)
-c isotropy:
-           zgam(ii,jj)=xp/xq*mask_tmp(ii,jj)
-c angle theta:
-           zthe(ii,jj)=57.29577951*atan2(xm,xl)/2.*mask_tmp(ii,jj)
-           ! GK211005 (CG) ne pas forcement lisser la topo
-           ! zphi(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
-           zphi(ii,jj)=zmea0(ii,jj)*mask_tmp(ii,jj)
-           !
-           zmea(ii,jj)=zmea(ii,jj)*mask_tmp(ii,jj)
-           zpic(ii,jj)=zpic(ii,jj)*mask_tmp(ii,jj)
-           zval(ii,jj)=zval(ii,jj)*mask_tmp(ii,jj)
-           zstd(ii,jj)=zstd(ii,jj)*mask_tmp(ii,jj)
-c          print 101,ii,jj,
-c    *           zmea(ii,jj),zstd(ii,jj),zsig(ii,jj),zgam(ii,jj),
-c    *           zthe(ii,jj)
-c101  format(1x,2(1x,i2),2(1x,f7.1),1x,f7.4,2x,f4.2,1x,f5.1)     
-         ELSE
-c           PRINT*, 'probleme,ii,jj=', ii,jj
-         ENDIF
-      zllmmea=AMAX1(zmea(ii,jj),zllmmea)
-      zllmstd=AMAX1(zstd(ii,jj),zllmstd)
-      zllmsig=AMAX1(zsig(ii,jj),zllmsig)
-      zllmgam=AMAX1(zgam(ii,jj),zllmgam)
-      zllmthe=AMAX1(zthe(ii,jj),zllmthe)
-      zminthe=amin1(zthe(ii,jj),zminthe)
-      zllmpic=AMAX1(zpic(ii,jj),zllmpic)
-      zllmval=AMAX1(zval(ii,jj),zllmval)
-       ENDDO
-       ENDDO
-      print *,'  MEAN ORO:',zllmmea
-      print *,'  ST. DEV.:',zllmstd
-      print *,'  PENTE:',zllmsig
-      print *,' ANISOTROP:',zllmgam
-      print *,'  ANGLE:',zminthe,zllmthe	
-      print *,'  pic:',zllmpic
-      print *,'  val:',zllmval
-      
-C
-c gamma and theta a 1. and 0. at poles
-c
-      DO jj=1,jmar
-      zmea(imar+1,jj)=zmea(1,jj)
-      zphi(imar+1,jj)=zphi(1,jj)
-      zpic(imar+1,jj)=zpic(1,jj)
-      zval(imar+1,jj)=zval(1,jj)
-      zstd(imar+1,jj)=zstd(1,jj)
-      zsig(imar+1,jj)=zsig(1,jj)
-      zgam(imar+1,jj)=zgam(1,jj)
-      zthe(imar+1,jj)=zthe(1,jj)
-      ENDDO
-
-
-      zmeanor=0.0
-      zmeasud=0.0
-      zstdnor=0.0
-      zstdsud=0.0
-      zsignor=0.0
-      zsigsud=0.0
-      zweinor=0.0
-      zweisud=0.0
-      zpicnor=0.0
-      zpicsud=0.0                                   
-      zvalnor=0.0
-      zvalsud=0.0 
-
-      DO ii=1,imar
-      zweinor=zweinor+              weight(ii,   1)
-      zweisud=zweisud+              weight(ii,jmar)
-      zmeanor=zmeanor+zmea(ii,   1)*weight(ii,   1)
-      zmeasud=zmeasud+zmea(ii,jmar)*weight(ii,jmar)
-      zstdnor=zstdnor+zstd(ii,   1)*weight(ii,   1)
-      zstdsud=zstdsud+zstd(ii,jmar)*weight(ii,jmar)
-      zsignor=zsignor+zsig(ii,   1)*weight(ii,   1)
-      zsigsud=zsigsud+zsig(ii,jmar)*weight(ii,jmar)
-      zpicnor=zpicnor+zpic(ii,   1)*weight(ii,   1)
-      zpicsud=zpicsud+zpic(ii,jmar)*weight(ii,jmar)
-      zvalnor=zvalnor+zval(ii,   1)*weight(ii,   1)
-      zvalsud=zvalsud+zval(ii,jmar)*weight(ii,jmar)
-      ENDDO
-
-      DO ii=1,imar+1
-      zmea(ii,   1)=zmeanor/zweinor
-      zmea(ii,jmar)=zmeasud/zweisud
-      zphi(ii,   1)=zmeanor/zweinor
-      zphi(ii,jmar)=zmeasud/zweisud
-      zpic(ii,   1)=zpicnor/zweinor
-      zpic(ii,jmar)=zpicsud/zweisud
-      zval(ii,   1)=zvalnor/zweinor
-      zval(ii,jmar)=zvalsud/zweisud
-      zstd(ii,   1)=zstdnor/zweinor
-      zstd(ii,jmar)=zstdsud/zweisud
-      zsig(ii,   1)=zsignor/zweinor
-      zsig(ii,jmar)=zsigsud/zweisud
-      zgam(ii,   1)=1.
-      zgam(ii,jmar)=1.
-      zthe(ii,   1)=0.
-      zthe(ii,jmar)=0.
-      ENDDO
-
-      RETURN
-      END
-
-      SUBROUTINE MVA9(X,IMAR,JMAR)
-
-C MAKE A MOVING AVERAGE OVER 9 GRIDPOINTS OF THE X FIELDS
-
-      REAL X(IMAR,JMAR),XF(IMAR,JMAR)
-      real WEIGHTpb(-1:1,-1:1)
-
-      
-      SUM=0.
-      DO IS=-1,1
-        DO JS=-1,1
-          WEIGHTpb(IS,JS)=1./REAL((1+IS**2)*(1+JS**2))
-          SUM=SUM+WEIGHTpb(IS,JS)
-        ENDDO
-      ENDDO
-      
-c     WRITE(*,*) 'MVA9 ', IMAR, JMAR
-c     WRITE(*,*) 'MVA9 ', WEIGHTpb
-c     WRITE(*,*) 'MVA9 SUM ', SUM
-      DO IS=-1,1
-        DO JS=-1,1
-          WEIGHTpb(IS,JS)=WEIGHTpb(IS,JS)/SUM
-        ENDDO
-      ENDDO
-
-      DO J=2,JMAR-1
-        DO I=2,IMAR-1
-          XF(I,J)=0.
-          DO IS=-1,1
-            DO JS=-1,1
-              XF(I,J)=XF(I,J)+X(I+IS,J+JS)*WEIGHTpb(IS,JS)
-            ENDDO
-          ENDDO
-        ENDDO
-      ENDDO
-
-      DO J=2,JMAR-1
-        XF(1,J)=0.
-        IS=IMAR-1
-        DO JS=-1,1 
-          XF(1,J)=XF(1,J)+X(IS,J+JS)*WEIGHTpb(-1,JS)
-        ENDDO
-        DO IS=0,1 
-          DO JS=-1,1 
-            XF(1,J)=XF(1,J)+X(1+IS,J+JS)*WEIGHTpb(IS,JS)
-          ENDDO
-        ENDDO
-        XF(IMAR,J)=XF(1,J)
-      ENDDO
-
-      DO I=1,IMAR
-        XF(I,1)=XF(I,2)
-        XF(I,JMAR)=XF(I,JMAR-1)
-      ENDDO
-
-      DO I=1,IMAR
-        DO J=1,JMAR
-          X(I,J)=XF(I,J)
-        ENDDO
-      ENDDO
-
-      RETURN
-      END
-
-
-
Index: trunk/LMDZ.COMMON/libf/dyn3dpar/startvar.F90
===================================================================
--- trunk/LMDZ.COMMON/libf/dyn3dpar/startvar.F90	(revision 780)
+++ 	(revision )
@@ -1,783 +1,0 @@
-!
-! $Id: startvar.F90 1425 2010-09-02 13:45:23Z lguez $
-!
-!*******************************************************************************
-!
-MODULE startvar
-!
-!*******************************************************************************
-!
-!-------------------------------------------------------------------------------
-! Purpose: Access data from the database of atmospheric to initialize the model.
-!-------------------------------------------------------------------------------
-! Comments:
-!
-!    *  This module is designed to work for Earth (and with ioipsl)
-!
-!    *  There are three ways to acces data, depending 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(...)
-!
-!  - A 1D variable on the physical grid :
-!    CALL startget_phys1d((varname, iml, jml,  lon_in,  lat_in,  nbindex,              &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
-!
-!  - A 2D variable on the dynamical grid :
-!    CALL startget_phys2d(varname, iml, jml,  lon_in,  lat_in,                        &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )             
-!
-!  - A 3D variable on the dynamical grid :
-!    CALL startget_dyn((varname, iml, jml,  lon_in,  lat_in,  lml, pls, workvar,    &
-!           champ, val_exp,      jml2, lon_in2, lat_in2, ibar )
-!
-!    *  Data needs to be in NetCDF format
-!
-!    *  Variables should have the following names in the files:
-!            '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
-!
-!   *   There is a big mess with the longitude size. 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 smaller fields when needed. This needs to be cleared once and for
-!  all in LMDZ. A convention is required.
-!-------------------------------------------------------------------------------
-#ifdef CPP_EARTH
-  USE ioipsl
-  IMPLICIT NONE
-
-  PRIVATE
-  PUBLIC startget_phys2d, startget_phys1d, startget_dyn
-!  INTERFACE startget
-!    MODULE PROCEDURE startget_phys1d, startget_phys2d, startget_dyn
-!  END INTERFACE
-
-  REAL,    SAVE :: deg2rad,  pi
-  INTEGER, SAVE ::           iml_rel,  jml_rel
-  INTEGER, SAVE :: fid_phys, iml_phys, jml_phys
-  INTEGER, SAVE :: fid_dyn,  iml_dyn,  jml_dyn,  llm_dyn,  ttm_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_phys, lon_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_phys, lat_dyn
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lon_rug, lon_alb, lon_rel
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: lat_rug, lat_alb, lat_rel
-  REAL, DIMENSION(:),     ALLOCATABLE, TARGET, SAVE :: levdyn_ini
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: relief, zstd, zsig, zgam
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: masque, zthe, zpic, zval
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: rugo, phis, tsol, qsol
-  REAL, DIMENSION(:,:),   ALLOCATABLE, TARGET, SAVE :: psol_dyn
-  REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET, SAVE :: var_ana3d
-
-   CONTAINS
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_phys1d(varname, iml, jml, lon_in, lat_in, nbindex, champ,  &
-                           val_exp ,jml2, lon_in2, lat_in2, ibar)
-!
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),         INTENT(IN)    :: varname
-  INTEGER,                  INTENT(IN)    :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)    :: lat_in
-  INTEGER,                  INTENT(IN)    :: nbindex
-  REAL, DIMENSION(nbindex), INTENT(INOUT) :: champ
-  REAL,                     INTENT(IN)    :: val_exp
-  INTEGER,                  INTENT(IN)    :: jml2
-  REAL, DIMENSION(iml),     INTENT(IN)    :: lon_in2
-  REAL, DIMENSION(jml2),    INTENT(IN)    :: lat_in2
-  LOGICAL,                  INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(:,:), POINTER :: v2d
-!-------------------------------------------------------------------------------
-  v2d=>NULL()
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
-    SELECT CASE(varname)
-      CASE('tsol')
-        IF(.NOT.ALLOCATED(tsol))                                               &
-         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('qsol')
-        IF(.NOT.ALLOCATED(qsol))                                               &
-         CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('psol')
-        IF(.NOT.ALLOCATED(psol_dyn))                                           &
-         CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('zmea','zstd','zsig','zgam','zthe','zpic','zval')
-        IF(.NOT.ALLOCATED(relief))                                             &
-         CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('rads','snow','tslab','seaice','rugmer','agsno')
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_phys1d'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                     //' from any data set'; STOP
-    END SELECT
-
-!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
-    SELECT CASE(varname)
-      CASE('rads','snow','tslab','seaice');  champ=0.0
-      CASE('rugmer');                        champ(:)=0.001
-      CASE('agsno');                         champ(:)=50.0
-      CASE DEFAULT
-        SELECT CASE(varname)
-          CASE('tsol'); v2d=>tsol
-          CASE('qsol'); v2d=>qsol
-          CASE('psol'); v2d=>psol_dyn
-          CASE('zmea'); v2d=>relief
-          CASE('zstd'); v2d=>zstd
-          CASE('zsig'); v2d=>zsig
-          CASE('zgam'); v2d=>zgam
-          CASE('zthe'); v2d=>zthe
-          CASE('zpic'); v2d=>zpic
-          CASE('zval'); v2d=>zval
-        END SELECT
-        IF(SIZE(v2d)/=SIZE(lon_in)*SIZE(lat_in)) THEN
-         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
-         STOP
-        END IF
-        CALL gr_dyn_fi(1,iml,jml,nbindex,v2d,champ)
-    END SELECT
-
-  ELSE
-
-!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
-    SELECT CASE(varname)
-      CASE('tsol')
-        IF(.NOT.ALLOCATED(tsol)) ALLOCATE(tsol(iml,jml))
-        CALL gr_fi_dyn(1,iml,jml,nbindex,champ,tsol)
-    END SELECT
-
-  END IF
-
-END SUBROUTINE  startget_phys1d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_phys2d(varname, iml, jml, lon_in, lat_in, champ, val_exp,  &
-                           jml2, lon_in2, lat_in2 , ibar, msk)
-!
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),         INTENT(IN)           :: varname
-  INTEGER,                  INTENT(IN)           :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
-  REAL, DIMENSION(iml,jml), INTENT(INOUT)        :: champ
-  REAL,                     INTENT(IN)           :: val_exp
-  INTEGER,                  INTENT(IN)           :: jml2
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in2
-  REAL, DIMENSION(jml2),    INTENT(IN)           :: lat_in2
-  LOGICAL,                  INTENT(IN)           :: ibar
-  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: msk
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(:,:), POINTER :: v2d=>NULL()
-  LOGICAL                       :: lrelief1, lrelief2
-!-------------------------------------------------------------------------------
-  v2d=>NULL()
-  lrelief1=(.NOT.ALLOCATED(relief).AND.     PRESENT(msk))
-  lrelief2=(.NOT.ALLOCATED(relief).AND..NOT.PRESENT(msk))
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-!--- CHECKING IF THE FIELD IS KNOWN ; READING UNALLOCATED FILES
-    SELECT CASE(varname)
-      CASE('psol')
-        IF(.NOT.ALLOCATED(psol_dyn))                                           &
-          CALL start_init_dyn (iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-      CASE('relief')
-        IF(lrelief1)             CALL start_init_orog(iml,jml,lon_in,lat_in,msk)
-        IF(lrelief2)             CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('surfgeo')
-        IF(.NOT.ALLOCATED(phis)) CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE('rugosite','masque')
-        IF(.NOT.ALLOCATED(rugo)) CALL start_init_orog(iml,jml,lon_in,lat_in)
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_phys2d'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                     //' from any data set'; STOP
-    END SELECT
-
-!--- SELECTING v2d FOR WANTED VARIABLE AND CHEKING ITS SIZE
-    SELECT CASE(varname)
-      CASE('psol');     v2d=>psol_dyn
-      CASE('relief');   v2d=>relief
-      CASE('rugosite'); v2d=>rugo
-      CASE('masque');   v2d=>masque
-      CASE('surfgeo');  v2d=>phis
-    END SELECT
-    IF(SIZE(champ)/=SIZE(v2d)) THEN
-      WRITE(lunout,*) 'STARTVAR module has been initialized to the wrong size'
-      STOP
-    END IF
-
-    champ(:,:)=v2d(:,:)
-
-  ELSE
-
-!--- SOME FIELDS ARE CAUGHT: MAY BE NEEDED FOR A 3D INTEPROLATION
-    SELECT CASE(varname)
-      CASE ('surfgeo')
-        IF(.NOT.ALLOCATED(phis)) ALLOCATE(phis(iml,jml))
-        IF(SIZE(phis)/=SIZE(champ)) THEN
-         WRITE(lunout,*)'STARTVAR module has been initialized to the wrong size'
-         STOP
-        END IF
-        phis(:,:)=champ(:,:)
-    END SELECT
-
-  END IF
-
-END SUBROUTINE startget_phys2d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE startget_dyn(varname,  lon_in,  lat_in, pls,workvar,&
-                     champ, val_exp, lon_in2, lat_in2, ibar)
-
-      use assert_eq_m, only: assert_eq
-
-
-!-------------------------------------------------------------------------------
-! Comment:
-!   This routine only works if the variable does not exist or is constant.
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN)    :: varname
-  REAL, INTENT(IN)    :: lon_in(:) ! dim(iml)
-  REAL, INTENT(IN)    :: lat_in(:) ! dim(jml)
-  REAL, INTENT(IN)    :: pls(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(IN)    :: workvar(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(INOUT) :: champ(:, :, :) ! dim(iml, jml, lml)
-  REAL, INTENT(IN)    :: val_exp
-  REAL, INTENT(IN)    :: lon_in2(:) ! dim(iml)
-  REAL, INTENT(IN)    :: lat_in2(:) ! dim(jml2)
-  LOGICAL,                      INTENT(IN)    :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-#include "dimensions.h"
-#include "comconst.h"
-#include "paramet.h"
-#include "comgeom2.h"
-  INTEGER    :: iml, jml
-  INTEGER    :: lml
-  INTEGER    :: jml2
-  REAL, DIMENSION(:,:,:), POINTER :: v3d=>NULL()
-  CHARACTER(LEN=10) :: vname
-  INTEGER :: il
-  REAL    :: xppn, xpps
-!-------------------------------------------------------------------------------
-  NULLIFY(v3d)
-  IF(MINVAL(champ)==MAXVAL(champ).AND.MINVAL(champ)==val_exp) THEN
-
-      iml = assert_eq((/size(lon_in), size(pls, 1), size(workvar, 1), &
-     &     size(champ, 1), size(lon_in2)/), "startget_dyn iml")
-      jml = assert_eq(size(lat_in), size(pls, 2), size(workvar, 2),   &
-     &     size(champ, 2), "startget_dyn jml")
-      lml = assert_eq(size(pls, 3), size(workvar, 3), size(champ, 3), &
-     &     "startget_dyn lml")
-      jml2 = size(lat_in2)
-
-!--- READING UNALLOCATED FILES
-    IF(.NOT.ALLOCATED(psol_dyn))                                               &
-      CALL start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-
-!--- CHECKING IF THE FIELD IS KNOWN AND INTERPOLATING 3D FIELDS
-    SELECT CASE(varname)
-      CASE('u');        vname='U'
-      CASE('v');        vname='V'
-      CASE('t','tpot'); vname='TEMP'
-      CASE('q');        vname='R'
-      CASE DEFAULT
-        WRITE(lunout,*)'startget_dyn'
-        WRITE(lunout,*)'No rule is present to extract variable '//TRIM(varname)&
-                //' from any data set'; STOP
-    END SELECT
-    CALL start_inter_3d(TRIM(vname), iml, jml, lml, lon_in, lat_in, jml2,      &
-                        lon_in2, lat_in2,  pls, champ,ibar )
-
-!--- COMPUTING THE REQUIRED FILED
-    SELECT CASE(varname)
-      CASE('u')                                            !--- Eastward wind
-        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
-        champ(iml,:,:)=champ(1,:,:)
-
-      CASE('v')                                            !--- Northward wind
-        DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
-        champ(iml,:,:)=champ(1,:,:)
-
-      CASE('tpot')                                         !--- Temperature
-        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
-          champ=champ*cpp/workvar
-          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
-          END DO
-        ELSE
-          WRITE(lunout,*)'Could not compute potential temperature as the'
-          WRITE(lunout,*)'Exner function is missing or constant.'; STOP
-        END IF
-
-      CASE('q')                                            !--- Relat. humidity
-        IF(MINVAL(workvar)/=MAXVAL(workvar)) THEN
-          champ=0.01*champ*workvar
-          WHERE(champ<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
-          END DO
-        ELSE
-          WRITE(lunout,*)'Could not compute specific humidity as the'
-          WRITE(lunout,*)'saturated humidity is missing or constant.'; STOP
-        END IF
-
-    END SELECT
-
-  END IF
-
-END SUBROUTINE startget_dyn
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_orog(iml,jml,lon_in,lat_in,masque_lu)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,                  INTENT(IN)           :: iml, jml
-  REAL, DIMENSION(iml),     INTENT(IN)           :: lon_in
-  REAL, DIMENSION(jml),     INTENT(IN)           :: lat_in
-  REAL, DIMENSION(iml,jml), INTENT(IN), OPTIONAL :: masque_lu
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: orofname
-  LOGICAL               :: check=.TRUE.
-  REAL,    DIMENSION(1) :: lev
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: fid, llm_tmp, ttm_tmp
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: relief_hi, tmp_var
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-!-------------------------------------------------------------------------------
-  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-
-  orofname = 'Relief.nc'; title='RELIEF'
-  IF(check) WRITE(lunout,*)'Reading the high resolution orography'
-  CALL flininfo(orofname, iml_rel, jml_rel, llm_tmp, ttm_tmp, fid)
-
-  ALLOCATE(lat_rel(iml_rel,jml_rel),lon_rel(iml_rel,jml_rel))
-  CALL flinopen(orofname, .FALSE., iml_rel, jml_rel, llm_tmp, lon_rel, lat_rel,&
-                lev, ttm_tmp, itau, date, dt, fid)
-  ALLOCATE(relief_hi(iml_rel,jml_rel))
-  CALL flinget(fid, title, iml_rel, jml_rel, llm_tmp, ttm_tmp, 1, 1, relief_hi)
-  CALL flinclo(fid)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_rel),lat_ini(jml_rel))
-  lon_ini(:)=lon_rel(:,1); IF(MAXVAL(lon_rel)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_rel(1,:); IF(MAXVAL(lat_rel)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_rel),lat_rad(jml_rel))
-  CALL conf_dat2d(title, iml_rel, jml_rel, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  relief_hi, .FALSE.)
-  DEALLOCATE(lon_ini,lat_ini)
-
-!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  IF(check) WRITE(lunout,*)'Computes all parameters needed for gravity wave dra&
-     &g code'
-
-  ALLOCATE(phis(iml,jml))      ! Geopotentiel au sol
-  ALLOCATE(zstd(iml,jml))      ! Deviation standard de l'orographie sous-maille
-  ALLOCATE(zsig(iml,jml))      ! Pente de l'orographie sous-maille 
-  ALLOCATE(zgam(iml,jml))      ! Anisotropie de l'orographie sous maille
-  ALLOCATE(zthe(iml,jml))      ! Orientation axe +grande pente d'oro sous maille
-  ALLOCATE(zpic(iml,jml))      ! Hauteur pics de la SSO
-  ALLOCATE(zval(iml,jml))      ! Hauteur vallees de la SSO
-  ALLOCATE(relief(iml,jml))    ! Orographie moyenne
-  ALLOCATE(masque(iml,jml))    ! Masque terre ocean
-  masque = -99999.
-  IF(PRESENT(masque_lu)) masque=masque_lu
-
-  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, masque)
-  phis = phis * 9.81
-
-!--- SURFACE ROUGHNESS COMPUTATION (UNUSED FOR THE MOMENT !!! )
-  IF(check) WRITE(lunout,*)'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)
-  rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
-  DEALLOCATE(relief_hi,tmp_var,lon_rad,lat_rad)
-  RETURN
-
-END SUBROUTINE start_init_orog
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,               INTENT(IN) :: iml, jml
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
-  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
-  INTEGER,               INTENT(IN) :: jml2
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
-  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
-  LOGICAL,               INTENT(IN) :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: physfname
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: llm_tmp, ttm_tmp
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),   ALLOCATABLE :: levphys_ini
-!-------------------------------------------------------------------------------
-  physfname = 'ECPHY.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(check) WRITE(lunout,*)'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))
-  ALLOCATE(levphys_ini(llm_tmp))
-  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_tmp, lon_phys,     &
-                lat_phys, levphys_ini, ttm_tmp, itau, date, dt, fid_phys)
-  DEALLOCATE(levphys_ini)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_phys),lat_ini(jml_phys))
-  lon_ini(:)=lon_phys(:,1); IF(MAXVAL(lon_phys)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_phys(1,:); IF(MAXVAL(lat_phys)>pi) lat_ini=lat_ini*deg2rad
-
-  ALLOCATE(var_ana(iml_phys,jml_phys),lon_rad(iml_phys),lat_rad(jml_phys))
-
-!--- SURFACE TEMPERATURE
-  title='ST'
-  ALLOCATE(tsol(iml,jml))
-  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
-  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
-                  var_ana , ibar  )
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
-      lon_in,   lat_in,   lon_in2, lat_in2, tsol)
-
-!--- SOIL MOISTURE
-  title='CDSW'
-  ALLOCATE(qsol(iml,jml))
-  CALL flinget(fid_phys,title,iml_phys,jml_phys,llm_tmp,ttm_tmp,1,1,var_ana)
-  CALL conf_dat2d(title,iml_phys, jml_phys, lon_ini, lat_ini, lon_rad, lat_rad,&
-                  var_ana, ibar  )
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_phys, jml_phys, lon_rad, lat_rad, var_ana, iml, jml, jml-1,          &
-      lon_in,   lat_in,   lon_in2, lat_in2, qsol)
-
-  CALL flinclo(fid_phys)
-
-  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
-
-END SUBROUTINE start_init_phys
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_dyn(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-!
-!-------------------------------------------------------------------------------
-! Arguments:
-  INTEGER,               INTENT(IN) :: iml, jml
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in
-  REAL, DIMENSION(jml),  INTENT(IN) :: lat_in
-  INTEGER,               INTENT(IN) :: jml2
-  REAL, DIMENSION(iml),  INTENT(IN) :: lon_in2
-  REAL, DIMENSION(jml2), INTENT(IN) :: lat_in2
-  LOGICAL,               INTENT(IN) :: ibar
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-#include "dimensions.h"
-#include "paramet.h"
-#include "comgeom2.h"
-  CHARACTER(LEN=25)     :: title
-  CHARACTER(LEN=120)    :: physfname
-  LOGICAL               :: check=.TRUE.
-  REAL                  :: date, dt
-  INTEGER, DIMENSION(1) :: itau
-  INTEGER               :: i, j
-  REAL,    DIMENSION(:,:), ALLOCATABLE :: var_ana, z
-  REAL,    DIMENSION(:),   ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL,    DIMENSION(:),   ALLOCATABLE :: xppn, xpps
-!-------------------------------------------------------------------------------
-
-!--- KINETIC ENERGY
-  physfname = 'ECDYN.nc'; pi=2.0*ASIN(1.0); deg2rad=pi/180.0
-  IF(check) WRITE(lunout,*) 'Opening the surface analysis'
-  CALL flininfo(physfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
-  IF(check) WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
-
-  ALLOCATE(lat_dyn(iml_dyn,jml_dyn))
-  ALLOCATE(lon_dyn(iml_dyn,jml_dyn))
-  ALLOCATE(levdyn_ini(llm_dyn))
-  CALL flinopen(physfname, .FALSE., iml_dyn, jml_dyn, llm_dyn, lon_dyn,lat_dyn,&
-                levdyn_ini, ttm_dyn, itau, date, dt, fid_dyn)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn),lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:,1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1,:); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
-
-  ALLOCATE(var_ana(iml_dyn,jml_dyn),lon_rad(iml_dyn),lat_rad(jml_dyn))
-
-!--- SURFACE GEOPOTENTIAL
-  title='Z'
-  ALLOCATE(z(iml,jml))
-  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
-  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
-      lon_in,  lat_in,  lon_in2, lat_in2, z)
-
-!--- SURFACE PRESSURE
-  title='SP'
-  ALLOCATE(psol_dyn(iml,jml))
-  CALL flinget(fid_dyn, title, iml_dyn, jml_dyn, 0, ttm_dyn, 1, 1, var_ana)
-  CALL conf_dat2d(title, iml_dyn, jml_dyn, lon_ini, lat_ini, lon_rad, lat_rad, &
-                  var_ana, ibar)
-  CALL interp_startvar(title, ibar, .TRUE.,                                    &
-      iml_dyn, jml_dyn, lon_rad, lat_rad, var_ana, iml, jml, jml-1,            &
-      lon_in,  lat_in,  lon_in2, lat_in2, psol_dyn)
-
-  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
-
-!--- ALLOCATION OF VARIABLES CREATED IN OR COMING FROM RESTART FILE
-  IF(.NOT.ALLOCATED(tsol)) THEN
-    CALL start_init_phys(iml,jml,lon_in,lat_in,jml2,lon_in2,lat_in2,ibar)
-  ELSE
-    IF(SIZE(tsol)/=SIZE(psol_dyn)) THEN
-      WRITE(lunout,*)'start_init_dyn :'
-      WRITE(lunout,*)'The temperature field we have does not have the right size'
-      STOP
-    END IF
-  END IF
-
-  IF(.NOT.ALLOCATED(phis)) THEN
-    CALL start_init_orog(iml,jml,lon_in,lat_in)
-  ELSE
-    IF(SIZE(phis)/=SIZE(psol_dyn)) THEN
-      WRITE(lunout,*)'start_init_dyn :'
-      WRITE(lunout,*)'The orography field we have does not have the right size'
-      STOP
-    END IF
-  END IF
-
-!--- 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))
-    END DO
-    psol_dyn(iml,j) = psol_dyn(1,j)
-  END DO
-  DEALLOCATE(z)
-
-  ALLOCATE(xppn(iml-1),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)
-  END DO
-  DO i = 1, iml
-    psol_dyn(i,1  ) = SUM(xppn)/apoln
-    psol_dyn(i,jml) = SUM(xpps)/apols
-  END DO
-  DEALLOCATE(xppn,xpps) 
-
-  RETURN
-
-END SUBROUTINE start_init_dyn
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_inter_3d(varname, iml, jml, lml, lon_in, lat_in, jml2, &
-     lon_in2, lat_in2, pls_in, var3d, ibar)
-
-  use pchsp_95_m, only: pchsp_95
-  use pchfe_95_m, only: pchfe_95
-
-! Arguments:
-  CHARACTER(LEN=*),             INTENT(IN)    :: varname
-  INTEGER,                      INTENT(IN)    :: iml, jml, lml
-  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in
-  REAL, DIMENSION(jml),         INTENT(IN)    :: lat_in
-  INTEGER,                      INTENT(IN)    :: jml2
-  REAL, DIMENSION(iml),         INTENT(IN)    :: lon_in2
-  REAL, DIMENSION(jml2),        INTENT(IN)    :: lat_in2
-  REAL, DIMENSION(iml, jml, lml), INTENT(IN)    :: pls_in
-  REAL, DIMENSION(iml, jml, lml), INTENT(OUT)   :: var3d
-  LOGICAL,                      INTENT(IN)    :: ibar
-!----------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  LOGICAL:: check=.TRUE., skip
-  REAL                  chmin, chmax
-  INTEGER ii, ij, il, ierr
-  integer n_extrap ! number of extrapolated points
-  REAL, DIMENSION(iml, jml, llm_dyn):: var_tmp3d
-  REAL,    DIMENSION(:),     ALLOCATABLE :: lon_rad, lat_rad, lon_ini, lat_ini
-  REAL, DIMENSION(llm_dyn):: lev_dyn, ax, ay, yder
-
-!---------------------------------------------------------------------------
-  IF(check) WRITE(lunout, *)'Going into flinget to extract the 3D  field.'
-  IF(check) WRITE(lunout, *) fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, &
-       ttm_dyn
-  IF(check) WRITE(lunout, *) 'Allocating space for interpolation', iml, jml, &
-       llm_dyn
-
-  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
-  CALL flinget(fid_dyn, varname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, 1, 1, &
-       var_ana3d)
-
-!--- IF ANGLES ARE IN DEGREES, THEY ARE CONVERTED INTO RADIANS
-  ALLOCATE(lon_ini(iml_dyn), lat_ini(jml_dyn))
-  lon_ini(:)=lon_dyn(:, 1); IF(MAXVAL(lon_dyn)>pi) lon_ini=lon_ini*deg2rad
-  lat_ini(:)=lat_dyn(1, :); IF(MAXVAL(lat_dyn)>pi) lat_ini=lat_ini*deg2rad
-
-!--- FIELDS ARE PROCESSED TO BE ON STANDARD ANGULAR DOMAINS
-  ALLOCATE(lon_rad(iml_dyn), lat_rad(jml_dyn))
-  CALL conf_dat3d (varname, iml_dyn, jml_dyn, llm_dyn, lon_ini, lat_ini,      &
-                   levdyn_ini, lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
-  DEALLOCATE(lon_ini, lat_ini)
-
-!--- COMPUTING THE REQUIRED FIELDS USING ROUTINE grid_noro
-  DO il=1, llm_dyn
-    CALL interp_startvar(varname, ibar, il==1, iml_dyn, jml_dyn, lon_rad, &
-         lat_rad, var_ana3d(:, :, il), iml, jml, jml2, lon_in, lat_in, &
-         lon_in2, lat_in2, var_tmp3d(:, :, il))
-  END DO
-  DEALLOCATE(lon_rad, lat_rad)
-
-!--- VERTICAL INTERPOLATION IS PERFORMED FROM TOP OF ATMOSPHERE TO GROUND
-  ax = lev_dyn(llm_dyn:1:-1) 
-  skip = .false.
-  n_extrap = 0
-  DO ij=1, jml
-    DO ii=1, iml-1
-      ay = var_tmp3d(ii, ij, llm_dyn:1:-1)
-      yder = pchsp_95(ax, ay, ibeg=2, iend=2, vc_beg=0., vc_end=0.)
-      CALL pchfe_95(ax, ay, yder, skip, pls_in(ii, ij, lml:1:-1), &
-           var3d(ii, ij, lml:1:-1), ierr)
-      if (ierr < 0) stop 1
-      n_extrap = n_extrap + ierr
-    END DO
-  END DO
-  if (n_extrap /= 0) then
-     print *, "start_inter_3d pchfe_95: n_extrap = ", n_extrap
-  end if
-  var3d(iml, :, :) = var3d(1, :, :) 
-
-  DO il=1, lml
-    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
-    WRITE(lunout, *)' '//TRIM(varname)//'  min max l ', il, chmin, chmax
-  END DO
-
-END SUBROUTINE start_inter_3d
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE interp_startvar(vname, ibar, ibeg, ii, jj,    lon,  lat,  vari,     &
-                                 i1, j1, j2, lon1, lat1, lon2, lat2, varo)
-!
-!-------------------------------------------------------------------------------
-
-  USE inter_barxy_m, only: inter_barxy
-
-! Arguments:
-  CHARACTER(LEN=*),       INTENT(IN)  :: vname
-  LOGICAL,                INTENT(IN)  :: ibar, ibeg
-  INTEGER,                INTENT(IN)  :: ii, jj
-  REAL, DIMENSION(ii),    INTENT(IN)  :: lon
-  REAL, DIMENSION(jj),    INTENT(IN)  :: lat
-  REAL, DIMENSION(ii,jj), INTENT(IN)  :: vari
-  INTEGER,                INTENT(IN)  :: i1, j1, j2
-  REAL, DIMENSION(i1),    INTENT(IN)  :: lon1
-  REAL, DIMENSION(j1),    INTENT(IN)  :: lat1
-  REAL, DIMENSION(i1),    INTENT(IN)  :: lon2
-  REAL, DIMENSION(j2),    INTENT(IN)  :: lat2
-  REAL, DIMENSION(i1,j1), INTENT(OUT) :: varo
-!-------------------------------------------------------------------------------
-! Local variables:
-#include "iniprint.h"
-  REAL, DIMENSION(i1-1,j1) :: vtmp
-!-------------------------------------------------------------------------------
-  IF(ibar) THEN
-    IF(ibeg) THEN
-      WRITE(lunout,*)                                                          &
-               '---------------------------------------------------------------'
-      WRITE(lunout,*)                                                          &
- '$$$ Utilisation de l interpolation barycentrique  pour  '//TRIM(vname)//' $$$'
-      WRITE(lunout,*)                                                          &
-               '---------------------------------------------------------------'
-    END IF
-    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2(:j2), vtmp)
-  ELSE
-    CALL grille_m   (ii, jj,   lon, lat, vari, i1-1, j1, lon1, lat1,     vtmp)
-  END IF
-  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
-
-END SUBROUTINE interp_startvar
-!
-!-------------------------------------------------------------------------------
-
-#endif
-! of #ifdef CPP_EARTH
-
-END MODULE startvar
-!
-!*******************************************************************************
