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

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

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