Context Navigation

gcm.F @ 1357

Last change on this file since 1357 was 1357, checked in by Ehouarn Millour, 14 years ago

Some cleanup and fixing the possibility to output fields in the dynamics, on the dynamical grids.

CLEANUPS:

arch-PW6_VARGAS.fcm : add potentially benefic compiling options
removed obsolete "control.h" in dyn3d/dyn3dpar (module control_mod.F90 is used instead)

OUTPUTS in the dynamics (3 sets of files, one for each grid: scalar, u, v):

removed "com_io_dyn.h" common; use module "com_io_dyn_mod.F90" instead
updated "initdynav.F","inithist.F","writehist.F" and "writedynav.F" in bibio: which field will be written is hard coded there.
flags "ok_dyn_ins" and "ok_dyn_ave" (loaded via conf_gcm.F) trigger output of fields in the dynamics: if ok_dyn_ins is true, then files "dyn_hist.nc", "dyn_histu.nc" and "dyn_histv.nc" are written (the frequency of the outputs is given by 'iecri' in run.def; values are written every 'iecri' dynamical step). if ok_dyn_ave is true then files "dyn_hist_ave.nc", "dyn_histu_ave.nc" and "dyn_histv_ave.nc" are written (the rate at which averages and made/written, in days, is given by 'periodav' in run.def).

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

File size: 16.6 KB

