source: LMDZ4/trunk/libf/dyn3d/gcm.F @ 937

Last change on this file since 937 was 863, checked in by lsce, 17 years ago

ACo : correction bug pour l'appel de INCA en mode sequentiel

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