source: LMDZ5/trunk/libf/dyn3d/gcm.F90 @ 2347

Last change on this file since 2347 was 2347, checked in by Ehouarn Millour, 9 years ago

Make iniphysiq a module.
Fix call to iniphysiq in lmdz1d (missing arguments and arrays of wrong sizes).
EM

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 15.0 KB
Line 
1!
2! $Id: gcm.F90 2347 2015-08-24 07:26:48Z 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
24
25#ifdef INCA
26  ! Only INCA needs these informations (from the Earth's physics)
27  USE indice_sol_mod
28  USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
29#endif
30
31!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
32  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
33  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
34  ! dynamique -> physique pour l'initialisation
35#ifdef CPP_PHYS
36  USE iniphysiq_mod, ONLY: iniphysiq
37#endif
38!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
39
40  IMPLICIT NONE
41
42  !      ......   Version  du 10/01/98    ..........
43
44  !             avec  coordonnees  verticales hybrides
45  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
46
47  !=======================================================================
48  !
49  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
50  !   -------
51  !
52  !   Objet:
53  !   ------
54  !
55  !   GCM LMD nouvelle grille
56  !
57  !=======================================================================
58  !
59  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
60  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
61  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
62  !  ... Possibilite de choisir le schema pour l'advection de
63  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
64  !
65  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
66  !      Pour Van-Leer iadv=10
67  !
68  !-----------------------------------------------------------------------
69  !   Declarations:
70  !   -------------
71
72  include "dimensions.h"
73  include "paramet.h"
74  include "comconst.h"
75  include "comdissnew.h"
76  include "comvert.h"
77  include "comgeom.h"
78  include "logic.h"
79  include "temps.h"
80!!!!!!!!!!!include "control.h"
81  include "ener.h"
82  include "description.h"
83  include "serre.h"
84  !include "com_io_dyn.h"
85  include "iniprint.h"
86  include "tracstoke.h"
87#ifdef INCA
88  ! Only INCA needs these informations (from the Earth's physics)
89  !include "indicesol.h"
90#endif
91
92  REAL zdtvr
93
94  !   variables dynamiques
95  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
96  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
97  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
98  REAL ps(ip1jmp1)                       ! pression  au sol
99  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
100  REAL masse(ip1jmp1,llm)                ! masse d'air
101  REAL phis(ip1jmp1)                     ! geopotentiel au sol
102  REAL phi(ip1jmp1,llm)                  ! geopotentiel
103  REAL w(ip1jmp1,llm)                    ! vitesse verticale
104
105  ! variables dynamiques intermediaire pour le transport
106
107  !   variables pour le fichier histoire
108  REAL dtav      ! intervalle de temps elementaire
109
110  REAL time_0
111
112  LOGICAL lafin
113  INTEGER ij,iq,l,i,j
114
115
116  real time_step, t_wrt, t_ops
117
118  LOGICAL first
119
120  !      LOGICAL call_iniphys
121  !      data call_iniphys/.true./
122
123  !+jld variables test conservation energie
124  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
125  !     Tendance de la temp. potentiel d (theta)/ d t due a la
126  !     tansformation d'energie cinetique en energie thermique
127  !     cree par la dissipation
128  REAL dhecdt(ip1jmp1,llm)
129  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
130  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
131  CHARACTER (len=15) :: ztit
132  !-jld
133
134
135  character (len=80) :: dynhist_file, dynhistave_file
136  character (len=20) :: modname
137  character (len=80) :: abort_message
138  ! locales pour gestion du temps
139  INTEGER :: an, mois, jour
140  REAL :: heure
141  logical use_filtre_fft
142
143  !-----------------------------------------------------------------------
144  !   Initialisations:
145  !   ----------------
146
147  abort_message = 'last timestep reached'
148  modname = 'gcm'
149  descript = 'Run GCM LMDZ'
150  lafin    = .FALSE.
151  dynhist_file = 'dyn_hist.nc'
152  dynhistave_file = 'dyn_hist_ave.nc'
153
154
155
156  !----------------------------------------------------------------------
157  !  lecture des fichiers gcm.def ou run.def
158  !  ---------------------------------------
159  !
160  CALL conf_gcm( 99, .TRUE.)
161
162  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
163       "iphysiq must be a multiple of iperiod", 1)
164
165  use_filtre_fft=.FALSE.
166  CALL getin('use_filtre_fft',use_filtre_fft)
167  IF (use_filtre_fft) call abort_gcm('FFT filter is not available in the ' &
168          // 'sequential version of the dynamics.', 1)
169
170!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
171  ! Initialisation de XIOS
172!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
173
174#ifdef CPP_XIOS
175  CALL wxios_init("LMDZ")
176#endif
177
178
179!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
180  ! FH 2008/05/02
181  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
182  ! dynamique -> physique pour l'initialisation
183#ifdef CPP_PHYS
184  CALL Init_Phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
185  !      call InitComgeomphy ! now done in iniphysiq
186#endif
187!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
188  !-----------------------------------------------------------------------
189  !   Choix du calendrier
190  !   -------------------
191
192  !      calend = 'earth_365d'
193
194#ifdef CPP_IOIPSL
195  if (calend == 'earth_360d') then
196     call ioconf_calendar('360d')
197     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
198  else if (calend == 'earth_365d') then
199     call ioconf_calendar('noleap')
200     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
201  else if (calend == 'gregorian') then
202     call ioconf_calendar('gregorian')
203     write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
204  else
205     abort_message = 'Mauvais choix de calendrier'
206     call abort_gcm(modname,abort_message,1)
207  endif
208#endif
209  !-----------------------------------------------------------------------
210
211  IF (type_trac == 'inca') THEN
212#ifdef INCA
213     call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,  &
214          nbsrf, is_oce,is_sic,is_ter,is_lic)
215     call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
216#endif
217  END IF
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     !        if (planet_type.eq."earth") then
249     ! Load an Earth-format start file
250     CALL dynetat0("start.nc",vcov,ucov, &
251          teta,q,masse,ps,phis, time_0)
252     !        endif ! of if (planet_type.eq."earth")
253
254     !       write(73,*) 'ucov',ucov
255     !       write(74,*) 'vcov',vcov
256     !       write(75,*) 'teta',teta
257     !       write(76,*) 'ps',ps
258     !       write(77,*) 'q',q
259
260  endif ! of if (read_start)
261
262  IF (type_trac == 'inca') THEN
263#ifdef INCA
264     call init_inca_dim(klon,llm,iim,jjm, &
265          rlonu,rlatu,rlonv,rlatv)
266#endif
267  END IF
268
269
270  ! le cas echeant, creation d un etat initial
271  IF (prt_level > 9) WRITE(lunout,*) &
272       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
273  if (.not.read_start) then
274     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
275  endif
276
277
278  !-----------------------------------------------------------------------
279  !   Lecture des parametres de controle pour la simulation :
280  !   -------------------------------------------------------
281  !  on recalcule eventuellement le pas de temps
282
283  IF(MOD(day_step,iperiod).NE.0) THEN
284     abort_message =  &
285          'Il faut choisir un nb de pas par jour multiple de iperiod'
286     call abort_gcm(modname,abort_message,1)
287  ENDIF
288
289  IF(MOD(day_step,iphysiq).NE.0) THEN
290     abort_message =  &
291          'Il faut choisir un nb de pas par jour multiple de iphysiq'
292     call abort_gcm(modname,abort_message,1)
293  ENDIF
294
295  zdtvr    = daysec/REAL(day_step)
296  IF(dtvr.NE.zdtvr) THEN
297     WRITE(lunout,*) &
298          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
299  ENDIF
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  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
375  ! Ehouarn: we still need to define JD_ref and JH_ref
376  ! and since we don't know how many days there are in a year
377  ! we set JD_ref to 0 (this should be improved ...)
378  jD_ref=0
379  jH_ref=0
380#endif
381
382
383  if (iflag_phys.eq.1) then
384     ! these initialisations have already been done (via iniacademic)
385     ! if running in SW or Newtonian mode
386     !-----------------------------------------------------------------------
387     !   Initialisation des constantes dynamiques :
388     !   ------------------------------------------
389     dtvr = zdtvr
390     CALL iniconst
391
392     !-----------------------------------------------------------------------
393     !   Initialisation de la geometrie :
394     !   --------------------------------
395     CALL inigeom
396
397     !-----------------------------------------------------------------------
398     !   Initialisation du filtre :
399     !   --------------------------
400     CALL inifilr
401  endif ! of if (iflag_phys.eq.1)
402  !
403  !-----------------------------------------------------------------------
404  !   Initialisation de la dissipation :
405  !   ----------------------------------
406
407  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
408       tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
409
410  !-----------------------------------------------------------------------
411  !   Initialisation de la physique :
412  !   -------------------------------
413
414  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
415     ! Physics:
416#ifdef CPP_PHYS
417     CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
418          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
419          iflag_phys)
420#endif
421  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
422
423  !  numero de stockage pour les fichiers de redemarrage:
424
425  !-----------------------------------------------------------------------
426  !   Initialisation des I/O :
427  !   ------------------------
428
429
430  if (nday>=0) then
431     day_end = day_ini + nday
432  else
433     day_end = day_ini - nday/day_step
434  endif
435  WRITE(lunout,300)day_ini,day_end
436300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
437
438#ifdef CPP_IOIPSL
439  call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
440  write (lunout,301)jour, mois, an
441  call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
442  write (lunout,302)jour, mois, an
443301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
444302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
445#endif
446
447  !      if (planet_type.eq."earth") then
448  ! Write an Earth-format restart file
449
450  CALL dynredem0("restart.nc", day_end, phis)
451  !      endif
452
453  ecripar = .TRUE.
454
455#ifdef CPP_IOIPSL
456  time_step = zdtvr
457  if (ok_dyn_ins) then
458     ! initialize output file for instantaneous outputs
459     ! t_ops = iecri * daysec ! do operations every t_ops
460     t_ops =((1.0*iecri)/day_step) * daysec 
461     t_wrt = daysec ! iecri * daysec ! write output every t_wrt
462     CALL inithist(day_ref,annee_ref,time_step, &
463          t_ops,t_wrt)
464  endif
465
466  IF (ok_dyn_ave) THEN
467     ! initialize output file for averaged outputs
468     t_ops = iperiod * time_step ! do operations every t_ops
469     t_wrt = periodav * daysec   ! write output every t_wrt
470     CALL initdynav(day_ref,annee_ref,time_step, &
471          t_ops,t_wrt)
472  END IF
473  dtav = iperiod*dtvr/daysec
474#endif
475  ! #endif of #ifdef CPP_IOIPSL
476
477  !  Choix des frequences de stokage pour le offline
478  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
479  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
480  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
481  istphy=istdyn/iphysiq     
482
483
484  !
485  !-----------------------------------------------------------------------
486  !   Integration temporelle du modele :
487  !   ----------------------------------
488
489  !       write(78,*) 'ucov',ucov
490  !       write(78,*) 'vcov',vcov
491  !       write(78,*) 'teta',teta
492  !       write(78,*) 'ps',ps
493  !       write(78,*) 'q',q
494
495
496  CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
497
498END PROGRAM gcm
Note: See TracBrowser for help on using the repository browser.