Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/etat0phys_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/etat0phys_netcdf.F90	(revision 2301)
+++ 	(revision )
@@ -1,558 +1,0 @@
-MODULE etat0phys
-!
-!*******************************************************************************
-! Purpose: Create physical 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)
-!
-!    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
-!         "start_init_phys" for variables contained in file "ECPHY.nc":
-!            'ST'     : Surface temperature
-!            'CDSW'   : Soil moisture
-!         "start_init_orog" for variables contained in file "Relief.nc":
-!            'RELIEF' : High resolution orography 
-!
-!    * The land mask and corresponding weights can be:
-!      1) 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"
-!      2) computed from topography file "Relief.nc" for forced runs.
-!
-!    * Allowed values for read_climoz flag are 0, 1 and 2:
-!      0: do not read an ozone climatology
-!      1: read a single ozone climatology that will be used day and night
-!      2: read two ozone climatologies, the average day and night climatology
-!         and the daylight climatology
-!-------------------------------------------------------------------------------
-!    * 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
-  USE assert_eq_m,        ONLY: assert_eq
-#ifdef CPP_EARTH
-  USE dimphy
-  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
-    rlon, solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
-    rlat, sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
-    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
-    zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
-    phys_state_var_init
-#endif
-
-  PRIVATE
-  PUBLIC :: etat0phys_netcdf
-
-  include "iniprint.h"
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-  include "comvert.h"
-  include "comconst.h"
-  include "dimsoil.h"
-  include "temps.h"
-  include "comdissnew.h"
-  include "serre.h"
-  include "clesphys.h"
-  REAL, SAVE :: deg2rad
-  REAL, SAVE, ALLOCATABLE :: tsol(:)
-!  REAL, SAVE, ALLOCATABLE :: rugo(:,:)  ! ??? COMPUTED BUT NOT USED ???
-  INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
-  REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
-  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc"
-  CHARACTER(LEN=256), PARAMETER :: title="RELIEF"
-
-
-CONTAINS
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE etat0phys_netcdf(ib, masque, phis)
-!
-!-------------------------------------------------------------------------------
-! Purpose: Creates initial states
-!-------------------------------------------------------------------------------
-! Note: This routine is designed to work for Earth
-!-------------------------------------------------------------------------------
-  USE control_mod
-#ifdef CPP_EARTH
-  USE infotrac
-  USE fonte_neige_mod
-  USE pbl_surface_mod
-  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
-  USE indice_sol_mod
-  USE conf_phys_m,    ONLY: conf_phys
-  USE exner_hyb_m,    ONLY: exner_hyb
-  USE exner_milieu_m, ONLY: exner_milieu
-  USE test_disvert_m, ONLY: test_disvert
-  USE grid_atob_m,    ONLY: grille_m
-#endif
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  LOGICAL, INTENT(IN)    :: ib          !--- Barycentric interpolation
-  REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
-  REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
-#ifndef CPP_EARTH
-  WRITE(lunout,*)'etat0phys_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
-  INTEGER            :: i, j, l, ji, iml, jml
-  REAL               :: phystep
-  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
-  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
-  REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
-  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
-
-!--- Local variables for sea-ice reading:
-  LOGICAL           :: read_mask
-  INTEGER           :: iml_lic, jml_lic, isst(klon-2)
-  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1)
-  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)
-  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)
-  REAL              :: date, lev(1), dummy
-  REAL              :: flic_tmp(SIZE(masque,1),SIZE(masque,2))
-
-!--- Arguments for conf_phys
-  LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
-  REAL    :: solarlong0, seuil_inversion, fact_cldcon, facttemps
-  LOGICAL :: ok_newmicro
-  INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
-  REAL    :: ratqsbas, ratqshaut, tau_ratqs
-  LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
-  INTEGER :: flag_aerosol
-  LOGICAL :: flag_aerosol_strat
-  LOGICAL :: new_aod
-  REAL    :: bl95_b0, bl95_b1
-  INTEGER :: read_climoz                        !--- Read ozone climatology
-  REAL    :: alp_offset
-
-  deg2rad= pi/180.0
-  iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
-
-! Grid construction and miscellanous initializations.
-!*******************************************************************************
-  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
-                   callstats,                                           &
-                   solarlong0,seuil_inversion,                          &
-                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
-                   iflag_cldcon,                                        &
-                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
-                   flag_aerosol, flag_aerosol_strat, new_aod,           &
-                   bl95_b0, bl95_b1,                                    &
-                   read_climoz,                                         &
-                   alp_offset)
-
-  CALL phys_state_var_init(read_climoz)
-
-  co2_ppm0 = co2_ppm  !--- Initial atmospheric CO2 conc. from .def file
-
-  rlat(1) = pi/2.
-  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
-  rlat(klon) = - pi/2.
-  rlat(:)=rlat(:)*(180.0/pi)
-
-  rlon(1) = 0.0
-  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
-  rlon(klon) = 0.0
-  rlon(:)=rlon(:)*(180.0/pi)
-
-! Compute ground geopotential, sub-cells quantities and possibly the mask.
-!*******************************************************************************
-  read_mask=ANY(masque/=-99999.); masque_tmp=masque
-  CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
-  WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
-  IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
-    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
- CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
-
-! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
-!*******************************************************************************
-  CALL start_init_phys(rlonv, rlatu, rlonu, rlatv, ib, phis)
-
-! Some initializations.
-!*******************************************************************************
-  sn    (:) = 0.0                               !--- Snow
-  radsol(:) = 0.0                               !--- Net radiation at ground
-  rugmer(:) = 0.001                             !--- Ocean rugosity
-  IF(read_climoz>=1) &                          !--- Ozone climatology
-    CALL regr_lat_time_climoz(read_climoz)
-
-! Sub-surfaces initialization
-!*******************************************************************************
-!--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
-  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
-  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
-  ALLOCATE( fraclic(iml_lic,jml_lic))
-  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
- &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
-  CALL flinclo(fid)
-  WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
-  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees
-  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
-  dlon_lic(:)=lon_lic(:,1)
-  dlat_lic(:)=lat_lic(1,:)
-  CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) )
-  flic_tmp(iml,:)=flic_tmp(1,:)
-
-!--- To the physical grid
-  pctsrf(:,:) = 0.
-  CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic))
-
-!--- Adequation with soil/sea mask
-  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
-  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
-  pctsrf(:,is_ter)=zmasq(:)
-  DO ji=1,klon
-    IF(zmasq(ji)>EPSFRA) THEN 
-      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
-        pctsrf(ji,is_lic)=zmasq(ji)
-        pctsrf(ji,is_ter)=0.
-      ELSE
-        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
-        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
-          pctsrf(ji,is_ter)=0.
-          pctsrf(ji,is_lic)=zmasq(ji)
-        END IF 
-      END IF 
-    END IF 
-  END DO 
-
-!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
-  pctsrf(:,is_oce)=(1.-zmasq(:))
-  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
-  IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:)
-  isst=0
-  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
-
-!--- It is checked that the sub-surfaces sum is equal to 1.
-  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
-  IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
-
-! Write physical initial state
-!*******************************************************************************
-  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
-  phystep = dtvr * FLOAT(iphysiq)
-  radpas  = NINT (86400./phystep/ FLOAT(nbapp_rad) )
-  WRITE(lunout,*)'phystep =', phystep, radpas
-
-! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
-!*******************************************************************************
-  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
-  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
-  falb_dir(:,is_ter,:) = 0.08
-  falb_dir(:,is_lic,:) = 0.6
-  falb_dir(:,is_oce,:) = 0.5
-  falb_dir(:,is_sic,:) = 0.6
-  fevap(:,:) = 0.
-  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
-  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
-  rain_fall  = 0.
-  snow_fall  = 0.
-  solsw      = 165.
-  sollw      = -53.
-  t_ancien   = 273.15
-  q_ancien   = 0.
-  agesno     = 0.
-
-  z0m(:,is_oce) = rugmer(:)
-  z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  z0m(:,is_sic) = 0.001
-  z0h(:,:)=z0m(:,:)
-
-  fder    = 0.0
-  clwcon  = 0.0
-  rnebcon = 0.0
-  ratqs   = 0.0
-  run_off_lic_0 = 0.0 
-  rugoro  = 0.0
-
-! Before phyredem calling, surface modules and values to be saved in startphy.nc
-! are initialized
-!*******************************************************************************
-  dummy            = 1.0
-  pbl_tke(:,:,:)   = 1.e-8 
-  zmax0(:)         = 40.
-  f0(:)            = 1.e-5
-  sig1(:,:)        = 0.
-  w01(:,:)         = 0.
-  wake_deltat(:,:) = 0.
-  wake_deltaq(:,:) = 0.
-  wake_s(:)        = 0.
-  wake_cstar(:)    = 0.
-  wake_fip(:)      = 0.
-  wake_pe          = 0.
-  fm_therm         = 0.
-  entr_therm       = 0.
-  detr_therm       = 0.
-
-  CALL fonte_neige_init(run_off_lic_0)
-  CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil )
-  CALL phyredem( "startphy.nc" )
-
-!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
-!  WRITE(lunout,*)'entree histclo'
-  CALL histclo()
-
-#endif 
-!#endif of #ifdef CPP_EARTH
-
-END SUBROUTINE etat0phys_netcdf
-!
-!-------------------------------------------------------------------------------
-
-
-#ifdef CPP_EARTH
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
-!
-!===============================================================================
-! Comment:
-!   This routine launch grid_noro, which computes parameters for SSO scheme as
-!   described in LOTT & MILLER (1997) and LOTT(1999).
-!===============================================================================
-  USE conf_dat_m,  ONLY: conf_dat2d
-!  USE grid_atob_m, ONLY: rugsoro
-  USE grid_noro_m, ONLY: grid_noro
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! 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_orog"
-  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(:,:), tmp_var  (:,:)
-  REAL, ALLOCATABLE  :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
-  REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
-!-------------------------------------------------------------------------------
-  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")
-
-!--- 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 THE REQUIRED FIELDS USING ROUTINE grid_noro
-  WRITE(lunout,*)
-  WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***'
-
-!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
-  ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
-  ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
-  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
-  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
-
-!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
-  CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0,     &
-                                      zsig0,zgam0,zthe0,zpic0,zval0,masque)
-  phis = phis * 9.81
-  phis(iml,:) = phis(1,:)
-
-!--- COMPUTE SURFACE ROUGHNESS
-!  WRITE(lunout,*)
-!  WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'
-!  ALLOCATE(tmp_var(iml-1,jml))
-!  CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)
-!  ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
-!  DEALLOCATE(tmp_var)
-  DEALLOCATE(relief_hi,lon_rad,lat_rad)
-
-!--- PUT QUANTITIES TO PHYSICAL GRID
-  CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
-  CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
-
-
-END SUBROUTINE start_init_orog
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis)
-!
-!===============================================================================
-! Purpose:   Compute tsol and qsol, knowing phis.
-!===============================================================================
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! 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) :: phis(:,:)                   ! dim (iml,jml)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc"
-  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  :: ts(:,:), qs(:,:)
-!-------------------------------------------------------------------------------
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")
-  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),              TRIM(modname)//" jml")
-  jml2=SIZE(lat_in2)
-
-  WRITE(lunout,*)'Opening the surface analysis'
-  CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
-  WRITE(lunout,*) 'Values read: ',   iml_phys, jml_phys, llm_phys, ttm_phys
-
-  ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
-  ALLOCATE(levphys_ini(llm_phys))
-  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
-                lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
-
-!--- 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))
-  CALL get_var_phys('ST'  ,ts)                   !--- SURFACE TEMPERATURE
-  CALL get_var_phys('CDSW',qs)                   !--- SOIL MOISTURE
-  CALL flinclo(fid_phys)
-  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
-
-!--- TSOL AND QSOL ON PHYSICAL GRID
-  ALLOCATE(tsol(klon))
-  CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
-  CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
-  DEALLOCATE(ts,qs)
-
-CONTAINS
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE get_var_phys(title,field)
-!
-!-------------------------------------------------------------------------------
-  USE conf_dat_m, ONLY: conf_dat2d
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),  INTENT(IN)    :: title
-  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
-!-------------------------------------------------------------------------------
-! Local variables:
-  INTEGER :: tllm
-!-------------------------------------------------------------------------------
-  SELECT CASE(title)
-    CASE('SP');        tllm=0
-    CASE('ST','CDSW'); tllm=llm_phys
-  END SELECT
-  IF(ALLOCATED(field)) RETURN
-  ALLOCATE(field(iml,jml)); field(:,:)=0.
-  CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,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)
-
-END SUBROUTINE get_var_phys
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE start_init_phys
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-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 
-!#endif of #ifdef CPP_EARTH
-
-END MODULE etat0phys
-!
-!*******************************************************************************
-
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 2302)
@@ -77,9 +77,7 @@
   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
