Ignore:
Timestamp:
Jul 24, 2024, 4:39:59 PM (2 months ago)
Author:
abarral
Message:

Replace iniprint.h by lmdz_iniprint.f90
(lint) along the way

File:
1 edited

Legend:

Unmodified
Added
Removed

  • LMDZ6/branches/Amaury_dev/libf/dyn3d/iniacademic.F90

    r5117 r5118  
    1 
    21! $Id$
    32
    4 SUBROUTINE iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
     3SUBROUTINE iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
    54
    65  USE lmdz_filtreg, ONLY: inifilr
    7   USE infotrac,    ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
    8   USE control_mod, ONLY: day_step,planet_type
     6  USE infotrac, ONLY: nqtot, niso, iqIsoPha, tracers, getKey, isoName
     7  USE control_mod, ONLY: day_step, planet_type
    98  USE exner_hyb_m, ONLY: exner_hyb
    109  USE exner_milieu_m, ONLY: exner_milieu
    … …  
    1514  USE comvert_mod, ONLY: ap, bp, preff, pa, presnivs, pressure_exner
    1615  USE temps_mod, ONLY: annee_ref, day_ini, day_ref
    17   USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
     16  USE ener_mod, ONLY: etot0, ptot0, ztot0, stot0, ang0
    1817  USE lmdz_readTracFiles, ONLY: addPhase
    19   USE netcdf, ONLY: nf90_nowrite,nf90_open,nf90_noerr,nf90_inq_varid,nf90_close,nf90_get_var
     18  USE netcdf, ONLY: nf90_nowrite, nf90_open, nf90_noerr, nf90_inq_varid, nf90_close, nf90_get_var
    2019  USE lmdz_ran1, ONLY: ran1
     20  USE lmdz_iniprint, ONLY: lunout, prt_level
    2121
    2222  !   Author:    Frederic Hourdin      original: 15/01/93
    … …  
    3333  include "comgeom.h"
    3434  include "academic.h"
    35   include "iniprint.h"
    3635
    3736  !   Arguments:
    3837  !   ----------
    3938
    40   REAL,INTENT(OUT) :: time_0
     39  REAL, INTENT(OUT) :: time_0
    4140
    4241  !   fields
    43   REAL,INTENT(OUT) :: vcov(ip1jm,llm) ! meridional covariant wind
    44   REAL,INTENT(OUT) :: ucov(ip1jmp1,llm) ! zonal covariant wind
    45   REAL,INTENT(OUT) :: teta(ip1jmp1,llm) ! potential temperature (K)
    46   REAL,INTENT(OUT) :: q(ip1jmp1,llm,nqtot) ! advected tracers (.../kg_of_air)
    47   REAL,INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
    48   REAL,INTENT(OUT) :: masse(ip1jmp1,llm) ! air mass in grid cell (kg)
    49   REAL,INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
     42  REAL, INTENT(OUT) :: vcov(ip1jm, llm) ! meridional covariant wind
     43  REAL, INTENT(OUT) :: ucov(ip1jmp1, llm) ! zonal covariant wind
     44  REAL, INTENT(OUT) :: teta(ip1jmp1, llm) ! potential temperature (K)
     45  REAL, INTENT(OUT) :: q(ip1jmp1, llm, nqtot) ! advected tracers (.../kg_of_air)
     46  REAL, INTENT(OUT) :: ps(ip1jmp1) ! surface pressure (Pa)
     47  REAL, INTENT(OUT) :: masse(ip1jmp1, llm) ! air mass in grid cell (kg)
     48  REAL, INTENT(OUT) :: phis(ip1jmp1) ! surface geopotential
    5049
    5150  !   Local:
    5251  !   ------
    5352
    54   REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
     53  REAL p (ip1jmp1, llmp1)               ! pression aux interfac.des couches
    5554  REAL pks(ip1jmp1)                      ! exner au  sol
    56   REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
    57   REAL phi(ip1jmp1,llm)                  ! geopotentiel
    58   REAL ddsin,zsig,tetapv,w_pv  ! variables auxiliaires
     55  REAL pk(ip1jmp1, llm)                   ! exner au milieu des couches
     56  REAL phi(ip1jmp1, llm)                  ! geopotentiel
     57  REAL ddsin, zsig, tetapv, w_pv  ! variables auxiliaires
    5958  REAL tetastrat ! potential temperature in the stratosphere, in K
    60   REAL tetajl(jjp1,llm)
    61   INTEGER i,j,l,lsup,ij, iq, iName, iPhase, iqParent
    62 
    63   INTEGER :: nid_relief,varid,ierr
    64   REAL, DIMENSION(iip1,jjp1) :: relief
    65 
    66   REAL teta0,ttp,delt_y,delt_z,eps ! Constantes pour profil de T
    67   REAL k_f,k_c_a,k_c_s         ! Constantes de rappel
     59  REAL tetajl(jjp1, llm)
     60  INTEGER i, j, l, lsup, ij, iq, iName, iPhase, iqParent
     61
     62  INTEGER :: nid_relief, varid, ierr
     63  REAL, DIMENSION(iip1, jjp1) :: relief
     64
     65  REAL teta0, ttp, delt_y, delt_z, eps ! Constantes pour profil de T
     66  REAL k_f, k_c_a, k_c_s         ! Constantes de rappel
    6867  LOGICAL ok_geost             ! Initialisation vent geost. ou nul
    6968  LOGICAL ok_pv                ! Polar Vortex
    70   REAL phi_pv,dphi_pv,gam_pv,tetanoise   ! Constantes pour polar vortex
     69  REAL phi_pv, dphi_pv, gam_pv, tetanoise   ! Constantes pour polar vortex
    7170
    7271  REAL zz
    … …  
    7473
    7574  REAL zdtvr, tnat, alpha_ideal
    76   LOGICAL,PARAMETER :: tnat1=.TRUE.
    77  
    78   CHARACTER(LEN=*),parameter :: modname="iniacademic"
    79   CHARACTER(LEN=80) :: abort_message
     75  LOGICAL, PARAMETER :: tnat1 = .TRUE.
     76
     77  CHARACTER(LEN = *), parameter :: modname = "iniacademic"
     78  CHARACTER(LEN = 80) :: abort_message
    8079
    8180  ! Sanity check: verify that options selected by user are not incompatible
    8281  IF ((iflag_phys==1).AND. .NOT. read_start) THEN
    83     WRITE(lunout,*) trim(modname)," error: if read_start is set to ", &
    84     " false then iflag_phys should not be 1"
    85     WRITE(lunout,*) "You most likely want an aquaplanet initialisation", &
    86     " (iflag_phys >= 100)"
    87     CALL abort_gcm(modname,"incompatible iflag_phys==1 and read_start==.FALSE.",1)
     82    WRITE(lunout, *) trim(modname), " error: if read_start is set to ", &
     83            " false then iflag_phys should not be 1"
     84    WRITE(lunout, *) "You most likely want an aquaplanet initialisation", &
     85            " (iflag_phys >= 100)"
     86    CALL abort_gcm(modname, "incompatible iflag_phys==1 and read_start==.FALSE.", 1)
    8887  ENDIF
    89  
     88
    9089  !-----------------------------------------------------------------------
    9190  ! 1. Initializations for Earth-like case
    … …  
    9594  CALL conf_planete
    9695
    97   time_0=0.
    98   day_ref=1
    99   annee_ref=0
    100 
    101   im         = iim
    102   jm         = jjm
    103   day_ini    = 1
    104   dtvr    = daysec/REAL(day_step)
    105   zdtvr=dtvr
    106   etot0      = 0.
    107   ptot0      = 0.
    108   ztot0      = 0.
    109   stot0      = 0.
    110   ang0       = 0.
     96  time_0 = 0.
     97  day_ref = 1
     98  annee_ref = 0
     99
     100  im = iim
     101  jm = jjm
     102  day_ini = 1
     103  dtvr = daysec / REAL(day_step)
     104  zdtvr = dtvr
     105  etot0 = 0.
     106  ptot0 = 0.
     107  ztot0 = 0.
     108  stot0 = 0.
     109  ang0 = 0.
    111110
    112111  IF (llm == 1) THEN
    113      ! specific initializations for the shallow water case
    114      kappa=1
     112    ! specific initializations for the shallow water case
     113    kappa = 1
    115114  ENDIF
    116115
    … …  
    126125  IF (.NOT. read_start) THEN
    127126
    128      !------------------------------------------------------------------
    129      ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
    130      ! du modele.
    131      ! On suppose que le fichier relief_in.nc est stoké sur une grille
    132      ! iim*jjp1
    133      ! Facile a créer à partir de la commande
    134      ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
    135      !------------------------------------------------------------------
    136 
    137      relief=0.
    138      ierr = nf90_open ('relief_in.nc', nf90_nowrite,nid_relief)
    139      IF (ierr==nf90_noerr) THEN
    140          ierr=nf90_inq_varid(nid_relief,'RELIEF',varid)
    141          IF (ierr==nf90_noerr) THEN
    142               ierr=nf90_get_var(nid_relief,varid,relief(1:iim,1:jjp1))
    143               relief(iip1,:)=relief(1,:)
    144          else
    145               CALL abort_gcm ('iniacademic','variable RELIEF pas la',1)
    146          endif
    147      endif
    148      ierr = nf90_close (nid_relief)
    149 
    150      !------------------------------------------------------------------
    151      ! Initialisation du geopotentiel au sol et de la pression
    152      !------------------------------------------------------------------
    153 
    154      PRINT*,'relief=',minval(relief),maxval(relief),'g=',g
    155      do j=1,jjp1
    156         do i=1,iip1
    157            phis((j-1)*iip1+i)=g*relief(i,j)
     127    !------------------------------------------------------------------
     128    ! Lecture eventuelle d'un fichier de relief interpollee sur la grille
     129    ! du modele.
     130    ! On suppose que le fichier relief_in.nc est stoké sur une grille
     131    ! iim*jjp1
     132    ! Facile a créer à partir de la commande
     133    ! cdo remapcon,fichier_output_phys.nc Relief.nc relief_in.nc
     134    !------------------------------------------------------------------
     135
     136    relief = 0.
     137    ierr = nf90_open ('relief_in.nc', nf90_nowrite, nid_relief)
     138    IF (ierr==nf90_noerr) THEN
     139      ierr = nf90_inq_varid(nid_relief, 'RELIEF', varid)
     140      IF (ierr==nf90_noerr) THEN
     141        ierr = nf90_get_var(nid_relief, varid, relief(1:iim, 1:jjp1))
     142        relief(iip1, :) = relief(1, :)
     143      else
     144        CALL abort_gcm ('iniacademic', 'variable RELIEF pas la', 1)
     145      endif
     146    endif
     147    ierr = nf90_close (nid_relief)
     148
     149    !------------------------------------------------------------------
     150    ! Initialisation du geopotentiel au sol et de la pression
     151    !------------------------------------------------------------------
     152
     153    PRINT*, 'relief=', minval(relief), maxval(relief), 'g=', g
     154    do j = 1, jjp1
     155      do i = 1, iip1
     156        phis((j - 1) * iip1 + i) = g * relief(i, j)
     157      enddo
     158    enddo
     159    PRINT*, 'phis=', minval(phis), maxval(phis), 'g=', g
     160
     161    ! ground geopotential
     162    !phis(:)=0.
     163    ps(:) = preff
     164    CALL pression (ip1jmp1, ap, bp, ps, p)
     165
     166    IF (pressure_exner) THEN
     167      CALL exner_hyb(ip1jmp1, ps, p, pks, pk)
     168    else
     169      CALL exner_milieu(ip1jmp1, ps, p, pks, pk)
     170    endif
     171    CALL massdair(p, masse)
     172  ENDIF
     173
     174  IF (llm == 1) THEN
     175    ! initialize fields for the shallow water case, if required
     176    IF (.NOT.read_start) THEN
     177      phis(:) = 0.
     178      q(:, :, :) = 0
     179      CALL sw_case_williamson91_6(vcov, ucov, teta, masse, ps)
     180    endif
     181  ENDIF
     182
     183  academic_case: if (iflag_phys >= 2) THEN
     184    ! initializations
     185
     186    ! 1. local parameters
     187    ! by convention, winter is in the southern hemisphere
     188    ! Geostrophic wind or no wind?
     189    ok_geost = .TRUE.
     190    CALL getin('ok_geost', ok_geost)
     191    ! Constants for Newtonian relaxation and friction
     192    k_f = 1.                !friction
     193    CALL getin('k_j', k_f)
     194    k_f = 1. / (daysec * k_f)
     195    k_c_s = 4.  !cooling surface
     196    CALL getin('k_c_s', k_c_s)
     197    k_c_s = 1. / (daysec * k_c_s)
     198    k_c_a = 40. !cooling free atm
     199    CALL getin('k_c_a', k_c_a)
     200    k_c_a = 1. / (daysec * k_c_a)
     201    ! Constants for Teta equilibrium profile
     202    teta0 = 315.     ! mean Teta (S.H. 315K)
     203    CALL getin('teta0', teta0)
     204    ttp = 200.       ! Tropopause temperature (S.H. 200K)
     205    CALL getin('ttp', ttp)
     206    eps = 0.         ! Deviation to N-S symmetry(~0-20K)
     207    CALL getin('eps', eps)
     208    delt_y = 60.     ! Merid Temp. Gradient (S.H. 60K)
     209    CALL getin('delt_y', delt_y)
     210    delt_z = 10.     ! Vertical Gradient (S.H. 10K)
     211    CALL getin('delt_z', delt_z)
     212    ! Polar vortex
     213    ok_pv = .FALSE.
     214    CALL getin('ok_pv', ok_pv)
     215    phi_pv = -50.            ! Latitude of edge of vortex
     216    CALL getin('phi_pv', phi_pv)
     217    phi_pv = phi_pv * pi / 180.
     218    dphi_pv = 5.             ! Width of the edge
     219    CALL getin('dphi_pv', dphi_pv)
     220    dphi_pv = dphi_pv * pi / 180.
     221    gam_pv = 4.              ! -dT/dz vortex (in K/km)
     222    CALL getin('gam_pv', gam_pv)
     223    tetanoise = 0.005
     224    CALL getin('tetanoise', tetanoise)
     225
     226    ! 2. Initialize fields towards which to relax
     227    ! Friction
     228    knewt_g = k_c_a
     229    DO l = 1, llm
     230      zsig = presnivs(l) / preff
     231      knewt_t(l) = (k_c_s - k_c_a) * MAX(0., (zsig - 0.7) / 0.3)
     232      kfrict(l) = k_f * MAX(0., (zsig - 0.7) / 0.3)
     233    ENDDO
     234    DO j = 1, jjp1
     235      clat4((j - 1) * iip1 + 1:j * iip1) = cos(rlatu(j))**4
     236    ENDDO
     237
     238    ! Potential temperature
     239    DO l = 1, llm
     240      zsig = presnivs(l) / preff
     241      tetastrat = ttp * zsig**(-kappa)
     242      tetapv = tetastrat
     243      IF ((ok_pv).AND.(zsig<0.1)) THEN
     244        tetapv = tetastrat * (zsig * 10.)**(kappa * cpp * gam_pv / 1000. / g)
     245      ENDIF
     246      DO j = 1, jjp1
     247        ! Troposphere
     248        ddsin = sin(rlatu(j))
     249        tetajl(j, l) = teta0 - delt_y * ddsin * ddsin + eps * ddsin &
     250                - delt_z * (1. - ddsin * ddsin) * log(zsig)
     251        IF (planet_type=="giant") THEN
     252          tetajl(j, l) = teta0 + (delt_y * &
     253                  ((sin(rlatu(j) * 3.14159 * eps + 0.0001))**2)   &
     254                  / ((rlatu(j) * 3.14159 * eps + 0.0001)**2))     &
     255                  - delt_z * log(zsig)
     256        endif
     257        ! Profil stratospherique isotherme (+vortex)
     258        w_pv = (1. - tanh((rlatu(j) - phi_pv) / dphi_pv)) / 2.
     259        tetastrat = tetastrat * (1. - w_pv) + tetapv * w_pv
     260        tetajl(j, l) = MAX(tetajl(j, l), tetastrat)
     261      ENDDO
     262    ENDDO
     263
     264    !          CALL writefield('theta_eq',tetajl)
     265
     266    do l = 1, llm
     267      do j = 1, jjp1
     268        do i = 1, iip1
     269          ij = (j - 1) * iip1 + i
     270          tetarappel(ij, l) = tetajl(j, l)
    158271        enddo
    159      enddo
    160      PRINT*,'phis=',minval(phis),maxval(phis),'g=',g
    161 
    162      ! ground geopotential
    163      !phis(:)=0.
    164      ps(:)=preff
    165      CALL pression ( ip1jmp1, ap, bp, ps, p       )
    166 
    167      IF (pressure_exner) THEN
    168        CALL exner_hyb( ip1jmp1, ps, p, pks, pk)
    169      else
    170        CALL exner_milieu(ip1jmp1,ps,p,pks,pk)
    171      endif
    172      CALL massdair(p,masse)
    173   ENDIF
    174 
    175   IF (llm == 1) THEN
    176      ! initialize fields for the shallow water case, if required
    177      IF (.NOT.read_start) THEN
    178         phis(:)=0.
    179         q(:,:,:)=0
    180         CALL sw_case_williamson91_6(vcov,ucov,teta,masse,ps)
    181      endif
    182   ENDIF
    183 
    184   academic_case: if (iflag_phys >= 2) THEN
    185      ! initializations
    186 
    187      ! 1. local parameters
    188      ! by convention, winter is in the southern hemisphere
    189      ! Geostrophic wind or no wind?
    190      ok_geost=.TRUE.
    191      CALL getin('ok_geost',ok_geost)
    192      ! Constants for Newtonian relaxation and friction
    193      k_f=1.                !friction
    194      CALL getin('k_j',k_f)
    195      k_f=1./(daysec*k_f)
    196      k_c_s=4.  !cooling surface
    197      CALL getin('k_c_s',k_c_s)
    198      k_c_s=1./(daysec*k_c_s)
    199      k_c_a=40. !cooling free atm
    200      CALL getin('k_c_a',k_c_a)
    201      k_c_a=1./(daysec*k_c_a)
    202      ! Constants for Teta equilibrium profile
    203      teta0=315.     ! mean Teta (S.H. 315K)
    204      CALL getin('teta0',teta0)
    205      ttp=200.       ! Tropopause temperature (S.H. 200K)
    206      CALL getin('ttp',ttp)
    207      eps=0.         ! Deviation to N-S symmetry(~0-20K)
    208      CALL getin('eps',eps)
    209      delt_y=60.     ! Merid Temp. Gradient (S.H. 60K)
    210      CALL getin('delt_y',delt_y)
    211      delt_z=10.     ! Vertical Gradient (S.H. 10K)
    212      CALL getin('delt_z',delt_z)
    213      ! Polar vortex
    214      ok_pv=.FALSE.
    215      CALL getin('ok_pv',ok_pv)
    216      phi_pv=-50.            ! Latitude of edge of vortex
    217      CALL getin('phi_pv',phi_pv)
    218      phi_pv=phi_pv*pi/180.
    219      dphi_pv=5.             ! Width of the edge
    220      CALL getin('dphi_pv',dphi_pv)
    221      dphi_pv=dphi_pv*pi/180.
    222      gam_pv=4.              ! -dT/dz vortex (in K/km)
    223      CALL getin('gam_pv',gam_pv)
    224      tetanoise=0.005
    225      CALL getin('tetanoise',tetanoise)
    226 
    227      ! 2. Initialize fields towards which to relax
    228      ! Friction
    229      knewt_g=k_c_a
    230      DO l=1,llm
    231         zsig=presnivs(l)/preff
    232         knewt_t(l)=(k_c_s-k_c_a)*MAX(0.,(zsig-0.7)/0.3)
    233         kfrict(l)=k_f*MAX(0.,(zsig-0.7)/0.3)
    234      ENDDO
    235      DO j=1,jjp1
    236         clat4((j-1)*iip1+1:j*iip1)=cos(rlatu(j))**4
    237      ENDDO
    238 
    239      ! Potential temperature
    240      DO l=1,llm
    241         zsig=presnivs(l)/preff
    242         tetastrat=ttp*zsig**(-kappa)
    243         tetapv=tetastrat
    244         IF ((ok_pv).AND.(zsig<0.1)) THEN
    245            tetapv=tetastrat*(zsig*10.)**(kappa*cpp*gam_pv/1000./g)
    246         ENDIF
    247         DO j=1,jjp1
    248            ! Troposphere
    249            ddsin=sin(rlatu(j))
    250            tetajl(j,l)=teta0-delt_y*ddsin*ddsin+eps*ddsin &
    251                 -delt_z*(1.-ddsin*ddsin)*log(zsig)
    252            IF (planet_type=="giant") THEN
    253              tetajl(j,l)=teta0+(delt_y*                   &
    254                 ((sin(rlatu(j)*3.14159*eps+0.0001))**2)   &
    255                 / ((rlatu(j)*3.14159*eps+0.0001)**2))     &
    256                 -delt_z*log(zsig)
    257            endif
    258            ! Profil stratospherique isotherme (+vortex)
    259            w_pv=(1.-tanh((rlatu(j)-phi_pv)/dphi_pv))/2.
    260            tetastrat=tetastrat*(1.-w_pv)+tetapv*w_pv             
    261            tetajl(j,l)=MAX(tetajl(j,l),tetastrat) 
    262         ENDDO
    263      ENDDO
    264 
    265      !          CALL writefield('theta_eq',tetajl)
    266 
    267      do l=1,llm
    268         do j=1,jjp1
    269            do i=1,iip1
    270               ij=(j-1)*iip1+i
    271               tetarappel(ij,l)=tetajl(j,l)
    272            enddo
     272      enddo
     273    enddo
     274
     275    ! 3. Initialize fields (if necessary)
     276    IF (.NOT. read_start) THEN
     277      ! bulk initialization of temperature
     278      IF (iflag_phys>10000) THEN
     279        ! Particular case to impose a constant temperature T0=0.01*iflag_physx
     280        teta(:, :) = 0.01 * iflag_phys / (pk(:, :) / cpp)
     281      ELSE
     282        teta(:, :) = tetarappel(:, :)
     283      ENDIF
     284      ! geopotential
     285      CALL geopot(ip1jmp1, teta, pk, pks, phis, phi)
     286
     287      DO l = 1, llm
     288        PRINT*, 'presnivs,play,l', presnivs(l), (pk(1, l) / cpp)**(1. / kappa) * preff
     289        !pks(ij) = (cpp/preff) * ps(ij)
     290        !pk(ij,1) = .5*pks(ij)
     291        ! pk = cpp * (p/preff)^kappa
     292      ENDDO
     293
     294      ! winds
     295      IF (ok_geost) THEN
     296        CALL ugeostr(phi, ucov)
     297      else
     298        ucov(:, :) = 0.
     299      endif
     300      vcov(:, :) = 0.
     301
     302      ! bulk initialization of tracers
     303      IF (planet_type=="earth") THEN
     304        ! Earth: first two tracers will be water
     305        do iq = 1, nqtot
     306          q(:, :, iq) = 0.
     307          IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:, :, iq) = 1.e-10
     308          IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:, :, iq) = 1.e-15
     309
     310          ! CRisi: init des isotopes
     311          ! distill de Rayleigh très simplifiée
     312          iName = tracers(iq)%iso_iName
     313          IF (niso <= 0 .OR. iName <= 0) CYCLE
     314          iPhase = tracers(iq)%iso_iPhase
     315          iqParent = tracers(iq)%iqParent
     316          IF(tracers(iq)%iso_iZone == 0) THEN
     317            IF (tnat1) THEN
     318              tnat = 1.0
     319              alpha_ideal = 1.0
     320              WRITE(*, *) 'Attention dans iniacademic: alpha_ideal=1'
     321            else
     322              IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
     323                      CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
     324            endif
     325            q(:, :, iq) = q(:, :, iqParent) * tnat * (q(:, :, iqParent) / 30.e-3)**(alpha_ideal - 1.)
     326          ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
     327            IF(tracers(iq)%iso_iZone == 1) THEN
     328              ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
     329              ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
     330              q(:, :, iq) = q(:, :, iqIsoPha(iName, iPhase))
     331            else !IF(tracers(iq)%iso_iZone == 1) THEN
     332              q(:, :, iq) = 0.
     333            endif !IF(tracers(iq)%iso_iZone == 1) THEN
     334          END IF !IF(tracers(iq)%iso_iZone == 0) THEN
    273335        enddo
    274      enddo
    275 
    276      ! 3. Initialize fields (if necessary)
    277      IF (.NOT. read_start) THEN
    278         ! bulk initialization of temperature
    279         IF (iflag_phys>10000) THEN
    280         ! Particular case to impose a constant temperature T0=0.01*iflag_physx
    281            teta(:,:)= 0.01*iflag_phys/(pk(:,:)/cpp)
    282         ELSE
    283            teta(:,:)=tetarappel(:,:)
    284         ENDIF
    285         ! geopotential
    286         CALL geopot(ip1jmp1,teta,pk,pks,phis,phi)
    287 
    288         DO l=1,llm
    289           PRINT*,'presnivs,play,l',presnivs(l),(pk(1,l)/cpp)**(1./kappa)*preff
    290          !pks(ij) = (cpp/preff) * ps(ij)
    291          !pk(ij,1) = .5*pks(ij)
    292          ! pk = cpp * (p/preff)^kappa
    293         ENDDO
    294 
    295         ! winds
    296         IF (ok_geost) THEN
    297            CALL ugeostr(phi,ucov)
    298         else
    299            ucov(:,:)=0.
    300         endif
    301         vcov(:,:)=0.
    302 
    303         ! bulk initialization of tracers
    304         IF (planet_type=="earth") THEN
    305            ! Earth: first two tracers will be water
    306            do iq=1,nqtot
    307               q(:,:,iq)=0.
    308               IF(tracers(iq)%name == addPhase('H2O', 'g')) q(:,:,iq)=1.e-10
    309               IF(tracers(iq)%name == addPhase('H2O', 'l')) q(:,:,iq)=1.e-15
    310 
    311               ! CRisi: init des isotopes
    312               ! distill de Rayleigh très simplifiée
    313               iName    = tracers(iq)%iso_iName
    314               IF (niso <= 0 .OR. iName <= 0) CYCLE
    315               iPhase   = tracers(iq)%iso_iPhase
    316               iqParent = tracers(iq)%iqParent
    317               IF(tracers(iq)%iso_iZone == 0) THEN
    318                  IF (tnat1) THEN
    319                          tnat=1.0
    320                          alpha_ideal=1.0
    321                          WRITE(*,*) 'Attention dans iniacademic: alpha_ideal=1'
    322                  else
    323                     IF(getKey('tnat', tnat, isoName(iName)) .OR. getKey('alpha', alpha_ideal, isoName(iName))) &
    324                     CALL abort_gcm(TRIM(modname), 'missing isotopic parameters', 1)
    325                  endif
    326                  q(:,:,iq) = q(:,:,iqParent)*tnat*(q(:,:,iqParent)/30.e-3)**(alpha_ideal-1.)
    327               ELSE !IF(tracers(iq)%iso_iZone == 0) THEN
    328                  IF(tracers(iq)%iso_iZone == 1) THEN
    329                     ! correction le 14 mai 2024 pour que tous les traceurs soient de la couleur 1.
    330                     ! Sinon, on va avoir des porblèmes de conservation de masse de traceurs.
    331                     q(:,:,iq) = q(:,:,iqIsoPha(iName,iPhase))
    332                  else !IF(tracers(iq)%iso_iZone == 1) THEN
    333                     q(:,:,iq) = 0.
    334                  endif !IF(tracers(iq)%iso_iZone == 1) THEN
    335               END IF !IF(tracers(iq)%iso_iZone == 0) THEN
    336            enddo
    337         else
    338            q(:,:,:)=0
    339         endif ! of if (planet_type=="earth")
    340 
    341         CALL check_isotopes_seq(q,1,ip1jmp1,'iniacademic_loc')
    342 
    343         ! add random perturbation to temperature
    344         idum  = -1
    345         zz = ran1(idum)
    346         idum  = 0
    347         do l=1,llm
    348            do ij=iip2,ip1jm
    349               teta(ij,l)=teta(ij,l)*(1.+tetanoise*ran1(idum))
    350            enddo
     336      else
     337        q(:, :, :) = 0
     338      endif ! of if (planet_type=="earth")
     339
     340      CALL check_isotopes_seq(q, 1, ip1jmp1, 'iniacademic_loc')
     341
     342      ! add random perturbation to temperature
     343      idum = -1
     344      zz = ran1(idum)
     345      idum = 0
     346      do l = 1, llm
     347        do ij = iip2, ip1jm
     348          teta(ij, l) = teta(ij, l) * (1. + tetanoise * ran1(idum))
    351349        enddo
    352 
    353         ! maintain periodicity in longitude
    354         do l=1,llm
    355            do ij=1,ip1jmp1,iip1
    356               teta(ij+iim,l)=teta(ij,l)
    357            enddo
     350      enddo
     351
     352      ! maintain periodicity in longitude
     353      do l = 1, llm
     354        do ij = 1, ip1jmp1, iip1
     355          teta(ij + iim, l) = teta(ij, l)
    358356        enddo
    359 
    360      ENDIF ! of IF (.NOT. read_start)
     357      enddo
     358
     359    ENDIF ! of IF (.NOT. read_start)
    361360  ENDIF academic_case
    362361
Note: See TracChangeset for help on using the changeset viewer.