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