source: LMDZ5/branches/testing/libf/dyn3dmem/gcm.F90 @ 2343

Last change on this file since 2343 was 2298, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes -r2237:2291 into testing branch

  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 15.1 KB
RevLine 
[2225]1! $Id: $
[1632]2
[2263]3PROGRAM gcm
4
[1632]5#ifdef CPP_IOIPSL
[2263]6  USE IOIPSL
[1632]7#endif
8
[2263]9  USE mod_const_mpi, ONLY: init_const_mpi
10  USE parallel_lmdz
11  USE infotrac
[2239]12#ifdef CPP_PHYS
[2263]13  USE mod_interface_dyn_phys, ONLY: init_interface_dyn_phys
[2239]14#endif
[2263]15  USE mod_hallo
16  USE Bands
17  USE filtreg_mod
18  USE control_mod
[1632]19
[1785]20#ifdef INCA
[2263]21  ! Only INCA needs these informations (from the Earth's physics)
22  USE indice_sol_mod
23  USE mod_phys_lmdz_omp_data, ONLY: klon_omp
[1785]24#endif
25
[1673]26#ifdef CPP_PHYS
[2263]27  !      USE mod_grid_phy_lmdz
28  !      USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
29  !      USE dimphy
30  !      USE comgeomphy
[1632]31#endif
[2263]32  IMPLICIT NONE
[1632]33
[2263]34  !      ......   Version  du 10/01/98    ..........
[1632]35
[2263]36  !             avec  coordonnees  verticales hybrides
37  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
[1632]38
[2263]39  !=======================================================================
40  !
41  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
42  !   -------
43  !
44  !   Objet:
45  !   ------
46  !
47  !   GCM LMD nouvelle grille
48  !
49  !=======================================================================
50  !
51  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
52  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
53  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
54  !  ... Possibilite de choisir le schema pour l'advection de
55  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
56  !
57  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
58  !      Pour Van-Leer iadv=10
59  !
60  !-----------------------------------------------------------------------
61  !   Declarations:
62  !   -------------
63  include "dimensions.h"
64  include "paramet.h"
65  include "comconst.h"
66  include "comdissnew.h"
67  include "comvert.h"
68  include "comgeom.h"
69  include "logic.h"
70  include "temps.h"
71  include "ener.h"
72  include "description.h"
73  include "serre.h"
74  !include "com_io_dyn.h"
75  include "iniprint.h"
76  include "tracstoke.h"
[1658]77
[1657]78#ifdef INCA
[2263]79  ! Only INCA needs these informations (from the Earth's physics)
80  !include "indicesol.h"
[1657]81#endif
[1632]82
[2263]83  REAL zdtvr
[1632]84
[2263]85  !   variables dynamiques
86  REAL,ALLOCATABLE,SAVE  :: vcov(:,:),ucov(:,:) ! vents covariants
87  REAL,ALLOCATABLE,SAVE  :: teta(:,:)     ! temperature potentielle
88  REAL, ALLOCATABLE,SAVE :: q(:,:,:)      ! champs advectes
89  REAL,ALLOCATABLE,SAVE  :: ps(:)         ! pression  au sol
90  !      REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
91  REAL,ALLOCATABLE,SAVE  :: masse(:,:)    ! masse d'air
92  REAL,ALLOCATABLE,SAVE  :: phis(:)       ! geopotentiel au sol
93  !      REAL phi(ip1jmp1,llm)                  ! geopotentiel
94  !      REAL w(ip1jmp1,llm)                    ! vitesse verticale
[1632]95
[2263]96  ! variables dynamiques intermediaire pour le transport
[1632]97
[2263]98  !   variables pour le fichier histoire
99  REAL dtav      ! intervalle de temps elementaire
[1632]100
[2263]101  REAL time_0
[1632]102
[2263]103  LOGICAL lafin
[1632]104
[2263]105  real time_step, t_wrt, t_ops
[1632]106
[2263]107  !+jld variables test conservation energie
108  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
109  !     Tendance de la temp. potentiel d (theta)/ d t due a la
110  !     tansformation d'energie cinetique en energie thermique
111  !     cree par la dissipation
112  !      REAL dhecdt(ip1jmp1,llm)
113  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
114  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
115  !      CHARACTER (len=15) :: ztit
116  !-jld
[1632]117
118
[2263]119  character (len=80) :: dynhist_file, dynhistave_file
120  character (len=20) :: modname
121  character (len=80) :: abort_message
122  ! locales pour gestion du temps
123  INTEGER :: an, mois, jour
124  REAL :: heure
[1632]125
126
[2263]127  !-----------------------------------------------------------------------
128  !   Initialisations:
129  !   ----------------
[1632]130
[2263]131  abort_message = 'last timestep reached'
132  modname = 'gcm'
133  descript = 'Run GCM LMDZ'
134  lafin    = .FALSE.
135  dynhist_file = 'dyn_hist'
136  dynhistave_file = 'dyn_hist_ave'
[1632]137
138
139
[2263]140  !----------------------------------------------------------------------
141  !  lecture des fichiers gcm.def ou run.def
142  !  ---------------------------------------
143  !
144  CALL conf_gcm( 99, .TRUE. )
145  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
146       "iphysiq must be a multiple of iperiod", 1)
147  !
148  !
149  !------------------------------------
150  !   Initialisation partie parallele
151  !------------------------------------
152  CALL init_const_mpi
[1632]153
[2263]154  call init_parallel
155  call Read_Distrib
[1673]156
157#ifdef CPP_PHYS
[2263]158  CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
159  !#endif
160  !      CALL set_bands
161  !#ifdef CPP_PHYS
162  CALL Init_interface_dyn_phys
[1632]163#endif
[2263]164  CALL barrier
[1632]165
[2263]166  CALL set_bands
167  if (mpi_rank==0) call WriteBands
168  call Set_Distrib(distrib_caldyn)
[1632]169
[2263]170  !$OMP PARALLEL
171  call Init_Mod_hallo
172  !$OMP END PARALLEL
[1632]173
[2263]174  !#ifdef CPP_PHYS
175  !c$OMP PARALLEL
176  !      call InitComgeomphy ! now done in iniphysiq
177  !c$OMP END PARALLEL
178  !#endif
[1632]179
[2263]180  !-----------------------------------------------------------------------
181  !   Choix du calendrier
182  !   -------------------
[1632]183
[2263]184  !      calend = 'earth_365d'
[1632]185
186#ifdef CPP_IOIPSL
[2263]187  if (calend == 'earth_360d') then
188     call ioconf_calendar('360d')
189     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
190  else if (calend == 'earth_365d') then
191     call ioconf_calendar('noleap')
192     write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
193  else if (calend == 'gregorian') then
194     call ioconf_calendar('gregorian')
195     write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
196  else
197     abort_message = 'Mauvais choix de calendrier'
198     call abort_gcm(modname,abort_message,1)
199  endif
[1632]200#endif
201
[2263]202  IF (type_trac == 'inca') THEN
[1632]203#ifdef INCA
[2263]204     call init_const_lmdz( &
205          nbtr,anneeref,dayref, &
206          iphysiq,day_step,nday,  &
207          nbsrf, is_oce,is_sic, &
208          is_ter,is_lic, calend)
[1632]209
[2263]210     call init_inca_para( &
211          iim,jjm+1,llm,klon_glo,mpi_size, &
212          distrib_phys,COMM_LMDZ)
[1632]213#endif
[2263]214  END IF
[1632]215
[2263]216  !-----------------------------------------------------------------------
217  !   Initialisation des traceurs
218  !   ---------------------------
219  !  Choix du nombre de traceurs et du schema pour l'advection
220  !  dans fichier traceur.def, par default ou via INCA
221  call infotrac_init
[1632]222
[2263]223  ! Allocation de la tableau q : champs advectes   
224  ALLOCATE(ucov(ijb_u:ije_u,llm))
225  ALLOCATE(vcov(ijb_v:ije_v,llm))
226  ALLOCATE(teta(ijb_u:ije_u,llm))
227  ALLOCATE(masse(ijb_u:ije_u,llm))
228  ALLOCATE(ps(ijb_u:ije_u))
229  ALLOCATE(phis(ijb_u:ije_u))
230  ALLOCATE(q(ijb_u:ije_u,llm,nqtot))
[1632]231
[2263]232  !-----------------------------------------------------------------------
233  !   Lecture de l'etat initial :
234  !   ---------------------------
[1632]235
[2263]236  !  lecture du fichier start.nc
237  if (read_start) then
238     ! we still need to run iniacademic to initialize some
239     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
240     if (iflag_phys.ne.1) then
241        CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
242     endif
[1657]243
[2263]244     !        if (planet_type.eq."earth") then
245     ! Load an Earth-format start file
246     CALL dynetat0_loc("start.nc",vcov,ucov, &
247          teta,q,masse,ps,phis, time_0)
248     !        endif ! of if (planet_type.eq."earth")
[1673]249
[2263]250     !       write(73,*) 'ucov',ucov
251     !       write(74,*) 'vcov',vcov
252     !       write(75,*) 'teta',teta
253     !       write(76,*) 'ps',ps
254     !       write(77,*) 'q',q
[1632]255
[2263]256  endif ! of if (read_start)
[1632]257
[2263]258  ! le cas echeant, creation d un etat initial
259  IF (prt_level > 9) WRITE(lunout,*) &
260       'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
261  if (.not.read_start) then
262     CALL iniacademic_loc(vcov,ucov,teta,q,masse,ps,phis,time_0)
263  endif
[1632]264
[2263]265  !-----------------------------------------------------------------------
266  !   Lecture des parametres de controle pour la simulation :
267  !   -------------------------------------------------------
268  !  on recalcule eventuellement le pas de temps
[1632]269
[2263]270  IF(MOD(day_step,iperiod).NE.0) THEN
271     abort_message =  &
272          'Il faut choisir un nb de pas par jour multiple de iperiod'
273     call abort_gcm(modname,abort_message,1)
274  ENDIF
[1632]275
[2263]276  IF(MOD(day_step,iphysiq).NE.0) THEN
277     abort_message =  &
278          'Il faut choisir un nb de pas par jour multiple de iphysiq'
279     call abort_gcm(modname,abort_message,1)
280  ENDIF
[1632]281
[2263]282  zdtvr    = daysec/REAL(day_step)
283  IF(dtvr.NE.zdtvr) THEN
284     WRITE(lunout,*) &
285          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
286  ENDIF
[1632]287
[2263]288  !
289  ! on remet le calendrier \`a zero si demande
290  !
291  IF (start_time /= starttime) then
292     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
293          ,' fichier restart ne correspond pas a celle lue dans le run.def'
294     IF (raz_date == 1) then
295        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
296        start_time = starttime
297     ELSE
298        WRITE(lunout,*)'Je m''arrete'
299        CALL abort
300     ENDIF
301  ENDIF
302  IF (raz_date == 1) THEN
303     annee_ref = anneeref
304     day_ref = dayref
305     day_ini = dayref
306     itau_dyn = 0
307     itau_phy = 0
308     time_0 = 0.
309     write(lunout,*) &
310          'GCM: On reinitialise a la date lue dans gcm.def'
311  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
312     write(lunout,*) &
313          'GCM: Attention les dates initiales lues dans le fichier'
314     write(lunout,*) &
315          ' restart ne correspondent pas a celles lues dans '
316     write(lunout,*)' gcm.def'
317     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
318     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
319     write(lunout,*)' Pas de remise a zero'
320  ENDIF
321  !      if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
322  !        write(lunout,*)
323  !     .  'GCM: Attention les dates initiales lues dans le fichier'
324  !        write(lunout,*)
325  !     .  ' restart ne correspondent pas a celles lues dans '
326  !        write(lunout,*)' gcm.def'
327  !        write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
328  !        write(lunout,*)' day_ref=',day_ref," dayref=",dayref
329  !        if (raz_date .ne. 1) then
330  !          write(lunout,*)
331  !     .    'GCM: On garde les dates du fichier restart'
332  !        else
333  !          annee_ref = anneeref
334  !          day_ref = dayref
335  !          day_ini = dayref
336  !          itau_dyn = 0
337  !          itau_phy = 0
338  !          time_0 = 0.
339  !          write(lunout,*)
340  !     .   'GCM: On reinitialise a la date lue dans gcm.def'
341  !        endif
342  !      ELSE
343  !        raz_date = 0
344  !      endif
[1632]345
346#ifdef CPP_IOIPSL
[2263]347  mois = 1
348  heure = 0.
349  call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
350  jH_ref = jD_ref - int(jD_ref)
351  jD_ref = int(jD_ref)
[1632]352
[2263]353  call ioconf_startdate(INT(jD_ref), jH_ref)
[1632]354
[2263]355  write(lunout,*)'DEBUG'
356  write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
357  write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
358  call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
359  write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
360  write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
[1632]361#else
[2263]362  ! Ehouarn: we still need to define JD_ref and JH_ref
363  ! and since we don't know how many days there are in a year
364  ! we set JD_ref to 0 (this should be improved ...)
365  jD_ref=0
366  jH_ref=0
[1632]367#endif
368
[2263]369  if (iflag_phys.eq.1) then
370     ! these initialisations have already been done (via iniacademic)
371     ! if running in SW or Newtonian mode
372     !-----------------------------------------------------------------------
373     !   Initialisation des constantes dynamiques :
374     !   ------------------------------------------
375     dtvr = zdtvr
376     CALL iniconst
[1632]377
[2263]378     !-----------------------------------------------------------------------
379     !   Initialisation de la geometrie :
380     !   --------------------------------
381     CALL inigeom
[1632]382
[2263]383     !-----------------------------------------------------------------------
384     !   Initialisation du filtre :
385     !   --------------------------
386     CALL inifilr
387  endif ! of if (iflag_phys.eq.1)
388  !
389  !-----------------------------------------------------------------------
390  !   Initialisation de la dissipation :
391  !   ----------------------------------
[1632]392
[2263]393  CALL inidissip( lstardis, nitergdiv, nitergrot, niterh   , &
394       tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
[1632]395
[2263]396  !-----------------------------------------------------------------------
397  !   Initialisation de la physique :
398  !   -------------------------------
399  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
400     ! Physics:
[1673]401#ifdef CPP_PHYS
[2263]402     CALL iniphysiq(iim,jjm,llm,daysec,day_ini,dtphys/nsplit_phys, &
403          rlatu,rlonv,aire,cu,cv,rad,g,r,cpp, &
404          iflag_phys)
[1632]405#endif
[2263]406  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
[1632]407
408
[2263]409  !-----------------------------------------------------------------------
410  !   Initialisation des dimensions d'INCA :
411  !   --------------------------------------
412  IF (type_trac == 'inca') THEN
413     !$OMP PARALLEL
[1632]414#ifdef INCA
[2263]415     CALL init_inca_dim(klon_omp,llm,iim,jjm, &
416          rlonu,rlatu,rlonv,rlatv)
[1632]417#endif
[2263]418     !$OMP END PARALLEL
419  END IF
[1632]420
[2263]421  !-----------------------------------------------------------------------
422  !   Initialisation des I/O :
423  !   ------------------------
[1632]424
425
[2263]426  if (nday>=0) then
427     day_end = day_ini + nday
428  else
429     day_end = day_ini - nday/day_step
430  endif
[1632]431
[2263]432  WRITE(lunout,300)day_ini,day_end
433300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
434
[1632]435#ifdef CPP_IOIPSL
[2263]436  call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
437  write (lunout,301)jour, mois, an
438  call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
439  write (lunout,302)jour, mois, an
440301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
441302 FORMAT('1'/,15x,'    au ', i2,'/',i2,'/',i4)
[1632]442#endif
443
[2263]444  !      if (planet_type.eq."earth") then
445  ! Write an Earth-format restart file
446  CALL dynredem0_loc("restart.nc", day_end, phis)
447  !      endif
[1632]448
[2263]449  ecripar = .TRUE.
[1632]450
451#ifdef CPP_IOIPSL
[2263]452  time_step = zdtvr
453  IF (mpi_rank==0) then
454     if (ok_dyn_ins) then
455        ! initialize output file for instantaneous outputs
456        ! t_ops = iecri * daysec ! do operations every t_ops
457        t_ops =((1.0*iecri)/day_step) * daysec 
458        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
459        t_wrt = daysec ! iecri * daysec ! write output every t_wrt
460        CALL inithist(day_ref,annee_ref,time_step, &
461             t_ops,t_wrt)
462     endif
[1632]463
[2263]464     IF (ok_dyn_ave) THEN
465        ! initialize output file for averaged outputs
466        t_ops = iperiod * time_step ! do operations every t_ops
467        t_wrt = periodav * daysec   ! write output every t_wrt
468        CALL initdynav_loc(day_ref,annee_ref,time_step,t_ops,t_wrt)
469     END IF
470  ENDIF
471  dtav = iperiod*dtvr/daysec
[1632]472#endif
[2263]473  ! #endif of #ifdef CPP_IOIPSL
[1632]474
[2263]475  !  Choix des frequences de stokage pour le offline
476  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
477  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
478  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
479  istphy=istdyn/iphysiq     
[1632]480
481
[2263]482  !
483  !-----------------------------------------------------------------------
484  !   Integration temporelle du modele :
485  !   ----------------------------------
[1632]486
[2263]487  !       write(78,*) 'ucov',ucov
488  !       write(78,*) 'vcov',vcov
489  !       write(78,*) 'teta',teta
490  !       write(78,*) 'ps',ps
491  !       write(78,*) 'q',q
[1632]492
[2263]493  !$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logici/,/logicl/)
494  CALL leapfrog_loc(ucov,vcov,teta,ps,masse,phis,q,time_0)
495  !$OMP END PARALLEL
[1632]496
[2263]497  !      OPEN(unit=5487,file='ok_lmdz',status='replace')
498  !      WRITE(5487,*) 'ok_lmdz'
499  !      CLOSE(5487)
500END PROGRAM gcm
Note: See TracBrowser for help on using the repository browser.