source: dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d/gcm.F90 @ 3817

Last change on this file since 3817 was 3817, checked in by millour, 10 years ago

Further cleanup and removal of references to iniprint.h.
Also added bench testcase 48x36x19.
EM

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