source: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F @ 492

Last change on this file since 492 was 492, checked in by emillour, 13 years ago

Common dynamics: updates to keep up with LMDZ5 Earth (rev 1605)
See file "DOC/chantiers/commit_importants.log" for details.
EM

File size: 17.8 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_PHYS
25      USE dimphy
26      USE comgeomphy
27#endif
28#ifdef CPP_EARTH
29      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
30#endif
31!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32
33      IMPLICIT NONE
34
35c      ......   Version  du 10/01/98    ..........
36
37c             avec  coordonnees  verticales hybrides
38c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
39
40c=======================================================================
41c
42c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
43c   -------
44c
45c   Objet:
46c   ------
47c
48c   GCM LMD nouvelle grille
49c
50c=======================================================================
51c
52c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
53c      et possibilite d'appeler une fonction f(y)  a derivee tangente
54c      hyperbolique a la  place de la fonction a derivee sinusoidale.
55c  ... Possibilite de choisir le schema pour l'advection de
56c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
57c
58c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
59c      Pour Van-Leer iadv=10
60c
61c-----------------------------------------------------------------------
62c   Declarations:
63c   -------------
64
65#include "dimensions.h"
66#include "paramet.h"
67#include "comconst.h"
68#include "comdissnew.h"
69#include "comvert.h"
70#include "comgeom.h"
71#include "logic.h"
72#include "temps.h"
73!!!!!!!!!!!#include "control.h"
74#include "ener.h"
75#include "description.h"
76#include "serre.h"
77!#include "com_io_dyn.h"
78#include "iniprint.h"
79#include "tracstoke.h"
80#ifdef INCA
81! Only INCA needs these informations (from the Earth's physics)
82#include "indicesol.h"
83#endif
84
85
86      REAL zdtvr
87
88c   variables dynamiques
89      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
90      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
91      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
92      REAL ps(ip1jmp1)                       ! pression  au sol
93      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
94      REAL pks(ip1jmp1)                      ! exner au  sol
95      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
96      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
97      REAL masse(ip1jmp1,llm)                ! masse d'air
98      REAL phis(ip1jmp1)                     ! geopotentiel au sol
99      REAL phi(ip1jmp1,llm)                  ! geopotentiel
100      REAL w(ip1jmp1,llm)                    ! vitesse verticale
101
102c variables dynamiques intermediaire pour le transport
103
104c   variables pour le fichier histoire
105      REAL dtav      ! intervalle de temps elementaire
106
107      REAL time_0
108
109      LOGICAL lafin
110      INTEGER ij,iq,l,i,j
111
112
113      real time_step, t_wrt, t_ops
114
115      LOGICAL first
116
117      LOGICAL call_iniphys
118      data call_iniphys/.true./
119
120      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
121c+jld variables test conservation energie
122c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
123C     Tendance de la temp. potentiel d (theta)/ d t due a la
124C     tansformation d'energie cinetique en energie thermique
125C     cree par la dissipation
126      REAL dhecdt(ip1jmp1,llm)
127c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
128c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
129      CHARACTER (len=15) :: ztit
130c-jld
131
132
133      character (len=80) :: dynhist_file, dynhistave_file
134      character (len=20) :: modname
135      character (len=80) :: abort_message
136! locales pour gestion du temps
137      INTEGER :: an, mois, jour
138      REAL :: heure
139
140
141c-----------------------------------------------------------------------
142c    variables pour l'initialisation de la physique :
143c    ------------------------------------------------
144      INTEGER ngridmx
145      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
146      REAL zcufi(ngridmx),zcvfi(ngridmx)
147      REAL latfi(ngridmx),lonfi(ngridmx)
148      REAL airefi(ngridmx)
149      SAVE latfi, lonfi, airefi
150
151c-----------------------------------------------------------------------
152c   Initialisations:
153c   ----------------
154
155      abort_message = 'last timestep reached'
156      modname = 'gcm'
157      descript = 'Run GCM LMDZ'
158      lafin    = .FALSE.
159      dynhist_file = 'dyn_hist.nc'
160      dynhistave_file = 'dyn_hist_ave.nc'
161
162
163
164c----------------------------------------------------------------------
165c  lecture des fichiers gcm.def ou run.def
166c  ---------------------------------------
167c
168! Ehouarn: dump possibility of using defrun
169!#ifdef CPP_IOIPSL
170      CALL conf_gcm( 99, .TRUE. )
171!#else
172!      CALL defrun( 99, .TRUE. , clesphy0 )
173!#endif
174
175!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
176! FH 2008/05/02
177! A nettoyer. On ne veut qu'une ou deux routines d'interface
178! dynamique -> physique pour l'initialisation
179! Ehouarn : temporarily (?) keep this only for Earth
180!      if (planet_type.eq."earth") then
181!#ifdef CPP_EARTH
182#ifdef CPP_PHYS
183      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
184      call initcomgeomphy
185#endif
186!      endif
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188c
189c Initialisations pour Cp(T) Venus
190      call ini_cpdet
191c
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 if (calend == 'titan') then
209!        call ioconf_calendar('titan')
210        write(lunout,*)'CALENDRIER CHOISI: Titan'
211        abort_message = 'A FAIRE...'
212        call abort_gcm(modname,abort_message,1)
213      else if (calend == 'venus') then
214!        call ioconf_calendar('venus')
215        write(lunout,*)'CALENDRIER CHOISI: Venus'
216        abort_message = 'A FAIRE...'
217        call abort_gcm(modname,abort_message,1)
218      else
219        abort_message = 'Mauvais choix de calendrier'
220        call abort_gcm(modname,abort_message,1)
221      endif
222#endif
223c-----------------------------------------------------------------------
224
225      IF (type_trac == 'inca') THEN
226#ifdef INCA
227      call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
228     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
229      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
230#endif
231      END IF
232c
233c
234c------------------------------------
235c   Initialisation partie parallele
236c------------------------------------
237
238c
239c
240c-----------------------------------------------------------------------
241c   Initialisation des traceurs
242c   ---------------------------
243c  Choix du nombre de traceurs et du schema pour l'advection
244c  dans fichier traceur.def, par default ou via INCA
245      call infotrac_init
246
247c Allocation de la tableau q : champs advectes   
248      allocate(q(ip1jmp1,llm,nqtot))
249
250c-----------------------------------------------------------------------
251c   Lecture de l'etat initial :
252c   ---------------------------
253
254c  lecture du fichier start.nc
255      if (read_start) then
256      ! we still need to run iniacademic to initialize some
257      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
258        if (iflag_phys.ne.1) then
259          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
260        endif
261
262        if (planet_type.eq."mars") then
263! POUR MARS, METTRE UNE FONCTION A PART, genre dynetat0_mars
264         abort_message = 'dynetat0_mars A FAIRE'
265         call abort_gcm(modname,abort_message,0)
266        else
267         CALL dynetat0("start.nc",vcov,ucov,
268     &              teta,q,masse,ps,phis, time_0)
269        endif ! of if (planet_type.eq."mars")
270       
271c       write(73,*) 'ucov',ucov
272c       write(74,*) 'vcov',vcov
273c       write(75,*) 'teta',teta
274c       write(76,*) 'ps',ps
275c       write(77,*) 'q',q
276
277      endif ! of if (read_start)
278
279      IF (type_trac == 'inca') THEN
280#ifdef INCA
281         call init_inca_dim(klon,llm,iim,jjm,
282     $        rlonu,rlatu,rlonv,rlatv)
283#endif
284      END IF
285
286
287c le cas echeant, creation d un etat initial
288      IF (prt_level > 9) WRITE(lunout,*)
289     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
290      if (.not.read_start) then
291         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
292      endif
293
294
295c-----------------------------------------------------------------------
296c   Lecture des parametres de controle pour la simulation :
297c   -------------------------------------------------------
298c  on recalcule eventuellement le pas de temps
299
300      IF(MOD(day_step,iperiod).NE.0) THEN
301        abort_message =
302     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
303        call abort_gcm(modname,abort_message,1)
304      ENDIF
305
306      IF(MOD(day_step,iphysiq).NE.0) THEN
307        abort_message =
308     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
309        call abort_gcm(modname,abort_message,1)
310      ENDIF
311
312      zdtvr    = daysec/REAL(day_step)
313        IF(dtvr.NE.zdtvr) THEN
314         WRITE(lunout,*)
315     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
316        ENDIF
317
318C
319C on remet le calendrier à zero si demande
320c
321      IF (start_time /= starttime) then
322        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
323     &,' fichier restart ne correspond pas à celle lue dans le run.def'
324        IF (raz_date == 1) then
325          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
326          start_time = starttime
327        ELSE
328          WRITE(lunout,*)'Je m''arrete'
329          CALL abort
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! Ce n'est defini pour l'instant que pour la Terre...
381      if (planet_type.eq.'earth') then
382      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
383      jH_ref = jD_ref - int(jD_ref)
384      jD_ref = int(jD_ref)
385
386      call ioconf_startdate(INT(jD_ref), jH_ref)
387
388      write(lunout,*)'DEBUG'
389      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
390      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
391      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
392      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
393      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
394      else
395! A voir pour Titan et Venus
396        jD_ref=0
397        jH_ref=0
398      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
399      write(lunout,*)jD_ref,jH_ref
400      endif ! planet_type
401#else
402! Ehouarn: we still need to define JD_ref and JH_ref
403! and since we don't know how many days there are in a year
404! we set JD_ref to 0 (this should be improved ...)
405      jD_ref=0
406      jH_ref=0
407#endif
408
409      if (iflag_phys.eq.1) then
410      ! these initialisations have already been done (via iniacademic)
411      ! if running in SW or Newtonian mode
412c-----------------------------------------------------------------------
413c   Initialisation des constantes dynamiques :
414c   ------------------------------------------
415        dtvr = zdtvr
416        CALL iniconst
417
418c-----------------------------------------------------------------------
419c   Initialisation de la geometrie :
420c   --------------------------------
421        CALL inigeom
422
423c-----------------------------------------------------------------------
424c   Initialisation du filtre :
425c   --------------------------
426        CALL inifilr
427      endif ! of if (iflag_phys.eq.1)
428c
429c-----------------------------------------------------------------------
430c   Initialisation de la dissipation :
431c   ----------------------------------
432
433      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
434     *                tetagdiv, tetagrot , tetatemp              )
435
436c-----------------------------------------------------------------------
437c   Initialisation de la physique :
438c   -------------------------------
439
440      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
441         latfi(1)=rlatu(1)
442         lonfi(1)=0.
443         zcufi(1) = cu(1)
444         zcvfi(1) = cv(1)
445         DO j=2,jjm
446            DO i=1,iim
447               latfi((j-2)*iim+1+i)= rlatu(j)
448               lonfi((j-2)*iim+1+i)= rlonv(i)
449               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
450               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
451            ENDDO
452         ENDDO
453         latfi(ngridmx)= rlatu(jjp1)
454         lonfi(ngridmx)= 0.
455         zcufi(ngridmx) = cu(ip1jm+1)
456         zcvfi(ngridmx) = cv(ip1jm-iim)
457         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
458         WRITE(lunout,*)
459     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
460
461! Initialisation de la physique: pose probleme quand on tourne
462! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
463! Il faut une cle CPP_PHYS
464#ifdef CPP_PHYS
465         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
466     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
467#endif ! CPP_PHYS
468         call_iniphys=.false.
469      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
470!#endif
471
472c  numero de stockage pour les fichiers de redemarrage:
473
474c-----------------------------------------------------------------------
475c   Initialisation des I/O :
476c   ------------------------
477
478
479      day_end = day_ini + nday
480      WRITE(lunout,300)day_ini,day_end
481 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
482
483#ifdef CPP_IOIPSL
484! Ce n'est defini pour l'instant que pour la Terre...
485      if (planet_type.eq.'earth') then
486      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
487      write (lunout,301)jour, mois, an
488      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
489      write (lunout,302)jour, mois, an
490      else
491! A voir pour Titan et Venus
492      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
493      endif ! planet_type
494
495 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
496 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
497#endif
498
499#ifdef CPP_EARTH
500! Create start file (startphy.nc) and boundary conditions (limit.nc)
501! for the Earth verstion
502       if (iflag_phys>=100) then
503          call iniaqua(ngridmx,latfi,lonfi,iflag_phys)
504       endif
505#endif
506
507      if (planet_type.eq."mars") then
508! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
509         abort_message = 'dynredem0_mars A FAIRE'
510         call abort_gcm(modname,abort_message,0)
511      else
512        CALL dynredem0("restart.nc", day_end, phis)
513      endif ! of if (planet_type.eq."mars")
514
515      ecripar = .TRUE.
516
517#ifdef CPP_IOIPSL
518      time_step = zdtvr
519      if (ok_dyn_ins) then
520        ! initialize output file for instantaneous outputs
521        ! t_ops = iecri * daysec ! do operations every t_ops
522        t_ops =((1.0*iecri)/day_step) * daysec 
523        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
524        CALL inithist(day_ref,annee_ref,time_step,
525     &              t_ops,t_wrt)
526      endif
527
528      IF (ok_dyn_ave) THEN
529        ! initialize output file for averaged outputs
530        t_ops = iperiod * time_step ! do operations every t_ops
531        t_wrt = periodav * daysec   ! write output every t_wrt
532        CALL initdynav(day_ref,annee_ref,time_step,
533     &       t_ops,t_wrt)
534      END IF
535      dtav = iperiod*dtvr/daysec
536#endif
537! #endif of #ifdef CPP_IOIPSL
538
539c  Choix des frequences de stokage pour le offline
540c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
541c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
542      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
543      istphy=istdyn/iphysiq     
544
545
546c
547c-----------------------------------------------------------------------
548c   Integration temporelle du modele :
549c   ----------------------------------
550
551c       write(78,*) 'ucov',ucov
552c       write(78,*) 'vcov',vcov
553c       write(78,*) 'teta',teta
554c       write(78,*) 'ps',ps
555c       write(78,*) 'q',q
556
557
558      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,
559     .              time_0)
560
561      END
562
Note: See TracBrowser for help on using the repository browser.