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

Last change on this file since 2360 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

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