source: LMDZ4/trunk/libf/dyn3dpar/gcm.F @ 801

Last change on this file since 801 was 791, checked in by lsce, 17 years ago

ACo : mise a jour pour inca

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