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

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