source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/dyn3d/gcm.F @ 222

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

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