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

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

Some partial cleanup on uses of "dimensions.h" in physics.
At this point 3D gcm compiles and bench seems to run fine :-)
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 de la physique :
404c   -------------------------------
405      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
406! Physics:
407#ifdef CPP_PHYS
408         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
409     &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,
410     &                iflag_phys)
411#endif
412      ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
413
414
415c-----------------------------------------------------------------------
416c   Initialisation des dimensions d'INCA :
417c   --------------------------------------
418      IF (type_trac == 'inca') THEN
419!$OMP PARALLEL
420#ifdef INCA
421         CALL init_inca_dim(klon_omp,llm,iim,jjm,
422     $        rlonu,rlatu,rlonv,rlatv)
423#endif
424!$OMP END PARALLEL
425      END IF
426
427c-----------------------------------------------------------------------
428c   Initialisation des I/O :
429c   ------------------------
430
431
432      if (nday>=0) then
433         day_end = day_ini + nday
434      else
435         day_end = day_ini - nday/day_step
436      endif
437 
438      WRITE(lunout,300)day_ini,day_end
439 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
440
441#ifdef CPP_IOIPSL
442      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
443      write (lunout,301)jour, mois, an
444      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
445      write (lunout,302)jour, mois, an
446 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
447 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
448#endif
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.