source: LMDZ5/trunk/libf/dyn3d/gcm.F90 @ 2597

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago

Cleanup in the dynamics: get rid of comconst.h, make it a module comconst_mod.
EM

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