source: trunk/LMDZ.COMMON/libf/dyn3dpar/gcm.F @ 965

Last change on this file since 965 was 965, checked in by emillour, 12 years ago

Common dynamics and generic/universal GCM:

  • LMDZ.COMMON: minor bug fix on the computation of physics mesh area in gcm.F
  • LMDZ.UNIVERSAL: missing clean initialization of tab_cntrl(:) array in phyredem.F90
  • LMDZ.GENERIC: minor bug fix in hydrol.F90, only output runoff if it is used. Update output routines so that all outputs files (stats, diagfi.nc, diagsoil.nc, diagspecIR.nc and diagspecVI.nc) can be generated when running LMDZ.UNIVERSAL in MPI mode.

EM

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