source: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F90 @ 1672

Last change on this file since 1672 was 1672, checked in by jvatant, 8 years ago

Re-plug chemistry of old Titan
+ minor correction on reading Titan's startfiles and haze routine
JVO

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