Context Navigation

gcm.F @ 1189

Last change on this file since 1189 was 1143, checked in by jghattas, 16 years ago

Recuperation des developpements fait uniquement sur la branche LMDZ4_V4_patches :

ajoute de la nouvelle flag ok_dynzon
ajoute du parametre aer_type
optimisation : isccp_cloud_types.F

+ bug pour le slab dans conf_phys.F90

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

File size: 13.4 KB

Line
1	!
2	! $Header$
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
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	c REAL rdayvrai,rdaym_ini,rday_ecri
116	c LOGICAL first
117
118	LOGICAL call_iniphys
119	data call_iniphys/.true./
120
121	c 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	c 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	c 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
138	C Calendrier
139	LOGICAL true_calendar
140	PARAMETER (true_calendar = .false.)
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	INTEGER :: ierr
153
154
155	c-----------------------------------------------------------------------
156	c Initialisations:
157	c ----------------
158
159	abort_message = 'last timestep reached'
160	modname = 'gcm'
161	descript = 'Run GCM LMDZ'
162	lafin = .FALSE.
163	dynhist_file = 'dyn_hist'
164	dynhistave_file = 'dyn_hist_ave'
165
166
167	c-----------------------------------------------------------------------
168	c Choix du calendrier
169	c -------------------
170
171	#ifdef CPP_IOIPSL
172	if (true_calendar) then
173	call ioconf_calendar('gregorian')
174	else
175	call ioconf_calendar('360d')
176	endif
177	#endif
178	c----------------------------------------------------------------------
179	c lecture des fichiers gcm.def ou run.def
180	c ---------------------------------------
181	c
182	! Ehouarn: dump possibility of using defrun
183	!#ifdef CPP_IOIPSL
184	CALL conf_gcm( 99, .TRUE. , clesphy0 )
185	!#else
186	! CALL defrun( 99, .TRUE. , clesphy0 )
187	!#endif
188	c
189	c
190	c------------------------------------
191	c Initialisation partie parallele
192	c------------------------------------
193	CALL init_const_mpi
194
195	call init_parallel
196	call Read_Distrib
197	! Ehouarn : temporarily (?) keep this only for Earth
198	if (planet_type.eq."earth") then
199	#ifdef CPP_EARTH
200	CALL Init_Phys_lmdz(iim,jjp1,llm,mpi_size,distrib_phys)
201	#endif
202	endif ! of if (planet_type.eq."earth")
203	CALL set_bands
204	CALL Init_interface_dyn_phys
205	CALL barrier
206
207	if (mpi_rank==0) call WriteBands
208	call SetDistrib(jj_Nb_Caldyn)
209
210	c$OMP PARALLEL
211	call Init_Mod_hallo
212	c$OMP END PARALLEL
213
214	! Ehouarn : temporarily (?) keep this only for Earth
215	if (planet_type.eq."earth") then
216	#ifdef CPP_EARTH
217	c$OMP PARALLEL
218	call InitComgeomphy
219	c$OMP END PARALLEL
220	#endif
221	endif ! of if (planet_type.eq."earth")
222
223	IF (config_inca /= 'none') THEN
224	#ifdef INCA
225	call init_const_lmdz(
226	$ nbtr,anneeref,dayref,
227	$ iphysiq,day_step,nday)
228
229	call init_inca_para(
230	$ iim,jjm+1,llm,klon_glo,mpi_size,
231	$ distrib_phys,COMM_LMDZ)
232	#endif
233	END IF
234
235	c-----------------------------------------------------------------------
236	c Initialisation des traceurs
237	c ---------------------------
238	c Choix du nombre de traceurs et du schema pour l'advection
239	c dans fichier traceur.def, par default ou via INCA
240	call infotrac_init
241
242	c Allocation de la tableau q : champs advectes
243	ALLOCATE(q(ip1jmp1,llm,nqtot))
244
245	c-----------------------------------------------------------------------
246	c Lecture de l'etat initial :
247	c ---------------------------
248
249	c lecture du fichier start.nc
250	if (read_start) then
251	! we still need to run iniacademic to initialize some
252	! constants & fields, if we run the 'newtonian' case:
253	if (iflag_phys.eq.2) then
254	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
255	endif
256	!#ifdef CPP_IOIPSL
257	if (planet_type.eq."earth") then
258	#ifdef CPP_EARTH
259	! Load an Earth-format start file
260	CALL dynetat0("start.nc",vcov,ucov,
261	. teta,q,masse,ps,phis, time_0)
262	#endif
263	endif ! of if (planet_type.eq."earth")
264	c write(73,*) 'ucov',ucov
265	c write(74,*) 'vcov',vcov
266	c write(75,*) 'teta',teta
267	c write(76,*) 'ps',ps
268	c write(77,*) 'q',q
269
270	endif ! of if (read_start)
271
272	c le cas echeant, creation d un etat initial
273	IF (prt_level > 9) WRITE(lunout,*)
274	. 'GCM: AVANT iniacademic AVANT AVANT AVANT AVANT'
275	if (.not.read_start) then
276	CALL iniacademic(vcov,ucov,teta,q,masse,ps,phis,time_0)
277	endif
278
279	c-----------------------------------------------------------------------
280	c Lecture des parametres de controle pour la simulation :
281	c -------------------------------------------------------
282	c on recalcule eventuellement le pas de temps
283
284	IF(MOD(day_step,iperiod).NE.0) THEN
285	abort_message =
286	. 'Il faut choisir un nb de pas par jour multiple de iperiod'
287	call abort_gcm(modname,abort_message,1)
288	ENDIF
289
290	IF(MOD(day_step,iphysiq).NE.0) THEN
291	abort_message =
292	* 'Il faut choisir un nb de pas par jour multiple de iphysiq'
293	call abort_gcm(modname,abort_message,1)
294	ENDIF
295
296	zdtvr = daysec/FLOAT(day_step)
297	IF(dtvr.NE.zdtvr) THEN
298	WRITE(lunout,*)
299	. 'WARNING!!! changement de pas de temps',dtvr,'>',zdtvr
300	ENDIF
301
302	C
303	C on remet le calendrier à zero si demande
304	c
305	if (annee_ref .ne. anneeref .or. day_ref .ne. dayref) then
306	write(lunout,*)
307	. ' Attention les dates initiales lues dans le fichier'
308	write(lunout,*)
309	. ' restart ne correspondent pas a celles lues dans '
310	write(lunout,*)' gcm.def'
311	if (raz_date .ne. 1) then
312	write(lunout,*)
313	. ' On garde les dates du fichier restart'
314	else
315	annee_ref = anneeref
316	day_ref = dayref
317	day_ini = dayref
318	itau_dyn = 0
319	itau_phy = 0
320	time_0 = 0.
321	write(lunout,*)
322	. ' On reinitialise a la date lue dans gcm.def'
323	endif
324	ELSE
325	raz_date = 0
326	endif
327
328
329	c nombre d'etats dans les fichiers demarrage et histoire
330	nbetatdem = nday / iecri
331	nbetatmoy = nday / periodav + 1
332
333	c-----------------------------------------------------------------------
334	c Initialisation des constantes dynamiques :
335	c ------------------------------------------
336	dtvr = zdtvr
337	CALL iniconst
338
339	c-----------------------------------------------------------------------
340	c Initialisation de la geometrie :
341	c --------------------------------
342	CALL inigeom
343
344	c-----------------------------------------------------------------------
345	c Initialisation du filtre :
346	c --------------------------
347	CALL inifilr
348	c
349	c-----------------------------------------------------------------------
350	c Initialisation de la dissipation :
351	c ----------------------------------
352
353	CALL inidissip( lstardis, nitergdiv, nitergrot, niterh ,
354	* tetagdiv, tetagrot , tetatemp )
355
356	c-----------------------------------------------------------------------
357	c Initialisation de la physique :
358	c -------------------------------
359	IF (call_iniphys.and.iflag_phys.eq.1) THEN
360	latfi(1)=rlatu(1)
361	lonfi(1)=0.
362	zcufi(1) = cu(1)
363	zcvfi(1) = cv(1)
364	DO j=2,jjm
365	DO i=1,iim
366	latfi((j-2)*iim+1+i)= rlatu(j)
367	lonfi((j-2)*iim+1+i)= rlonv(i)
368	zcufi((j-2)iim+1+i) = cu((j-1)iip1+i)
369	zcvfi((j-2)iim+1+i) = cv((j-1)iip1+i)
370	ENDDO
371	ENDDO
372	latfi(ngridmx)= rlatu(jjp1)
373	lonfi(ngridmx)= 0.
374	zcufi(ngridmx) = cu(ip1jm+1)
375	zcvfi(ngridmx) = cv(ip1jm-iim)
376	CALL gr_dyn_fi(1,iip1,jjp1,ngridmx,aire,airefi)
377
378	WRITE(lunout,*)
379	. 'GCM: WARNING!!! vitesse verticale nulle dans la physique'
380	! Earth:
381	if (planet_type.eq."earth") then
382	#ifdef CPP_EARTH
383	CALL iniphysiq(ngridmx,llm,daysec,day_ini,dtphys ,
384	, latfi,lonfi,airefi,zcufi,zcvfi,rad,g,r,cpp )
385	#endif
386	endif ! of if (planet_type.eq."earth")
387	call_iniphys=.false.
388	ENDIF ! of IF (call_iniphys.and.(iflag_phys.eq.1))
389
390
391	c-----------------------------------------------------------------------
392	c Initialisation des dimensions d'INCA :
393	c --------------------------------------
394	IF (config_inca /= 'none') THEN
395	!$OMP PARALLEL
396	#ifdef INCA
397	CALL init_inca_dim(klon_omp,llm,iim,jjm,
398	$ rlonu,rlatu,rlonv,rlatv)
399	#endif
400	!$OMP END PARALLEL
401	END IF
402
403	c-----------------------------------------------------------------------
404	c Initialisation des I/O :
405	c ------------------------
406
407
408	day_end = day_ini + nday
409	WRITE(lunout,300)day_ini,day_end
410	300 FORMAT('1'/,15x,'run du jour',i7,2x,'au jour',i7//)
411
412	!#ifdef CPP_IOIPSL
413	if (planet_type.eq."earth") then
414	#ifdef CPP_EARTH
415	CALL dynredem0_p("restart.nc", day_end, phis)
416	#endif
417	endif
418
419	ecripar = .TRUE.
420
421	#ifdef CPP_IOIPSL
422	if ( 1.eq.1) then
423	time_step = zdtvr
424	t_ops = iecri * daysec
425	t_wrt = iecri * daysec
426	CALL inithist_p(dynhist_file,day_ref,annee_ref,time_step,
427	. t_ops, t_wrt, histid, histvid)
428
429	IF (ok_dynzon) THEN
430	t_ops = iperiod * time_step
431	t_wrt = periodav * daysec
432	CALL initdynav_p(dynhistave_file,day_ref,annee_ref,time_step,
433	. t_ops, t_wrt, histaveid)
434	END IF
435	dtav = iperiod*dtvr/daysec
436	endif
437
438
439	#endif
440	! #endif of #ifdef CPP_IOIPSL
441
442	c Choix des frequences de stokage pour le offline
443	c istdyn=day_step/4 ! stockage toutes les 6h=1jour/4
444	c istdyn=day_step/12 ! stockage toutes les 2h=1jour/12
445	istdyn=day_step/4 ! stockage toutes les 6h=1jour/12
446	istphy=istdyn/iphysiq
447
448
449	c
450	c-----------------------------------------------------------------------
451	c Integration temporelle du modele :
452	c ----------------------------------
453
454	c write(78,*) 'ucov',ucov
455	c write(78,*) 'vcov',vcov
456	c write(78,*) 'teta',teta
457	c write(78,*) 'ps',ps
458	c write(78,*) 'q',q
459
460	c$OMP PARALLEL DEFAULT(SHARED) COPYIN(/temps/,/logic/)
461	CALL leapfrog_p(ucov,vcov,teta,ps,masse,phis,q,clesphy0,
462	. time_0)
463	c$OMP END PARALLEL
464
465
466	END
467

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

Download in other formats:

Original Format