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

Last change on this file since 2513 was 2507, checked in by romain.vande, 4 years ago

For LMDZ MARS: Update of day_ini, time and hour_ini in restart and retstartfi.
Hour_ini is obsolete. If we write one restart file: day_ini is the last day
of the simulation and the remaining time is in Time (often=0), if we write
multiple restart nothing changes

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