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

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

Première étape de l'implémentation de XIOS. Modifications isolées dans des flags CPP_XIOS. Sorties opérationnelles (sauf stations et régionalisation) en modes séquentiel et omp, pas mpi.
UG
...........................................
First step of the XIOS implementation. Modifications are confined into CPP_XIOS flags. Output is operationnal (except for stations and regionalization) in sequential and omp modes (not mpi).
UG

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