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

Last change on this file since 2436 was 2272, checked in by emillour, 5 years ago

Common dynamics:
Remove unnecessary call to ini_getparam() in gcm.F, which moreover would create an (empty) out.def file at runtime.
EM

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