Index: LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90
===================================================================
--- LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2248)
+++ LMDZ5/trunk/libf/dynlonlat_phylonlat/phylmd/ce0l.F90	(revision 2293)
@@ -1,8 +1,4 @@
+PROGRAM ce0l
 !
-! $Id$
-!
-!-------------------------------------------------------------------------------
-!
-PROGRAM ce0l
 !-------------------------------------------------------------------------------
 ! Purpose: Calls etat0, creates initial states and limit_netcdf
@@ -15,68 +11,151 @@
 !-------------------------------------------------------------------------------
   USE control_mod
-#ifdef CPP_EARTH
-! This prog. is designed to work for Earth
+
+!#ifdef CPP_EARTH
+!! This prog. is designed to work for Earth
+  USE etat0dyn
+  USE netcdf, ONLY: NF90_OPEN, NF90_NOWRITE, NF90_CLOSE, NF90_NOERR
+  USE ioipsl, ONLY: ioconf_calendar, getin, flininfo, flinopen, flinget, flinclo
+
+#ifdef CPP_PHYS
+  USE etat0phys
   USE dimphy
   USE comgeomphy
   USE infotrac
   USE indice_sol_mod
+  USE test_disvert_m, ONLY: test_disvert
+#endif
+! of #ifdef CPP_PHYS
+
+!#endif
+!! of #ifdef CPP_EARTH
+
+  IMPLICIT NONE
+
+!-------------------------------------------------------------------------------
+! Local variables:
+!#ifndef CPP_EARTH
+!  include "iniprint.h"
+!!-------------------------------------------------------------------------------
+!  WRITE(lunout,*)'ce0l: Earth-specific routine, needs Earth physics'
+!#else
+  include "dimensions.h"
+  include "paramet.h"
+  include "comconst.h"
+  include "comvert.h"
+  include "iniprint.h"
+  include "temps.h"
+  include "logic.h"
+  REAL               :: masque(iip1,jjp1)             !--- CONTINENTAL MASK
+  REAL               :: phis  (iip1,jjp1)             !--- GROUND GEOPOTENTIAL
+  CHARACTER(LEN=256) :: modname, fmt, calnd           !--- CALENDAR TYPE
+  LOGICAL            :: use_filtre_fft
+  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
+
+!--- Local variables for ocean mask reading:
+  INTEGER            :: nid_o2a, iml_omask, jml_omask, j
+  INTEGER            :: fid, iret, llm_tmp, ttm_tmp, itaul(1)
+  REAL, ALLOCATABLE  :: lon_omask(:,:), dlon_omask(:), ocemask(:,:)
+  REAL, ALLOCATABLE  :: lat_omask(:,:), dlat_omask(:), ocetmp (:,:)
+  REAL               :: date, lev(1)
+!-------------------------------------------------------------------------------
+  modname="ce0l"
+
+!--- Constants
+  pi     = 4. * ATAN(1.)
+  rad    = 6371229.
+  daysec = 86400.
+  omeg   = 2.*pi/daysec
+  g      = 9.8
+  kappa  = 0.2857143
+  cpp    = 1004.70885
+  jmp1   = jjm + 1
+  preff   = 101325.
+  pa      = 50000.
+
+  CALL conf_gcm( 99, .TRUE. )
+
+  dtvr = daysec/FLOAT(day_step)
+  WRITE(lunout,*)'dtvr',dtvr
+
+  CALL iniconst()
+#ifdef CPP_PHYS
+  IF(pressure_exner) CALL test_disvert
+#endif
+  CALL inigeom()
 
 #ifdef CPP_IOIPSL
