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 #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 vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
91  REAL teta(ip1jmp1,llm) ! temperature potentielle
92  REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
93  REAL ps(ip1jmp1) ! 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 masse(ip1jmp1,llm) ! masse d'air
99  REAL phis(ip1jmp1) ! 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 setdistrib(jj_nb_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(q(ip1jmp1,llm,nqtot))
256 
257 c-----------------------------------------------------------------------
258 c Lecture de l'etat initial :
259 c ---------------------------
260 
261 c lecture du fichier start.nc
262  if (read_start) then
263  ! we still need to run iniacademic to initialize some
264  ! constants & fields, if we run the 'newtonian' or 'SW' cases:
265  if (iflag_phys.ne.1) then
266  CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
267  endif
268 
269 ! if (planet_type.eq."earth") then
270 ! Load an Earth-format start file
271  CALL dynetat0("start.nc",vcov,ucov,
272  & teta,q,masse,ps,phis, time_0)
273 ! endif ! of if (planet_type.eq."earth")
274 
275 c write(73,*) 'ucov',ucov
276 c write(74,*) 'vcov',vcov
277 c write(75,*) 'teta',teta
278 c write(76,*) 'ps',ps
279 c write(77,*) 'q',q
280 
281  endif ! of if (read_start)
282 
283 c le cas echeant, creation d un etat initial
284  IF (prt_level > 9) WRITE(lunout,*)
285  . 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
286  if (.not.read_start) then
287  CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
288  endif
289 
290 c-----------------------------------------------------------------------
291 c Lecture des parametres de controle pour la simulation :
292 c -------------------------------------------------------
293 c on recalcule eventuellement le pas de temps
294 
295  IF(mod(day_step,iperiod).NE.0) THEN
296  abort_message =
297  . 'Il faut choisir un nb de pas par jour multiple de iperiod'
298  call abort_gcm(modname,abort_message,1)
299  ENDIF
300 
301  IF(mod(day_step,iphysiq).NE.0) THEN
302  abort_message =
303  * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
304  call abort_gcm(modname,abort_message,1)
305  ENDIF
306 
307  zdtvr = daysec/REAL(day_step)
308  IF(dtvr.NE.zdtvr) THEN
309  WRITE(lunout,*)
310  . 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
311  ENDIF
312 
313 C
314 C on remet le calendrier à zero si demande
315 c
316  IF (start_time /= starttime) then
317  WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
318  &,' fichier restart ne correspond pas à celle lue dans le run.def'
319  IF (raz_date == 1) then
320  WRITE(lunout,*)'Je prends l''heure lue dans run.def'
321  start_time = starttime
322  ELSE
323  WRITE(lunout,*)'Je m''arrete'
324  CALL abort
325  ENDIF
326  ENDIF
327  IF (raz_date == 1) THEN
328  annee_ref = anneeref
329  day_ref = dayref
330  day_ini = dayref
331  itau_dyn = 0
332  itau_phy = 0
333  time_0 = 0.
334  write(lunout,*)
335  . 'GCM: On reinitialise a la date lue dans gcm.def'
336  ELSE 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  write(lunout,*)' Pas de remise a zero'
345  ENDIF
346 c if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
347 c write(lunout,*)
348 c . 'GCM: Attention les dates initiales lues dans le fichier'
349 c write(lunout,*)
350 c . ' restart ne correspondent pas a celles lues dans '
351 c write(lunout,*)' gcm.def'
352 c write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
353 c write(lunout,*)' day_ref=',day_ref," dayref=",dayref
354 c if (raz_date .ne. 1) then
355 c write(lunout,*)
356 c . 'GCM: On garde les dates du fichier restart'
357 c else
358 c annee_ref = anneeref
359 c day_ref = dayref
360 c day_ini = dayref
361 c itau_dyn = 0
362 c itau_phy = 0
363 c time_0 = 0.
364 c write(lunout,*)
365 c . 'GCM: On reinitialise a la date lue dans gcm.def'
366 c endif
367 c ELSE
368 c raz_date = 0
369 c endif
370 
371 #ifdef CPP_IOIPSL
372  mois = 1
373  heure = 0.
374  call ymds2ju(annee_ref, mois, day_ref, heure, jd_ref)
375  jh_ref = jd_ref - int(jd_ref)
376  jd_ref = int(jd_ref)
377 
378  call ioconf_startdate(int(jd_ref), jh_ref)
379 
380  write(lunout,*)'DEBUG'
381  write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
382  write(lunout,*)annee_ref, mois, day_ref, heure, jd_ref
383  call ju2ymds(jd_ref+jh_ref,an, mois, jour, heure)
384  write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
385  write(lunout,*)jd_ref+jh_ref,an, mois, jour, heure
386 #else
387 ! Ehouarn: we still need to define JD_ref and JH_ref
388 ! and since we don't know how many days there are in a year
389 ! we set JD_ref to 0 (this should be improved ...)
390  jd_ref=0
391  jh_ref=0
392 #endif
393 
394 
395  if (iflag_phys.eq.1) then
396  ! these initialisations have already been done (via iniacademic)
397  ! if running in SW or Newtonian mode
398 c-----------------------------------------------------------------------
399 c Initialisation des constantes dynamiques :
400 c ------------------------------------------
401  dtvr = zdtvr
402  CALL iniconst
403 
404 c-----------------------------------------------------------------------
405 c Initialisation de la geometrie :
406 c --------------------------------
407  CALL inigeom
408 
409 c-----------------------------------------------------------------------
410 c Initialisation du filtre :
411 c --------------------------
412  CALL inifilr
413  endif ! of if (iflag_phys.eq.1)
414 c
415 c-----------------------------------------------------------------------
416 c Initialisation de la dissipation :
417 c ----------------------------------
418 
419  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
420  * tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
421 
422 c-----------------------------------------------------------------------
423 c Initialisation de la physique :
424 c -------------------------------
425  IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
426  latfi(1)=rlatu(1)
427  lonfi(1)=0.
428  zcufi(1) = cu(1)
429  zcvfi(1) = cv(1)
430  DO j=2,jjm
431  DO i=1,iim
432  latfi((j-2)*iim+1+i)= rlatu(j)
433  lonfi((j-2)*iim+1+i)= rlonv(i)
434  zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
435  zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
436  ENDDO
437  ENDDO
438  latfi(ngridmx)= rlatu(jjp1)
439  lonfi(ngridmx)= 0.
440  zcufi(ngridmx) = cu(ip1jm+1)
441  zcvfi(ngridmx) = cv(ip1jm-iim)
442  CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
443 
444  WRITE(lunout,*)
445  . 'GCM: WARNING!!! vitesse verticale nulle dans la physique'
446 ! Physics:
447 #ifdef CPP_PHYS
448  CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
449  & latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
450  & iflag_phys)
451 #endif
452  call_iniphys=.false.
453  ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
454 
455 
456 c-----------------------------------------------------------------------
457 c Initialisation des dimensions d'INCA :
458 c --------------------------------------
459  IF (type_trac == 'inca') THEN
460 !$OMP PARALLEL
461 #ifdef INCA
462  CALL init_inca_dim(klon_omp,llm,iim,jjm,
463  $ rlonu,rlatu,rlonv,rlatv)
464 #endif
465 !$OMP END PARALLEL
466  END IF
467 
468 c-----------------------------------------------------------------------
469 c Initialisation des I/O :
470 c ------------------------
471 
472 
473  day_end = day_ini + nday
474  WRITE(lunout,300)day_ini,day_end
475  300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
476 
477 #ifdef CPP_IOIPSL
478  call ju2ymds(jd_ref + day_ini - day_ref, an, mois, jour, heure)
479  write (lunout,301)jour, mois, an
480  call ju2ymds(jd_ref + day_end - day_ref, an, mois, jour, heure)
481  write (lunout,302)jour, mois, an
482  301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
483  302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4)
484 #endif
485 
486 ! if (planet_type.eq."earth") then
487 ! Write an Earth-format restart file
488  CALL dynredem0_p("restart.nc", day_end, phis)
489 ! endif
490 
491  ecripar = .true.
492 
493 #ifdef CPP_IOIPSL
494  time_step = zdtvr
495  IF (mpi_rank==0) then
496  if (ok_dyn_ins) then
497  ! initialize output file for instantaneous outputs
498  ! t_ops = iecri * daysec ! do operations every t_ops
499  t_ops =((1.0*iecri)/day_step) * daysec
500  t_wrt = daysec ! iecri * daysec ! write output every t_wrt
501  t_wrt = daysec ! iecri * daysec ! write output every t_wrt
502  CALL inithist(day_ref,annee_ref,time_step,
503  & t_ops,t_wrt)
504  endif
505 
506  IF (ok_dyn_ave) THEN
507  ! initialize output file for averaged outputs
508  t_ops = iperiod * time_step ! do operations every t_ops
509  t_wrt = periodav * daysec ! write output every t_wrt
510  CALL initdynav(day_ref,annee_ref,time_step,
511  & t_ops,t_wrt)
512 ! CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
513 ! . t_ops, t_wrt, histaveid)
514  END IF
515  ENDIF
516  dtav = iperiod*dtvr/daysec
517 #endif
518 ! #endif of #ifdef CPP_IOIPSL
519 
520 c Choix des frequences de stokage pour le offline
521 c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
522 c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12
523  istdyn=day_step/4 ! stockage toutes les 6h=1jour/12
524  istphy=istdyn/iphysiq
525 
526 
527 c
528 c-----------------------------------------------------------------------
529 c Integration temporelle du modele :
530 c ----------------------------------
531 
532 c write(78,*) 'ucov',ucov
533 c write(78,*) 'vcov',vcov
534 c write(78,*) 'teta',teta
535 c write(78,*) 'ps',ps
536 c write(78,*) 'q',q
537 
538 c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
539  CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
540  . time_0)
541 c$OMP END PARALLEL
542 
543 
544  END
545