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

Last change on this file since 1032 was 1017, checked in by lsce, 16 years ago
  • Copure des lignes trop longues

JG

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