source: LMDZ4/branches/LMDZ4V5.0-LF/libf/dyn3d/gcm.F @ 3603

Last change on this file since 3603 was 1316, checked in by acozic, 15 years ago

acozic: add argument in init_const_lmdz call for initialisation of INCA model

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