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

Last change on this file since 1243 was 1107, checked in by emillour, 11 years ago

Common dynamics: Updates and modifications to enable running Mars physics with

LMDZ.COMMON dynamics:

  • For compilation: adapted makelmdz, create_make_gcm and makelmdz_fcm, bld.cfg to compile aeronomy routines in "aerono$physique" if it exists, and added "-P -traditional" preprocessing flags in "arch-linux-ifort*"
  • Added function "cbrt.F" (cubic root) in 'bibio'
  • Adapted the reading/writing of dynamics (re)start.nc files for Mars. The main issue is that different information (on time, reference and current) is stored and used differently, hence a few if (planet_type =="mars") here and there. Moreover in the martian case there is the possibility to store fields over multiple times. Some Mars-specific variables (ecritphy,ecritstart,timestart) added in control_mod.F and (hour_ini) in temps.h

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