Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/callphysiq_mod.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/callphysiq_mod.F90	(revision 5267)
+++ 	(revision )
@@ -1,109 +1,0 @@
-!
-! $Id: $
-!
-MODULE callphysiq_mod
-
-IMPLICIT NONE
-
-CONTAINS
-
-SUBROUTINE call_physiq(klon,llm,nqtot,tname,                              &
-                       debut_split,lafin_split,                           &
-                       jD_cur,jH_cur_split,zdt_split,                     &
-                       zplev_omp,zplay_omp,                               &
-                       zpk_omp,zphi_omp,zphis_omp,                        &
-                       presnivs_omp,                                      &
-                       zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,      &
-                       flxwfi_omp,pducov,                                 &
-                       zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,zdpsrf_omp)
-
-
-  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
-  USE control_mod, ONLY: planet_type, ok_dyn_xios
-  USE physiq_mod, ONLY: physiq
-  USE lmdz_xios, ONLY : xios_set_current_context, xios_get_current_context, xios_context
-  IMPLICIT NONE
-
-  INTEGER,INTENT(IN) :: klon ! (local) number of atmospheric columns
-  INTEGER,INTENT(IN) :: llm  ! number of atmospheric layers
-  INTEGER,INTENT(IN) :: nqtot ! number of tracers
-  CHARACTER(len=*),INTENT(IN) :: tname(nqtot) ! tracer names
-  LOGICAL,INTENT(IN) :: debut_split ! .true. if very first call to physics
-  LOGICAL,INTENT(IN) :: lafin_split ! .true. if last call to physics
-  REAL,INTENT(IN) :: JD_cur ! Julian day
-  REAL,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day)
-  REAL,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated
-  REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
-  REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
-  REAL,INTENT(IN) :: zpk_omp(klon,llm) ! Exner function
-  REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
-  REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential
-  REAL,INTENT(IN) :: presnivs_omp(llm) ! approximate pressure of atm. layers
-  REAL,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s)
-  REAL,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s)
-  REAL,INTENT(IN) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1
-  REAL,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K)
-  REAL,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air)
-  REAL,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) 
-  REAL,INTENT(IN) :: pducov(nbp_lon+1,nbp_lat,llm) ! dynamical tendency on ucov
-  ! tendencies (in */s) from the physics:
-  REAL,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds
-  REAL,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds
-  REAL,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature
-  REAL,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers
-  REAL,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure
-  TYPE(xios_context) :: dyn3d_ctx_handle
-
-  ! Local variables
-  CHARACTER(len=11) :: modname="call_physiq"
-  LOGICAL,SAVE :: firstcall=.true.
-!$OMP THREADPRIVATE(firstcall)
-
-! Sanity check on physics package type
-  IF (firstcall) THEN
-    IF (planet_type.ne."earth") THEN
-      CALL abort_gcm(modname,"wrong planet_type for this physics package",1)
-    ENDIF
-    firstcall=.false.
-  ENDIF
-
-  !$OMP MASTER
-  if (ok_dyn_xios) then
-     CALL xios_get_current_context(dyn3d_ctx_handle)
-  endif
-  !$OMP END MASTER
-
-! Call physics package with required inputs/outputs
-  CALL physiq(klon,           &
-              llm,            &
-              debut_split,    &
-              lafin_split,    &
-              zdt_split,      &
-              zplev_omp,      &
-              zplay_omp,      &
-              zphi_omp,       &
-              zphis_omp,      &
-              presnivs_omp,   &
-              zufi_omp,       &
-              zvfi_omp,       &
-              zrfi_omp,       &
-              ztfi_omp,       &
-              zqfi_omp,       &
-              flxwfi_omp,     &
-              zdufi_omp,      &
-              zdvfi_omp,      &
-              zdtfi_omp,      &
-              zdqfi_omp,      &
-              zdpsrf_omp)
-
-! switching back to LMDZDYN context
-!$OMP MASTER
-  if (ok_dyn_xios) then
-     CALL xios_set_current_context(dyn3d_ctx_handle)
-  endif
-!$OMP END MASTER
-
-
-END SUBROUTINE call_physiq
-
-END MODULE callphysiq_mod
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/callphysiq_mod.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/callphysiq_mod.f90	(revision 5268)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/callphysiq_mod.f90	(revision 5268)
@@ -0,0 +1,109 @@
+!
+! $Id: $
+!
+MODULE callphysiq_mod
+
+IMPLICIT NONE
+
+CONTAINS
+
+SUBROUTINE call_physiq(klon,llm,nqtot,tname,                              &
+                       debut_split,lafin_split,                           &
+                       jD_cur,jH_cur_split,zdt_split,                     &
+                       zplev_omp,zplay_omp,                               &
+                       zpk_omp,zphi_omp,zphis_omp,                        &
+                       presnivs_omp,                                      &
+                       zufi_omp,zvfi_omp,zrfi_omp,ztfi_omp,zqfi_omp,      &
+                       flxwfi_omp,pducov,                                 &
+                       zdufi_omp,zdvfi_omp,zdtfi_omp,zdqfi_omp,zdpsrf_omp)
+
+
+  USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat
+  USE control_mod, ONLY: planet_type, ok_dyn_xios
+  USE physiq_mod, ONLY: physiq
+  USE lmdz_xios, ONLY : xios_set_current_context, xios_get_current_context, xios_context
+  IMPLICIT NONE
+
+  INTEGER,INTENT(IN) :: klon ! (local) number of atmospheric columns
+  INTEGER,INTENT(IN) :: llm  ! number of atmospheric layers
+  INTEGER,INTENT(IN) :: nqtot ! number of tracers
+  CHARACTER(len=*),INTENT(IN) :: tname(nqtot) ! tracer names
+  LOGICAL,INTENT(IN) :: debut_split ! .true. if very first call to physics
+  LOGICAL,INTENT(IN) :: lafin_split ! .true. if last call to physics
+  REAL,INTENT(IN) :: JD_cur ! Julian day
+  REAL,INTENT(IN) :: JH_cur_split ! Julian hour (fraction of day)
+  REAL,INTENT(IN) :: zdt_split ! time step over which the physics are evaluated
+  REAL,INTENT(IN) :: zplev_omp(klon,llm+1) ! interlayer pressure (Pa)
+  REAL,INTENT(IN) :: zplay_omp(klon,llm) ! mid-layer pressure (Pa)
+  REAL,INTENT(IN) :: zpk_omp(klon,llm) ! Exner function
+  REAL,INTENT(IN) :: zphi_omp(klon,llm) ! geopotential at midlayer
+  REAL,INTENT(IN) :: zphis_omp(klon) ! surface geopotential
+  REAL,INTENT(IN) :: presnivs_omp(llm) ! approximate pressure of atm. layers
+  REAL,INTENT(IN) :: zufi_omp(klon,llm) ! zonal wind (m/s)
+  REAL,INTENT(IN) :: zvfi_omp(klon,llm) ! meridional wind (m/s)
+  REAL,INTENT(IN) :: zrfi_omp(klon,llm) ! relative wind vorticity, in s-1
+  REAL,INTENT(IN) :: ztfi_omp(klon,llm) ! temperature (K)
+  REAL,INTENT(IN) :: zqfi_omp(klon,llm,nqtot) ! tracers (*/kg of air)
+  REAL,INTENT(IN) :: flxwfi_omp(klon,llm) ! Vertical mass flux on lower mesh interfaces (kg/s) 
+  REAL,INTENT(IN) :: pducov(nbp_lon+1,nbp_lat,llm) ! dynamical tendency on ucov
+  ! tendencies (in */s) from the physics:
+  REAL,INTENT(OUT) :: zdufi_omp(klon,llm) ! tendency on zonal winds
+  REAL,INTENT(OUT) :: zdvfi_omp(klon,llm) ! tendency on meridional winds
+  REAL,INTENT(OUT) :: zdtfi_omp(klon,llm) ! tendency on temperature
+  REAL,INTENT(OUT) :: zdqfi_omp(klon,llm,nqtot) ! tendency on tracers
+  REAL,INTENT(OUT) :: zdpsrf_omp(klon) ! tendency on surface pressure
+  TYPE(xios_context) :: dyn3d_ctx_handle
+
+  ! Local variables
+  CHARACTER(len=11) :: modname="call_physiq"
+  LOGICAL,SAVE :: firstcall=.true.
+!$OMP THREADPRIVATE(firstcall)
+
+! Sanity check on physics package type
+  IF (firstcall) THEN
+    IF (planet_type.ne."earth") THEN
+      CALL abort_gcm(modname,"wrong planet_type for this physics package",1)
+    ENDIF
+    firstcall=.false.
+  ENDIF
+
+  !$OMP MASTER
+  if (ok_dyn_xios) then
+     CALL xios_get_current_context(dyn3d_ctx_handle)
+  endif
+  !$OMP END MASTER
+
+! Call physics package with required inputs/outputs
+  CALL physiq(klon,           &
+              llm,            &
+              debut_split,    &
+              lafin_split,    &
+              zdt_split,      &
+              zplev_omp,      &
+              zplay_omp,      &
+              zphi_omp,       &
+              zphis_omp,      &
+              presnivs_omp,   &
+              zufi_omp,       &
+              zvfi_omp,       &
+              zrfi_omp,       &
+              ztfi_omp,       &
+              zqfi_omp,       &
+              flxwfi_omp,     &
+              zdufi_omp,      &
+              zdvfi_omp,      &
+              zdtfi_omp,      &
+              zdqfi_omp,      &
+              zdpsrf_omp)
+
+! switching back to LMDZDYN context
+!$OMP MASTER
+  if (ok_dyn_xios) then
+     CALL xios_set_current_context(dyn3d_ctx_handle)
+  endif
+!$OMP END MASTER
+
+
+END SUBROUTINE call_physiq
+
+END MODULE callphysiq_mod
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.F90	(revision 5267)
+++ 	(revision )
@@ -1,576 +1,0 @@
-!
-! $Id$
-!
-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,             ONLY: klon
-  USE conf_dat_m,         ONLY: conf_dat2d
-  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
-          solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
-          sollw,sollwdown, 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, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke,  &
-    phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
-    prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
-    ale_bl, ale_bl_trig, alp_bl, &
-    ale_wake, ale_bl_stat, AWAKE_S, &
-    cf_ancien, rvc_ancien
-
-  USE comconst_mod, ONLY: pi, dtvr
-
-  PRIVATE
-  PUBLIC :: etat0phys_netcdf
-
-  include "iniprint.h"
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-  include "dimsoil.h"
-  include "clesphys.h"
-  REAL, SAVE :: deg2rad
-  REAL, SAVE, ALLOCATABLE :: tsol(:)
-  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 :: oroparam="oro_params.nc"
-  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
-  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
-  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
-
-
-CONTAINS
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE etat0phys_netcdf(masque, phis)
-!
-!-------------------------------------------------------------------------------
-! Purpose: Creates initial states
-!-------------------------------------------------------------------------------
-! Notes:  1) This routine is designed to work for Earth
-!         2) If masque(:,:)/=-99999., masque and phis are already known.
-!         Otherwise: compute it.
-!-------------------------------------------------------------------------------
-  USE control_mod
-  USE fonte_neige_mod
-  USE pbl_surface_mod
-  USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
-  USE indice_sol_mod
-  USE conf_phys_m, ONLY: conf_phys
-  USE init_ssrf_m, ONLY: start_init_subsurf
-  USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
-       ratqs_inter_
-  !use ioipsl
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  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
-  LOGICAL            :: read_mask
-  REAL               :: phystep, dummy
-  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
-  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
-  REAL, DIMENSION(klon,nbsrf)         :: qsurf, snsrf
-  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
-
-!--- 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_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
-  INTEGER :: flag_aerosol
-  INTEGER :: flag_aerosol_strat
-  INTEGER :: flag_volc_surfstrat
-  LOGICAL :: flag_aer_feedback
-  LOGICAL :: flag_bc_internal_mixture
-  REAL    :: bl95_b0, bl95_b1
-  INTEGER :: read_climoz                        !--- Read ozone climatology
-  REAL    :: alp_offset
-  LOGICAL :: filtre_oro=.false.
-
-  INCLUDE "compbl.h"
-  INCLUDE "alpale.h"
-  
-  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")
-
-! Physics configuration
-!*******************************************************************************
-  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,                                        &
-                   ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat,     &
-                   aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat,  & 
-                   flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1,       &
-                   read_climoz, alp_offset)
-  CALL phys_state_var_init(read_climoz)
-
-!--- Initial atmospheric CO2 conc. from .def file
-  co2_ppm0 = co2_ppm
-
-! 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)
-
-  CALL getin('filtre_oro',filtre_oro)
-  IF (filtre_oro) CALL filtreoro(size(phis,1),size(phis,2),phis,masque_tmp,rlatu)
-
-  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(rlonu, rlatv, phis)
-
-! Some initializations.
-!*******************************************************************************
-  sn    (:) = 0.0                               !--- Snow
-  radsol(:) = 0.0                               !--- Net radiation at ground
-  rugmer(:) = 0.001                             !--- Ocean rugosity
-  !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
-  IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
-
-! Sub-surfaces initialization.
-!*******************************************************************************
-  CALL start_init_subsurf(read_mask)
-
-! 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, qsurf, 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
-
-!ym warning missing init for falb_dif => set to 0
-  falb_dif(:,:,:)=0
-
-  u10m(:,:)=0  
-  v10m(:,:)=0  
-  treedrg(:,:,:)=0
-
-  fevap(:,:) = 0.
-  qsurf = 0.
-  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
-  rain_fall  = 0.
-  snow_fall  = 0.
-  solsw      = 165.
-  solswfdiff = 1. 
-  sollw      = -53.
-!ym warning missing init for sollwdown => set to 0
-  sollwdown  = 0.
-  t_ancien   = 273.15
-  q_ancien   = 0.
-  ql_ancien = 0.
-  qs_ancien = 0.
-  prlw_ancien = 0.
-  prsw_ancien = 0.
-  prw_ancien = 0.
-  agesno     = 0.
-  
-  u_ancien = 0.
-  v_ancien = 0.
-  wake_delta_pbl_TKE(:,:,:)=0
-  wake_dens(:)=0
-  awake_dens = 0.
-  cv_gen = 0.
-  ale_bl = 0.
-  ale_bl_trig =0.
-  alp_bl=0.
-  ale_wake=0.
-  ale_bl_stat=0.
-
-  cf_ancien = 0.
-  rvc_ancien = 0.
-  
-  z0m(:,:)=0 ! ym missing 5th subsurface initialization
- 
-  z0m(:,is_oce) = rugmer(:)
-  z0m(:,is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
-  z0m(:,is_lic) = 0.001 !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.
-  awake_s = 0.
-
-  CALL fonte_neige_init(run_off_lic_0)
-  CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
-
-  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then
-     delta_tsurf = 0.
-     beta_aridity = 0.
-  end IF
-
-  ratqs_inter_ = 0.002
-  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).
-!   In case the file oroparam is present and the key read_orop is activated,
-!   grid_noro is bypassed and sub-cell parameters are read from the file.
-!===============================================================================
-  USE grid_noro_m, ONLY: grid_noro, read_noro
-  USE logic_mod,   ONLY: read_orop
-  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
-  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
-  INTEGER            :: ierr
-  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(:,:)
-!-------------------------------------------------------------------------------
-  modname="start_init_orog"
-  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, orogvar, 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(orogvar, 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
-  zsig0(:,:)=0   !ym uninitialized variable
-  zgam0(:,:)=0   !ym uninitialized variable
-  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
-  zthe0(:,:)=0   !ym uninitialized variable
-  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
-
-!--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
-  OPEN(UNIT=66,FILE=oroparam,STATUS='OLD',IOSTAT=ierr)
-  IF(ierr==0.AND.read_orop) THEN
-    CLOSE(UNIT=66)
-    CALL read_noro(lon_in,lat_in,oroparam,                                     &
-                   phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
-  ELSE
-!--- 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)
-  END IF
-  phis = phis * 9.81
-  phis(iml,:) = phis(1,:)
-  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,phis)
-!
-!===============================================================================
-! Purpose:   Compute tsol and qsol, knowing phis.
-!===============================================================================
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
-  REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname
-  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(:,:)
-!-------------------------------------------------------------------------------
-  modname="start_init_phys"
-  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
-  jml=SIZE(phis,2); jml2=SIZE(lat_in)
-
-  WRITE(lunout,*)'Opening the surface analysis'
-  CALL flininfo(phyfname, 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(phyfname, .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(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
-  CALL get_var_phys(qsolvar,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)
-!
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*),  INTENT(IN)    :: title
-  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
-!-------------------------------------------------------------------------------
-! Local variables:
-  INTEGER :: tllm
-!-------------------------------------------------------------------------------
-  SELECT CASE(title)
-    CASE(psrfvar);         tllm=0
-    CASE(tsrfvar,qsolvar); 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, .TRUE.)
-  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
-                                      lon_in,  lat_in,  field)
-
-END SUBROUTINE get_var_phys
-!
-!-------------------------------------------------------------------------------
-!
-END SUBROUTINE start_init_phys
-!
-!-------------------------------------------------------------------------------
-
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
-!
-!-------------------------------------------------------------------------------
-  USE inter_barxy_m, ONLY: inter_barxy
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  CHARACTER(LEN=*), INTENT(IN)  :: nam
-  LOGICAL,          INTENT(IN)  :: ibeg
-  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
-  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
-  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
-  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
-!-------------------------------------------------------------------------------
-! Local variables:
-  CHARACTER(LEN=256) :: modname
-  INTEGER            :: ii, jj, i1, j1, j2
-  REAL, ALLOCATABLE  :: vtmp(:,:)
-!-------------------------------------------------------------------------------
-  modname="interp_startvar"
-  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(lon2),SIZE(varo,1),TRIM(modname)//" i1")
-  j1=SIZE(varo,2); j2=SIZE(lat2)
-  ALLOCATE(vtmp(i1-1,j1))
-  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)
-  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
-
-END SUBROUTINE interp_startvar
-!
-!-------------------------------------------------------------------------------
-!
-!*******************************************************************************
-
-SUBROUTINE filtreoro(imp1,jmp1,phis,masque,rlatu)
-
-IMPLICIT NONE
-
-  INTEGER imp1,jmp1
-  REAL, DIMENSION(imp1,jmp1) :: phis,masque
-  REAL, DIMENSION(jmp1) :: rlatu
-  REAL, DIMENSION(imp1) :: wwf
-  REAL, DIMENSION(imp1,jmp1) :: phiso
-  INTEGER :: ifiltre,ifi,ii,i,j
-  REAL :: coslat0,ssz
-
-  coslat0=0.5
-  phiso=phis
-  do j=2,jmp1-1
-     print*,'avant if ',cos(rlatu(j)),coslat0
-     if (cos(rlatu(j))<coslat0) then
-         ! nb de pts affectes par le filtrage de part et d'autre du pt
-         ifiltre=(coslat0/cos(rlatu(j))-1.)/2.
-         wwf=0.
-         do i=1,ifiltre
-            wwf(i)=1.
-         enddo
-         wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre
-         do i=1,imp1-1
-            if (masque(i,j)>0.9) then
-               ssz=phis(i,j)
-               do ifi=1,ifiltre+1
-                  ii=i+ifi
-                  if (ii>imp1-1) ii=ii-imp1+1
-                  ssz=ssz+wwf(ifi)*phis(ii,j)
-                  ii=i-ifi
-                  if (ii<1) ii=ii+imp1-1
-                  ssz=ssz+wwf(ifi)*phis(ii,j)
-               enddo
-               phis(i,j)=ssz*cos(rlatu(j))/coslat0
-            endif
-         enddo
-         print*,'j=',j,coslat0/cos(rlatu(j)), (1.+2.*sum(wwf))*cos(rlatu(j))/coslat0
-     endif
-  enddo
-  call dump2d(imp1,jmp1,phis,'phis ')
-  call dump2d(imp1,jmp1,masque,'masque ')
-  call dump2d(imp1,jmp1,phis-phiso,'dphis ')
-
-END SUBROUTINE filtreoro
-
-
-END MODULE etat0phys
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90	(revision 5268)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/etat0phys_netcdf.f90	(revision 5268)
@@ -0,0 +1,576 @@
+!
+! $Id$
+!
+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,             ONLY: klon
+  USE conf_dat_m,         ONLY: conf_dat2d
+  USE phys_state_var_mod, ONLY: zmea, zstd, zsig, zgam, zthe, zpic, zval, z0m, &
+          solsw, solswfdiff, radsol, t_ancien, wake_deltat, wake_s,  rain_fall, qsol, z0h, &
+          sollw,sollwdown, 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, falb_dif, wake_fip,    agesno,  detr_therm, pbl_tke,  &
+    phys_state_var_init, ql_ancien, qs_ancien, prlw_ancien, prsw_ancien, &
+    prw_ancien, u10m,v10m, treedrg, u_ancien, v_ancien, wake_delta_pbl_TKE, wake_dens, &
+    ale_bl, ale_bl_trig, alp_bl, &
+    ale_wake, ale_bl_stat, AWAKE_S, &
+    cf_ancien, rvc_ancien
+
+  USE comconst_mod, ONLY: pi, dtvr
+
+  PRIVATE
+  PUBLIC :: etat0phys_netcdf
+
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+  include "dimsoil.h"
+  include "clesphys.h"
+  REAL, SAVE :: deg2rad
+  REAL, SAVE, ALLOCATABLE :: tsol(:)
+  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 :: oroparam="oro_params.nc"
+  CHARACTER(LEN=256), PARAMETER :: orofname="Relief.nc", orogvar="RELIEF"
+  CHARACTER(LEN=256), PARAMETER :: phyfname="ECPHY.nc",  psrfvar="SP"
+  CHARACTER(LEN=256), PARAMETER :: qsolvar="CDSW",       tsrfvar="ST"
+
+
+CONTAINS
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE etat0phys_netcdf(masque, phis)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Creates initial states
+!-------------------------------------------------------------------------------
+! Notes:  1) This routine is designed to work for Earth
+!         2) If masque(:,:)/=-99999., masque and phis are already known.
+!         Otherwise: compute it.
+!-------------------------------------------------------------------------------
+  USE control_mod
+  USE fonte_neige_mod
+  USE pbl_surface_mod
+  USE regr_horiz_time_climoz_m, ONLY: regr_horiz_time_climoz
+  USE indice_sol_mod
+  USE conf_phys_m, ONLY: conf_phys
+  USE init_ssrf_m, ONLY: start_init_subsurf
+  USE phys_state_var_mod, ONLY: beta_aridity, delta_tsurf, awake_dens, cv_gen, &
+       ratqs_inter_
+  !use ioipsl
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  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
+  LOGICAL            :: read_mask
+  REAL               :: phystep, dummy
+  REAL, DIMENSION(SIZE(masque,1),SIZE(masque,2)) :: masque_tmp,phiso
+  REAL, DIMENSION(klon)               :: sn, rugmer, run_off_lic_0, fder
+  REAL, DIMENSION(klon,nbsrf)         :: qsurf, snsrf
+  REAL, DIMENSION(klon,nsoilmx,nbsrf) :: tsoil
+
+!--- 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_volcan, ok_alw, ok_cdnc, aerosol_couple, chemistry_couple
+  INTEGER :: flag_aerosol
+  INTEGER :: flag_aerosol_strat
+  INTEGER :: flag_volc_surfstrat
+  LOGICAL :: flag_aer_feedback
+  LOGICAL :: flag_bc_internal_mixture
+  REAL    :: bl95_b0, bl95_b1
+  INTEGER :: read_climoz                        !--- Read ozone climatology
+  REAL    :: alp_offset
+  LOGICAL :: filtre_oro=.false.
+
+  INCLUDE "compbl.h"
+  INCLUDE "alpale.h"
+  
+  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")
+
+! Physics configuration
+!*******************************************************************************
+  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,                                        &
+                   ratqsbas,ratqshaut,tau_ratqs,            &
+                   ok_ade, ok_aie, ok_alw, ok_cdnc, ok_volcan, flag_volc_surfstrat,     &
+                   aerosol_couple, chemistry_couple, flag_aerosol, flag_aerosol_strat,  & 
+                   flag_aer_feedback, flag_bc_internal_mixture, bl95_b0, bl95_b1,       &
+                   read_climoz, alp_offset)
+  CALL phys_state_var_init(read_climoz)
+
+!--- Initial atmospheric CO2 conc. from .def file
+  co2_ppm0 = co2_ppm
+
+! 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)
+
+  CALL getin('filtre_oro',filtre_oro)
+  IF (filtre_oro) CALL filtreoro(size(phis,1),size(phis,2),phis,masque_tmp,rlatu)
+
+  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(rlonu, rlatv, phis)
+
+! Some initializations.
+!*******************************************************************************
+  sn    (:) = 0.0                               !--- Snow
+  radsol(:) = 0.0                               !--- Net radiation at ground
+  rugmer(:) = 0.001                             !--- Ocean rugosity
+  !--- Ozone (zonal or 3D) interpolation in space and time (if 2nd arg is TRUE)
+  IF(read_climoz>=1) CALL regr_horiz_time_climoz(read_climoz,ok_daily_climoz)
+
+! Sub-surfaces initialization.
+!*******************************************************************************
+  CALL start_init_subsurf(read_mask)
+
+! 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, qsurf, 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
+
+!ym warning missing init for falb_dif => set to 0
+  falb_dif(:,:,:)=0
+
+  u10m(:,:)=0  
+  v10m(:,:)=0  
+  treedrg(:,:,:)=0
+
+  fevap(:,:) = 0.
+  qsurf = 0.
+  DO i=1,nbsrf; DO j=1,nsoilmx; tsoil(:,j,i) = tsol; END DO; END DO
+  rain_fall  = 0.
+  snow_fall  = 0.
+  solsw      = 165.
+  solswfdiff = 1. 
+  sollw      = -53.
+!ym warning missing init for sollwdown => set to 0
+  sollwdown  = 0.
+  t_ancien   = 273.15
+  q_ancien   = 0.
+  ql_ancien = 0.
+  qs_ancien = 0.
+  prlw_ancien = 0.
+  prsw_ancien = 0.
+  prw_ancien = 0.
+  agesno     = 0.
+  
+  u_ancien = 0.
+  v_ancien = 0.
+  wake_delta_pbl_TKE(:,:,:)=0
+  wake_dens(:)=0
+  awake_dens = 0.
+  cv_gen = 0.
+  ale_bl = 0.
+  ale_bl_trig =0.
+  alp_bl=0.
+  ale_wake=0.
+  ale_bl_stat=0.
+
+  cf_ancien = 0.
+  rvc_ancien = 0.
+  
+  z0m(:,:)=0 ! ym missing 5th subsurface initialization
+ 
+  z0m(:,is_oce) = rugmer(:)
+  z0m(:,is_ter) = 0.01 !MAX(1.0e-05,zstd(:)*zsig(:)/2.0)
+  z0m(:,is_lic) = 0.001 !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.
+  awake_s = 0.
+
+  CALL fonte_neige_init(run_off_lic_0)
+  CALL pbl_surface_init( fder, snsrf, qsurf, tsoil )
+
+  IF (iflag_pbl>1 .AND. iflag_wake>=1  .AND. iflag_pbl_split >=1) then
+     delta_tsurf = 0.
+     beta_aridity = 0.
+  end IF
+
+  ratqs_inter_ = 0.002
+  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).
+!   In case the file oroparam is present and the key read_orop is activated,
+!   grid_noro is bypassed and sub-cell parameters are read from the file.
+!===============================================================================
+  USE grid_noro_m, ONLY: grid_noro, read_noro
+  USE logic_mod,   ONLY: read_orop
+  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
+  INTEGER            :: fid, llm_tmp,ttm_tmp, iml,jml, iml_rel,jml_rel, itau(1)
+  INTEGER            :: ierr
+  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(:,:)
+!-------------------------------------------------------------------------------
+  modname="start_init_orog"
+  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, orogvar, 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(orogvar, 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
+  zsig0(:,:)=0   !ym uninitialized variable
+  zgam0(:,:)=0   !ym uninitialized variable
+  ALLOCATE(zthe0(iml,jml))                !--- Highest slope orientation
+  zthe0(:,:)=0   !ym uninitialized variable
+  ALLOCATE(zpic0(iml,jml),zval0(iml,jml)) !--- Peaks and valley heights
+
+!--- READ SUB-CELL SCALES PARAMETERS FROM A FILE (AT RIGHT RESOLUTION)
+  OPEN(UNIT=66,FILE=oroparam,STATUS='OLD',IOSTAT=ierr)
+  IF(ierr==0.AND.read_orop) THEN
+    CLOSE(UNIT=66)
+    CALL read_noro(lon_in,lat_in,oroparam,                                     &
+                   phis,zmea0,zstd0,zsig0,zgam0,zthe0,zpic0,zval0,masque)
+  ELSE
+!--- 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)
+  END IF
+  phis = phis * 9.81
+  phis(iml,:) = phis(1,:)
+  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,phis)
+!
+!===============================================================================
+! Purpose:   Compute tsol and qsol, knowing phis.
+!===============================================================================
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  REAL,    INTENT(IN) :: lon_in(:),  lat_in(:)       ! dim (iml) (jml2)
+  REAL,    INTENT(IN) :: phis(:,:)                   ! dim (iml,jml)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname
+  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(:,:)
+!-------------------------------------------------------------------------------
+  modname="start_init_phys"
+  iml=assert_eq(SIZE(lon_in),SIZE(phis,1),TRIM(modname)//" iml")
+  jml=SIZE(phis,2); jml2=SIZE(lat_in)
+
+  WRITE(lunout,*)'Opening the surface analysis'
+  CALL flininfo(phyfname, 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(phyfname, .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(tsrfvar,ts)                   !--- SURFACE TEMPERATURE
+  CALL get_var_phys(qsolvar,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)
+!
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*),  INTENT(IN)    :: title
+  REAL, ALLOCATABLE, INTENT(INOUT) :: field(:,:)
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER :: tllm
+!-------------------------------------------------------------------------------
+  SELECT CASE(title)
+    CASE(psrfvar);         tllm=0
+    CASE(tsrfvar,qsolvar); 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, .TRUE.)
+  CALL interp_startvar(title, .TRUE., lon_rad, lat_rad, var_ana,               &
+                                      lon_in,  lat_in,  field)
+
+END SUBROUTINE get_var_phys
+!
+!-------------------------------------------------------------------------------
+!
+END SUBROUTINE start_init_phys
+!
+!-------------------------------------------------------------------------------
+
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE interp_startvar(nam,ibeg,lon,lat,vari,lon2,lat2,varo)
+!
+!-------------------------------------------------------------------------------
+  USE inter_barxy_m, ONLY: inter_barxy
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  CHARACTER(LEN=*), INTENT(IN)  :: nam
+  LOGICAL,          INTENT(IN)  :: ibeg
+  REAL,             INTENT(IN)  :: lon(:), lat(:)   ! dim (ii) (jj)
+  REAL,             INTENT(IN)  :: vari(:,:)        ! dim (ii,jj)
+  REAL,             INTENT(IN)  :: lon2(:), lat2(:) ! dim (i1) (j2)
+  REAL,             INTENT(OUT) :: varo(:,:)        ! dim (i1) (j1)
+!-------------------------------------------------------------------------------
+! Local variables:
+  CHARACTER(LEN=256) :: modname
+  INTEGER            :: ii, jj, i1, j1, j2
+  REAL, ALLOCATABLE  :: vtmp(:,:)
+!-------------------------------------------------------------------------------
+  modname="interp_startvar"
+  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(lon2),SIZE(varo,1),TRIM(modname)//" i1")
+  j1=SIZE(varo,2); j2=SIZE(lat2)
+  ALLOCATE(vtmp(i1-1,j1))
+  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)
+  CALL gr_int_dyn(vtmp, varo, i1-1, j1)
+
+END SUBROUTINE interp_startvar
+!
+!-------------------------------------------------------------------------------
+!
+!*******************************************************************************
+
+SUBROUTINE filtreoro(imp1,jmp1,phis,masque,rlatu)
+
+IMPLICIT NONE
+
+  INTEGER imp1,jmp1
+  REAL, DIMENSION(imp1,jmp1) :: phis,masque
+  REAL, DIMENSION(jmp1) :: rlatu
+  REAL, DIMENSION(imp1) :: wwf
+  REAL, DIMENSION(imp1,jmp1) :: phiso
+  INTEGER :: ifiltre,ifi,ii,i,j
+  REAL :: coslat0,ssz
+
+  coslat0=0.5
+  phiso=phis
+  do j=2,jmp1-1
+     print*,'avant if ',cos(rlatu(j)),coslat0
+     if (cos(rlatu(j))<coslat0) then
+         ! nb de pts affectes par le filtrage de part et d'autre du pt
+         ifiltre=(coslat0/cos(rlatu(j))-1.)/2.
+         wwf=0.
+         do i=1,ifiltre
+            wwf(i)=1.
+         enddo
+         wwf(ifiltre+1)=(coslat0/cos(rlatu(j))-1.)/2.-ifiltre
+         do i=1,imp1-1
+            if (masque(i,j)>0.9) then
+               ssz=phis(i,j)
+               do ifi=1,ifiltre+1
+                  ii=i+ifi
+                  if (ii>imp1-1) ii=ii-imp1+1
+                  ssz=ssz+wwf(ifi)*phis(ii,j)
+                  ii=i-ifi
+                  if (ii<1) ii=ii+imp1-1
+                  ssz=ssz+wwf(ifi)*phis(ii,j)
+               enddo
+               phis(i,j)=ssz*cos(rlatu(j))/coslat0
+            endif
+         enddo
+         print*,'j=',j,coslat0/cos(rlatu(j)), (1.+2.*sum(wwf))*cos(rlatu(j))/coslat0
+     endif
+  enddo
+  call dump2d(imp1,jmp1,phis,'phis ')
+  call dump2d(imp1,jmp1,masque,'masque ')
+  call dump2d(imp1,jmp1,phis-phiso,'dphis ')
+
+END SUBROUTINE filtreoro
+
+
+END MODULE etat0phys
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/init_ssrf_m.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/init_ssrf_m.F90	(revision 5267)
+++ 	(revision )
@@ -1,148 +1,0 @@
-MODULE init_ssrf_m
-!
-!*******************************************************************************
-
-  USE indice_sol_mod,     ONLY: is_ter, is_oce, is_oce, is_lic, epsfra
-  USE dimphy,             ONLY: klon, zmasq
-  USE phys_state_var_mod, ONLY: pctsrf
-  USE geometry_mod,       ONLY : longitude_deg, latitude_deg
-  USE grid_atob_m,        ONLY: grille_m
-  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
-  USE ioipsl_getin_p_mod, ONLY: getin_p
-  USE comconst_mod,       ONLY: im, pi
-  USE surface_data,       ONLY: landice_opt
-
-  CHARACTER(LEN=256), PARAMETER :: icefname="landiceref.nc", icevar="landice"
-  PRIVATE
-  PUBLIC :: start_init_subsurf
-  include "iniprint.h"
-  include "dimensions.h"
-  include "paramet.h"
-  include "comgeom2.h"
-
-CONTAINS
-
-!-------------------------------------------------------------------------------
-!
-SUBROUTINE start_init_subsurf(known_mask)
-!
-!-------------------------------------------------------------------------------
-! Purpose: Subsurfaces initialization.
-!-------------------------------------------------------------------------------
-! Comment: Called by etat0phys_netcdf ; also called by limit_netcdf in case
-!          no starting states are required (ok_etat0==.FALSE.).
-!-------------------------------------------------------------------------------
-  IMPLICIT NONE
-!-------------------------------------------------------------------------------
-! Arguments:
-  LOGICAL, INTENT(IN) :: known_mask
-!-------------------------------------------------------------------------------
-! Local variables:
-  INTEGER           :: iml_lic, jml_lic
-  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1), ji, j
-  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic (:,:)
-  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:), flic_tmp(:,:), vtmp(:,:)
-  REAL              :: date, lev(1), dt, deg2rad
-  LOGICAL           :: no_ter_antartique   ! If true, no land points are allowed at Antartic
-!-------------------------------------------------------------------------------
-  deg2rad= pi/180.0
-
-!--- Physical grid points coordinates
-  DO j=2,jjm; latitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);    END DO
-  DO j=2,jjm; longitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:im); END DO
-  latitude_deg(1) = pi/2.; latitude_deg(klon) = - pi/2.
-  latitude_deg(:)=latitude_deg(:)/deg2rad
-  longitude_deg(1) = 0.0;   longitude_deg(klon) = 0.0;
-  longitude_deg(:)=longitude_deg(:)/deg2rad
-
-! Compute ground geopotential, sub-cells quantities and possibly the mask.
-! Sub-surfaces initialization
-!*******************************************************************************
-  IF (landice_opt .LT. 2) THEN
-     ! Continue with reading landice.nc file
-     WRITE(lunout,*)"Read landice.nc file to attribute is_lic fraction"
-
-     !--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
-     CALL flininfo(icefname, iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
-     ALLOCATE(lat_lic(iml_lic,jml_lic),lon_lic(iml_lic,jml_lic))
-     ALLOCATE(fraclic(iml_lic,jml_lic))
-     CALL flinopen(icefname, .FALSE., iml_lic, jml_lic, llm_tmp,  &
-          &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
-     CALL flinget(fid, icevar, 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
-     
-     ALLOCATE(dlon_lic(iml_lic),dlat_lic(jml_lic))
-     dlon_lic(:)=lon_lic(:,1); IF(MAXVAL(dlon_lic)>pi) dlon_lic=dlon_lic*pi/180.
-     dlat_lic(:)=lat_lic(1,:); IF(MAXVAL(dlat_lic)>pi) dlat_lic=dlat_lic*pi/180.
-     DEALLOCATE(lon_lic,lat_lic); ALLOCATE(flic_tmp(iip1,jjp1))
-     CALL grille_m(dlon_lic,dlat_lic,fraclic,rlonv(1:iim),rlatu,flic_tmp(1:iim,:))
-     flic_tmp(iip1,:)=flic_tmp(1,:)
-     
-     !--- To the physical grid
-     pctsrf(:,:) = 0.
-     CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
-     DEALLOCATE(flic_tmp)
-     
-     !--- 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
-     
-     
-     !--- Option no_ter_antartique removes all land fractions souther than 60S.
-     !--- Land ice is set instead of the land fractions on these latitudes.
-     !--- The ocean and sea-ice fractions are not changed.
-     no_ter_antartique=.FALSE.
-     CALL getin_p('no_ter_antartique',no_ter_antartique)
-     WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
-     IF (no_ter_antartique) THEN
-        ! Remove all land fractions souther than 60S and set land-ice instead
-        WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
-        WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
-        DO ji=1, klon
-           IF (latitude_deg(ji)<-60.0) THEN
-              pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
-              pctsrf(ji,is_ter) = 0
-           END IF
-        END DO
-     END IF
-     
-  ELSE
-     ! landice_opt=2 and higher
-     WRITE(lunout,*) 'No landice is attributed is_lic sub-surface because landice_opt=2 or higher.'
-     WRITE(lunout,*) 'This means that the land model will handel land ice as well as all other land areas.'
-     pctsrf(:,is_ter) = zmasq(:)
-     pctsrf(:,is_lic) = 0.0
-  END IF
-
-!--- 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(known_mask) pctsrf(:,is_oce)=1-zmasq(:)
-
-!--- 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'
-
-END SUBROUTINE start_init_subsurf
-!
-!-------------------------------------------------------------------------------
-
-END MODULE init_ssrf_m
-!
-!*******************************************************************************
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/init_ssrf_m.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/init_ssrf_m.f90	(revision 5268)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/init_ssrf_m.f90	(revision 5268)
@@ -0,0 +1,148 @@
+MODULE init_ssrf_m
+!
+!*******************************************************************************
+
+  USE indice_sol_mod,     ONLY: is_ter, is_oce, is_oce, is_lic, epsfra
+  USE dimphy,             ONLY: klon, zmasq
+  USE phys_state_var_mod, ONLY: pctsrf
+  USE geometry_mod,       ONLY : longitude_deg, latitude_deg
+  USE grid_atob_m,        ONLY: grille_m
+  USE ioipsl,             ONLY: flininfo, flinopen, flinget, flinclo
+  USE ioipsl_getin_p_mod, ONLY: getin_p
+  USE comconst_mod,       ONLY: im, pi
+  USE surface_data,       ONLY: landice_opt
+
+  CHARACTER(LEN=256), PARAMETER :: icefname="landiceref.nc", icevar="landice"
+  PRIVATE
+  PUBLIC :: start_init_subsurf
+  include "iniprint.h"
+  include "dimensions.h"
+  include "paramet.h"
+  include "comgeom2.h"
+
+CONTAINS
+
+!-------------------------------------------------------------------------------
+!
+SUBROUTINE start_init_subsurf(known_mask)
+!
+!-------------------------------------------------------------------------------
+! Purpose: Subsurfaces initialization.
+!-------------------------------------------------------------------------------
+! Comment: Called by etat0phys_netcdf ; also called by limit_netcdf in case
+!          no starting states are required (ok_etat0==.FALSE.).
+!-------------------------------------------------------------------------------
+  IMPLICIT NONE
+!-------------------------------------------------------------------------------
+! Arguments:
+  LOGICAL, INTENT(IN) :: known_mask
+!-------------------------------------------------------------------------------
+! Local variables:
+  INTEGER           :: iml_lic, jml_lic
+  INTEGER           :: fid, llm_tmp, ttm_tmp, itaul(1), ji, j
+  REAL, ALLOCATABLE :: dlon_lic(:), lon_lic(:,:), fraclic (:,:)
+  REAL, ALLOCATABLE :: dlat_lic(:), lat_lic(:,:), flic_tmp(:,:), vtmp(:,:)
+  REAL              :: date, lev(1), dt, deg2rad
+  LOGICAL           :: no_ter_antartique   ! If true, no land points are allowed at Antartic
+!-------------------------------------------------------------------------------
+  deg2rad= pi/180.0
+
+!--- Physical grid points coordinates
+  DO j=2,jjm; latitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlatu(j);    END DO
+  DO j=2,jjm; longitude_deg((j-2)*iim+2:(j-1)*iim+1)=rlonv(1:im); END DO
+  latitude_deg(1) = pi/2.; latitude_deg(klon) = - pi/2.
+  latitude_deg(:)=latitude_deg(:)/deg2rad
+  longitude_deg(1) = 0.0;   longitude_deg(klon) = 0.0;
+  longitude_deg(:)=longitude_deg(:)/deg2rad
+
+! Compute ground geopotential, sub-cells quantities and possibly the mask.
+! Sub-surfaces initialization
+!*******************************************************************************
+  IF (landice_opt .LT. 2) THEN
+     ! Continue with reading landice.nc file
+     WRITE(lunout,*)"Read landice.nc file to attribute is_lic fraction"
+
+     !--- Read and interpolate on model T-grid soil fraction and soil ice fraction.
+     CALL flininfo(icefname, iml_lic, jml_lic, llm_tmp, ttm_tmp, fid)
+     ALLOCATE(lat_lic(iml_lic,jml_lic),lon_lic(iml_lic,jml_lic))
+     ALLOCATE(fraclic(iml_lic,jml_lic))
+     CALL flinopen(icefname, .FALSE., iml_lic, jml_lic, llm_tmp,  &
+          &               lon_lic, lat_lic, lev, ttm_tmp, itaul, date, dt, fid)
+     CALL flinget(fid, icevar, 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
+     
+     ALLOCATE(dlon_lic(iml_lic),dlat_lic(jml_lic))
+     dlon_lic(:)=lon_lic(:,1); IF(MAXVAL(dlon_lic)>pi) dlon_lic=dlon_lic*pi/180.
+     dlat_lic(:)=lat_lic(1,:); IF(MAXVAL(dlat_lic)>pi) dlat_lic=dlat_lic*pi/180.
+     DEALLOCATE(lon_lic,lat_lic); ALLOCATE(flic_tmp(iip1,jjp1))
+     CALL grille_m(dlon_lic,dlat_lic,fraclic,rlonv(1:iim),rlatu,flic_tmp(1:iim,:))
+     flic_tmp(iip1,:)=flic_tmp(1,:)
+     
+     !--- To the physical grid
+     pctsrf(:,:) = 0.
+     CALL gr_dyn_fi(1, iip1, jjp1, klon, flic_tmp, pctsrf(:,is_lic))
+     DEALLOCATE(flic_tmp)
+     
+     !--- 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
+     
+     
+     !--- Option no_ter_antartique removes all land fractions souther than 60S.
+     !--- Land ice is set instead of the land fractions on these latitudes.
+     !--- The ocean and sea-ice fractions are not changed.
+     no_ter_antartique=.FALSE.
+     CALL getin_p('no_ter_antartique',no_ter_antartique)
+     WRITE(lunout,*)"no_ter_antartique=",no_ter_antartique
+     IF (no_ter_antartique) THEN
+        ! Remove all land fractions souther than 60S and set land-ice instead
+        WRITE(lunout,*) "Remove land fractions souther than 60deg south by increasing"
+        WRITE(lunout,*) "the continental ice fractions. No land can now be found at Antartic."
+        DO ji=1, klon
+           IF (latitude_deg(ji)<-60.0) THEN
+              pctsrf(ji,is_lic) = pctsrf(ji,is_lic) + pctsrf(ji,is_ter)
+              pctsrf(ji,is_ter) = 0
+           END IF
+        END DO
+     END IF
+     
+  ELSE
+     ! landice_opt=2 and higher
+     WRITE(lunout,*) 'No landice is attributed is_lic sub-surface because landice_opt=2 or higher.'
+     WRITE(lunout,*) 'This means that the land model will handel land ice as well as all other land areas.'
+     pctsrf(:,is_ter) = zmasq(:)
+     pctsrf(:,is_lic) = 0.0
+  END IF
+
+!--- 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(known_mask) pctsrf(:,is_oce)=1-zmasq(:)
+
+!--- 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'
+
+END SUBROUTINE start_init_subsurf
+!
+!-------------------------------------------------------------------------------
+
+END MODULE init_ssrf_m
+!
+!*******************************************************************************
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/test_disvert_m.F90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/test_disvert_m.F90	(revision 5267)
+++ 	(revision )
@@ -1,66 +1,0 @@
-module test_disvert_m
-
-  implicit none
-
-contains
-
-  subroutine test_disvert
-
-    ! Author: Lionel GUEZ
-
-    ! This procedure tests the order of pressure values at half-levels
-    ! and full levels. We arbitrarily choose to test ngrid values of
-    ! the surface pressure, which sample possible values on Earth.
-
-    use exner_hyb_m, only: exner_hyb
-    use vertical_layers_mod, only: ap,bp,preff
-    use comconst_mod, only: kappa, cpp
-
-    ! For llm:
-    include "dimensions.h"
-
-    ! Local:
-    integer l, i
-    integer, parameter:: ngrid = 7
-    real p(ngrid, llm + 1) ! pressure at half-level, in Pa
-    real pks(ngrid) ! exner function at the surface, in J K-1 kg-1 
-    real pk(ngrid, llm) ! exner function at full level, in J K-1 kg-1 
-    real ps(ngrid) ! surface pressure, in Pa
-    real p_lay(ngrid, llm) ! pressure at full level, in Pa
-    real delta_ps ! in Pa
-
-    !---------------------
-
-    print *, "Call sequence information: test_disvert"
-
-    delta_ps = 6e4 / (ngrid - 1)
-    ps = (/(5e4 + delta_ps * i, i = 0, ngrid - 1)/)
-    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
-    call exner_hyb(ngrid, ps, p, pks, pk)
-    p_lay = preff * (pk / cpp)**(1. / kappa)
-
-    ! Are pressure values in the right order?
-    if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) then
-       ! List details and stop:
-       do l = 1, llm
-          do i = 1, ngrid
-             if (p(i, l) <= p_lay(i, l)) then
-                print 1000, "ps = ", ps(i) / 100., "hPa, p(level ",  l, &
-                     ") = ", p(i, l) / 100., " hPa <= p(layer ", l, ") = ", &
-                     p_lay(i, l) / 100., " hPa"
-             end if
-             if (p_lay(i, l) <= p(i, l + 1)) then
-                print 1000, "ps = ", ps(i) / 100., "hPa, p(layer ", l, ") = ", &
-                     p_lay(i, l) / 100., " hPa <= p(level ", l + 1, ") = ", &
-                     p(i, l + 1) / 100., " hPa"
-             end if
-          end do
-       end do
-       call abort_physic("test_disvert", "bad order of pressure values", 1)
-    end if
-
-1000 format (3(a, g10.4, a, i0))
-
-  end subroutine test_disvert
-
-end module test_disvert_m
Index: LMDZ6/trunk/libf/dynphy_lonlat/phylmd/test_disvert_m.f90
===================================================================
--- LMDZ6/trunk/libf/dynphy_lonlat/phylmd/test_disvert_m.f90	(revision 5268)
+++ LMDZ6/trunk/libf/dynphy_lonlat/phylmd/test_disvert_m.f90	(revision 5268)
@@ -0,0 +1,66 @@
+module test_disvert_m
+
+  implicit none
+
+contains
+
+  subroutine test_disvert
+
+    ! Author: Lionel GUEZ
+
+    ! This procedure tests the order of pressure values at half-levels
+    ! and full levels. We arbitrarily choose to test ngrid values of
+    ! the surface pressure, which sample possible values on Earth.
+
+    use exner_hyb_m, only: exner_hyb
+    use vertical_layers_mod, only: ap,bp,preff
+    use comconst_mod, only: kappa, cpp
+
+    ! For llm:
+    include "dimensions.h"
+
+    ! Local:
+    integer l, i
+    integer, parameter:: ngrid = 7
+    real p(ngrid, llm + 1) ! pressure at half-level, in Pa
+    real pks(ngrid) ! exner function at the surface, in J K-1 kg-1 
+    real pk(ngrid, llm) ! exner function at full level, in J K-1 kg-1 
+    real ps(ngrid) ! surface pressure, in Pa
+    real p_lay(ngrid, llm) ! pressure at full level, in Pa
+    real delta_ps ! in Pa
+
+    !---------------------
+
+    print *, "Call sequence information: test_disvert"
+
+    delta_ps = 6e4 / (ngrid - 1)
+    ps = (/(5e4 + delta_ps * i, i = 0, ngrid - 1)/)
+    forall (l = 1: llm + 1) p(:, l) = ap(l) + bp(l) * ps
+    call exner_hyb(ngrid, ps, p, pks, pk)
+    p_lay = preff * (pk / cpp)**(1. / kappa)
+
+    ! Are pressure values in the right order?
+    if (any(p(:, :llm) <= p_lay .or. p_lay <= p(:, 2:))) then
+       ! List details and stop:
+       do l = 1, llm
+          do i = 1, ngrid
+             if (p(i, l) <= p_lay(i, l)) then
+                print 1000, "ps = ", ps(i) / 100., "hPa, p(level ",  l, &
+                     ") = ", p(i, l) / 100., " hPa <= p(layer ", l, ") = ", &
+                     p_lay(i, l) / 100., " hPa"
+             end if
+             if (p_lay(i, l) <= p(i, l + 1)) then
+                print 1000, "ps = ", ps(i) / 100., "hPa, p(layer ", l, ") = ", &
+                     p_lay(i, l) / 100., " hPa <= p(level ", l + 1, ") = ", &
+                     p(i, l + 1) / 100., " hPa"
+             end if
+          end do
+       end do
+       call abort_physic("test_disvert", "bad order of pressure values", 1)
+    end if
+
+1000 format (3(a, g10.4, a, i0))
+
+  end subroutine test_disvert
+
+end module test_disvert_m
