source: dynamico_lmdz/aquaplanet/LMDZ5/libf/dyn3d/gcm.F90 @ 3831

Last change on this file since 3831 was 3831, checked in by ymipsl, 10 years ago

module reorganisation for a cleaner dyn-phys interface
YM

File size: 15.2 KB
Line 
1!
2! $Id: gcm.F90 2247 2015-03-25 17:24:15Z lguez $
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  USE mod_const_mpi, ONLY: COMM_LMDZ
25
26#ifdef INCA
27  ! Only INCA needs these informations (from the Earth's physics)
28  USE indice_sol_mod
29  USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
30#endif
31
32!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33  ! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
34  ! A nettoyer. On ne veut qu'une ou deux routines d'interface
35  ! dynamique -> physique pour l'initialisation
36#ifdef CPP_PHYS
37  !      USE dimphy
38  !      USE comgeomphy
39#endif
40!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
41
42  IMPLICIT NONE
43
44  !      ......   Version  du 10/01/98    ..........
45
46  !             avec  coordonnees  verticales hybrides
47  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
48
49  !=======================================================================
50  !
51  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
52  !   -------
53  !
54  !   Objet:
55  !   ------
56  !
57  !   GCM LMD nouvelle grille
58  !
59  !=======================================================================
60  !
61  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
62  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
63  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
64  !  ... Possibilite de choisir le schema pour l'advection de
65  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
66  !
67  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
68  !      Pour Van-Leer iadv=10
69  !
70  !-----------------------------------------------------------------------
71  !   Declarations:
72  !   -------------
73
74  include "dimensions.h"
75  include "paramet.h"
76  include "comconst.h"
77  include "comdissnew.h"
78  include "comvert.h"
79  include "comgeom.h"
80  include "logic.h"
81  include "temps.h"
82!!!!!!!!!!!include "control.h"
83  include "ener.h"
84  include "description.h"
85  include "serre.h"
86  !include "com_io_dyn.h"
87  include "iniprint.h"
88  include "tracstoke.h"
89#ifdef INCA
90  ! Only INCA needs these informations (from the Earth's physics)
91  !include "indicesol.h"
92#endif
93
94  REAL zdtvr
95
96  !   variables dynamiques
97  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
98  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
99  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
100  REAL ps(ip1jmp1)                       ! pression  au sol
101!  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
102  REAL masse(ip1jmp1,llm)                ! masse d'air
103  REAL phis(ip1jmp1)                     ! geopotentiel au sol
104!  REAL phi(ip1jmp1,llm)                  ! geopotentiel
105!  REAL w(ip1jmp1,llm)                    ! vitesse verticale
106
107  ! variables dynamiques intermediaire pour le transport
108
109  !   variables pour le fichier histoire
110  REAL dtav      ! intervalle de temps elementaire
111
112  REAL time_0
113
114  LOGICAL lafin
115  INTEGER ij,iq,l,i,j
116
117
118  real time_step, t_wrt, t_ops
119
120  LOGICAL first
121
122  !      LOGICAL call_iniphys
123  !      data call_iniphys/.true./
124
125  !+jld variables test conservation energie
126  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
127  !     Tendance de la temp. potentiel d (theta)/ d t due a la
128  !     tansformation d'energie cinetique en energie thermique
129  !     cree par la dissipation
130  REAL dhecdt(ip1jmp1,llm)
131  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
132  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
133  CHARACTER (len=15) :: ztit
134  !-jld
135
136
137  character (len=80) :: dynhist_file, dynhistave_file
138  character (len=20) :: modname
139  character (len=80) :: abort_message
140  ! locales pour gestion du temps
141  INTEGER :: an, mois, jour
142  REAL :: heure
143  logical use_filtre_fft
144
145  !-----------------------------------------------------------------------
146  !   Initialisations:
147  !   ----------------
148
149  abort_message = 'last timestep reached'
150  modname = 'gcm'
151  descript = 'Run GCM LMDZ'
152  lafin    = .FALSE.
153  dynhist_file = 'dyn_hist.nc'
154  dynhistave_file = 'dyn_hist_ave.nc'
155
156
157
158  !----------------------------------------------------------------------
159  !  lecture des fichiers gcm.def ou run.def
160  !  ---------------------------------------
161  !
162  CALL conf_gcm( 99, .TRUE.)
163
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('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/),COMM_LMDZ)
187!  !      call InitComgeomphy ! now done in iniphysiq
188!#endif
189!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
190  !-----------------------------------------------------------------------
191  !   Choix du calendrier
192  !   -------------------
193
194  !      calend = 'earth_365d'
195
196#ifdef CPP_IOIPSL
197  if (calend == 'earth_360d') then
198     call ioconf_calendar('360d')
199     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
200  else if (calend == 'earth_365d') then
201     call ioconf_calendar('noleap')
202     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
203  else if (calend == 'gregorian') then
204     call ioconf_calendar('gregorian')
205     write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
206  else
207     abort_message = 'Mauvais choix de calendrier'
208     call abort_gcm(modname,abort_message,1)
209  endif
210#endif
211  !-----------------------------------------------------------------------
212
213  IF (type_trac == 'inca') THEN
214#ifdef INCA
215     call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,  &
216          nbsrf, is_oce,is_sic,is_ter,is_lic)
217     call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
218#endif
219  END IF
220  !
221  !
222  !------------------------------------
223  !   Initialisation partie parallele
224  !------------------------------------
225
226  !
227  !
228  !-----------------------------------------------------------------------
229  !   Initialisation des traceurs
230  !   ---------------------------
231  !  Choix du nombre de traceurs et du schema pour l'advection
232  !  dans fichier traceur.def, par default ou via INCA
233  call infotrac_init
234
235  ! Allocation de la tableau q : champs advectes   
236  allocate(q(ip1jmp1,llm,nqtot))
237
238  !-----------------------------------------------------------------------
239  !   Lecture de l'etat initial :
240  !   ---------------------------
241
242  !  lecture du fichier start.nc
243  if (read_start) then
244     ! we still need to run iniacademic to initialize some
245     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
246     if (iflag_phys.ne.1) then
247        CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
248     endif
249
250     !        if (planet_type.eq."earth") then
251     ! Load an Earth-format start file
252     CALL dynetat0("start.nc",vcov,ucov, &
253          teta,q,masse,ps,phis, time_0)
254     !        endif ! of if (planet_type.eq."earth")
255
256     !       write(73,*) 'ucov',ucov
257     !       write(74,*) 'vcov',vcov
258     !       write(75,*) 'teta',teta
259     !       write(76,*) 'ps',ps
260     !       write(77,*) 'q',q
261
262  endif ! of if (read_start)
263
264  IF (type_trac == 'inca') THEN
265#ifdef INCA
266     call init_inca_dim(klon,llm,iim,jjm, &
267          rlonu,rlatu,rlonv,rlatv)
268#endif
269  END IF
270
271
272  ! le cas echeant, creation d un etat initial
273  IF (prt_level > 9) WRITE(lunout,*) &
274       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
275  if (.not.read_start) then
276     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
277  endif
278
279
280  !-----------------------------------------------------------------------
281  !   Lecture des parametres de controle pour la simulation :
282  !   -------------------------------------------------------
283  !  on recalcule eventuellement le pas de temps
284
285  IF(MOD(day_step,iperiod).NE.0) THEN
286     abort_message =  &
287          'Il faut choisir un nb de pas par jour multiple de iperiod'
288     call abort_gcm(modname,abort_message,1)
289  ENDIF
290
291  IF(MOD(day_step,iphysiq).NE.0) THEN
292     abort_message =  &
293          'Il faut choisir un nb de pas par jour multiple de iphysiq'
294     call abort_gcm(modname,abort_message,1)
295  ENDIF
296
297  zdtvr    = daysec/REAL(day_step)
298  IF(dtvr.NE.zdtvr) THEN
299     WRITE(lunout,*) &
300          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
301  ENDIF
302
303  !
304  ! on remet le calendrier \`a zero si demande
305  !
306  IF (start_time /= starttime) then
307     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
308          ,' fichier restart ne correspond pas a celle lue dans le run.def'
309     IF (raz_date == 1) then
310        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
311        start_time = starttime
312     ELSE
313        call abort_gcm("gcm", "'Je m''arrete'", 1)
314     ENDIF
315  ENDIF
316  IF (raz_date == 1) THEN
317     annee_ref = anneeref
318     day_ref = dayref
319     day_ini = dayref
320     itau_dyn = 0
321     itau_phy = 0
322     time_0 = 0.
323     write(lunout,*) &
324          'GCM: On reinitialise a la date lue dans gcm.def'
325  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
326     write(lunout,*) &
327          'GCM: Attention les dates initiales lues dans le fichier'
328     write(lunout,*) &
329          ' restart ne correspondent pas a celles lues dans '
330     write(lunout,*)' gcm.def'
331     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
332     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
333     write(lunout,*)' Pas de remise a zero'
334  ENDIF
335
336  !      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
337  !        write(lunout,*)
338  !     .  'GCM: Attention les dates initiales lues dans le fichier'
339  !        write(lunout,*)
340  !     .  ' restart ne correspondent pas a celles lues dans '
341  !        write(lunout,*)' gcm.def'
342  !        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
343  !        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
344  !        if (raz_date .ne. 1) then
345  !          write(lunout,*)
346  !     .    'GCM: On garde les dates du fichier restart'
347  !        else
348  !          annee_ref = anneeref
349  !          day_ref = dayref
350  !          day_ini = dayref
351  !          itau_dyn = 0
352  !          itau_phy = 0
353  !          time_0 = 0.
354  !          write(lunout,*)
355  !     .   'GCM: On reinitialise a la date lue dans gcm.def'
356  !        endif
357  !      ELSE
358  !        raz_date = 0
359  !      endif
360
361#ifdef CPP_IOIPSL
362  mois = 1
363  heure = 0.
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
377  ! Ehouarn: we still need to define JD_ref and JH_ref
378  ! and since we don't know how many days there are in a year
379  ! we set JD_ref to 0 (this should be improved ...)
380  jD_ref=0
381  jH_ref=0
382#endif
383
384
385  if (iflag_phys.eq.1) then
386     ! these initialisations have already been done (via iniacademic)
387     ! if running in SW or Newtonian mode
388     !-----------------------------------------------------------------------
389     !   Initialisation des constantes dynamiques :
390     !   ------------------------------------------
391     dtvr = zdtvr
392     CALL iniconst
393
394     !-----------------------------------------------------------------------
395     !   Initialisation de la geometrie :
396     !   --------------------------------
397     CALL inigeom
398
399     !-----------------------------------------------------------------------
400     !   Initialisation du filtre :
401     !   --------------------------
402     CALL inifilr
403  endif ! of if (iflag_phys.eq.1)
404  !
405  !-----------------------------------------------------------------------
406  !   Initialisation de la dissipation :
407  !   ----------------------------------
408
409  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
410       tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
411
412  !  numero de stockage pour les fichiers de redemarrage:
413
414  !-----------------------------------------------------------------------
415  !   Initialisation des I/O :
416  !   ------------------------
417
418
419  if (nday>=0) then
420     day_end = day_ini + nday
421  else
422     day_end = day_ini - nday/day_step
423  endif
424  WRITE(lunout,300)day_ini,day_end
425300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
426
427#ifdef CPP_IOIPSL
428  call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
429  write (lunout,301)jour, mois, an
430  call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
431  write (lunout,302)jour, mois, an
432301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
433302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
434#endif
435
436  !-----------------------------------------------------------------------
437  !   Initialisation de la physique :
438  !   -------------------------------
439
440  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
441     ! Physics:
442#ifdef CPP_PHYS
443     CALL iniphysiq(iim,jjm,(jjm-1)*iim+2-1/iim,comm_lmdz, &
444          llm,daysec, dtphys/nsplit_phys, &
445          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
446          iflag_phys)
447#endif
448  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
449
450  !      if (planet_type.eq."earth") then
451  ! Write an Earth-format restart file
452
453  CALL dynredem0("restart.nc", day_end, phis)
454  !      endif
455
456  ecripar = .TRUE.
457
458#ifdef CPP_IOIPSL
459  time_step = zdtvr
460  if (ok_dyn_ins) then
461     ! initialize output file for instantaneous outputs
462     ! t_ops = iecri * daysec ! do operations every t_ops
463     t_ops =((1.0*iecri)/day_step) * daysec 
464     t_wrt = daysec ! iecri * daysec ! write output every t_wrt
465     CALL inithist(day_ref,annee_ref,time_step, &
466          t_ops,t_wrt)
467  endif
468
469  IF (ok_dyn_ave) THEN
470     ! initialize output file for averaged outputs
471     t_ops = iperiod * time_step ! do operations every t_ops
472     t_wrt = periodav * daysec   ! write output every t_wrt
473     CALL initdynav(day_ref,annee_ref,time_step, &
474          t_ops,t_wrt)
475  END IF
476  dtav = iperiod*dtvr/daysec
477#endif
478  ! #endif of #ifdef CPP_IOIPSL
479
480  !  Choix des frequences de stokage pour le offline
481  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
482  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
483  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
484  istphy=istdyn/iphysiq     
485
486
487  !
488  !-----------------------------------------------------------------------
489  !   Integration temporelle du modele :
490  !   ----------------------------------
491
492  !       write(78,*) 'ucov',ucov
493  !       write(78,*) 'vcov',vcov
494  !       write(78,*) 'teta',teta
495  !       write(78,*) 'ps',ps
496  !       write(78,*) 'q',q
497
498
499  CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,time_0)
500
501END PROGRAM gcm
Note: See TracBrowser for help on using the repository browser.