source: LMDZ5/trunk/libf/dyn3dpar/gcm.F @ 2346

Last change on this file since 2346 was 2346, checked in by Ehouarn Millour, 9 years ago

Physics/dynamics separation:

  • remove all references to dimensions.h from physics. nbp_lon (==iim) , nbp_lat (==jjm+1) and nbp_lev (==llm) from mod_grid_phy_lmdz should be used instead.
  • added module regular_lonlat_mod in phy_common to store information about the global (lon-lat) grid cell boundaries and centers.

EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.4 KB
Line 
1!
2! $Id: gcm.F 2346 2015-08-21 15:13:46Z 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
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!#endif
167!      CALL set_bands
168!#ifdef CPP_PHYS
169      CALL Init_interface_dyn_phys
170#endif
171      CALL barrier
172
173      CALL set_bands
174      if (mpi_rank==0) call WriteBands
175      call SetDistrib(jj_Nb_Caldyn)
176
177c$OMP PARALLEL
178      call Init_Mod_hallo
179c$OMP END PARALLEL
180
181!#ifdef CPP_PHYS
182!c$OMP PARALLEL
183!      call InitComgeomphy ! now done in iniphysiq
184!c$OMP END PARALLEL
185!#endif
186
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188! Initialisation de XIOS
189!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190
191
192c-----------------------------------------------------------------------
193c   Choix du calendrier
194c   -------------------
195
196c      calend = 'earth_365d'
197
198#ifdef CPP_IOIPSL
199      if (calend == 'earth_360d') then
200        call ioconf_calendar('360d')
201        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
202      else if (calend == 'earth_365d') then
203        call ioconf_calendar('noleap')
204        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
205      else if (calend == 'earth_366d') then
206        call ioconf_calendar('gregorian')
207        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
208      else
209        abort_message = 'Mauvais choix de calendrier'
210        call abort_gcm(modname,abort_message,1)
211      endif
212#endif
213
214      IF (type_trac == 'inca') THEN
215#ifdef INCA
216         call init_const_lmdz(
217     $        nbtr,anneeref,dayref,
218     $        iphysiq,day_step,nday,
219     $        nbsrf, is_oce,is_sic,
220     $        is_ter,is_lic, calend)
221
222         call init_inca_para(
223     $        iim,jjm+1,llm,klon_glo,mpi_size,
224     $        distrib_phys,COMM_LMDZ)
225#endif
226      END IF
227
228c-----------------------------------------------------------------------
229c   Initialisation des traceurs
230c   ---------------------------
231c  Choix du nombre de traceurs et du schema pour l'advection
232c  dans fichier traceur.def, par default ou via INCA
233      call infotrac_init
234
235c Allocation de la tableau q : champs advectes   
236      ALLOCATE(q(ip1jmp1,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(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("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(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 à 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 à 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          call abort_gcm("gcm", "'Je m''arrete'", 1)
305        ENDIF
306      ENDIF
307      IF (raz_date == 1) THEN
308        annee_ref = anneeref
309        day_ref = dayref
310        day_ini = dayref
311        itau_dyn = 0
312        itau_phy = 0
313        time_0 = 0.
314        write(lunout,*)
315     .   'GCM: On reinitialise a la date lue dans gcm.def'
316      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
317        write(lunout,*)
318     .  'GCM: Attention les dates initiales lues dans le fichier'
319        write(lunout,*)
320     .  ' restart ne correspondent pas a celles lues dans '
321        write(lunout,*)' gcm.def'
322        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
323        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
324        write(lunout,*)' Pas de remise a zero'
325      ENDIF
326c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
327c        write(lunout,*)
328c     .  'GCM: Attention les dates initiales lues dans le fichier'
329c        write(lunout,*)
330c     .  ' restart ne correspondent pas a celles lues dans '
331c        write(lunout,*)' gcm.def'
332c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
333c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
334c        if (raz_date .ne. 1) then
335c          write(lunout,*)
336c     .    'GCM: On garde les dates du fichier restart'
337c        else
338c          annee_ref = anneeref
339c          day_ref = dayref
340c          day_ini = dayref
341c          itau_dyn = 0
342c          itau_phy = 0
343c          time_0 = 0.
344c          write(lunout,*)
345c     .   'GCM: On reinitialise a la date lue dans gcm.def'
346c        endif
347c      ELSE
348c        raz_date = 0
349c      endif
350
351#ifdef CPP_IOIPSL
352      mois = 1
353      heure = 0.
354      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
355      jH_ref = jD_ref - int(jD_ref)
356      jD_ref = int(jD_ref)
357
358      call ioconf_startdate(INT(jD_ref), jH_ref)
359
360      write(lunout,*)'DEBUG'
361      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
362      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
363      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
364      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
365      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
366#else
367! Ehouarn: we still need to define JD_ref and JH_ref
368! and since we don't know how many days there are in a year
369! we set JD_ref to 0 (this should be improved ...)
370      jD_ref=0
371      jH_ref=0
372#endif
373
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_p("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(day_ref,annee_ref,time_step,
475     &                   t_ops,t_wrt)
476!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
477!     .        t_ops, t_wrt, histaveid)
478        END IF
479      ENDIF
480      dtav = iperiod*dtvr/daysec
481#endif
482! #endif of #ifdef CPP_IOIPSL
483
484c  Choix des frequences de stokage pour le offline
485c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
486c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
487      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
488      istphy=istdyn/iphysiq     
489
490
491c
492c-----------------------------------------------------------------------
493c   Integration temporelle du modele :
494c   ----------------------------------
495
496c       write(78,*) 'ucov',ucov
497c       write(78,*) 'vcov',vcov
498c       write(78,*) 'teta',teta
499c       write(78,*) 'ps',ps
500c       write(78,*) 'q',q
501
502c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
503      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,time_0)
504c$OMP END PARALLEL
505
506
507      END
508
Note: See TracBrowser for help on using the repository browser.