Context Navigation

gcm.F @ 1202

Last change on this file since 1202 was 1201, checked in by Laurent Fairhead, 15 years ago
Modifications nécessaires a l'inclusion d'un calendrier réaliste. La date courante est calculée dans leapfrog.F et exprimée en Jour Julien (modifié). On en a profité pour faire un peu de ménage dans la gestion des dates du modèle. Dans la physique, on utilise les routines de passages entre calendrier Julien et Gregorien incluses dans IOIPSL pour calculer le nombre de jours écoulés depuis le 1er janvier (pour les conditions aux limites) ou l'equinoxe (pour le calcul de la longitude solaire). Le calcul de l'orbite reprend celui du gcm planétaire (codé par FH) On décide du calendrier à utiliser à l'aide du paramètre calend du run.def. Par défaut celui-ci est à earth_360d LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 14.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format