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

Last change on this file since 3785 was 3615, checked in by emillour, 4 months ago

Venus PCM: Corrections to enable 1+1=2

  • store correctly the time_of_day in restart.nc to enable proper restart
  • enforce recomputation of CP in the physics at all time steps (otherwise when without thermosphere the value was only computed at first step and kept unchanged afterwards).

EM

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