source: LMDZ5/trunk/libf/dyn3d/gcm.F @ 1972

Last change on this file since 1972 was 1930, checked in by lguez, 11 years ago

abort, dfloat and pause are not in the Fortran standard. Replaced
abort by abort_gcm and dfloat by dble. Note: I modified dyn3dpar files
that were identical to dyn3d modified files.

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