source: LMDZ6/trunk/libf/dyn3d/gcm.F90 @ 5258

Last change on this file since 5258 was 5250, checked in by abarral, 5 weeks ago

Wrap uses of cpp key CPP_PHYS

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