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

Last change on this file since 2513 was 2507, checked in by romain.vande, 4 years ago

For LMDZ MARS: Update of day_ini, time and hour_ini in restart and retstartfi.
Hour_ini is obsolete. If we write one restart file: day_ini is the last day
of the simulation and the remaining time is in Time (often=0), if we write
multiple restart nothing changes

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