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

Last change on this file since 966 was 965, checked in by emillour, 12 years ago

Common dynamics and generic/universal GCM:

  • LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
  • LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
  • LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

EM

File size: 17.6 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: the following are needed with (parallel) physics:
24#ifdef CPP_PHYS
25      USE dimphy
26      USE comgeomphy
27#endif
28#ifdef INCA
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#ifdef CPP_PHYS
180      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
181      call initcomgeomphy
182#endif
183!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
184c
185c Initialisations pour Cp(T) Venus
186      call ini_cpdet
187c
188c-----------------------------------------------------------------------
189c   Choix du calendrier
190c   -------------------
191
192c      calend = 'earth_365d'
193
194#ifdef CPP_IOIPSL
195      if (calend == 'earth_360d') then
196        call ioconf_calendar('360d')
197        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
198      else if (calend == 'earth_365d') then
199        call ioconf_calendar('noleap')
200        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
201      else if (calend == 'earth_366d') then
202        call ioconf_calendar('gregorian')
203        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
204      else if (calend == 'titan') then
205!        call ioconf_calendar('titan')
206        write(lunout,*)'CALENDRIER CHOISI: Titan'
207        abort_message = 'A FAIRE...'
208        call abort_gcm(modname,abort_message,1)
209      else if (calend == 'venus') then
210!        call ioconf_calendar('venus')
211        write(lunout,*)'CALENDRIER CHOISI: Venus'
212        abort_message = 'A FAIRE...'
213        call abort_gcm(modname,abort_message,1)
214      else
215        abort_message = 'Mauvais choix de calendrier'
216        call abort_gcm(modname,abort_message,1)
217      endif
218#endif
219c-----------------------------------------------------------------------
220
221      IF (type_trac == 'inca') THEN
222#ifdef INCA
223      call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
224     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
225      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
226#endif
227      END IF
228c
229c
230c------------------------------------
231c   Initialisation partie parallele
232c------------------------------------
233
234c
235c
236c-----------------------------------------------------------------------
237c   Initialisation des traceurs
238c   ---------------------------
239c  Choix du nombre de traceurs et du schema pour l'advection
240c  dans fichier traceur.def, par default ou via INCA
241      call infotrac_init
242
243c Allocation de la tableau q : champs advectes   
244      allocate(q(ip1jmp1,llm,nqtot))
245
246c-----------------------------------------------------------------------
247c   Lecture de l'etat initial :
248c   ---------------------------
249
250c  lecture du fichier start.nc
251      if (read_start) then
252      ! we still need to run iniacademic to initialize some
253      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
254        if (iflag_phys.ne.1) then
255          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
256        endif
257
258        if (planet_type.eq."mars") then
259! POUR MARS, METTRE UNE FONCTION A PART, genre dynetat0_mars
260         abort_message = 'dynetat0_mars A FAIRE'
261         call abort_gcm(modname,abort_message,0)
262        else
263         CALL dynetat0("start.nc",vcov,ucov,
264     &              teta,q,masse,ps,phis, time_0)
265        endif ! of if (planet_type.eq."mars")
266       
267c       write(73,*) 'ucov',ucov
268c       write(74,*) 'vcov',vcov
269c       write(75,*) 'teta',teta
270c       write(76,*) 'ps',ps
271c       write(77,*) 'q',q
272
273      endif ! of if (read_start)
274
275      IF (type_trac == 'inca') THEN
276#ifdef INCA
277         call init_inca_dim(klon,llm,iim,jjm,
278     $        rlonu,rlatu,rlonv,rlatv)
279#endif
280      END IF
281
282
283c le cas echeant, creation d un etat initial
284      IF (prt_level > 9) WRITE(lunout,*)
285     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
286      if (.not.read_start) then
287         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
288      endif
289
290
291c-----------------------------------------------------------------------
292c   Lecture des parametres de controle pour la simulation :
293c   -------------------------------------------------------
294c  on recalcule eventuellement le pas de temps
295
296      IF(MOD(day_step,iperiod).NE.0) THEN
297        abort_message =
298     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
299        call abort_gcm(modname,abort_message,1)
300      ENDIF
301
302      IF(MOD(day_step,iphysiq).NE.0) THEN
303        abort_message =
304     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
305        call abort_gcm(modname,abort_message,1)
306      ENDIF
307
308      zdtvr    = daysec/REAL(day_step)
309        IF(dtvr.NE.zdtvr) THEN
310         WRITE(lunout,*)
311     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
312        ENDIF
313
314C
315C on remet le calendrier à zero si demande
316c
317      IF (start_time /= starttime) then
318        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
319     &,' fichier restart ne correspond pas à celle lue dans le run.def'
320        IF (raz_date == 1) then
321          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
322          start_time = starttime
323        ELSE
324          WRITE(lunout,*)'Je m''arrete'
325          CALL abort
326        ENDIF
327      ENDIF
328      IF (raz_date == 1) THEN
329        annee_ref = anneeref
330        day_ref = dayref
331        day_ini = dayref
332        itau_dyn = 0
333        itau_phy = 0
334        time_0 = 0.
335        write(lunout,*)
336     .   'GCM: On reinitialise a la date lue dans gcm.def'
337      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
338        write(lunout,*)
339     .  'GCM: Attention les dates initiales lues dans le fichier'
340        write(lunout,*)
341     .  ' restart ne correspondent pas a celles lues dans '
342        write(lunout,*)' gcm.def'
343        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
344        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
345        write(lunout,*)' Pas de remise a zero'
346      ENDIF
347
348c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
349c        write(lunout,*)
350c     .  'GCM: Attention les dates initiales lues dans le fichier'
351c        write(lunout,*)
352c     .  ' restart ne correspondent pas a celles lues dans '
353c        write(lunout,*)' gcm.def'
354c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
355c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
356c        if (raz_date .ne. 1) then
357c          write(lunout,*)
358c     .    'GCM: On garde les dates du fichier restart'
359c        else
360c          annee_ref = anneeref
361c          day_ref = dayref
362c          day_ini = dayref
363c          itau_dyn = 0
364c          itau_phy = 0
365c          time_0 = 0.
366c          write(lunout,*)
367c     .   'GCM: On reinitialise a la date lue dans gcm.def'
368c        endif
369c      ELSE
370c        raz_date = 0
371c      endif
372
373#ifdef CPP_IOIPSL
374      mois = 1
375      heure = 0.
376! Ce n'est defini pour l'instant que pour la Terre...
377      if (planet_type.eq.'earth') then
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! A voir pour Titan et Venus
392        jD_ref=0
393        jH_ref=0
394      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
395      write(lunout,*)jD_ref,jH_ref
396      endif ! planet_type
397#else
398! Ehouarn: we still need to define JD_ref and JH_ref
399! and since we don't know how many days there are in a year
400! we set JD_ref to 0 (this should be improved ...)
401      jD_ref=0
402      jH_ref=0
403#endif
404
405      if (iflag_phys.eq.1) then
406      ! these initialisations have already been done (via iniacademic)
407      ! if running in SW or Newtonian mode
408c-----------------------------------------------------------------------
409c   Initialisation des constantes dynamiques :
410c   ------------------------------------------
411        dtvr = zdtvr
412        CALL iniconst
413
414c-----------------------------------------------------------------------
415c   Initialisation de la geometrie :
416c   --------------------------------
417        CALL inigeom
418
419c-----------------------------------------------------------------------
420c   Initialisation du filtre :
421c   --------------------------
422        CALL inifilr
423      endif ! of if (iflag_phys.eq.1)
424c
425c-----------------------------------------------------------------------
426c   Initialisation de la dissipation :
427c   ----------------------------------
428
429      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
430     *                tetagdiv, tetagrot , tetatemp              )
431
432c-----------------------------------------------------------------------
433c   Initialisation de la physique :
434c   -------------------------------
435
436      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
437         latfi(1)=rlatu(1)
438         lonfi(1)=0.
439         zcufi(1) = cu(1)
440         zcvfi(1) = cv(1)
441         DO j=2,jjm
442            DO i=1,iim
443               latfi((j-2)*iim+1+i)= rlatu(j)
444               lonfi((j-2)*iim+1+i)= rlonv(i)
445               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
446               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
447            ENDDO
448         ENDDO
449         latfi(ngridmx)= rlatu(jjp1)
450         lonfi(ngridmx)= 0.
451         zcufi(ngridmx) = cu(ip1jm+1)
452         zcvfi(ngridmx) = cv(ip1jm-iim)
453
454         ! build airefi(), mesh area on physics grid
455         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
456         ! Poles are single points on physics grid
457         airefi(1)=airefi(1)*iim
458         airefi(ngridmx)=airefi(ngridmx)*iim
459
460! Initialisation de la physique: pose probleme quand on tourne
461! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
462! Il faut une cle CPP_PHYS
463#ifdef CPP_PHYS
464         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
465     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
466     &                iflag_phys)
467#endif
468         call_iniphys=.false.
469      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
470
471c  numero de stockage pour les fichiers de redemarrage:
472
473c-----------------------------------------------------------------------
474c   Initialisation des I/O :
475c   ------------------------
476
477
478      day_end = day_ini + nday
479      WRITE(lunout,300)day_ini,day_end
480 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
481
482#ifdef CPP_IOIPSL
483! Ce n'est defini pour l'instant que pour la Terre...
484      if (planet_type.eq.'earth') then
485      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
486      write (lunout,301)jour, mois, an
487      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
488      write (lunout,302)jour, mois, an
489      else
490! A voir pour Titan et Venus
491      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
492      endif ! planet_type
493
494 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
495 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
496#endif
497
498      if (planet_type.eq."mars") then
499! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
500         abort_message = 'dynredem0_mars A FAIRE'
501         call abort_gcm(modname,abort_message,0)
502      else
503        CALL dynredem0("restart.nc", day_end, phis)
504      endif ! of if (planet_type.eq."mars")
505
506      ecripar = .TRUE.
507
508#ifdef CPP_IOIPSL
509      time_step = zdtvr
510      if (ok_dyn_ins) then
511        ! initialize output file for instantaneous outputs
512        ! t_ops = iecri * daysec ! do operations every t_ops
513        t_ops =((1.0*iecri)/day_step) * daysec 
514        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
515        CALL inithist(day_ref,annee_ref,time_step,
516     &              t_ops,t_wrt)
517      endif
518
519      IF (ok_dyn_ave) THEN
520        ! initialize output file for averaged outputs
521        t_ops = iperiod * time_step ! do operations every t_ops
522        t_wrt = periodav * daysec   ! write output every t_wrt
523        CALL initdynav(day_ref,annee_ref,time_step,
524     &       t_ops,t_wrt)
525      END IF
526      dtav = iperiod*dtvr/daysec
527#endif
528! #endif of #ifdef CPP_IOIPSL
529
530c  Choix des frequences de stokage pour le offline
531c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
532c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
533      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
534      istphy=istdyn/iphysiq     
535
536
537c
538c-----------------------------------------------------------------------
539c   Integration temporelle du modele :
540c   ----------------------------------
541
542c       write(78,*) 'ucov',ucov
543c       write(78,*) 'vcov',vcov
544c       write(78,*) 'teta',teta
545c       write(78,*) 'ps',ps
546c       write(78,*) 'q',q
547
548
549      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,
550     .              time_0)
551
552      END
553
Note: See TracBrowser for help on using the repository browser.