Index: /LMDZ5/trunk/libf/phylmd/iniphysiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/iniphysiq.F	(revision 1850)
+++ /LMDZ5/trunk/libf/phylmd/iniphysiq.F	(revision 1851)
@@ -14,5 +14,5 @@
      &                               klon_omp_end,klon_mpi_begin
       USE comgeomphy, only : airephy,cuphy,cvphy,rlond,rlatd
-
+      USE phyaqua_mod, only: iniaqua
       IMPLICIT NONE
 c
Index: DZ5/trunk/libf/phylmd/phyaqua.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyaqua.F	(revision 1850)
+++ 	(revision )
@@ -1,909 +1,0 @@
-! Routines complementaires pour la physique planetaire.
-
-
-      subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  Creation d'un etat initial et de conditions aux limites 
-!  (resp startphy.nc et limit.nc) pour des configurations idealisees 
-! du modele LMDZ dans sa version terrestre.
-!  iflag_phys est un parametre qui controle
-!  iflag_phys = N  
-!    de 100 a 199 : aqua planetes avec SST forcees
-!                 N-100 determine le type de SSTs
-!    de 200 a 299 : terra planetes avec Ts calcule
-!        
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      USE comgeomphy, only : rlatd,rlond
-      USE dimphy, only : klon
-      USE surface_data, only : type_ocean,ok_veget
-      USE pbl_surface_mod, only : pbl_surface_init
-      USE fonte_neige_mod, only : fonte_neige_init
-      USE phys_state_var_mod
-      USE control_mod, only : dayref,nday,iphysiq
-      USE indice_sol_mod
-
-      USE IOIPSL 
-      IMPLICIT NONE
-
-#include "dimensions.h"
-!   #include "dimphy.h"
-!   #include "YOMCST.h"
-#include "comconst.h"
-#include "clesphys.h"
-#include "dimsoil.h"
-#include "temps.h"
-
-      integer,intent(in) :: nlon,iflag_phys
-cIM ajout latfi, lonfi
-      real,intent(in) :: lonfi(nlon),latfi(nlon)
-
-      INTEGER type_profil,type_aqua
-
-c  Ajouts initialisation des surfaces
-      REAL :: run_off_lic_0(nlon)
-      REAL :: qsolsrf(nlon,nbsrf),snsrf(nlon,nbsrf)
-      REAL :: frugs(nlon,nbsrf)
-      REAL :: agesno(nlon,nbsrf)
-      REAL :: tsoil(nlon,nsoilmx,nbsrf)
-      REAL :: tslab(nlon), seaice(nlon)
-      REAL evap(nlon,nbsrf),fder(nlon)
-
-
-
-c    Arguments :
-c    -----------
-
-!      integer radpas  
-      integer it,unit,i,k,itap
-
-      real airefi,zcufi,zcvfi
-
-      real rugos,albedo
-      REAL tsurf
-      REAL time,timestep,day,day0
-      real qsol_f,qsol(nlon)
-      real rugsrel(nlon)
-!      real zmea(nlon),zstd(nlon),zsig(nlon)
-!      real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
-!      real rlon(nlon),rlat(nlon)
-      logical alb_ocean
-!      integer demih_pas
-
-      CHARACTER*80 ans,file_forctl, file_fordat, file_start
-      character*100 file,var
-      character*2 cnbl
-
-      REAL phy_nat(nlon,360)
-      REAL phy_alb(nlon,360)
-      REAL phy_sst(nlon,360)
-      REAL phy_bil(nlon,360)
-      REAL phy_rug(nlon,360)
-      REAL phy_ice(nlon,360)
-      REAL phy_fter(nlon,360)
-      REAL phy_foce(nlon,360)
-      REAL phy_fsic(nlon,360)
-      REAL phy_flic(nlon,360)
-
-      integer, save::  read_climoz=0 ! read ozone climatology
-
-! intermediate variables to use getin
-      integer :: nbapp_rad_omp
-      real :: co2_ppm_omp,solaire_omp
-      logical :: alb_ocean_omp
-      real :: rugos_omp
-!-------------------------------------------------------------------------
-!  declaration pour l'appel a phyredem
-!-------------------------------------------------------------------------
-
-!      real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
-      real falbe(nlon,nbsrf),falblw(nlon,nbsrf)
-!      real pbl_tke(nlon,llm,nbsrf)
-!      real rain_fall(nlon),snow_fall(nlon)
-!      real solsw(nlon), sollw(nlon),radsol(nlon)
-!      real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
-!      real ratqs(nlon,llm)
-!      real clwcon(nlon,llm)
-
-      INTEGER        longcles
-      PARAMETER    ( longcles = 20 )
-      REAL clesphy0( longcles )
-
-
-c-----------------------------------------------------------------------
-c   dynamial tendencies :
-c   ---------------------
-
-      INTEGER l,ierr,aslun
-
-      REAL longitude,latitude
-      REAL paire 
-
-      DATA latitude,longitude/48.,0./
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! INITIALISATIONS
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-!-----------------------------------------------------------------------
-!    Initialisations  des constantes
-!    -------------------------------
-
-
-      type_aqua=iflag_phys/100
-      type_profil=iflag_phys-type_aqua*100
-      print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
-
-      if (klon.ne.nlon) then
-        write(*,*)"iniaqua: klon=",klon," nlon=",nlon
-        stop'probleme de dimensions dans iniaqua'
-      endif
-      call phys_state_var_init(read_climoz)
-
-
-      read_climoz=0
-      day0=217.
-      day=day0
-      it=0
-      time=0.
-
-cIM ajout latfi, lonfi
-      rlatd=latfi
-      rlond=lonfi
-      rlat=rlatd*180./pi
-      rlon=rlond*180./pi
-
-!-----------------------------------------------------------------------
-!  initialisations de la physique
-!-----------------------------------------------------------------------
-
-         day_ini=dayref
-         day_end=day_ini+nday
-         airefi=1.
-         zcufi=1.
-         zcvfi=1.
-!$OMP MASTER
-      nbapp_rad_omp=24
-      CALL getin('nbapp_rad',nbapp_rad_omp)
-!$OMP END MASTER
-!$OMP BARRIER
-      nbapp_rad=nbapp_rad_omp
-
-!---------------------------------------------------------------------
-c Creation des conditions aux limites:
-c ------------------------------------
-! Initialisations des constantes
-! Ajouter les manquants dans planete.def... (albedo etc)
-!$OMP MASTER
-      co2_ppm_omp=348.
-      CALL getin('co2_ppm',co2_ppm_omp)
-      solaire_omp=1365.
-      CALL getin('solaire',solaire_omp)
-!      CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua
-      alb_ocean_omp=.true.
-      CALL getin('alb_ocean',alb_ocean_omp)
-!$OMP END MASTER
-!$OMP BARRIER
-      co2_ppm=co2_ppm_omp
-      solaire=solaire_omp
-      alb_ocean=alb_ocean_omp
-
-      radsol=0.
-      qsol_f=10.
-
-c  Conditions aux limites:
-c  -----------------------
-
-      qsol(:)    = qsol_f
-      rugsrel = 0.0    ! (rugsrel = rugoro)
-      rugoro = 0.0
-      u_ancien = 0.0
-      v_ancien = 0.0
-      agesno  = 50.0
-! Relief plat
-      zmea = 0.
-      zstd = 0.
-      zsig = 0.
-      zgam = 0.
-      zthe = 0.
-      zpic = 0.
-      zval = 0.
-
-! Une seule surface
-      pctsrf=0.
-      if (type_aqua==1) then
-         rugos=1.e-4
-         albedo=0.19
-         pctsrf(:,is_oce)=1.
-      else if (type_aqua==2) then
-         rugos=0.03
-         albedo=0.1
-         pctsrf(:,is_ter)=1.
-      endif
-
-!$OMP MASTER
-      rugos_omp=rugos
-      CALL getin('rugos',rugos_omp)
-!$OMP END MASTER
-!$OMP BARRIER
-      rugos=rugos_omp
-      zmasq(:)=pctsrf(:,is_oce)
-
-! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
-! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
-
-! Si alb_ocean on calcule un albedo oceanique moyen
-!  if (alb_ocean) then
-! Voir pourquoi on avait ca.
-!          CALL ini_alb_oce(phy_alb)
-!      else 
-      phy_alb(:,:) = albedo ! albedo land only (old value condsurf_jyg=0.3)
-!      endif !alb_ocean
-      
-      do i=1,360
-cIM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
-cIM ajout calcul profil sst selon le cas considere (cf. FBr)
-
-      phy_nat(:,i) = 1.0    ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
-      phy_bil(:,i) = 1.0    ! ne sert que pour les slab_ocean
-      phy_rug(:,i) = rugos  ! longueur rugosite utilisee sur land only 
-      phy_ice(:,i) = 0.0    ! fraction de glace (?)
-      phy_fter(:,i) = pctsrf(:,is_ter)  ! fraction de glace (?)
-      phy_foce(:,i) = pctsrf(:,is_oce)  ! fraction de glace (?)
-      phy_fsic(:,i) = pctsrf(:,is_sic)  ! fraction de glace (?)
-      phy_flic(:,i) = pctsrf(:,is_lic)  ! fraction de glace (?)
-      enddo
-cIM calcul profil sst
-      call profil_sst(nlon, rlatd, type_profil, phy_sst)
-
-      call writelim
-     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
-     s    phy_fter,phy_foce,phy_flic,phy_fsic)
-
-
-!---------------------------------------------------------------------
-c Ecriture de l'etat initial:
-c ---------------------------
-
-C
-C Ecriture etat initial physique
-C
-      timestep   = dtvr * FLOAT(iphysiq)
-      radpas    = NINT (daysec/timestep/ FLOAT(nbapp_rad) )
-
-      DO i = 1, longcles
-       clesphy0(i) = 0.
-      ENDDO
-      clesphy0(1) = FLOAT( iflag_con )
-      clesphy0(2) = FLOAT( nbapp_rad )
-c     IF( cycle_diurne  ) clesphy0(3) =  1. 
-      clesphy0(3)=1. ! cycle_diurne
-      clesphy0(4)=1. ! soil_model
-      clesphy0(5)=1. ! new_oliq
-      clesphy0(6)=0. ! ok_orodr
-      clesphy0(7)=0. ! ok_orolf
-      clesphy0(8)=0. ! ok_limitvrai
-
-
-c=======================================================================
-c  Profils initiaux
-c=======================================================================
-
-! On initialise les temperatures de surfaces comme les sst
-      do i=1,nlon
-         ftsol(i,:)=phy_sst(i,1)
-         tsoil(i,:,:)=phy_sst(i,1)
-         tslab(i)=phy_sst(i,1)
-      enddo
-
-      falbe(:,:)=albedo
-      falblw(:,:)=albedo
-      rain_fall(:)=0.
-      snow_fall(:)=0.
-      solsw(:)=0.
-      sollw(:)=0.
-      radsol(:)=0.
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!  intialisation bidon mais pas grave
-      t_ancien(:,:)=0.
-      q_ancien(:,:)=0.
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-      rnebcon=0.
-      ratqs=0.
-      clwcon=0.
-      pbl_tke=1.e-8
-
-! variables supplementaires pour appel a plb_surface_init
-      fder(:)=0.
-      seaice(:)=0.
-      run_off_lic_0=0.
-      evap=0.
-
-
-! Initialisations necessaires avant phyredem
-      type_ocean = "force"
-      call fonte_neige_init(run_off_lic_0)
-      qsolsrf(:,:)=qsol(1) ! humidite du sol des sous surface
-      snsrf(:,:)=0.        ! couverture de neige des sous surface
-      frugs(:,:)=rugos        ! couverture de neige des sous surface
-
-
-      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
-     .     evap, frugs, agesno, tsoil)
-
-        print*,'iniaqua: before phyredem'
-
-      falb1=albedo
-      falb2=albedo
-      zmax0=0.
-      f0=0.
-      sig1=0.
-      w01=0.
-      wake_deltat=0.
-      wake_deltaq=0.
-      wake_s=0.
-      wake_cstar=0.
-      wake_pe=0.
-      wake_fip=0.
-      fm_therm=0.
-      entr_therm=0.
-      detr_therm=0.
-
-
-      CALL phyredem ("startphy.nc")
-
-        print*,'iniaqua: after phyredem'
-      call phys_state_var_end
-
-      return
-      end
-
-
-c====================================================================
-c====================================================================
-      SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
-      USE dimphy
-      IMPLICIT none
-c====================================================================
-c=============================================================
-c         CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
-c Auteur : A. Campoy et F. Hourdin
-c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
-c          et l'ensoleillement moyen entre gmtime1 et gmtime2 
-c          connaissant la declinaison, la latitude et la longitude.
-c 
-c   Dans cette version particuliere, on calcule le rayonnement
-c  moyen sur l'année à chaque latitude.
-c   angle zenithal calculé pour obtenir un 
-c   Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
-c   en moyenne annuelle.
-c   Spécifique de la terre. Utilisé pour les aqua planetes.
-c
-c Rque   : Different de la routine angle en ce sens que zenang 
-c          fournit des moyennes de pmu0 et non des valeurs 
-c          instantanees, du coup frac prend toutes les valeurs 
-c          entre 0 et 1.
-c Date   : premiere version le 13 decembre 1994
-c          revu pour  GCM  le 30 septembre 1996
-c===============================================================
-c longi----INPUT : la longitude vraie de la terre dans son plan
-c                  solaire a partir de l'equinoxe de printemps (degre)
-c gmtime---INPUT : temps universel en fraction de jour
-c pdtrad---INPUT : pas de temps du rayonnement (secondes)
-c lat------INPUT : latitude en degres
-c long-----INPUT : longitude en degres
-c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
-c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
-c================================================================
-#include "YOMCST.h"
-c================================================================
-      logical cycle_diurne
-      real  gmtime
-      real rlat(klon), rlon(klon), rmu0(klon), fract(klon)
-c================================================================
-      integer i
-      real gmtime1, gmtime2
-      real pi_local
-
-
-      real rmu0m(klon),rmu0a(klon)
-c
-
-      pi_local = 4.0 * ATAN(1.0)
-
-c================================================================
-c  Calcul de l'angle zenithal moyen sur la journee
-c================================================================
-
-      DO i=1,klon
-        fract(i)=1. 
-!  Calcule du flux moyen
-        IF (abs(rlat(i)).LE.28.75) THEN
-           rmu0m(i)=(210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
-        ELSEIF (abs(rlat(i)).LE.43.75) THEN
-          rmu0m(i)=(187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
-        ELSEIF (abs(rlat(i)).LE.71.25) THEN
-          rmu0m(i)=(162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
-        ELSE
-          rmu0m(i)=(172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
-        ENDIF
-      ENDDO
-
-c================================================================
-!  Avec ou sans cycle diurne
-c================================================================
-
-      IF (cycle_diurne) THEN
-
-!  On redecompose flux  au sommet suivant un cycle diurne idealise
-!  identique a toutes les latitudes.
-
-         DO i=1,klon
-           rmu0a(i)=2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
-           rmu0(i)=rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*
-     &      rlon(i)/360.))-rmu0a(i)/sqrt(2.)
-         ENDDO
-
-         DO i=1,klon
-           IF (rmu0(i).LE.0.) THEN
-              rmu0(i)=0.
-              fract(i)=0.
-           ELSE
-              fract(i)=1.
-           ENDIF
-         ENDDO
-
-!  Affichage de l'angel zenitale
-!     print*,'************************************'
-!     print*,'************************************'
-!     print*,'************************************'
-!     print*,'latitude=',rlat(i),'longitude=',rlon(i)
-!     print*,'rmu0m=',rmu0m(i)
-!     print*,'rmu0a=',rmu0a(i)
-!     print*,'rmu0=',rmu0(i)
-                                                               
-      ELSE
-
-        DO i=1,klon
-           fract(i)=0.5
-           rmu0(i)=rmu0m(i)*2.
-        ENDDO
-
-      ENDIF
-
-      RETURN
-      END
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      subroutine writelim
-     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
-     s    phy_fter,phy_foce,phy_flic,phy_fsic)
-c
-      use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root
-      use mod_grid_phy_lmdz, only : klon_glo
-      use mod_phys_lmdz_transfert_para, only : gather
-!#include "dimensions.h"
-!#include "dimphy.h"
-#include "netcdf.inc"
- 
-      integer,intent(in) :: klon
-      real,intent(in) :: phy_nat(klon,360)
-      real,intent(in) :: phy_alb(klon,360)
-      real,intent(in) :: phy_sst(klon,360)
-      real,intent(in) :: phy_bil(klon,360)
-      real,intent(in) :: phy_rug(klon,360)
-      real,intent(in) :: phy_ice(klon,360)
-      real,intent(in) :: phy_fter(klon,360)
-      real,intent(in) :: phy_foce(klon,360)
-      real,intent(in) :: phy_flic(klon,360)
-      real,intent(in) :: phy_fsic(klon,360)
-
-      real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:)
-                                    ! on the whole physics grid 
-      INTEGER ierr
-      INTEGER dimfirst(3)
-      INTEGER dimlast(3)
-c
-      INTEGER nid, ndim, ntim
-      INTEGER dims(2), debut(2), epais(2)
-      INTEGER id_tim
-      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
-      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
- 
-      if (is_mpi_root.and.is_omp_root) then
-      
-        PRINT*, 'writelim: Ecriture du fichier limit'
-c
-        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
-c
-        ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
-     .                       "Fichier conditions aux limites")
-!!        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
-        ierr = NF_DEF_DIM (nid, "points_physiques", klon_glo, ndim)
-        ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
-c
-        dims(1) = ndim
-        dims(2) = ntim
-c
-ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
-        ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
-        ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
-     .                        "Jour dans l annee")
-ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
-        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
-        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
-     .                        "Nature du sol (0,1,2,3)")
-ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
-        ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
-        ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
-     .                        "Temperature superficielle de la mer")
-ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
-        ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
-        ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
-     .                        "Reference flux de chaleur au sol")
-ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
-        ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
-        ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
-     .                        "Albedo a la surface")
-ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
-        ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
-        ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
-     .                        "Rugosite")
-
-        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
-        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
-        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
-        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
-        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
-        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
-        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
-        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
-c
-        ierr = NF_ENDDEF(nid)
-c
-
-! write the 'times'
-        do k=1,360
-#ifdef NC_DOUBLE
-          ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
-#else
-          ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
-#endif
-        enddo
-
-      endif ! of if (is_mpi_root.and.is_omp_root)
-
-! write the fields, after having collected them on master
-
-      call gather(phy_nat,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_nat"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_sst,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_sst"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_bil,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_bil"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_alb,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_alb"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_rug,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_rug"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_fter,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_fter"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_foce,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_foce"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_fsic,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_fsic"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-      call gather(phy_flic,phy_glo)
-      if (is_mpi_root.and.is_omp_root) then
-#ifdef NC_DOUBLE
-        ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo)
-#else
-        ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo)
-#endif
-        if(ierr.ne.NF_NOERR) then
-          write(*,*) "writelim error with phy_flic"
-          write(*,*) NF_STRERROR(ierr)
-        endif
-      endif
-
-!  close file:
-      if (is_mpi_root.and.is_omp_root) then
-        ierr = NF_CLOSE(nid)
-      endif
-
-      end
-
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-
-      SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
-      use dimphy
-      IMPLICIT none
-c
-      INTEGER nlon, type_profil, i, k, j
-      REAL :: rlatd(nlon), phy_sst(nlon, 360) 
-      INTEGER imn, imx, amn, amx, kmn, kmx
-      INTEGER p, pplus, nlat_max
-      parameter (nlat_max=72)
-      REAL x_anom_sst(nlat_max)
-c
-      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
-      do i=1,360
-c      phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
-
-c Rajout fbrlmd
-
-      if(type_profil.EQ.1)then
-c     Méthode 1 "Control" faible plateau à l'Equateur
-      do j=1,klon
-       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**2)
-c        PI/3=1.047197551
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-      if(type_profil.EQ.2)then
-c     Méthode 2 "Flat" fort plateau à l'Equateur
-      do j=1,klon
-       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-
-      if (type_profil.EQ.3) then
-c     Méthode 3 "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.4) then
-c     Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.5) then
-c     Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.+2.
-        endif
-
-      enddo
-      endif
-
-      if(type_profil.EQ.6)then
-c     Méthode 6 "cst" valeur constante de SST
-      do j=1,klon
-       phy_sst(j,i)=288.
-      enddo
-      endif
-
-
-        if(type_profil.EQ.7)then
-c     Méthode 7 "cst" valeur constante de SST +2
-      do j=1,klon
-       phy_sst(j,i)=288.+2.
-      enddo
-      endif
-
-        p=0
-        if(type_profil.EQ.8)then
-c     Méthode 8 profil anomalies SST du modèle couplé AR4
-       do j=1,klon
-         if (rlatd(j).EQ.rlatd(j-1)) then
-       phy_sst(j,i)=273.+x_anom_sst(pplus)
-     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
-          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-            phy_sst(j,i)=273.+x_anom_sst(pplus)
-          endif
-        else
-          p=p+1
-          pplus=73-p
-        phy_sst(j,i)=273.+x_anom_sst(pplus)
-     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
-          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-            phy_sst(j,i)=273.+x_anom_sst(pplus)
-          endif
-          write (*,*) rlatd(j),x_anom_sst(pplus),phy_sst(j,i)
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.9) then
-c     Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.-2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.-2.
-        endif
-      enddo
-      endif
-
-
-      if (type_profil.EQ.10) then
-c     Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.+4.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.11) then
-c     Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.12) then
-c     Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
-      do j=1,klon
-        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.+4.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.13) then
-c     Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
-      do j=1,klon
-        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-      if (type_profil.EQ.14) then
-c     Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
-      do j=1,klon
-        phy_sst(j,i)=273.+2.+0.5*29.*(2-sin(1.5*rlatd(j))**2
-     &                 -sin(1.5*rlatd(j))**4)
-      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
-         phy_sst(j,i)=273.
-        endif
-      enddo
-      endif
-
-      enddo
-
-cIM beg : verif profil SST: phy_sst
-       amn=MIN(phy_sst(1,1),1000.)
-       amx=MAX(phy_sst(1,1),-1000.)
-       DO k=1, 360
-       DO i=2, nlon
-        IF(phy_sst(i,k).LT.amn) THEN
-         amn=phy_sst(i,k)
-         imn=i
-         kmn=k
-        ENDIF
-        IF(phy_sst(i,k).GT.amx) THEN
-         amx=phy_sst(i,k)
-         imx=i
-         kmx=k
-        ENDIF
-       ENDDO
-       ENDDO
-c
-       PRINT*,' debut avant writelim min max phy_sst',imn,kmn,amn,
-     & imx,kmx,amx
-cIM end : verif profil SST: phy_sst
-
-       return 
-       end
Index: /LMDZ5/trunk/libf/phylmd/phyaqua_mod.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/phyaqua_mod.F	(revision 1851)
+++ /LMDZ5/trunk/libf/phylmd/phyaqua_mod.F	(revision 1851)
@@ -0,0 +1,930 @@
+      module phyaqua_mod
+! Routines complementaires pour la physique planetaire.
+      implicit none
+
+      contains
+
+      subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  Creation d'un etat initial et de conditions aux limites 
+!  (resp startphy.nc et limit.nc) pour des configurations idealisees 
+! du modele LMDZ dans sa version terrestre.
+!  iflag_phys est un parametre qui controle
+!  iflag_phys = N  
+!    de 100 a 199 : aqua planetes avec SST forcees
+!                 N-100 determine le type de SSTs
+!    de 200 a 299 : terra planetes avec Ts calcule
+!        
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      USE comgeomphy, only : rlatd,rlond
+      USE dimphy, only : klon
+      USE surface_data, only : type_ocean,ok_veget
+      USE pbl_surface_mod, only : pbl_surface_init
+      USE fonte_neige_mod, only : fonte_neige_init
+      USE phys_state_var_mod
+      USE control_mod, only : dayref,nday,iphysiq
+      USE indice_sol_mod
+
+      USE IOIPSL 
+      IMPLICIT NONE
+
+#include "dimensions.h"
+!   #include "dimphy.h"
+!   #include "YOMCST.h"
+#include "comconst.h"
+#include "clesphys.h"
+#include "dimsoil.h"
+#include "temps.h"
+
+      integer,intent(in) :: nlon,iflag_phys
+cIM ajout latfi, lonfi
+      real,intent(in) :: lonfi(nlon),latfi(nlon)
+
+      INTEGER type_profil,type_aqua
+
+c  Ajouts initialisation des surfaces
+      REAL :: run_off_lic_0(nlon)
+      REAL :: qsolsrf(nlon,nbsrf),snsrf(nlon,nbsrf)
+      REAL :: frugs(nlon,nbsrf)
+      REAL :: agesno(nlon,nbsrf)
+      REAL :: tsoil(nlon,nsoilmx,nbsrf)
+      REAL :: tslab(nlon), seaice(nlon)
+      REAL evap(nlon,nbsrf),fder(nlon)
+
+
+
+c    Arguments :
+c    -----------
+
+!      integer radpas  
+      integer it,unit,i,k,itap
+
+      real airefi,zcufi,zcvfi
+
+      real rugos,albedo
+      REAL tsurf
+      REAL time,timestep,day,day0
+      real qsol_f,qsol(nlon)
+      real rugsrel(nlon)
+!      real zmea(nlon),zstd(nlon),zsig(nlon)
+!      real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
+!      real rlon(nlon),rlat(nlon)
+      logical alb_ocean
+!      integer demih_pas
+
+      CHARACTER*80 ans,file_forctl, file_fordat, file_start
+      character*100 file,var
+      character*2 cnbl
+
+      REAL phy_nat(nlon,360)
+      REAL phy_alb(nlon,360)
+      REAL phy_sst(nlon,360)
+      REAL phy_bil(nlon,360)
+      REAL phy_rug(nlon,360)
+      REAL phy_ice(nlon,360)
+      REAL phy_fter(nlon,360)
+      REAL phy_foce(nlon,360)
+      REAL phy_fsic(nlon,360)
+      REAL phy_flic(nlon,360)
+
+      integer, save::  read_climoz=0 ! read ozone climatology
+
+! intermediate variables to use getin (need to be "save" to be shared by all threads)
+      integer,save :: nbapp_rad_omp
+      real,save :: co2_ppm_omp,solaire_omp
+      logical,save :: alb_ocean_omp
+      real,save :: rugos_omp
+!-------------------------------------------------------------------------
+!  declaration pour l'appel a phyredem
+!-------------------------------------------------------------------------
+
+!      real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
+      real falbe(nlon,nbsrf),falblw(nlon,nbsrf)
+!      real pbl_tke(nlon,llm,nbsrf)
+!      real rain_fall(nlon),snow_fall(nlon)
+!      real solsw(nlon), sollw(nlon),radsol(nlon)
+!      real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
+!      real ratqs(nlon,llm)
+!      real clwcon(nlon,llm)
+
+      INTEGER        longcles
+      PARAMETER    ( longcles = 20 )
+      REAL clesphy0( longcles )
+
+
+c-----------------------------------------------------------------------
+c   dynamial tendencies :
+c   ---------------------
+
+      INTEGER l,ierr,aslun
+
+      REAL longitude,latitude
+      REAL paire 
+
+      DATA latitude,longitude/48.,0./
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! INITIALISATIONS
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+!-----------------------------------------------------------------------
+!    Initialisations  des constantes
+!    -------------------------------
+
+
+      type_aqua=iflag_phys/100
+      type_profil=iflag_phys-type_aqua*100
+      print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
+
+      if (klon.ne.nlon) then
+        write(*,*)"iniaqua: klon=",klon," nlon=",nlon
+        stop'probleme de dimensions dans iniaqua'
+      endif
+      call phys_state_var_init(read_climoz)
+
+
+      read_climoz=0
+      day0=217.
+      day=day0
+      it=0
+      time=0.
+
+cIM ajout latfi, lonfi
+      rlatd=latfi
+      rlond=lonfi
+      rlat=rlatd*180./pi
+      rlon=rlond*180./pi
+
+!-----------------------------------------------------------------------
+!  initialisations de la physique
+!-----------------------------------------------------------------------
+
+         day_ini=dayref
+         day_end=day_ini+nday
+         airefi=1.
+         zcufi=1.
+         zcvfi=1.
+!$OMP MASTER
+      nbapp_rad_omp=24
+      CALL getin('nbapp_rad',nbapp_rad_omp)
+!$OMP END MASTER
+!$OMP BARRIER
+      nbapp_rad=nbapp_rad_omp
+
+!---------------------------------------------------------------------
+c Creation des conditions aux limites:
+c ------------------------------------
+! Initialisations des constantes
+! Ajouter les manquants dans planete.def... (albedo etc)
+!$OMP MASTER
+      co2_ppm_omp=348.
+      CALL getin('co2_ppm',co2_ppm_omp)
+      solaire_omp=1365.
+      CALL getin('solaire',solaire_omp)
+!      CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua
+      alb_ocean_omp=.true.
+      CALL getin('alb_ocean',alb_ocean_omp)
+!$OMP END MASTER
+!$OMP BARRIER
+      co2_ppm=co2_ppm_omp
+      write(*,*)"iniaqua: co2_ppm=",co2_ppm
+      solaire=solaire_omp
+      write(*,*)"iniaqua: solaire=",solaire
+      alb_ocean=alb_ocean_omp
+      write(*,*)"iniaqua: alb_ocean=",alb_ocean
+
+      radsol=0.
+      qsol_f=10.
+
+c  Conditions aux limites:
+c  -----------------------
+
+      qsol(:)    = qsol_f
+      rugsrel = 0.0    ! (rugsrel = rugoro)
+      rugoro = 0.0
+      u_ancien = 0.0
+      v_ancien = 0.0
+      agesno  = 50.0
+! Relief plat
+      zmea = 0.
+      zstd = 0.
+      zsig = 0.
+      zgam = 0.
+      zthe = 0.
+      zpic = 0.
+      zval = 0.
+
+! Une seule surface
+      pctsrf=0.
+      if (type_aqua==1) then
+         rugos=1.e-4
+         albedo=0.19
+         pctsrf(:,is_oce)=1.
+      else if (type_aqua==2) then
+         rugos=0.03
+         albedo=0.1
+         pctsrf(:,is_ter)=1.
+      endif
+
+!$OMP MASTER
+      rugos_omp=rugos
+      CALL getin('rugos',rugos_omp)
+!$OMP END MASTER
+!$OMP BARRIER
+      rugos=rugos_omp
+      write(*,*) "iniaqua: rugos=",rugos
+      zmasq(:)=pctsrf(:,is_oce)
+
+! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
+! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
+
+! Si alb_ocean on calcule un albedo oceanique moyen
+!  if (alb_ocean) then
+! Voir pourquoi on avait ca.
+!          CALL ini_alb_oce(phy_alb)
+!      else 
+      phy_alb(:,:) = albedo ! albedo land only (old value condsurf_jyg=0.3)
+!      endif !alb_ocean
+      
+      do i=1,360
+cIM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
+cIM ajout calcul profil sst selon le cas considere (cf. FBr)
+
+      phy_nat(:,i) = 1.0    ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
+      phy_bil(:,i) = 1.0    ! ne sert que pour les slab_ocean
+      phy_rug(:,i) = rugos  ! longueur rugosite utilisee sur land only 
+      phy_ice(:,i) = 0.0    ! fraction de glace (?)
+      phy_fter(:,i) = pctsrf(:,is_ter)  ! fraction de glace (?)
+      phy_foce(:,i) = pctsrf(:,is_oce)  ! fraction de glace (?)
+      phy_fsic(:,i) = pctsrf(:,is_sic)  ! fraction de glace (?)
+      phy_flic(:,i) = pctsrf(:,is_lic)  ! fraction de glace (?)
+      enddo
+cIM calcul profil sst
+      call profil_sst(nlon, rlatd, type_profil, phy_sst)
+
+      call writelim
+     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
+     s    phy_fter,phy_foce,phy_flic,phy_fsic)
+
+
+!---------------------------------------------------------------------
+c Ecriture de l'etat initial:
+c ---------------------------
+
+C
+C Ecriture etat initial physique
+C
+      timestep   = dtvr * FLOAT(iphysiq)
+      radpas    = NINT (daysec/timestep/ FLOAT(nbapp_rad) )
+
+      DO i = 1, longcles
+       clesphy0(i) = 0.
+      ENDDO
+      clesphy0(1) = FLOAT( iflag_con )
+      clesphy0(2) = FLOAT( nbapp_rad )
+c     IF( cycle_diurne  ) clesphy0(3) =  1. 
+      clesphy0(3)=1. ! cycle_diurne
+      clesphy0(4)=1. ! soil_model
+      clesphy0(5)=1. ! new_oliq
+      clesphy0(6)=0. ! ok_orodr
+      clesphy0(7)=0. ! ok_orolf
+      clesphy0(8)=0. ! ok_limitvrai
+
+
+c=======================================================================
+c  Profils initiaux
+c=======================================================================
+
+! On initialise les temperatures de surfaces comme les sst
+      do i=1,nlon
+         ftsol(i,:)=phy_sst(i,1)
+         tsoil(i,:,:)=phy_sst(i,1)
+         tslab(i)=phy_sst(i,1)
+      enddo
+
+      falbe(:,:)=albedo
+      falblw(:,:)=albedo
+      rain_fall(:)=0.
+      snow_fall(:)=0.
+      solsw(:)=0.
+      sollw(:)=0.
+      radsol(:)=0.
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!  intialisation bidon mais pas grave
+      t_ancien(:,:)=0.
+      q_ancien(:,:)=0.
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+      rnebcon=0.
+      ratqs=0.
+      clwcon=0.
+      pbl_tke=1.e-8
+
+! variables supplementaires pour appel a plb_surface_init
+      fder(:)=0.
+      seaice(:)=0.
+      run_off_lic_0=0.
+      evap=0.
+
+
+! Initialisations necessaires avant phyredem
+      type_ocean = "force"
+      call fonte_neige_init(run_off_lic_0)
+      qsolsrf(:,:)=qsol(1) ! humidite du sol des sous surface
+      snsrf(:,:)=0.        ! couverture de neige des sous surface
+      frugs(:,:)=rugos        ! couverture de neige des sous surface
+
+
+      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,
+     .     evap, frugs, agesno, tsoil)
+
+        print*,'iniaqua: before phyredem'
+
+      falb1=albedo
+      falb2=albedo
+      zmax0=0.
+      f0=0.
+      sig1=0.
+      w01=0.
+      wake_deltat=0.
+      wake_deltaq=0.
+      wake_s=0.
+      wake_cstar=0.
+      wake_pe=0.
+      wake_fip=0.
+      fm_therm=0.
+      entr_therm=0.
+      detr_therm=0.
+
+
+      CALL phyredem ("startphy.nc")
+
+        print*,'iniaqua: after phyredem'
+      call phys_state_var_end
+
+      return
+      end subroutine iniaqua
+
+
+c====================================================================
+c====================================================================
+      SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
+      USE dimphy
+      IMPLICIT none
+c====================================================================
+c=============================================================
+c         CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
+c Auteur : A. Campoy et F. Hourdin
+c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
+c          et l'ensoleillement moyen entre gmtime1 et gmtime2 
+c          connaissant la declinaison, la latitude et la longitude.
+c 
+c   Dans cette version particuliere, on calcule le rayonnement
+c  moyen sur l'année à chaque latitude.
+c   angle zenithal calculé pour obtenir un 
+c   Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
+c   en moyenne annuelle.
+c   Spécifique de la terre. Utilisé pour les aqua planetes.
+c
+c Rque   : Different de la routine angle en ce sens que zenang 
+c          fournit des moyennes de pmu0 et non des valeurs 
+c          instantanees, du coup frac prend toutes les valeurs 
+c          entre 0 et 1.
+c Date   : premiere version le 13 decembre 1994
+c          revu pour  GCM  le 30 septembre 1996
+c===============================================================
+c longi----INPUT : la longitude vraie de la terre dans son plan
+c                  solaire a partir de l'equinoxe de printemps (degre)
+c gmtime---INPUT : temps universel en fraction de jour
+c pdtrad---INPUT : pas de temps du rayonnement (secondes)
+c lat------INPUT : latitude en degres
+c long-----INPUT : longitude en degres
+c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
+c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
+c================================================================
+#include "YOMCST.h"
+c================================================================
+      logical cycle_diurne
+      real  gmtime
+      real rlat(klon), rlon(klon), rmu0(klon), fract(klon)
+c================================================================
+      integer i
+      real gmtime1, gmtime2
+      real pi_local
+
+
+      real rmu0m(klon),rmu0a(klon)
+c
+
+      pi_local = 4.0 * ATAN(1.0)
+
+c================================================================
+c  Calcul de l'angle zenithal moyen sur la journee
+c================================================================
+
+      DO i=1,klon
+        fract(i)=1. 
+!  Calcule du flux moyen
+        IF (abs(rlat(i)).LE.28.75) THEN
+           rmu0m(i)=(210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
+        ELSEIF (abs(rlat(i)).LE.43.75) THEN
+          rmu0m(i)=(187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
+        ELSEIF (abs(rlat(i)).LE.71.25) THEN
+          rmu0m(i)=(162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
+        ELSE
+          rmu0m(i)=(172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
+        ENDIF
+      ENDDO
+
+c================================================================
+!  Avec ou sans cycle diurne
+c================================================================
+
+      IF (cycle_diurne) THEN
+
+!  On redecompose flux  au sommet suivant un cycle diurne idealise
+!  identique a toutes les latitudes.
+
+         DO i=1,klon
+           rmu0a(i)=2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
+           rmu0(i)=rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*
+     &      rlon(i)/360.))-rmu0a(i)/sqrt(2.)
+         ENDDO
+
+         DO i=1,klon
+           IF (rmu0(i).LE.0.) THEN
+              rmu0(i)=0.
+              fract(i)=0.
+           ELSE
+              fract(i)=1.
+           ENDIF
+         ENDDO
+
+!  Affichage de l'angel zenitale
+!     print*,'************************************'
+!     print*,'************************************'
+!     print*,'************************************'
+!     print*,'latitude=',rlat(i),'longitude=',rlon(i)
+!     print*,'rmu0m=',rmu0m(i)
+!     print*,'rmu0a=',rmu0a(i)
+!     print*,'rmu0=',rmu0(i)
+                                                               
+      ELSE
+
+        DO i=1,klon
+           fract(i)=0.5
+           rmu0(i)=rmu0m(i)*2.
+        ENDDO
+
+      ENDIF
+
+      RETURN
+      END subroutine zenang_an
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      subroutine writelim
+     s   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,
+     s    phy_fter,phy_foce,phy_flic,phy_fsic)
+c
+      use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root
+      use mod_grid_phy_lmdz, only : klon_glo
+      use mod_phys_lmdz_transfert_para, only : gather
+      implicit none
+!#include "dimensions.h"
+!#include "dimphy.h"
+#include "netcdf.inc"
+ 
+      integer,intent(in) :: klon
+      real,intent(in) :: phy_nat(klon,360)
+      real,intent(in) :: phy_alb(klon,360)
+      real,intent(in) :: phy_sst(klon,360)
+      real,intent(in) :: phy_bil(klon,360)
+      real,intent(in) :: phy_rug(klon,360)
+      real,intent(in) :: phy_ice(klon,360)
+      real,intent(in) :: phy_fter(klon,360)
+      real,intent(in) :: phy_foce(klon,360)
+      real,intent(in) :: phy_flic(klon,360)
+      real,intent(in) :: phy_fsic(klon,360)
+
+      real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:)
+                                    ! on the whole physics grid 
+      integer :: k
+      INTEGER ierr
+      INTEGER dimfirst(3)
+      INTEGER dimlast(3)
+c
+      INTEGER nid, ndim, ntim
+      INTEGER dims(2), debut(2), epais(2)
+      INTEGER id_tim
+      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
+      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
+ 
+      if (is_mpi_root.and.is_omp_root) then
+      
+        PRINT*, 'writelim: Ecriture du fichier limit'
+c
+        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
+c
+        ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,
+     .                       "Fichier conditions aux limites")
+!!        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
+        ierr = NF_DEF_DIM (nid, "points_physiques", klon_glo, ndim)
+        ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
+c
+        dims(1) = ndim
+        dims(2) = ntim
+c
+ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
+        ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
+        ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,
+     .                        "Jour dans l annee")
+ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
+        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
+        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,
+     .                        "Nature du sol (0,1,2,3)")
+ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
+        ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
+        ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,
+     .                        "Temperature superficielle de la mer")
+ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
+        ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
+        ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,
+     .                        "Reference flux de chaleur au sol")
+ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
+        ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
+        ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,
+     .                        "Albedo a la surface")
+ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
+        ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
+        ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,
+     .                        "Rugosite")
+
+        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
+        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
+        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
+c
+        ierr = NF_ENDDEF(nid)
+        if (ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error: failed to end define mode"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+c
+
+! write the 'times'
+        do k=1,360
+#ifdef NC_DOUBLE
+          ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
+#else
+          ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
+#endif
+         if (ierr.ne.NF_NOERR) then
+           write(*,*) "writelim error with temps(k),k=",k
+           write(*,*) NF_STRERROR(ierr)
+         endif
+        enddo
+
+      endif ! of if (is_mpi_root.and.is_omp_root)
+
+! write the fields, after having collected them on master
+
+      call gather(phy_nat,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_nat"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_sst,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_sst"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_bil,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_bil"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_alb,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_alb"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_rug,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_rug"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_fter,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_fter"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_foce,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_foce"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_fsic,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_fsic"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+      call gather(phy_flic,phy_glo)
+      if (is_mpi_root.and.is_omp_root) then
+#ifdef NC_DOUBLE
+        ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo)
+#else
+        ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo)
+#endif
+        if(ierr.ne.NF_NOERR) then
+          write(*,*) "writelim error with phy_flic"
+          write(*,*) NF_STRERROR(ierr)
+        endif
+      endif
+
+!  close file:
+      if (is_mpi_root.and.is_omp_root) then
+        ierr = NF_CLOSE(nid)
+      endif
+
+      end subroutine writelim 
+
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+
+      SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
+      use dimphy
+      IMPLICIT none
+c
+      INTEGER nlon, type_profil, i, k, j
+      REAL :: rlatd(nlon), phy_sst(nlon, 360) 
+      INTEGER imn, imx, amn, amx, kmn, kmx
+      INTEGER p, pplus, nlat_max
+      parameter (nlat_max=72)
+      REAL x_anom_sst(nlat_max)
+c
+      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
+      write(*,*)" profil_sst: type_profil=",type_profil
+      do i=1,360
+c      phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
+
+c Rajout fbrlmd
+
+      if(type_profil.EQ.1)then
+c     Méthode 1 "Control" faible plateau à l'Equateur
+      do j=1,klon
+       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**2)
+c        PI/3=1.047197551
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+      if(type_profil.EQ.2)then
+c     Méthode 2 "Flat" fort plateau à l'Equateur
+      do j=1,klon
+       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+
+      if (type_profil.EQ.3) then
+c     Méthode 3 "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.4) then
+c     Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.5) then
+c     Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.+2.
+        endif
+
+      enddo
+      endif
+
+      if(type_profil.EQ.6)then
+c     Méthode 6 "cst" valeur constante de SST
+      do j=1,klon
+       phy_sst(j,i)=288.
+      enddo
+      endif
+
+
+        if(type_profil.EQ.7)then
+c     Méthode 7 "cst" valeur constante de SST +2
+      do j=1,klon
+       phy_sst(j,i)=288.+2.
+      enddo
+      endif
+
+        p=0
+        if(type_profil.EQ.8)then
+c     Méthode 8 profil anomalies SST du modèle couplé AR4
+       do j=1,klon
+         if (rlatd(j).EQ.rlatd(j-1)) then
+       phy_sst(j,i)=273.+x_anom_sst(pplus)
+     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
+          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+            phy_sst(j,i)=273.+x_anom_sst(pplus)
+          endif
+        else
+          p=p+1
+          pplus=73-p
+        phy_sst(j,i)=273.+x_anom_sst(pplus)
+     &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
+          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+            phy_sst(j,i)=273.+x_anom_sst(pplus)
+          endif
+          write (*,*) rlatd(j),x_anom_sst(pplus),phy_sst(j,i)
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.9) then
+c     Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.-2.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.-2.
+        endif
+      enddo
+      endif
+
+
+      if (type_profil.EQ.10) then
+c     Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.+4.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.11) then
+c     Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.12) then
+c     Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
+      do j=1,klon
+        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.+4.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.13) then
+c     Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
+      do j=1,klon
+        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+      if (type_profil.EQ.14) then
+c     Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
+      do j=1,klon
+        phy_sst(j,i)=273.+2.+0.5*29.*(2-sin(1.5*rlatd(j))**2
+     &                 -sin(1.5*rlatd(j))**4)
+      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
+         phy_sst(j,i)=273.
+        endif
+      enddo
+      endif
+
+      enddo
+
+cIM beg : verif profil SST: phy_sst
+       amn=MIN(phy_sst(1,1),1000.)
+       amx=MAX(phy_sst(1,1),-1000.)
+       imn=1 ; kmn=1 ; imx=1 ; kmx=1
+       DO k=1, 360
+       DO i=2, nlon
+        IF(phy_sst(i,k).LT.amn) THEN
+         amn=phy_sst(i,k)
+         imn=i
+         kmn=k
+        ENDIF
+        IF(phy_sst(i,k).GT.amx) THEN
+         amx=phy_sst(i,k)
+         imx=i
+         kmx=k
+        ENDIF
+       ENDDO
+       ENDDO
+c
+       PRINT*,'profil_sst: imn, kmn, phy_sst(imn,kmn) ',imn,kmn,amn
+       PRINT*,'profil_sst: imx, kmx, phy_sst(imx,kmx) ',imx,kmx,amx
+cIM end : verif profil SST: phy_sst
+
+       return 
+       end subroutine profil_sst
+     
+      end module phyaqua_mod
Index: /LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1850)
+++ /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1851)
@@ -45,4 +45,5 @@
       use conf_phys_m, only: conf_phys
       use radlwsw_m, only: radlwsw
+      use phyaqua_mod, only: zenang_an
       USE control_mod
 #ifdef REPROBUS
