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