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