source: LMDZ5/branches/testing/libf/dyn3d/gcm.F @ 2258

Last change on this file since 2258 was 2258, checked in by Laurent Fairhead, 10 years ago

Merged trunk changes 2216:2237 into testing branch

  • 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.1 KB
Line 
1!
2! $Id: gcm.F 2258 2015-04-13 08:21:09Z fairhead $
3!
4c
5c
6      PROGRAM gcm
7
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      USE ioipsl_getincom
13#endif
14
15
16#ifdef CPP_XIOS
17    ! ug Pour les sorties XIOS
18        USE wxios
19#endif
20
21      USE filtreg_mod
22      USE infotrac
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_para, ONLY : klon_mpi_para_nb
29#endif
30
31!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
33! A nettoyer. On ne veut qu'une ou deux routines d'interface
34! dynamique -> physique pour l'initialisation
35#ifdef CPP_PHYS
36!      USE dimphy
37!      USE comgeomphy
38#endif
39!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
40
41      IMPLICIT NONE
42
43c      ......   Version  du 10/01/98    ..........
44
45c             avec  coordonnees  verticales hybrides
46c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
47
48c=======================================================================
49c
50c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
51c   -------
52c
53c   Objet:
54c   ------
55c
56c   GCM LMD nouvelle grille
57c
58c=======================================================================
59c
60c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
61c      et possibilite d'appeler une fonction f(y)  a derivee tangente
62c      hyperbolique a la  place de la fonction a derivee sinusoidale.
63c  ... Possibilite de choisir le schema pour l'advection de
64c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
65c
66c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
67c      Pour Van-Leer iadv=10
68c
69c-----------------------------------------------------------------------
70c   Declarations:
71c   -------------
72
73#include "dimensions.h"
74#include "paramet.h"
75#include "comconst.h"
76#include "comdissnew.h"
77#include "comvert.h"
78#include "comgeom.h"
79#include "logic.h"
80#include "temps.h"
81!!!!!!!!!!!#include "control.h"
82#include "ener.h"
83#include "description.h"
84#include "serre.h"
85!#include "com_io_dyn.h"
86#include "iniprint.h"
87#include "tracstoke.h"
88#ifdef INCA
89! Only INCA needs these informations (from the Earth's physics)
90!#include "indicesol.h"
91#endif
92
93      REAL zdtvr
94
95c   variables dynamiques
96      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
97      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
98      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
99      REAL ps(ip1jmp1)                       ! pression  au sol
100      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
101      REAL masse(ip1jmp1,llm)                ! masse d'air
102      REAL phis(ip1jmp1)                     ! geopotentiel au sol
103      REAL phi(ip1jmp1,llm)                  ! geopotentiel
104      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
114      INTEGER ij,iq,l,i,j
115
116
117      real time_step, t_wrt, t_ops
118
119      LOGICAL first
120
121!      LOGICAL call_iniphys
122!      data call_iniphys/.true./
123
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
129      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
132      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
143c-----------------------------------------------------------------------
144c   Initialisations:
145c   ----------------
146
147      abort_message = 'last timestep reached'
148      modname = 'gcm'
149      descript = 'Run GCM LMDZ'
150      lafin    = .FALSE.
151      dynhist_file = 'dyn_hist.nc'
152      dynhistave_file = 'dyn_hist_ave.nc'
153
154
155
156c----------------------------------------------------------------------
157c  lecture des fichiers gcm.def ou run.def
158c  ---------------------------------------
159c
160      CALL conf_gcm( 99, .TRUE.)
161      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
162     s "iphysiq must be a multiple of iperiod", 1)
163
164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165! Initialisation de XIOS
166!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167
168#ifdef CPP_XIOS
169        CALL wxios_init("LMDZ")
170#endif
171
172
173!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174! FH 2008/05/02
175! A nettoyer. On ne veut qu'une ou deux routines d'interface
176! dynamique -> physique pour l'initialisation
177#ifdef CPP_PHYS
178      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
179!      call InitComgeomphy ! now done in iniphysiq
180#endif
181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182c-----------------------------------------------------------------------
183c   Choix du calendrier
184c   -------------------
185
186c      calend = 'earth_365d'
187
188#ifdef CPP_IOIPSL
189      if (calend == 'earth_360d') then
190        call ioconf_calendar('360d')
191        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
192      else if (calend == 'earth_365d') then
193        call ioconf_calendar('noleap')
194        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
195      else if (calend == 'gregorian') then
196        call ioconf_calendar('gregorian')
197        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
198      else
199        abort_message = 'Mauvais choix de calendrier'
200        call abort_gcm(modname,abort_message,1)
201      endif
202#endif
203c-----------------------------------------------------------------------
204
205      IF (type_trac == 'inca') THEN
206#ifdef INCA
207      call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
208     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
209      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
210#endif
211      END IF
212c
213c
214c------------------------------------
215c   Initialisation partie parallele
216c------------------------------------
217
218c
219c
220c-----------------------------------------------------------------------
221c   Initialisation des traceurs
222c   ---------------------------
223c  Choix du nombre de traceurs et du schema pour l'advection
224c  dans fichier traceur.def, par default ou via INCA
225      call infotrac_init
226
227c Allocation de la tableau q : champs advectes   
228      allocate(q(ip1jmp1,llm,nqtot))
229
230c-----------------------------------------------------------------------
231c   Lecture de l'etat initial :
232c   ---------------------------
233
234c  lecture du fichier start.nc
235      if (read_start) then
236      ! we still need to run iniacademic to initialize some
237      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
238        if (iflag_phys.ne.1) then
239          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
240        endif
241
242!        if (planet_type.eq."earth") then
243! Load an Earth-format start file
244         CALL dynetat0("start.nc",vcov,ucov,
245     &              teta,q,masse,ps,phis, time_0)
246!        endif ! of if (planet_type.eq."earth")
247       
248c       write(73,*) 'ucov',ucov
249c       write(74,*) 'vcov',vcov
250c       write(75,*) 'teta',teta
251c       write(76,*) 'ps',ps
252c       write(77,*) 'q',q
253
254      endif ! of if (read_start)
255
256      IF (type_trac == 'inca') THEN
257#ifdef INCA
258         call init_inca_dim(klon,llm,iim,jjm,
259     $        rlonu,rlatu,rlonv,rlatv)
260#endif
261      END IF
262
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
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 \`a 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 a 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
327
328c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
329c        write(lunout,*)
330c     .  'GCM: Attention les dates initiales lues dans le fichier'
331c        write(lunout,*)
332c     .  ' restart ne correspondent pas a celles lues dans '
333c        write(lunout,*)' gcm.def'
334c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
335c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
336c        if (raz_date .ne. 1) then
337c          write(lunout,*)
338c     .    'GCM: On garde les dates du fichier restart'
339c        else
340c          annee_ref = anneeref
341c          day_ref = dayref
342c          day_ini = dayref
343c          itau_dyn = 0
344c          itau_phy = 0
345c          time_0 = 0.
346c          write(lunout,*)
347c     .   'GCM: On reinitialise a la date lue dans gcm.def'
348c        endif
349c      ELSE
350c        raz_date = 0
351c      endif
352
353#ifdef CPP_IOIPSL
354      mois = 1
355      heure = 0.
356      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
357      jH_ref = jD_ref - int(jD_ref)
358      jD_ref = int(jD_ref)
359
360      call ioconf_startdate(INT(jD_ref), jH_ref)
361
362      write(lunout,*)'DEBUG'
363      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
364      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
365      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
366      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
367      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
368#else
369! Ehouarn: we still need to define JD_ref and JH_ref
370! and since we don't know how many days there are in a year
371! we set JD_ref to 0 (this should be improved ...)
372      jD_ref=0
373      jH_ref=0
374#endif
375
376
377      if (iflag_phys.eq.1) then
378      ! these initialisations have already been done (via iniacademic)
379      ! if running in SW or Newtonian mode
380c-----------------------------------------------------------------------
381c   Initialisation des constantes dynamiques :
382c   ------------------------------------------
383        dtvr = zdtvr
384        CALL iniconst
385
386c-----------------------------------------------------------------------
387c   Initialisation de la geometrie :
388c   --------------------------------
389        CALL inigeom
390
391c-----------------------------------------------------------------------
392c   Initialisation du filtre :
393c   --------------------------
394        CALL inifilr
395      endif ! of if (iflag_phys.eq.1)
396c
397c-----------------------------------------------------------------------
398c   Initialisation de la dissipation :
399c   ----------------------------------
400
401      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
402     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
403
404c-----------------------------------------------------------------------
405c   Initialisation de la physique :
406c   -------------------------------
407
408      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
409! Physics:
410#ifdef CPP_PHYS
411         CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys,
412     &                rlatu,rlonv,aire,cu,cv,rad,g,r,cpp,
413     &                iflag_phys)
414#endif
415      ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
416
417c  numero de stockage pour les fichiers de redemarrage:
418
419c-----------------------------------------------------------------------
420c   Initialisation des I/O :
421c   ------------------------
422
423
424      if (nday>=0) then
425         day_end = day_ini + nday
426      else
427         day_end = day_ini - nday/day_step
428      endif
429      WRITE(lunout,300)day_ini,day_end
430 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
431
432#ifdef CPP_IOIPSL
433      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
434      write (lunout,301)jour, mois, an
435      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
436      write (lunout,302)jour, mois, an
437 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
438 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
439#endif
440
441!      if (planet_type.eq."earth") then
442! Write an Earth-format restart file
443
444        CALL dynredem0("restart.nc", day_end, phis)
445!      endif
446
447      ecripar = .TRUE.
448
449#ifdef CPP_IOIPSL
450      time_step = zdtvr
451      if (ok_dyn_ins) then
452        ! initialize output file for instantaneous outputs
453        ! t_ops = iecri * daysec ! do operations every t_ops
454        t_ops =((1.0*iecri)/day_step) * daysec 
455        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
456        CALL inithist(day_ref,annee_ref,time_step,
457     &              t_ops,t_wrt)
458      endif
459
460      IF (ok_dyn_ave) THEN
461        ! initialize output file for averaged outputs
462        t_ops = iperiod * time_step ! do operations every t_ops
463        t_wrt = periodav * daysec   ! write output every t_wrt
464        CALL initdynav(day_ref,annee_ref,time_step,
465     &       t_ops,t_wrt)
466      END IF
467      dtav = iperiod*dtvr/daysec
468#endif
469! #endif of #ifdef CPP_IOIPSL
470
471c  Choix des frequences de stokage pour le offline
472c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
473c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
474      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
475      istphy=istdyn/iphysiq     
476
477
478c
479c-----------------------------------------------------------------------
480c   Integration temporelle du modele :
481c   ----------------------------------
482
483c       write(78,*) 'ucov',ucov
484c       write(78,*) 'vcov',vcov
485c       write(78,*) 'teta',teta
486c       write(78,*) 'ps',ps
487c       write(78,*) 'q',q
488
489
490      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
491
492      END
493
Note: See TracBrowser for help on using the repository browser.