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

Last change on this file since 3604 was 3604, checked in by jbclement, 9 days ago

Mars PCM:
Bug correction introduced in r3574: the program has to be in a fixed-form Fortran format.
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 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!!! 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.