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

Last change on this file since 2353 was 2351, checked in by Ehouarn Millour, 9 years ago

More on physics/dynamics separation and cleanup:

  • Set things up so that all physics-related initializations are done via iniphysiq.
  • Created a "geometry_mod.F90" module in phy_common to store information on the loacl grid (i.e. replaces comgeomphy) and moreover give these variables more obvious names (e.g.: rlond => longitude, rlatd => latitude, airephy => cell_area).
  • removed obsolete comgeomphy.h and comgeomphy.F90

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