source: LMDZ4/branches/unlabeled-1.1.1/libf/dyn3d/gcm.F @ 3722

Last change on this file since 3722 was 524, checked in by lmdzadmin, 20 years ago

Initial revision

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 11.2 KB
Line 
1!
2! $Header$
3!
4c
5c
6      PROGRAM gcm
7
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#endif
11
12      IMPLICIT NONE
13
14c      ......   Version  du 10/01/98    ..........
15
16c             avec  coordonnees  verticales hybrides
17c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
18
19c=======================================================================
20c
21c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
22c   -------
23c
24c   Objet:
25c   ------
26c
27c   GCM LMD nouvelle grille
28c
29c=======================================================================
30c
31c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
32c      et possibilite d'appeler une fonction f(y)  a derivee tangente
33c      hyperbolique a la  place de la fonction a derivee sinusoidale.
34c  ... Possibilite de choisir le schema pour l'advection de
35c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
36c
37c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
38c      Pour Van-Leer iadv=10
39c
40c-----------------------------------------------------------------------
41c   Declarations:
42c   -------------
43
44#include "dimensions.h"
45#include "paramet.h"
46#include "comconst.h"
47#include "comdissnew.h"
48#include "comvert.h"
49#include "comgeom.h"
50#include "logic.h"
51#include "temps.h"
52#include "control.h"
53#include "ener.h"
54#include "description.h"
55#include "serre.h"
56#include "com_io_dyn.h"
57#include "iniprint.h"
58
59c#include "tracstoke.h"
60
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
115      INTEGER   ip_ebil_dyn  ! PRINT level for energy conserv. diag.
116      SAVE      ip_ebil_dyn
117      DATA      ip_ebil_dyn/0/
118c-jld
119
120
121      LOGICAL offline  ! Controle du stockage ds "fluxmass"
122      PARAMETER (offline=.false.)
123
124      character*80 dynhist_file, dynhistave_file
125      character*20 modname
126      character*80 abort_message
127
128C Calendrier
129      LOGICAL true_calendar
130      PARAMETER (true_calendar = .false.)
131
132c-----------------------------------------------------------------------
133c    variables pour l'initialisation de la physique :
134c    ------------------------------------------------
135      INTEGER ngridmx,nq
136      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
137      REAL zcufi(ngridmx),zcvfi(ngridmx)
138      REAL latfi(ngridmx),lonfi(ngridmx)
139      REAL airefi(ngridmx)
140      SAVE latfi, lonfi, airefi
141
142c-----------------------------------------------------------------------
143c   Initialisations:
144c   ----------------
145
146      abort_message = 'last timestep reached'
147      modname = 'gcm'
148      descript = 'Run GCM LMDZ'
149      lafin    = .FALSE.
150      dynhist_file = 'dyn_hist.nc'
151      dynhistave_file = 'dyn_hist_ave.nc'
152
153c--------------------------------------------------------------------------
154c   Iflag_phys controle l'appel a la physique :
155c   -------------------------------------------
156c      0 : pas de physique
157c      1 : Normale (appel a phylmd, phymars ...)
158c      2 : rappel Newtonien pour la temperature + friction au sol
159      iflag_phys=1
160
161c--------------------------------------------------------------------------
162c   Lecture de l'etat initial :
163c   ---------------------------
164c     T : on lit start.nc
165c     F : le modele s'autoinitialise avec un cas academique (iniacademic)
166#ifdef CPP_IOIPSL
167      read_start=.true.
168#else
169      read_start=.false.
170#endif
171
172c-----------------------------------------------------------------------
173c   Choix du calendrier
174c   -------------------
175
176#ifdef CPP_IOIPSL
177      if (true_calendar) then
178        call ioconf_calendar('gregorian')
179      else
180        call ioconf_calendar('360d')
181      endif
182#endif
183c----------------------------------------------------------------------
184c  lecture des fichiers gcm.def ou run.def
185c  ---------------------------------------
186c
187#ifdef CPP_IOIPSL
188      CALL conf_gcm( 99, .TRUE. , clesphy0 )
189#else
190      CALL defrun( 99, .TRUE. , clesphy0 )
191#endif
192c
193c
194c-----------------------------------------------------------------------
195c   Initialisation des traceurs
196c   ---------------------------
197c  Choix du schema pour l'advection
198c  dans fichier trac.def ou via INCA
199
200       call iniadvtrac(nq)
201c
202c-----------------------------------------------------------------------
203c   Lecture de l'etat initial :
204c   ---------------------------
205
206c  lecture du fichier start.nc
207      if (read_start) then
208#ifdef CPP_IOIPSL
209         CALL dynetat0("start.nc",nqmx,vcov,ucov,
210     .              teta,q,masse,ps,phis, time_0)
211c       write(73,*) 'ucov',ucov
212c       write(74,*) 'vcov',vcov
213c       write(75,*) 'teta',teta
214c       write(76,*) 'ps',ps
215c       write(77,*) 'q',q
216
217#endif
218      endif
219
220
221
222c le cas echeant, creation d un etat initial
223      IF (prt_level > 9) WRITE(lunout,*)
224     .                 'AVANT iniacademic AVANT AVANT AVANT AVANT'
225      if (.not.read_start) then
226         CALL iniacademic(nqmx,vcov,ucov,teta,q,masse,ps,phis,time_0)
227      endif
228
229
230c-----------------------------------------------------------------------
231c   Lecture des parametres de controle pour la simulation :
232c   -------------------------------------------------------
233c  on recalcule eventuellement le pas de temps
234
235      IF(MOD(day_step,iperiod).NE.0) THEN
236        abort_message =
237     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
238        call abort_gcm(modname,abort_message,1)
239      ENDIF
240
241      IF(MOD(day_step,iphysiq).NE.0) THEN
242        abort_message =
243     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
244        call abort_gcm(modname,abort_message,1)
245      ENDIF
246
247      zdtvr    = daysec/FLOAT(day_step)
248        IF(dtvr.NE.zdtvr) THEN
249         WRITE(lunout,*)
250     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
251        ENDIF
252
253C
254C on remet le calendrier à zero si demande
255c
256      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
257        write(lunout,*)
258     .  ' Attention les dates initiales lues dans le fichier'
259        write(lunout,*)
260     .  ' restart ne correspondent pas a celles lues dans '
261        write(lunout,*)' gcm.def'
262        if (raz_date .ne. 1) then
263          write(lunout,*)
264     .    ' On garde les dates du fichier restart'
265        else
266          annee_ref = anneeref
267          day_ref = dayref
268          day_ini = dayref
269          itau_dyn = 0
270          itau_phy = 0
271          time_0 = 0.
272          write(lunout,*)
273     .   ' On reinitialise a la date lue dans gcm.def'
274        endif
275      ELSE
276        raz_date = 0
277      endif
278
279
280c  nombre d'etats dans les fichiers demarrage et histoire
281      nbetatdem = nday / iecri
282      nbetatmoy = nday / periodav + 1
283
284c-----------------------------------------------------------------------
285c   Initialisation des constantes dynamiques :
286c   ------------------------------------------
287      dtvr = zdtvr
288      CALL iniconst
289
290c-----------------------------------------------------------------------
291c   Initialisation de la geometrie :
292c   --------------------------------
293      CALL inigeom
294
295c-----------------------------------------------------------------------
296c   Initialisation du filtre :
297c   --------------------------
298      CALL inifilr
299c
300c-----------------------------------------------------------------------
301c   Initialisation de la dissipation :
302c   ----------------------------------
303
304      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
305     *                tetagdiv, tetagrot , tetatemp              )
306
307c-----------------------------------------------------------------------
308c   Initialisation de la physique :
309c   -------------------------------
310#ifdef CPP_PHYS
311      IF (call_iniphys.and.iflag_phys.eq.1) THEN
312         latfi(1)=rlatu(1)
313         lonfi(1)=0.
314         zcufi(1) = cu(1)
315         zcvfi(1) = cv(1)
316         DO j=2,jjm
317            DO i=1,iim
318               latfi((j-2)*iim+1+i)= rlatu(j)
319               lonfi((j-2)*iim+1+i)= rlonv(i)
320               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
321               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
322            ENDDO
323         ENDDO
324         latfi(ngridmx)= rlatu(jjp1)
325         lonfi(ngridmx)= 0.
326         zcufi(ngridmx) = cu(ip1jm+1)
327         zcvfi(ngridmx) = cv(ip1jm-iim)
328         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
329         WRITE(lunout,*)
330     .           'WARNING!!! vitesse verticale nulle dans la physique'
331         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
332     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
333         call_iniphys=.false.
334      ENDIF
335#endif
336
337c  numero de stockage pour les fichiers de redemarrage:
338
339c-----------------------------------------------------------------------
340c   Initialisation des I/O :
341c   ------------------------
342
343
344      day_end = day_ini + nday
345      WRITE(lunout,300)day_ini,day_end
346
347#ifdef CPP_IOIPSL
348      CALL dynredem0("restart.nc", day_end, phis, nqmx)
349
350      ecripar = .TRUE.
351
352      if ( 1.eq.1) then
353      time_step = zdtvr
354      t_ops = iecri * daysec
355      t_wrt = iecri * daysec
356      CALL inithist(dynhist_file,day_ref,annee_ref,time_step,
357     .              t_ops, t_wrt, nqmx, histid, histvid)
358
359      t_ops = iperiod * time_step
360      t_wrt = periodav * daysec
361      CALL initdynav(dynhistave_file,day_ref,annee_ref,time_step,
362     .              t_ops, t_wrt, nqmx, histaveid)
363
364      dtav = iperiod*dtvr/daysec
365      endif
366
367
368#endif
369
370c
371c-----------------------------------------------------------------------
372c   Integration temporelle du modele :
373c   ----------------------------------
374
375c       write(78,*) 'ucov',ucov
376c       write(78,*) 'vcov',vcov
377c       write(78,*) 'teta',teta
378c       write(78,*) 'ps',ps
379c       write(78,*) 'q',q
380
381
382      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0)
383
384
385
386 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
387     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
388      END
389
Note: See TracBrowser for help on using the repository browser.