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

Last change on this file since 1428 was 1422, checked in by milmd, 10 years ago

In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).

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