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

Last change on this file since 1032 was 962, checked in by Laurent Fairhead, 16 years ago

Un peu de menage pour avoir moins de sorties a la compilation sur Nec
LF

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