source: LMDZ5/branches/LMDZ6_rc0/libf/dyn3d/gcm.F @ 2942

Last change on this file since 2942 was 2275, checked in by Laurent Fairhead, 9 years ago

Backport of trunk r2229 in LMDZ6_rc0 branch

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