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

Last change on this file since 2495 was 2478, checked in by emillour, 4 years ago

Common dynamics
Bug fix for the (rare) case when running without physics (iflag_phys=0),
where the dynamical time steps was not correctly initialized, as it
was overwritten by the value in controle() from start.nc instead of being
computed from daysec and day_step.
EM

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  dtvr=zdtvr
299
300  !
301  ! on remet le calendrier à zero si demande
302  !
303  IF (start_time /= starttime) then
304     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
305     ,' fichier restart ne correspond pas à celle lue dans le run.def'
306     IF (raz_date == 1) then
307        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
308        start_time = starttime
309     ELSE
310        call abort_gcm("gcm", "'Je m''arrete'", 1)
311     ENDIF
312  ENDIF
313  IF (raz_date == 1) THEN
314     annee_ref = anneeref
315     day_ref = dayref
316     day_ini = dayref
317     itau_dyn = 0
318     itau_phy = 0
319     time_0 = 0.
320     write(lunout,*) &
321         'GCM: On reinitialise a la date lue dans gcm.def'
322  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
323     write(lunout,*) &
324        'GCM: Attention les dates initiales lues dans le fichier'
325     write(lunout,*) &
326        ' restart ne correspondent pas a celles lues dans '
327     write(lunout,*)' gcm.def'
328     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
329     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
330     write(lunout,*)' Pas de remise a zero'
331  ENDIF
332
333!      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
334!        write(lunout,*)
335!     .  'GCM: Attention les dates initiales lues dans le fichier'
336!        write(lunout,*)
337!     .  ' restart ne correspondent pas a celles lues dans '
338!        write(lunout,*)' gcm.def'
339!        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
340!        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
341!        if (raz_date .ne. 1) then
342!          write(lunout,*)
343!     .    'GCM: On garde les dates du fichier restart'
344!        else
345!          annee_ref = anneeref
346!          day_ref = dayref
347!          day_ini = dayref
348!          itau_dyn = 0
349!          itau_phy = 0
350!          time_0 = 0.
351!          write(lunout,*)
352!     .   'GCM: On reinitialise a la date lue dans gcm.def'
353!        endif
354!      ELSE
355!        raz_date = 0
356!      endif
357
358#ifdef CPP_IOIPSL
359  mois = 1
360  heure = 0.
361! Ce n'est defini pour l'instant que pour la Terre...
362  if (planet_type.eq.'earth') then
363    call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
364    jH_ref = jD_ref - int(jD_ref)
365    jD_ref = int(jD_ref)
366
367    call ioconf_startdate(INT(jD_ref), jH_ref)
368
369    write(lunout,*)'DEBUG'
370    write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
371    write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
372    call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
373    write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
374    write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
375!  else if (planet_type.eq.'titan') then
376!    jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
377  else
378! A voir pour Titan et Venus
379    jD_ref=0
380    jH_ref=0
381    write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
382    write(lunout,*)jD_ref,jH_ref
383  endif ! planet_type
384#else
385  ! Ehouarn: we still need to define JD_ref and JH_ref
386  ! and since we don't know how many days there are in a year
387  ! we set JD_ref to 0 (this should be improved ...)
388  jD_ref=0
389  jH_ref=0
390#endif
391
392  if (iflag_phys.eq.1) then
393     ! these initialisations have already been done (via iniacademic)
394     ! if running in SW or Newtonian mode
395     !-----------------------------------------------------------------------
396     !   Initialisation des constantes dynamiques :
397     !   ------------------------------------------
398     dtvr = zdtvr
399     CALL iniconst
400
401     !-----------------------------------------------------------------------
402     !   Initialisation de la geometrie :
403     !   --------------------------------
404     CALL inigeom
405
406     !-----------------------------------------------------------------------
407     !   Initialisation du filtre :
408     !   --------------------------
409     CALL inifilr
410  endif ! of if (iflag_phys.eq.1)
411  !
412  !-----------------------------------------------------------------------
413  !   Initialisation de la dissipation :
414  !   ----------------------------------
415
416  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
417                  tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
418
419  !-----------------------------------------------------------------------
420  !   Initialisation des I/O :
421  !   ------------------------
422
423  if (nday>=0) then ! standard case
424     day_end=day_ini+nday
425  else ! special case when nday <0, run -nday dynamical steps
426     day_end=day_ini-nday/day_step
427  endif
428  if (less1day) then
429     day_end=day_ini+floor(time_0+fractday)
430  endif
431  if (ndynstep.gt.0) then
432     day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
433  endif
434     
435  WRITE(lunout,'(a,i7,a,i7)') &
436               "run from day ",day_ini,"  to day",day_end
437
438#ifdef CPP_IOIPSL
439  ! Ce n'est defini pour l'instant que pour la Terre...
440  if (planet_type.eq.'earth') then
441    call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
442    write (lunout,301)jour, mois, an
443    call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
444    write (lunout,302)jour, mois, an
445  else
446  ! A voir pour Titan et Venus
447    write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
448  endif ! planet_type
449301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
450302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
451#endif
452
453
454  !-----------------------------------------------------------------------
455  !   Initialisation de la physique :
456  !   -------------------------------
457
458  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
459     ! Physics:
460#ifdef CPP_PHYS
461     CALL iniphysiq(iim,jjm,llm, &
462          (jjm-1)*iim+2,comm_lmdz, &
463          daysec,day_ini,dtphys/nsplit_phys, &
464          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
465          iflag_phys)
466#endif
467  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
468
469
470  if (planet_type=="mars") then
471    ! For Mars we transmit day_ini
472    CALL dynredem0("restart.nc", day_ini, phis)
473  else
474    CALL dynredem0("restart.nc", day_end, phis)
475  endif
476  ecripar = .TRUE.
477
478#ifdef CPP_IOIPSL
479  time_step = zdtvr
480  if (ok_dyn_ins) then
481    ! initialize output file for instantaneous outputs
482    ! t_ops = iecri * daysec ! do operations every t_ops
483    t_ops =((1.0*iecri)/day_step) * daysec 
484    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
485    CALL inithist(day_ref,annee_ref,time_step, &
486                   t_ops,t_wrt)
487  endif
488
489  IF (ok_dyn_ave) THEN
490    ! initialize output file for averaged outputs
491    t_ops = iperiod * time_step ! do operations every t_ops
492    t_wrt = periodav * daysec   ! write output every t_wrt
493    CALL initdynav(day_ref,annee_ref,time_step, &
494            t_ops,t_wrt)
495  END IF
496  dtav = iperiod*dtvr/daysec
497#endif
498! #endif of #ifdef CPP_IOIPSL
499
500  !  Choix des frequences de stokage pour le offline
501  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
502  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
503  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
504  istphy=istdyn/iphysiq     
505
506
507  !
508  !-----------------------------------------------------------------------
509  !   Integration temporelle du modele :
510  !   ----------------------------------
511
512  !       write(78,*) 'ucov',ucov
513  !       write(78,*) 'vcov',vcov
514  !       write(78,*) 'teta',teta
515  !       write(78,*) 'ps',ps
516  !       write(78,*) 'q',q
517
518
519  CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
520
521END PROGRAM gcm
522
Note: See TracBrowser for help on using the repository browser.