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

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

Phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ4 (HEAD)
Validation effectuee par comparaison des fichiers de sorties debug (u, v, t, q, masse, etc ...) d'une simulation sans physique
faite avec la version du modele donnee par Y. Meurdesoif et la version phasee avec la r1428 (fin du tronc LMDZ4)


Phasing of the localised (low memory) parallel dynamics package with the LMDZ4 trunk 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 r1428 (HEAD of the LMDZ4 trunk)

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