source: LMDZ4/branches/LMDZ4-dev/libf/dyn3d/gcm.F @ 1086

Last change on this file since 1086 was 1086, checked in by yann meurdesoif, 16 years ago

Modifications Othman Bouzi : optimisation du filtre (remplacement opération matrice/vecteurs par matrice/matrice/matrice - BLAS), l'allocation des tableaux du filtre se fait maintenant dynamiquement (plus d'intervention manuelle dans parafiltre.h)

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