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

Last change on this file since 2438 was 2438, checked in by lguez, 8 years ago

Bug fix

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