source: LMDZ5/branches/LMDZ5V1.0-dev/libf/dyn3dpar/gcm.F @ 3998

Last change on this file since 3998 was 1446, checked in by Ehouarn Millour, 14 years ago

Implemented modifications to enable running with only one tracer for planet types different from "earth". Rem: If flag 'planet_type' is set to "earth" (default behaviour) then there must be at least 2 tracers for the dynamics to function properly.

These updates do not induce any changes in model outputs with respect to previous revisions.

EM

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.7 KB
Line 
1!
2! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $
3!
4c
5c
6      PROGRAM gcm
7
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#endif
11
12      USE mod_const_mpi, ONLY: init_const_mpi
13      USE parallel
14      USE infotrac
15      USE mod_interface_dyn_phys
16      USE mod_hallo
17      USE Bands
18      USE getparam
19      USE filtreg_mod
20      USE control_mod
21
22! Ehouarn: for now these only apply to Earth:
23#ifdef CPP_EARTH
24      USE mod_grid_phy_lmdz
25      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
26      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
27      USE dimphy
28      USE comgeomphy
29#endif
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#include "dimensions.h"
62#include "paramet.h"
63#include "comconst.h"
64#include "comdissnew.h"
65#include "comvert.h"
66#include "comgeom.h"
67#include "logic.h"
68#include "temps.h"
69#include "ener.h"
70#include "description.h"
71#include "serre.h"
72!#include "com_io_dyn.h"
73#include "iniprint.h"
74#include "tracstoke.h"
75
76#ifdef INCA
77! Only INCA needs these informations (from the Earth's physics)
78#include "indicesol.h"
79#endif
80
81      INTEGER         longcles
82      PARAMETER     ( longcles = 20 )
83      REAL  clesphy0( longcles )
84      SAVE  clesphy0
85
86
87
88      REAL zdtvr
89c      INTEGER nbetatmoy, nbetatdem,nbetat
90      INTEGER nbetatmoy, nbetatdem
91
92c   variables dynamiques
93      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
94      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
95      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
96      REAL ps(ip1jmp1)                       ! pression  au sol
97c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
98c      REAL pks(ip1jmp1)                      ! exner au  sol
99c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
100c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
101      REAL masse(ip1jmp1,llm)                ! masse d'air
102      REAL phis(ip1jmp1)                     ! geopotentiel au sol
103c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
104c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
105
106c variables dynamiques intermediaire pour le transport
107
108c   variables pour le fichier histoire
109      REAL dtav      ! intervalle de temps elementaire
110
111      REAL time_0
112
113      LOGICAL lafin
114c      INTEGER ij,iq,l,i,j
115      INTEGER i,j
116
117
118      real time_step, t_wrt, t_ops
119
120
121      LOGICAL call_iniphys
122      data call_iniphys/.true./
123
124c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
125c+jld variables test conservation energie
126c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
127C     Tendance de la temp. potentiel d (theta)/ d t due a la
128C     tansformation d'energie cinetique en energie thermique
129C     cree par la dissipation
130c      REAL dhecdt(ip1jmp1,llm)
131c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
132c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
133c      CHARACTER (len=15) :: ztit
134c-jld
135
136
137      character (len=80) :: dynhist_file, dynhistave_file
138      character (len=20) :: modname
139      character (len=80) :: abort_message
140! locales pour gestion du temps
141      INTEGER :: an, mois, jour
142      REAL :: heure
143
144
145c-----------------------------------------------------------------------
146c    variables pour l'initialisation de la physique :
147c    ------------------------------------------------
148      INTEGER ngridmx
149      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
150      REAL zcufi(ngridmx),zcvfi(ngridmx)
151      REAL latfi(ngridmx),lonfi(ngridmx)
152      REAL airefi(ngridmx)
153      SAVE latfi, lonfi, airefi
154     
155      INTEGER :: ierr
156
157
158c-----------------------------------------------------------------------
159c   Initialisations:
160c   ----------------
161
162      abort_message = 'last timestep reached'
163      modname = 'gcm'
164      descript = 'Run GCM LMDZ'
165      lafin    = .FALSE.
166      dynhist_file = 'dyn_hist'
167      dynhistave_file = 'dyn_hist_ave'
168
169
170
171c----------------------------------------------------------------------
172c  lecture des fichiers gcm.def ou run.def
173c  ---------------------------------------
174c
175! Ehouarn: dump possibility of using defrun
176!#ifdef CPP_IOIPSL
177      CALL conf_gcm( 99, .TRUE. , clesphy0 )
178!#else
179!      CALL defrun( 99, .TRUE. , clesphy0 )
180!#endif
181c
182c
183c------------------------------------
184c   Initialisation partie parallele
185c------------------------------------
186      CALL init_const_mpi
187
188      call init_parallel
189      call ini_getparam("out.def")
190      call Read_Distrib
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,mpi_size,distrib_phys)
195#endif
196      endif ! of if (planet_type.eq."earth")
197      CALL set_bands
198#ifdef CPP_EARTH
199! Ehouarn: For now only Earth physics is parallel
200      CALL Init_interface_dyn_phys
201#endif
202      CALL barrier
203
204      if (mpi_rank==0) call WriteBands
205      call SetDistrib(jj_Nb_Caldyn)
206
207c$OMP PARALLEL
208      call Init_Mod_hallo
209c$OMP END PARALLEL
210
211! Ehouarn : temporarily (?) keep this only for Earth
212      if (planet_type.eq."earth") then
213#ifdef CPP_EARTH
214c$OMP PARALLEL
215      call InitComgeomphy
216c$OMP END PARALLEL
217#endif
218      endif ! of if (planet_type.eq."earth")
219
220c-----------------------------------------------------------------------
221c   Choix du calendrier
222c   -------------------
223
224c      calend = 'earth_365d'
225
226#ifdef CPP_IOIPSL
227      if (calend == 'earth_360d') then
228        call ioconf_calendar('360d')
229        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
230      else if (calend == 'earth_365d') then
231        call ioconf_calendar('noleap')
232        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
233      else if (calend == 'earth_366d') then
234        call ioconf_calendar('gregorian')
235        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
236      else
237        abort_message = 'Mauvais choix de calendrier'
238        call abort_gcm(modname,abort_message,1)
239      endif
240#endif
241
242      IF (config_inca /= 'none') THEN
243#ifdef INCA
244         call init_const_lmdz(
245     $        nbtr,anneeref,dayref,
246     $        iphysiq,day_step,nday,
247     $        nbsrf, is_oce,is_sic,
248     $        is_ter,is_lic)
249
250         call init_inca_para(
251     $        iim,jjm+1,llm,klon_glo,mpi_size,
252     $        distrib_phys,COMM_LMDZ)
253#endif
254      END IF
255
256c-----------------------------------------------------------------------
257c   Initialisation des traceurs
258c   ---------------------------
259c  Choix du nombre de traceurs et du schema pour l'advection
260c  dans fichier traceur.def, par default ou via INCA
261      call infotrac_init
262
263c Allocation de la tableau q : champs advectes   
264      ALLOCATE(q(ip1jmp1,llm,nqtot))
265
266c-----------------------------------------------------------------------
267c   Lecture de l'etat initial :
268c   ---------------------------
269
270c  lecture du fichier start.nc
271      if (read_start) then
272      ! we still need to run iniacademic to initialize some
273      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
274        if (iflag_phys.ne.1) then
275          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
276        endif
277
278!        if (planet_type.eq."earth") then
279! Load an Earth-format start file
280         CALL dynetat0("start.nc",vcov,ucov,
281     &              teta,q,masse,ps,phis, time_0)
282!        endif ! of if (planet_type.eq."earth")
283
284c       write(73,*) 'ucov',ucov
285c       write(74,*) 'vcov',vcov
286c       write(75,*) 'teta',teta
287c       write(76,*) 'ps',ps
288c       write(77,*) 'q',q
289
290      endif ! of if (read_start)
291
292c le cas echeant, creation d un etat initial
293      IF (prt_level > 9) WRITE(lunout,*)
294     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
295      if (.not.read_start) then
296         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
297      endif
298
299c-----------------------------------------------------------------------
300c   Lecture des parametres de controle pour la simulation :
301c   -------------------------------------------------------
302c  on recalcule eventuellement le pas de temps
303
304      IF(MOD(day_step,iperiod).NE.0) THEN
305        abort_message =
306     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
307        call abort_gcm(modname,abort_message,1)
308      ENDIF
309
310      IF(MOD(day_step,iphysiq).NE.0) THEN
311        abort_message =
312     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
313        call abort_gcm(modname,abort_message,1)
314      ENDIF
315
316      zdtvr    = daysec/REAL(day_step)
317        IF(dtvr.NE.zdtvr) THEN
318         WRITE(lunout,*)
319     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
320        ENDIF
321
322C
323C on remet le calendrier à zero si demande
324c
325      IF (raz_date == 1) THEN
326        annee_ref = anneeref
327        day_ref = dayref
328        day_ini = dayref
329        itau_dyn = 0
330        itau_phy = 0
331        time_0 = 0.
332        write(lunout,*)
333     .   'GCM: On reinitialise a la date lue dans gcm.def'
334      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
335        write(lunout,*)
336     .  'GCM: Attention les dates initiales lues dans le fichier'
337        write(lunout,*)
338     .  ' restart ne correspondent pas a celles lues dans '
339        write(lunout,*)' gcm.def'
340        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
341        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
342        write(lunout,*)' Pas de remise a zero'
343      ENDIF
344c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
345c        write(lunout,*)
346c     .  'GCM: Attention les dates initiales lues dans le fichier'
347c        write(lunout,*)
348c     .  ' restart ne correspondent pas a celles lues dans '
349c        write(lunout,*)' gcm.def'
350c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
351c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
352c        if (raz_date .ne. 1) then
353c          write(lunout,*)
354c     .    'GCM: On garde les dates du fichier restart'
355c        else
356c          annee_ref = anneeref
357c          day_ref = dayref
358c          day_ini = dayref
359c          itau_dyn = 0
360c          itau_phy = 0
361c          time_0 = 0.
362c          write(lunout,*)
363c     .   'GCM: On reinitialise a la date lue dans gcm.def'
364c        endif
365c      ELSE
366c        raz_date = 0
367c      endif
368
369#ifdef CPP_IOIPSL
370      mois = 1
371      heure = 0.
372      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
373      jH_ref = jD_ref - int(jD_ref)
374      jD_ref = int(jD_ref)
375
376      call ioconf_startdate(INT(jD_ref), jH_ref)
377
378      write(lunout,*)'DEBUG'
379      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
380      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
381      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
382      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
383      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
384#else
385! Ehouarn: we still need to define JD_ref and JH_ref
386! and since we don't know how many days there are in a year
387! we set JD_ref to 0 (this should be improved ...)
388      jD_ref=0
389      jH_ref=0
390#endif
391
392c  nombre d'etats dans les fichiers demarrage et histoire
393      nbetatdem = nday / iecri
394      nbetatmoy = nday / periodav + 1
395
396      if (iflag_phys.eq.1) then
397      ! these initialisations have already been done (via iniacademic)
398      ! if running in SW or Newtonian mode
399c-----------------------------------------------------------------------
400c   Initialisation des constantes dynamiques :
401c   ------------------------------------------
402        dtvr = zdtvr
403        CALL iniconst
404
405c-----------------------------------------------------------------------
406c   Initialisation de la geometrie :
407c   --------------------------------
408        CALL inigeom
409
410c-----------------------------------------------------------------------
411c   Initialisation du filtre :
412c   --------------------------
413        CALL inifilr
414      endif ! of if (iflag_phys.eq.1)
415c
416c-----------------------------------------------------------------------
417c   Initialisation de la dissipation :
418c   ----------------------------------
419
420      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
421     *                tetagdiv, tetagrot , tetatemp              )
422
423c-----------------------------------------------------------------------
424c   Initialisation de la physique :
425c   -------------------------------
426      IF (call_iniphys.and.iflag_phys.eq.1) THEN
427         latfi(1)=rlatu(1)
428         lonfi(1)=0.
429         zcufi(1) = cu(1)
430         zcvfi(1) = cv(1)
431         DO j=2,jjm
432            DO i=1,iim
433               latfi((j-2)*iim+1+i)= rlatu(j)
434               lonfi((j-2)*iim+1+i)= rlonv(i)
435               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
436               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
437            ENDDO
438         ENDDO
439         latfi(ngridmx)= rlatu(jjp1)
440         lonfi(ngridmx)= 0.
441         zcufi(ngridmx) = cu(ip1jm+1)
442         zcvfi(ngridmx) = cv(ip1jm-iim)
443         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
444
445         WRITE(lunout,*)
446     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
447! Earth:
448         if (planet_type.eq."earth") then
449#ifdef CPP_EARTH
450         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
451     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
452#endif
453         endif ! of if (planet_type.eq."earth")
454         call_iniphys=.false.
455      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
456
457
458c-----------------------------------------------------------------------
459c   Initialisation des dimensions d'INCA :
460c   --------------------------------------
461      IF (config_inca /= 'none') THEN
462!$OMP PARALLEL
463#ifdef INCA
464         CALL init_inca_dim(klon_omp,llm,iim,jjm,
465     $        rlonu,rlatu,rlonv,rlatv)
466#endif
467!$OMP END PARALLEL
468      END IF
469
470c-----------------------------------------------------------------------
471c   Initialisation des I/O :
472c   ------------------------
473
474
475      day_end = day_ini + nday
476      WRITE(lunout,300)day_ini,day_end
477 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
478
479#ifdef CPP_IOIPSL
480      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
481      write (lunout,301)jour, mois, an
482      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
483      write (lunout,302)jour, mois, an
484 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
485 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
486#endif
487
488!      if (planet_type.eq."earth") then
489! Write an Earth-format restart file
490        CALL dynredem0_p("restart.nc", day_end, phis)
491!      endif
492
493      ecripar = .TRUE.
494
495#ifdef CPP_IOIPSL
496      time_step = zdtvr
497      IF (mpi_rank==0) then
498        if (ok_dyn_ins) then
499          ! initialize output file for instantaneous outputs
500          ! t_ops = iecri * daysec ! do operations every t_ops
501          t_ops =((1.0*iecri)/day_step) * daysec 
502          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
503          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
504          CALL inithist(day_ref,annee_ref,time_step,
505     &                  t_ops,t_wrt)
506        endif
507
508        IF (ok_dyn_ave) THEN
509          ! initialize output file for averaged outputs
510          t_ops = iperiod * time_step ! do operations every t_ops
511          t_wrt = periodav * daysec   ! write output every t_wrt
512          CALL initdynav(day_ref,annee_ref,time_step,
513     &                   t_ops,t_wrt)
514!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
515!     .        t_ops, t_wrt, histaveid)
516        END IF
517      ENDIF
518      dtav = iperiod*dtvr/daysec
519#endif
520! #endif of #ifdef CPP_IOIPSL
521
522c  Choix des frequences de stokage pour le offline
523c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
524c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
525      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
526      istphy=istdyn/iphysiq     
527
528
529c
530c-----------------------------------------------------------------------
531c   Integration temporelle du modele :
532c   ----------------------------------
533
534c       write(78,*) 'ucov',ucov
535c       write(78,*) 'vcov',vcov
536c       write(78,*) 'teta',teta
537c       write(78,*) 'ps',ps
538c       write(78,*) 'q',q
539
540c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
541      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
542     .              time_0)
543c$OMP END PARALLEL
544
545
546      END
547
Note: See TracBrowser for help on using the repository browser.