source: trunk/LMDZ.COMMON/libf/dyn3d/gcm.F @ 1076

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

Common dynamics: Small fix: nsplit_phys needs to be know in gcm if calling iniphysiq...
EM

File size: 18.7 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#else
11! if not using IOIPSL, we still need to use (a local version of) getin
12      USE ioipsl_getincom
13#endif
14
15
16#ifdef CPP_XIOS
17    ! ug Pour les sorties XIOS
18        USE wxios
19#endif
20
21      USE filtreg_mod
22      USE infotrac
23      USE control_mod, only: planet_type,nday,day_step,iperiod,iphysiq,
24     &                       raz_date,anneeref,starttime,dayref,
25     &                       ok_dyn_ins,ok_dyn_ave,iecri,periodav,
26     &                       less1day,fractday,ndynstep,nsplit_phys
27      use cpdet_mod, only: ini_cpdet
28
29#ifdef INCA
30! Only INCA needs these informations (from the Earth's physics)
31      USE indice_sol_mod
32#endif
33
34!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
35! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
36! A nettoyer. On ne veut qu'une ou deux routines d'interface
37! dynamique -> physique pour l'initialisation
38! Ehouarn: the following are needed with (parallel) physics:
39#ifdef CPP_PHYS
40      USE dimphy
41      USE comgeomphy
42#endif
43#ifdef INCA
44      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
45#endif
46!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
47
48      IMPLICIT NONE
49
50c      ......   Version  du 10/01/98    ..........
51
52c             avec  coordonnees  verticales hybrides
53c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
54
55c=======================================================================
56c
57c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
58c   -------
59c
60c   Objet:
61c   ------
62c
63c   GCM LMD nouvelle grille
64c
65c=======================================================================
66c
67c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
68c      et possibilite d'appeler une fonction f(y)  a derivee tangente
69c      hyperbolique a la  place de la fonction a derivee sinusoidale.
70c  ... Possibilite de choisir le schema pour l'advection de
71c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
72c
73c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
74c      Pour Van-Leer iadv=10
75c
76c-----------------------------------------------------------------------
77c   Declarations:
78c   -------------
79
80#include "dimensions.h"
81#include "paramet.h"
82#include "comconst.h"
83#include "comdissnew.h"
84#include "comvert.h"
85#include "comgeom.h"
86#include "logic.h"
87#include "temps.h"
88!!!!!!!!!!!#include "control.h"
89#include "ener.h"
90#include "description.h"
91#include "serre.h"
92!#include "com_io_dyn.h"
93#include "iniprint.h"
94#include "tracstoke.h"
95#ifdef INCA
96! Only INCA needs these informations (from the Earth's physics)
97!#include "indicesol.h"
98#endif
99
100
101      REAL zdtvr
102
103c   variables dynamiques
104      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
105      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
106      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
107      REAL ps(ip1jmp1)                       ! pression  au sol
108      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
109      REAL pks(ip1jmp1)                      ! exner au  sol
110      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
111      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
112      REAL masse(ip1jmp1,llm)                ! masse d'air
113      REAL phis(ip1jmp1)                     ! geopotentiel au sol
114      REAL phi(ip1jmp1,llm)                  ! geopotentiel
115      REAL w(ip1jmp1,llm)                    ! vitesse verticale
116
117c variables dynamiques intermediaire pour le transport
118
119c   variables pour le fichier histoire
120      REAL dtav      ! intervalle de temps elementaire
121
122      REAL time_0
123
124      LOGICAL lafin
125      INTEGER ij,iq,l,i,j
126
127
128      real time_step, t_wrt, t_ops
129
130      LOGICAL first
131
132      LOGICAL call_iniphys
133      data call_iniphys/.true./
134
135      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
136c+jld variables test conservation energie
137c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
138C     Tendance de la temp. potentiel d (theta)/ d t due a la
139C     tansformation d'energie cinetique en energie thermique
140C     cree par la dissipation
141      REAL dhecdt(ip1jmp1,llm)
142c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
143c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
144      CHARACTER (len=15) :: ztit
145c-jld
146
147
148      character (len=80) :: dynhist_file, dynhistave_file
149      character (len=20) :: modname
150      character (len=80) :: abort_message
151! locales pour gestion du temps
152      INTEGER :: an, mois, jour
153      REAL :: heure
154
155
156c-----------------------------------------------------------------------
157c    variables pour l'initialisation de la physique :
158c    ------------------------------------------------
159      INTEGER ngridmx
160      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
161      REAL zcufi(ngridmx),zcvfi(ngridmx)
162      REAL latfi(ngridmx),lonfi(ngridmx)
163      REAL airefi(ngridmx)
164      SAVE latfi, lonfi, airefi
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.nc'
175      dynhistave_file = 'dyn_hist_ave.nc'
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. )
186!#else
187!      CALL defrun( 99, .TRUE. , clesphy0 )
188!#endif
189
190!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
191! Initialisation de XIOS
192!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
193
194#ifdef CPP_XIOS
195        CALL wxios_init("LMDZ")
196#endif
197
198
199!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
200! FH 2008/05/02
201! A nettoyer. On ne veut qu'une ou deux routines d'interface
202! dynamique -> physique pour l'initialisation
203#ifdef CPP_PHYS
204      CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
205      call initcomgeomphy
206#endif
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(nbtr,anneeref,dayref,iphysiq,day_step,nday,
248     $        nbsrf, is_oce,is_sic,is_ter,is_lic)
249      call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
250#endif
251      END IF
252c
253c
254c------------------------------------
255c   Initialisation partie parallele
256c------------------------------------
257
258c
259c
260c-----------------------------------------------------------------------
261c   Initialisation des traceurs
262c   ---------------------------
263c  Choix du nombre de traceurs et du schema pour l'advection
264c  dans fichier traceur.def, par default ou via INCA
265      call infotrac_init
266
267c Allocation de la tableau q : champs advectes   
268      allocate(q(ip1jmp1,llm,nqtot))
269
270c-----------------------------------------------------------------------
271c   Lecture de l'etat initial :
272c   ---------------------------
273
274c  lecture du fichier start.nc
275      if (read_start) then
276      ! we still need to run iniacademic to initialize some
277      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
278        if (iflag_phys.ne.1) then
279          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
280        endif
281
282        if (planet_type.eq."mars") then
283! POUR MARS, METTRE UNE FONCTION A PART, genre dynetat0_mars
284         abort_message = 'dynetat0_mars A FAIRE'
285         call abort_gcm(modname,abort_message,0)
286        else
287         CALL dynetat0("start.nc",vcov,ucov,
288     &              teta,q,masse,ps,phis, time_0)
289        endif ! of if (planet_type.eq."mars")
290       
291        ! Load relaxation fields (simple nudging). AS 09/2013
292        ! ---------------------------------------------------
293        if (planet_type.eq."generic") then
294         if (ok_guide) then
295           CALL relaxetat0("relax.nc")
296         endif
297        endif
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
307      IF (type_trac == 'inca') THEN
308#ifdef INCA
309         call init_inca_dim(klon,llm,iim,jjm,
310     $        rlonu,rlatu,rlonv,rlatv)
311#endif
312      END IF
313
314
315c le cas echeant, creation d un etat initial
316      IF (prt_level > 9) WRITE(lunout,*)
317     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
318      if (.not.read_start) then
319         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
320      endif
321
322
323c-----------------------------------------------------------------------
324c   Lecture des parametres de controle pour la simulation :
325c   -------------------------------------------------------
326c  on recalcule eventuellement le pas de temps
327
328      IF(MOD(day_step,iperiod).NE.0) THEN
329        abort_message =
330     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
331        call abort_gcm(modname,abort_message,1)
332      ENDIF
333
334      IF(MOD(day_step,iphysiq).NE.0) THEN
335        abort_message =
336     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
337        call abort_gcm(modname,abort_message,1)
338      ENDIF
339
340      zdtvr    = daysec/REAL(day_step)
341        IF(dtvr.NE.zdtvr) THEN
342         WRITE(lunout,*)
343     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
344        ENDIF
345
346C
347C on remet le calendrier à zero si demande
348c
349      IF (start_time /= starttime) then
350        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
351     &,' fichier restart ne correspond pas à celle lue dans le run.def'
352        IF (raz_date == 1) then
353          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
354          start_time = starttime
355        ELSE
356          WRITE(lunout,*)'Je m''arrete'
357          CALL abort
358        ENDIF
359      ENDIF
360      IF (raz_date == 1) THEN
361        annee_ref = anneeref
362        day_ref = dayref
363        day_ini = dayref
364        itau_dyn = 0
365        itau_phy = 0
366        time_0 = 0.
367        write(lunout,*)
368     .   'GCM: On reinitialise a la date lue dans gcm.def'
369      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
370        write(lunout,*)
371     .  'GCM: Attention les dates initiales lues dans le fichier'
372        write(lunout,*)
373     .  ' restart ne correspondent pas a celles lues dans '
374        write(lunout,*)' gcm.def'
375        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
376        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
377        write(lunout,*)' Pas de remise a zero'
378      ENDIF
379
380c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
381c        write(lunout,*)
382c     .  'GCM: Attention les dates initiales lues dans le fichier'
383c        write(lunout,*)
384c     .  ' restart ne correspondent pas a celles lues dans '
385c        write(lunout,*)' gcm.def'
386c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
387c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
388c        if (raz_date .ne. 1) then
389c          write(lunout,*)
390c     .    'GCM: On garde les dates du fichier restart'
391c        else
392c          annee_ref = anneeref
393c          day_ref = dayref
394c          day_ini = dayref
395c          itau_dyn = 0
396c          itau_phy = 0
397c          time_0 = 0.
398c          write(lunout,*)
399c     .   'GCM: On reinitialise a la date lue dans gcm.def'
400c        endif
401c      ELSE
402c        raz_date = 0
403c      endif
404
405#ifdef CPP_IOIPSL
406      mois = 1
407      heure = 0.
408! Ce n'est defini pour l'instant que pour la Terre...
409      if (planet_type.eq.'earth') then
410      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
411      jH_ref = jD_ref - int(jD_ref)
412      jD_ref = int(jD_ref)
413
414      call ioconf_startdate(INT(jD_ref), jH_ref)
415
416      write(lunout,*)'DEBUG'
417      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
418      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
419      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
420      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
421      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
422      else
423! A voir pour Titan et Venus
424        jD_ref=0
425        jH_ref=0
426      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
427      write(lunout,*)jD_ref,jH_ref
428      endif ! planet_type
429#else
430! Ehouarn: we still need to define JD_ref and JH_ref
431! and since we don't know how many days there are in a year
432! we set JD_ref to 0 (this should be improved ...)
433      jD_ref=0
434      jH_ref=0
435#endif
436
437      if (iflag_phys.eq.1) then
438      ! these initialisations have already been done (via iniacademic)
439      ! if running in SW or Newtonian mode
440c-----------------------------------------------------------------------
441c   Initialisation des constantes dynamiques :
442c   ------------------------------------------
443        dtvr = zdtvr
444        CALL iniconst
445
446c-----------------------------------------------------------------------
447c   Initialisation de la geometrie :
448c   --------------------------------
449        CALL inigeom
450
451c-----------------------------------------------------------------------
452c   Initialisation du filtre :
453c   --------------------------
454        CALL inifilr
455      endif ! of if (iflag_phys.eq.1)
456c
457c-----------------------------------------------------------------------
458c   Initialisation de la dissipation :
459c   ----------------------------------
460
461      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
462     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
463
464c-----------------------------------------------------------------------
465c   Initialisation de la physique :
466c   -------------------------------
467
468      IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
469         latfi(1)=rlatu(1)
470         lonfi(1)=0.
471         zcufi(1) = cu(1)
472         zcvfi(1) = cv(1)
473         DO j=2,jjm
474            DO i=1,iim
475               latfi((j-2)*iim+1+i)= rlatu(j)
476               lonfi((j-2)*iim+1+i)= rlonv(i)
477               zcufi((j-2)*iim+1+i) = cu((j-1)*iip1+i)
478               zcvfi((j-2)*iim+1+i) = cv((j-1)*iip1+i)
479            ENDDO
480         ENDDO
481         latfi(ngridmx)= rlatu(jjp1)
482         lonfi(ngridmx)= 0.
483         zcufi(ngridmx) = cu(ip1jm+1)
484         zcvfi(ngridmx) = cv(ip1jm-iim)
485
486         ! build airefi(), mesh area on physics grid
487         CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
488         ! Poles are single points on physics grid
489         airefi(1)=airefi(1)*iim
490         airefi(ngridmx)=airefi(ngridmx)*iim
491
492! Initialisation de la physique: pose probleme quand on tourne
493! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
494! Il faut une cle CPP_PHYS
495#ifdef CPP_PHYS
496         CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
497     &                latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
498     &                iflag_phys)
499#endif
500         call_iniphys=.false.
501      ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
502
503c  numero de stockage pour les fichiers de redemarrage:
504
505c-----------------------------------------------------------------------
506c   Initialisation des I/O :
507c   ------------------------
508
509      day_end=day_ini+nday
510      if (less1day) then
511        day_end=day_ini+floor(time_0+fractday)
512      endif
513      if (ndynstep.gt.0) then
514        day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
515      endif
516     
517      WRITE(lunout,'(a,i7,a,i7)')
518     &             "run from day ",day_ini,"  to day",day_end
519
520#ifdef CPP_IOIPSL
521! Ce n'est defini pour l'instant que pour la Terre...
522      if (planet_type.eq.'earth') then
523      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
524      write (lunout,301)jour, mois, an
525      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
526      write (lunout,302)jour, mois, an
527      else
528! A voir pour Titan et Venus
529      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
530      endif ! planet_type
531
532 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
533 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
534#endif
535
536      if (planet_type.eq."mars") then
537! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
538         abort_message = 'dynredem0_mars A FAIRE'
539         call abort_gcm(modname,abort_message,0)
540      else
541        CALL dynredem0("restart.nc", day_end, phis)
542      endif ! of if (planet_type.eq."mars")
543
544      ecripar = .TRUE.
545
546#ifdef CPP_IOIPSL
547      time_step = zdtvr
548      if (ok_dyn_ins) then
549        ! initialize output file for instantaneous outputs
550        ! t_ops = iecri * daysec ! do operations every t_ops
551        t_ops =((1.0*iecri)/day_step) * daysec 
552        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
553        CALL inithist(day_ref,annee_ref,time_step,
554     &              t_ops,t_wrt)
555      endif
556
557      IF (ok_dyn_ave) THEN
558        ! initialize output file for averaged outputs
559        t_ops = iperiod * time_step ! do operations every t_ops
560        t_wrt = periodav * daysec   ! write output every t_wrt
561        CALL initdynav(day_ref,annee_ref,time_step,
562     &       t_ops,t_wrt)
563      END IF
564      dtav = iperiod*dtvr/daysec
565#endif
566! #endif of #ifdef CPP_IOIPSL
567
568c  Choix des frequences de stokage pour le offline
569c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
570c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
571      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
572      istphy=istdyn/iphysiq     
573
574
575c
576c-----------------------------------------------------------------------
577c   Integration temporelle du modele :
578c   ----------------------------------
579
580c       write(78,*) 'ucov',ucov
581c       write(78,*) 'vcov',vcov
582c       write(78,*) 'teta',teta
583c       write(78,*) 'ps',ps
584c       write(78,*) 'q',q
585
586
587      CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,
588     .              time_0)
589
590      END
591
Note: See TracBrowser for help on using the repository browser.