@@ -93,7 +91,4 @@
   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:
@@ -176,5 +171,4 @@
 
 #ifdef CPP_PHYS
-#ifdef CPP_EARTH
 ! Parameterization of ozone chemistry:
 !*******************************************************************************
@@ -188,5 +182,4 @@
   END IF
 
-#endif
 #endif
   q3d(iip1,:,:,:)=q3d(1,:,:,:)
@@ -232,7 +225,4 @@
   CALL histclo()
 
-!#endif 
-! of ifdef CPP_EARTH
-
 END SUBROUTINE etat0dyn_netcdf
 !
@@ -240,5 +230,4 @@
 
 
-!#ifdef CPP_EARTH
 !-------------------------------------------------------------------------------
 !
@@ -763,7 +752,4 @@
 !-------------------------------------------------------------------------------
 
-!#endif
-! of ifdef CPP_EARTH
-
 END MODULE etat0dyn
 !
Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2302)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/etat0phys_netcdf.F90	(revision 2302)
@@ -0,0 +1,545 @@
+MODULE etat0phys
+!
+!*******************************************************************************
+! Purpose: Create physical 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)
+!
+!    *  etat0phys_netcdf routine can access to NetCDF data through subroutines:
+!         "start_init_phys" for variables contained in file "ECPHY.nc":
+!            'ST'     : Surface temperature
+!            'CDSW'   : Soil moisture
+!         "start_init_orog" for variables contained in file "Relief.nc":
+!            'RELIEF' : High resolution orography 
+!
+!    * The land mask and corresponding weights can be:
+!      1) 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"
+!      2) computed from topography file "Relief.nc" for forced runs.
+!
+!    * Allowed values for read_climoz flag are 0, 1 and 2:
+!      0: do not read an ozone climatology
+!      1: read a single ozone climatology that will be used day and night
+!      2: read two ozone climatologies, the average day and night climatology
+!         and the daylight climatology
+!-------------------------------------------------------------------------------
+!    * 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
+  USE assert_eq_m,        ONLY: assert_eq
+  USE dimphy
+  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
+    rlon, solsw, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
+    rlat, sollw, rugoro, q_ancien, wake_deltaq, wake_pe, snow_fall, ratqs,w01, &
+    sig1, ftsol, clwcon, fm_therm, wake_Cstar,  pctsrf,  entr_therm,radpas, f0,&
+    zmax0,fevap, rnebcon,falb_dir, wake_fip,    agesno,  detr_therm, pbl_tke,  &
+    phys_state_var_init
+
+  PRIVATE
+  PUBLIC :: etat0phys_netcdf
+
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "comvert.h"
+  include "comconst.h"
+  include "dimsoil.h"
+  include "temps.h"
+  include "comdissnew.h"
+  include "serre.h"
+  include "clesphys.h"
+  REAL, SAVE :: deg2rad
+  REAL, SAVE, ALLOCATABLE :: tsol(:)
+!  REAL, SAVE, ALLOCATABLE :: rugo(:,:)  ! ??? COMPUTED BUT NOT USED ???
+  INTEGER,            SAVE      :: iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys
+  REAL, ALLOCATABLE,  SAVE      :: lon_phys(:,:), lat_phys(:,:), levphys_ini(:)
+  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc"
+  CHARACTER(LEN=256), PARAMETER :: title="RELIEF"
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE etat0phys_netcdf(ib, masque, phis)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Creates initial states
+!-------------------------------------------------------------------------------
+! Note: This routine is designed to work for Earth
+!-------------------------------------------------------------------------------
+  USE control_mod
+  USE infotrac
+  USE fonte_neige_mod
+  USE pbl_surface_mod
+  USE regr_lat_time_climoz_m, ONLY: regr_lat_time_climoz
+  USE indice_sol_mod
+  USE conf_phys_m,    ONLY: conf_phys
+  USE exner_hyb_m,    ONLY: exner_hyb
+  USE exner_milieu_m, ONLY: exner_milieu
+  USE test_disvert_m, ONLY: test_disvert
+  USE grid_atob_m,    ONLY: grille_m
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL, INTENT(IN)    :: ib          !--- Barycentric interpolation
+  REAL,    INTENT(INOUT) :: masque(:,:) !--- Land mask           dim(iip1,jjp1)
+  REAL,    INTENT(INOUT) :: phis  (:,:) !--- Ground geopotential dim(iip1,jjp1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="etat0phys_netcdf", fmt
+  INTEGER            :: i, j, l, ji, iml, jml
+  REAL               :: phystep
+  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp
+  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
+  REAL, DIMENSION(klon,nbsrf)         :: qsolsrf, snsrf
+  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
+
+!--- Local variables for sea-ice reading:
+  LOGICAL           :: read_mask
+  INTEGER           :: iml_lic, jml_lic, isst(klon-2)
+  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1)
+  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic(:,:)
+  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:)
+  REAL              :: date, lev(1), dummy
+  REAL              :: flic_tmp(SIZE(masque,1),SIZE(masque,2))
+
+!--- Arguments for conf_phys
+  LOGICAL :: ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES, callstats
+  REAL    :: solarlong0, seuil_inversion, fact_cldcon, facttemps
+  LOGICAL :: ok_newmicro
+  INTEGER :: iflag_radia, iflag_cldcon, iflag_ratqs
+  REAL    :: ratqsbas, ratqshaut, tau_ratqs
+  LOGICAL :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
+  INTEGER :: flag_aerosol
+  LOGICAL :: flag_aerosol_strat
+  LOGICAL :: new_aod
+  REAL    :: bl95_b0, bl95_b1
+  INTEGER :: read_climoz                        !--- Read ozone climatology
+  REAL    :: alp_offset
+
+  deg2rad= pi/180.0
+  iml=assert_eq(SIZE(masque,1),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(masque,2),SIZE(phis,2),TRIM(modname)//" jml")
+
+! Grid construction and miscellanous initializations.
+!*******************************************************************************
+  CALL conf_phys(  ok_journe, ok_mensuel, ok_instan, ok_hf, ok_LES,     &
+                   callstats,                                           &
+                   solarlong0,seuil_inversion,                          &
+                   fact_cldcon, facttemps,ok_newmicro,iflag_radia,      &
+                   iflag_cldcon,                                        &
+                   iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
+                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
+                   flag_aerosol, flag_aerosol_strat, new_aod,           &
+                   bl95_b0, bl95_b1,                                    &
+                   read_climoz,                                         &
+                   alp_offset)
+
+  CALL phys_state_var_init(read_climoz)
+
+  co2_ppm0 = co2_ppm  !--- Initial atmospheric CO2 conc. from .def file
+
+  rlat(1) = pi/2.
+  DO j=2,jjm; rlat((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);     END DO
+  rlat(klon) = - pi/2.
+  rlat(:)=rlat(:)*(180.0/pi)
+
+  rlon(1) = 0.0
+  DO j=2,jjm; rlon((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:iim); END DO
+  rlon(klon) = 0.0
+  rlon(:)=rlon(:)*(180.0/pi)
+
+! Compute ground geopotential, sub-cells quantities and possibly the mask.
+!*******************************************************************************
+  read_mask=ANY(masque/=-99999.); masque_tmp=masque
+  CALL start_init_orog(rlonv, rlatu, phis, masque_tmp)
+  WRITE(fmt,"(i4,'i1)')")iml ; fmt='('//ADJUSTL(fmt)
+  IF(.NOT.read_mask) THEN                       !--- Keep mask form orography
+    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
+ CALL gr_dyn_fi(1,iml,jml,klon,masque,zmasq) !--- Land mask to physical grid
+
+! Compute tsol and qsol on physical grid, knowing phis on 2D grid.
+!*******************************************************************************
+  CALL start_init_phys(rlonv, rlatu, rlonu, rlatv, ib, phis)
+
+! Some initializations.
+!*******************************************************************************
+  sn    (:) = 0.0                               !--- Snow
+  radsol(:) = 0.0                               !--- Net radiation at ground
+  rugmer(:) = 0.001                             !--- Ocean rugosity
+  IF(read_climoz>=1) &                          !--- Ozone climatology
+    CALL regr_lat_time_climoz(read_climoz)
+
+! Sub-surfaces initialization
+!*******************************************************************************
+!--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
+  CALL flininfo("landiceref.nc", iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
+  ALLOCATE( lat_lic(iml_lic,jml_lic),lon_lic(iml_lic, jml_lic))
+  ALLOCATE(dlat_lic(jml_lic),       dlon_lic(iml_lic))
+  ALLOCATE( fraclic(iml_lic,jml_lic))
+  CALL flinopen("landiceref.nc", .FALSE., iml_lic, jml_lic, llm_tmp,  &
+ &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
+  CALL flinget(fid, 'landice', iml_lic, jml_lic, llm_tmp, ttm_tmp, 1,1, fraclic)
+  CALL flinclo(fid)
+  WRITE(lunout,*)'landice dimensions: iml_lic, jml_lic : ',iml_lic,jml_lic
+  IF(MAXVAL(lon_lic)>pi) lon_lic=lon_lic*pi/180. !--- Conversion to degrees
+  IF(MAXVAL(lat_lic)>pi) lat_lic=lat_lic*pi/180.
+  dlon_lic(:)=lon_lic(:,1)
+  dlat_lic(:)=lat_lic(1,:)
+  CALL grille_m(dlon_lic, dlat_lic, fraclic, rlonv(1:iim), rlatu, flic_tmp(1:iim,:) )
+  flic_tmp(iml,:)=flic_tmp(1,:)
+
+!--- To the physical grid
+  pctsrf(:,:) = 0.
+  CALL gr_dyn_fi(1, iml, jml, klon, flic_tmp, pctsrf(:,is_lic))
+
+!--- Adequation with soil/sea mask
+  WHERE(pctsrf(:,is_lic)<EPSFRA) pctsrf(:,is_lic)=0. 
+  WHERE(zmasq(:)<EPSFRA)         pctsrf(:,is_lic)=0.
+  pctsrf(:,is_ter)=zmasq(:)
+  DO ji=1,klon
+    IF(zmasq(ji)>EPSFRA) THEN 
+      IF(pctsrf(ji,is_lic)>=zmasq(ji)) THEN
+        pctsrf(ji,is_lic)=zmasq(ji)
+        pctsrf(ji,is_ter)=0.
+      ELSE
+        pctsrf(ji,is_ter)=zmasq(ji)-pctsrf(ji,is_lic)
+        IF(pctsrf(ji,is_ter)<EPSFRA) THEN
+          pctsrf(ji,is_ter)=0.
+          pctsrf(ji,is_lic)=zmasq(ji)
+        END IF 
+      END IF 
+    END IF 
+  END DO 
+
+!--- Sub-surface ocean and sea ice (sea ice set to zero for start).
+  pctsrf(:,is_oce)=(1.-zmasq(:))
+  WHERE(pctsrf(:,is_oce)<EPSFRA) pctsrf(:,is_oce)=0.
+  IF(read_mask) pctsrf(:,is_oce)=1-zmasq(:)
+  isst=0
+  WHERE(pctsrf(2:klon-1,is_oce)>0.) isst=1
+
+!--- It is checked that the sub-surfaces sum is equal to 1.
+  ji=COUNT((ABS(SUM(pctsrf(:,:),dim=2))-1.0)>EPSFRA)
+  IF(ji/=0) WRITE(lunout,*) 'Sub-cell distribution problem for ',ji,' points'
+
+! Write physical initial state
+!*******************************************************************************
+  WRITE(lunout,*)'phystep ',dtvr,iphysiq,nbapp_rad
+  phystep = dtvr * FLOAT(iphysiq)
+  radpas  = NINT (86400./phystep/ FLOAT(nbapp_rad) )
+  WRITE(lunout,*)'phystep =', phystep, radpas
+
+! Init: ftsol, snsrf, qsolsrf, tsoil, rain_fall, snow_fall, solsw, sollw, z0
+!*******************************************************************************
+  DO i=1,nbsrf; ftsol(:,i) = tsol; END DO
+  DO i=1,nbsrf; snsrf(:,i) = sn;   END DO
+  falb_dir(:,is_ter,:) = 0.08
+  falb_dir(:,is_lic,:) = 0.6
+  falb_dir(:,is_oce,:) = 0.5
+  falb_dir(:,is_sic,:) = 0.6
+  fevap(:,:) = 0.
+  DO i=1,nbsrf; qsolsrf(:,i)=150.; END DO
+  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
+  rain_fall  = 0.
+  snow_fall  = 0.
+  solsw      = 165.
+  sollw      = -53.
+  t_ancien   = 273.15
+  q_ancien   = 0.
+  agesno     = 0.
+
+  z0m(:,is_oce) = rugmer(:)
+  z0m(:,is_ter) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+  z0m(:,is_lic) = MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+  z0m(:,is_sic) = 0.001
+  z0h(:,:)=z0m(:,:)
+
+  fder    = 0.0
+  clwcon  = 0.0
+  rnebcon = 0.0
+  ratqs   = 0.0
+  run_off_lic_0 = 0.0 
+  rugoro  = 0.0
+
+! Before phyredem calling, surface modules and values to be saved in startphy.nc
+! are initialized
+!*******************************************************************************
+  dummy            = 1.0
+  pbl_tke(:,:,:)   = 1.e-8 
+  zmax0(:)         = 40.
+  f0(:)            = 1.e-5
+  sig1(:,:)        = 0.
+  w01(:,:)         = 0.
+  wake_deltat(:,:) = 0.
+  wake_deltaq(:,:) = 0.
+  wake_s(:)        = 0.
+  wake_cstar(:)    = 0.
+  wake_fip(:)      = 0.
+  wake_pe          = 0.
+  fm_therm         = 0.
+  entr_therm       = 0.
+  detr_therm       = 0.
+
+  CALL fonte_neige_init(run_off_lic_0)
+  CALL pbl_surface_init( fder, snsrf, qsolsrf, tsoil )
+  CALL phyredem( "startphy.nc" )
+
+!  WRITE(lunout,*)'CCCCCCCCCCCCCCCCCC REACTIVER SORTIE VISU DANS ETAT0'
+!  WRITE(lunout,*)'entree histclo'
+  CALL histclo()
+
+END SUBROUTINE etat0phys_netcdf
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_orog(lon_in,lat_in,phis,masque)
+!
+!===============================================================================
+! Comment:
+!   This routine launch grid_noro, which computes parameters for SSO scheme as
+!   described in LOTT & MILLER (1997) and LOTT(1999).
+!===============================================================================
+  USE conf_dat_m,  ONLY: conf_dat2d
+!  USE grid_atob_m, ONLY: rugsoro
+  USE grid_noro_m, ONLY: grid_noro
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! 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_orog"
+  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(:,:), tmp_var  (:,:)
+  REAL, ALLOCATABLE  :: zmea0(:,:), zstd0(:,:), zsig0(:,:)
+  REAL, ALLOCATABLE  :: zgam0(:,:), zthe0(:,:), zpic0(:,:), zval0(:,:)
+!-------------------------------------------------------------------------------
+  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")
+
+!--- 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 THE REQUIRED FIELDS USING ROUTINE grid_noro
+  WRITE(lunout,*)
+  WRITE(lunout,*)'*** Compute parameters needed for gravity wave drag code ***'
+
+!--- ALLOCATIONS OF SUB-CELL SCALES QUANTITIES
+  ALLOCATE(zmea0(iml,jml),zstd0(iml,jml)) !--- Mean orography and std deviation
+  ALLOCATE(zsig0(iml,jml),zgam0(iml,jml)) !--- Slope and nisotropy
+  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
+  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
+
+!--- CALL OROGRAPHY MODULE TO COMPUTE FIELDS
+  CALL grid_noro(lon_rad,lat_rad,relief_hi,lon_in,lat_in,phis,zmea0,zstd0,     &
+                                      zsig0,zgam0,zthe0,zpic0,zval0,masque)
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+
+!--- COMPUTE SURFACE ROUGHNESS
+!  WRITE(lunout,*)
+!  WRITE(lunout,*)'*** Compute surface roughness induced by the orography ***'
+!  ALLOCATE(tmp_var(iml-1,jml))
+!  CALL rugsoro(lon_rad, lat_rad, relief_hi, lon_in(1:iml-1), lat_in, tmp_var)
+!  ALLOCATE(rugo(iml,jml)); rugo(1:iml-1,:)=tmp_var; rugo(iml,:)=tmp_var(1,:)
+!  DEALLOCATE(tmp_var)
+  DEALLOCATE(relief_hi,lon_rad,lat_rad)
+
+!--- PUT QUANTITIES TO PHYSICAL GRID
+  CALL gr_dyn_fi(1,iml,jml,klon,zmea0,zmea); DEALLOCATE(zmea0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zstd0,zstd); DEALLOCATE(zstd0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zsig0,zsig); DEALLOCATE(zsig0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zgam0,zgam); DEALLOCATE(zgam0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zthe0,zthe); DEALLOCATE(zthe0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zpic0,zpic); DEALLOCATE(zpic0)
+  CALL gr_dyn_fi(1,iml,jml,klon,zval0,zval); DEALLOCATE(zval0)
+
+
+END SUBROUTINE start_init_orog
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_phys(lon_in,lat_in,lon_in2,lat_in2,ibar,phis)
+!
+!===============================================================================
+! Purpose:   Compute tsol and qsol, knowing phis.
+!===============================================================================
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! 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) :: phis(:,:)                   ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname="start_init_phys", physfname="ECPHY.nc"
+  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  :: ts(:,:), qs(:,:)
+!-------------------------------------------------------------------------------
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),SIZE(lon_in2),TRIM(modname)//" iml")
+  jml=assert_eq(SIZE(lat_in),SIZE(phis,2),              TRIM(modname)//" jml")
+  jml2=SIZE(lat_in2)
+
+  WRITE(lunout,*)'Opening the surface analysis'
+  CALL flininfo(physfname, iml_phys, jml_phys, llm_phys, ttm_phys, fid_phys)
+  WRITE(lunout,*) 'Values read: ',   iml_phys, jml_phys, llm_phys, ttm_phys
+
+  ALLOCATE(lat_phys(iml_phys,jml_phys),lon_phys(iml_phys,jml_phys))
+  ALLOCATE(levphys_ini(llm_phys))
+  CALL flinopen(physfname, .FALSE., iml_phys, jml_phys, llm_phys,              &
+                lon_phys,lat_phys,levphys_ini,ttm_phys,itau,date,dt,fid_phys)
+
+!--- 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))
+  CALL get_var_phys('ST'  ,ts)                   !--- SURFACE TEMPERATURE
+  CALL get_var_phys('CDSW',qs)                   !--- SOIL MOISTURE
+  CALL flinclo(fid_phys)
+  DEALLOCATE(var_ana,lon_rad,lat_rad,lon_ini,lat_ini)
+
+!--- TSOL AND QSOL ON PHYSICAL GRID
+  ALLOCATE(tsol(klon))
+  CALL gr_dyn_fi(1,iml,jml,klon,ts,tsol)
+  CALL gr_dyn_fi(1,iml,jml,klon,qs,qsol)
+  DEALLOCATE(ts,qs)
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE get_var_phys(title,field)
+!
+!-------------------------------------------------------------------------------
+  USE conf_dat_m, ONLY: conf_dat2d
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)    :: title
+  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: tllm
+!-------------------------------------------------------------------------------
+  SELECT CASE(title)
+    CASE('SP');        tllm=0
+    CASE('ST','CDSW'); tllm=llm_phys
+  END SELECT
+  IF(ALLOCATED(field)) RETURN
+  ALLOCATE(field(iml,jml)); field(:,:)=0.
+  CALL flinget(fid_phys,title,iml_phys,jml_phys,tllm,ttm_phys,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)
+
+END SUBROUTINE get_var_phys
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE start_init_phys
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+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
+!
+!-------------------------------------------------------------------------------
+
+
+END MODULE etat0phys
+!
+!*******************************************************************************
+
