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

Last change on this file since 5268 was 5268, checked in by abarral, 3 days ago

.f90 <-> .F90 depending on cpp key use

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