Index: LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90
===================================================================
--- LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5117)
+++ LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90	(revision 5118)
@@ -1,10 +1,9 @@
-
 ! $Id$
 
-SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
+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 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
@@ -15,8 +14,9 @@
   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 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 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
 
   !   Author:    Frederic Hourdin      original: 15/01/93
@@ -33,40 +33,39 @@
   include "comgeom.h"
   include "academic.h"
-  include "iniprint.h"
 
   !   Arguments:
   !   ----------
 
-  REAL,INTENT(OUT) :: time_0
+  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
+  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 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 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
+  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 phi_pv, dphi_pv, gam_pv, tetanoise   ! Constantes pour polar vortex
 
   REAL zz
@@ -74,18 +73,18 @@
 
   REAL zdtvr, tnat, alpha_ideal
-  LOGICAL,PARAMETER :: tnat1=.TRUE.
-  
-  CHARACTER(LEN=*),parameter :: modname="iniacademic"
-  CHARACTER(LEN=80) :: abort_message
+  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)
+    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
@@ -95,22 +94,22 @@
   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.
+  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
+    ! specific initializations for the shallow water case
+    kappa = 1
   ENDIF
 
@@ -126,237 +125,237 @@
   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)
+    !------------------------------------------------------------------
+    ! 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
-     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
-     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
+      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
-
-        ! maintain periodicity in longitude
-        do l=1,llm
-           do ij=1,ip1jmp1,iip1
-              teta(ij+iim,l)=teta(ij,l)
-           enddo
+      enddo
+
+      ! maintain periodicity in longitude
+      do l = 1, llm
+        do ij = 1, ip1jmp1, iip1
+          teta(ij + iim, l) = teta(ij, l)
         enddo
-
-     ENDIF ! of IF (.NOT. read_start)
+      enddo
+
+    ENDIF ! of IF (.NOT. read_start)
   ENDIF academic_case
 