-  USE ioipsl, ONLY: ioconf_calendar, getin
-#else
-  ! if not using IOIPSL, we still need to use (a local version of) getin
-  use ioipsl_getincom, only: getin
-#endif
-
-#endif
-  IMPLICIT NONE
-#ifndef CPP_EARTH
-#include "iniprint.h"
-  WRITE(lunout,*)'limit_netcdf: Earth-specific routine, needs Earth physics'
-#else
-!-------------------------------------------------------------------------------
-! Local variables:
-  LOGICAL, PARAMETER :: interbar=.TRUE., extrap=.FALSE., oldice=.FALSE.
-#include "dimensions.h"
-#include "paramet.h"
-!#include "indicesol.h"
-#include "iniprint.h"
-#include "temps.h"
-#include "logic.h"
-  REAL,    DIMENSION(iip1,jjp1) :: masque
-  CHARACTER(LEN=15)             :: calnd
-  REAL,    DIMENSION(iip1,jjp1) :: phis ! geopotentiel au sol
-  logical use_filtre_fft
-!-------------------------------------------------------------------------------
-  CALL conf_gcm( 99, .TRUE. )
+  calnd='gregorian'
+  SELECT CASE(calend)
+    CASE('earth_360d');CALL ioconf_calendar('360d');   calnd='with 360 days/year'
+    CASE('earth_365d');CALL ioconf_calendar('noleap'); calnd='with no leap year'
+    CASE('earth_366d');CALL ioconf_calendar('366d');   calnd='with leap years only'
+    CASE('gregorian'); CALL ioconf_calendar('gregorian')
+    CASE('standard');  CALL ioconf_calendar('gregorian')
+    CASE('julian');    CALL ioconf_calendar('julian'); calnd='julian'
+    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
+  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
+    CASE DEFAULT
+      CALL abort_gcm('ce0l','Bad choice for calendar',1)
+  END SELECT
+  WRITE(lunout,*)'CHOSEN CALENDAR: Earth '//TRIM(calnd)
+#endif
 
   use_filtre_fft=.FALSE.
   CALL getin('use_filtre_fft',use_filtre_fft)
-  IF (use_filtre_fft) THEN
-     write(lunout, fmt = *) 'FFT filter is not available in the ' &
-          // 'sequential version of the dynamics.'
-     write(lunout, fmt = *) &
-          "Your setting of variable use_filtre_fft is not used."
+  IF(use_filtre_fft) THEN
+     WRITE(lunout,*)"FFT filter not available for sequential dynamics."
+     WRITE(lunout,*)"Your setting of variable use_filtre_fft is not used."
   ENDIF
 
+!--- LAND MASK. TWO CASES:
+!   1) read from ocean model    file "o2a.nc"    (coupled runs)
+!   2) computed from topography file "Relief.nc" (masque(:,:)=-99999.)
+! Coupled simulations (case 1) use the ocean model mask to compute the
+! weights to ensure ocean fractions are the same for atmosphere and ocean.
+!*******************************************************************************
+  IF(NF90_OPEN("o2a.nc", NF90_NOWRITE, nid_o2a)/=NF90_NOERR) THEN
+    WRITE(lunout,*)'BEWARE !! No ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Forced run.'
+    masque(:,:)=-99999.
+  ELSE
+    iret=NF90_CLOSE(nid_o2a)
+    WRITE(lunout,*)'BEWARE !! Ocean mask "o2a.nc" file found'
+    WRITE(lunout,*)'Coupled run.'
+    CALL flininfo("o2a.nc", iml_omask, jml_omask, llm_tmp, ttm_tmp, nid_o2a)
+    IF(iml_omask/=iim .OR.jml_omask/=jjp1) THEN
+      WRITE(lunout,*)'Mismatching dimensions for ocean mask'
+      WRITE(lunout,*)'iim  = ',iim ,' iml_omask = ',iml_omask
+      WRITE(lunout,*)'jjp1 = ',jjp1,' jml_omask = ',jml_omask
+      CALL abort_gcm(modname,'',1)
+    END IF
+    ALLOCATE(ocemask(iim,jjp1),lon_omask(iim,jjp1),dlon_omask(iim ))
+    ALLOCATE(ocetmp (iim,jjp1),lat_omask(iim,jjp1),dlat_omask(jjp1))
+    CALL flinopen("o2a.nc", .FALSE.,iim,jjp1,llm_tmp,lon_omask,lat_omask,lev,  &
+     &                                       ttm_tmp,itaul,date,dt,fid)
+    CALL flinget(fid, "OceMask",    iim,jjp1,llm_tmp,ttm_tmp,1,1,ocetmp)
+    CALL flinclo(fid)
+    dlon_omask(1:iim ) = lon_omask(1:iim,1)
+    dlat_omask(1:jjp1) = lat_omask(1,1:jjp1)
+    ocemask = ocetmp
+    IF(dlat_omask(1)<dlat_omask(jml_omask)) THEN
+      DO j=1,jjp1; ocemask(:,j) = ocetmp(:,jjp1-j+1); END DO
+    END IF
+    DEALLOCATE(ocetmp,lon_omask,lat_omask,dlon_omask,dlat_omask)
+    IF(prt_level>=1) THEN
+      WRITE(fmt,"(i4,'i1)')")iim ; fmt='('//ADJUSTL(fmt)
+      WRITE(lunout,*)'OCEAN MASK :'
+      WRITE(lunout,fmt) NINT(ocemask)
+    END IF
+    masque(1:iim,:)=1.-ocemask(:,:)
+    masque(iip1 ,:)=masque(1,:)
+    DEALLOCATE(ocemask)
+  END IF
+  phis(:,:)=-99999.
+
+#ifdef CPP_PHYS
   CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
   WRITE(lunout,*)'---> klon=',klon
   CALL InitComgeomphy
