source: trunk/LMDZ.PLUTO/libf/dynphy_lonlat/phypluto/nogcm_pluto.F90 @ 3595

Last change on this file since 3595 was 3548, checked in by tbertrand, 11 months ago

LMDZ.PLUTO
Adding alternative routines for nogcm that better conserve the mass (but overall equivalent to nogcm). We want to keep both options (nogcm and nogcm_pluto) for now.
TB

File size: 13.1 KB
RevLine 
[3548]1!
2! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $
3!
4! En cours de modification par T. Bertrand
5!
6PROGRAM nogcm_pluto
7
8#ifdef CPP_IOIPSL
9  USE IOIPSL
10#else
11  ! if not using IOIPSL, we still need to use (a local version of) getin
12  USE ioipsl_getincom
13#endif
14
15
16#ifdef CPP_XIOS
17  ! ug Pour les sorties XIOS
18  USE wxios
19#endif
20
21  USE filtreg_mod
22  USE infotrac
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 mod_const_mpi, ONLY: COMM_LMDZ
28  use cpdet_mod, only: ini_cpdet
29  USE temps_mod, ONLY: calend,start_time,annee_ref,day_ref, &
30                itau_dyn,itau_phy,day_ini,jD_ref,jH_ref,day_end
31
32
33!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
35! A nettoyer. On ne veut qu'une ou deux routines d'interface
36! dynamique -> physique pour l'initialisation
37! Ehouarn: the following are needed with (parallel) physics:
38#ifdef CPP_PHYS
39  USE iniphysiq_mod, ONLY: iniphysiq
40#endif
41
42  USE comconst_mod, ONLY: daysec,dtvr,dtphys,rad,g,r,cpp
43  USE logic_mod, ONLY: read_start,iflag_phys,ok_guide,ecripar
44
45!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
46
47  IMPLICIT NONE
48
49  !      ......   Version  du 10/01/98    ..........
50
51  !             avec  coordonnees  verticales hybrides
52  !   avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
53
54  !=======================================================================
55  !
56  !   Auteur:  P. Le Van /L. Fairhead/F.Hourdin
57  !   -------
58  !
59  !   Objet:
60  !   ------
61  !
62  !   GCM LMD nouvelle grille
63  !
64  !=======================================================================
65  !
66  !  ... Dans inigeom , nouveaux calculs pour les elongations  cu , cv
67  !      et possibilite d'appeler une fonction f(y)  a derivee tangente
68  !      hyperbolique a la  place de la fonction a derivee sinusoidale.
69  !  ... Possibilite de choisir le schema pour l'advection de
70  !        q  , en modifiant iadv dans traceur.def  (MAF,10/02) .
71  !
72  !      Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
73  !      Pour Van-Leer iadv=10
74  !
75  !-----------------------------------------------------------------------
76  !   Declarations:
77  !   -------------
78
79  include "dimensions.h"
80  include "paramet.h"
81  include "comdissnew.h"
82  include "comgeom.h"
83!!!!!!!!!!!#include "control.h"
84!#include "com_io_dyn.h"
85  include "iniprint.h"
86  include "tracstoke.h"
87
88  REAL zdtvr
89
90  !   variables dynamiques
91  REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
92  REAL teta(ip1jmp1,llm)                 ! temperature potentielle
93  REAL, ALLOCATABLE, DIMENSION(:,:,:):: q! champs advectes
94  REAL ps(ip1jmp1)                       ! pression  au sol
95  REAL p (ip1jmp1,llmp1  )               ! pression aux interfac.des couches
96  REAL masse(ip1jmp1,llm)                ! masse d'air
97  REAL phis(ip1jmp1)                     ! geopotentiel au sol
98  REAL phi(ip1jmp1,llm)                  ! geopotentiel
99  REAL w(ip1jmp1,llm)                    ! vitesse verticale
100
101  ! variables dynamiques intermediaire pour le transport
102
103  !   variables pour le fichier histoire
104  REAL dtav      ! intervalle de temps elementaire
105
106  REAL time_0
107
108  LOGICAL lafin
109  INTEGER ij,iq,l,i,j
110
111
112  real time_step, t_wrt, t_ops
113
114  LOGICAL first
115
116  !      LOGICAL call_iniphys
117  !      data call_iniphys/.true./
118
119  !+jld variables test conservation energie
120  !      REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
121  !     Tendance de la temp. potentiel d (theta)/ d t due a la
122  !     tansformation d'energie cinetique en energie thermique
123  !     cree par la dissipation
124  REAL dhecdt(ip1jmp1,llm)
125  !      REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
126  !      REAL      d_h_vcol, d_qt, d_qw, d_ql, d_ec
127  CHARACTER (len=15) :: ztit
128  !-jld
129
130
131  character (len=80) :: dynhist_file, dynhistave_file
132  character (len=20) :: modname
133  character (len=80) :: abort_message
134  ! locales pour gestion du temps
135  INTEGER :: an, mois, jour
136  REAL :: heure
137  logical use_filtre_fft
138
139!-----------------------------------------------------------------------
140!   Initialisations:
141!   ----------------
142
143  abort_message = 'last timestep reached'
144  modname = 'gcm'
145  lafin    = .FALSE.
146  dynhist_file = 'dyn_hist.nc'
147  dynhistave_file = 'dyn_hist_ave.nc'
148
149
150
151!----------------------------------------------------------------------
152!  lecture des fichiers gcm.def ou run.def
153!  ---------------------------------------
154!
155  CALL conf_gcm( 99, .TRUE. )
156  if (mod(iphysiq, iperiod) /= 0) call abort_gcm("conf_gcm", &
157       "iphysiq must be a multiple of iperiod", 1)
158
159  use_filtre_fft=.FALSE.
160  CALL getin('use_filtre_fft',use_filtre_fft)
161  IF (use_filtre_fft) call abort_gcm("gcm",'FFT filter is not available in the ' &
162          // 'sequential version of the dynamics.', 1)
163
164!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
165! Initialisation de XIOS
166!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
167
168#ifdef CPP_XIOS
169  CALL wxios_init("LMDZ")
170#endif
171
172
173!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
174! FH 2008/05/02
175! A nettoyer. On ne veut qu'une ou deux routines d'interface
176! dynamique -> physique pour l'initialisation
177!#ifdef CPP_PHYS
178!  CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
179!!      call initcomgeomphy ! now done in iniphysiq
180!#endif
181!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
182!
183! Initialisations pour Cp(T) Venus
184  call ini_cpdet
185!
186!-----------------------------------------------------------------------
187!   Choix du calendrier
188!   -------------------
189
190!      calend = 'earth_365d'
191
192#ifdef CPP_IOIPSL
193  if (calend == 'earth_360d') then
194    call ioconf_calendar('360d')
195    write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
196  else if (calend == 'earth_365d') then
197    call ioconf_calendar('noleap')
198    write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
199  else if (calend == 'earth_366d') then
200    call ioconf_calendar('gregorian')
201    write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
202  else
203    abort_message = 'Mauvais choix de calendrier'
204    call abort_gcm(modname,abort_message,1)
205  endif
206#endif
207!-----------------------------------------------------------------------
208  !
209  !
210  !------------------------------------
211  !   Initialisation partie parallele
212  !------------------------------------
213
214  !
215  !
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
222
223  ! Allocation de la tableau q : champs advectes   
224  allocate(q(ip1jmp1,llm,nqtot))
225
226  !-----------------------------------------------------------------------
227  !   Lecture de l'etat initial :
228  !   ---------------------------
229
230  !  lecture du fichier start.nc
231  if (read_start) then
232     ! we still need to run iniacademic to initialize some
233     ! constants & fields, if we run the 'newtonian' or 'SW' cases:
234     if (iflag_phys.ne.1) then
235        CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
236     endif
237
238     CALL dynetat0("start.nc",vcov,ucov, &
239                    teta,q,masse,ps,phis, time_0)
240       
241     !       write(73,*) 'ucov',ucov
242     !       write(74,*) 'vcov',vcov
243     !       write(75,*) 'teta',teta
244     !       write(76,*) 'ps',ps
245     !       write(77,*) 'q',q
246
247  endif ! of if (read_start)
248
249
250  ! le cas echeant, creation d un etat initial
251  IF (prt_level > 9) WRITE(lunout,*) &
252       'NOGCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
253  if (.not.read_start) then
254     CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
255  endif
256
257
258  !-----------------------------------------------------------------------
259  !   Lecture des parametres de controle pour la simulation :
260  !   -------------------------------------------------------
261  !  on recalcule eventuellement le pas de temps
262
263  IF(MOD(day_step,iperiod).NE.0) THEN
264     abort_message = &
265       'Il faut choisir un nb de pas par jour multiple de iperiod'
266     call abort_gcm(modname,abort_message,1)
267  ENDIF
268
269!  IF(MOD(day_step,iphysiq).NE.0) THEN
270!     abort_message = &
271!       'Il faut choisir un nb de pas par jour multiple de iphysiq'
272!     call abort_gcm(modname,abort_message,1)
273!  ENDIF
274
275  zdtvr    = daysec/REAL(day_step)
276  IF(dtvr.NE.zdtvr) THEN
277     WRITE(lunout,*) &
278          'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
279  ENDIF
280
281  !
282  ! on remet le calendrier à zero si demande
283  !
284  IF (start_time /= starttime) then
285     WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le' &
286     ,' fichier restart ne correspond pas à celle lue dans le run.def'
287     IF (raz_date == 1) then
288        WRITE(lunout,*)'Je prends l''heure lue dans run.def'
289        start_time = starttime
290     ELSE
291        call abort_gcm("gcm", "'Je m''arrete'", 1)
292     ENDIF
293  ENDIF
294  IF (raz_date == 1) THEN
295     annee_ref = anneeref
296     day_ref = dayref
297     day_ini = dayref
298     itau_dyn = 0
299     itau_phy = 0
300     time_0 = 0.
301     write(lunout,*) &
302         'GCM: On reinitialise a la date lue dans gcm.def'
303  ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
304     write(lunout,*) &
305        'GCM: Attention les dates initiales lues dans le fichier'
306     write(lunout,*) &
307        ' restart ne correspondent pas a celles lues dans '
308     write(lunout,*)' gcm.def'
309     write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
310     write(lunout,*)' day_ref=',day_ref," dayref=",dayref
311     write(lunout,*)' Pas de remise a zero'
312  ENDIF
313
314#ifdef CPP_IOIPSL
315  mois = 1
316  heure = 0.
317  ! Ehouarn: we still need to define JD_ref and JH_ref
318  ! and since we don't know how many days there are in a year
319  ! we set JD_ref to 0 (this should be improved ...)
320  jD_ref=0
321  jH_ref=0
322#endif
323
324  if (iflag_phys.eq.1) then
325     ! these initialisations have already been done (via iniacademic)
326     ! if running in SW or Newtonian mode
327     !-----------------------------------------------------------------------
328     !   Initialisation des constantes dynamiques :
329     !   ------------------------------------------
330     dtvr = zdtvr
331     CALL iniconst
332
333     !-----------------------------------------------------------------------
334     !   Initialisation de la geometrie :
335     !   --------------------------------
336     CALL inigeom
337
338     !-----------------------------------------------------------------------
339     !   Initialisation du filtre :
340     !   --------------------------
341     CALL inifilr
342  endif ! of if (iflag_phys.eq.1)
343  !
344  !-----------------------------------------------------------------------
345  !   Initialisation des I/O :
346  !   ------------------------
347
348  if (nday>=0) then ! standard case
349     day_end=day_ini+nday
350  else ! special case when nday <0, run -nday dynamical steps
351     day_end=day_ini-nday/day_step
352  endif
353  if (less1day) then
354     day_end=day_ini+floor(time_0+fractday)
355  endif
356  if (ndynstep.gt.0) then
357     day_end=day_ini+floor(time_0+float(ndynstep)/float(day_step))
358  endif
359     
360  WRITE(lunout,'(a,i7,a,i7)') &
361               "run from day ",day_ini,"  to day",day_end
362
363  !-----------------------------------------------------------------------
364  !   Initialisation de la physique :
365  !   -------------------------------
366
367  IF ((iflag_phys==1).or.(iflag_phys>=100)) THEN
368     ! Physics:
369#ifdef CPP_PHYS
370     CALL iniphysiq(iim,jjm,llm, &
371          (jjm-1)*iim+2,comm_lmdz, &
372          daysec,day_ini,dtphys/nsplit_phys, &
373          rlatu,rlatv,rlonu,rlonv,aire,cu,cv,rad,g,r,cpp, &
374          iflag_phys)
375#endif
376  ENDIF ! of IF ((iflag_phys==1).or.(iflag_phys>=100))
377
378
379  CALL dynredem0("restart.nc", day_end, phis)
380  ecripar = .TRUE.
381
382#ifdef CPP_IOIPSL
383  time_step = zdtvr
384  if (ok_dyn_ins) then
385    ! initialize output file for instantaneous outputs
386    ! t_ops = iecri * daysec ! do operations every t_ops
387    t_ops =((1.0*iecri)/day_step) * daysec 
388    t_wrt = daysec ! iecri * daysec ! write output every t_wrt
389    CALL inithist(day_ref,annee_ref,time_step, &
390                   t_ops,t_wrt)
391  endif
392
393  IF (ok_dyn_ave) THEN
394    ! initialize output file for averaged outputs
395    t_ops = iperiod * time_step ! do operations every t_ops
396    t_wrt = periodav * daysec   ! write output every t_wrt
397    CALL initdynav(day_ref,annee_ref,time_step, &
398            t_ops,t_wrt)
399  END IF
400  dtav = iperiod*dtvr/daysec
401#endif
402! #endif of #ifdef CPP_IOIPSL
403
404  !  Choix des frequences de stokage pour le offline
405  !      istdyn=day_step/4     ! stockage toutes les 6h=1jour/4
406  !      istdyn=day_step/12     ! stockage toutes les 2h=1jour/12
407  istdyn=day_step/4     ! stockage toutes les 6h=1jour/12
408  istphy=istdyn/iphysiq     
409
410
411  !
412  !-----------------------------------------------------------------------
413  !   Integration temporelle du modele :
414  !   ----------------------------------
415
416  !       write(78,*) 'ucov',ucov
417  !       write(78,*) 'vcov',vcov
418  !       write(78,*) 'teta',teta
419  !       write(78,*) 'ps',ps
420  !       write(78,*) 'q',q
421
422
423  CALL leapfrog_nogcm_pluto(ucov,vcov,teta,ps,masse,phis,q,time_0)
424
425END PROGRAM nogcm_pluto
426
427
Note: See TracBrowser for help on using the repository browser.