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

Last change on this file since 3594 was 3574, checked in by jbclement, 10 days ago

COMMON:
The compilation generates a Fortran subroutine to track compilation and version (SVN or Git) details through the executable file. Put the argument 'version' as an option when executing the code to display these information and create a file "version_diff.txt" containing the diff result.
It can work with every program but it has been implemented only for: 'gcm', 'parallel gcm', 'pem', '1D Mars PCM', 'Mars newstart', 'Mars start2archive' and '1D Generic PCM'.
JBC

File size: 18.8 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 '
169     &//'is 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!!! changement de pas de temps',dtvr,'>',zdtvr
337        ENDIF
338      dtvr=zdtvr
339
340C
341C on remet le calendrier a zero si demande
342c
343      IF (start_time /= starttime) then
344        WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
345     &,' fichier restart ne correspond pas a celle lue dans le run.def'
346        IF (raz_date == 1) then
347          WRITE(lunout,*)'Je prends l''heure lue dans run.def'
348          start_time = starttime
349        ELSE
350          call abort_gcm("gcm", "'Je m''arrete'", 1)
351        ENDIF
352      ENDIF
353      IF (raz_date == 1) THEN
354        annee_ref = anneeref
355        day_ref = dayref
356        day_ini = dayref
357        itau_dyn = 0
358        itau_phy = 0
359        time_0 = 0.
360        write(lunout,*)
361     .   'GCM: On reinitialise a la date lue dans gcm.def'
362      ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
363        write(lunout,*)
364     .  'GCM: Attention les dates initiales lues dans le fichier'
365        write(lunout,*)
366     .  ' restart ne correspondent pas a celles lues dans '
367        write(lunout,*)' gcm.def'
368        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
369        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
370        write(lunout,*)' Pas de remise a zero'
371      ENDIF
372
373c      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
374c        write(lunout,*)
375c     .  'GCM: Attention les dates initiales lues dans le fichier'
376c        write(lunout,*)
377c     .  ' restart ne correspondent pas a celles lues dans '
378c        write(lunout,*)' gcm.def'
379c        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
380c        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
381c        if (raz_date .ne. 1) then
382c          write(lunout,*)
383c     .    'GCM: On garde les dates du fichier restart'
384c        else
385c          annee_ref = anneeref
386c          day_ref = dayref
387c          day_ini = dayref
388c          itau_dyn = 0
389c          itau_phy = 0
390c          time_0 = 0.
391c          write(lunout,*)
392c     .   'GCM: On reinitialise a la date lue dans gcm.def'
393c        endif
394c      ELSE
395c        raz_date = 0
396c      endif
397
398#ifdef CPP_IOIPSL
399      mois = 1
400      heure = 0.
401! Ce n'est defini pour l'instant que pour la Terre...
402      if (planet_type.eq.'earth') then
403      call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
404      jH_ref = jD_ref - int(jD_ref)
405      jD_ref = int(jD_ref)
406
407      call ioconf_startdate(INT(jD_ref), jH_ref)
408
409      write(lunout,*)'DEBUG'
410      write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
411      write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
412      call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
413      write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
414      write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
415!      else if (planet_type.eq.'titan') then
416!        jD_ref=1 ! Only if we use old Titan starts (J.V.O 2016)
417      else
418! A voir pour Titan et Venus
419        jD_ref=0
420        jH_ref=0
421      write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
422      write(lunout,*)jD_ref,jH_ref
423      endif ! planet_type
424#else
425! Ehouarn: we still need to define JD_ref and JH_ref
426! and since we don't know how many days there are in a year
427! we set JD_ref to 0 (this should be improved ...)
428      jD_ref=0
429      jH_ref=0
430#endif
431
432      if (iflag_phys.eq.1) then
433      ! these initialisations have already been done (via iniacademic)
434      ! if running in SW or Newtonian mode
435c-----------------------------------------------------------------------
436c   Initialisation des constantes dynamiques :
437c   ------------------------------------------
438        dtvr = zdtvr
439        CALL iniconst
440
441c-----------------------------------------------------------------------
442c   Initialisation de la geometrie :
443c   --------------------------------
444        CALL inigeom
445
446c-----------------------------------------------------------------------
447c   Initialisation du filtre :
448c   --------------------------
449        CALL inifilr
450      endif ! of if (iflag_phys.eq.1)
451c
452c-----------------------------------------------------------------------
453c   Initialisation de la dissipation :
454c   ----------------------------------
455
456      CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   ,
457     *                tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
458
459c-----------------------------------------------------------------------
460c   Initialisation des I/O :
461c   ------------------------
462
463
464      if (nday>=0) then ! standard case
465        day_end=day_ini+nday
466      else ! special case when nday <0, run -nday dynamical steps
467        day_end=day_ini-nday/day_step
468      endif
469      if (less1day) then
470        day_end=day_ini+floor(time_0+fractday)
471      endif
472      if (ndynstep.gt.0) then
473        day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
474      endif
475     
476      WRITE(lunout,'(a,i7,a,i7)')
477     &             "run from day ",day_ini,"  to day",day_end
478
479#ifdef CPP_IOIPSL
480! Ce n'est defini pour l'instant que pour la Terre...
481      if (planet_type.eq.'earth') then
482      call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
483      write (lunout,301)jour, mois, an
484      call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
485      write (lunout,302)jour, mois, an
486      else
487! A voir pour Titan et Venus
488      write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
489      endif ! planet_type
490
491 301  FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
492 302  FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
493#endif
494
495
496c-----------------------------------------------------------------------
497c   Initialisation de la physique :
498c   -------------------------------
499
500      IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
501! Physics
502#ifdef CPP_PHYS
503         CALL iniphysiq(iim,jjm,llm,
504     &                distrib_phys(mpi_rank),comm_lmdz,
505     &                daysec,day_ini,dtphys/nsplit_phys,
506     &                rlatu,rlatv,rlonu,rlonv,aire,cu,cv,
507     &                rad,g,r,cpp,iflag_phys)
508#endif
509!         call_iniphys=.false.
510      ENDIF ! of IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100))
511
512
513      if (planet_type=="mars") then
514         ! For Mars we transmit day_ini
515        !CALL dynredem0_p("restart.nc", day_ini, phis)
516        if (ecritstart.GT.0) then
517         CALL dynredem0_p("restart.nc", day_ini, phis)
518        else
519         CALL dynredem0_p("restart.nc", day_end, phis)
520        endif
521      else
522        CALL dynredem0_p("restart.nc", day_end, phis)
523      endif
524      ecripar = .TRUE.
525
526#ifdef CPP_IOIPSL
527      time_step = zdtvr
528      IF (mpi_rank==0) then
529        if (ok_dyn_ins) then
530        ! initialize output file for instantaneous outputs
531        ! t_ops = iecri * daysec ! do operations every t_ops
532        t_ops =((1.0*iecri)/day_step) * daysec 
533        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
534        CALL inithist(day_ref,annee_ref,time_step,
535     &                  t_ops,t_wrt)
536        endif
537
538        IF (ok_dyn_ave) THEN
539          ! initialize output file for averaged outputs
540          t_ops = iperiod * time_step ! do operations every t_ops
541          t_wrt = periodav * daysec   ! write output every t_wrt
542          CALL initdynav(day_ref,annee_ref,time_step,
543     &                   t_ops,t_wrt)
544!         CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
545!     .        t_ops, t_wrt, histaveid)
546        END IF
547      ENDIF
548      dtav = iperiod*dtvr/daysec
549#endif
550! #endif of #ifdef CPP_IOIPSL
551
552c  Choix des frequences de stokage pour le offline
553c      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
554c      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
555      istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
556      istphy=istdyn/iphysiq     
557
558
559c
560c-----------------------------------------------------------------------
561c   Integration temporelle du modele :
562c   ----------------------------------
563
564c       write(78,*) 'ucov',ucov
565c       write(78,*) 'vcov',vcov
566c       write(78,*) 'teta',teta
567c       write(78,*) 'ps',ps
568c       write(78,*) 'q',q
569
570!c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
571!variable temps no longer exists
572c$OMP PARALLEL DEFAULT(SHARED)
573c       Copy all threadprivate variables from temps_mod
574c$OMP1 COPYIN(dt,jD_ref,jH_ref,start_time,hour_ini,day_ini,day_end)
575c$OMP1 COPYIN(annee_ref,day_ref,itau_dyn,itau_phy,itaufin,calend)
576c       Copy all threadprivate variables from logic_mod
577c$OMP1 COPYIN(purmats,forward,leapf,apphys,statcl,conser,apdiss,apdelq)
578c$OMP1 COPYIN(saison,ecripar,fxyhypb,ysinus,read_start,ok_guide)
579c$OMP1 COPYIN(ok_strato,tidal,ok_gradsfile,ok_limit,ok_etat0)
580c$OMP1 COPYIN(iflag_phys,iflag_trac)
581
582
583      CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,
584     .              time_0)
585c$OMP END PARALLEL
586
587
588      END
589
Note: See TracBrowser for help on using the repository browser.