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