source: LMDZ6/branches/Amaury_dev/libf/dyn3d/gcm.F90 @ 5136

Last change on this file since 5136 was 5136, checked in by abarral, 3 months ago

Put comgeom.h, comgeom2.h into modules

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 14.5 KB
Line 
1! $Id: gcm.F90 5136 2024-07-28 14:17:54Z abarral $
2
3
4
5PROGRAM gcm
6
7  USE IOIPSL
8  USE lmdz_wxios  ! ug Pour les sorties XIOS
9
10  USE lmdz_filtreg, ONLY: inifilr
11  USE infotrac, ONLY: nqtot, init_infotrac
12  USE control_mod
13  USE mod_const_mpi, ONLY: COMM_LMDZ
14  USE temps_mod, ONLY: calend, start_time, annee_ref, day_ref, &
15          itau_dyn, itau_phy, day_ini, jD_ref, jH_ref, day_end
16  USE comconst_mod, ONLY: cpp, daysec, dtphys, dtvr, g, r, rad
17  USE logic_mod, ONLY: ecripar, iflag_phys, read_start
18
19  USE lmdz_cppkeys_wrapper, ONLY: CPPKEY_PHYS
20  USE lmdz_description, ONLY: descript
21
22  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
23  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
24  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
25  ! dynamique -> physique pour l'initialisation
26  ! AB 2024/07/20: remplace CPP key by fortran logical, but ^ still relevant, see later use of iniphys later on
27  USE iniphysiq_mod, ONLY: iniphysiq
28  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
29  USE lmdz_iniprint, ONLY: lunout, prt_level
30  USE lmdz_comdissnew, ONLY: lstardis, nitergdiv, nitergrot, niterh, tetagdiv, &
31          tetagrot, tetatemp, coefdis, vert_prof_dissip
32  USE lmdz_comgeom
33
34  IMPLICIT NONE
35
36  !      ......   Version  du 10/01/98    ..........
37
38  !             avec  coordonnees  verticales hybrides
39  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
40
41  !=======================================================================
42
43  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
44  !   -------
45
46  !   Objet:
47  !   ------
48
49  !   GCM LMD nouvelle grille
50
51  !=======================================================================
52
53  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
54  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
55  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
56  !  ... Possibilite de choisir le schema pour l'advection de
57  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
58
59  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
60  !      Pour Van-Leer iadv=10
61
62  !-----------------------------------------------------------------------
63  !   Declarations:
64  !   -------------
65
66  INCLUDE "dimensions.h"
67  INCLUDE "paramet.h"
68  INCLUDE "tracstoke.h"
69
70  REAL zdtvr
71
72  !   variables dynamiques
73  REAL vcov(ip1jm, llm), ucov(ip1jmp1, llm) ! vents covariants
74  REAL teta(ip1jmp1, llm)                 ! temperature potentielle
75  REAL, ALLOCATABLE, DIMENSION(:, :, :) :: q! champs advectes
76  REAL ps(ip1jmp1)                       ! pression  au sol
77  !  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
78  REAL masse(ip1jmp1, llm)                ! masse d'air
79  REAL phis(ip1jmp1)                     ! geopotentiel au sol
80  !  REAL phi(ip1jmp1,llm)                  ! geopotentiel
81  !  REAL w(ip1jmp1,llm)                    ! vitesse verticale
82
83  ! variables dynamiques intermediaire pour le transport
84
85  !   variables pour le fichier histoire
86  REAL dtav      ! intervalle de temps elementaire
87
88  REAL time_0
89
90  LOGICAL lafin
91
92  REAL time_step, t_wrt, t_ops
93
94  !      LOGICAL call_iniphys
95  !      data call_iniphys/.TRUE./
96
97  !+jld variables test conservation energie
98  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
99  !     Tendance de la temp. potentiel d (theta)/ d t due a la
100  !     tansformation d'energie cinetique en energie thermique
101  !     cree par la dissipation
102  !  REAL dhecdt(ip1jmp1,llm)
103  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
104  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
105  !  CHARACTER (len=15) :: ztit
106  !-jld
107
108  CHARACTER (LEN = 80) :: dynhist_file, dynhistave_file
109  CHARACTER (LEN = 20) :: modname
110  CHARACTER (LEN = 80) :: abort_message
111  ! locales pour gestion du temps
112  INTEGER :: an, mois, jour
113  REAL :: heure
114  LOGICAL use_filtre_fft
115
116  !-----------------------------------------------------------------------
117  !   Initialisations:
118  !   ----------------
119
120  abort_message = 'last timestep reached'
121  modname = 'gcm'
122  descript = 'Run GCM LMDZ'
123  lafin = .FALSE.
124  dynhist_file = 'dyn_hist.nc'
125  dynhistave_file = 'dyn_hist_ave.nc'
126
127
128
129  !----------------------------------------------------------------------
130  !  lecture des fichiers gcm.def ou run.def
131  !  ---------------------------------------
132
133  CALL conf_gcm(99, .TRUE.)
134
135  IF (mod(iphysiq, iperiod) /= 0) CALL abort_gcm("conf_gcm", &
136          "iphysiq must be a multiple of iperiod", 1)
137
138  use_filtre_fft = .FALSE.
139  CALL getin('use_filtre_fft', use_filtre_fft)
140  IF (use_filtre_fft) CALL abort_gcm("gcm", 'FFT filter is not available in ' &
141          // 'the sequential version of the dynamics.', 1)
142
143  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
144  ! Initialisation de XIOS
145  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
146
147  IF (using_xios) THEN
148    CALL wxios_init("LMDZ")
149  ENDIF
150
151
152  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
153  ! FH 2008/05/02
154  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
155  ! dynamique -> physique pour l'initialisation
156  !#ifdef CPP_PHYS
157  !  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
158  !  !      CALL InitComgeomphy ! now done in iniphysiq
159  !#endif
160  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
161  !-----------------------------------------------------------------------
162  !   Choix du calendrier
163  !   -------------------
164
165  !      calend = 'earth_365d'
166
167  IF (calend == 'earth_360d') THEN
168    CALL ioconf_calendar('360_day')
169    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
170  ELSE IF (calend == 'earth_365d') THEN
171    CALL ioconf_calendar('noleap')
172    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
173  ELSE IF (calend == 'gregorian') THEN
174    CALL ioconf_calendar('gregorian')
175    WRITE(lunout, *)'CALENDRIER CHOISI: Terrestre bissextile'
176  else
177    abort_message = 'Mauvais choix de calendrier'
178    CALL abort_gcm(modname, abort_message, 1)
179  ENDIF
180  !-----------------------------------------------------------------------
181
182
183  !------------------------------------
184  !   Initialisation partie parallele
185  !------------------------------------
186
187
188  !-----------------------------------------------------------------------
189  !   Initialisation des traceurs
190  !   ---------------------------
191  !  Choix du nombre de traceurs et du schema pour l'advection
192  !  dans fichier traceur.def, par default ou via INCA
193  CALL init_infotrac
194
195  ! Allocation de la tableau q : champs advectes   
196  allocate(q(ip1jmp1, llm, nqtot))
197
198  !-----------------------------------------------------------------------
199  !   Lecture de l'etat initial :
200  !   ---------------------------
201
202  !  lecture du fichier start.nc
203  IF (read_start) THEN
204    ! we still need to run iniacademic to initialize some
205    ! constants & fields, if we run the 'newtonian' or 'SW' cases:
206    IF (iflag_phys/=1) THEN
207      CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
208    endif
209
210    !        if (planet_type.EQ."earth") THEN
211    ! Load an Earth-format start file
212    CALL dynetat0("start.nc", vcov, ucov, &
213            teta, q, masse, ps, phis, time_0)
214    !        endif ! of if (planet_type.EQ."earth")
215
216    !       WRITE(73,*) 'ucov',ucov
217    !       WRITE(74,*) 'vcov',vcov
218    !       WRITE(75,*) 'teta',teta
219    !       WRITE(76,*) 'ps',ps
220    !       WRITE(77,*) 'q',q
221
222  ENDIF ! of if (read_start)
223
224
225  ! le cas echeant, creation d un etat initial
226  IF (prt_level > 9) WRITE(lunout, *) &
227          'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
228  IF (.NOT.read_start) THEN
229    start_time = 0.
230    annee_ref = anneeref
231    CALL iniacademic(vcov, ucov, teta, q, masse, ps, phis, time_0)
232  ENDIF
233
234
235  !-----------------------------------------------------------------------
236  !   Lecture des parametres de controle pour la simulation :
237  !   -------------------------------------------------------
238  !  on recalcule eventuellement le pas de temps
239
240  IF(MOD(day_step, iperiod)/=0) THEN
241    abort_message = &
242            'Il faut choisir un nb de pas par jour multiple de iperiod'
243    CALL abort_gcm(modname, abort_message, 1)
244  ENDIF
245
246  IF(MOD(day_step, iphysiq)/=0) THEN
247    abort_message = &
248            'Il faut choisir un nb de pas par jour multiple de iphysiq'
249    CALL abort_gcm(modname, abort_message, 1)
250  ENDIF
251
252  zdtvr = daysec / REAL(day_step)
253  IF(dtvr/=zdtvr) THEN
254    WRITE(lunout, *) &
255            'WARNING!!! changement de pas de temps', dtvr, '>', zdtvr
256  ENDIF
257
258  ! on remet le calendrier \`a zero si demande
259
260  IF (start_time /= starttime) THEN
261    WRITE(lunout, *)' GCM: Attention l''heure de depart lue dans le' &
262            , ' fichier restart ne correspond pas a celle lue dans le run.def'
263    IF (raz_date == 1) THEN
264      WRITE(lunout, *)'Je prends l''heure lue dans run.def'
265      start_time = starttime
266    ELSE
267      CALL abort_gcm("gcm", "'Je m''arrete'", 1)
268    ENDIF
269  ENDIF
270  IF (raz_date == 1) THEN
271    annee_ref = anneeref
272    day_ref = dayref
273    day_ini = dayref
274    itau_dyn = 0
275    itau_phy = 0
276    time_0 = 0.
277    WRITE(lunout, *) &
278            'GCM: On reinitialise a la date lue dans gcm.def'
279  ELSE IF (annee_ref /= anneeref .OR. day_ref /= dayref) THEN
280    WRITE(lunout, *) &
281            'GCM: Attention les dates initiales lues dans le fichier'
282    WRITE(lunout, *) &
283            ' restart ne correspondent pas a celles lues dans '
284    WRITE(lunout, *)' gcm.def'
285    WRITE(lunout, *)' annee_ref=', annee_ref, " anneeref=", anneeref
286    WRITE(lunout, *)' day_ref=', day_ref, " dayref=", dayref
287    WRITE(lunout, *)' Pas de remise a zero'
288  ENDIF
289
290  !      if (annee_ref .NE. anneeref .OR. day_ref .NE. dayref) THEN
291  !        WRITE(lunout,*)
292  !     .  'GCM: Attention les dates initiales lues dans le fichier'
293  !        WRITE(lunout,*)
294  !     .  ' restart ne correspondent pas a celles lues dans '
295  !        WRITE(lunout,*)' gcm.def'
296  !        WRITE(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
297  !        WRITE(lunout,*)' day_ref=',day_ref," dayref=",dayref
298  !        if (raz_date .NE. 1) THEN
299  !          WRITE(lunout,*)
300  !     .    'GCM: On garde les dates du fichier restart'
301  !        else
302  !          annee_ref = anneeref
303  !          day_ref = dayref
304  !          day_ini = dayref
305  !          itau_dyn = 0
306  !          itau_phy = 0
307  !          time_0 = 0.
308  !          WRITE(lunout,*)
309  !     .   'GCM: On reinitialise a la date lue dans gcm.def'
310  !        endif
311  !      ELSE
312  !        raz_date = 0
313  !      endif
314
315  mois = 1
316  heure = 0.
317  CALL ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
318  jH_ref = jD_ref - int(jD_ref)
319  jD_ref = int(jD_ref)
320
321  CALL ioconf_startdate(INT(jD_ref), jH_ref)
322
323  WRITE(lunout, *)'DEBUG'
324  WRITE(lunout, *)'annee_ref, mois, day_ref, heure, jD_ref'
325  WRITE(lunout, *)annee_ref, mois, day_ref, heure, jD_ref
326  CALL ju2ymds(jD_ref + jH_ref, an, mois, jour, heure)
327  WRITE(lunout, *)'jD_ref+jH_ref,an, mois, jour, heure'
328  WRITE(lunout, *)jD_ref + jH_ref, an, mois, jour, heure
329
330  IF (iflag_phys==1) THEN
331    ! these initialisations have already been done (via iniacademic)
332    ! if running in SW or Newtonian mode
333    !-----------------------------------------------------------------------
334    !   Initialisation des constantes dynamiques :
335    !   ------------------------------------------
336    dtvr = zdtvr
337    CALL iniconst
338
339    !-----------------------------------------------------------------------
340    !   Initialisation de la geometrie :
341    !   --------------------------------
342    CALL inigeom
343
344    !-----------------------------------------------------------------------
345    !   Initialisation du filtre :
346    !   --------------------------
347    CALL inifilr
348  ENDIF ! of if (iflag_phys.EQ.1)
349
350  !-----------------------------------------------------------------------
351  !   Initialisation de la dissipation :
352  !   ----------------------------------
353
354  CALL inidissip(lstardis, nitergdiv, nitergrot, niterh, &
355          tetagdiv, tetagrot, tetatemp, vert_prof_dissip)
356
357  !  numero de stockage pour les fichiers de redemarrage:
358
359  !-----------------------------------------------------------------------
360  !   Initialisation des I/O :
361  !   ------------------------
362
363  IF (nday>=0) THEN
364    day_end = day_ini + nday
365  else
366    day_end = day_ini - nday / day_step
367  ENDIF
368  WRITE(lunout, 300)day_ini, day_end
369  300 FORMAT('1'/, 15x, 'run du jour', i7, 2x, 'au jour', i7//)
370
371  CALL ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
372  write (lunout, 301)jour, mois, an
373  CALL ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
374  write (lunout, 302)jour, mois, an
375  301 FORMAT('1'/, 15x, 'run du ', i2, '/', i2, '/', i4)
376  302 FORMAT('1'/, 15x, '    au ', i2, '/', i2, '/', i4)
377
378  !-----------------------------------------------------------------------
379  !   Initialisation de la physique :
380  !   -------------------------------
381
382  IF ((iflag_phys==1).OR.(iflag_phys>=100)) THEN
383    ! Physics:
384    IF (CPPKEY_PHYS) THEN
385      CALL iniphysiq(iim, jjm, llm, &
386              (jjm - 1) * iim + 2, comm_lmdz, &
387              daysec, day_ini, dtphys / nsplit_phys, &
388              rlatu, rlatv, rlonu, rlonv, aire, cu, cv, rad, g, r, cpp, &
389              iflag_phys)
390    END IF
391  ENDIF ! of IF ((iflag_phys==1).OR.(iflag_phys>=100))
392
393  !      if (planet_type.EQ."earth") THEN
394  ! Write an Earth-format restart file
395
396  CALL dynredem0("restart.nc", day_end, phis)
397  !      endif
398
399  ecripar = .TRUE.
400
401  time_step = zdtvr
402  IF (ok_dyn_ins) THEN
403    ! initialize output file for instantaneous outputs
404    ! t_ops = iecri * daysec ! do operations every t_ops
405    t_ops = ((1.0 * iecri) / day_step) * daysec
406    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
407    CALL inithist(day_ref, annee_ref, time_step, &
408            t_ops, t_wrt)
409  ENDIF
410
411  IF (ok_dyn_ave) THEN
412    ! initialize output file for averaged outputs
413    t_ops = iperiod * time_step ! do operations every t_ops
414    t_wrt = periodav * daysec   ! write output every t_wrt
415    CALL initdynav(day_ref, annee_ref, time_step, &
416            t_ops, t_wrt)
417  END IF
418  dtav = iperiod * dtvr / daysec
419
420  !  Choix des frequences de stokage pour le offline
421  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
422  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
423  istdyn = day_step / 4     ! stockage toutes les 6h=1jour/12
424  istphy = istdyn / iphysiq
425
426  !-----------------------------------------------------------------------
427  !   Integration temporelle du modele :
428  !   ----------------------------------
429
430  !       WRITE(78,*) 'ucov',ucov
431  !       WRITE(78,*) 'vcov',vcov
432  !       WRITE(78,*) 'teta',teta
433  !       WRITE(78,*) 'ps',ps
434  !       WRITE(78,*) 'q',q
435
436  CALL leapfrog(ucov, vcov, teta, ps, masse, phis, q, time_0)
437
438END PROGRAM gcm
Note: See TracBrowser for help on using the repository browser.