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

Last change on this file since 3093 was 2564, checked in by emillour, 3 years ago

Common dynamics:
Check that there is indeed a run.def file to get options, flags, etc. from.
Stop if there isn't, as most likely the defaults won't do.
EM

File size: 18.3 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      logical :: present
141
142c-----------------------------------------------------------------------
143c    variables pour l'initialisation de la physique :
144c    ------------------------------------------------
145!      INTEGER ngridmx
146!      PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm   )
147!      REAL zcufi(ngridmx),zcvfi(ngridmx)
148!      REAL latfi(ngridmx),lonfi(ngridmx)
149!      REAL airefi(ngridmx)
150!      SAVE latfi, lonfi, airefi
151     
152      INTEGER :: ierr
153
154
155c-----------------------------------------------------------------------
156c   Initialisations:
157c   ----------------
158
159      abort_message = 'last timestep reached'
160      modname = 'gcm'
161      lafin    = .FALSE.
162      dynhist_file = 'dyn_hist'
163      dynhistave_file = 'dyn_hist_ave'
164
165
166
167c----------------------------------------------------------------------
168c  Load flags and options from run.def
169c  ---------------------------------------
170c
171      ! Start by checking that there indeed is a run.def file
172      INQUIRE(file="run.def", EXIST=present)
173      IF(.not. present) then
174        CALL abort_gcm("gcm", "file run.def not found, Aborting!",1)
175      ENDIF
176
177      CALL conf_gcm( 99, .TRUE. )
178      if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm",
179     s "iphysiq must be a multiple of iperiod", 1)
180c
181c
182c------------------------------------
183c   Initialisation partie parallele
184c------------------------------------
185
186      CALL init_const_mpi
187      call init_parallel
188      call Read_Distrib
189
190!#ifdef CPP_PHYS
191!        CALL init_phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
192!#endif
193!      CALL set_bands
194!#ifdef CPP_PHYS
195!      CALL Init_interface_dyn_phys
196!#endif
197      CALL barrier
198
199      CALL set_bands
200      if (mpi_rank==0) call WriteBands
201      call SetDistrib(jj_Nb_Caldyn)
202
203c$OMP PARALLEL
204      call Init_Mod_hallo
205c$OMP END PARALLEL
206
207#ifdef CPP_PHYS
208!c$OMP PARALLEL
209!      call initcomgeomphy
210!c$OMP END PARALLEL
211#endif
212
213!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
214! Initialisation de XIOS
215!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
216
217
218c
219c Initialisations pour Cp(T) Venus
220      call ini_cpdet
221c
222c-----------------------------------------------------------------------
223c   Choix du calendrier
224c   -------------------
225
226c      calend = 'earth_365d'
227
228#ifdef CPP_IOIPSL
229      if (calend == 'earth_360d') then
230        call ioconf_calendar('360d')
231        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
232      else if (calend == 'earth_365d') then
233        call ioconf_calendar('noleap')
234        write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
235      else if (calend == 'earth_366d') then
236        call ioconf_calendar('gregorian')
237        write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
238      else if (calend == 'titan') then
239!        call ioconf_calendar('titan')
240        write(lunout,*)'CALENDRIER CHOISI: Titan'
241        abort_message = 'A FAIRE...'
242        call abort_gcm(modname,abort_message,1)
243      else if (calend == 'venus') then
244!        call ioconf_calendar('venus')
245        write(lunout,*)'CALENDRIER CHOISI: Venus'
246        abort_message = 'A FAIRE...'
247        call abort_gcm(modname,abort_message,1)
248      else
249        abort_message = 'Mauvais choix de calendrier'
250        call abort_gcm(modname,abort_message,1)
251      endif
252#endif
253c-----------------------------------------------------------------------
254c   Initialisation des traceurs
255c   ---------------------------
256c  Choix du nombre de traceurs et du schema pour l'advection
257c  dans fichier traceur.def, par default ou via INCA
258      call infotrac_init
259
260c Allocation de la tableau q : champs advectes   
261      ALLOCATE(q(ip1jmp1,llm,nqtot))
262
263c-----------------------------------------------------------------------
264c   Lecture de l'etat initial :
265c   ---------------------------
266
267c  lecture du fichier start.nc
268      if (read_start) then
269      ! we still need to run iniacademic to initialize some
270      ! constants & fields, if we run the 'newtonian' or 'SW' cases:
271        if (iflag_phys.ne.1) then
272          CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
273        endif
274
275        CALL dynetat0("start.nc",vcov,ucov,
276     &              teta,q,masse,ps,phis, time_0)
277       
278        ! Load relaxation fields (simple nudging). AS 09/2013
279        ! ---------------------------------------------------
280        if (planet_type.eq."generic") then
281         if (ok_guide) then
282           CALL relaxetat0("relax.nc")
283         endif
284        endif
285 
286c       write(73,*) 'ucov',ucov
287c       write(74,*) 'vcov',vcov
288c       write(75,*) 'teta',teta
289c       write(76,*) 'ps',ps
290c       write(77,*) 'q',q
291
292      endif ! of if (read_start)
293
294c le cas echeant, creation d un etat initial
295      IF (prt_level > 9) WRITE(lunout,*)
296     .              'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
297      if (.not.read_start) then
298         CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
299      endif
300
301
302c-----------------------------------------------------------------------
303c   Lecture des parametres de controle pour la simulation :
304c   -------------------------------------------------------
305c  on recalcule eventuellement le pas de temps
306
307      IF(MOD(day_step,iperiod).NE.0) THEN
308        abort_message =
309     .  'Il faut choisir un nb de pas par jour multiple de iperiod'
310        call abort_gcm(modname,abort_message,1)
311      ENDIF
312
313      IF(MOD(day_step,iphysiq).NE.0) THEN
314        abort_message =
315     * 'Il faut choisir un nb de pas par jour multiple de iphysiq'
316        call abort_gcm(modname,abort_message,1)
317      ENDIF
318
319      zdtvr    = daysec/REAL(day_step)
320        IF(dtvr.NE.zdtvr) THEN
321         WRITE(lunout,*)
322     .    'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
323        ENDIF
324      dtvr=zdtvr
325
326C
327C on remet le calendrier a zero si demande
328c
329      IF (start_time /= starttime) then
330        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
331     &,' fichier restart ne correspond pas a celle lue dans le run.def'
332        IF (raz_date == 1) then
333          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
334          start_time = starttime
335        ELSE
336          call abort_gcm("gcm", "'Je m''arrete'", 1)
337        ENDIF
338      ENDIF
339      IF (raz_date == 1) THEN
340        annee_ref = anneeref
341        day_ref = dayref
342        day_ini = dayref
343        itau_dyn = 0
344        itau_phy = 0
345        time_0 = 0.
346        write(lunout,*)
347     .   'GCM: On reinitialise a la date lue dans gcm.def'
348      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
349        write(lunout,*)
350     .  'GCM: Attention les dates initiales lues dans le fichier'
351        write(lunout,*)
352     .  ' restart ne correspondent pas a celles lues dans '
353        write(lunout,*)' gcm.def'
354        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
355        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
356        write(lunout,*)' Pas de remise a zero'
357      ENDIF
358
359c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
360c        write(lunout,*)
361c     .  'GCM: Attention les dates initiales lues dans le fichier'
362c        write(lunout,*)
363c     .  ' restart ne correspondent pas a celles lues dans '
364c        write(lunout,*)' gcm.def'
365c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
366c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
367c        if (raz_date .ne. 1) then
368c          write(lunout,*)
369c     .    'GCM: On garde les dates du fichier restart'
370c        else
371c          annee_ref = anneeref
372c          day_ref = dayref
373c          day_ini = dayref
374c          itau_dyn = 0
375c          itau_phy = 0
376c          time_0 = 0.
377c          write(lunout,*)
378c     .   'GCM: On reinitialise a la date lue dans gcm.def'
379c        endif
380c      ELSE
381c        raz_date = 0
382c      endif
383
384#ifdef CPP_IOIPSL
385      mois = 1
386      heure = 0.
387! Ce n'est defini pour l'instant que pour la Terre...
388      if (planet_type.eq.'earth') then
389      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
390      jH_ref = jD_ref - int(jD_ref)
391      jD_ref = int(jD_ref)
392
393      call ioconf_startdate(INT(jD_ref), jH_ref)
394
395      write(lunout,*)'DEBUG'
396      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
397      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
398      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
399      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
400      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
401!      else if (planet_type.eq.'titan') then
402!        jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
403      else
404! A voir pour Titan et Venus
405        jD_ref=0
406        jH_ref=0
407      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
408      write(lunout,*)jD_ref,jH_ref
409      endif ! planet_type
410#else
411! Ehouarn: we still need to define JD_ref and JH_ref
412! and since we don't know how many days there are in a year
413! we set JD_ref to 0 (this should be improved ...)
414      jD_ref=0
415      jH_ref=0
416#endif
417
418      if (iflag_phys.eq.1) then
419      ! these initialisations have already been done (via iniacademic)
420      ! if running in SW or Newtonian mode
421c-----------------------------------------------------------------------
422c   Initialisation des constantes dynamiques :
423c   ------------------------------------------
424        dtvr = zdtvr
425        CALL iniconst
426
427c-----------------------------------------------------------------------
428c   Initialisation de la geometrie :
429c   --------------------------------
430        CALL inigeom
431
432c-----------------------------------------------------------------------
433c   Initialisation du filtre :
434c   --------------------------
435        CALL inifilr
436      endif ! of if (iflag_phys.eq.1)
437c
438c-----------------------------------------------------------------------
439c   Initialisation de la dissipation :
440c   ----------------------------------
441
442      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
443     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
444
445c-----------------------------------------------------------------------
446c   Initialisation des I/O :
447c   ------------------------
448
449
450      if (nday>=0) then ! standard case
451        day_end=day_ini+nday
452      else ! special case when nday <0, run -nday dynamical steps
453        day_end=day_ini-nday/day_step
454      endif
455      if (less1day) then
456        day_end=day_ini+floor(time_0+fractday)
457      endif
458      if (ndynstep.gt.0) then
459        day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
460      endif
461     
462      WRITE(lunout,'(a,i7,a,i7)')
463     &             "run from day ",day_ini,"  to day",day_end
464
465#ifdef CPP_IOIPSL
466! Ce n'est defini pour l'instant que pour la Terre...
467      if (planet_type.eq.'earth') then
468      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
469      write (lunout,301)jour, mois, an
470      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
471      write (lunout,302)jour, mois, an
472      else
473! A voir pour Titan et Venus
474      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
475      endif ! planet_type
476
477 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
478 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
479#endif
480
481
482c-----------------------------------------------------------------------
483c   Initialisation de la physique :
484c   -------------------------------
485
486      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
487! Physics
488#ifdef CPP_PHYS
489         CALL iniphysiq(iim,jjm,llm,
490     &                distrib_phys(mpi_rank),comm_lmdz,
491     &                daysec,day_ini,dtphys/nsplit_phys,
492     &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp,
493     &                iflag_phys)
494#endif
495!         call_iniphys=.false.
496      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
497
498
499      if (planet_type=="mars") then
500         ! For Mars we transmit day_ini
501        !CALL dynredem0_p("restart.nc", day_ini, phis)
502        if (ecritstart.GT.0) then
503         CALL dynredem0_p("restart.nc", day_ini, phis)
504        else
505         CALL dynredem0_p("restart.nc", day_end, phis)
506        endif
507      else
508        CALL dynredem0_p("restart.nc", day_end, phis)
509      endif
510      ecripar = .TRUE.
511
512#ifdef CPP_IOIPSL
513      time_step = zdtvr
514      IF (mpi_rank==0) then
515        if (ok_dyn_ins) then
516        ! initialize output file for instantaneous outputs
517        ! t_ops = iecri * daysec ! do operations every t_ops
518        t_ops =((1.0*iecri)/day_step) * daysec 
519        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
520        CALL inithist(day_ref,annee_ref,time_step,
521     &                  t_ops,t_wrt)
522        endif
523
524        IF (ok_dyn_ave) THEN
525          ! initialize output file for averaged outputs
526          t_ops = iperiod * time_step ! do operations every t_ops
527          t_wrt = periodav * daysec   ! write output every t_wrt
528          CALL initdynav(day_ref,annee_ref,time_step,
529     &                   t_ops,t_wrt)
530!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
531!     .        t_ops, t_wrt, histaveid)
532        END IF
533      ENDIF
534      dtav = iperiod*dtvr/daysec
535#endif
536! #endif of #ifdef CPP_IOIPSL
537
538c  Choix des frequences de stokage pour le offline
539c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
540c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
541      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
542      istphy=istdyn/iphysiq     
543
544
545c
546c-----------------------------------------------------------------------
547c   Integration temporelle du modele :
548c   ----------------------------------
549
550c       write(78,*) 'ucov',ucov
551c       write(78,*) 'vcov',vcov
552c       write(78,*) 'teta',teta
553c       write(78,*) 'ps',ps
554c       write(78,*) 'q',q
555
556!c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
557!variable temps no longer exists
558c$OMP PARALLEL DEFAULT(SHARED)
559c       Copy all threadprivate variables from temps_mod
560c$OMP1 COPYIN(dt,jD_ref,jH_ref,start_time,hour_ini,day_ini,day_end)
561c$OMP1 COPYIN(annee_ref,day_ref,itau_dyn,itau_phy,itaufin,calend)
562c       Copy all threadprivate variables from logic_mod
563c$OMP1 COPYIN(purmats,forward,leapf,apphys,statcl,conser,apdiss,apdelq)
564c$OMP1 COPYIN(saison,ecripar,fxyhypb,ysinus,read_start,ok_guide)
565c$OMP1 COPYIN(ok_strato,tidal,ok_gradsfile,ok_limit,ok_etat0)
566c$OMP1 COPYIN(iflag_phys,iflag_trac)
567
568
569      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,
570     .              time_0)
571c$OMP END PARALLEL
572
573
574      END
575
Note: See TracBrowser for help on using the repository browser.