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

Last change on this file since 1151 was 1143, checked in by jghattas, 16 years ago

Recuperation des developpements fait uniquement sur la branche LMDZ4_V4_patches :

  • ajoute de la nouvelle flag ok_dynzon
  • ajoute du parametre aer_type
  • optimisation : isccp_cloud_types.F

+ bug pour le slab dans conf_phys.F90

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