source: LMDZ4/branches/LMDZ4-dev/libf/dyn3dpar/gcm.F @ 1195

Last change on this file since 1195 was 1195, checked in by yann meurdesoif, 15 years ago

Correction : getparam n'etait pas initialise, chaque processeur ecrivait dans le meme fichier 99.fort, le code plantait sur mercure en fin de run lors du close(99) dans leapfrog_p.F

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 13.5 KB
Line 
1!
2! $Header$
3!
4c
5c
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 mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
15      USE infotrac
16      USE mod_interface_dyn_phys
17      USE mod_hallo
18      USE Bands
19      USE getparam
20      USE filtreg_mod
21
22! Ehouarn: for now these only apply to Earth:
23#ifdef CPP_EARTH
24      USE mod_grid_phy_lmdz
25      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
26      USE dimphy
27      USE comgeomphy
28#endif
29      IMPLICIT NONE
30
31c      ......   Version  du 10/01/98    ..........
32
33c             avec  coordonnees  verticales hybrides
34c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
35
36c=======================================================================
37c
38c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
39c   -------
40c
41c   Objet:
42c   ------
43c
44c   GCM LMD nouvelle grille
45c
46c=======================================================================
47c
48c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
49c      et possibilite d'appeler une fonction f(y)  a derivee tangente
50c      hyperbolique a la  place de la fonction a derivee sinusoidale.
51c  ... Possibilite de choisir le schema pour l'advection de
52c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
53c
54c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
55c      Pour Van-Leer iadv=10
56c
57c-----------------------------------------------------------------------
58c   Declarations:
59c   -------------
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 "control.h"
69#include "ener.h"
70#include "description.h"
71#include "serre.h"
72#include "com_io_dyn.h"
73#include "iniprint.h"
74#include "tracstoke.h"
75
76      INTEGER         longcles
77      PARAMETER     ( longcles = 20 )
78      REAL  clesphy0( longcles )
79      SAVE  clesphy0
80
81
82
83      REAL zdtvr
84c      INTEGER nbetatmoy, nbetatdem,nbetat
85      INTEGER nbetatmoy, nbetatdem
86
87c   variables dynamiques
88      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
89      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
90      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
91      REAL ps(ip1jmp1)                       ! pression  au sol
92c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
93c      REAL pks(ip1jmp1)                      ! exner au  sol
94c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
95c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
96      REAL masse(ip1jmp1,llm)                ! masse d'air
97      REAL phis(ip1jmp1)                     ! geopotentiel au sol
98c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
99c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
100
101c variables dynamiques intermediaire pour le transport
102
103c   variables pour le fichier histoire
104      REAL dtav      ! intervalle de temps elementaire
105
106      REAL time_0
107
108      LOGICAL lafin
109c      INTEGER ij,iq,l,i,j
110      INTEGER i,j
111
112
113      real time_step, t_wrt, t_ops
114
115c      REAL rdayvrai,rdaym_ini,rday_ecri
116c      LOGICAL first
117
118      LOGICAL call_iniphys
119      data call_iniphys/.true./
120
121c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
122c+jld variables test conservation energie
123c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
124C     Tendance de la temp. potentiel d (theta)/ d t due a la
125C     tansformation d'energie cinetique en energie thermique
126C     cree par la dissipation
127c      REAL dhecdt(ip1jmp1,llm)
128c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
129c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
130c      CHARACTER (len=15) :: ztit
131c-jld
132
133
134      character (len=80) :: dynhist_file, dynhistave_file
135      character (len=20) ::modname
136      character (len=80) ::abort_message
137
138C Calendrier
139      LOGICAL true_calendar
140      PARAMETER (true_calendar = .false.)
141
142c-----------------------------------------------------------------------
143c    variables pour l'initialisation de la physique :
144c    ------------------------------------------------
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
155c-----------------------------------------------------------------------
156c   Initialisations:
157c   ----------------
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
167c-----------------------------------------------------------------------
168c   Choix du calendrier
169c   -------------------
170
171#ifdef CPP_IOIPSL
172      if (true_calendar) then
173        call ioconf_calendar('gregorian')
174      else
175        call ioconf_calendar('360d')
176      endif
177#endif
178c----------------------------------------------------------------------
179c  lecture des fichiers gcm.def ou run.def
180c  ---------------------------------------
181c
182! Ehouarn: dump possibility of using defrun
183!#ifdef CPP_IOIPSL
184      CALL conf_gcm( 99, .TRUE. , clesphy0 )
185!#else
186!      CALL defrun( 99, .TRUE. , clesphy0 )
187!#endif
188c
189c
190c------------------------------------
191c   Initialisation partie parallele
192c------------------------------------
193      CALL init_const_mpi
194
195      call init_parallel
196      call ini_getparam("out.def")
197      call Read_Distrib
198! Ehouarn : temporarily (?) keep this only for Earth
199      if (planet_type.eq."earth") then
200#ifdef CPP_EARTH
201        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
202#endif
203      endif ! of if (planet_type.eq."earth")
204      CALL set_bands
205      CALL Init_interface_dyn_phys
206      CALL barrier
207
208      if (mpi_rank==0) call WriteBands
209      call SetDistrib(jj_Nb_Caldyn)
210
211c$OMP PARALLEL
212      call Init_Mod_hallo
213c$OMP END PARALLEL
214
215! Ehouarn : temporarily (?) keep this only for Earth
216      if (planet_type.eq."earth") then
217#ifdef CPP_EARTH
218c$OMP PARALLEL
219      call InitComgeomphy
220c$OMP END PARALLEL
221#endif
222      endif ! of if (planet_type.eq."earth")
223
224      IF (config_inca /= 'none') THEN
225#ifdef INCA
226         call init_const_lmdz(
227     $        nbtr,anneeref,dayref,
228     $        iphysiq,day_step,nday)
229
230         call init_inca_para(
231     $        iim,jjm+1,llm,klon_glo,mpi_size,
232     $        distrib_phys,COMM_LMDZ)
233#endif
234      END IF
235
236c-----------------------------------------------------------------------
237c   Initialisation des traceurs
238c   ---------------------------
239c  Choix du nombre de traceurs et du schema pour l'advection
240c  dans fichier traceur.def, par default ou via INCA
241      call infotrac_init
242
243c Allocation de la tableau q : champs advectes   
244      ALLOCATE(q(ip1jmp1,llm,nqtot))
245
246c-----------------------------------------------------------------------
247c   Lecture de l'etat initial :
248c   ---------------------------
249
250c  lecture du fichier start.nc
251      if (read_start) then
252      ! we still need to run iniacademic to initialize some
253      ! constants & fields, if we run the 'newtonian' case:
254        if (iflag_phys.eq.2) then
255          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
256        endif
257!#ifdef CPP_IOIPSL
258        if (planet_type.eq."earth") then
259#ifdef CPP_EARTH
260! Load an Earth-format start file
261         CALL dynetat0("start.nc",vcov,ucov,
262     .              teta,q,masse,ps,phis, time_0)
263#endif
264        endif ! of if (planet_type.eq."earth")
265c       write(73,*) 'ucov',ucov
266c       write(74,*) 'vcov',vcov
267c       write(75,*) 'teta',teta
268c       write(76,*) 'ps',ps
269c       write(77,*) 'q',q
270
271      endif ! of if (read_start)
272
273c le cas echeant, creation d un etat initial
274      IF (prt_level > 9) WRITE(lunout,*)
275     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
276      if (.not.read_start) then
277         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
278      endif
279
280c-----------------------------------------------------------------------
281c   Lecture des parametres de controle pour la simulation :
282c   -------------------------------------------------------
283c  on recalcule eventuellement le pas de temps
284
285      IF(MOD(day_step,iperiod).NE.0) THEN
286        abort_message =
287     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
288        call abort_gcm(modname,abort_message,1)
289      ENDIF
290
291      IF(MOD(day_step,iphysiq).NE.0) THEN
292        abort_message =
293     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
294        call abort_gcm(modname,abort_message,1)
295      ENDIF
296
297      zdtvr    = daysec/FLOAT(day_step)
298        IF(dtvr.NE.zdtvr) THEN
299         WRITE(lunout,*)
300     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
301        ENDIF
302
303C
304C on remet le calendrier à zero si demande
305c
306      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
307        write(lunout,*)
308     .  ' Attention les dates initiales lues dans le fichier'
309        write(lunout,*)
310     .  ' restart ne correspondent pas a celles lues dans '
311        write(lunout,*)' gcm.def'
312        if (raz_date .ne. 1) then
313          write(lunout,*)
314     .    ' On garde les dates du fichier restart'
315        else
316          annee_ref = anneeref
317          day_ref = dayref
318          day_ini = dayref
319          itau_dyn = 0
320          itau_phy = 0
321          time_0 = 0.
322          write(lunout,*)
323     .   ' On reinitialise a la date lue dans gcm.def'
324        endif
325      ELSE
326        raz_date = 0
327      endif
328
329
330c  nombre d'etats dans les fichiers demarrage et histoire
331      nbetatdem = nday / iecri
332      nbetatmoy = nday / periodav + 1
333
334c-----------------------------------------------------------------------
335c   Initialisation des constantes dynamiques :
336c   ------------------------------------------
337      dtvr = zdtvr
338      CALL iniconst
339
340c-----------------------------------------------------------------------
341c   Initialisation de la geometrie :
342c   --------------------------------
343      CALL inigeom
344
345c-----------------------------------------------------------------------
346c   Initialisation du filtre :
347c   --------------------------
348      CALL inifilr
349c
350c-----------------------------------------------------------------------
351c   Initialisation de la dissipation :
352c   ----------------------------------
353
354      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
355     *                tetagdiv, tetagrot , tetatemp              )
356
357c-----------------------------------------------------------------------
358c   Initialisation de la physique :
359c   -------------------------------
360      IF (call_iniphys.and.iflag_phys.eq.1) THEN
361         latfi(1)=rlatu(1)
362         lonfi(1)=0.
363         zcufi(1) = cu(1)
364         zcvfi(1) = cv(1)
365         DO j=2,jjm
366            DO i=1,iim
367               latfi((j-2)*iim+1+i)= rlatu(j)
368               lonfi((j-2)*iim+1+i)= rlonv(i)
369               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
370               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
371            ENDDO
372         ENDDO
373         latfi(ngridmx)= rlatu(jjp1)
374         lonfi(ngridmx)= 0.
375         zcufi(ngridmx) = cu(ip1jm+1)
376         zcvfi(ngridmx) = cv(ip1jm-iim)
377         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
378
379         WRITE(lunout,*)
380     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
381! Earth:
382         if (planet_type.eq."earth") then
383#ifdef CPP_EARTH
384         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
385     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
386#endif
387         endif ! of if (planet_type.eq."earth")
388         call_iniphys=.false.
389      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
390
391
392c-----------------------------------------------------------------------
393c   Initialisation des dimensions d'INCA :
394c   --------------------------------------
395      IF (config_inca /= 'none') THEN
396!$OMP PARALLEL
397#ifdef INCA
398         CALL init_inca_dim(klon_omp,llm,iim,jjm,
399     $        rlonu,rlatu,rlonv,rlatv)
400#endif
401!$OMP END PARALLEL
402      END IF
403
404c-----------------------------------------------------------------------
405c   Initialisation des I/O :
406c   ------------------------
407
408
409      day_end = day_ini + nday
410      WRITE(lunout,300)day_ini,day_end
411 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
412
413!#ifdef CPP_IOIPSL
414      if (planet_type.eq."earth") then
415#ifdef CPP_EARTH
416      CALL dynredem0_p("restart.nc", day_end, phis)
417#endif
418      endif
419
420      ecripar = .TRUE.
421
422#ifdef CPP_IOIPSL
423      if ( 1.eq.1) then
424      time_step = zdtvr
425      t_ops = iecri * daysec
426      t_wrt = iecri * daysec
427      CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
428     .              t_ops, t_wrt, histid, histvid)
429
430      IF (ok_dynzon) THEN
431         t_ops = iperiod * time_step
432         t_wrt = periodav * daysec
433         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
434     .        t_ops, t_wrt, histaveid)
435      END IF
436      dtav = iperiod*dtvr/daysec
437      endif
438
439
440#endif
441! #endif of #ifdef CPP_IOIPSL
442
443c  Choix des frequences de stokage pour le offline
444c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
445c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
446      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
447      istphy=istdyn/iphysiq     
448
449
450c
451c-----------------------------------------------------------------------
452c   Integration temporelle du modele :
453c   ----------------------------------
454
455c       write(78,*) 'ucov',ucov
456c       write(78,*) 'vcov',vcov
457c       write(78,*) 'teta',teta
458c       write(78,*) 'ps',ps
459c       write(78,*) 'q',q
460
461c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
462      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
463     .              time_0)
464c$OMP END PARALLEL
465
466
467      END
468
Note: See TracBrowser for help on using the repository browser.