source: LMDZ5/trunk/libf/dyn3dmem/gcm.F @ 1657

Last change on this file since 1657 was 1657, checked in by Laurent Fairhead, 12 years ago

Phasage de la dynamique parallèle localisée (petite mémoire) avec la branche LMDZ4V5.0-dev (fin de la branche)
Validation effectuée par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modèle donnée paY. Meurdesoif et la version phasée avec la r1399 (fin de la branche LMDZ4V5.0-dev)


Phasing of the localised (low memory) parallel dynamics package with the LMDZ4V5.0-dev version of LMDZ
Validation consisted in comparing output debug files (u, v, t, q, masse, etc... ) of a no physics simulation
run with the version of the code given by Y. Meurdesoif and this version phased with r1399 (end of the LMDZ4V5.0-dev branch)

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