source: LMDZ5/branches/testing/libf/dyn3dmem/gcm.F @ 1707

Last change on this file since 1707 was 1707, checked in by Laurent Fairhead, 12 years ago

Version testing basée sur la r1706


Testing release based on r1706

File size: 16.6 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
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 CPP_PHYS
23      USE mod_grid_phy_lmdz
24      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
25      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
26      USE dimphy
27      USE comgeomphy
28#endif
29      IMPLICIT NONE
30
31c      ......   Version  du 10/01/98    ..........
32
33c             avec  coordonnees  verticales hybrides
34c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
35
36c=======================================================================
37c
38c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
39c   -------
40c
41c   Objet:
42c   ------
43c
44c   GCM LMD nouvelle grille
45c
46c=======================================================================
47c
48c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
49c      et possibilite d'appeler une fonction f(y)  a derivee tangente
50c      hyperbolique a la  place de la fonction a derivee sinusoidale.
51c  ... Possibilite de choisir le schema pour l'advection de
52c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
53c
54c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
55c      Pour Van-Leer iadv=10
56c
57c-----------------------------------------------------------------------
58c   Declarations:
59c   -------------
60#include "dimensions.h"
61#include "paramet.h"
62#include "comconst.h"
63#include "comdissnew.h"
64#include "comvert.h"
65#include "comgeom.h"
66#include "logic.h"
67#include "temps.h"
68#include "ener.h"
69#include "description.h"
70#include "serre.h"
71!#include "com_io_dyn.h"
72#include "iniprint.h"
73#include "tracstoke.h"
74
75#ifdef INCA
76! Only INCA needs these informations (from the Earth's physics)
77#include "indicesol.h"
78#endif
79
80      INTEGER         longcles
81      PARAMETER     ( longcles = 20 )
82      REAL  clesphy0( longcles )
83      SAVE  clesphy0
84
85
86
87      REAL zdtvr
88
89c   variables dynamiques
90      REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
91      REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle
92      REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
93      REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
94c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
95c      REAL pks(ip1jmp1)                      ! exner au  sol
96c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
97c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
98      REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
99      REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
100c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
101c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
102
103c variables dynamiques intermediaire pour le transport
104
105c   variables pour le fichier histoire
106      REAL dtav      ! intervalle de temps elementaire
107
108      REAL time_0
109
110      LOGICAL lafin
111c      INTEGER ij,iq,l,i,j
112      INTEGER i,j
113
114
115      real time_step, t_wrt, t_ops
116
117
118      LOGICAL call_iniphys
119      data call_iniphys/.true./
120
121c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
122c+jld variables test conservation energie
123c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
124C     Tendance de la temp. potentiel d (theta)/ d t due a la
125C     tansformation d'energie cinetique en energie thermique
126C     cree par la dissipation
127c      REAL dhecdt(ip1jmp1,llm)
128c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
129c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
130c      CHARACTER (len=15) :: ztit
131c-jld
132
133
134      character (len=80) :: dynhist_file, dynhistave_file
135      character (len=20) :: modname
136      character (len=80) :: abort_message
137! locales pour gestion du temps
138      INTEGER :: an, mois, jour
139      REAL :: heure
140
141
142c-----------------------------------------------------------------------
143c    variables pour l'initialisation de la physique :
144c    ------------------------------------------------
145      INTEGER ngridmx
146      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
147      REAL zcufi(ngridmx),zcvfi(ngridmx)
148      REAL latfi(ngridmx),lonfi(ngridmx)
149      REAL airefi(ngridmx)
150      SAVE latfi, lonfi, airefi
151     
152      INTEGER :: ierr
153
154
155c-----------------------------------------------------------------------
156c   Initialisations:
157c   ----------------
158
159      abort_message = 'last timestep reached'
160      modname = 'gcm'
161      descript = 'Run GCM LMDZ'
162      lafin    = .FALSE.
163      dynhist_file = 'dyn_hist'
164      dynhistave_file = 'dyn_hist_ave'
165
166
167
168c----------------------------------------------------------------------
169c  lecture des fichiers gcm.def ou run.def
170c  ---------------------------------------
171c
172! Ehouarn: dump possibility of using defrun
173!#ifdef CPP_IOIPSL
174      CALL conf_gcm( 99, .TRUE. , clesphy0 )
175!#else
176!      CALL defrun( 99, .TRUE. , clesphy0 )
177!#endif
178c
179c
180c------------------------------------
181c   Initialisation partie parallele
182c------------------------------------
183      CALL init_const_mpi
184
185      call init_parallel
186      call ini_getparam("out.def")
187      call Read_Distrib
188
189#ifdef CPP_PHYS
190        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
191#endif
192      CALL set_bands
193#ifdef CPP_PHYS
194      CALL Init_interface_dyn_phys
195#endif
196      CALL barrier
197
198      if (mpi_rank==0) call WriteBands
199      call Set_Distrib(distrib_caldyn)
200
201c$OMP PARALLEL
202      call Init_Mod_hallo
203c$OMP END PARALLEL
204
205#ifdef CPP_PHYS
206c$OMP PARALLEL
207      call InitComgeomphy
208c$OMP END PARALLEL
209#endif
210
211c-----------------------------------------------------------------------
212c   Choix du calendrier
213c   -------------------
214
215c      calend = 'earth_365d'
216
217#ifdef CPP_IOIPSL
218      if (calend == 'earth_360d') then
219        call ioconf_calendar('360d')
220        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
221      else if (calend == 'earth_365d') then
222        call ioconf_calendar('noleap')
223        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
224      else if (calend == 'earth_366d') then
225        call ioconf_calendar('gregorian')
226        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
227      else
228        abort_message = 'Mauvais choix de calendrier'
229        call abort_gcm(modname,abort_message,1)
230      endif
231#endif
232
233      IF (type_trac == 'inca') THEN
234#ifdef INCA
235         call init_const_lmdz(
236     $        nbtr,anneeref,dayref,
237     $        iphysiq,day_step,nday,
238     $        nbsrf, is_oce,is_sic,
239     $        is_ter,is_lic)
240
241         call init_inca_para(
242     $        iim,jjm+1,llm,klon_glo,mpi_size,
243     $        distrib_phys,COMM_LMDZ)
244#endif
245      END IF
246
247c-----------------------------------------------------------------------
248c   Initialisation des traceurs
249c   ---------------------------
250c  Choix du nombre de traceurs et du schema pour l'advection
251c  dans fichier traceur.def, par default ou via INCA
252      call infotrac_init
253
254c Allocation de la tableau q : champs advectes   
255      ALLOCATE(ucov(ijb_u:ije_u,llm))
256      ALLOCATE(vcov(ijb_v:ije_v,llm))
257      ALLOCATE(teta(ijb_u:ije_u,llm))
258      ALLOCATE(masse(ijb_u:ije_u,llm))
259      ALLOCATE(ps(ijb_u:ije_u))
260      ALLOCATE(phis(ijb_u:ije_u))
261      ALLOCATE(q(ijb_u:ije_u,llm,nqtot))
262
263c-----------------------------------------------------------------------
264c   Lecture de l'etat initial :
265c   ---------------------------
266
267c  lecture du fichier start.nc
268      if (read_start) then
269      ! we still need to run iniacademic to initialize some
270      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
271        if (iflag_phys.ne.1) then
272          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
273        endif
274
275!        if (planet_type.eq."earth") then
276! Load an Earth-format start file
277         CALL dynetat0_loc("start.nc",vcov,ucov,
278     &              teta,q,masse,ps,phis, time_0)
279!        endif ! of if (planet_type.eq."earth")
280
281c       write(73,*) 'ucov',ucov
282c       write(74,*) 'vcov',vcov
283c       write(75,*) 'teta',teta
284c       write(76,*) 'ps',ps
285c       write(77,*) 'q',q
286
287      endif ! of if (read_start)
288
289c le cas echeant, creation d un etat initial
290      IF (prt_level > 9) WRITE(lunout,*)
291     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
292      if (.not.read_start) then
293         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
294      endif
295
296c-----------------------------------------------------------------------
297c   Lecture des parametres de controle pour la simulation :
298c   -------------------------------------------------------
299c  on recalcule eventuellement le pas de temps
300
301      IF(MOD(day_step,iperiod).NE.0) THEN
302        abort_message =
303     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
304        call abort_gcm(modname,abort_message,1)
305      ENDIF
306
307      IF(MOD(day_step,iphysiq).NE.0) THEN
308        abort_message =
309     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
310        call abort_gcm(modname,abort_message,1)
311      ENDIF
312
313      zdtvr    = daysec/REAL(day_step)
314        IF(dtvr.NE.zdtvr) THEN
315         WRITE(lunout,*)
316     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
317        ENDIF
318
319C
320C on remet le calendrier à zero si demande
321c
322      IF (start_time /= starttime) then
323        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
324     &,' fichier restart ne correspond pas à celle lue dans le run.def'
325        IF (raz_date == 1) then
326          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
327          start_time = starttime
328        ELSE
329          WRITE(lunout,*)'Je m''arrete'
330          CALL abort
331        ENDIF
332      ENDIF
333      IF (raz_date == 1) THEN
334        annee_ref = anneeref
335        day_ref = dayref
336        day_ini = dayref
337        itau_dyn = 0
338        itau_phy = 0
339        time_0 = 0.
340        write(lunout,*)
341     .   'GCM: On reinitialise a la date lue dans gcm.def'
342      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
343        write(lunout,*)
344     .  'GCM: Attention les dates initiales lues dans le fichier'
345        write(lunout,*)
346     .  ' restart ne correspondent pas a celles lues dans '
347        write(lunout,*)' gcm.def'
348        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
349        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
350        write(lunout,*)' Pas de remise a zero'
351      ENDIF
352c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
353c        write(lunout,*)
354c     .  'GCM: Attention les dates initiales lues dans le fichier'
355c        write(lunout,*)
356c     .  ' restart ne correspondent pas a celles lues dans '
357c        write(lunout,*)' gcm.def'
358c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
359c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
360c        if (raz_date .ne. 1) then
361c          write(lunout,*)
362c     .    'GCM: On garde les dates du fichier restart'
363c        else
364c          annee_ref = anneeref
365c          day_ref = dayref
366c          day_ini = dayref
367c          itau_dyn = 0
368c          itau_phy = 0
369c          time_0 = 0.
370c          write(lunout,*)
371c     .   'GCM: On reinitialise a la date lue dans gcm.def'
372c        endif
373c      ELSE
374c        raz_date = 0
375c      endif
376
377#ifdef CPP_IOIPSL
378      mois = 1
379      heure = 0.
380      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
381      jH_ref = jD_ref - int(jD_ref)
382      jD_ref = int(jD_ref)
383
384      call ioconf_startdate(INT(jD_ref), jH_ref)
385
386      write(lunout,*)'DEBUG'
387      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
388      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
389      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
390      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
391      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
392#else
393! Ehouarn: we still need to define JD_ref and JH_ref
394! and since we don't know how many days there are in a year
395! we set JD_ref to 0 (this should be improved ...)
396      jD_ref=0
397      jH_ref=0
398#endif
399
400
401c-----------------------------------------------------------------------
402c   Initialisation des constantes dynamiques :
403c   ------------------------------------------
404        dtvr = zdtvr
405        CALL iniconst
406
407c-----------------------------------------------------------------------
408c   Initialisation de la geometrie :
409c   --------------------------------
410        CALL inigeom
411
412c-----------------------------------------------------------------------
413c   Initialisation du filtre :
414c   --------------------------
415        CALL inifilr
416c
417c-----------------------------------------------------------------------
418c   Initialisation de la dissipation :
419c   ----------------------------------
420
421      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
422     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
423
424c-----------------------------------------------------------------------
425c   Initialisation de la physique :
426c   -------------------------------
427      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
428         latfi(1)=rlatu(1)
429         lonfi(1)=0.
430         zcufi(1) = cu(1)
431         zcvfi(1) = cv(1)
432         DO j=2,jjm
433            DO i=1,iim
434               latfi((j-2)*iim+1+i)= rlatu(j)
435               lonfi((j-2)*iim+1+i)= rlonv(i)
436               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
437               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
438            ENDDO
439         ENDDO
440         latfi(ngridmx)= rlatu(jjp1)
441         lonfi(ngridmx)= 0.
442         zcufi(ngridmx) = cu(ip1jm+1)
443         zcvfi(ngridmx) = cv(ip1jm-iim)
444         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
445
446         WRITE(lunout,*)
447     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
448! Physics:
449#ifdef CPP_PHYS
450         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
451     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
452     &                iflag_phys)
453#endif
454         call_iniphys=.false.
455      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
456
457
458c-----------------------------------------------------------------------
459c   Initialisation des dimensions d'INCA :
460c   --------------------------------------
461      IF (type_trac == 'inca') 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_loc("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_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
513        END IF
514      ENDIF
515      dtav = iperiod*dtvr/daysec
516#endif
517! #endif of #ifdef CPP_IOIPSL
518
519c  Choix des frequences de stokage pour le offline
520c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
521c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
522      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
523      istphy=istdyn/iphysiq     
524
525
526c
527c-----------------------------------------------------------------------
528c   Integration temporelle du modele :
529c   ----------------------------------
530
531c       write(78,*) 'ucov',ucov
532c       write(78,*) 'vcov',vcov
533c       write(78,*) 'teta',teta
534c       write(78,*) 'ps',ps
535c       write(78,*) 'q',q
536
537c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
538      CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
539     .              time_0)
540c$OMP END PARALLEL
541
542      OPEN(unit=5487,file='ok_lmdz',status='replace')
543      WRITE(5487,*) 'ok_lmdz'
544      CLOSE(5487)
545      END
546
Note: See TracBrowser for help on using the repository browser.