! $Id: iniacademic.F90 5134 2024-07-26 15:56:37Z abarral $

SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)

  USE lmdz_filtreg, ONLY: inifilr
  USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
  USE control_mod, ONLY: day_step, planet_type
  USE exner_hyb_m, ONLY: exner_hyb
  USE exner_milieu_m, ONLY: exner_milieu
  USE IOIPSL, ONLY: getin
  USE lmdz_write_field
  USE comconst_mod, ONLY: cpp, kappa, g, daysec, dtvr, pi, im, jm
  USE logic_mod, ONLY: iflag_phys, read_start
  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
  USE lmdz_readTracFiles, ONLY: addPhase
  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var
  USE lmdz_ran1, ONLY: ran1
  USE lmdz_iniprint, ONLY: lunout, prt_level
  USE lmdz_academic, ONLY: tetarappel, knewt_t, kfrict, knewt_g, clat4

  !   Author:    Frederic Hourdin      original: 15/01/93
  ! The forcing defined here is from Held and Suarez, 1994, Bulletin
  ! of the American Meteorological Society, 75, 1825.

  IMPLICIT NONE

  !   Declararations:
  !   ---------------

  INCLUDE "dimensions.h"
  INCLUDE "paramet.h"
  INCLUDE "comgeom.h"

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

  REAL, INTENT(OUT) :: time_0

  !   fields
  REAL, INTENT(OUT) :: vcov(ip1jm, llm) ! meridional covariant wind
  REAL, INTENT(OUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
  REAL, INTENT(OUT) :: teta(ip1jmp1, llm) ! potential temperature (K)
  REAL, INTENT(OUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers (.../kg_of_air)
  REAL, INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
  REAL, INTENT(OUT) :: masse(ip1jmp1, llm) ! air mass in grid cell (kg)
  REAL, INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential

  !   Local:
  !   ------

  REAL p (ip1jmp1, llmp1)               ! pression aux interfac.des couches
  REAL pks(ip1jmp1)                      ! exner au  sol
  REAL pk(ip1jmp1, llm)                   ! exner au milieu des couches
  REAL phi(ip1jmp1, llm)                  ! geopotentiel
  REAL ddsin, zsig, tetapv, w_pv  ! variables auxiliaires
  REAL tetastrat ! potential temperature in the stratosphere, in K
  REAL tetajl(jjp1, llm)
  INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent

  INTEGER :: nid_relief, varid, ierr
  REAL, DIMENSION(iip1, jjp1) :: relief

  REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T
  REAL k_f, k_c_a, k_c_s         ! Constantes de rappel
  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
  LOGICAL ok_pv                ! Polar Vortex
  REAL phi_pv, dphi_pv, gam_pv, tetanoise   ! Constantes pour polar vortex

  REAL zz
  INTEGER idum

  REAL zdtvr, tnat, alpha_ideal
  LOGICAL, PARAMETER :: tnat1 = .TRUE.

  CHARACTER(LEN = *), parameter :: modname = "iniacademic"
  CHARACTER(LEN = 80) :: abort_message

  ! Sanity check: verify that options selected by user are not incompatible
  IF ((iflag_phys==1).AND. .NOT. read_start) THEN
    WRITE(lunout, *) trim(modname), " error: if read_start is set to ", &
            " false then iflag_phys should not be 1"
    WRITE(lunout, *) "You most likely want an aquaplanet initialisation", &
            " (iflag_phys >= 100)"
    CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1)
  ENDIF

  !-----------------------------------------------------------------------
  ! 1. Initializations for Earth-like case
  ! --------------------------------------

  ! initialize planet radius, rotation rate,...
  CALL conf_planete

  time_0 = 0.
  day_ref = 1
  annee_ref = 0

  im = iim
  jm = jjm
  day_ini = 1
  dtvr = daysec / REAL(day_step)
  zdtvr = dtvr
  etot0 = 0.
  ptot0 = 0.
  ztot0 = 0.
  stot0 = 0.
  ang0 = 0.

  IF (llm == 1) THEN
    ! specific initializations for the shallow water case
    kappa = 1
  ENDIF

  CALL iniconst
  CALL inigeom
  CALL inifilr


  !------------------------------------------------------------------
  ! Initialize pressure and mass field if read_start=.FALSE.
  !------------------------------------------------------------------

  IF (.NOT. read_start) THEN

    !------------------------------------------------------------------
    ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
    ! du modele.
    ! On suppose que le fichier relief_in.nc est stoké sur une grille
    ! iim*jjp1
    ! Facile a créer à partir de la commande
    ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
    !------------------------------------------------------------------

    relief = 0.
    ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief)
    IF (ierr==nf90_noerr) THEN
      ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid)
      IF (ierr==nf90_noerr) THEN
        ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1))
        relief(iip1, :) = relief(1, :)
      else
        CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1)
      endif
    endif
    ierr = nf90_close (nid_relief)

    !------------------------------------------------------------------
    ! Initialisation du geopotentiel au sol et de la pression
    !------------------------------------------------------------------

    PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g
    do j = 1, jjp1
      do i = 1, iip1
        phis((j - 1) * iip1 + i) = g * relief(i, j)
      enddo
    enddo
    PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g

    ! ground geopotential
    !phis(:)=0.
    ps(:) = preff
    CALL pression (ip1jmp1, ap, bp, ps, p)

    IF (pressure_exner) THEN
      CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
    else
      CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
    endif
    CALL massdair(p, masse)
  ENDIF

  IF (llm == 1) THEN
    ! initialize fields for the shallow water case, if required
    IF (.NOT.read_start) THEN
      phis(:) = 0.
      q(:, :, :) = 0
      CALL sw_case_williamson91_6(vcov, ucov, teta, masse, ps)
    endif
  ENDIF

  academic_case: if (iflag_phys >= 2) THEN
    ! initializations

    ! 1. local parameters
    ! by convention, winter is in the southern hemisphere
    ! Geostrophic wind or no wind?
    ok_geost = .TRUE.
    CALL getin('ok_geost', ok_geost)
    ! Constants for Newtonian relaxation and friction
    k_f = 1.                !friction
    CALL getin('k_j', k_f)
    k_f = 1. / (daysec * k_f)
    k_c_s = 4.  !cooling surface
    CALL getin('k_c_s', k_c_s)
    k_c_s = 1. / (daysec * k_c_s)
    k_c_a = 40. !cooling free atm
    CALL getin('k_c_a', k_c_a)
    k_c_a = 1. / (daysec * k_c_a)
    ! Constants for Teta equilibrium profile
    teta0 = 315.     ! mean Teta (S.H. 315K)
    CALL getin('teta0', teta0)
    ttp = 200.       ! Tropopause temperature (S.H. 200K)
    CALL getin('ttp', ttp)
    eps = 0.         ! Deviation to N-S symmetry(~0-20K)
    CALL getin('eps', eps)
    delt_y = 60.     ! Merid Temp. Gradient (S.H. 60K)
    CALL getin('delt_y', delt_y)
    delt_z = 10.     ! Vertical Gradient (S.H. 10K)
    CALL getin('delt_z', delt_z)
    ! Polar vortex
    ok_pv = .FALSE.
    CALL getin('ok_pv', ok_pv)
    phi_pv = -50.            ! Latitude of edge of vortex
    CALL getin('phi_pv', phi_pv)
    phi_pv = phi_pv * pi / 180.
    dphi_pv = 5.             ! Width of the edge
    CALL getin('dphi_pv', dphi_pv)
    dphi_pv = dphi_pv * pi / 180.
    gam_pv = 4.              ! -dT/dz vortex (in K/km)
    CALL getin('gam_pv', gam_pv)
    tetanoise = 0.005
    CALL getin('tetanoise', tetanoise)

    ! 2. Initialize fields towards which to relax
    ! Friction
    knewt_g = k_c_a
    DO l = 1, llm
      zsig = presnivs(l) / preff
      knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3)
      kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3)
    ENDDO
    DO j = 1, jjp1
      clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4
    ENDDO

    ! Potential temperature
    DO l = 1, llm
      zsig = presnivs(l) / preff
      tetastrat = ttp * zsig**(-kappa)
      tetapv = tetastrat
      IF ((ok_pv).AND.(zsig<0.1)) THEN
        tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g)
      ENDIF
      DO j = 1, jjp1
        ! Troposphere
        ddsin = sin(rlatu(j))
        tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin &
                - delt_z * (1. - ddsin * ddsin) * log(zsig)
        IF (planet_type=="giant") THEN
          tetajl(j, l) = teta0 + (delt_y * &
                  ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2)   &
                  / ((rlatu(j) * 3.14159 * eps + 0.0001)**2))     &
                  - delt_z * log(zsig)
        endif
        ! Profil stratospherique isotherme (+vortex)
        w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2.
        tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv
        tetajl(j, l) = MAX(tetajl(j, l), tetastrat)
      ENDDO
    ENDDO

    !          CALL writefield('theta_eq',tetajl)

    do l = 1, llm
      do j = 1, jjp1
        do i = 1, iip1
          ij = (j - 1) * iip1 + i
          tetarappel(ij, l) = tetajl(j, l)
        enddo
      enddo
    enddo

    ! 3. Initialize fields (if necessary)
    IF (.NOT. read_start) THEN
      ! bulk initialization of temperature
      IF (iflag_phys>10000) THEN
        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
        teta(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp)
      ELSE
        teta(:, :) = tetarappel(:, :)
      ENDIF
      ! geopotential
      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)

      DO l = 1, llm
        PRINT*, 'presnivs,play,l', presnivs(l), (pk(1, l) / cpp)**(1. / kappa) * preff
        !pks(ij) = (cpp/preff) * ps(ij)
        !pk(ij,1) = .5*pks(ij)
        ! pk = cpp * (p/preff)^kappa
      ENDDO

      ! winds
      IF (ok_geost) THEN
        CALL ugeostr(phi, ucov)
      else
        ucov(:, :) = 0.
      endif
      vcov(:, :) = 0.

      ! bulk initialization of tracers
      IF (planet_type=="earth") THEN
        ! Earth: first two tracers will be water
        do iq = 1, nqtot
          q(:, :, iq) = 0.
          IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:, :, iq) = 1.e-10
          IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:, :, iq) = 1.e-15

          ! CRisi: init des isotopes
          ! distill de Rayleigh très simplifiée
          iName = tracers(iq)%iso_iName
          IF (niso <= 0 .OR. iName <= 0) CYCLE
          iPhase = tracers(iq)%iso_iPhase
          iqParent = tracers(iq)%iqParent
          IF(tracers(iq)%iso_iZone == 0) THEN
            IF (tnat1) THEN
              tnat = 1.0
              alpha_ideal = 1.0
              WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1'
            else
              IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
                      CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
            endif
            q(:, :, iq) = q(:, :, iqParent) * tnat * (q(:, :, iqParent) / 30.e-3)**(alpha_ideal - 1.)
          ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
            IF(tracers(iq)%iso_iZone == 1) THEN
              ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
              ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
              q(:, :, iq) = q(:, :, iqIsoPha(iName, iPhase))
            else !IF(tracers(iq)%iso_iZone == 1) THEN
              q(:, :, iq) = 0.
            endif !IF(tracers(iq)%iso_iZone == 1) THEN
          END IF !IF(tracers(iq)%iso_iZone == 0) THEN
        enddo
      else
        q(:, :, :) = 0
      endif ! of if (planet_type=="earth")

      CALL check_isotopes_seq(q, 1, ip1jmp1, 'iniacademic_loc')

      ! add random perturbation to temperature
      idum = -1
      zz = ran1(idum)
      idum = 0
      do l = 1, llm
        do ij = iip2, ip1jm
          teta(ij, l) = teta(ij, l) * (1. + tetanoise * ran1(idum))
        enddo
      enddo

      ! maintain periodicity in longitude
      do l = 1, llm
        do ij = 1, ip1jmp1, iip1
          teta(ij + iim, l) = teta(ij, l)
        enddo
      enddo

    ENDIF ! of IF (.NOT. read_start)
  ENDIF academic_case

END SUBROUTINE iniacademic