source: LMDZ5/trunk/libf/dyn3dmem/gcm.F90 @ 2339

Last change on this file since 2339 was 2263, checked in by lguez, 10 years ago

In the sequential version, ini_getparam and fin_getparam were not
called so guide_init created a file "fort.99". Added call to
ini_getparam and fin_getparam in guide_init. The created file is now
"nudging_parameters_out.txt".

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