source: trunk/libf/dyn3dpar/gcm.F @ 98

Last change on this file since 98 was 97, checked in by slebonnois, 14 years ago

Serie de modifs SL pour homogeneisation des phytitan et phyvenus
Ca touche aussi aux liens phy/dyn (surtout a propos de clesphy0),
a verifier avec les autres, donc...

File size: 18.1 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: for now these only apply to Earth:
23#ifdef CPP_EARTH
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
84c      INTEGER nbetatmoy, nbetatdem,nbetat
85      INTEGER nbetatmoy, nbetatdem
86
87c   variables dynamiques
88      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
89      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
90      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
91      REAL ps(ip1jmp1)                       ! pression  au sol
92c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
93c      REAL pks(ip1jmp1)                      ! exner au  sol
94c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
95c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
96      REAL masse(ip1jmp1,llm)                ! masse d'air
97      REAL phis(ip1jmp1)                     ! geopotentiel au sol
98c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
99c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
100
101c variables dynamiques intermediaire pour le transport
102
103c   variables pour le fichier histoire
104      REAL dtav      ! intervalle de temps elementaire
105
106      REAL time_0
107
108      LOGICAL lafin
109c      INTEGER ij,iq,l,i,j
110      INTEGER i,j
111
112
113      real time_step, t_wrt, t_ops
114
115
116      LOGICAL call_iniphys
117      data call_iniphys/.true./
118
119c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
120c+jld variables test conservation energie
121c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
122C     Tendance de la temp. potentiel d (theta)/ d t due a la
123C     tansformation d'energie cinetique en energie thermique
124C     cree par la dissipation
125c      REAL dhecdt(ip1jmp1,llm)
126c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
127c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
128c      CHARACTER (len=15) :: ztit
129c-jld
130
131
132      character (len=80) :: dynhist_file, dynhistave_file
133      character (len=20) :: modname
134      character (len=80) :: abort_message
135! locales pour gestion du temps
136      INTEGER :: an, mois, jour
137      REAL :: heure
138
139
140c-----------------------------------------------------------------------
141c    variables pour l'initialisation de la physique :
142c    ------------------------------------------------
143      INTEGER ngridmx
144      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
145      REAL zcufi(ngridmx),zcvfi(ngridmx)
146      REAL latfi(ngridmx),lonfi(ngridmx)
147      REAL airefi(ngridmx)
148      SAVE latfi, lonfi, airefi
149     
150      INTEGER :: ierr
151
152
153c-----------------------------------------------------------------------
154c   Initialisations:
155c   ----------------
156
157      abort_message = 'last timestep reached'
158      modname = 'gcm'
159      descript = 'Run GCM LMDZ'
160      lafin    = .FALSE.
161      dynhist_file = 'dyn_hist'
162      dynhistave_file = 'dyn_hist_ave'
163
164
165
166c----------------------------------------------------------------------
167c  lecture des fichiers gcm.def ou run.def
168c  ---------------------------------------
169c
170! Ehouarn: dump possibility of using defrun
171!#ifdef CPP_IOIPSL
172      CALL conf_gcm( 99, .TRUE. )
173!#else
174!      CALL defrun( 99, .TRUE. , clesphy0 )
175!#endif
176c
177c
178c------------------------------------
179c   Initialisation partie parallele
180c------------------------------------
181      CALL init_const_mpi
182
183      call init_parallel
184      call ini_getparam("out.def")
185      call Read_Distrib
186! Ehouarn : temporarily (?) keep this only for Earth
187!      if (planet_type.eq."earth") then
188!#ifdef CPP_EARTH
189#ifdef CPP_PHYS
190        CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
191#endif
192!      endif ! of if (planet_type.eq."earth")
193      CALL set_bands
194#ifdef CPP_PHYS
195! Ehouarn: NB: For now only Earth physics is parallel
196      CALL Init_interface_dyn_phys
197#endif
198      CALL barrier
199
200      if (mpi_rank==0) call WriteBands
201      call SetDistrib(jj_Nb_Caldyn)
202
203c$OMP PARALLEL
204      call Init_Mod_hallo
205c$OMP END PARALLEL
206
207! Ehouarn : temporarily (?) keep this only for Earth
208!      if (planet_type.eq."earth") then
209!#ifdef CPP_EARTH
210#ifdef CPP_PHYS
211c$OMP PARALLEL
212      call InitComgeomphy
213c$OMP END PARALLEL
214#endif
215!      endif ! of if (planet_type.eq."earth")
216!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
217c
218c Initialisations pour Cp(T) Venus
219      call ini_cpdet
220c
221c-----------------------------------------------------------------------
222c   Choix du calendrier
223c   -------------------
224
225c      calend = 'earth_365d'
226
227#ifdef CPP_IOIPSL
228      if (calend == 'earth_360d') then
229        call ioconf_calendar('360d')
230        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
231      else if (calend == 'earth_365d') then
232        call ioconf_calendar('noleap')
233        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
234      else if (calend == 'earth_366d') then
235        call ioconf_calendar('gregorian')
236        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
237      else if (calend == 'titan') then
238!        call ioconf_calendar('titan')
239        write(lunout,*)'CALENDRIER CHOISI: Titan'
240        abort_message = 'A FAIRE...'
241        call abort_gcm(modname,abort_message,1)
242      else if (calend == 'venus') then
243!        call ioconf_calendar('venus')
244        write(lunout,*)'CALENDRIER CHOISI: Venus'
245        abort_message = 'A FAIRE...'
246        call abort_gcm(modname,abort_message,1)
247      else
248        abort_message = 'Mauvais choix de calendrier'
249        call abort_gcm(modname,abort_message,1)
250      endif
251#endif
252c-----------------------------------------------------------------------
253
254      IF (config_inca /= 'none') THEN
255#ifdef INCA
256         call init_const_lmdz(
257     $        nbtr,anneeref,dayref,
258     $        iphysiq,day_step,nday,
259     $        nbsrf, is_oce,is_sic,
260     $        is_ter,is_lic)
261
262         call init_inca_para(
263     $        iim,jjm+1,llm,klon_glo,mpi_size,
264     $        distrib_phys,COMM_LMDZ)
265#endif
266      END IF
267
268c-----------------------------------------------------------------------
269c   Initialisation des traceurs
270c   ---------------------------
271c  Choix du nombre de traceurs et du schema pour l'advection
272c  dans fichier traceur.def, par default ou via INCA
273      call infotrac_init
274
275c Allocation de la tableau q : champs advectes   
276      ALLOCATE(q(ip1jmp1,llm,nqtot))
277
278c-----------------------------------------------------------------------
279c   Lecture de l'etat initial :
280c   ---------------------------
281
282c  lecture du fichier start.nc
283      if (read_start) then
284      ! we still need to run iniacademic to initialize some
285      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
286        if (iflag_phys.ne.1) then
287          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
288        endif
289
290        if (planet_type.eq."mars") then
291! POUR MARS, METTRE UNE FONCTION A PART, genre dynetat0_mars
292         abort_message = 'dynetat0_mars A FAIRE'
293         call abort_gcm(modname,abort_message,0)
294        else
295         CALL dynetat0("start.nc",vcov,ucov,
296     &              teta,q,masse,ps,phis, time_0)
297        endif ! of if (planet_type.eq."mars")
298       
299c       write(73,*) 'ucov',ucov
300c       write(74,*) 'vcov',vcov
301c       write(75,*) 'teta',teta
302c       write(76,*) 'ps',ps
303c       write(77,*) 'q',q
304
305      endif ! of if (read_start)
306
307c le cas echeant, creation d un etat initial
308      IF (prt_level > 9) WRITE(lunout,*)
309     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
310      if (.not.read_start) then
311         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
312      endif
313
314
315c-----------------------------------------------------------------------
316c   Lecture des parametres de controle pour la simulation :
317c   -------------------------------------------------------
318c  on recalcule eventuellement le pas de temps
319
320      IF(MOD(day_step,iperiod).NE.0) THEN
321        abort_message =
322     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
323        call abort_gcm(modname,abort_message,1)
324      ENDIF
325
326      IF(MOD(day_step,iphysiq).NE.0) THEN
327        abort_message =
328     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
329        call abort_gcm(modname,abort_message,1)
330      ENDIF
331
332      zdtvr    = daysec/REAL(day_step)
333        IF(dtvr.NE.zdtvr) THEN
334         WRITE(lunout,*)
335     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
336        ENDIF
337
338C
339C on remet le calendrier à zero si demande
340c
341      IF (raz_date == 1) THEN
342        annee_ref = anneeref
343        day_ref = dayref
344        day_ini = dayref
345        itau_dyn = 0
346        itau_phy = 0
347        time_0 = 0.
348        write(lunout,*)
349     .   'GCM: On reinitialise a la date lue dans gcm.def'
350      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
351        write(lunout,*)
352     .  'GCM: Attention les dates initiales lues dans le fichier'
353        write(lunout,*)
354     .  ' restart ne correspondent pas a celles lues dans '
355        write(lunout,*)' gcm.def'
356        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
357        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
358        write(lunout,*)' Pas de remise a zero'
359      ENDIF
360
361c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
362c        write(lunout,*)
363c     .  'GCM: Attention les dates initiales lues dans le fichier'
364c        write(lunout,*)
365c     .  ' restart ne correspondent pas a celles lues dans '
366c        write(lunout,*)' gcm.def'
367c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
368c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
369c        if (raz_date .ne. 1) then
370c          write(lunout,*)
371c     .    'GCM: On garde les dates du fichier restart'
372c        else
373c          annee_ref = anneeref
374c          day_ref = dayref
375c          day_ini = dayref
376c          itau_dyn = 0
377c          itau_phy = 0
378c          time_0 = 0.
379c          write(lunout,*)
380c     .   'GCM: On reinitialise a la date lue dans gcm.def'
381c        endif
382c      ELSE
383c        raz_date = 0
384c      endif
385
386#ifdef CPP_IOIPSL
387      mois = 1
388      heure = 0.
389! Ce n'est defini pour l'instant que pour la Terre...
390      if (planet_type.eq.'earth') then
391      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
392      jH_ref = jD_ref - int(jD_ref)
393      jD_ref = int(jD_ref)
394
395      call ioconf_startdate(INT(jD_ref), jH_ref)
396
397      write(lunout,*)'DEBUG'
398      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
399      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
400      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
401      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
402      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
403      else
404! A voir pour Titan et Venus
405        jD_ref=0
406        jH_ref=0
407      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
408      write(lunout,*)jD_ref,jH_ref
409      endif ! planet_type
410#else
411! Ehouarn: we still need to define JD_ref and JH_ref
412! and since we don't know how many days there are in a year
413! we set JD_ref to 0 (this should be improved ...)
414      jD_ref=0
415      jH_ref=0
416#endif
417
418c  nombre d'etats dans les fichiers demarrage et histoire
419      nbetatdem = nday / iecri
420      nbetatmoy = nday / periodav + 1
421
422      if (iflag_phys.eq.1) then
423      ! these initialisations have already been done (via iniacademic)
424      ! if running in SW or Newtonian mode
425c-----------------------------------------------------------------------
426c   Initialisation des constantes dynamiques :
427c   ------------------------------------------
428        dtvr = zdtvr
429        CALL iniconst
430
431c-----------------------------------------------------------------------
432c   Initialisation de la geometrie :
433c   --------------------------------
434        CALL inigeom
435
436c-----------------------------------------------------------------------
437c   Initialisation du filtre :
438c   --------------------------
439        CALL inifilr
440      endif ! of if (iflag_phys.eq.1)
441c
442c-----------------------------------------------------------------------
443c   Initialisation de la dissipation :
444c   ----------------------------------
445
446      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
447     *                tetagdiv, tetagrot , tetatemp              )
448
449c-----------------------------------------------------------------------
450c   Initialisation de la physique :
451c   -------------------------------
452
453      IF (call_iniphys.and.(iflag_phys.eq.1)) THEN
454         latfi(1)=rlatu(1)
455         lonfi(1)=0.
456         zcufi(1) = cu(1)
457         zcvfi(1) = cv(1)
458         DO j=2,jjm
459            DO i=1,iim
460               latfi((j-2)*iim+1+i)= rlatu(j)
461               lonfi((j-2)*iim+1+i)= rlonv(i)
462               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
463               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
464            ENDDO
465         ENDDO
466         latfi(ngridmx)= rlatu(jjp1)
467         lonfi(ngridmx)= 0.
468         zcufi(ngridmx) = cu(ip1jm+1)
469         zcvfi(ngridmx) = cv(ip1jm-iim)
470         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
471         WRITE(lunout,*)
472     .       'GCM: WARNING!!! vitesse verticale nulle dans la physique'
473
474! Initialisation de la physique: pose probleme quand on tourne
475! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
476! Il faut une cle CPP_PHYS
477#ifdef CPP_PHYS
478         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
479     ,                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp     )
480#endif ! CPP_PHYS
481         call_iniphys=.false.
482      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
483!#endif
484
485c-----------------------------------------------------------------------
486c   Initialisation des dimensions d'INCA :
487c   --------------------------------------
488      IF (config_inca /= 'none') THEN
489!$OMP PARALLEL
490#ifdef INCA
491         CALL init_inca_dim(klon_omp,llm,iim,jjm,
492     $        rlonu,rlatu,rlonv,rlatv)
493#endif
494!$OMP END PARALLEL
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 ET 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/,/logic/)
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.