source: LMDZ5/trunk/libf/dyn3dmem/gcm.F90 @ 2597

Last change on this file since 2597 was 2597, checked in by Ehouarn Millour, 8 years ago

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