source: LMDZ6/trunk/libf/dyn3dmem/gcm.F90 @ 5229

Last change on this file since 5229 was 4996, checked in by evignon, 12 months ago

ajout d'un flag pour le calcul de qsat dans la condtion de "francis"
pour l'advection de l'humidite (q<qsat_aval). En activant ce flag,
on calcule qsat /liquide quelque soit la temperature et on peut donc
ainsi autoriser l'advection de sursaturations / glace à T<0oC.

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