Context Navigation

gcm.F @ 1017

Last change on this file since 1017 was 1017, checked in by emillour, 11 years ago

Common dynamics: (and collateral adaptations in Venus physics)
Improved cpdet routines in and additional sponge mode:

Additionnal sponge mode (trigered with "callsponge" flag), in line with the one used in the Generic and Martian GCM. This sponge is called whenever there is a dissipation step.
Improvement of the cpdet routines : created routines tpot2t_glo_p and t2tpot_glo_p which handle fields on the whole dynamics (scaler) grid, which are more efficient than calling tpot2t_p or t2tpot_p with slabs of data (generated use of intermediate copies of these chunks of data at every call)
Turned cpdet.F into a module cpdet_mod.F90 (and correspondingly adapted all routines in the Venus physics).

File size: 17.6 KB

Line
1	!
2	! $Id: gcm.F 1446 2010-10-22 09:27:25Z emillour $
3	!
4	c
5	c
6	PROGRAM gcm
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	USE filtreg_mod
16	USE infotrac
17	USE control_mod
18	use cpdet_mod, only: ini_cpdet
19
20	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
21	! FH 2008/05/09 On elimine toutes les clefs physiques dans la dynamique
22	! A nettoyer. On ne veut qu'une ou deux routines d'interface
23	! dynamique -> physique pour l'initialisation
24	! Ehouarn: the following are needed with (parallel) physics:
25	#ifdef CPP_PHYS
26	USE dimphy
27	USE comgeomphy
28	#endif
29	#ifdef INCA
30	USE mod_phys_lmdz_para, ONLY : klon_mpi_para_nb
31	#endif
32	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
33
34	IMPLICIT NONE
35
36	c ...... Version du 10/01/98 ..........
37
38	c avec coordonnees verticales hybrides
39	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
40
41	c=======================================================================
42	c
43	c Auteur: P. Le Van /L. Fairhead/F.Hourdin
44	c -------
45	c
46	c Objet:
47	c ------
48	c
49	c GCM LMD nouvelle grille
50	c
51	c=======================================================================
52	c
53	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
54	c et possibilite d'appeler une fonction f(y) a derivee tangente
55	c hyperbolique a la place de la fonction a derivee sinusoidale.
56	c ... Possibilite de choisir le schema pour l'advection de
57	c q , en modifiant iadv dans traceur.def (MAF,10/02) .
58	c
59	c Pour Van-Leer + Vapeur d'eau saturee, iadv(1)=4. (F.Codron,10/99)
60	c Pour Van-Leer iadv=10
61	c
62	c-----------------------------------------------------------------------
63	c Declarations:
64	c -------------
65
66	#include "dimensions.h"
67	#include "paramet.h"
68	#include "comconst.h"
69	#include "comdissnew.h"
70	#include "comvert.h"
71	#include "comgeom.h"
72	#include "logic.h"
73	#include "temps.h"
74	!!!!!!!!!!!#include "control.h"
75	#include "ener.h"
76	#include "description.h"
77	#include "serre.h"
78	!#include "com_io_dyn.h"
79	#include "iniprint.h"
80	#include "tracstoke.h"
81	#ifdef INCA
82	! Only INCA needs these informations (from the Earth's physics)
83	#include "indicesol.h"
84	#endif
85
86
87	REAL zdtvr
88
89	c 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
94	REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
95	REAL pks(ip1jmp1) ! exner au sol
96	REAL pk(ip1jmp1,llm) ! exner au milieu des couches
97	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
100	REAL phi(ip1jmp1,llm) ! geopotentiel
101	REAL w(ip1jmp1,llm) ! vitesse verticale
102
103	c variables dynamiques intermediaire pour le transport
104
105	c variables pour le fichier histoire
106	REAL dtav ! intervalle de temps elementaire
107
108	REAL time_0
109
110	LOGICAL lafin
111	INTEGER ij,iq,l,i,j
112
113
114	real time_step, t_wrt, t_ops
115
116	LOGICAL first
117
118	LOGICAL call_iniphys
119	data call_iniphys/.true./
120
121	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
122	c+jld variables test conservation energie
123	c REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
124	C Tendance de la temp. potentiel d (theta)/ d t due a la
125	C tansformation d'energie cinetique en energie thermique
126	C cree par la dissipation
127	REAL dhecdt(ip1jmp1,llm)
128	c REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
129	c REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
130	CHARACTER (len=15) :: ztit
131	c-jld
132
133
134	character (len=80) :: dynhist_file, dynhistave_file
135	character (len=20) :: modname
136	character (len=80) :: abort_message
137	! locales pour gestion du temps
138	INTEGER :: an, mois, jour
139	REAL :: heure
140
141
142	c-----------------------------------------------------------------------
143	c variables pour l'initialisation de la physique :
144	c ------------------------------------------------
145	INTEGER ngridmx
146	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
147	REAL zcufi(ngridmx),zcvfi(ngridmx)
148	REAL latfi(ngridmx),lonfi(ngridmx)
149	REAL airefi(ngridmx)
150	SAVE latfi, lonfi, airefi
151
152	c-----------------------------------------------------------------------
153	c Initialisations:
154	c ----------------
155
156	abort_message = 'last timestep reached'
157	modname = 'gcm'
158	descript = 'Run GCM LMDZ'
159	lafin = .FALSE.
160	dynhist_file = 'dyn_hist.nc'
161	dynhistave_file = 'dyn_hist_ave.nc'
162
163
164
165	c----------------------------------------------------------------------
166	c lecture des fichiers gcm.def ou run.def
167	c ---------------------------------------
168	c
169	! Ehouarn: dump possibility of using defrun
170	!#ifdef CPP_IOIPSL
171	CALL conf_gcm( 99, .TRUE. )
172	!#else
173	! CALL defrun( 99, .TRUE. , clesphy0 )
174	!#endif
175
176	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
177	! FH 2008/05/02
178	! A nettoyer. On ne veut qu'une ou deux routines d'interface
179	! dynamique -> physique pour l'initialisation
180	#ifdef CPP_PHYS
181	CALL init_phys_lmdz(iim,jjp1,llm,1,(/(jjm-1)*iim+2/))
182	call initcomgeomphy
183	#endif
184	!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
185	c
186	c Initialisations pour Cp(T) Venus
187	call ini_cpdet
188	c
189	c-----------------------------------------------------------------------
190	c Choix du calendrier
191	c -------------------
192
193	c calend = 'earth_365d'
194
195	#ifdef CPP_IOIPSL
196	if (calend == 'earth_360d') then
197	call ioconf_calendar('360d')
198	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
199	else if (calend == 'earth_365d') then
200	call ioconf_calendar('noleap')
201	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
202	else if (calend == 'earth_366d') then
203	call ioconf_calendar('gregorian')
204	write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
205	else if (calend == 'titan') then
206	! call ioconf_calendar('titan')
207	write(lunout,*)'CALENDRIER CHOISI: Titan'
208	abort_message = 'A FAIRE...'
209	call abort_gcm(modname,abort_message,1)
210	else if (calend == 'venus') then
211	! call ioconf_calendar('venus')
212	write(lunout,*)'CALENDRIER CHOISI: Venus'
213	abort_message = 'A FAIRE...'
214	call abort_gcm(modname,abort_message,1)
215	else
216	abort_message = 'Mauvais choix de calendrier'
217	call abort_gcm(modname,abort_message,1)
218	endif
219	#endif
220	c-----------------------------------------------------------------------
221
222	IF (type_trac == 'inca') THEN
223	#ifdef INCA
224	call init_const_lmdz(nbtr,anneeref,dayref,iphysiq,day_step,nday,
225	$ nbsrf, is_oce,is_sic,is_ter,is_lic)
226	call init_inca_para(iim,jjm+1,klon,1,klon_mpi_para_nb,0)
227	#endif
228	END IF
229	c
230	c
231	c------------------------------------
232	c Initialisation partie parallele
233	c------------------------------------
234
235	c
236	c
237	c-----------------------------------------------------------------------
238	c Initialisation des traceurs
239	c ---------------------------
240	c Choix du nombre de traceurs et du schema pour l'advection
241	c dans fichier traceur.def, par default ou via INCA
242	call infotrac_init
243
244	c Allocation de la tableau q : champs advectes
245	allocate(q(ip1jmp1,llm,nqtot))
246
247	c-----------------------------------------------------------------------
248	c Lecture de l'etat initial :
249	c ---------------------------
250
251	c lecture du fichier start.nc
252	if (read_start) then
253	! we still need to run iniacademic to initialize some
254	! constants & fields, if we run the 'newtonian' or 'SW' cases:
255	if (iflag_phys.ne.1) then
256	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
257	endif
258
259	if (planet_type.eq."mars") then
260	! POUR MARS, METTRE UNE FONCTION A PART, genre dynetat0_mars
261	abort_message = 'dynetat0_mars A FAIRE'
262	call abort_gcm(modname,abort_message,0)
263	else
264	CALL dynetat0("start.nc",vcov,ucov,
265	& teta,q,masse,ps,phis, time_0)
266	endif ! of if (planet_type.eq."mars")
267
268	c write(73,*) 'ucov',ucov
269	c write(74,*) 'vcov',vcov
270	c write(75,*) 'teta',teta
271	c write(76,*) 'ps',ps
272	c write(77,*) 'q',q
273
274	endif ! of if (read_start)
275
276	IF (type_trac == 'inca') THEN
277	#ifdef INCA
278	call init_inca_dim(klon,llm,iim,jjm,
279	$ rlonu,rlatu,rlonv,rlatv)
280	#endif
281	END IF
282
283
284	c le cas echeant, creation d un etat initial
285	IF (prt_level > 9) WRITE(lunout,*)
286	. 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
287	if (.not.read_start) then
288	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
289	endif
290
291
292	c-----------------------------------------------------------------------
293	c Lecture des parametres de controle pour la simulation :
294	c -------------------------------------------------------
295	c on recalcule eventuellement le pas de temps
296
297	IF(MOD(day_step,iperiod).NE.0) THEN
298	abort_message =
299	. 'Il faut choisir un nb de pas par jour multiple de iperiod'
300	call abort_gcm(modname,abort_message,1)
301	ENDIF
302
303	IF(MOD(day_step,iphysiq).NE.0) THEN
304	abort_message =
305	* 'Il faut choisir un nb de pas par jour multiple de iphysiq'
306	call abort_gcm(modname,abort_message,1)
307	ENDIF
308
309	zdtvr = daysec/REAL(day_step)
310	IF(dtvr.NE.zdtvr) THEN
311	WRITE(lunout,*)
312	. 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
313	ENDIF
314
315	C
316	C on remet le calendrier à zero si demande
317	c
318	IF (start_time /= starttime) then
319	WRITE(lunout,*)' GCM: Attention l''heure de depart lue dans le'
320	&,' fichier restart ne correspond pas à celle lue dans le run.def'
321	IF (raz_date == 1) then
322	WRITE(lunout,*)'Je prends l''heure lue dans run.def'
323	start_time = starttime
324	ELSE
325	WRITE(lunout,*)'Je m''arrete'
326	CALL abort
327	ENDIF
328	ENDIF
329	IF (raz_date == 1) THEN
330	annee_ref = anneeref
331	day_ref = dayref
332	day_ini = dayref
333	itau_dyn = 0
334	itau_phy = 0
335	time_0 = 0.
336	write(lunout,*)
337	. 'GCM: On reinitialise a la date lue dans gcm.def'
338	ELSE IF (annee_ref .ne. anneeref .or. day_ref .ne. dayref) THEN
339	write(lunout,*)
340	. 'GCM: Attention les dates initiales lues dans le fichier'
341	write(lunout,*)
342	. ' restart ne correspondent pas a celles lues dans '
343	write(lunout,*)' gcm.def'
344	write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
345	write(lunout,*)' day_ref=',day_ref," dayref=",dayref
346	write(lunout,*)' Pas de remise a zero'
347	ENDIF
348
349	c if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
350	c write(lunout,*)
351	c . 'GCM: Attention les dates initiales lues dans le fichier'
352	c write(lunout,*)
353	c . ' restart ne correspondent pas a celles lues dans '
354	c write(lunout,*)' gcm.def'
355	c write(lunout,*)' annee_ref=',annee_ref," anneeref=",anneeref
356	c write(lunout,*)' day_ref=',day_ref," dayref=",dayref
357	c if (raz_date .ne. 1) then
358	c write(lunout,*)
359	c . 'GCM: On garde les dates du fichier restart'
360	c else
361	c annee_ref = anneeref
362	c day_ref = dayref
363	c day_ini = dayref
364	c itau_dyn = 0
365	c itau_phy = 0
366	c time_0 = 0.
367	c write(lunout,*)
368	c . 'GCM: On reinitialise a la date lue dans gcm.def'
369	c endif
370	c ELSE
371	c raz_date = 0
372	c endif
373
374	#ifdef CPP_IOIPSL
375	mois = 1
376	heure = 0.
377	! Ce n'est defini pour l'instant que pour la Terre...
378	if (planet_type.eq.'earth') then
379	call ymds2ju(annee_ref, mois, day_ref, heure, jD_ref)
380	jH_ref = jD_ref - int(jD_ref)
381	jD_ref = int(jD_ref)
382
383	call ioconf_startdate(INT(jD_ref), jH_ref)
384
385	write(lunout,*)'DEBUG'
386	write(lunout,*)'annee_ref, mois, day_ref, heure, jD_ref'
387	write(lunout,*)annee_ref, mois, day_ref, heure, jD_ref
388	call ju2ymds(jD_ref+jH_ref,an, mois, jour, heure)
389	write(lunout,*)'jD_ref+jH_ref,an, mois, jour, heure'
390	write(lunout,*)jD_ref+jH_ref,an, mois, jour, heure
391	else
392	! A voir pour Titan et Venus
393	jD_ref=0
394	jH_ref=0
395	write(lunout,*)'A VOIR POUR VENUS ET TITAN: jD_ref, jH_ref'
396	write(lunout,*)jD_ref,jH_ref
397	endif ! planet_type
398	#else
399	! Ehouarn: we still need to define JD_ref and JH_ref
400	! and since we don't know how many days there are in a year
401	! we set JD_ref to 0 (this should be improved ...)
402	jD_ref=0
403	jH_ref=0
404	#endif
405
406	if (iflag_phys.eq.1) then
407	! these initialisations have already been done (via iniacademic)
408	! if running in SW or Newtonian mode
409	c-----------------------------------------------------------------------
410	c Initialisation des constantes dynamiques :
411	c ------------------------------------------
412	dtvr = zdtvr
413	CALL iniconst
414
415	c-----------------------------------------------------------------------
416	c Initialisation de la geometrie :
417	c --------------------------------
418	CALL inigeom
419
420	c-----------------------------------------------------------------------
421	c Initialisation du filtre :
422	c --------------------------
423	CALL inifilr
424	endif ! of if (iflag_phys.eq.1)
425	c
426	c-----------------------------------------------------------------------
427	c Initialisation de la dissipation :
428	c ----------------------------------
429
430	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
431	* tetagdiv, tetagrot , tetatemp, vert_prof_dissip)
432
433	c-----------------------------------------------------------------------
434	c Initialisation de la physique :
435	c -------------------------------
436
437	IF (call_iniphys.and.(iflag_phys==1.or.iflag_phys>=100)) THEN
438	latfi(1)=rlatu(1)
439	lonfi(1)=0.
440	zcufi(1) = cu(1)
441	zcvfi(1) = cv(1)
442	DO j=2,jjm
443	DO i=1,iim
444	latfi((j-2)*iim+1+i)= rlatu(j)
445	lonfi((j-2)*iim+1+i)= rlonv(i)
446	zcufi((j-2)iim+1+i) = cu((j-1)iip1+i)
447	zcvfi((j-2)iim+1+i) = cv((j-1)iip1+i)
448	ENDDO
449	ENDDO
450	latfi(ngridmx)= rlatu(jjp1)
451	lonfi(ngridmx)= 0.
452	zcufi(ngridmx) = cu(ip1jm+1)
453	zcvfi(ngridmx) = cv(ip1jm-iim)
454
455	! build airefi(), mesh area on physics grid
456	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
457	! Poles are single points on physics grid
458	airefi(1)=airefi(1)*iim
459	airefi(ngridmx)=airefi(ngridmx)*iim
460
461	! Initialisation de la physique: pose probleme quand on tourne
462	! SANS physique, car iniphysiq.F est dans le repertoire phy[]...
463	! Il faut une cle CPP_PHYS
464	#ifdef CPP_PHYS
465	CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys,
466	& latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp,
467	& iflag_phys)
468	#endif
469	call_iniphys=.false.
470	ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
471
472	c numero de stockage pour les fichiers de redemarrage:
473
474	c-----------------------------------------------------------------------
475	c Initialisation des I/O :
476	c ------------------------
477
478
479	day_end = day_ini + nday
480	WRITE(lunout,300)day_ini,day_end
481	300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
482
483	#ifdef CPP_IOIPSL
484	! Ce n'est defini pour l'instant que pour la Terre...
485	if (planet_type.eq.'earth') then
486	call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
487	write (lunout,301)jour, mois, an
488	call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
489	write (lunout,302)jour, mois, an
490	else
491	! A voir pour Titan et Venus
492	write(lunout,*)'A VOIR POUR VENUS/TITAN: separation en annees...'
493	endif ! planet_type
494
495	301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
496	302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4)
497	#endif
498
499	if (planet_type.eq."mars") then
500	! POUR MARS, METTRE UNE FONCTION A PART, genre dynredem0_mars
501	abort_message = 'dynredem0_mars A FAIRE'
502	call abort_gcm(modname,abort_message,0)
503	else
504	CALL dynredem0("restart.nc", day_end, phis)
505	endif ! of if (planet_type.eq."mars")
506
507	ecripar = .TRUE.
508
509	#ifdef CPP_IOIPSL
510	time_step = zdtvr
511	if (ok_dyn_ins) then
512	! initialize output file for instantaneous outputs
513	! t_ops = iecri * daysec ! do operations every t_ops
514	t_ops =((1.0iecri)/day_step) daysec
515	t_wrt = daysec ! iecri * daysec ! write output every t_wrt
516	CALL inithist(day_ref,annee_ref,time_step,
517	& t_ops,t_wrt)
518	endif
519
520	IF (ok_dyn_ave) THEN
521	! initialize output file for averaged outputs
522	t_ops = iperiod * time_step ! do operations every t_ops
523	t_wrt = periodav * daysec ! write output every t_wrt
524	CALL initdynav(day_ref,annee_ref,time_step,
525	& t_ops,t_wrt)
526	END IF
527	dtav = iperiod*dtvr/daysec
528	#endif
529	! #endif of #ifdef CPP_IOIPSL
530
531	c Choix des frequences de stokage pour le offline
532	c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
533	c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12
534	istdyn=day_step/4 ! stockage toutes les 6h=1jour/12
535	istphy=istdyn/iphysiq
536
537
538	c
539	c-----------------------------------------------------------------------
540	c Integration temporelle du modele :
541	c ----------------------------------
542
543	c write(78,*) 'ucov',ucov
544	c write(78,*) 'vcov',vcov
545	c write(78,*) 'teta',teta
546	c write(78,*) 'ps',ps
547	c write(78,*) 'q',q
548
549
550	CALL leapfrog(ucov,vcov,teta,ps,masse,phis,q,
551	. time_0)
552
553	END
554

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

Download in other formats:

Original Format