Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2301)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0dyn_netcdf.F90	(revision 2301)
@@ -0,0 +1,770 @@
+MODULE etat0dyn
+!
+!*******************************************************************************
+! Purpose: Create dynamical initial state using atmospheric fields from a
+!          database of atmospheric to initialize the model.
+!-------------------------------------------------------------------------------
+! Comments:
+!
+!    *  This module is designed to work for Earth (and with ioipsl)
+!
+!    *  etat0dyn_netcdf routine can access to NetCDF data through the following
+!  routine (to be called after restget):
+!    CALL startget_dyn3d(varname, lon_in,  lat_in, pls, workvar,&
+!                          champ, lon_in2, lat_in2, ibar)
+!
+!    *  Variables should have the following names in the NetCDF files:
+!            'U'      : East ward wind              (in "ECDYN.nc")
+!            'V'      : Northward wind              (in "ECDYN.nc")
+!            'TEMP'   : Temperature                 (in "ECDYN.nc")
+!            'R'      : Relative humidity           (in "ECDYN.nc")
+!            'RELIEF' : High resolution orography   (in "Relief.nc") 
+!
+!    * The land mask and corresponding weights can be:
+!      1) already known (in particular if etat0dyn has been called before) ;
+!         in this case, ANY(masque(:,:)/=-99999.) = .TRUE.
+!      2) computed using the ocean mask from the ocean model (to ensure ocean
+!         fractions are the same for atmosphere and ocean) for coupled runs.
+!         File name: "o2a.nc"  ;  variable name: "OceMask"
+!      3) computed from topography file "Relief.nc" for forced runs.
+!
+!   *   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.
+!-------------------------------------------------------------------------------
+  USE ioipsl,         ONLY: flininfo, flinopen, flinget, flinclo, histclo
+  USE assert_eq_m,    ONLY: assert_eq
+#ifdef CPP_PHYS
+  USE indice_sol_mod, ONLY: epsfra
+#endif
+  IMPLICIT NONE
+
+  PRIVATE
+  PUBLIC :: etat0dyn_netcdf
+
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "temps.h"
+  include "comdissnew.h"
+  include "serre.h"
+  REAL, SAVE :: deg2rad
+#ifndef CPP_PHYS
+  REAL, SAVE :: epsfra= 1.E-5
+#endif
+  INTEGER,            SAVE      :: iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn
+  REAL, ALLOCATABLE,  SAVE      :: lon_dyn(:,:), lat_dyn(:,:), levdyn_ini(:)
+  CHARACTER(LEN=120), PARAMETER :: dynfname='ECDYN.nc'
+  CHARACTER(LEN=120), PARAMETER :: orofname='Relief.nc'
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE etat0dyn_netcdf(ib, masque, phis)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Create dynamical initial states.
+!-------------------------------------------------------------------------------
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: read it from ocean model file (coupled run) or compute it.
+!-------------------------------------------------------------------------------
+  USE control_mod
+#ifdef CPP_PHYS
+#ifdef CPP_EARTH
+  USE regr_lat_time_coefoz_m, ONLY: regr_lat_time_coefoz
+  USE regr_pr_o3_m,   ONLY: regr_pr_o3
+  USE press_coefoz_m, ONLY: press_coefoz
+#endif
+#endif
+  USE exner_hyb_m,    ONLY: exner_hyb
+  USE exner_milieu_m, ONLY: exner_milieu
+  USE infotrac
+  USE filtreg_mod
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL, INTENT(IN)    :: ib                  !--- Barycentric interpolation
+  REAL,    INTENT(INOUT) :: masque(iip1,jjp1)   !--- Land-ocean mask
+  REAL,    INTENT(INOUT) :: phis  (iip1,jjp1)   !--- Ground geopotential
+!#ifndef CPP_EARTH
+!  WRITE(lunout,*)'etat0dyn_netcdf: Earth-specific routine, needs Earth physics'
+!#else
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname, fmt
+  INTEGER            :: i, j, l, ji, itau, iday
+  REAL               :: xpn, xps, time, phystep
+  REAL, DIMENSION(iip1,jjp1)       :: psol, masque_tmp
+  REAL, DIMENSION(iip1,jjp1,llm+1) :: p3d
+  REAL, DIMENSION(iip1,jjp1,llm)   :: uvent, t3d, tpot, qsat, qd
+  REAL, DIMENSION(iip1,jjp1,llm)   :: pk, pls, y, masse
+  REAL, DIMENSION(iip1,jjm ,llm)   :: vvent
+  REAL, DIMENSION(ip1jm    ,llm)   :: pbarv
+  REAL, DIMENSION(ip1jmp1  ,llm)   :: pbaru, phi, w
+  REAL, DIMENSION(ip1jmp1)         :: pks
+  REAL, DIMENSION(iim)             :: xppn, xpps
+  REAL, ALLOCATABLE                :: q3d(:,:,:,:)
+!-------------------------------------------------------------------------------
+  modname='etat0dyn_netcdf'
+
+  deg2rad = pi/180.0
+
+! Initializations for tracers and filter
+!*******************************************************************************
+  CALL infotrac_init
+  CALL inifilr()
+
+! Compute ground geopotential and possibly the mask.
+!*******************************************************************************
+  masque_tmp(:,:)=masque(:,:)
+  CALL start_init_orog0(rlonv, rlatu, phis, masque_tmp)
+  WRITE(fmt,"(i4,'i1)')")iip1 ; fmt='('//ADJUSTL(fmt)
+  IF(ALL(masque==-99999.)) THEN                         !--- KEEP NEW MASK
+    masque=masque_tmp
+    IF(prt_level>=1) THEN
+      WRITE(lunout,*)'BUILT MASK :'
+      WRITE(lunout,fmt) NINT(masque)
+    END IF
+    WHERE(   masque(:,:)<EPSFRA) masque(:,:)=0.
+    WHERE(1.-masque(:,:)<EPSFRA) masque(:,:)=1.
+  END IF
+
+! Compute psol AND tsol, knowing phis.
+!*******************************************************************************
+  CALL start_init_dyn(rlonv, rlatu, rlonu, rlatv, ib, phis, psol)
+
+! Mid-levels pressure computation
+!*******************************************************************************
+  CALL pression(ip1jmp1, ap, bp, psol, p3d)             !--- Update p3d
+  IF(pressure_exner) THEN                               !--- Update pk, pks
+    CALL exner_hyb   (ip1jmp1,psol,p3d,pks,pk)
+  ELSE
+    CALL exner_milieu(ip1jmp1,psol,p3d,pks,pk)
+  END IF
+  pls(:,:,:)=preff*(pk(:,:,:)/cpp)**(1./kappa)          !--- Update pls
+
+! Update uvent, vvent, t3d and tpot
+!*******************************************************************************
+  uvent(:,:,:) = 0.0 ; vvent(:,:,:) = 0.0 ; t3d (:,:,:) = 0.0
+  CALL startget_dyn3d('u'   ,rlonu,rlatu,pls,y ,uvent,rlonv,rlatv,ib)
+  CALL startget_dyn3d('v'   ,rlonv,rlatv,pls(:,:jjm,:),y(:,:jjm,:),vvent,      &
+ &                           rlonu,rlatu(:jjm),ib)
+  CALL startget_dyn3d('t'   ,rlonv,rlatu,pls,y ,t3d ,rlonu,rlatv,ib)
+  tpot(:,:,:)=t3d(:,:,:)
+  CALL startget_dyn3d('tpot',rlonv,rlatu,pls,pk,tpot,rlonu,rlatv,ib)
+
+  WRITE(lunout,*) 'T3D min,max:',MINVAL(t3d(:,:,:)),MAXVAL(t3d(:,:,:))
+  WRITE(lunout,*) 'PLS min,max:',MINVAL(pls(:,:,:)),MAXVAL(pls(:,:,:))
+
+! Humidity at saturation computation
+!*******************************************************************************
+  WRITE(lunout,*) 'avant q_sat'
+  CALL q_sat(llm*jjp1*iip1, t3d, pls, qsat)
+  WRITE(lunout,*) 'apres q_sat'
+  WRITE(lunout,*) 'QSAT min,max:',MINVAL(qsat(:,:,:)),MAXVAL(qsat(:,:,:))
+!  WRITE(lunout,*) 'QSAT :',qsat(10,20,:)
+  qd (:,:,:) = 0.0
+  CALL startget_dyn3d('q',rlonv,rlatu,pls,qsat,qd,rlonu,rlatv,ib)
+  ALLOCATE(q3d(iip1,jjp1,llm,nqtot)); q3d(:,:,:,:)=0.0 ; q3d(:,:,:,1)=qd(:,:,:)
+  CALL flinclo(fid_dyn)
+
+#ifdef CPP_PHYS
+#ifdef CPP_EARTH
+! Parameterization of ozone chemistry:
+!*******************************************************************************
+! Look for ozone tracer:
+  DO i=1,nqtot; IF(ANY(["O3","o3"]==tname(i))) EXIT; END DO
+  IF(i/=nqtot+1) THEN
+    CALL regr_lat_time_coefoz
+    CALL press_coefoz
+    CALL regr_pr_o3(p3d, q3d(:,:,:,i))
+    q3d(:,:,:,i)=q3d(:,:,:,i)*48./ 29.                  !--- Mole->mass fraction         
+  END IF
+
+#endif
+#endif
+  q3d(iip1,:,:,:)=q3d(1,:,:,:)
+
+! Intermediate computation
+!*******************************************************************************
+  CALL massdair(p3d,masse)
+  WRITE(lunout,*)' ALPHAX ',alphax
+  DO l=1,llm
+    xppn(:)=aire(1:iim,1   )*masse(1:iim,1   ,l)
+    xpps(:)=aire(1:iim,jjp1)*masse(1:iim,jjp1,l)
+    xpn=SUM(xppn)/apoln
+    xps=SUM(xpps)/apols
+    masse(:,1   ,l)=xpn
+    masse(:,jjp1,l)=xps
+  END DO
+
+! Writing
+!*******************************************************************************
+  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, tetagdiv, tetagrot,   &
+       tetatemp, vert_prof_dissip)
+  WRITE(lunout,*)'sortie inidissip'
+  itau=0
+  itau_dyn=0
+  itau_phy=0
+  iday=dayref+itau/day_step
+  time=FLOAT(itau-(iday-dayref)*day_step)/day_step
+  IF(time>1.) THEN
+   time=time-1
+   iday=iday+1
+  END IF
+  day_ref=dayref
+  annee_ref=anneeref
+  CALL geopot( ip1jmp1, tpot, pk, pks, phis, phi )
+  WRITE(lunout,*)'sortie geopot'
+  CALL caldyn0( itau, uvent, vvent, tpot, psol, masse, pk, phis,               &
+                phi,  w, pbaru, pbarv, time+iday-dayref)
+  WRITE(lunout,*)'sortie caldyn0'     
+  CALL dynredem0( "start.nc", dayref, phis)
+  WRITE(lunout,*)'sortie dynredem0'
+  CALL dynredem1( "start.nc", 0.0, vvent, uvent, tpot, q3d, masse, psol)
+  WRITE(lunout,*)'sortie dynredem1' 
+  CALL histclo()
+
+!#endif 
+! of ifdef CPP_EARTH
+
+END SUBROUTINE etat0dyn_netcdf
+!
+!-------------------------------------------------------------------------------
+
+
+!#ifdef CPP_EARTH
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE startget_dyn3d(var, lon_in,  lat_in,  pls,  workvar,&
+                        champ, lon_in2, lat_in2, ibar)
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose: Compute some quantities (u,v,t,q,tpot) using variables U,V,TEMP and R
+!     (3D fields) of file dynfname.
+!-------------------------------------------------------------------------------
+! Note: An input auxilliary field "workvar" has to be specified in two cases:
+!     * for "q":    the saturated humidity.
+!     * for "tpot": the Exner function.
+!===============================================================================
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)    :: var
+  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)    :: lon_in2(:)       ! dim (iml)
+  REAL,             INTENT(IN)    :: lat_in2(:)       ! dim (jml2)
+  LOGICAL,          INTENT(IN)    :: ibar
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=10)  :: vname
+  CHARACTER(LEN=256) :: msg, modname="startget_dyn3d"
+  INTEGER            :: iml, jml, jml2, lml, il
+  REAL               :: xppn, xpps
+!-------------------------------------------------------------------------------
+  iml=assert_eq([SIZE(lon_in),SIZE(pls,1),SIZE(workvar,1),SIZE(champ,1),       &
+     &                                    SIZE(lon_in2)], TRIM(modname)//" iml")
+  jml=assert_eq( SIZE(lat_in),SIZE(pls,2),SIZE(workvar,2),SIZE(champ,2),       &
+     &                                                    TRIM(modname)//" jml")
+  lml=assert_eq(              SIZE(pls,3),SIZE(workvar,3),SIZE(champ,3),       &
+     &                                                    TRIM(modname)//" lml")
+  jml2=SIZE(lat_in2)
+
+!--- CHECK IF THE FIELD IS KNOWN
+   SELECT CASE(var)
+    CASE('u');    vname='U'
+    CASE('v');    vname='V'
+    CASE('t');    vname='TEMP'
+    CASE('q');    vname='R';    msg='humidity as the saturated humidity'
+    CASE('tpot'); msg='potential temperature as the Exner function'
+    CASE DEFAULT; msg='No rule to extract variable '//TRIM(var)
+      CALL abort_gcm(modname,TRIM(msg)//' from any data set',1)
+  END SELECT
+
+!--- CHECK IF SOMETHING IS MISSING
+  IF((var=='tpot'.OR.var=='q').AND.MINVAL(workvar)==MAXVAL(workvar)) THEN
+    msg='Could not compute '//TRIM(msg)//' is missing or constant.'
+    CALL abort_gcm(modname,TRIM(msg),1)
+  END IF
+
+!--- INTERPOLATE 3D FIELD IF NEEDED
+  IF(var/='tpot') CALL start_inter_3d(TRIM(vname),lon_in,lat_in,lon_in2,      &
+                                                  lat_in2,pls,champ,ibar)
+
+!--- COMPUTE THE REQUIRED FILED
+  SELECT CASE(var)
+    CASE('u'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cu(:,1:jml); END DO
+      champ(iml,:,:)=champ(1,:,:)                   !--- Eastward wind
+
+    CASE('v'); DO il=1,lml; champ(:,:,il)=champ(:,:,il)*cv(:,1:jml); END DO
+      champ(iml,:,:)=champ(1,:,:)                   !--- Northward wind
+
+    CASE('tpot','q')
+      IF(var=='tpot') THEN; champ=champ*cpp/workvar !--- Potential temperature
+      ELSE;                 champ=champ*.01*workvar !--- Relative humidity
+        WHERE(champ<0.) champ=1.0E-10
+      END IF
+      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
+  END SELECT
+
+END SUBROUTINE startget_dyn3d
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_orog0(lon_in,lat_in,phis,masque)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat2d
+  IMPLICIT NONE
+!===============================================================================
+! Purpose:  Compute "phis" just like it would be in start_init_orog.
+!===============================================================================
+! Arguments:
+  REAL, INTENT(IN)    :: lon_in(:), lat_in(:)   ! dim (iml) (jml)
+  REAL, INTENT(INOUT) :: phis(:,:), masque(:,:) ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="start_init_orog0"
+  CHARACTER(LEN=256) :: title="RELIEF"
+  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
+  REAL               :: lev(1), date, dt
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), lon_rel(:,:), relief_hi(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:), lat_rel(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(masque,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),SIZE(masque,2),TRIM(modname)//" jml")
+  IF(iml/=iip1) CALL abort_gcm(TRIM(modname),'iml/=iip1',1)
+  IF(jml/=jjp1) CALL abort_gcm(TRIM(modname),'jml/=jjp1',1)
+  pi=2.0*ASIN(1.0); deg2rad=pi/180.0
+  IF(ANY(phis/=-99999.)) RETURN                  !--- phis ALREADY KNOWN
+
+!--- 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, lon_ini, lat_ini, lon_rad, lat_rad, relief_hi, .FALSE.)
+  DEALLOCATE(lon_ini,lat_ini)
+
+!--- COMPUTING SURFACE GEOPOTENTIAL USING ROUTINE grid_noro0
+  WRITE(lunout,*)
+  WRITE(lunout,*)'*** Compute surface geopotential ***'
+
+!--- CALL OROGRAPHY MODULE (REDUCED VERSION) TO COMPUTE FIELDS
+  CALL grid_noro0(lon_rad, lat_rad, relief_hi, lon_in, lat_in, phis, masque)
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+  DEALLOCATE(relief_hi,lon_rad,lat_rad)
+
+END SUBROUTINE start_init_orog0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE grid_noro0(xd,yd,zd,x,y,zphi,mask)
+!
+!===============================================================================
+! Purpose: Extracted from grid_noro to provide geopotential height for dynamics
+!          without any call to physics subroutines.
+!===============================================================================
+  IMPLICIT NONE 
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL, INTENT(IN)   :: xd(:), yd(:) !--- INPUT  COORDINATES     (imdp) (jmdp)
+  REAL, INTENT(IN)   :: zd(:,:)      !--- INPUT  FIELD           (imdp,jmdp)
+  REAL, INTENT(IN)   :: x(:), y(:)   !--- OUTPUT COORDINATES     (imar+1) (jmar)
+  REAL, INTENT(OUT)  :: zphi(:,:)    !--- GEOPOTENTIAL           (imar+1,jmar)
+  REAL, INTENT(INOUT):: mask(:,:)    !--- MASK                   (imar+1,jmar)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="grid_noro0"
+  REAL, ALLOCATABLE :: xusn(:), yusn(:)           ! dim (imdp+2*iext) (jmdp+2)
+  REAL, ALLOCATABLE :: zusn(:,:)                  ! dim (imdp+2*iext,jmdp+2)
+  REAL, ALLOCATABLE :: weight(:,:)                ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: mask_tmp(:,:), zmea(:,:)   ! dim (imar+1,jmar)
+  REAL, ALLOCATABLE :: num_tot(:,:), num_lan(:,:) ! dim (imax,jmax)
+  REAL, ALLOCATABLE :: a(:), b(:)                 ! dim (imax)
+  REAL, ALLOCATABLE :: c(:), d(:)                 ! dim (jmax)
+  LOGICAL :: masque_lu
+  INTEGER :: i, ii, imdp, imar, iext
+  INTEGER :: j, jj, jmdp, jmar, nn
+  REAL    :: xpi, zlenx, weighx, xincr,  zbordnor, zmeanor, zweinor, zbordest
+  REAL    :: rad, zleny, weighy, masque, zbordsud, zmeasud, zweisud, zbordoue
+!-------------------------------------------------------------------------------
+  imdp=assert_eq(SIZE(xd),SIZE(zd,1),TRIM(modname)//" imdp")
+  jmdp=assert_eq(SIZE(yd),SIZE(zd,2),TRIM(modname)//" jmdp")
+  imar=assert_eq(SIZE(x),SIZE(zphi,1),SIZE(mask,1),TRIM(modname)//" imar")-1
+  jmar=assert_eq(SIZE(y),SIZE(zphi,2),SIZE(mask,2),TRIM(modname)//" jmar")
+  IF(imar/=iim)   CALL abort_gcm(TRIM(modname),'imar/=iim'  ,1)
+  IF(jmar/=jjm+1) CALL abort_gcm(TRIM(modname),'jmar/=jjm+1',1)
+  iext=imdp/10
+  xpi = ACOS(-1.)
+  rad = 6371229.
+
+!--- ARE WE USING A READ MASK ?
+  masque_lu=ANY(mask/=-99999.); IF(.NOT.masque_lu) mask=0.0
+  WRITE(lunout,*)'Masque lu: ',masque_lu
+
+!--- EXTENSION OF THE INPUT DATABASE TO PROCEED COMPUTATIONS AT BOUNDARIES:
+  ALLOCATE(xusn(imdp+2*iext))
+  xusn(1     +iext:imdp  +iext)=xd(:)
+  xusn(1          :       iext)=xd(1+imdp-iext:imdp)-2.*xpi
+  xusn(1+imdp+iext:imdp+2*iext)=xd(1          :iext)+2.*xpi
+
+  ALLOCATE(yusn(jmdp+2))
+  yusn(1       )=yd(1)   +(yd(1)   -yd(2))
+  yusn(2:jmdp+1)=yd(:)
+  yusn(  jmdp+2)=yd(jmdp)+(yd(jmdp)-yd(jmdp-1))
+
+  ALLOCATE(zusn(imdp+2*iext,jmdp+2))
+  zusn(1       +iext:imdp  +iext,2:jmdp+1)=zd  (:                   ,     :)
+  zusn(1            :       iext,2:jmdp+1)=zd  (imdp-iext+1:imdp    ,     :)
+  zusn(1+imdp  +iext:imdp+2*iext,2:jmdp+1)=zd  (1:iext              ,     :)
+  zusn(1            :imdp/2+iext,       1)=zusn(1+imdp/2:imdp  +iext,     2)
+  zusn(1+imdp/2+iext:imdp+2*iext,       1)=zusn(1       :imdp/2+iext,     2)
+  zusn(1            :imdp/2+iext,  jmdp+2)=zusn(1+imdp/2:imdp  +iext,jmdp+1)
+  zusn(1+imdp/2+iext:imdp+2*iext,  jmdp+2)=zusn(1       :imdp/2+iext,jmdp+1)
+
+!--- COMPUTE LIMITS OF MODEL GRIDPOINT AREA (REGULAR GRID)
+  ALLOCATE(a(imar+1),b(imar+1))
+  b(1:imar)=(x(1:imar  )+ x(2:imar+1))/2.0
+  b(imar+1)= x(  imar+1)+(x(  imar+1)-x(imar))/2.0
+  a(1)=x(1)-(x(2)-x(1))/2.0
+  a(2:imar+1)= b(1:imar)
+
+  ALLOCATE(c(jmar),d(jmar))
+  d(1:jmar-1)=(y(1:jmar-1)+ y(2:jmar))/2.0
+  d(  jmar  )= y(  jmar  )+(y(  jmar)-y(jmar-1))/2.0
+  c(1)=y(1)-(y(2)-y(1))/2.0
+  c(2:jmar)=d(1:jmar-1)
+
+!--- INITIALIZATIONS:
+  ALLOCATE(weight(imar+1,jmar)); weight(:,:)= 0.0
+  ALLOCATE(zmea  (imar+1,jmar)); zmea  (:,:)= 0.0
+
+!--- SUMMATION OVER GRIDPOINT AREA
+  zleny=xpi/REAL(jmdp)*rad
+  xincr=xpi/REAL(jmdp)/2.
+  ALLOCATE(num_tot(imar+1,jmar)); num_tot(:,:)=0.
+  ALLOCATE(num_lan(imar+1,jmar)); num_lan(:,:)=0.
+  DO ii = 1, imar+1
+    DO jj = 1, jmar
+      DO j = 2,jmdp+1 
+        zlenx  =zleny  *COS(yusn(j))
+        zbordnor=(xincr+c(jj)-yusn(j))*rad
+        zbordsud=(xincr-d(jj)+yusn(j))*rad
+        weighy=AMAX1(0.,AMIN1(zbordnor,zbordsud,zleny))
+        IF(weighy/=0) THEN
+          DO i = 2, imdp+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/=0)THEN
+              num_tot(ii,jj)=num_tot(ii,jj)+1.0
+              IF(zusn(i,j)>=1.)num_lan(ii,jj)=num_lan(ii,jj)+1.0
+              weight(ii,jj)=weight(ii,jj)+weighx*weighy
+              zmea  (ii,jj)=zmea  (ii,jj)+zusn(i,j)*weighx*weighy !--- MEAN
+            END IF
+          END DO
+        END IF
+      END DO
+    END DO
+  END DO
+
+!--- COMPUTE PARAMETERS NEEDED BY LOTT & MILLER (1997) AND LOTT (1999) SSO SCHEME
+  IF(.NOT.masque_lu) THEN
+    WHERE(weight(:,1:jmar-1)/=0.0) mask=num_lan(:,:)/num_tot(:,:)
+  END IF
+  nn=COUNT(weight(:,1:jmar-1)==0.0)
+  IF(nn/=0) WRITE(lunout,*)'Problem with weight ; vanishing occurrences: ',nn
+  WHERE(weight/=0.0) zmea(:,:)=zmea(:,:)/weight(:,:)
+
+!--- MASK BASED ON GROUND MAXIMUM, 10% THRESHOLD (<10%: SURF PARAMS MEANINGLESS)
+  ALLOCATE(mask_tmp(imar+1,jmar)); mask_tmp(:,:)=0.0
+  WHERE(mask>=0.1) mask_tmp = 1.
+  WHERE(weight(:,:)/=0.0)
+    zphi(:,:)=mask_tmp(:,:)*zmea(:,:)
+    zmea(:,:)=mask_tmp(:,:)*zmea(:,:)
+  END WHERE
+  WRITE(lunout,*)'  MEAN ORO:' ,MAXVAL(zmea)
+
+!--- Values at poles
+  zphi(imar+1,:)=zphi(1,:)
+
+  zweinor=SUM(weight(1:imar,   1),DIM=1)
+  zweisud=SUM(weight(1:imar,jmar),DIM=1)
+  zmeanor=SUM(weight(1:imar,   1)*zmea(1:imar,   1),DIM=1)
+  zmeasud=SUM(weight(1:imar,jmar)*zmea(1:imar,jmar),DIM=1)
+  zphi(:,1)=zmeanor/zweinor; zphi(:,jmar)=zmeasud/zweisud
+
+END SUBROUTINE grid_noro0
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_dyn(lon_in,lat_in,lon_in2,lat_in2,ibar,zs,psol)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!===============================================================================
+! Purpose:   Compute psol, knowing phis.
+!===============================================================================
+! Arguments:
+  REAL,    INTENT(IN)  :: lon_in (:),  lat_in (:)    ! dim (iml) (jml)
+  REAL,    INTENT(IN)  :: lon_in2(:),  lat_in2(:)    ! dim (iml) (jml2)
+  LOGICAL, INTENT(IN)  :: ibar
+  REAL,    INTENT(IN)  :: zs  (:,:)                  ! dim (iml,jml)
+  REAL,    INTENT(OUT) :: psol(:,:)                  ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname='start_init_dyn'
+  REAL               :: date, dt
+  INTEGER            :: iml, jml, jml2, itau(1)
+  REAL, ALLOCATABLE  :: lon_rad(:), lon_ini(:), var_ana(:,:)
+  REAL, ALLOCATABLE  :: lat_rad(:), lat_ini(:)
+  REAL, ALLOCATABLE  :: z(:,:), ps(:,:), ts(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(zs,1),SIZE(psol,1),SIZE(lon_in2),            &
+      &                                              TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(zs,2),SIZE(psol,2),TRIM(modname)//" jml")
+  jml2=SIZE(lat_in2)
+
+  WRITE(lunout,*) 'Opening the surface analysis'
+  CALL flininfo(dynfname, iml_dyn, jml_dyn, llm_dyn, ttm_dyn, fid_dyn)
+  WRITE(lunout,*) 'Values read: ', iml_dyn, jml_dyn, llm_dyn, ttm_dyn
+
+  ALLOCATE(lon_dyn(iml_dyn,jml_dyn), lat_dyn(iml_dyn,jml_dyn))
+  ALLOCATE(levdyn_ini(llm_dyn))
+  CALL flinopen(dynfname, .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))
+  CALL get_var_dyn('Z',z)                        !--- SURFACE GEOPOTENTIAL
+  CALL get_var_dyn('SP',ps)                      !--- SURFACE PRESSURE
+  CALL get_var_dyn('ST',ts)                      !--- SURFACE TEMPERATURE
+!  CALL flinclo(fid_dyn)
+  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
+
+!--- PSOL IS COMPUTED IN PASCALS
+  psol(:iml-1,:) = ps(:iml-1,:)*(1.0+(z(:iml-1,:)-zs(:iml-1,:))/287.0          &
+     &            /ts(:iml-1,:))
+  psol(iml,:)=psol(1,:)
+  DEALLOCATE(z,ps,ts)
+  psol(:,1  )=SUM(aire(1:iml-1,1  )*psol(1:iml-1,1  ))/apoln  !--- NORTH POLE
+  psol(:,jml)=SUM(aire(1:iml-1,jml)*psol(1:iml-1,jml))/apols  !--- SOUTH POLE
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_var_dyn(title,field)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat2d
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)    :: title
+  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: msg
+  INTEGER :: tllm
+!-------------------------------------------------------------------------------
+  SELECT CASE(title)
+    CASE('Z');     tllm=0;       msg='geopotential'
+    CASE('SP');    tllm=0;       msg='surface pressure'
+    CASE('ST');    tllm=llm_dyn; msg='temperature'
+  END SELECT
+  IF(.NOT.ALLOCATED(field)) THEN
+    ALLOCATE(field(iml,jml))
+    CALL flinget(fid_dyn, title, iml_dyn,jml_dyn, tllm, ttm_dyn, 1, 1, var_ana)
+    CALL conf_dat2d(title, lon_ini, lat_ini, lon_rad, lat_rad, var_ana, ibar)
+    CALL interp_startvar(title, ibar, .TRUE., lon_rad, lat_rad, var_ana,       &
+                              lon_in, lat_in, lon_in2, lat_in2, field)
+  ELSE IF(SIZE(field)/=SIZE(z)) THEN
+    msg='The '//TRIM(msg)//' field we have does not have the right size'
+    CALL abort_gcm(TRIM(modname),msg,1)
+  END IF
+
+END SUBROUTINE get_var_dyn
+!
+!-------------------------------------------------------------------------------
+
+END SUBROUTINE start_init_dyn
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_inter_3d(var,lon_in,lat_in,lon_in2,lat_in2,pls_in,var3d,ibar)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat3d
+  USE pchsp_95_m, ONLY: pchsp_95
+  USE pchfe_95_m, ONLY: pchfe_95
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN) :: var
+  REAL,    INTENT(IN)  :: lon_in(:),  lat_in(:)   ! dim (iml) (jml)
+  REAL,    INTENT(IN)  :: lon_in2(:), lat_in2(:)  ! dim (iml) (jml2)
+  REAL,    INTENT(IN)  :: pls_in(:,:,:)           ! dim (iml,jml,lml)
+  REAL,    INTENT(OUT) :: var3d (:,:,:)           ! dim (iml,jml,lml)
+  LOGICAL, INTENT(IN)  :: ibar
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname='start_inter_3d'
+  LOGICAL :: skip
+  REAL    :: chmin, chmax
+  INTEGER :: iml, jml, lml, jml2, ii, ij, il, ierr
+  INTEGER :: n_extrap                             ! Extrapolated points number
+  REAL, ALLOCATABLE :: ax(:), lon_rad(:), lon_ini(:), lev_dyn(:), yder(:)
+  REAL, ALLOCATABLE :: ay(:), lat_rad(:), lat_ini(:), var_tmp3d(:,:,:)
+  REAL, ALLOCATABLE, SAVE :: var_ana3d(:,:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(lon_in2),SIZE(pls_in,1),SIZE(var3d,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),              SIZE(pls_in,2),SIZE(var3d,2),TRIM(modname)//" jml")
+  lml=assert_eq(SIZE(pls_in,3),SIZE(var3d,3),TRIM(modname)//" lml"); jml2=SIZE(lat_in2)
+
+  WRITE(lunout, *)'Going into flinget to extract the 3D field.'
+  IF(.NOT.ALLOCATED(var_ana3d)) ALLOCATE(var_ana3d(iml_dyn, jml_dyn, llm_dyn))
+  CALL flinget(fid_dyn,var,iml_dyn,jml_dyn,llm_dyn,ttm_dyn,1,1,var_ana3d)
+
+!--- ANGLES IN DEGREES 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), lev_dyn(llm_dyn))
+  CALL conf_dat3d(var, lon_ini, lat_ini, levdyn_ini,                           &
+                       lon_rad, lat_rad, lev_dyn, var_ana3d, ibar)
+  DEALLOCATE(lon_ini, lat_ini)
+
+!--- COMPUTE THE REQUIRED FIELDS USING ROUTINE grid_noro
+  ALLOCATE(var_tmp3d(iml,jml,llm_dyn))
+  DO il = 1,llm_dyn
+    CALL interp_startvar(var,ibar,il==1,lon_rad,lat_rad,var_ana3d(:,:,il),     &
+                          lon_in,lat_in,lon_in2,lat_in2,var_tmp3d(:,:,il))
+  END DO
+  DEALLOCATE(lon_rad, lat_rad)
+
+!--- VERTICAL INTERPOLATION FROM TOP OF ATMOSPHERE TO GROUND
+  ALLOCATE(ax(llm_dyn),ay(llm_dyn),yder(llm_dyn))
+  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) CALL abort_gcm(TRIM(modname),'error in pchfe_95',1)
+      n_extrap = n_extrap + ierr
+    END DO
+  END DO
+  IF(n_extrap/=0) WRITE(lunout,*)TRIM(modname)//" pchfe_95: n_extrap=", n_extrap
+  var3d(iml, :, :) = var3d(1, :, :) 
+
+  DO il=1, lml
+    CALL minmax(iml*jml, var3d(1, 1, il), chmin, chmax)
+    WRITE(lunout, *)' '//TRIM(var)//'  min max l ', il, chmin, chmax
+  END DO
+
+END SUBROUTINE start_inter_3d
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE interp_startvar(nam,ibar,ibeg,lon,lat,vari,lon1,lat1,lon2,lat2,varo)
+!
+!-------------------------------------------------------------------------------
+  USE inter_barxy_m, ONLY: inter_barxy
+  USE grid_atob_m,   ONLY: grille_m
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)  :: nam
+  LOGICAL,          INTENT(IN)  :: ibar, ibeg
+  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
+  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
+  REAL,             INTENT(IN)  :: lon1(:), lat1(:) ! dim (i1) (j1)
+  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
+  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="interp_startvar"
+  INTEGER            :: ii, jj, i1, j1, j2
+  REAL, ALLOCATABLE  :: vtmp(:,:)
+!-------------------------------------------------------------------------------
+  ii=assert_eq(SIZE(lon),            SIZE(vari,1),TRIM(modname)//" ii")
+  jj=assert_eq(SIZE(lat),            SIZE(vari,2),TRIM(modname)//" jj")
+  i1=assert_eq(SIZE(lon1),SIZE(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=assert_eq(SIZE(lat1),           SIZE(varo,2),TRIM(modname)//" j1")
+  j2=SIZE(lat2)
+  ALLOCATE(vtmp(i1-1,j1))
+  IF(ibar) THEN
+    IF(ibeg.AND.prt_level>1) THEN
+      WRITE(lunout,*)"---------------------------------------------------------"
+      WRITE(lunout,*)"$$$ Interpolation barycentrique pour "//TRIM(nam)//" $$$"
+      WRITE(lunout,*)"---------------------------------------------------------"
+    END IF
+    CALL inter_barxy(lon, lat(:jj-1), vari, lon2(:i1-1), lat2, vtmp)
+  ELSE
+    CALL grille_m   (lon, lat,        vari, 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 etat0dyn
+!
+!*******************************************************************************
