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

Last change on this file since 3553 was 3316, checked in by jbclement, 8 months ago

Mars PCM:
Reversion of r3305 and r3307.
JBC

File size: 17.0 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  logical :: present
140
141!-----------------------------------------------------------------------
142!   Initialisations:
143!   ----------------
144
145  abort_message = 'last timestep reached'
146  modname = 'gcm'
147  lafin    = .FALSE.
148  dynhist_file = 'dyn_hist.nc'
149  dynhistave_file = 'dyn_hist_ave.nc'
150
151
152
153!----------------------------------------------------------------------
154!  Load flags and options from run.def
155!  ---------------------------------------
156!
157! Start by checking that there indeed is a run.def file
158  INQUIRE(file="run.def", EXIST=present)
159  IF(.not. present) then
160    CALL abort_gcm("gcm", "file run.def not found, Aborting!",1)
161  ENDIF
162
163  CALL conf_gcm( 99, .TRUE. )
164  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
165       "iphysiq must be a multiple of iperiod", 1)
166
167  use_filtre_fft=.FALSE.
168  CALL getin('use_filtre_fft',use_filtre_fft)
169  IF (use_filtre_fft) call abort_gcm("gcm",'FFT filter is not available in the ' &
170          // 'sequential version of the dynamics.', 1)
171
172!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173! Initialisation de XIOS
174!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
175
176#ifdef CPP_XIOS
177  CALL wxios_init("LMDZ")
178#endif
179
180
181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182! FH 2008/05/02
183! A nettoyer. On ne veut qu'une ou deux routines d'interface
184! dynamique -> physique pour l'initialisation
185!#ifdef CPP_PHYS
186!  CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
187!!      call initcomgeomphy ! now done in iniphysiq
188!#endif
189!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190!
191! Initialisations pour Cp(T) Venus
192  call ini_cpdet
193!
194!-----------------------------------------------------------------------
195!   Choix du calendrier
196!   -------------------
197
198!      calend = 'earth_365d'
199
200#ifdef CPP_IOIPSL
201  if (calend == 'earth_360d') then
202    call ioconf_calendar('360d')
203    write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
204  else if (calend == 'earth_365d') then
205    call ioconf_calendar('noleap')
206    write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
207  else if (calend == 'earth_366d') then
208    call ioconf_calendar('gregorian')
209    write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
210  else if (calend == 'titan') then
211!        call ioconf_calendar('titan')
212    write(lunout,*)'CALENDRIER CHOISI: Titan'
213    abort_message = 'A FAIRE...'
214    call abort_gcm(modname,abort_message,1)
215  else if (calend == 'venus') then
216!        call ioconf_calendar('venus')
217    write(lunout,*)'CALENDRIER CHOISI: Venus'
218    abort_message = 'A FAIRE...'
219    call abort_gcm(modname,abort_message,1)
220  else
221    abort_message = 'Mauvais choix de calendrier'
222    call abort_gcm(modname,abort_message,1)
223  endif
224#endif
225!-----------------------------------------------------------------------
226  !
227  !
228  !------------------------------------
229  !   Initialisation partie parallele
230  !------------------------------------
231
232  !
233  !
234  !-----------------------------------------------------------------------
235  !   Initialisation des traceurs
236  !   ---------------------------
237  !  Choix du nombre de traceurs et du schema pour l'advection
238  !  dans fichier traceur.def, par default ou via INCA
239  call infotrac_init
240
241  ! Allocation de la tableau q : champs advectes   
242  allocate(q(ip1jmp1,llm,nqtot))
243
244  !-----------------------------------------------------------------------
245  !   Lecture de l'etat initial :
246  !   ---------------------------
247
248  !  lecture du fichier start.nc
249  if (read_start) then
250     ! we still need to run iniacademic to initialize some
251     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
252     if (iflag_phys.ne.1) then
253        CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
254     endif
255
256     CALL dynetat0("start.nc",vcov,ucov, &
257                    teta,q,masse,ps,phis, time_0)
258       
259     ! Load relaxation fields (simple nudging). AS 09/2013
260     ! ---------------------------------------------------
261     if (planet_type.eq."generic") then
262       if (ok_guide) then
263         CALL relaxetat0("relax.nc")
264       endif
265     endif
266 
267     !       write(73,*) 'ucov',ucov
268     !       write(74,*) 'vcov',vcov
269     !       write(75,*) 'teta',teta
270     !       write(76,*) 'ps',ps
271     !       write(77,*) 'q',q
272
273  endif ! of if (read_start)
274
275
276  ! le cas echeant, creation d un etat initial
277  IF (prt_level > 9) WRITE(lunout,*) &
278       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
279  if (.not.read_start) then
280     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
281  endif
282
283
284  !-----------------------------------------------------------------------
285  !   Lecture des parametres de controle pour la simulation :
286  !   -------------------------------------------------------
287  !  on recalcule eventuellement le pas de temps
288
289  IF(MOD(day_step,iperiod).NE.0) THEN
290     abort_message = &
291       'Il faut choisir un nb de pas par jour multiple de iperiod'
292     call abort_gcm(modname,abort_message,1)
293  ENDIF
294
295  IF(MOD(day_step,iphysiq).NE.0) THEN
296     abort_message = &
297       'Il faut choisir un nb de pas par jour multiple de iphysiq'
298     call abort_gcm(modname,abort_message,1)
299  ENDIF
300
301  zdtvr    = daysec/REAL(day_step)
302  IF(dtvr.NE.zdtvr) THEN
303     WRITE(lunout,*) &
304          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
305  ENDIF
306  dtvr=zdtvr
307
308  !
309  ! on remet le calendrier a zero si demande
310  !
311  IF (start_time /= starttime) then
312     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
313     ,' fichier restart ne correspond pas a celle lue dans le run.def'
314     IF (raz_date == 1) then
315        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
316        start_time = starttime
317     ELSE
318        call abort_gcm("gcm", "'Je m''arrete'", 1)
319     ENDIF
320  ENDIF
321  IF (raz_date == 1) THEN
322     annee_ref = anneeref
323     day_ref = dayref
324     day_ini = dayref
325     itau_dyn = 0
326     itau_phy = 0
327     time_0 = 0.
328     write(lunout,*) &
329         'GCM: On reinitialise a la date lue dans gcm.def'
330  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
331     write(lunout,*) &
332        'GCM: Attention les dates initiales lues dans le fichier'
333     write(lunout,*) &
334        ' restart ne correspondent pas a celles lues dans '
335     write(lunout,*)' gcm.def'
336     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
337     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
338     write(lunout,*)' Pas de remise a zero'
339  ENDIF
340
341!      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
342!        write(lunout,*)
343!     .  'GCM: Attention les dates initiales lues dans le fichier'
344!        write(lunout,*)
345!     .  ' restart ne correspondent pas a celles lues dans '
346!        write(lunout,*)' gcm.def'
347!        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
348!        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
349!        if (raz_date .ne. 1) then
350!          write(lunout,*)
351!     .    'GCM: On garde les dates du fichier restart'
352!        else
353!          annee_ref = anneeref
354!          day_ref = dayref
355!          day_ini = dayref
356!          itau_dyn = 0
357!          itau_phy = 0
358!          time_0 = 0.
359!          write(lunout,*)
360!     .   'GCM: On reinitialise a la date lue dans gcm.def'
361!        endif
362!      ELSE
363!        raz_date = 0
364!      endif
365
366#ifdef CPP_IOIPSL
367  mois = 1
368  heure = 0.
369! Ce n'est defini pour l'instant que pour la Terre...
370  if (planet_type.eq.'earth') then
371    call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
372    jH_ref = jD_ref - int(jD_ref)
373    jD_ref = int(jD_ref)
374
375    call ioconf_startdate(INT(jD_ref), jH_ref)
376
377    write(lunout,*)'DEBUG'
378    write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
379    write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
380    call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
381    write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
382    write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
383!  else if (planet_type.eq.'titan') then
384!    jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
385  else
386! A voir pour Titan et Venus
387    jD_ref=0
388    jH_ref=0
389    write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
390    write(lunout,*)jD_ref,jH_ref
391  endif ! planet_type
392#else
393  ! Ehouarn: we still need to define JD_ref and JH_ref
394  ! and since we don't know how many days there are in a year
395  ! we set JD_ref to 0 (this should be improved ...)
396  jD_ref=0
397  jH_ref=0
398#endif
399
400  if (iflag_phys.eq.1) then
401     ! these initialisations have already been done (via iniacademic)
402     ! if running in SW or Newtonian mode
403     !-----------------------------------------------------------------------
404     !   Initialisation des constantes dynamiques :
405     !   ------------------------------------------
406     dtvr = zdtvr
407     CALL iniconst
408
409     !-----------------------------------------------------------------------
410     !   Initialisation de la geometrie :
411     !   --------------------------------
412     CALL inigeom
413
414     !-----------------------------------------------------------------------
415     !   Initialisation du filtre :
416     !   --------------------------
417     CALL inifilr
418  endif ! of if (iflag_phys.eq.1)
419  !
420  !-----------------------------------------------------------------------
421  !   Initialisation de la dissipation :
422  !   ----------------------------------
423
424  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
425                  tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
426
427  !-----------------------------------------------------------------------
428  !   Initialisation des I/O :
429  !   ------------------------
430
431  if (nday>=0) then ! standard case
432     day_end=day_ini+nday
433  else ! special case when nday <0, run -nday dynamical steps
434     day_end=day_ini-nday/day_step
435  endif
436  if (less1day) then
437     day_end=day_ini+floor(time_0+fractday)
438  endif
439  if (ndynstep.gt.0) then
440     day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
441  endif
442     
443  WRITE(lunout,'(a,i7,a,i7)') &
444               "run from day ",day_ini,"  to day",day_end
445
446#ifdef CPP_IOIPSL
447  ! Ce n'est defini pour l'instant que pour la Terre...
448  if (planet_type.eq.'earth') then
449    call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
450    write (lunout,301)jour, mois, an
451    call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
452    write (lunout,302)jour, mois, an
453  else
454  ! A voir pour Titan et Venus
455    write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
456  endif ! planet_type
457301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
458302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
459#endif
460
461
462  !-----------------------------------------------------------------------
463  !   Initialisation de la physique :
464  !   -------------------------------
465
466  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
467     ! Physics:
468#ifdef CPP_PHYS
469     CALL iniphysiq(iim,jjm,llm,                        &
470                    (jjm-1)*iim+2,comm_lmdz,            &
471                    daysec,day_ini,dtphys/nsplit_phys,  &
472                    rlatu,rlatv,rlonu,rlonv,aire,cu,cv, &
473                    rad,g,r,cpp,iflag_phys)
474#endif
475  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
476
477
478! If we want to save multiple restart at different time (ecritstart.GT.0)
479! We will write in restart the initial day of the simulation and put all
480! the different time in the Time variable
481! However, if we write only one restart, we save the last day of the
482! simulation and put the remaining time (which is 0 if we do fulls days)
483! in the Time variable
484
485  if (planet_type=="mars") then
486    if (ecritstart.GT.0) then
487    ! For Mars we transmit day_ini
488      CALL dynredem0("restart.nc", day_ini, phis)
489    else
490      CALL dynredem0("restart.nc", day_end, phis)
491    endif
492  else
493    CALL dynredem0("restart.nc", day_end, phis)
494  endif
495  ecripar = .TRUE.
496
497#ifdef CPP_IOIPSL
498  time_step = zdtvr
499  if (ok_dyn_ins) then
500    ! initialize output file for instantaneous outputs
501    ! t_ops = iecri * daysec ! do operations every t_ops
502    t_ops =((1.0*iecri)/day_step) * daysec 
503    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
504    CALL inithist(day_ref,annee_ref,time_step,t_ops,t_wrt)
505  endif
506
507  IF (ok_dyn_ave) THEN
508    ! initialize output file for averaged outputs
509    t_ops = iperiod * time_step ! do operations every t_ops
510    t_wrt = periodav * daysec   ! write output every t_wrt
511    CALL initdynav(day_ref,annee_ref,time_step,t_ops,t_wrt)
512  END IF
513  dtav = iperiod*dtvr/daysec
514#endif
515! #endif of #ifdef CPP_IOIPSL
516
517  !  Choix des frequences de stokage pour le offline
518  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
519  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
520  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
521  istphy=istdyn/iphysiq     
522
523
524  !
525  !-----------------------------------------------------------------------
526  !   Integration temporelle du modele :
527  !   ----------------------------------
528
529  !       write(78,*) 'ucov',ucov
530  !       write(78,*) 'vcov',vcov
531  !       write(78,*) 'teta',teta
532  !       write(78,*) 'ps',ps
533  !       write(78,*) 'q',q
534
535
536  CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
537
538END PROGRAM gcm
539
Note: See TracBrowser for help on using the repository browser.