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

Last change on this file since 2223 was 2222, checked in by Ehouarn Millour, 10 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
File size: 16.7 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      REAL zdtvr
86
87c   variables dynamiques
88      REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
89      REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle
90      REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
91      REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
92c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
93      REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
94      REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
95c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
96c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
97
98c variables dynamiques intermediaire pour le transport
99
100c   variables pour le fichier histoire
101      REAL dtav      ! intervalle de temps elementaire
102
103      REAL time_0
104
105      LOGICAL lafin
106c      INTEGER ij,iq,l,i,j
107      INTEGER i,j
108
109
110      real time_step, t_wrt, t_ops
111
112
113      LOGICAL call_iniphys
114      data call_iniphys/.true./
115
116c+jld variables test conservation energie
117c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
118C     Tendance de la temp. potentiel d (theta)/ d t due a la
119C     tansformation d'energie cinetique en energie thermique
120C     cree par la dissipation
121c      REAL dhecdt(ip1jmp1,llm)
122c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
123c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
124c      CHARACTER (len=15) :: ztit
125c-jld
126
127
128      character (len=80) :: dynhist_file, dynhistave_file
129      character (len=20) :: modname
130      character (len=80) :: abort_message
131! locales pour gestion du temps
132      INTEGER :: an, mois, jour
133      REAL :: heure
134
135
136c-----------------------------------------------------------------------
137c    variables pour l'initialisation de la physique :
138c    ------------------------------------------------
139      INTEGER ngridmx
140      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
141      REAL zcufi(ngridmx),zcvfi(ngridmx)
142      REAL latfi(ngridmx),lonfi(ngridmx)
143      REAL airefi(ngridmx)
144      SAVE latfi, lonfi, airefi
145     
146      INTEGER :: ierr
147
148
149c-----------------------------------------------------------------------
150c   Initialisations:
151c   ----------------
152
153      abort_message = 'last timestep reached'
154      modname = 'gcm'
155      descript = 'Run GCM LMDZ'
156      lafin    = .FALSE.
157      dynhist_file = 'dyn_hist'
158      dynhistave_file = 'dyn_hist_ave'
159
160
161
162c----------------------------------------------------------------------
163c  lecture des fichiers gcm.def ou run.def
164c  ---------------------------------------
165c
166! Ehouarn: dump possibility of using defrun
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      CALL init_const_mpi
176
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 Set_Distrib(distrib_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
203c-----------------------------------------------------------------------
204c   Choix du calendrier
205c   -------------------
206
207c      calend = 'earth_365d'
208
209#ifdef CPP_IOIPSL
210      if (calend == 'earth_360d') then
211        call ioconf_calendar('360d')
212        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
213      else if (calend == 'earth_365d') then
214        call ioconf_calendar('noleap')
215        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
216      else if (calend == 'earth_366d') then
217        call ioconf_calendar('gregorian')
218        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
219      else
220        abort_message = 'Mauvais choix de calendrier'
221        call abort_gcm(modname,abort_message,1)
222      endif
223#endif
224
225      IF (type_trac == 'inca') THEN
226#ifdef INCA
227         call init_const_lmdz(
228     $        nbtr,anneeref,dayref,
229     $        iphysiq,day_step,nday,
230     $        nbsrf, is_oce,is_sic,
231     $        is_ter,is_lic, calend)
232
233         call init_inca_para(
234     $        iim,jjm+1,llm,klon_glo,mpi_size,
235     $        distrib_phys,COMM_LMDZ)
236#endif
237      END IF
238
239c-----------------------------------------------------------------------
240c   Initialisation des traceurs
241c   ---------------------------
242c  Choix du nombre de traceurs et du schema pour l'advection
243c  dans fichier traceur.def, par default ou via INCA
244      call infotrac_init
245
246c Allocation de la tableau q : champs advectes   
247      ALLOCATE(ucov(ijb_u:ije_u,llm))
248      ALLOCATE(vcov(ijb_v:ije_v,llm))
249      ALLOCATE(teta(ijb_u:ije_u,llm))
250      ALLOCATE(masse(ijb_u:ije_u,llm))
251      ALLOCATE(ps(ijb_u:ije_u))
252      ALLOCATE(phis(ijb_u:ije_u))
253      ALLOCATE(q(ijb_u:ije_u,llm,nqtot))
254
255c-----------------------------------------------------------------------
256c   Lecture de l'etat initial :
257c   ---------------------------
258
259c  lecture du fichier start.nc
260      if (read_start) then
261      ! we still need to run iniacademic to initialize some
262      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
263        if (iflag_phys.ne.1) then
264          CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
265        endif
266
267!        if (planet_type.eq."earth") then
268! Load an Earth-format start file
269         CALL dynetat0_loc("start.nc",vcov,ucov,
270     &              teta,q,masse,ps,phis, time_0)
271!        endif ! of if (planet_type.eq."earth")
272
273c       write(73,*) 'ucov',ucov
274c       write(74,*) 'vcov',vcov
275c       write(75,*) 'teta',teta
276c       write(76,*) 'ps',ps
277c       write(77,*) 'q',q
278
279      endif ! of if (read_start)
280
281c le cas echeant, creation d un etat initial
282      IF (prt_level > 9) WRITE(lunout,*)
283     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
284      if (.not.read_start) then
285         CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
286      endif
287
288c-----------------------------------------------------------------------
289c   Lecture des parametres de controle pour la simulation :
290c   -------------------------------------------------------
291c  on recalcule eventuellement le pas de temps
292
293      IF(MOD(day_step,iperiod).NE.0) THEN
294        abort_message =
295     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
296        call abort_gcm(modname,abort_message,1)
297      ENDIF
298
299      IF(MOD(day_step,iphysiq).NE.0) THEN
300        abort_message =
301     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
302        call abort_gcm(modname,abort_message,1)
303      ENDIF
304
305      zdtvr    = daysec/REAL(day_step)
306        IF(dtvr.NE.zdtvr) THEN
307         WRITE(lunout,*)
308     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
309        ENDIF
310
311C
312C on remet le calendrier à zero si demande
313c
314      IF (start_time /= starttime) then
315        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
316     &,' fichier restart ne correspond pas à celle lue dans le run.def'
317        IF (raz_date == 1) then
318          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
319          start_time = starttime
320        ELSE
321          WRITE(lunout,*)'Je m''arrete'
322          CALL abort
323        ENDIF
324      ENDIF
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
392      if (iflag_phys.eq.1) then
393      ! these initialisations have already been done (via iniacademic)
394      ! if running in SW or Newtonian mode
395c-----------------------------------------------------------------------
396c   Initialisation des constantes dynamiques :
397c   ------------------------------------------
398        dtvr = zdtvr
399        CALL iniconst
400
401c-----------------------------------------------------------------------
402c   Initialisation de la geometrie :
403c   --------------------------------
404        CALL inigeom
405
406c-----------------------------------------------------------------------
407c   Initialisation du filtre :
408c   --------------------------
409        CALL inifilr
410      endif ! of if (iflag_phys.eq.1)
411c
412c-----------------------------------------------------------------------
413c   Initialisation de la dissipation :
414c   ----------------------------------
415
416      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
417     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
418
419c-----------------------------------------------------------------------
420c   Initialisation de la physique :
421c   -------------------------------
422      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
423         latfi(1)=rlatu(1)
424         lonfi(1)=0.
425         zcufi(1) = cu(1)
426         zcvfi(1) = cv(1)
427         DO j=2,jjm
428            DO i=1,iim
429               latfi((j-2)*iim+1+i)= rlatu(j)
430               lonfi((j-2)*iim+1+i)= rlonv(i)
431               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
432               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
433            ENDDO
434         ENDDO
435         latfi(ngridmx)= rlatu(jjp1)
436         lonfi(ngridmx)= 0.
437         zcufi(ngridmx) = cu(ip1jm+1)
438         zcvfi(ngridmx) = cv(ip1jm-iim)
439         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
440         ! Poles are single points on physics grid
441         airefi(1)=sum(aire(1:iim))
442         airefi(ngridmx)=sum(aire(ip1jmp1-(iim+1)+1:ip1jmp1-1))
443!         WRITE(lunout,*)
444!     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
445! Physics:
446#ifdef CPP_PHYS
447         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
448     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
449     &                iflag_phys)
450#endif
451         call_iniphys=.false.
452      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
453
454
455c-----------------------------------------------------------------------
456c   Initialisation des dimensions d'INCA :
457c   --------------------------------------
458      IF (type_trac == 'inca') THEN
459!$OMP PARALLEL
460#ifdef INCA
461         CALL init_inca_dim(klon_omp,llm,iim,jjm,
462     $        rlonu,rlatu,rlonv,rlatv)
463#endif
464!$OMP END PARALLEL
465      END IF
466
467c-----------------------------------------------------------------------
468c   Initialisation des I/O :
469c   ------------------------
470
471
472      if (nday>=0) then
473         day_end = day_ini + nday
474      else
475         day_end = day_ini - nday/day_step
476      endif
477 
478      WRITE(lunout,300)day_ini,day_end
479 300  FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
480
481#ifdef CPP_IOIPSL
482      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
483      write (lunout,301)jour, mois, an
484      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
485      write (lunout,302)jour, mois, an
486 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
487 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
488#endif
489
490!      if (planet_type.eq."earth") then
491! Write an Earth-format restart file
492        CALL dynredem0_loc("restart.nc", day_end, phis)
493!      endif
494
495      ecripar = .TRUE.
496
497#ifdef CPP_IOIPSL
498      time_step = zdtvr
499      IF (mpi_rank==0) then
500        if (ok_dyn_ins) then
501          ! initialize output file for instantaneous outputs
502          ! t_ops = iecri * daysec ! do operations every t_ops
503          t_ops =((1.0*iecri)/day_step) * daysec 
504          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
505          t_wrt = daysec ! iecri * daysec ! write output every t_wrt
506          CALL inithist(day_ref,annee_ref,time_step,
507     &                  t_ops,t_wrt)
508        endif
509
510      IF (ok_dyn_ave) THEN
511         ! initialize output file for averaged outputs
512         t_ops = iperiod * time_step ! do operations every t_ops
513         t_wrt = periodav * daysec   ! write output every t_wrt
514         CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
515        END IF
516      ENDIF
517      dtav = iperiod*dtvr/daysec
518#endif
519! #endif of #ifdef CPP_IOIPSL
520
521c  Choix des frequences de stokage pour le offline
522c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
523c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
524      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
525      istphy=istdyn/iphysiq     
526
527
528c
529c-----------------------------------------------------------------------
530c   Integration temporelle du modele :
531c   ----------------------------------
532
533c       write(78,*) 'ucov',ucov
534c       write(78,*) 'vcov',vcov
535c       write(78,*) 'teta',teta
536c       write(78,*) 'ps',ps
537c       write(78,*) 'q',q
538
539c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
540      CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
541c$OMP END PARALLEL
542
543!      OPEN(unit=5487,file='ok_lmdz',status='replace')
544!      WRITE(5487,*) 'ok_lmdz'
545!      CLOSE(5487)
546      END
547
Note: See TracBrowser for help on using the repository browser.