source: LMDZ6/trunk/libf/dyn3dmem/gcm.F90 @ 4379

Last change on this file since 4379 was 4361, checked in by lguez, 2 years ago

Change calendar attribute "360d" to "360_day"

"360_day" is the correct attribute according to CF convention. "360d"
leads to an error when opening a history file with xarray.

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