Line
1	!
2	! $Id: gcm.F 1357 2010-04-14 14:03:19Z 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
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	#ifdef INCA
76	! Only INCA needs these informations (from the Earth's physics)
77	#include "indicesol.h"
78	#endif
79
80	INTEGER longcles
81	PARAMETER ( longcles = 20 )
82	REAL clesphy0( longcles )
83	SAVE clesphy0
84
85
86
87	REAL zdtvr
88	c INTEGER nbetatmoy, nbetatdem,nbetat
89	INTEGER nbetatmoy, nbetatdem
90
91	c variables dynamiques
92	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
93	REAL teta(ip1jmp1,llm) ! temperature potentielle
94	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: q ! champs advectes
95	REAL ps(ip1jmp1) ! pression au sol
96	c REAL p (ip1jmp1,llmp1 ) ! pression aux interfac.des couches
97	c REAL pks(ip1jmp1) ! exner au sol
98	c REAL pk(ip1jmp1,llm) ! exner au milieu des couches
99	c REAL pkf(ip1jmp1,llm) ! exner filt.au milieu des couches
100	REAL masse(ip1jmp1,llm) ! masse d'air
101	REAL phis(ip1jmp1) ! geopotentiel au sol
102	c REAL phi(ip1jmp1,llm) ! geopotentiel
103	c REAL w(ip1jmp1,llm) ! vitesse verticale
104
105	c variables dynamiques intermediaire pour le transport
106
107	c variables pour le fichier histoire
108	REAL dtav ! intervalle de temps elementaire
109
110	REAL time_0
111
112	LOGICAL lafin
113	c INTEGER ij,iq,l,i,j
114	INTEGER i,j
115
116
117	real time_step, t_wrt, t_ops
118
119
120	LOGICAL call_iniphys
121	data call_iniphys/.true./
122
123	c REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
124	c+jld variables test conservation energie
125	c REAL ecin(ip1jmp1,llm),ecin0(ip1jmp1,llm)
126	C Tendance de la temp. potentiel d (theta)/ d t due a la
127	C tansformation d'energie cinetique en energie thermique
128	C cree par la dissipation
129	c REAL dhecdt(ip1jmp1,llm)
130	c REAL vcont(ip1jm,llm),ucont(ip1jmp1,llm)
131	c REAL d_h_vcol, d_qt, d_qw, d_ql, d_ec
132	c CHARACTER (len=15) :: ztit
133	c-jld
134
135
136	character (len=80) :: dynhist_file, dynhistave_file
137	character (len=20) :: modname
138	character (len=80) :: abort_message
139	! locales pour gestion du temps
140	INTEGER :: an, mois, jour
141	REAL :: heure
142
143
144	c-----------------------------------------------------------------------
145	c variables pour l'initialisation de la physique :
146	c ------------------------------------------------
147	INTEGER ngridmx
148	PARAMETER( ngridmx = 2+(jjm-1)*iim - 1/jjm )
149	REAL zcufi(ngridmx),zcvfi(ngridmx)
150	REAL latfi(ngridmx),lonfi(ngridmx)
151	REAL airefi(ngridmx)
152	SAVE latfi, lonfi, airefi
153
154	INTEGER :: ierr
155
156
157	c-----------------------------------------------------------------------
158	c Initialisations:
159	c ----------------
160
161	abort_message = 'last timestep reached'
162	modname = 'gcm'
163	descript = 'Run GCM LMDZ'
164	lafin = .FALSE.
165	dynhist_file = 'dyn_hist'
166	dynhistave_file = 'dyn_hist_ave'
167
168
169
170	c----------------------------------------------------------------------
171	c lecture des fichiers gcm.def ou run.def
172	c ---------------------------------------
173	c
174	! Ehouarn: dump possibility of using defrun
175	!#ifdef CPP_IOIPSL
176	CALL conf_gcm( 99, .TRUE. , clesphy0 )
177	!#else
178	! CALL defrun( 99, .TRUE. , clesphy0 )
179	!#endif
180	c
181	c
182	c------------------------------------
183	c Initialisation partie parallele
184	c------------------------------------
185	CALL init_const_mpi
186
187	call init_parallel
188	call ini_getparam("out.def")
189	call Read_Distrib
190	! Ehouarn : temporarily (?) keep this only for Earth
191	if (planet_type.eq."earth") then
192	#ifdef CPP_EARTH
193	CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
194	#endif
195	endif ! of if (planet_type.eq."earth")
196	CALL set_bands
197	#ifdef CPP_EARTH
198	! Ehouarn: For now only Earth physics is parallel
199	CALL Init_interface_dyn_phys
200	#endif
201	CALL barrier
202
203	if (mpi_rank==0) call WriteBands
204	call SetDistrib(jj_Nb_Caldyn)
205
206	c$OMP PARALLEL
207	call Init_Mod_hallo
208	c$OMP END PARALLEL
209
210	! Ehouarn : temporarily (?) keep this only for Earth
211	if (planet_type.eq."earth") then
212	#ifdef CPP_EARTH
213	c$OMP PARALLEL
214	call InitComgeomphy
215	c$OMP END PARALLEL
216	#endif
217	endif ! of if (planet_type.eq."earth")
218
219	c-----------------------------------------------------------------------
220	c Choix du calendrier
221	c -------------------
222
223	c calend = 'earth_365d'
224
225	#ifdef CPP_IOIPSL
226	if (calend == 'earth_360d') then
227	call ioconf_calendar('360d')
228	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 360 jours/an'
229	else if (calend == 'earth_365d') then
230	call ioconf_calendar('noleap')
231	write(lunout,*)'CALENDRIER CHOISI: Terrestre a 365 jours/an'
232	else if (calend == 'earth_366d') then
233	call ioconf_calendar('gregorian')
234	write(lunout,*)'CALENDRIER CHOISI: Terrestre bissextile'
235	else
236	abort_message = 'Mauvais choix de calendrier'
237	call abort_gcm(modname,abort_message,1)
238	endif
239	#endif
240
241	IF (config_inca /= 'none') THEN
242	#ifdef INCA
243	call init_const_lmdz(
244	$ nbtr,anneeref,dayref,
245	$ iphysiq,day_step,nday,
246	$ nbsrf, is_oce,is_sic,
247	$ is_ter,is_lic)
248
249	call init_inca_para(
250	$ iim,jjm+1,llm,klon_glo,mpi_size,
251	$ distrib_phys,COMM_LMDZ)
252	#endif
253	END IF
254
255	c-----------------------------------------------------------------------
256	c Initialisation des traceurs
257	c ---------------------------
258	c Choix du nombre de traceurs et du schema pour l'advection
259	c dans fichier traceur.def, par default ou via INCA
260	call infotrac_init
261
262	c Allocation de la tableau q : champs advectes
263	ALLOCATE(q(ip1jmp1,llm,nqtot))
264
265	c-----------------------------------------------------------------------
266	c Lecture de l'etat initial :
267	c ---------------------------
268
269	c lecture du fichier start.nc
270	if (read_start) then
271	! we still need to run iniacademic to initialize some
272	! constants & fields, if we run the 'newtonian' case:
273	if (iflag_phys.eq.2) then
274	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
275	endif
276	!#ifdef CPP_IOIPSL
277	if (planet_type.eq."earth") then
278	#ifdef CPP_EARTH
279	! Load an Earth-format start file
280	CALL dynetat0("start.nc",vcov,ucov,
281	. teta,q,masse,ps,phis, time_0)
282	#endif
283	endif ! of if (planet_type.eq."earth")
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	c-----------------------------------------------------------------------
397	c Initialisation des constantes dynamiques :
398	c ------------------------------------------
399	dtvr = zdtvr
400	CALL iniconst
401
402	c-----------------------------------------------------------------------
403	c Initialisation de la geometrie :
404	c --------------------------------
405	CALL inigeom
406
407	c-----------------------------------------------------------------------
408	c Initialisation du filtre :
409	c --------------------------
410	CALL inifilr
411	c
412	c-----------------------------------------------------------------------
413	c Initialisation de la dissipation :
414	c ----------------------------------
415
416	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
417	* tetagdiv, tetagrot , tetatemp )
418
419	c-----------------------------------------------------------------------
420	c Initialisation de la physique :
421	c -------------------------------
422	IF (call_iniphys.and.iflag_phys.eq.1) THEN
423	latfi(1)=rlatu(1)
424	lonfi(1)=0.
425	zcufi(1) = cu(1)
426	zcvfi(1) = cv(1)
427	DO j=2,jjm
428	DO i=1,iim
429	latfi((j-2)*iim+1+i)= rlatu(j)
430	lonfi((j-2)*iim+1+i)= rlonv(i)
431	zcufi((j-2)iim+1+i) = cu((j-1)iip1+i)
432	zcvfi((j-2)iim+1+i) = cv((j-1)iip1+i)
433	ENDDO
434	ENDDO
435	latfi(ngridmx)= rlatu(jjp1)
436	lonfi(ngridmx)= 0.
437	zcufi(ngridmx) = cu(ip1jm+1)
438	zcvfi(ngridmx) = cv(ip1jm-iim)
439	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
440
441	WRITE(lunout,*)
442	. 'GCM: WARNING!!! vitesse verticale nulle dans la physique'
443	! Earth:
444	if (planet_type.eq."earth") then
445	#ifdef CPP_EARTH
446	CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys/nsplit_phys ,
447	, latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp )
448	#endif
449	endif ! of if (planet_type.eq."earth")
450	call_iniphys=.false.
451	ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
452
453
454	c-----------------------------------------------------------------------
455	c Initialisation des dimensions d'INCA :
456	c --------------------------------------
457	IF (config_inca /= 'none') THEN
458	!$OMP PARALLEL
459	#ifdef INCA
460	CALL init_inca_dim(klon_omp,llm,iim,jjm,
461	$ rlonu,rlatu,rlonv,rlatv)
462	#endif
463	!$OMP END PARALLEL
464	END IF
465
466	c-----------------------------------------------------------------------
467	c Initialisation des I/O :
468	c ------------------------
469
470
471	day_end = day_ini + nday
472	WRITE(lunout,300)day_ini,day_end
473	300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
474
475	#ifdef CPP_IOIPSL
476	call ju2ymds(jD_ref + day_ini - day_ref, an, mois, jour, heure)
477	write (lunout,301)jour, mois, an
478	call ju2ymds(jD_ref + day_end - day_ref, an, mois, jour, heure)
479	write (lunout,302)jour, mois, an
480	301 FORMAT('1'/,15x,'run du ', i2,'/',i2,'/',i4)
481	302 FORMAT('1'/,15x,' au ', i2,'/',i2,'/',i4)
482	#endif
483
484	if (planet_type.eq."earth") then
485	CALL dynredem0_p("restart.nc", day_end, phis)
486	endif
487
488	ecripar = .TRUE.
489
490	#ifdef CPP_IOIPSL
491	time_step = zdtvr
492	if (ok_dyn_ins) then
493	! initialize output file for instantaneous outputs
494	! t_ops = iecri * daysec ! do operations every t_ops
495	t_ops =((1.0iecri)/day_step) daysec
496	t_wrt = daysec ! iecri * daysec ! write output every t_wrt
497	t_wrt = daysec ! iecri * daysec ! write output every t_wrt
498	CALL inithist(day_ref,annee_ref,time_step,
499	& t_ops,t_wrt)
500	! CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
501	! . t_ops, t_wrt, histid, histvid)
502	endif
503
504	IF (ok_dyn_ave) THEN
505	! initialize output file for averaged outputs
506	t_ops = iperiod * time_step ! do operations every t_ops
507	t_wrt = periodav * daysec ! write output every t_wrt
508	CALL initdynav(day_ref,annee_ref,time_step,
509	& t_ops,t_wrt)
510	! CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
511	! . t_ops, t_wrt, histaveid)
512	END IF
513	dtav = iperiod*dtvr/daysec
514	#endif
515	! #endif of #ifdef CPP_IOIPSL
516
517	c Choix des frequences de stokage pour le offline
518	c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
519	c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12
520	istdyn=day_step/4 ! stockage toutes les 6h=1jour/12
521	istphy=istdyn/iphysiq
522
523
524	c
525	c-----------------------------------------------------------------------
526	c Integration temporelle du modele :
527	c ----------------------------------
528
529	c write(78,*) 'ucov',ucov
530	c write(78,*) 'vcov',vcov
531	c write(78,*) 'teta',teta
532	c write(78,*) 'ps',ps
533	c write(78,*) 'q',q
534
535	c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
536	CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
537	. time_0)
538	c$OMP END PARALLEL
539
540
541	END
542

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

Download in other formats:

Original Format