source: LMDZ4/branches/LMDZ4_V3_patches/libf/dyn3d/gcm.F @ 5026

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