source: LMDZ5/trunk/libf/dyn3dmem/gcm.F @ 2174

Last change on this file since 2174 was 2171, checked in by acozic, 10 years ago

There are some commits that we must not do just before holiday .... so be back to rev 2168

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 16.8 KB
Line 
1!
2! $Id$
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_lmdz
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#ifdef INCA
23! Only INCA needs these informations (from the Earth's physics)
24      USE indice_sol_mod
25#endif
26
27#ifdef CPP_PHYS
28      USE mod_grid_phy_lmdz
29      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
30      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
31      USE dimphy
32      USE comgeomphy
33#endif
34      IMPLICIT NONE
35
36c      ......   Version  du 10/01/98    ..........
37
38c             avec  coordonnees  verticales hybrides
39c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
40
41c=======================================================================
42c
43c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
44c   -------
45c
46c   Objet:
47c   ------
48c
49c   GCM LMD nouvelle grille
50c
51c=======================================================================
52c
53c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
54c      et possibilite d'appeler une fonction f(y)  a derivee tangente
55c      hyperbolique a la  place de la fonction a derivee sinusoidale.
56c  ... Possibilite de choisir le schema pour l'advection de
57c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
58c
59c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
60c      Pour Van-Leer iadv=10
61c
62c-----------------------------------------------------------------------
63c   Declarations:
64c   -------------
65#include "dimensions.h"
66#include "paramet.h"
67#include "comconst.h"
68#include "comdissnew.h"
69#include "comvert.h"
70#include "comgeom.h"
71#include "logic.h"
72#include "temps.h"
73#include "ener.h"
74#include "description.h"
75#include "serre.h"
76!#include "com_io_dyn.h"
77#include "iniprint.h"
78#include "tracstoke.h"
79
80#ifdef INCA
81! Only INCA needs these informations (from the Earth's physics)
82!#include "indicesol.h"
83#endif
84
85      INTEGER         longcles
86      PARAMETER     ( longcles = 20 )
87      REAL  clesphy0( longcles )
88      SAVE  clesphy0
89
90
91
92      REAL zdtvr
93
94c   variables dynamiques
95      REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
96      REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle
97      REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
98      REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
99c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
100      REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
101      REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
102c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
103c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
104
105c variables dynamiques intermediaire pour le transport
106
107c   variables pour le fichier histoire
108      REAL dtav      ! intervalle de temps elementaire
109
110      REAL time_0
111
112      LOGICAL lafin
113c      INTEGER ij,iq,l,i,j
114      INTEGER i,j
115
116
117      real time_step, t_wrt, t_ops
118
119
120      LOGICAL call_iniphys
121      data call_iniphys/.true./
122
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
128c      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
131c      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     
153      INTEGER :: ierr
154
155
156c-----------------------------------------------------------------------
157c   Initialisations:
158c   ----------------
159
160      abort_message = 'last timestep reached'
161      modname = 'gcm'
162      descript = 'Run GCM LMDZ'
163      lafin    = .FALSE.
164      dynhist_file = 'dyn_hist'
165      dynhistave_file = 'dyn_hist_ave'
166
167
168
169c----------------------------------------------------------------------
170c  lecture des fichiers gcm.def ou run.def
171c  ---------------------------------------
172c
173! Ehouarn: dump possibility of using defrun
174!#ifdef CPP_IOIPSL
175      CALL conf_gcm( 99, .TRUE. , clesphy0 )
176      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
177     s "iphysiq must be a multiple of iperiod", 1)
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
192#ifdef CPP_PHYS
193        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
194#endif
195      CALL set_bands
196#ifdef CPP_PHYS
197      CALL Init_interface_dyn_phys
198#endif
199      CALL barrier
200
201      if (mpi_rank==0) call WriteBands
202      call Set_Distrib(distrib_caldyn)
203
204c$OMP PARALLEL
205      call Init_Mod_hallo
206c$OMP END PARALLEL
207
208#ifdef CPP_PHYS
209c$OMP PARALLEL
210      call InitComgeomphy
211c$OMP END PARALLEL
212#endif
213
214c-----------------------------------------------------------------------
215c   Choix du calendrier
216c   -------------------
217
218c      calend = 'earth_365d'
219
220#ifdef CPP_IOIPSL
221      if (calend == 'earth_360d') then
222        call ioconf_calendar('360d')
223        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
224      else if (calend == 'earth_365d') then
225        call ioconf_calendar('noleap')
226        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
227      else if (calend == 'earth_366d') then
228        call ioconf_calendar('gregorian')
229        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
230      else
231        abort_message = 'Mauvais choix de calendrier'
232        call abort_gcm(modname,abort_message,1)
233      endif
234#endif
235
236      IF (type_trac == 'inca') THEN
237#ifdef INCA
238         call init_const_lmdz(
239     $        nbtr,anneeref,dayref,
240     $        iphysiq,day_step,nday,
241     $        nbsrf, is_oce,is_sic,
242     $        is_ter,is_lic)
243
244         call init_inca_para(
245     $        iim,jjm+1,llm,klon_glo,mpi_size,
246     $        distrib_phys,COMM_LMDZ)
247#endif
248      END IF
249
250c-----------------------------------------------------------------------
251c   Initialisation des traceurs
252c   ---------------------------
253c  Choix du nombre de traceurs et du schema pour l'advection
254c  dans fichier traceur.def, par default ou via INCA
255      call infotrac_init
256
257c Allocation de la tableau q : champs advectes   
258      ALLOCATE(ucov(ijb_u:ije_u,llm))
259      ALLOCATE(vcov(ijb_v:ije_v,llm))
260      ALLOCATE(teta(ijb_u:ije_u,llm))
261      ALLOCATE(masse(ijb_u:ije_u,llm))
262      ALLOCATE(ps(ijb_u:ije_u))
263      ALLOCATE(phis(ijb_u:ije_u))
264      ALLOCATE(q(ijb_u:ije_u,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_loc(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_loc("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_loc(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 (start_time /= starttime) then
326        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
327     &,' fichier restart ne correspond pas à celle lue dans le run.def'
328        IF (raz_date == 1) then
329          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
330          start_time = starttime
331        ELSE
332          WRITE(lunout,*)'Je m''arrete'
333          CALL abort
334        ENDIF
335      ENDIF
336      IF (raz_date == 1) THEN
337        annee_ref = anneeref
338        day_ref = dayref
339        day_ini = dayref
340        itau_dyn = 0
341        itau_phy = 0
342        time_0 = 0.
343        write(lunout,*)
344     .   'GCM: On reinitialise a la date lue dans gcm.def'
345      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
346        write(lunout,*)
347     .  'GCM: Attention les dates initiales lues dans le fichier'
348        write(lunout,*)
349     .  ' restart ne correspondent pas a celles lues dans '
350        write(lunout,*)' gcm.def'
351        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
352        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
353        write(lunout,*)' Pas de remise a zero'
354      ENDIF
355c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
356c        write(lunout,*)
357c     .  'GCM: Attention les dates initiales lues dans le fichier'
358c        write(lunout,*)
359c     .  ' restart ne correspondent pas a celles lues dans '
360c        write(lunout,*)' gcm.def'
361c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
362c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
363c        if (raz_date .ne. 1) then
364c          write(lunout,*)
365c     .    'GCM: On garde les dates du fichier restart'
366c        else
367c          annee_ref = anneeref
368c          day_ref = dayref
369c          day_ini = dayref
370c          itau_dyn = 0
371c          itau_phy = 0
372c          time_0 = 0.
373c          write(lunout,*)
374c     .   'GCM: On reinitialise a la date lue dans gcm.def'
375c        endif
376c      ELSE
377c        raz_date = 0
378c      endif
379
380#ifdef CPP_IOIPSL
381      mois = 1
382      heure = 0.
383      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
384      jH_ref = jD_ref - int(jD_ref)
385      jD_ref = int(jD_ref)
386
387      call ioconf_startdate(INT(jD_ref), jH_ref)
388
389      write(lunout,*)'DEBUG'
390      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
391      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
392      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
393      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
394      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
395#else
396! Ehouarn: we still need to define JD_ref and JH_ref
397! and since we don't know how many days there are in a year
398! we set JD_ref to 0 (this should be improved ...)
399      jD_ref=0
400      jH_ref=0
401#endif
402
403      if (iflag_phys.eq.1) then
404      ! these initialisations have already been done (via iniacademic)
405      ! if running in SW or Newtonian mode
406c-----------------------------------------------------------------------
407c   Initialisation des constantes dynamiques :
408c   ------------------------------------------
409        dtvr = zdtvr
410        CALL iniconst
411
412c-----------------------------------------------------------------------
413c   Initialisation de la geometrie :
414c   --------------------------------
415        CALL inigeom
416
417c-----------------------------------------------------------------------
418c   Initialisation du filtre :
419c   --------------------------
420        CALL inifilr
421      endif ! of if (iflag_phys.eq.1)
422c
423c-----------------------------------------------------------------------
424c   Initialisation de la dissipation :
425c   ----------------------------------
426
427      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
428     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
429
430c-----------------------------------------------------------------------
431c   Initialisation de la physique :
432c   -------------------------------
433      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
434         latfi(1)=rlatu(1)
435         lonfi(1)=0.
436         zcufi(1) = cu(1)
437         zcvfi(1) = cv(1)
438         DO j=2,jjm
439            DO i=1,iim
440               latfi((j-2)*iim+1+i)= rlatu(j)
441               lonfi((j-2)*iim+1+i)= rlonv(i)
442               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
443               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
444            ENDDO
445         ENDDO
446         latfi(ngridmx)= rlatu(jjp1)
447         lonfi(ngridmx)= 0.
448         zcufi(ngridmx) = cu(ip1jm+1)
449         zcvfi(ngridmx) = cv(ip1jm-iim)
450         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
451
452         WRITE(lunout,*)
453     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
454! Physics:
455#ifdef CPP_PHYS
456         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
457     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
458     &                iflag_phys)
459#endif
460         call_iniphys=.false.
461      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
462
463
464c-----------------------------------------------------------------------
465c   Initialisation des dimensions d'INCA :
466c   --------------------------------------
467      IF (type_trac == 'inca') THEN
468!$OMP PARALLEL
469#ifdef INCA
470         CALL init_inca_dim(klon_omp,llm,iim,jjm,
471     $        rlonu,rlatu,rlonv,rlatv)
472#endif
473!$OMP END PARALLEL
474      END IF
475
476c-----------------------------------------------------------------------
477c   Initialisation des I/O :
478c   ------------------------
479
480
481      if (nday>=0) then
482         day_end = day_ini + nday
483      else
484         day_end = day_ini - nday/day_step
485      endif
486 
487      WRITE(lunout,300)day_ini,day_end
488 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
489
490#ifdef CPP_IOIPSL
491      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
492      write (lunout,301)jour, mois, an
493      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
494      write (lunout,302)jour, mois, an
495 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
496 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
497#endif
498
499!      if (planet_type.eq."earth") then
500! Write an Earth-format restart file
501        CALL dynredem0_loc("restart.nc", day_end, phis)
502!      endif
503
504      ecripar = .TRUE.
505
506#ifdef CPP_IOIPSL
507      time_step = zdtvr
508      IF (mpi_rank==0) then
509        if (ok_dyn_ins) then
510          ! initialize output file for instantaneous outputs
511          ! t_ops = iecri * daysec ! do operations every t_ops
512          t_ops =((1.0*iecri)/day_step) * daysec 
513          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
514          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
515          CALL inithist(day_ref,annee_ref,time_step,
516     &                  t_ops,t_wrt)
517        endif
518
519      IF (ok_dyn_ave) THEN
520         ! initialize output file for averaged outputs
521         t_ops = iperiod * time_step ! do operations every t_ops
522         t_wrt = periodav * daysec   ! write output every t_wrt
523         CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
524        END IF
525      ENDIF
526      dtav = iperiod*dtvr/daysec
527#endif
528! #endif of #ifdef CPP_IOIPSL
529
530c  Choix des frequences de stokage pour le offline
531c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
532c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
533      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
534      istphy=istdyn/iphysiq     
535
536
537c
538c-----------------------------------------------------------------------
539c   Integration temporelle du modele :
540c   ----------------------------------
541
542c       write(78,*) 'ucov',ucov
543c       write(78,*) 'vcov',vcov
544c       write(78,*) 'teta',teta
545c       write(78,*) 'ps',ps
546c       write(78,*) 'q',q
547
548c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
549      CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
550     .              time_0)
551c$OMP END PARALLEL
552
553!      OPEN(unit=5487,file='ok_lmdz',status='replace')
554!      WRITE(5487,*) 'ok_lmdz'
555!      CLOSE(5487)
556      END
557
Note: See TracBrowser for help on using the repository browser.