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

Last change on this file since 3615 was 3615, checked in by emillour, 29 hours ago

Venus PCM: Corrections to enable 1+1=2

  • store correctly the time_of_day in restart.nc to enable proper restart
  • enforce recomputation of CP in the physics at all time steps (otherwise when without thermosphere the value was only computed at first step and kept unchanged afterwards).

EM

File size: 18.5 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      USE mod_const_mpi, ONLY: init_const_mpi
14      USE parallel_lmdz
15      USE infotrac
16!#ifdef CPP_PHYS
17!      USE mod_interface_dyn_phys
18!#endif
19      USE mod_hallo
20      USE Bands
21      USE getparam
22      USE filtreg_mod
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     &                       ecritstart
28      use cpdet_mod, only: ini_cpdet
29      use version_info_mod, only: print_version_info
30
31
32! Ehouarn: the following are needed with (parallel) physics:
33#ifdef CPP_PHYS
34      USE iniphysiq_mod, ONLY: iniphysiq
35      USE mod_grid_phy_lmdz
36!      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
37      USE mod_phys_lmdz_omp_data, ONLY: klon_omp
38      USE dimphy
39#endif
40      USE comconst_mod, ONLY: daysec,dtvr,dtphys,rad,g,r,cpp
41      USE logic_mod
42      USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref,
43     .          itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end,
44     .          dt,hour_ini,itaufin
45      IMPLICIT NONE
46
47c      ......   Version  du 10/01/98    ..........
48
49c             avec  coordonnees  verticales hybrides
50c   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
51
52c=======================================================================
53c
54c   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
55c   -------
56c
57c   Objet:
58c   ------
59c
60c   GCM LMD nouvelle grille
61c
62c=======================================================================
63c
64c  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
65c      et possibilite d'appeler une fonction f(y)  a derivee tangente
66c      hyperbolique a la  place de la fonction a derivee sinusoidale.
67c  ... Possibilite de choisir le schema pour l'advection de
68c        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
69c
70c      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
71c      Pour Van-Leer iadv=10
72c
73c-----------------------------------------------------------------------
74c   Declarations:
75c   -------------
76
77#include "dimensions.h"
78#include "paramet.h"
79#include "comdissnew.h"
80#include "comgeom.h"
81!!!!!!!!!!!#include "control.h"
82!#include "com_io_dyn.h"
83#include "iniprint.h"
84#include "tracstoke.h"
85
86
87      REAL zdtvr
88
89c   variables dynamiques
90      REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
91      REAL teta(ip1jmp1,llm)                 ! temperature potentielle
92      REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
93      REAL ps(ip1jmp1)                       ! pression  au sol
94c      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
95c      REAL pks(ip1jmp1)                      ! exner au  sol
96c      REAL pk(ip1jmp1,llm)                   ! exner au milieu des couches
97c      REAL pkf(ip1jmp1,llm)                  ! exner filt.au milieu des couches
98      REAL masse(ip1jmp1,llm)                ! masse d'air
99      REAL phis(ip1jmp1)                     ! geopotentiel au sol
100c      REAL phi(ip1jmp1,llm)                  ! geopotentiel
101c      REAL w(ip1jmp1,llm)                    ! vitesse verticale
102
103c variables dynamiques intermediaire pour le transport
104
105c   variables pour le fichier histoire
106      REAL dtav      ! intervalle de temps elementaire
107
108      REAL time_0
109
110      LOGICAL lafin
111c      INTEGER ij,iq,l,i,j
112      INTEGER i,j
113
114
115      real time_step, t_wrt, t_ops
116
117
118!      LOGICAL call_iniphys
119!      data call_iniphys/.true./
120
121c      REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
122c+jld variables test conservation energie
123c      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
124C     Tendance de la temp. potentiel d (theta)/ d t due a la
125C     tansformation d'energie cinetique en energie thermique
126C     cree par la dissipation
127c      REAL dhecdt(ip1jmp1,llm)
128c      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
129c      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
130c      CHARACTER (len=15) :: ztit
131c-jld
132
133
134      character (len=80) :: dynhist_file, dynhistave_file
135      character (len=20) :: modname
136      character (len=80) :: abort_message
137      character(100)     :: arg ! To read command-line arguments
138! locales pour gestion du temps
139      INTEGER :: an, mois, jour
140      REAL :: heure
141
142      logical :: present
143
144c-----------------------------------------------------------------------
145c    variables pour l'initialisation de la physique :
146c    ------------------------------------------------
147!      INTEGER ngridmx
148!      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
149!      REAL zcufi(ngridmx),zcvfi(ngridmx)
150!      REAL latfi(ngridmx),lonfi(ngridmx)
151!      REAL airefi(ngridmx)
152!      SAVE latfi, lonfi, airefi
153     
154      INTEGER :: ierr
155
156
157c-----------------------------------------------------------------------
158c   Initialisations:
159c   ----------------
160
161      if (command_argument_count() > 0) then ! Get the number of command-line arguments
162          call get_command_argument(1,arg) ! Read the argument given to the program
163          select case (trim(adjustl(arg)))
164              case('version')
165                  call print_version_info()
166                  stop
167              case default
168                  error stop 'The argument given to the program is '
169     &//'unknown!'
170          end select
171      endif
172
173      abort_message = 'last timestep reached'
174      modname = 'gcm'
175      lafin    = .FALSE.
176      dynhist_file = 'dyn_hist'
177      dynhistave_file = 'dyn_hist_ave'
178
179
180
181c----------------------------------------------------------------------
182c  Load flags and options from run.def
183c  ---------------------------------------
184c
185      ! Start by checking that there indeed is a run.def file
186      INQUIRE(file="run.def", EXIST=present)
187      IF(.not. present) then
188        CALL abort_gcm("gcm", "file run.def not found, Aborting!",1)
189      ENDIF
190
191      CALL conf_gcm( 99, .TRUE. )
192      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
193     s "iphysiq must be a multiple of iperiod", 1)
194c
195c
196c------------------------------------
197c   Initialisation partie parallele
198c------------------------------------
199
200      CALL init_const_mpi
201      call init_parallel
202      call Read_Distrib
203
204!#ifdef CPP_PHYS
205!        CALL init_phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
206!#endif
207!      CALL set_bands
208!#ifdef CPP_PHYS
209!      CALL Init_interface_dyn_phys
210!#endif
211      CALL barrier
212
213      CALL set_bands
214      if (mpi_rank==0) call WriteBands
215      call SetDistrib(jj_Nb_Caldyn)
216
217c$OMP PARALLEL
218      call Init_Mod_hallo
219c$OMP END PARALLEL
220
221#ifdef CPP_PHYS
222!c$OMP PARALLEL
223!      call initcomgeomphy
224!c$OMP END PARALLEL
225#endif
226
227!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
228! Initialisation de XIOS
229!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
230
231
232c
233c Initialisations pour Cp(T) Venus
234      call ini_cpdet
235c
236c-----------------------------------------------------------------------
237c   Choix du calendrier
238c   -------------------
239
240c      calend = 'earth_365d'
241
242#ifdef CPP_IOIPSL
243      if (calend == 'earth_360d') then
244        call ioconf_calendar('360d')
245        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
246      else if (calend == 'earth_365d') then
247        call ioconf_calendar('noleap')
248        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
249      else if (calend == 'earth_366d') then
250        call ioconf_calendar('gregorian')
251        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
252      else if (calend == 'titan') then
253!        call ioconf_calendar('titan')
254        write(lunout,*)'CALENDRIER CHOISI: Titan'
255        abort_message = 'A FAIRE...'
256        call abort_gcm(modname,abort_message,1)
257      else if (calend == 'venus') then
258!        call ioconf_calendar('venus')
259        write(lunout,*)'CALENDRIER CHOISI: Venus'
260        abort_message = 'A FAIRE...'
261        call abort_gcm(modname,abort_message,1)
262      else
263        abort_message = 'Mauvais choix de calendrier'
264        call abort_gcm(modname,abort_message,1)
265      endif
266#endif
267c-----------------------------------------------------------------------
268c   Initialisation des traceurs
269c   ---------------------------
270c  Choix du nombre de traceurs et du schema pour l'advection
271c  dans fichier traceur.def, par default ou via INCA
272      call infotrac_init
273
274c Allocation de la tableau q : champs advectes   
275      ALLOCATE(q(ip1jmp1,llm,nqtot))
276
277c-----------------------------------------------------------------------
278c   Lecture de l'etat initial :
279c   ---------------------------
280
281c  lecture du fichier start.nc
282      if (read_start) then
283      ! we still need to run iniacademic to initialize some
284      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
285        if (iflag_phys.ne.1) then
286          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
287        endif
288
289        CALL dynetat0("start.nc",vcov,ucov,
290     &              teta,q,masse,ps,phis, time_0)
291       
292        ! Load relaxation fields (simple nudging). AS 09/2013
293        ! ---------------------------------------------------
294        if (planet_type.eq."generic") then
295         if (ok_guide) then
296           CALL relaxetat0("relax.nc")
297         endif
298        endif
299 
300c       write(73,*) 'ucov',ucov
301c       write(74,*) 'vcov',vcov
302c       write(75,*) 'teta',teta
303c       write(76,*) 'ps',ps
304c       write(77,*) 'q',q
305
306      endif ! of if (read_start)
307
308c le cas echeant, creation d un etat initial
309      IF (prt_level > 9) WRITE(lunout,*)
310     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
311      if (.not.read_start) then
312         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
313      endif
314
315
316c-----------------------------------------------------------------------
317c   Lecture des parametres de controle pour la simulation :
318c   -------------------------------------------------------
319c  on recalcule eventuellement le pas de temps
320
321      IF(MOD(day_step,iperiod).NE.0) THEN
322        abort_message =
323     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
324        call abort_gcm(modname,abort_message,1)
325      ENDIF
326
327      IF(MOD(day_step,iphysiq).NE.0) THEN
328        abort_message =
329     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
330        call abort_gcm(modname,abort_message,1)
331      ENDIF
332
333      zdtvr    = daysec/REAL(day_step)
334        IF(dtvr.NE.zdtvr) THEN
335         WRITE(lunout,*)
336     .    'WARNING!!! time step change: ',dtvr,'>',zdtvr
337        ENDIF
338      dtvr=zdtvr
339
340C
341C reset calendar if requested
342c
343      IF (raz_date == 1) then
344        WRITE(lunout,*)'Reinitializing to start time from run.def'
345        start_time = starttime
346      ENDIF
347
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 if (planet_type.eq.'titan') then
411!        jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
412      else
413! A voir pour Titan et Venus
414        jD_ref=0
415        jH_ref=0
416      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
417      write(lunout,*)jD_ref,jH_ref
418      endif ! planet_type
419#else
420! Ehouarn: we still need to define JD_ref and JH_ref
421! and since we don't know how many days there are in a year
422! we set JD_ref to 0 (this should be improved ...)
423      jD_ref=0
424      jH_ref=0
425#endif
426
427      if (iflag_phys.eq.1) then
428      ! these initialisations have already been done (via iniacademic)
429      ! if running in SW or Newtonian mode
430c-----------------------------------------------------------------------
431c   Initialisation des constantes dynamiques :
432c   ------------------------------------------
433        dtvr = zdtvr
434        CALL iniconst
435
436c-----------------------------------------------------------------------
437c   Initialisation de la geometrie :
438c   --------------------------------
439        CALL inigeom
440
441c-----------------------------------------------------------------------
442c   Initialisation du filtre :
443c   --------------------------
444        CALL inifilr
445      endif ! of if (iflag_phys.eq.1)
446c
447c-----------------------------------------------------------------------
448c   Initialisation de la dissipation :
449c   ----------------------------------
450
451      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
452     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
453
454c-----------------------------------------------------------------------
455c   Initialisation des I/O :
456c   ------------------------
457
458
459      if (nday>=0) then ! standard case
460        day_end=day_ini+nday
461      else ! special case when nday <0, run -nday dynamical steps
462        day_end=day_ini-nday/day_step
463      endif
464      if (less1day) then
465        day_end=day_ini+floor(time_0+fractday)
466      endif
467      if (ndynstep.gt.0) then
468        day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
469      endif
470     
471      WRITE(lunout,'(a,i7,a,i7)')
472     &             "run from day ",day_ini,"  to day",day_end
473
474#ifdef CPP_IOIPSL
475! Ce n'est defini pour l'instant que pour la Terre...
476      if (planet_type.eq.'earth') then
477      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
478      write (lunout,301)jour, mois, an
479      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
480      write (lunout,302)jour, mois, an
481      else
482! A voir pour Titan et Venus
483      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
484      endif ! planet_type
485
486 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
487 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
488#endif
489
490
491c-----------------------------------------------------------------------
492c   Initialisation de la physique :
493c   -------------------------------
494
495      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
496! Physics
497#ifdef CPP_PHYS
498         CALL iniphysiq(iim,jjm,llm,
499     &                distrib_phys(mpi_rank),comm_lmdz,
500     &                daysec,day_ini,dtphys/nsplit_phys,
501     &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,
502     &                rad,g,r,cpp,iflag_phys)
503#endif
504!         call_iniphys=.false.
505      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
506
507
508      if (planet_type=="mars") then
509         ! For Mars we transmit day_ini
510        !CALL dynredem0_p("restart.nc", day_ini, phis)
511        if (ecritstart.GT.0) then
512         CALL dynredem0_p("restart.nc", day_ini, phis)
513        else
514         CALL dynredem0_p("restart.nc", day_end, phis)
515        endif
516      else
517        CALL dynredem0_p("restart.nc", day_end, phis)
518      endif
519      ecripar = .TRUE.
520
521#ifdef CPP_IOIPSL
522      time_step = zdtvr
523      IF (mpi_rank==0) then
524        if (ok_dyn_ins) then
525        ! initialize output file for instantaneous outputs
526        ! t_ops = iecri * daysec ! do operations every t_ops
527        t_ops =((1.0*iecri)/day_step) * daysec 
528        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
529        CALL inithist(day_ref,annee_ref,time_step,
530     &                  t_ops,t_wrt)
531        endif
532
533        IF (ok_dyn_ave) THEN
534          ! initialize output file for averaged outputs
535          t_ops = iperiod * time_step ! do operations every t_ops
536          t_wrt = periodav * daysec   ! write output every t_wrt
537          CALL initdynav(day_ref,annee_ref,time_step,
538     &                   t_ops,t_wrt)
539!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
540!     .        t_ops, t_wrt, histaveid)
541        END IF
542      ENDIF
543      dtav = iperiod*dtvr/daysec
544#endif
545! #endif of #ifdef CPP_IOIPSL
546
547c  Choix des frequences de stokage pour le offline
548c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
549c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
550      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
551      istphy=istdyn/iphysiq     
552
553
554c
555c-----------------------------------------------------------------------
556c   Integration temporelle du modele :
557c   ----------------------------------
558
559c       write(78,*) 'ucov',ucov
560c       write(78,*) 'vcov',vcov
561c       write(78,*) 'teta',teta
562c       write(78,*) 'ps',ps
563c       write(78,*) 'q',q
564
565!c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
566!variable temps no longer exists
567c$OMP PARALLEL DEFAULT(SHARED)
568c       Copy all threadprivate variables from temps_mod
569c$OMP1 COPYIN(dt,jD_ref,jH_ref,start_time,hour_ini,day_ini,day_end)
570c$OMP1 COPYIN(annee_ref,day_ref,itau_dyn,itau_phy,itaufin,calend)
571c       Copy all threadprivate variables from logic_mod
572c$OMP1 COPYIN(purmats,forward,leapf,apphys,statcl,conser,apdiss,apdelq)
573c$OMP1 COPYIN(saison,ecripar,fxyhypb,ysinus,read_start,ok_guide)
574c$OMP1 COPYIN(ok_strato,tidal,ok_gradsfile,ok_limit,ok_etat0)
575c$OMP1 COPYIN(iflag_phys,iflag_trac)
576
577
578      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,
579     .              time_0)
580c$OMP END PARALLEL
581
582
583      END
584
Note: See TracBrowser for help on using the repository browser.