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

Last change on this file since 1235 was 1107, checked in by emillour, 11 years ago

Common dynamics: Updates and modifications to enable running Mars physics with

LMDZ.COMMON dynamics:

  • For compilation: adapted makelmdz, create_make_gcm and makelmdz_fcm, bld.cfg to compile aeronomy routines in "aerono$physique" if it exists, and added "-P -traditional" preprocessing flags in "arch-linux-ifort*"
  • Added function "cbrt.F" (cubic root) in 'bibio'
  • Adapted the reading/writing of dynamics (re)start.nc files for Mars. The main issue is that different information (on time, reference and current) is stored and used differently, hence a few if (planet_type =="mars") here and there. Moreover in the martian case there is the possibility to store fields over multiple times. Some Mars-specific variables (ecritphy,ecritstart,timestart) added in control_mod.F and (hour_ini) in temps.h

EM

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