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

Last change on this file since 1085 was 1085, checked in by jghattas, 15 years ago

Bug fix pour OpenMP.
Anne Cozic

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