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