source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3d/gcm.F @ 1446

Last change on this file since 1446 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.1 KB
Line 
1!
2! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $
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      USE filtreg_mod
16      USE infotrac
17      USE control_mod
18
19!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
21! A nettoyer. On ne veut qu'une ou deux routines d'interface
22! dynamique -> physique pour l'initialisation
23! Ehouarn: for now these only apply to Earth:
24#ifdef CPP_EARTH
25      USE dimphy
26      USE comgeomphy
27      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
28#endif
29!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
30
31      IMPLICIT NONE
32
33c      ......   Version  du 10/01/98    ..........
34
35c             avec  coordonnees  verticales hybrides
36c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
37
38c=======================================================================
39c
40c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
41c   -------
42c
43c   Objet:
44c   ------
45c
46c   GCM LMD nouvelle grille
47c
48c=======================================================================
49c
50c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
51c      et possibilite d'appeler une fonction f(y)  a derivee tangente
52c      hyperbolique a la  place de la fonction a derivee sinusoidale.
53c  ... Possibilite de choisir le schema pour l'advection de
54c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
55c
56c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
57c      Pour Van-Leer iadv=10
58c
59c-----------------------------------------------------------------------
60c   Declarations:
61c   -------------
62
63#include "dimensions.h"
64#include "paramet.h"
65#include "comconst.h"
66#include "comdissnew.h"
67#include "comvert.h"
68#include "comgeom.h"
69#include "logic.h"
70#include "temps.h"
71!!!!!!!!!!!#include "control.h"
72#include "ener.h"
73#include "description.h"
74#include "serre.h"
75!#include "com_io_dyn.h"
76#include "iniprint.h"
77#include "tracstoke.h"
78#ifdef INCA
79! Only INCA needs these informations (from the Earth's physics)
80#include "indicesol.h"
81#endif
82      INTEGER         longcles
83      PARAMETER     ( longcles = 20 )
84      REAL  clesphy0( longcles )
85      SAVE  clesphy0
86
87
88
89      REAL zdtvr
90      INTEGER nbetatmoy, nbetatdem,nbetat
91
92c   variables dynamiques
93      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
94      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
95      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
96      REAL ps(ip1jmp1)                       ! pression  au sol
97      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
98      REAL pks(ip1jmp1)                      ! exner au  sol
99      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
100      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu 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
124      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
125c+jld variables test conservation energie
126c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
127C     Tendance de la temp. potentiel d (theta)/ d t due a la
128C     tansformation d'energie cinetique en energie thermique
129C     cree par la dissipation
130      REAL dhecdt(ip1jmp1,llm)
131c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
132c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
133      CHARACTER (len=15) :: ztit
134c-jld
135
136
137      character (len=80) :: dynhist_file, dynhistave_file
138      character (len=20) :: modname
139      character (len=80) :: abort_message
140! locales pour gestion du temps
141      INTEGER :: an, mois, jour
142      REAL :: heure
143
144
145c-----------------------------------------------------------------------
146c    variables pour l'initialisation de la physique :
147c    ------------------------------------------------
148      INTEGER ngridmx
149      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
150      REAL zcufi(ngridmx),zcvfi(ngridmx)
151      REAL latfi(ngridmx),lonfi(ngridmx)
152      REAL airefi(ngridmx)
153      SAVE latfi, lonfi, airefi
154
155c-----------------------------------------------------------------------
156c   Initialisations:
157c   ----------------
158
159      abort_message = 'last timestep reached'
160      modname = 'gcm'
161      descript = 'Run GCM LMDZ'
162      lafin    = .FALSE.
163      dynhist_file = 'dyn_hist.nc'
164      dynhistave_file = 'dyn_hist_ave.nc'
165
166
167
168c----------------------------------------------------------------------
169c  lecture des fichiers gcm.def ou run.def
170c  ---------------------------------------
171c
172! Ehouarn: dump possibility of using defrun
173!#ifdef CPP_IOIPSL
174      CALL conf_gcm( 99, .TRUE. , clesphy0 )
175!#else
176!      CALL defrun( 99, .TRUE. , clesphy0 )
177!#endif
178
179!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180! FH 2008/05/02
181! A nettoyer. On ne veut qu'une ou deux routines d'interface
182! dynamique -> physique pour l'initialisation
183! Ehouarn : temporarily (?) keep this only for Earth
184      if (planet_type.eq."earth") then
185#ifdef CPP_EARTH
186      CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
187      call InitComgeomphy
188#endif
189      endif
190!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191c-----------------------------------------------------------------------
192c   Choix du calendrier
193c   -------------------
194
195c      calend = 'earth_365d'
196
197#ifdef CPP_IOIPSL
198      if (calend == 'earth_360d') then
199        call ioconf_calendar('360d')
200        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
201      else if (calend == 'earth_365d') then
202        call ioconf_calendar('noleap')
203        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
204      else if (calend == 'earth_366d') then
205        call ioconf_calendar('gregorian')
206        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
207      else
208        abort_message = 'Mauvais choix de calendrier'
209        call abort_gcm(modname,abort_message,1)
210      endif
211#endif
212c-----------------------------------------------------------------------
213
214      IF (config_inca /= 'none') THEN
215#ifdef INCA
216      call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
217     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
218      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
219#endif
220      END IF
221c
222c
223c------------------------------------
224c   Initialisation partie parallele
225c------------------------------------
226
227c
228c
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
265      IF (config_inca /= 'none') THEN
266#ifdef INCA
267         call init_inca_dim(klon,llm,iim,jjm,
268     $        rlonu,rlatu,rlonv,rlatv)
269#endif
270      END IF
271
272
273c le cas echeant, creation d un etat initial
274      IF (prt_level > 9) WRITE(lunout,*)
275     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
276      if (.not.read_start) then
277         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
278      endif
279
280
281c-----------------------------------------------------------------------
282c   Lecture des parametres de controle pour la simulation :
283c   -------------------------------------------------------
284c  on recalcule eventuellement le pas de temps
285
286      IF(MOD(day_step,iperiod).NE.0) THEN
287        abort_message =
288     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
289        call abort_gcm(modname,abort_message,1)
290      ENDIF
291
292      IF(MOD(day_step,iphysiq).NE.0) THEN
293        abort_message =
294     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
295        call abort_gcm(modname,abort_message,1)
296      ENDIF
297
298      zdtvr    = daysec/REAL(day_step)
299        IF(dtvr.NE.zdtvr) THEN
300         WRITE(lunout,*)
301     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
302        ENDIF
303
304C
305C on remet le calendrier à zero si demande
306c
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
326
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
375c  nombre d'etats dans les fichiers demarrage et histoire
376      nbetatdem = nday / iecri
377      nbetatmoy = nday / periodav + 1
378
379      if (iflag_phys.eq.1) then
380      ! these initialisations have already been done (via iniacademic)
381      ! if running in SW or Newtonian mode
382c-----------------------------------------------------------------------
383c   Initialisation des constantes dynamiques :
384c   ------------------------------------------
385        dtvr = zdtvr
386        CALL iniconst
387
388c-----------------------------------------------------------------------
389c   Initialisation de la geometrie :
390c   --------------------------------
391        CALL inigeom
392
393c-----------------------------------------------------------------------
394c   Initialisation du filtre :
395c   --------------------------
396        CALL inifilr
397      endif ! of if (iflag_phys.eq.1)
398c
399c-----------------------------------------------------------------------
400c   Initialisation de la dissipation :
401c   ----------------------------------
402
403      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
404     *                tetagdiv, tetagrot , tetatemp              )
405
406c-----------------------------------------------------------------------
407c   Initialisation de la physique :
408c   -------------------------------
409
410      IF (call_iniphys.and.(iflag_phys.eq.1)) THEN
411         latfi(1)=rlatu(1)
412         lonfi(1)=0.
413         zcufi(1) = cu(1)
414         zcvfi(1) = cv(1)
415         DO j=2,jjm
416            DO i=1,iim
417               latfi((j-2)*iim+1+i)= rlatu(j)
418               lonfi((j-2)*iim+1+i)= rlonv(i)
419               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
420               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
421            ENDDO
422         ENDDO
423         latfi(ngridmx)= rlatu(jjp1)
424         lonfi(ngridmx)= 0.
425         zcufi(ngridmx) = cu(ip1jm+1)
426         zcvfi(ngridmx) = cv(ip1jm-iim)
427         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
428         WRITE(lunout,*)
429     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
430! Earth:
431         if (planet_type.eq."earth") then
432#ifdef CPP_EARTH
433         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
434     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
435#endif
436         endif ! of if (planet_type.eq."earth")
437         call_iniphys=.false.
438      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
439!#endif
440
441c  numero de stockage pour les fichiers de redemarrage:
442
443c-----------------------------------------------------------------------
444c   Initialisation des I/O :
445c   ------------------------
446
447
448      day_end = day_ini + nday
449      WRITE(lunout,300)day_ini,day_end
450 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
451
452#ifdef CPP_IOIPSL
453      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
454      write (lunout,301)jour, mois, an
455      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
456      write (lunout,302)jour, mois, an
457 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
458 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
459#endif
460
461!      if (planet_type.eq."earth") then
462! Write an Earth-format restart file
463        CALL dynredem0("restart.nc", day_end, phis)
464!      endif
465
466      ecripar = .TRUE.
467
468#ifdef CPP_IOIPSL
469      time_step = zdtvr
470      if (ok_dyn_ins) then
471        ! initialize output file for instantaneous outputs
472        ! t_ops = iecri * daysec ! do operations every t_ops
473        t_ops =((1.0*iecri)/day_step) * daysec 
474        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
475        CALL inithist(day_ref,annee_ref,time_step,
476     &              t_ops,t_wrt)
477      endif
478
479      IF (ok_dyn_ave) THEN
480        ! initialize output file for averaged outputs
481        t_ops = iperiod * time_step ! do operations every t_ops
482        t_wrt = periodav * daysec   ! write output every t_wrt
483        CALL initdynav(day_ref,annee_ref,time_step,
484     &       t_ops,t_wrt)
485      END IF
486      dtav = iperiod*dtvr/daysec
487#endif
488! #endif of #ifdef CPP_IOIPSL
489
490c  Choix des frequences de stokage pour le offline
491c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
492c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
493      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
494      istphy=istdyn/iphysiq     
495
496
497c
498c-----------------------------------------------------------------------
499c   Integration temporelle du modele :
500c   ----------------------------------
501
502c       write(78,*) 'ucov',ucov
503c       write(78,*) 'vcov',vcov
504c       write(78,*) 'teta',teta
505c       write(78,*) 'ps',ps
506c       write(78,*) 'q',q
507
508
509      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
510     .              time_0)
511
512      END
513
Note: See TracBrowser for help on using the repository browser.