source: LMDZ5/trunk/libf/dyn3dpar/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
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 16.8 KB
Line 
1!
2! $Id: gcm.F 2171 2014-12-19 15:21:08Z jescribano $
3!
4c
5c
6      PROGRAM gcm
7
8#ifdef CPP_IOIPSL
9      USE IOIPSL
10#endif
11
12
13      USE mod_const_mpi, ONLY: init_const_mpi
14      USE parallel_lmdz
15      USE infotrac
16      USE mod_interface_dyn_phys
17      USE mod_hallo
18      USE Bands
19      USE getparam
20      USE filtreg_mod
21      USE control_mod
22
23#ifdef INCA
24! Only INCA needs these informations (from the Earth's physics)
25      USE indice_sol_mod
26#endif
27
28#ifdef CPP_PHYS
29      USE mod_grid_phy_lmdz
30      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
31      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
32      USE dimphy
33      USE comgeomphy
34#endif
35      IMPLICIT NONE
36
37c      ......   Version  du 10/01/98    ..........
38
39c             avec  coordonnees  verticales hybrides
40c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
41
42c=======================================================================
43c
44c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
45c   -------
46c
47c   Objet:
48c   ------
49c
50c   GCM LMD nouvelle grille
51c
52c=======================================================================
53c
54c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
55c      et possibilite d'appeler une fonction f(y)  a derivee tangente
56c      hyperbolique a la  place de la fonction a derivee sinusoidale.
57c  ... Possibilite de choisir le schema pour l'advection de
58c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
59c
60c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
61c      Pour Van-Leer iadv=10
62c
63c-----------------------------------------------------------------------
64c   Declarations:
65c   -------------
66#include "dimensions.h"
67#include "paramet.h"
68#include "comconst.h"
69#include "comdissnew.h"
70#include "comvert.h"
71#include "comgeom.h"
72#include "logic.h"
73#include "temps.h"
74#include "ener.h"
75#include "description.h"
76#include "serre.h"
77!#include "com_io_dyn.h"
78#include "iniprint.h"
79#include "tracstoke.h"
80
81#ifdef INCA
82! Only INCA needs these informations (from the Earth's physics)
83!#include "indicesol.h"
84#endif
85
86      INTEGER         longcles
87      PARAMETER     ( longcles = 20 )
88      REAL  clesphy0( longcles )
89      SAVE  clesphy0
90
91
92
93      REAL zdtvr
94
95c   variables dynamiques
96      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
97      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
98      REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
99      REAL ps(ip1jmp1)                       ! pression  au sol
100c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.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+jld variables test conservation energie
125c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
126C     Tendance de la temp. potentiel d (theta)/ d t due a la
127C     tansformation d'energie cinetique en energie thermique
128C     cree par la dissipation
129c      REAL dhecdt(ip1jmp1,llm)
130c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
131c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
132c      CHARACTER (len=15) :: ztit
133c-jld
134
135
136      character (len=80) :: dynhist_file, dynhistave_file
137      character (len=20) :: modname
138      character (len=80) :: abort_message
139! locales pour gestion du temps
140      INTEGER :: an, mois, jour
141      REAL :: heure
142
143
144c-----------------------------------------------------------------------
145c    variables pour l'initialisation de la physique :
146c    ------------------------------------------------
147      INTEGER ngridmx
148      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
149      REAL zcufi(ngridmx),zcvfi(ngridmx)
150      REAL latfi(ngridmx),lonfi(ngridmx)
151      REAL airefi(ngridmx)
152      SAVE latfi, lonfi, airefi
153     
154      INTEGER :: ierr
155
156
157c-----------------------------------------------------------------------
158c   Initialisations:
159c   ----------------
160
161      abort_message = 'last timestep reached'
162      modname = 'gcm'
163      descript = 'Run GCM LMDZ'
164      lafin    = .FALSE.
165      dynhist_file = 'dyn_hist'
166      dynhistave_file = 'dyn_hist_ave'
167
168
169
170c----------------------------------------------------------------------
171c  lecture des fichiers gcm.def ou run.def
172c  ---------------------------------------
173c
174! Ehouarn: dump possibility of using defrun
175!#ifdef CPP_IOIPSL
176      CALL conf_gcm( 99, .TRUE. , clesphy0 )
177      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
178     s "iphysiq must be a multiple of iperiod", 1)
179!#else
180!      CALL defrun( 99, .TRUE. , clesphy0 )
181!#endif
182c
183c
184c------------------------------------
185c   Initialisation partie parallele
186c------------------------------------
187
188      CALL init_const_mpi
189      call init_parallel
190      call ini_getparam("out.def")
191      call Read_Distrib
192
193#ifdef CPP_PHYS
194        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
195#endif
196      CALL set_bands
197#ifdef CPP_PHYS
198      CALL Init_interface_dyn_phys
199#endif
200      CALL barrier
201
202      if (mpi_rank==0) call WriteBands
203      call SetDistrib(jj_Nb_Caldyn)
204
205c$OMP PARALLEL
206      call Init_Mod_hallo
207c$OMP END PARALLEL
208
209#ifdef CPP_PHYS
210c$OMP PARALLEL
211      call InitComgeomphy
212c$OMP END PARALLEL
213#endif
214
215!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
216! Initialisation de XIOS
217!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
218
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 (type_trac == 'inca') 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 (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          call abort_gcm("gcm", "'Je m''arrete'", 1)
333        ENDIF
334      ENDIF
335      IF (raz_date == 1) THEN
336        annee_ref = anneeref
337        day_ref = dayref
338        day_ini = dayref
339        itau_dyn = 0
340        itau_phy = 0
341        time_0 = 0.
342        write(lunout,*)
343     .   'GCM: On reinitialise a la date lue dans gcm.def'
344      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
345        write(lunout,*)
346     .  'GCM: Attention les dates initiales lues dans le fichier'
347        write(lunout,*)
348     .  ' restart ne correspondent pas a celles lues dans '
349        write(lunout,*)' gcm.def'
350        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
351        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
352        write(lunout,*)' Pas de remise a zero'
353      ENDIF
354c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
355c        write(lunout,*)
356c     .  'GCM: Attention les dates initiales lues dans le fichier'
357c        write(lunout,*)
358c     .  ' restart ne correspondent pas a celles lues dans '
359c        write(lunout,*)' gcm.def'
360c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
361c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
362c        if (raz_date .ne. 1) then
363c          write(lunout,*)
364c     .    'GCM: On garde les dates du fichier restart'
365c        else
366c          annee_ref = anneeref
367c          day_ref = dayref
368c          day_ini = dayref
369c          itau_dyn = 0
370c          itau_phy = 0
371c          time_0 = 0.
372c          write(lunout,*)
373c     .   'GCM: On reinitialise a la date lue dans gcm.def'
374c        endif
375c      ELSE
376c        raz_date = 0
377c      endif
378
379#ifdef CPP_IOIPSL
380      mois = 1
381      heure = 0.
382      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
383      jH_ref = jD_ref - int(jD_ref)
384      jD_ref = int(jD_ref)
385
386      call ioconf_startdate(INT(jD_ref), jH_ref)
387
388      write(lunout,*)'DEBUG'
389      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
390      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
391      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
392      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
393      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
394#else
395! Ehouarn: we still need to define JD_ref and JH_ref
396! and since we don't know how many days there are in a year
397! we set JD_ref to 0 (this should be improved ...)
398      jD_ref=0
399      jH_ref=0
400#endif
401
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_p("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(day_ref,annee_ref,time_step,
524     &                   t_ops,t_wrt)
525!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
526!     .        t_ops, t_wrt, histaveid)
527        END IF
528      ENDIF
529      dtav = iperiod*dtvr/daysec
530#endif
531! #endif of #ifdef CPP_IOIPSL
532
533c  Choix des frequences de stokage pour le offline
534c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
535c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
536      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
537      istphy=istdyn/iphysiq     
538
539
540c
541c-----------------------------------------------------------------------
542c   Integration temporelle du modele :
543c   ----------------------------------
544
545c       write(78,*) 'ucov',ucov
546c       write(78,*) 'vcov',vcov
547c       write(78,*) 'teta',teta
548c       write(78,*) 'ps',ps
549c       write(78,*) 'q',q
550
551c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
552      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
553     .              time_0)
554c$OMP END PARALLEL
555
556
557      END
558
Note: See TracBrowser for help on using the repository browser.