source: LMDZ4/branches/V3_test/libf/dyn3d/gcm.F @ 5002

Last change on this file since 5002 was 726, checked in by Laurent Fairhead, 18 years ago

Modifications pour rendre INCA plus independant de LMDZ ACo
LF

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 12.0 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      IMPLICIT NONE
14
15c      ......   Version  du 10/01/98    ..........
16
17c             avec  coordonnees  verticales hybrides
18c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
19
20c=======================================================================
21c
22c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
23c   -------
24c
25c   Objet:
26c   ------
27c
28c   GCM LMD nouvelle grille
29c
30c=======================================================================
31c
32c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
33c      et possibilite d'appeler une fonction f(y)  a derivee tangente
34c      hyperbolique a la  place de la fonction a derivee sinusoidale.
35c  ... Possibilite de choisir le schema pour l'advection de
36c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
37c
38c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
39c      Pour Van-Leer iadv=10
40c
41c-----------------------------------------------------------------------
42c   Declarations:
43c   -------------
44
45#include "dimensions.h"
46#include "paramet.h"
47#include "comconst.h"
48#include "comdissnew.h"
49#include "comvert.h"
50#include "comgeom.h"
51#include "logic.h"
52#include "temps.h"
53#include "control.h"
54#include "ener.h"
55#include "description.h"
56#include "serre.h"
57#include "com_io_dyn.h"
58#include "iniprint.h"
59#include "tracstoke.h"
60#include "advtrac.h"
61
62      INTEGER         longcles
63      PARAMETER     ( longcles = 20 )
64      REAL  clesphy0( longcles )
65      SAVE  clesphy0
66
67
68
69      REAL zdtvr
70      INTEGER nbetatmoy, nbetatdem,nbetat
71
72c   variables dynamiques
73      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
74      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
75      REAL q(ip1jmp1,llm,nqmx)               ! champs advectes
76      REAL ps(ip1jmp1)                       ! pression  au sol
77      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
78      REAL pks(ip1jmp1)                      ! exner au  sol
79      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
80      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
81      REAL masse(ip1jmp1,llm)                ! masse d'air
82      REAL phis(ip1jmp1)                     ! geopotentiel au sol
83      REAL phi(ip1jmp1,llm)                  ! geopotentiel
84      REAL w(ip1jmp1,llm)                    ! vitesse verticale
85
86c variables dynamiques intermediaire pour le transport
87
88c   variables pour le fichier histoire
89      REAL dtav      ! intervalle de temps elementaire
90
91      REAL time_0
92
93      LOGICAL lafin
94      INTEGER ij,iq,l,i,j
95
96
97      real time_step, t_wrt, t_ops
98
99      REAL rdayvrai,rdaym_ini,rday_ecri
100      LOGICAL first
101
102      LOGICAL call_iniphys
103      data call_iniphys/.true./
104
105      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
106c+jld variables test conservation energie
107      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
108C     Tendance de la temp. potentiel d (theta)/ d t due a la
109C     tansformation d'energie cinetique en energie thermique
110C     cree par la dissipation
111      REAL dhecdt(ip1jmp1,llm)
112      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
113      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
114      CHARACTER*15 ztit
115c-jld
116
117
118      character*80 dynhist_file, dynhistave_file
119      character*20 modname
120      character*80 abort_message
121
122C Calendrier
123      LOGICAL true_calendar
124      PARAMETER (true_calendar = .false.)
125
126c-----------------------------------------------------------------------
127c    variables pour l'initialisation de la physique :
128c    ------------------------------------------------
129      INTEGER ngridmx,nq
130      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
131      REAL zcufi(ngridmx),zcvfi(ngridmx)
132      REAL latfi(ngridmx),lonfi(ngridmx)
133      REAL airefi(ngridmx)
134      SAVE latfi, lonfi, airefi
135
136c-----------------------------------------------------------------------
137c   Initialisations:
138c   ----------------
139
140      abort_message = 'last timestep reached'
141      modname = 'gcm'
142      descript = 'Run GCM LMDZ'
143      lafin    = .FALSE.
144      dynhist_file = 'dyn_hist.nc'
145      dynhistave_file = 'dyn_hist_ave.nc'
146
147c initialisation Anne
148      hadv_flg(:) = 0.
149      vadv_flg(:) = 0.
150      conv_flg(:) = 0.
151      pbl_flg(:)  = 0.
152      tracnam(:)  = '        '
153      nprath = 1
154      nbtrac = 0
155      mmt_adj(:,:,:,:) = 1
156
157
158c--------------------------------------------------------------------------
159c   Iflag_phys controle l'appel a la physique :
160c   -------------------------------------------
161c      0 : pas de physique
162c      1 : Normale (appel a phylmd, phymars ...)
163c      2 : rappel Newtonien pour la temperature + friction au sol
164      iflag_phys=1
165
166c--------------------------------------------------------------------------
167c   Lecture de l'etat initial :
168c   ---------------------------
169c     T : on lit start.nc
170c     F : le modele s'autoinitialise avec un cas academique (iniacademic)
171#ifdef CPP_IOIPSL
172      read_start=.true.
173#else
174      read_start=.false.
175#endif
176
177c-----------------------------------------------------------------------
178c   Choix du calendrier
179c   -------------------
180
181#ifdef CPP_IOIPSL
182      if (true_calendar) then
183        call ioconf_calendar('gregorian')
184      else
185        call ioconf_calendar('360d')
186      endif
187#endif
188c----------------------------------------------------------------------
189c  lecture des fichiers gcm.def ou run.def
190c  ---------------------------------------
191c
192#ifdef CPP_IOIPSL
193      CALL conf_gcm( 99, .TRUE. , clesphy0 )
194#else
195      CALL defrun( 99, .TRUE. , clesphy0 )
196#endif
197
198      call InitDimphy
199      call init_phys_openmp
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,klon2,phy_size,klon_para_nb)
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.