source: LMDZ6/trunk/libf/dyn3d/gcm.f90

Last change on this file was 5272, checked in by abarral, 2 days ago

Turn paramet.h into a module

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