source: dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3dmem/gcm.F @ 3822

Last change on this file since 3822 was 3822, checked in by millour, 10 years ago

A couple of bug fixes.
Now the bench (in debug mode) yields identical results in seq/mpi/omp/mpi_omp, and also identical restart files to rev 5 (ie before any modifications to LMDZ5 source files).
EM

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