Context Navigation

gcm.F @ 1565

Last change on this file since 1565 was 1563, checked in by jghattas, 13 years ago

Added reading of paramter type_trac from *.def file : type_trac=lmdz(default), inca or repr(soon). While running with INCA, 2 parameters are now necessare in .def : type_trac=inca and config_inca=aero/chem. If type_trac=lmdz or repr, config_inca will not be used.

Removed print of ecrit_mth, ecrit_day etc in physiq.F. Removed the variable ecrit_hf2mth which is no longer used.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 16.7 KB

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

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

Download in other formats:

Original Format