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

Last change on this file since 1856 was 1795, checked in by Ehouarn Millour, 11 years ago

Version testing basee sur la r1794


Testing release based on r1794

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