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

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

Download in other formats:

Original Format