source: LMDZ5/trunk/libf/dyn3dpar/gcm.F @ 2222

Last change on this file since 2222 was 2222, checked in by Ehouarn Millour, 9 years ago

Bug fix: Poles are single points on physics grid. True mesh area there is thus the sum of corresponding dynamics "polar meshes" areas.
Note that this also implies that some extra work is needed in physics to generate correct areas at "polar points" for outputs (especially for zoomed grids).
EM

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