-
-#ifdef CPP_IOIPSL
-  SELECT CASE(calend)
-    CASE('earth_360d');CALL ioconf_calendar('360d');      calnd='a 360 jours/an'
-    CASE('earth_365d');CALL ioconf_calendar('noleap');    calnd='a 365 jours/an'
-    CASE('earth_366d');CALL ioconf_calendar('366d');      calnd='bissextile'
-    CASE('gregorian'); CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('standard');  CALL ioconf_calendar('gregorian'); calnd='gregorien'
-    CASE('julian');    CALL ioconf_calendar('julian');    calnd='julien'
-    CASE('proleptic_gregorian'); CALL ioconf_calendar('gregorian')
-  !--- DC Bof...  => IOIPSL a mettre a jour: proleptic_gregorian /= gregorian
-    CASE DEFAULT
-      CALL abort_gcm('ce0l','Mauvais choix de calendrier',1)
-  END SELECT
-  WRITE(lunout,*)'CALENDRIER CHOISI: Terrestre '//TRIM(calnd)
-#endif
 
   IF (type_trac == 'inca') THEN
@@ -86,24 +165,39 @@
     WRITE(lunout,*)'nbtr =' , nbtr 
 #endif
-  END IF
-
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  etat0_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) ' interbar = ',interbar
-  CALL etat0_netcdf(interbar,masque,phis,ok_etat0)
-
+! of #ifdef INCA
+  END IF
+  IF(ok_etat0) THEN
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0phy_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) ' interbar = ',interbar
+    CALL etat0phys_netcdf(interbar,masque,phis)
+  END IF
+#endif
+! of #ifdef CPP_PHYS
+
+  IF(ok_etat0) THEN
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,*) '  ***  etat0dyn_netcdf ***  '
+    WRITE(lunout,*) '  ************************  '
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) ' interbar = ',interbar
+    CALL etat0dyn_netcdf(interbar,masque,phis)
+  END IF
+
+#ifdef CPP_PHYS
   IF(ok_limit) THEN
-  WRITE(lunout,'(//)')
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,*) '  ***  Limit_netcdf ***  '
-  WRITE(lunout,*) '  *********************  '
-  WRITE(lunout,'(//)')
-  CALL limit_netcdf(interbar,extrap,oldice,masque)
-  END IF
-
-  
+    WRITE(lunout,'(//)')
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,*) '  ***  Limit_netcdf ***  '
+    WRITE(lunout,*) '  *********************  '
+    WRITE(lunout,'(//)')
+    CALL limit_netcdf(interbar,extrap,oldice,masque)
+  END IF
+#endif
+
   WRITE(lunout,'(//)')
   WRITE(lunout,*) '  ***************************  '
@@ -113,5 +207,5 @@
   CALL grilles_gcm_netcdf_sub(masque,phis)
 
-#endif
+!#endif
 ! of #ifndef CPP_EARTH #else
 
