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

Last change on this file since 1086 was 1086, checked in by yann meurdesoif, 15 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: 13.1 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
21      USE filtre
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#include "dimensions.h"
55#include "paramet.h"
56#include "comconst.h"
57#include "comdissnew.h"
58#include "comvert.h"
59#include "comgeom.h"
60#include "logic.h"
61#include "temps.h"
62#include "control.h"
63#include "ener.h"
64#include "description.h"
65#include "serre.h"
66#include "com_io_dyn.h"
67#include "iniprint.h"
68#include "tracstoke.h"
69#include "advtrac.h"
70
71      INTEGER         longcles
72      PARAMETER     ( longcles = 20 )
73      REAL  clesphy0( longcles )
74      SAVE  clesphy0
75
76
77
78      REAL zdtvr
79c      INTEGER nbetatmoy, nbetatdem,nbetat
80      INTEGER nbetatmoy, nbetatdem
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
87c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
88c      REAL pks(ip1jmp1)                      ! exner au  sol
89c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
90c      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
93c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
94c      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
104c      INTEGER ij,iq,l,i,j
105      INTEGER i,j
106
107
108      real time_step, t_wrt, t_ops
109
110c      REAL rdayvrai,rdaym_ini,rday_ecri
111c      LOGICAL first
112
113      LOGICAL call_iniphys
114      data call_iniphys/.true./
115
116c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
117c+jld variables test conservation energie
118c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
119C     Tendance de la temp. potentiel d (theta)/ d t due a la
120C     tansformation d'energie cinetique en energie thermique
121C     cree par la dissipation
122c      REAL dhecdt(ip1jmp1,llm)
123c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
124c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
125c      CHARACTER (len=15) :: ztit
126c-jld
127
128
129      character (len=80) :: dynhist_file, dynhistave_file
130      character (len=20) ::modname
131      character (len=80) ::abort_message
132
133C Calendrier
134      LOGICAL true_calendar
135      PARAMETER (true_calendar = .false.)
136
137c-----------------------------------------------------------------------
138c    variables pour l'initialisation de la physique :
139c    ------------------------------------------------
140      INTEGER ngridmx,nq
141      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
142      REAL zcufi(ngridmx),zcvfi(ngridmx)
143      REAL latfi(ngridmx),lonfi(ngridmx)
144      REAL airefi(ngridmx)
145      SAVE latfi, lonfi, airefi
146     
147      INTEGER :: ierr
148
149
150c-----------------------------------------------------------------------
151c   Initialisations:
152c   ----------------
153
154      abort_message = 'last timestep reached'
155      modname = 'gcm'
156      descript = 'Run GCM LMDZ'
157      lafin    = .FALSE.
158      dynhist_file = 'dyn_hist'
159      dynhistave_file = 'dyn_hist_ave'
160
161
162c initialisation Anne
163      hadv_flg(:) = 0.
164      vadv_flg(:) = 0.
165      conv_flg(:) = 0.
166      pbl_flg(:)  = 0.
167      tracnam(:)  = '        '
168      nprath = 1
169      nbtrac = 0
170      mmt_adj(:,:,:,:) = 1
171
172
173c--------------------------------------------------------------------------
174c   Iflag_phys controle l'appel a la physique :
175c   -------------------------------------------
176c      0 : pas de physique
177c      1 : Normale (appel a phylmd, phymars ...)
178c      2 : rappel Newtonien pour la temperature + friction au sol
179      iflag_phys=1
180
181c--------------------------------------------------------------------------
182c   Lecture de l'etat initial :
183c   ---------------------------
184c     T : on lit start.nc
185c     F : le modele s'autoinitialise avec un cas academique (iniacademic)
186#ifdef CPP_IOIPSL
187      read_start=.true.
188#else
189      read_start=.false.
190#endif
191
192c-----------------------------------------------------------------------
193c   Choix du calendrier
194c   -------------------
195
196#ifdef CPP_IOIPSL
197      if (true_calendar) then
198        call ioconf_calendar('gregorian')
199      else
200        call ioconf_calendar('360d')
201      endif
202#endif
203c----------------------------------------------------------------------
204c  lecture des fichiers gcm.def ou run.def
205c  ---------------------------------------
206c
207#ifdef CPP_IOIPSL
208      CALL conf_gcm( 99, .TRUE. , clesphy0 )
209#else
210      CALL defrun( 99, .TRUE. , clesphy0 )
211#endif
212c
213c
214c------------------------------------
215c   Initialisation partie parallele
216c------------------------------------
217      CALL init_const_mpi
218
219      call init_parallel
220      call Read_Distrib
221      CALL Init_Phys_lmdz(iim,jjp1,llm,nqmx-2,mpi_size,distrib_phys)
222      CALL set_bands
223      CALL Init_interface_dyn_phys
224      CALL barrier
225
226      if (mpi_rank==0) call WriteBands
227      call SetDistrib(jj_Nb_Caldyn)
228
229c$OMP PARALLEL
230      call Init_Mod_hallo
231c$OMP END PARALLEL
232
233c$OMP PARALLEL
234      call InitComgeomphy
235c$OMP END PARALLEL
236
237      IF (config_inca /= 'none') THEN
238#ifdef INCA
239         call init_const_lmdz(
240     $        nbtrac,anneeref,dayref,
241     $        iphysiq,day_step,nday)
242
243         call init_inca_para(
244     $        iim,jjm+1,llm,klon_glo,mpi_size,
245     $        distrib_phys,COMM_LMDZ)
246#endif
247      END IF
248
249c-----------------------------------------------------------------------
250c   Initialisation des traceurs
251c   ---------------------------
252c  Choix du schema pour l'advection
253c  dans fichier trac.def ou via INCA
254
255       call iniadvtrac(nq)
256c
257c-----------------------------------------------------------------------
258c   Lecture de l'etat initial :
259c   ---------------------------
260
261c  lecture du fichier start.nc
262      if (read_start) then
263#ifdef CPP_IOIPSL
264         CALL dynetat0("start.nc",nqmx,vcov,ucov,
265     .              teta,q,masse,ps,phis, time_0)
266c       write(73,*) 'ucov',ucov
267c       write(74,*) 'vcov',vcov
268c       write(75,*) 'teta',teta
269c       write(76,*) 'ps',ps
270c       write(77,*) 'q',q
271
272#endif
273      endif
274
275c le cas echeant, creation d un etat initial
276      IF (prt_level > 9) WRITE(lunout,*)
277     .                 'AVANT iniacademic AVANT AVANT AVANT AVANT'
278      if (.not.read_start) then
279         CALL iniacademic(nqmx,vcov,ucov,teta,q,masse,ps,phis,time_0)
280      endif
281
282
283c-----------------------------------------------------------------------
284c   Lecture des parametres de controle pour la simulation :
285c   -------------------------------------------------------
286c  on recalcule eventuellement le pas de temps
287
288      IF(MOD(day_step,iperiod).NE.0) THEN
289        abort_message =
290     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
291        call abort_gcm(modname,abort_message,1)
292      ENDIF
293
294      IF(MOD(day_step,iphysiq).NE.0) THEN
295        abort_message =
296     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
297        call abort_gcm(modname,abort_message,1)
298      ENDIF
299
300      zdtvr    = daysec/FLOAT(day_step)
301        IF(dtvr.NE.zdtvr) THEN
302         WRITE(lunout,*)
303     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
304        ENDIF
305
306C
307C on remet le calendrier à zero si demande
308c
309      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
310        write(lunout,*)
311     .  ' Attention les dates initiales lues dans le fichier'
312        write(lunout,*)
313     .  ' restart ne correspondent pas a celles lues dans '
314        write(lunout,*)' gcm.def'
315        if (raz_date .ne. 1) then
316          write(lunout,*)
317     .    ' On garde les dates du fichier restart'
318        else
319          annee_ref = anneeref
320          day_ref = dayref
321          day_ini = dayref
322          itau_dyn = 0
323          itau_phy = 0
324          time_0 = 0.
325          write(lunout,*)
326     .   ' On reinitialise a la date lue dans gcm.def'
327        endif
328      ELSE
329        raz_date = 0
330      endif
331
332
333c  nombre d'etats dans les fichiers demarrage et histoire
334      nbetatdem = nday / iecri
335      nbetatmoy = nday / periodav + 1
336
337c-----------------------------------------------------------------------
338c   Initialisation des constantes dynamiques :
339c   ------------------------------------------
340      dtvr = zdtvr
341      CALL iniconst
342
343c-----------------------------------------------------------------------
344c   Initialisation de la geometrie :
345c   --------------------------------
346      CALL inigeom
347
348c-----------------------------------------------------------------------
349c   Initialisation du filtre :
350c   --------------------------
351      CALL inifilr
352c
353c-----------------------------------------------------------------------
354c   Initialisation de la dissipation :
355c   ----------------------------------
356
357      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
358     *                tetagdiv, tetagrot , tetatemp              )
359
360c-----------------------------------------------------------------------
361c   Initialisation de la physique :
362c   -------------------------------
363#ifdef CPP_PHYS
364      IF (call_iniphys.and.iflag_phys.eq.1) THEN
365         latfi(1)=rlatu(1)
366         lonfi(1)=0.
367         zcufi(1) = cu(1)
368         zcvfi(1) = cv(1)
369         DO j=2,jjm
370            DO i=1,iim
371               latfi((j-2)*iim+1+i)= rlatu(j)
372               lonfi((j-2)*iim+1+i)= rlonv(i)
373               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
374               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
375            ENDDO
376         ENDDO
377         latfi(ngridmx)= rlatu(jjp1)
378         lonfi(ngridmx)= 0.
379         zcufi(ngridmx) = cu(ip1jm+1)
380         zcvfi(ngridmx) = cv(ip1jm-iim)
381         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
382
383         WRITE(lunout,*)
384     .           'WARNING!!! vitesse verticale nulle dans la physique'
385
386         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
387     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
388
389         call_iniphys=.false.
390
391      ENDIF
392#endif
393
394
395c-----------------------------------------------------------------------
396c   Initialisation des dimensions d'INCA :
397c   --------------------------------------
398      IF (config_inca /= 'none') THEN
399!$OMP PARALLEL
400#ifdef INCA
401         CALL init_inca_dim(klon_omp,llm,iim,jjm,
402     $        rlonu,rlatu,rlonv,rlatv)
403#endif
404!$OMP END PARALLEL
405      END IF
406
407c-----------------------------------------------------------------------
408c   Initialisation des I/O :
409c   ------------------------
410
411
412      day_end = day_ini + nday
413      WRITE(lunout,300)day_ini,day_end
414
415#ifdef CPP_IOIPSL
416      CALL dynredem0_p("restart.nc", day_end, phis, nqmx)
417
418      ecripar = .TRUE.
419
420      if ( 1.eq.1) then
421      time_step = zdtvr
422      t_ops = iecri * daysec
423      t_wrt = iecri * daysec
424      CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
425     .              t_ops, t_wrt, nqmx, histid, histvid)
426
427      t_ops = iperiod * time_step
428      t_wrt = periodav * daysec
429      CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
430     .              t_ops, t_wrt, nqmx, histaveid)
431
432      dtav = iperiod*dtvr/daysec
433      endif
434
435
436#endif
437
438c  Choix des frequences de stokage pour le offline
439c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
440c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
441      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
442      istphy=istdyn/iphysiq     
443
444
445c
446c-----------------------------------------------------------------------
447c   Integration temporelle du modele :
448c   ----------------------------------
449
450c       write(78,*) 'ucov',ucov
451c       write(78,*) 'vcov',vcov
452c       write(78,*) 'teta',teta
453c       write(78,*) 'ps',ps
454c       write(78,*) 'q',q
455
456c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
457      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,nq,q,clesphy0,
458     .              time_0)
459c$OMP END PARALLEL
460
461
462 300  FORMAT('1'/,15x,'run du pas',i7,2x,'au pas',i7,2x,
463     . 'c''est a dire du jour',i7,3x,'au jour',i7//)
464      END
465
Note: See TracBrowser for help on using the repository browser.