source: LMDZ6/branches/Amaury_dev/libf/phylmd/phyaqua_mod.F90 @ 5411

! $Id: phyaqua_mod.F90 5158 2024-08-02 12:12:03Z evignon $

MODULE phyaqua_mod
  ! Routines complementaires pour la physique planetaire.
  USE lmdz_abort_physic, ONLY: abort_physic
  IMPLICIT NONE

CONTAINS

  SUBROUTINE iniaqua(nlon, year_len, 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 dimphy, ONLY: klon
    USE lmdz_geometry, ONLY: latitude
    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 time_phylmdz_mod, ONLY: day_ref, ndays, pdtphys, &
            day_ini, day_end
    USE indice_sol_mod
    USE lmdz_physical_constants, ONLY: pi
    !    USE ioipsl
    USE lmdz_phys_para, ONLY: is_master
    USE lmdz_phys_transfert_para, ONLY: bcast
    USE lmdz_grid_phy
    USE lmdz_ioipsl_getin_p, ONLY: getin_p
    USE phys_cal_mod, ONLY: calend, year_len_phy => year_len
    USE lmdz_clesphys
    USE lmdz_yomcst

    IMPLICIT NONE

    include "dimsoil.h"

    INTEGER, INTENT (IN) :: nlon, year_len, iflag_phys
    ! IM ajout latfi, lonfi
    !    REAL, INTENT (IN) :: lonfi(nlon), latfi(nlon)

    INTEGER type_profil, type_aqua

    ! Ajouts initialisation des surfaces
    REAL :: run_off_lic_0(nlon)
    REAL :: qsolsrf(nlon, nbsrf), snsrf(nlon, nbsrf)
    REAL :: tsoil(nlon, nsoilmx, nbsrf)
    REAL :: tslab(nlon), seaice(nlon)
    REAL fder(nlon)



    ! Arguments :
    ! -----------

    ! integer radpas
    INTEGER it, unit, i, k, itap

    REAL rugos, albedo
    REAL tsurf
    REAL time, timestep, day, day0
    REAL qsol_f
    REAL rugsrel(nlon)
    LOGICAL alb_ocean

    CHARACTER *80 ans, file_forctl, file_fordat, file_start
    CHARACTER *100 file, var
    CHARACTER *2 cnbl

    REAL phy_nat(nlon, year_len)
    REAL phy_alb(nlon, year_len)
    REAL phy_sst(nlon, year_len)
    REAL phy_bil(nlon, year_len)
    REAL phy_rug(nlon, year_len)
    REAL phy_ice(nlon, year_len)
    REAL phy_fter(nlon, year_len)
    REAL phy_foce(nlon, year_len)
    REAL phy_fsic(nlon, year_len)
    REAL phy_flic(nlon, year_len)

    INTEGER, SAVE :: read_climoz = 0 ! read ozone climatology
    !$OMP THREADPRIVATE(read_climoz)

    ! -------------------------------------------------------------------------
    ! 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)


    ! -----------------------------------------------------------------------
    ! dynamial tendencies :
    ! ---------------------

    INTEGER l, ierr, aslun

    REAL paire

    ! Local
    CHARACTER (LEN = 20) :: modname = 'phyaqua'
    CHARACTER (LEN = 80) :: abort_message


    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    ! INITIALISATIONS
    ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! -----------------------------------------------------------------------
    ! Initialisations  des constantes
    ! -------------------------------

    !IF (calend .EQ. "earth_360d") Then
    year_len_phy = year_len
    !END IF

    IF (year_len/=360) THEN
      write (*, *) year_len
      CALL abort_physic("iniaqua", 'iniaqua: 360 day calendar is required !', 1)
    endif

    type_aqua = iflag_phys / 100
    type_profil = iflag_phys - type_aqua * 100
    PRINT *, 'iniaqua:type_aqua, type_profil', type_aqua, type_profil

    IF (klon/=nlon) THEN
      WRITE (*, *) 'iniaqua: klon=', klon, ' nlon=', nlon
      abort_message = 'probleme de dimensions dans iniaqua'
      CALL abort_physic(modname, abort_message, 1)
    END IF
    CALL phys_state_var_init(read_climoz)

    read_climoz = 0
    day0 = 217.
    day = day0
    it = 0
    time = 0.

    ! -----------------------------------------------------------------------
    ! initialisations de la physique
    ! -----------------------------------------------------------------------

    day_ini = day_ref
    day_end = day_ini + ndays

    nbapp_rad = 24
    CALL getin_p('nbapp_rad', nbapp_rad)

    ! ---------------------------------------------------------------------
    ! Creation des conditions aux limites:
    ! ------------------------------------
    ! Initialisations des constantes
    ! Ajouter les manquants dans planete.def... (albedo etc)
    co2_ppm = 348.
    CALL getin_p('co2_ppm', co2_ppm)

    solaire = 1365.
    CALL getin_p('solaire', solaire)

    ! CALL getin('albedo',albedo) ! albedo is set below, depending on
    ! type_aqua
    alb_ocean = .TRUE.
    CALL getin_p('alb_ocean', alb_ocean)

    WRITE (*, *) 'iniaqua: co2_ppm=', co2_ppm
    WRITE (*, *) 'iniaqua: solaire=', solaire
    WRITE (*, *) 'iniaqua: alb_ocean=', alb_ocean

    radsol = 0.
    qsol_f = 10.

    ! Conditions aux limites:
    ! -----------------------

    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.
    END IF

    CALL getin_p('rugos', rugos)

    WRITE (*, *) 'iniaqua: rugos=', rugos
    zmasq(:) = pctsrf(:, is_ter)

    ! 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)
    ! END IF !alb_ocean

    DO i = 1, year_len
      ! IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
      ! IM 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 (?)
    END DO
    ! IM calcul profil sst
    CALL profil_sst(nlon, latitude, type_profil, phy_sst)

    IF (grid_type==unstructured) THEN
      CALL writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
              phy_fter, phy_foce, phy_flic, phy_fsic)
    ELSE

      CALL writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, phy_ice, &
              phy_fter, phy_foce, phy_flic, phy_fsic)
    ENDIF

    ! ---------------------------------------------------------------------
    ! Ecriture de l'etat initial:
    ! ---------------------------


    ! Ecriture etat initial physique

    timestep = pdtphys
    radpas = nint(rday / timestep / float(nbapp_rad))

    DO i = 1, longcles
      clesphy0(i) = 0.
    END DO
    clesphy0(1) = float(iflag_con)
    clesphy0(2) = float(nbapp_rad)
    ! 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


    ! =======================================================================
    ! Profils initiaux
    ! =======================================================================

    ! 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)
    END DO

    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.
    fevap = 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
    z0m(:, :) = rugos ! couverture de neige des sous surface
    z0h = z0m

    CALL pbl_surface_init(fder, snsrf, qsolsrf, tsoil)

    PRINT *, 'iniaqua: before phyredem'

    pbl_tke(:, :, :) = 1.e-8
    falb1 = albedo
    falb2 = albedo
    zmax0 = 0.
    f0 = 0.
    sig1 = 0.
    w01 = 0.
    wake_deltat = 0.
    wake_deltaq = 0.
    wake_s = 0.
    wake_dens = 0.
    wake_cstar = 0.
    wake_pe = 0.
    wake_fip = 0.
    fm_therm = 0.
    entr_therm = 0.
    detr_therm = 0.
    ale_bl = 0.
    ale_bl_trig = 0.
    alp_bl = 0.
    treedrg(:, :, :) = 0.

    u10m = 0.
    v10m = 0.

    ql_ancien = 0.
    qs_ancien = 0.
    qbs_ancien = 0.
    u_ancien = 0.
    v_ancien = 0.
    prw_ancien = 0.
    prlw_ancien = 0.
    prsw_ancien = 0.
    prbsw_ancien = 0.

    ale_wake = 0.
    ale_bl_stat = 0.


    !ym error : the sub surface dimension is the third not second : forgotten for iniaqua
    !    falb_dir(:,is_ter,:)=0.08; falb_dir(:,is_lic,:)=0.6
    !    falb_dir(:,is_oce,:)=0.5;  falb_dir(:,is_sic,:)=0.6
    falb_dir(:, :, is_ter) = 0.08; falb_dir(:, :, is_lic) = 0.6
    falb_dir(:, :, is_oce) = 0.5;  falb_dir(:, :, is_sic) = 0.6

    !ym falb_dif has been forgotten, initialize with defaukt value found in phyetat0 or 0 ?
    !ym probably the uninitialized value was 0 for standard (regular grid) case
    falb_dif(:, :, :) = 0

    CALL phyredem('startphy.nc')

    PRINT *, 'iniaqua: after phyredem'
    CALL phys_state_var_end

  END SUBROUTINE iniaqua


  ! ====================================================================
  ! ====================================================================
  SUBROUTINE zenang_an(cycle_diurne, gmtime, rlat, rlon, rmu0, fract)
    USE dimphy
    USE lmdz_yomcst

    IMPLICIT NONE
    ! ====================================================================
    ! =============================================================
    ! CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
    ! Auteur : A. Campoy et F. Hourdin
    ! Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
    ! et l'ensoleillement moyen entre gmtime1 et gmtime2
    ! connaissant la declinaison, la latitude et la longitude.

    ! Dans cette version particuliere, on calcule le rayonnement
    ! moyen sur l'année à chaque latitude.
    ! angle zenithal calculé pour obtenir un
    ! Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
    ! en moyenne annuelle.
    ! Spécifique de la terre. Utilisé pour les aqua planetes.

    ! Rque   : Different de la routine angle en ce sens que zenang
    ! fournit des moyennes de pmu0 et non des valeurs
    ! instantanees, du coup frac prend toutes les valeurs
    ! entre 0 et 1.
    ! Date   : premiere version le 13 decembre 1994
    ! revu pour  GCM  le 30 septembre 1996
    ! ===============================================================
    ! longi----INPUT : la longitude vraie de la terre dans son plan
    ! solaire a partir de l'equinoxe de printemps (degre)
    ! gmtime---INPUT : temps universel en fraction de jour
    ! pdtrad---INPUT : pas de temps du rayonnement (secondes)
    ! lat------INPUT : latitude en degres
    ! long-----INPUT : longitude en degres
    ! pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
    ! frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
    ! ================================================================
    LOGICAL cycle_diurne
    REAL gmtime
    REAL rlat(klon), rlon(klon), rmu0(klon), fract(klon)
    ! ================================================================
    INTEGER i
    REAL gmtime1, gmtime2
    REAL pi_local

    REAL rmu0m(klon), rmu0a(klon)

    pi_local = 4.0 * atan(1.0)

    ! ================================================================
    ! Calcul de l'angle zenithal moyen sur la journee
    ! ================================================================

    DO i = 1, klon
      fract(i) = 1.
      ! Calcule du flux moyen
      IF (abs(rlat(i))<=28.75) THEN
        rmu0m(i) = (210.1924 + 206.6059 * cos(0.0174533 * rlat(i))**2) / 1365.
      ELSE IF (abs(rlat(i))<=43.75) THEN
        rmu0m(i) = (187.4562 + 236.1853 * cos(0.0174533 * rlat(i))**2) / 1365.
      ELSE IF (abs(rlat(i))<=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.
      END IF
    END DO

    ! ================================================================
    ! Avec ou sans cycle diurne
    ! ================================================================

    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.)
      END DO

      DO i = 1, klon
        IF (rmu0(i)<=0.) THEN
          rmu0(i) = 0.
          fract(i) = 0.
        ELSE
          fract(i) = 1.
        END IF
      END DO

      ! 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.
      END DO

    END IF

  END SUBROUTINE zenang_an

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE writelim_unstruct(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)

    USE lmdz_phys_para, ONLY: is_omp_master, klon_mpi
    USE lmdz_phys_transfert_para, ONLY: gather_omp
    USE lmdz_xios
    IMPLICIT NONE

    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_mpi(klon_mpi, 360) ! temporary variable, to store phy_***(:)
    ! on the whole physics grid

    IF (using_xios) THEN
      PRINT *, 'writelim: Ecriture du fichier limit'

      CALL gather_omp(phy_foce, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('foce_limout', phy_mpi)

      CALL gather_omp(phy_fsic, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('fsic_limout', phy_mpi)

      CALL gather_omp(phy_fter, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('fter_limout', phy_mpi)

      CALL gather_omp(phy_flic, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('flic_limout', phy_mpi)

      CALL gather_omp(phy_sst, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('sst_limout', phy_mpi)

      CALL gather_omp(phy_bil, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('bils_limout', phy_mpi)

      CALL gather_omp(phy_alb, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('alb_limout', phy_mpi)

      CALL gather_omp(phy_rug, phy_mpi)
      IF (is_omp_master) CALL xios_send_field('rug_limout', phy_mpi)
    ENDIF
  END SUBROUTINE writelim_unstruct


  SUBROUTINE writelim(klon, phy_nat, phy_alb, phy_sst, phy_bil, phy_rug, &
          phy_ice, phy_fter, phy_foce, phy_flic, phy_fsic)

    USE lmdz_phys_para, ONLY: is_master
    USE lmdz_grid_phy, ONLY: klon_glo
    USE lmdz_phys_transfert_para, ONLY: gather
    USE phys_cal_mod, ONLY: year_len
    USE netcdf, ONLY: nf90_def_var, nf90_put_var, nf90_get_var, nf90_strerror, nf90_close, &
            nf90_enddef, nf90_put_att, nf90_unlimited, nf90_noerr, nf90_global, nf90_clobber, &
            nf90_64bit_offset, nf90_def_dim, nf90_create
    USE lmdz_cppkeys_wrapper, ONLY: nf90_format
    IMPLICIT NONE

    INTEGER, INTENT (IN) :: klon
    REAL, INTENT (IN) :: phy_nat(klon, year_len)
    REAL, INTENT (IN) :: phy_alb(klon, year_len)
    REAL, INTENT (IN) :: phy_sst(klon, year_len)
    REAL, INTENT (IN) :: phy_bil(klon, year_len)
    REAL, INTENT (IN) :: phy_rug(klon, year_len)
    REAL, INTENT (IN) :: phy_ice(klon, year_len)
    REAL, INTENT (IN) :: phy_fter(klon, year_len)
    REAL, INTENT (IN) :: phy_foce(klon, year_len)
    REAL, INTENT (IN) :: phy_flic(klon, year_len)
    REAL, INTENT (IN) :: phy_fsic(klon, year_len)

    REAL :: phy_glo(klon_glo, year_len) ! temporary variable, to store phy_***(:)
    ! on the whole physics grid
    INTEGER :: k
    INTEGER ierr
    INTEGER dimfirst(3)
    INTEGER dimlast(3)

    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_master) THEN

      PRINT *, 'writelim: Ecriture du fichier limit'

      ierr = nf90_create('limit.nc', IOR(nf90_clobber, nf90_64bit_offset), nid)

      ierr = nf90_put_att(nid, nf90_global, 'title', &
              'Fichier conditions aux limites')
      !        ierr = nf90_def_dim (nid, "points_physiques", klon, ndim)
      ierr = nf90_def_dim(nid, 'points_physiques', klon_glo, ndim)
      ierr = nf90_def_dim(nid, 'time', nf90_unlimited, ntim)

      dims(1) = ndim
      dims(2) = ntim

      ierr = nf90_def_var(nid, 'TEMPS', nf90_format, [ntim], id_tim)
      ierr = nf90_put_att(nid, id_tim, 'title', 'Jour dans l annee')

      ierr = nf90_def_var(nid, 'NAT', nf90_format, dims, id_nat)
      ierr = nf90_put_att(nid, id_nat, 'title', &
              'Nature du sol (0,1,2,3)')

      ierr = nf90_def_var(nid, 'SST', nf90_format, dims, id_sst)
      ierr = nf90_put_att(nid, id_sst, 'title', &
              'Temperature superficielle de la mer')

      ierr = nf90_def_var(nid, 'BILS', nf90_format, dims, id_bils)
      ierr = nf90_put_att(nid, id_bils, 'title', &
              'Reference flux de chaleur au sol')

      ierr = nf90_def_var(nid, 'ALB', nf90_format, dims, id_alb)
      ierr = nf90_put_att(nid, id_alb, 'title', 'Albedo a la surface')

      ierr = nf90_def_var(nid, 'RUG', nf90_format, dims, id_rug)
      ierr = nf90_put_att(nid, id_rug, 'title', 'Rugosite')

      ierr = nf90_def_var(nid, 'FTER', nf90_format, dims, id_fter)
      ierr = nf90_put_att(nid, id_fter, 'title', 'Frac. Land')
      ierr = nf90_def_var(nid, 'FOCE', nf90_format, dims, id_foce)
      ierr = nf90_put_att(nid, id_foce, 'title', 'Frac. Ocean')
      ierr = nf90_def_var(nid, 'FSIC', nf90_format, dims, id_fsic)
      ierr = nf90_put_att(nid, id_fsic, 'title', 'Frac. Sea Ice')
      ierr = nf90_def_var(nid, 'FLIC', nf90_format, dims, id_flic)
      ierr = nf90_put_att(nid, id_flic, 'title', 'Frac. Land Ice')

      ierr = nf90_enddef(nid)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error: failed to end define mode'
        WRITE (*, *) nf90_strerror(ierr)
      END IF


      ! write the 'times'
      DO k = 1, year_len
        ierr = nf90_put_var(nid, id_tim, k, [k])
        IF (ierr/=nf90_noerr) THEN
          WRITE (*, *) 'writelim error with temps(k),k=', k
          WRITE (*, *) nf90_strerror(ierr)
        END IF
      END DO

    END IF ! of if (is_master)

    ! write the fields, after having collected them on master

    CALL gather(phy_nat, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_nat, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_nat'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_sst, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_sst, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_sst'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_bil, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_bils, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_bil'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_alb, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_alb, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_alb'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_rug, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_rug, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_rug'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_fter, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_fter, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_fter'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_foce, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_foce, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_foce'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_fsic, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_fsic, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_fsic'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    CALL gather(phy_flic, phy_glo)
    IF (is_master) THEN
      ierr = nf90_put_var(nid, id_flic, phy_glo)
      IF (ierr/=nf90_noerr) THEN
        WRITE (*, *) 'writelim error with phy_flic'
        WRITE (*, *) nf90_strerror(ierr)
      END IF
    END IF

    ! close file:
    IF (is_master) THEN
      ierr = nf90_close(nid)
    END IF

  END SUBROUTINE writelim

  ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
    USE dimphy
    USE phys_cal_mod, ONLY: year_len
    IMPLICIT NONE

    INTEGER nlon, type_profil, i, k, j
    REAL :: rlatd(nlon), phy_sst(nlon, year_len)
    INTEGER imn, imx, amn, amx, kmn, kmx
    INTEGER p, pplus, nlat_max
    PARAMETER (nlat_max = 72)
    REAL x_anom_sst(nlat_max)
    CHARACTER (LEN = 20) :: modname = 'profil_sst'
    CHARACTER (LEN = 80) :: abort_message

    IF (klon/=nlon) THEN
      abort_message = 'probleme de dimensions dans profil_sst'
      CALL abort_physic(modname, abort_message, 1)
    ENDIF
    WRITE (*, *) ' profil_sst: type_profil=', type_profil
    DO i = 1, year_len
      ! phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2

      ! Rajout fbrlmd

      IF (type_profil==1) THEN
        ! 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)
          ! PI/3=1.047197551
          IF ((rlatd(j)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF
      IF (type_profil==2) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==3) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==4) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==5) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273. + 2.
          END IF

        END DO
      END IF

      IF (type_profil==6) THEN
        ! Méthode 6 "cst" valeur constante de SST
        DO j = 1, klon
          phy_sst(j, i) = 288.
        END DO
      END IF

      IF (type_profil==7) THEN
        ! Méthode 7 "cst" valeur constante de SST +2
        DO j = 1, klon
          phy_sst(j, i) = 288. + 2.
        END DO
      END IF

      p = 0
      IF (type_profil==8) THEN
        ! Méthode 8 profil anomalies SST du modèle couplé AR4
        DO j = 1, klon
          IF (rlatd(j)==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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
              phy_sst(j, i) = 273. + x_anom_sst(pplus)
            END IF
          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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
              phy_sst(j, i) = 273. + x_anom_sst(pplus)
            END IF
            WRITE (*, *) rlatd(j), x_anom_sst(pplus), phy_sst(j, i)
          END IF
        END DO
      END IF

      IF (type_profil==9) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273. - 2.
          END IF
        END DO
      END IF

      IF (type_profil==10) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273. + 4.
          END IF
        END DO
      END IF

      IF (type_profil==11) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==12) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273. + 4.
          END IF
        END DO
      END IF

      IF (type_profil==13) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==14) THEN
        ! 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)>1.0471975) .OR. (rlatd(j)<-1.0471975)) THEN
            phy_sst(j, i) = 273.
          END IF
        END DO
      END IF

      IF (type_profil==20) THEN
        PRINT*, 'Profile SST 20'
        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K

        DO j = 1, klon
          phy_sst(j, i) = 248. + 55. * (1 - sin(rlatd(j))**2)
        enddo
      endif

      IF (type_profil==21) THEN
        PRINT*, 'Profile SST 21'
        !     Méthode 13 "Qmax2K" plateau réel é|  l'Equateur augmenté +2K
        DO j = 1, klon
          phy_sst(j, i) = 252. + 55. * (1 - sin(rlatd(j))**2)
        enddo
      endif

    END DO

    ! IM 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, year_len
      DO i = 2, nlon
        IF (phy_sst(i, k)<amn) THEN
          amn = phy_sst(i, k)
          imn = i
          kmn = k
        END IF
        IF (phy_sst(i, k)>amx) THEN
          amx = phy_sst(i, k)
          imx = i
          kmx = k
        END IF
      END DO
    END DO

    PRINT *, 'profil_sst: imn, kmn, phy_sst(imn,kmn) ', imn, kmn, amn
    PRINT *, 'profil_sst: imx, kmx, phy_sst(imx,kmx) ', imx, kmx, amx
    ! IM end : verif profil SST: phy_sst

  END SUBROUTINE profil_sst

END MODULE phyaqua_mod