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

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

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