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

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

Make iniphysiq a module.
Fix call to iniphysiq in lmdz1d (missing arguments and arrays of wrong sizes).
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.0 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,daysec,day_ini,dtphys/nsplit_phys, &
400          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
401          iflag_phys)
402#endif
403  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
404
405
406  !-----------------------------------------------------------------------
407  !   Initialisation des dimensions d'INCA :
408  !   --------------------------------------
409  IF (type_trac == 'inca') THEN
410     !$OMP PARALLEL
411#ifdef INCA
412     CALL init_inca_dim(klon_omp,llm,iim,jjm, &
413          rlonu,rlatu,rlonv,rlatv)
414#endif
415     !$OMP END PARALLEL
416  END IF
417
418  !-----------------------------------------------------------------------
419  !   Initialisation des I/O :
420  !   ------------------------
421
422
423  if (nday>=0) then
424     day_end = day_ini + nday
425  else
426     day_end = day_ini - nday/day_step
427  endif
428
429  WRITE(lunout,300)day_ini,day_end
430300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
431
432#ifdef CPP_IOIPSL
433  call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
434  write (lunout,301)jour, mois, an
435  call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
436  write (lunout,302)jour, mois, an
437301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
438302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
439#endif
440
441  !      if (planet_type.eq."earth") then
442  ! Write an Earth-format restart file
443  CALL dynredem0_loc("restart.nc", day_end, phis)
444  !      endif
445
446  ecripar = .TRUE.
447
448#ifdef CPP_IOIPSL
449  time_step = zdtvr
450  IF (mpi_rank==0) then
451     if (ok_dyn_ins) then
452        ! initialize output file for instantaneous outputs
453        ! t_ops = iecri * daysec ! do operations every t_ops
454        t_ops =((1.0*iecri)/day_step) * daysec 
455        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
456        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
457        CALL inithist(day_ref,annee_ref,time_step, &
458             t_ops,t_wrt)
459     endif
460
461     IF (ok_dyn_ave) THEN
462        ! initialize output file for averaged outputs
463        t_ops = iperiod * time_step ! do operations every t_ops
464        t_wrt = periodav * daysec   ! write output every t_wrt
465        CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
466     END IF
467  ENDIF
468  dtav = iperiod*dtvr/daysec
469#endif
470  ! #endif of #ifdef CPP_IOIPSL
471
472  !  Choix des frequences de stokage pour le offline
473  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
474  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
475  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
476  istphy=istdyn/iphysiq     
477
478
479  !
480  !-----------------------------------------------------------------------
481  !   Integration temporelle du modele :
482  !   ----------------------------------
483
484  !       write(78,*) 'ucov',ucov
485  !       write(78,*) 'vcov',vcov
486  !       write(78,*) 'teta',teta
487  !       write(78,*) 'ps',ps
488  !       write(78,*) 'q',q
489
490  !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
491  CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
492  !$OMP END PARALLEL
493
494  !      OPEN(unit=5487,file='ok_lmdz',status='replace')
495  !      WRITE(5487,*) 'ok_lmdz'
496  !      CLOSE(5487)
497END PROGRAM gcm
Note: See TracBrowser for help on using the repository browser.