Context Navigation

guide_mod.F90 @ 1508

Last change on this file since 1508 was 1508, checked in by emillour, 9 years ago

Common dynamics:
Updates in the dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2325):
IMPORTANT: Modifications for isotopes are only done in dyn3d, not in dyn3dpar

as in LMDZ5 these modifications were done in dyn3dmem.
Related LMDZ5 revisions are r2270 and r2281

in dynlonlat_phylonlat:
add module "grid_atob_m.F90" (a regridding utility so far only used by phylmd/ce0l.F90, used to be dyn3d_common/grid_atob.F)

in misc:
follow up updates on wxios.F (add missing_val module variable)

in dyn3d_common:
pression.F => pression.F90
misc_mod.F90: moved from misc to dyn3d_common
added new iso_verif_dyn.F
covcont.F => covcont.F90
infotrac.F90 : add handling of isotopes (reading of corresponding traceur.def for planets not implemented)
dynetat0.F => dynetat0.F90 with some code factorization
dynredem.F => dynredem.F90 with some code factorization
added dynredem_mod.F90: routines used by dynredem
iniacademic.F90 : added isotopes-related initialization for Earth case

in dyn3d:
added check_isotopes.F
modified (isotopes) advtrac.F90, caladvtrac.F
guide_mod.F90: ported updates
leapfrog.F : (isotopes) updates (NB: call integrd with nqtot tracers)
qminimium.F : adaptations for isotopes (copied over, except that #include comvert.h is not needed).
vlsplt.F: adaptations for isotopes (copied over, except than #include logic.h, comvert.h not needed, and replace "include comconst.h" with use comconst_mod, ONLY: pi)
vlspltqs.F : same as vlsplt.F, but also keeping added modification for CP(T)

in dyn3dpar:
leapfrog_p.F: remove unecessary #ifdef CPP_EARTH cpp flag. and call integrd_p with nqtot tracers (only important for Earth)
dynredem_p.F => dynredem_p.F90 and some code factorization
and no isotopes-relates changes in dyn3dpar (since these changes have been made in LMDZ5 dyn3dmem).

File size: 58.5 KB

Line
1	!
2	! $Id$
3	!
4	MODULE guide_mod
5
6	!=======================================================================
7	! Auteur: F.Hourdin
8	! F. Codron 01/09
9	!=======================================================================
10
11	USE getparam
12	USE Write_Field
13	use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
14	use pres2lev_mod
15	USE serre_mod, ONLY: grossismx
16
17	IMPLICIT NONE
18
19	! ---------------------------------------------
20	! Declarations des cles logiques et parametres
21	! ---------------------------------------------
22	INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
23	INTEGER, PRIVATE, SAVE :: nlevnc
24	LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P
25	LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta
26	LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
27	LOGICAL, PRIVATE, SAVE :: guide_modele,invert_p,invert_y,ini_anal
28	LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav
29
30	REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
31	REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
32	REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
33	REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
34	REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
35
36	REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
37	REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
38	REAL, PRIVATE, SAVE :: tau_lon,tau_lat
39
40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
41	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q
42	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
43
44	! ---------------------------------------------
45	! Variables de guidage
46	! ---------------------------------------------
47	! Variables des fichiers de guidage
48	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
49	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
50	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
51	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
52	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2
53	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc
54	! Variables aux dimensions du modele
55	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2
56	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2
57	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2
58	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2
59	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2
60
61	CONTAINS
62	!=======================================================================
63
64	SUBROUTINE guide_init
65
66	USE control_mod
67
68	IMPLICIT NONE
69
70	INCLUDE "dimensions.h"
71	INCLUDE "paramet.h"
72	INCLUDE "netcdf.inc"
73
74	INTEGER :: error,ncidpl,rid,rcod
75	CHARACTER (len = 80) :: abort_message
76	CHARACTER (len = 20) :: modname = 'guide_init'
77
78	! ---------------------------------------------
79	! Lecture des parametres:
80	! ---------------------------------------------
81	call ini_getparam("nudging_parameters_out.txt")
82	! Variables guidees
83	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
84	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
85	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
86	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
87	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
88	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
89	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
90
91	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
92	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
93	if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
94	call abort_gcm("guide_init", &
95	"zonal nudging requires grid regular in longitude", 1)
96
97	! Constantes de rappel. Unite : fraction de jour
98	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
99	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
100	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
101	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
102	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
103	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
104	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
105	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
106	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
107	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
108	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
109	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
110
111	! Sauvegarde du for�age
112	CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
113	CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
114	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
115	IF (iguide_sav.GT.0) THEN
116	iguide_sav=day_step/iguide_sav
117	ELSE if (iguide_sav == 0) then
118	iguide_sav = huge(0)
119	ELSE
120	iguide_sav=day_step*iguide_sav
121	ENDIF
122
123	! Guidage regional seulement (sinon constant ou suivant le zoom)
124	CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
125	CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
126	CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
127	CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
128	CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
129	CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
130	CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
131
132	! Parametres pour lecture des fichiers
133	CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
134	CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
135	IF (iguide_int.EQ.0) THEN
136	iguide_int=1
137	ELSEIF (iguide_int.GT.0) THEN
138	iguide_int=day_step/iguide_int
139	ELSE
140	iguide_int=day_step*iguide_int
141	ENDIF
142	CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele')
143	CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
144	CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
145	CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
146	CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
147
148	call fin_getparam
149
150	! ---------------------------------------------
151	! Determination du nombre de niveaux verticaux
152	! des fichiers guidage
153	! ---------------------------------------------
154	ncidpl=-99
155	if (guide_modele) then
156	if (ncidpl.eq.-99) then
157	rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
158	if (rcod.NE.NF_NOERR) THEN
159	CALL abort_gcm(modname, &
160	'Guide: probleme -> pas de fichier apbp.nc',1)
161	endif
162	endif
163	else
164	if (guide_u) then
165	if (ncidpl.eq.-99) then
166	rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
167	if (rcod.NE.NF_NOERR) THEN
168	CALL abort_gcm(modname, &
169	'Guide: probleme -> pas de fichier u.nc',1)
170	endif
171	endif
172	elseif (guide_v) then
173	if (ncidpl.eq.-99) then
174	rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
175	if (rcod.NE.NF_NOERR) THEN
176	CALL abort_gcm(modname, &
177	'Guide: probleme -> pas de fichier v.nc',1)
178	endif
179	endif
180	elseif (guide_T) then
181	if (ncidpl.eq.-99) then
182	rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
183	if (rcod.NE.NF_NOERR) THEN
184	CALL abort_gcm(modname, &
185	'Guide: probleme -> pas de fichier T.nc',1)
186	endif
187	endif
188	elseif (guide_Q) then
189	if (ncidpl.eq.-99) then
190	rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
191	if (rcod.NE.NF_NOERR) THEN
192	CALL abort_gcm(modname, &
193	'Guide: probleme -> pas de fichier hur.nc',1)
194	endif
195	endif
196	endif
197	endif
198	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
199	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
200	IF (error.NE.NF_NOERR) THEN
201	CALL abort_gcm(modname,'Guide: probleme lecture niveaux pression',1)
202	ENDIF
203	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
204	print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
205	rcod = nf90_close(ncidpl)
206
207	! ---------------------------------------------
208	! Allocation des variables
209	! ---------------------------------------------
210	abort_message='pb in allocation guide'
211
212	ALLOCATE(apnc(nlevnc), stat = error)
213	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
214	ALLOCATE(bpnc(nlevnc), stat = error)
215	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
216	apnc=0.;bpnc=0.
217
218	ALLOCATE(alpha_pcor(llm), stat = error)
219	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
220	ALLOCATE(alpha_u(ip1jmp1), stat = error)
221	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
222	ALLOCATE(alpha_v(ip1jm), stat = error)
223	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
224	ALLOCATE(alpha_T(ip1jmp1), stat = error)
225	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
226	ALLOCATE(alpha_Q(ip1jmp1), stat = error)
227	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
228	ALLOCATE(alpha_P(ip1jmp1), stat = error)
229	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
230	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
231
232	IF (guide_u) THEN
233	ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
234	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
235	ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
236	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237	ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
238	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239	ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
240	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
241	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
242	ENDIF
243
244	IF (guide_T) THEN
245	ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
246	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
247	ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
248	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
249	ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
250	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
251	ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
252	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
253	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
254	ENDIF
255
256	IF (guide_Q) THEN
257	ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
258	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
259	ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
260	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
261	ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
262	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
263	ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
264	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
265	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
266	ENDIF
267
268	IF (guide_v) THEN
269	ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
270	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
271	ALLOCATE(vgui1(ip1jm,llm), stat = error)
272	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
273	ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
274	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
275	ALLOCATE(vgui2(ip1jm,llm), stat = error)
276	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
277	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
278	ENDIF
279
280	IF (guide_P.OR.guide_modele) THEN
281	ALLOCATE(psnat1(iip1,jjp1), stat = error)
282	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
283	ALLOCATE(psnat2(iip1,jjp1), stat = error)
284	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
285	psnat1=0.;psnat2=0.;
286	ENDIF
287	IF (guide_P) THEN
288	ALLOCATE(psgui2(ip1jmp1), stat = error)
289	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
290	ALLOCATE(psgui1(ip1jmp1), stat = error)
291	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
292	psgui1=0.;psgui2=0.
293	ENDIF
294
295	! ---------------------------------------------
296	! Lecture du premier etat de guidage.
297	! ---------------------------------------------
298	IF (guide_2D) THEN
299	CALL guide_read2D(1)
300	ELSE
301	CALL guide_read(1)
302	ENDIF
303	IF (guide_v) vnat1=vnat2
304	IF (guide_u) unat1=unat2
305	IF (guide_T) tnat1=tnat2
306	IF (guide_Q) qnat1=qnat2
307	IF (guide_P.OR.guide_modele) psnat1=psnat2
308
309	END SUBROUTINE guide_init
310
311	!=======================================================================
312	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
313
314	USE control_mod
315	USE comvert_mod, ONLY: ap,bp,preff,presnivs
316	USE comconst_mod, ONLY: daysec,dtvr
317
318	IMPLICIT NONE
319
320	INCLUDE "dimensions.h"
321	INCLUDE "paramet.h"
322
323	! Variables entree
324	INTEGER, INTENT(IN) :: itau !pas de temps
325	REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
326	REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
327	REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
328
329	! Variables locales
330	LOGICAL, SAVE :: first=.TRUE.
331	LOGICAL :: f_out ! sortie guidage
332	REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
333	REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P
334	! Compteurs temps:
335	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
336	REAL :: ditau, dday_step
337	REAL :: tau,reste ! position entre 2 etats de guidage
338	REAL, SAVE :: factt ! pas de temps en fraction de jour
339
340	INTEGER :: l
341
342	!-----------------------------------------------------------------------
343	! Initialisations au premier passage
344	!-----------------------------------------------------------------------
345	IF (first) THEN
346	first=.FALSE.
347	CALL guide_init
348	itau_test=1001
349	step_rea=1
350	count_no_rea=0
351	! Calcul des constantes de rappel
352	factt=dtvr*iperiod/daysec
353	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
354	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
355	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
356	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
357	call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
358	! correction de rappel dans couche limite
359	if (guide_BL) then
360	alpha_pcor(:)=1.
361	else
362	do l=1,llm
363	alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
364	enddo
365	endif
366	! ini_anal: etat initial egal au guidage
367	IF (ini_anal) THEN
368	CALL guide_interp(ps,teta)
369	IF (guide_u) ucov=ugui2
370	IF (guide_v) vcov=ugui2
371	IF (guide_T) teta=tgui2
372	IF (guide_Q) q=qgui2
373	IF (guide_P) THEN
374	ps=psgui2
375	CALL pression(ip1jmp1,ap,bp,ps,p)
376	CALL massdair(p,masse)
377	ENDIF
378	RETURN
379	ENDIF
380	! Verification structure guidage
381	IF (guide_u) THEN
382	CALL writefield('unat',unat1)
383	CALL writefield('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
384	ENDIF
385	IF (guide_T) THEN
386	CALL writefield('tnat',tnat1)
387	CALL writefield('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
388	ENDIF
389
390	ENDIF !first
391
392	!-----------------------------------------------------------------------
393	! Lecture des fichiers de guidage ?
394	!-----------------------------------------------------------------------
395	IF (iguide_read.NE.0) THEN
396	ditau=real(itau)
397	dday_step=real(day_step)
398	IF (iguide_read.LT.0) THEN
399	tau=ditau/dday_step/REAL(iguide_read)
400	ELSE
401	tau=REAL(iguide_read)*ditau/dday_step
402	ENDIF
403	reste=tau-AINT(tau)
404	IF (reste.EQ.0.) THEN
405	IF (itau_test.EQ.itau) THEN
406	write(,)'deuxieme passage de advreel a itau=',itau
407	stop
408	ELSE
409	IF (guide_v) vnat1=vnat2
410	IF (guide_u) unat1=unat2
411	IF (guide_T) tnat1=tnat2
412	IF (guide_Q) qnat1=qnat2
413	IF (guide_P.OR.guide_modele) psnat1=psnat2
414	step_rea=step_rea+1
415	itau_test=itau
416	print*,'Lecture fichiers guidage, pas ',step_rea, &
417	'apres ',count_no_rea,' non lectures'
418	IF (guide_2D) THEN
419	CALL guide_read2D(step_rea)
420	ELSE
421	CALL guide_read(step_rea)
422	ENDIF
423	count_no_rea=0
424	ENDIF
425	ELSE
426	count_no_rea=count_no_rea+1
427
428	ENDIF
429	ENDIF !iguide_read=0
430
431	!-----------------------------------------------------------------------
432	! Interpolation et conversion des champs de guidage
433	!-----------------------------------------------------------------------
434	IF (MOD(itau,iguide_int).EQ.0) THEN
435	CALL guide_interp(ps,teta)
436	ENDIF
437	! Repartition entre 2 etats de guidage
438	IF (iguide_read.NE.0) THEN
439	tau=reste
440	ELSE
441	tau=1.
442	ENDIF
443
444	!-----------------------------------------------------------------------
445	! Ajout des champs de guidage
446	!-----------------------------------------------------------------------
447	! Sauvegarde du guidage?
448	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
449	IF (f_out) CALL guide_out("SP",jjp1,1,ps)
450
451	if (guide_u) then
452	if (guide_add) then
453	f_add=(1.-tau)ugui1+tauugui2
454	else
455	f_add=(1.-tau)ugui1+tauugui2-ucov
456	endif
457	if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
458	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
459	IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)ugui1+tauugui2)
460	IF (f_out) CALL guide_out("u",jjp1,llm,ucov)
461	IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add/factt)
462	ucov=ucov+f_add
463	endif
464
465	if (guide_T) then
466	if (guide_add) then
467	f_add=(1.-tau)tgui1+tautgui2
468	else
469	f_add=(1.-tau)tgui1+tautgui2-teta
470	endif
471	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
472	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
473	IF (f_out) CALL guide_out("teta",jjp1,llm,f_add/factt)
474	teta=teta+f_add
475	endif
476
477	if (guide_P) then
478	if (guide_add) then
479	f_add(1:ip1jmp1,1)=(1.-tau)psgui1+taupsgui2
480	else
481	f_add(1:ip1jmp1,1)=(1.-tau)psgui1+taupsgui2-ps
482	endif
483	if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
484	CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
485	IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt)
486	ps=ps+f_add(1:ip1jmp1,1)
487	CALL pression(ip1jmp1,ap,bp,ps,p)
488	CALL massdair(p,masse)
489	endif
490
491	if (guide_Q) then
492	if (guide_add) then
493	f_add=(1.-tau)qgui1+tauqgui2
494	else
495	f_add=(1.-tau)qgui1+tauqgui2-q
496	endif
497	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
498	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
499	IF (f_out) CALL guide_out("q",jjp1,llm,f_add/factt)
500	q=q+f_add
501	endif
502
503	if (guide_v) then
504	if (guide_add) then
505	f_add(1:ip1jm,:)=(1.-tau)vgui1+tauvgui2
506	else
507	f_add(1:ip1jm,:)=(1.-tau)vgui1+tauvgui2-vcov
508	endif
509	if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
510	CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
511	IF (f_out) CALL guide_out("v",jjm,llm,vcov)
512	IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)vgui1+tauvgui2)
513	IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt)
514	vcov=vcov+f_add(1:ip1jm,:)
515	endif
516	END SUBROUTINE guide_main
517
518	!=======================================================================
519	SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
520	! field1=afield1+alphafield2
521
522	IMPLICIT NONE
523
524	! input variables
525	INTEGER, INTENT(IN) :: hsize
526	INTEGER, INTENT(IN) :: vsize
527	REAL, DIMENSION(hsize), INTENT(IN) :: alpha
528	REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
529
530	! Local variables
531	INTEGER :: l
532
533	do l=1,vsize
534	field(:,l)=alphafield(:,l)alpha_pcor(l)
535	enddo
536
537	END SUBROUTINE guide_addfield
538
539	!=======================================================================
540	SUBROUTINE guide_zonave(typ,hsize,vsize,field)
541
542	USE comconst_mod, ONLY: pi
543
544	IMPLICIT NONE
545
546	INCLUDE "dimensions.h"
547	INCLUDE "paramet.h"
548	INCLUDE "comgeom.h"
549
550	! input/output variables
551	INTEGER, INTENT(IN) :: typ
552	INTEGER, INTENT(IN) :: vsize
553	INTEGER, INTENT(IN) :: hsize
554	REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
555
556	! Local variables
557	LOGICAL, SAVE :: first=.TRUE.
558	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
559	INTEGER :: i,j,l,ij
560	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
561	REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
562
563	IF (first) THEN
564	first=.FALSE.
565	!Compute domain for averaging
566	lond=rlonu*180./pi
567	imin(1)=1;imax(1)=iip1;
568	imin(2)=1;imax(2)=iip1;
569	IF (guide_reg) THEN
570	DO i=1,iim
571	IF (lond(i).LT.lon_min_g) imin(1)=i
572	IF (lond(i).LE.lon_max_g) imax(1)=i
573	ENDDO
574	lond=rlonv*180./pi
575	DO i=1,iim
576	IF (lond(i).LT.lon_min_g) imin(2)=i
577	IF (lond(i).LE.lon_max_g) imax(2)=i
578	ENDDO
579	ENDIF
580	ENDIF
581
582	fieldm=0.
583	DO l=1,vsize
584	! Compute zonal average
585	DO j=1,hsize
586	DO i=imin(typ),imax(typ)
587	ij=(j-1)*iip1+i
588	fieldm(j,l)=fieldm(j,l)+field(ij,l)
589	ENDDO
590	ENDDO
591	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
592	! Compute forcing
593	DO j=1,hsize
594	DO i=1,iip1
595	ij=(j-1)*iip1+i
596	field(ij,l)=fieldm(j,l)
597	ENDDO
598	ENDDO
599	ENDDO
600
601	END SUBROUTINE guide_zonave
602
603	!=======================================================================
604	SUBROUTINE guide_interp(psi,teta)
605
606	use exner_hyb_m, only: exner_hyb
607	use exner_milieu_m, only: exner_milieu
608	USE comvert_mod, ONLY: ap,bp,preff,pressure_exner
609	USE comconst_mod, ONLY: cpp,kappa
610
611	IMPLICIT NONE
612
613	include "dimensions.h"
614	include "paramet.h"
615	include "comgeom2.h"
616
617	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
618	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
619
620	LOGICAL, SAVE :: first=.TRUE.
621	! Variables pour niveaux pression:
622	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
623	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
624	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
625	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
626	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
627	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
628	REAL, DIMENSION (iip1,jjm,llm) :: pbary
629	! Variables pour fonction Exner (P milieu couche)
630	REAL, DIMENSION (iip1,jjp1,llm) :: pk
631	REAL, DIMENSION (iip1,jjp1) :: pks
632	REAL :: prefkap,unskap
633	! Pression de vapeur saturante
634	REAL, DIMENSION (ip1jmp1,llm) :: qsat
635	!Variables intermediaires interpolation
636	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
637	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
638
639	INTEGER :: i,j,l,ij
640
641	print *,'Guide: conversion variables guidage'
642	! -----------------------------------------------------------------
643	! Calcul des niveaux de pression champs guidage
644	! -----------------------------------------------------------------
645	if (guide_modele) then
646	do i=1,iip1
647	do j=1,jjp1
648	do l=1,nlevnc
649	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
650	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
651	enddo
652	enddo
653	enddo
654	else
655	do i=1,iip1
656	do j=1,jjp1
657	do l=1,nlevnc
658	plnc2(i,j,l)=apnc(l)
659	plnc1(i,j,l)=apnc(l)
660	enddo
661	enddo
662	enddo
663
664	endif
665	if (first) then
666	first=.FALSE.
667	print*,'Guide: verification ordre niveaux verticaux'
668	print*,'LMDZ :'
669	do l=1,llm
670	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
671	+psi(1,jjp1)*(bp(l)+bp(l+1))/2.
672	enddo
673	print*,'Fichiers guidage'
674	do l=1,nlevnc
675	print*,'PL(',l,')=',plnc2(1,1,l)
676	enddo
677	print *,'inversion de l''ordre: invert_p=',invert_p
678	if (guide_u) then
679	do l=1,nlevnc
680	print*,'U(',l,')=',unat2(1,1,l)
681	enddo
682	endif
683	if (guide_T) then
684	do l=1,nlevnc
685	print*,'T(',l,')=',tnat2(1,1,l)
686	enddo
687	endif
688	endif
689
690	! -----------------------------------------------------------------
691	! Calcul niveaux pression modele
692	! -----------------------------------------------------------------
693	CALL pression( ip1jmp1, ap, bp, psi, p )
694	if (pressure_exner) then
695	CALL exner_hyb(ip1jmp1,psi,p,pks,pk)
696	else
697	CALL exner_milieu(ip1jmp1,psi,p,pks,pk)
698	endif
699	! .... Calcul de pls , pression au milieu des couches ,en Pascals
700	unskap=1./kappa
701	prefkap = preff ** kappa
702	DO l = 1, llm
703	DO j=1,jjp1
704	DO i =1, iip1
705	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
706	ENDDO
707	ENDDO
708	ENDDO
709
710	! calcul des pressions pour les grilles u et v
711	do l=1,llm
712	do j=1,jjp1
713	do i=1,iip1
714	pext(i,j,l)=pls(i,j,l)*aire(i,j)
715	enddo
716	enddo
717	enddo
718	call massbar(pext, pbarx, pbary )
719	do l=1,llm
720	do j=1,jjp1
721	do i=1,iip1
722	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
723	plsnc(i,j,l)=pls(i,j,l)
724	enddo
725	enddo
726	enddo
727	do l=1,llm
728	do j=1,jjm
729	do i=1,iip1
730	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
731	enddo
732	enddo
733	enddo
734
735	! -----------------------------------------------------------------
736	! Interpolation champs guidage sur niveaux modele (+inversion N/S)
737	! Conversion en variables gcm (ucov, vcov...)
738	! -----------------------------------------------------------------
739	if (guide_P) then
740	do j=1,jjp1
741	do i=1,iim
742	ij=(j-1)*iip1+i
743	psgui1(ij)=psnat1(i,j)
744	psgui2(ij)=psnat2(i,j)
745	enddo
746	psgui1(iip1*j)=psnat1(1,j)
747	psgui2(iip1*j)=psnat2(1,j)
748	enddo
749	endif
750
751	IF (guide_u) THEN
752	CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
753	CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
754	do l=1,llm
755	do j=1,jjp1
756	do i=1,iim
757	ij=(j-1)*iip1+i
758	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
759	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
760	enddo
761	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
762	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
763	enddo
764	do i=1,iip1
765	ugui1(i,l)=0.
766	ugui1(ip1jm+i,l)=0.
767	ugui2(i,l)=0.
768	ugui2(ip1jm+i,l)=0.
769	enddo
770	enddo
771	ENDIF
772
773	IF (guide_T) THEN
774	CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
775	CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
776	do l=1,llm
777	do j=1,jjp1
778	IF (guide_teta) THEN
779	do i=1,iim
780	ij=(j-1)*iip1+i
781	tgui1(ij,l)=zu1(i,j,l)
782	tgui2(ij,l)=zu2(i,j,l)
783	enddo
784	ELSE
785	do i=1,iim
786	ij=(j-1)*iip1+i
787	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
788	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
789	enddo
790	ENDIF
791	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
792	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
793	enddo
794	do i=1,iip1
795	tgui1(i,l)=tgui1(1,l)
796	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
797	tgui2(i,l)=tgui2(1,l)
798	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
799	enddo
800	enddo
801	ENDIF
802
803	IF (guide_v) THEN
804
805	CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
806	CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
807
808	do l=1,llm
809	do j=1,jjm
810	do i=1,iim
811	ij=(j-1)*iip1+i
812	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
813	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
814	enddo
815	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
816	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
817	enddo
818	enddo
819	ENDIF
820
821	IF (guide_Q) THEN
822	! On suppose qu'on a la bonne variable dans le fichier de guidage:
823	! Hum.Rel si guide_hr, Hum.Spec. sinon.
824	CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
825	CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
826	do l=1,llm
827	do j=1,jjp1
828	do i=1,iim
829	ij=(j-1)*iip1+i
830	qgui1(ij,l)=zu1(i,j,l)
831	qgui2(ij,l)=zu2(i,j,l)
832	enddo
833	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
834	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
835	enddo
836	do i=1,iip1
837	qgui1(i,l)=qgui1(1,l)
838	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
839	qgui2(i,l)=qgui2(1,l)
840	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
841	enddo
842	enddo
843	IF (guide_hr) THEN
844	CALL q_sat(iip1jjp1llm,teta*pk/cpp,plsnc,qsat)
845	qgui1=qgui1qsat0.01 !hum. rel. en %
846	qgui2=qgui2qsat0.01
847	ENDIF
848	ENDIF
849
850	END SUBROUTINE guide_interp
851
852	!=======================================================================
853	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
854
855	! Calcul des constantes de rappel alpha (=1/tau)
856
857	USE comconst_mod, ONLY: pi
858	USE serre_mod, ONLY: clon,clat,grossismy
859
860	implicit none
861
862	include "dimensions.h"
863	include "paramet.h"
864	include "comgeom2.h"
865
866	! input arguments :
867	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
868	INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
869	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
870	REAL, INTENT(IN) :: taumin,taumax
871	! output arguments:
872	REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
873
874	! local variables:
875	LOGICAL, SAVE :: first=.TRUE.
876	REAL, SAVE :: gamma,dxdy_min,dxdy_max
877	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
878	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
879	REAL, DIMENSION (iip1,jjm) :: dxdyv
880	real dxdy_
881	real zlat,zlon
882	real alphamin,alphamax,xi
883	integer i,j,ilon,ilat
884
885
886	alphamin=factt/taumax
887	alphamax=factt/taumin
888	IF (guide_reg.OR.guide_add) THEN
889	alpha=alphamax
890	!-----------------------------------------------------------------------
891	! guide_reg: alpha=alpha_min dans region, 0. sinon.
892	!-----------------------------------------------------------------------
893	IF (guide_reg) THEN
894	do j=1,pjm
895	do i=1,pim
896	if (typ.eq.2) then
897	zlat=rlatu(j)*180./pi
898	zlon=rlonu(i)*180./pi
899	elseif (typ.eq.1) then
900	zlat=rlatu(j)*180./pi
901	zlon=rlonv(i)*180./pi
902	elseif (typ.eq.3) then
903	zlat=rlatv(j)*180./pi
904	zlon=rlonv(i)*180./pi
905	endif
906	alpha(i,j)=alphamax/16.* &
907	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
908	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
909	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
910	(1.+tanh((lon_max_g-zlon)/tau_lon))
911	enddo
912	enddo
913	ENDIF
914	ELSE
915	!-----------------------------------------------------------------------
916	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
917	!-----------------------------------------------------------------------
918	!Calcul de l'aire des mailles
919	do j=2,jjm
920	do i=2,iip1
921	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
922	enddo
923	zdx(1,j)=zdx(iip1,j)
924	enddo
925	do j=2,jjm
926	do i=1,iip1
927	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
928	enddo
929	enddo
930	do i=1,iip1
931	zdx(i,1)=zdx(i,2)
932	zdx(i,jjp1)=zdx(i,jjm)
933	zdy(i,1)=zdy(i,2)
934	zdy(i,jjp1)=zdy(i,jjm)
935	enddo
936	do j=1,jjp1
937	do i=1,iip1
938	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
939	enddo
940	enddo
941	IF (typ.EQ.2) THEN
942	do j=1,jjp1
943	do i=1,iim
944	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
945	enddo
946	dxdyu(iip1,j)=dxdyu(1,j)
947	enddo
948	ENDIF
949	IF (typ.EQ.3) THEN
950	do j=1,jjm
951	do i=1,iip1
952	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
953	enddo
954	enddo
955	ENDIF
956	! Premier appel: calcul des aires min et max et de gamma.
957	IF (first) THEN
958	first=.FALSE.
959	! coordonnees du centre du zoom
960	CALL coordij(clon,clat,ilon,ilat)
961	! aire de la maille au centre du zoom
962	dxdy_min=dxdys(ilon,ilat)
963	! dxdy maximale de la maille
964	dxdy_max=0.
965	do j=1,jjp1
966	do i=1,iip1
967	dxdy_max=max(dxdy_max,dxdys(i,j))
968	enddo
969	enddo
970	! Calcul de gamma
971	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
972	print*,'ATTENTION modele peu zoome'
973	print*,'ATTENTION on prend une constante de guidage cste'
974	gamma=0.
975	else
976	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
977	print*,'gamma=',gamma
978	if (gamma.lt.1.e-5) then
979	print*,'gamma =',gamma,'<1e-5'
980	stop
981	endif
982	gamma=log(0.5)/log(gamma)
983	if (gamma4) then
984	gamma=min(gamma,4.)
985	endif
986	print*,'gamma=',gamma
987	endif
988	ENDIF !first
989
990	do j=1,pjm
991	do i=1,pim
992	if (typ.eq.1) then
993	dxdy_=dxdys(i,j)
994	zlat=rlatu(j)*180./pi
995	elseif (typ.eq.2) then
996	dxdy_=dxdyu(i,j)
997	zlat=rlatu(j)*180./pi
998	elseif (typ.eq.3) then
999	dxdy_=dxdyv(i,j)
1000	zlat=rlatv(j)*180./pi
1001	endif
1002	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1003	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
1004	alpha(i,j)=alphamin
1005	else
1006	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
1007	xi=min(xi,1.)
1008	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
1009	alpha(i,j)=xialphamin+(1.-xi)alphamax
1010	else
1011	alpha(i,j)=0.
1012	endif
1013	endif
1014	enddo
1015	enddo
1016	ENDIF ! guide_reg
1017
1018	if (.not. guide_add) alpha = 1. - exp(- alpha)
1019
1020	END SUBROUTINE tau2alpha
1021
1022	!=======================================================================
1023	SUBROUTINE guide_read(timestep)
1024
1025	IMPLICIT NONE
1026
1027	#include "netcdf.inc"
1028	#include "dimensions.h"
1029	#include "paramet.h"
1030
1031	INTEGER, INTENT(IN) :: timestep
1032
1033	LOGICAL, SAVE :: first=.TRUE.
1034	! Identification fichiers et variables NetCDF:
1035	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1036	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1037	INTEGER :: ncidpl,varidpl,varidap,varidbp
1038	! Variables auxiliaires NetCDF:
1039	INTEGER, DIMENSION(4) :: start,count
1040	INTEGER :: status,rcode
1041
1042	CHARACTER (len = 80) :: abort_message
1043	CHARACTER (len = 20) :: modname = 'guide_read'
1044	! -----------------------------------------------------------------
1045	! Premier appel: initialisation de la lecture des fichiers
1046	! -----------------------------------------------------------------
1047	if (first) then
1048	ncidpl=-99
1049	print*,'Guide: ouverture des fichiers guidage '
1050	! Niveaux de pression si non constants
1051	if (guide_modele) then
1052	print *,'Lecture du guidage sur niveaux modele'
1053	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1054	IF (rcode.NE.NF_NOERR) THEN
1055	print *,'Guide: probleme -> pas de fichier apbp.nc'
1056	CALL abort_gcm(modname,abort_message,1)
1057	ENDIF
1058	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1059	IF (rcode.NE.NF_NOERR) THEN
1060	print *,'Guide: probleme -> pas de variable AP, fichier apbp.nc'
1061	CALL abort_gcm(modname,abort_message,1)
1062	ENDIF
1063	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1064	IF (rcode.NE.NF_NOERR) THEN
1065	print *,'Guide: probleme -> pas de variable BP, fichier apbp.nc'
1066	CALL abort_gcm(modname,abort_message,1)
1067	ENDIF
1068	print*,'ncidpl,varidap',ncidpl,varidap
1069	endif
1070	! Vent zonal
1071	if (guide_u) then
1072	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1073	IF (rcode.NE.NF_NOERR) THEN
1074	print *,'Guide: probleme -> pas de fichier u.nc'
1075	CALL abort_gcm(modname,abort_message,1)
1076	ENDIF
1077	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1078	IF (rcode.NE.NF_NOERR) THEN
1079	print *,'Guide: probleme -> pas de variable UWND, fichier u.nc'
1080	CALL abort_gcm(modname,abort_message,1)
1081	ENDIF
1082	print*,'ncidu,varidu',ncidu,varidu
1083	if (ncidpl.eq.-99) ncidpl=ncidu
1084	endif
1085	! Vent meridien
1086	if (guide_v) then
1087	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1088	IF (rcode.NE.NF_NOERR) THEN
1089	print *,'Guide: probleme -> pas de fichier v.nc'
1090	CALL abort_gcm(modname,abort_message,1)
1091	ENDIF
1092	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1093	IF (rcode.NE.NF_NOERR) THEN
1094	print *,'Guide: probleme -> pas de variable VWND, fichier v.nc'
1095	CALL abort_gcm(modname,abort_message,1)
1096	ENDIF
1097	print*,'ncidv,varidv',ncidv,varidv
1098	if (ncidpl.eq.-99) ncidpl=ncidv
1099	endif
1100	! Temperature
1101	if (guide_T) then
1102	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1103	IF (rcode.NE.NF_NOERR) THEN
1104	print *,'Guide: probleme -> pas de fichier T.nc'
1105	CALL abort_gcm(modname,abort_message,1)
1106	ENDIF
1107	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1108	IF (rcode.NE.NF_NOERR) THEN
1109	print *,'Guide: probleme -> pas de variable AIR, fichier T.nc'
1110	CALL abort_gcm(modname,abort_message,1)
1111	ENDIF
1112	print*,'ncidT,varidT',ncidt,varidt
1113	if (ncidpl.eq.-99) ncidpl=ncidt
1114	endif
1115	! Humidite
1116	if (guide_Q) then
1117	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1118	IF (rcode.NE.NF_NOERR) THEN
1119	print *,'Guide: probleme -> pas de fichier hur.nc'
1120	CALL abort_gcm(modname,abort_message,1)
1121	ENDIF
1122	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1123	IF (rcode.NE.NF_NOERR) THEN
1124	print *,'Guide: probleme -> pas de variable RH, fichier hur.nc'
1125	CALL abort_gcm(modname,abort_message,1)
1126	ENDIF
1127	print*,'ncidQ,varidQ',ncidQ,varidQ
1128	if (ncidpl.eq.-99) ncidpl=ncidQ
1129	endif
1130	! Pression de surface
1131	if ((guide_P).OR.(guide_modele)) then
1132	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1133	IF (rcode.NE.NF_NOERR) THEN
1134	print *,'Guide: probleme -> pas de fichier ps.nc'
1135	CALL abort_gcm(modname,abort_message,1)
1136	ENDIF
1137	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1138	IF (rcode.NE.NF_NOERR) THEN
1139	print *,'Guide: probleme -> pas de variable SP, fichier ps.nc'
1140	CALL abort_gcm(modname,abort_message,1)
1141	ENDIF
1142	print*,'ncidps,varidps',ncidps,varidps
1143	endif
1144	! Coordonnee verticale
1145	if (.not.guide_modele) then
1146	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1147	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1148	print*,'ncidpl,varidpl',ncidpl,varidpl
1149	endif
1150	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1151	if (guide_modele) then
1152	#ifdef NC_DOUBLE
1153	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1154	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1155	#else
1156	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1157	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1158	#endif
1159	else
1160	#ifdef NC_DOUBLE
1161	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1162	#else
1163	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1164	#endif
1165	apnc=apnc*100.! conversion en Pascals
1166	bpnc(:)=0.
1167	endif
1168	first=.FALSE.
1169	endif ! (first)
1170
1171	! -----------------------------------------------------------------
1172	! lecture des champs u, v, T, Q, ps
1173	! -----------------------------------------------------------------
1174
1175	! dimensions pour les champs scalaires et le vent zonal
1176	start(1)=1
1177	start(2)=1
1178	start(3)=1
1179	start(4)=timestep
1180
1181	count(1)=iip1
1182	count(2)=jjp1
1183	count(3)=nlevnc
1184	count(4)=1
1185
1186	! Vent zonal
1187	if (guide_u) then
1188	#ifdef NC_DOUBLE
1189	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
1190	#else
1191	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
1192	#endif
1193	IF (invert_y) THEN
1194	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1195	ENDIF
1196	endif
1197
1198	! Temperature
1199	if (guide_T) then
1200	#ifdef NC_DOUBLE
1201	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
1202	#else
1203	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
1204	#endif
1205	IF (invert_y) THEN
1206	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1207	ENDIF
1208	endif
1209
1210	! Humidite
1211	if (guide_Q) then
1212	#ifdef NC_DOUBLE
1213	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
1214	#else
1215	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
1216	#endif
1217	IF (invert_y) THEN
1218	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1219	ENDIF
1220
1221	endif
1222
1223	! Vent meridien
1224	if (guide_v) then
1225	count(2)=jjm
1226	#ifdef NC_DOUBLE
1227	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
1228	#else
1229	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
1230	#endif
1231	IF (invert_y) THEN
1232	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1233	ENDIF
1234	endif
1235
1236	! Pression de surface
1237	if ((guide_P).OR.(guide_modele)) then
1238	start(3)=timestep
1239	start(4)=0
1240	count(2)=jjp1
1241	count(3)=1
1242	count(4)=0
1243	#ifdef NC_DOUBLE
1244	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
1245	#else
1246	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
1247	#endif
1248	IF (invert_y) THEN
1249	CALL invert_lat(iip1,jjp1,1,psnat2)
1250	ENDIF
1251	endif
1252
1253	END SUBROUTINE guide_read
1254
1255	!=======================================================================
1256	SUBROUTINE guide_read2D(timestep)
1257
1258	IMPLICIT NONE
1259
1260	#include "netcdf.inc"
1261	#include "dimensions.h"
1262	#include "paramet.h"
1263
1264	INTEGER, INTENT(IN) :: timestep
1265
1266	LOGICAL, SAVE :: first=.TRUE.
1267	! Identification fichiers et variables NetCDF:
1268	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1269	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1270	INTEGER :: ncidpl,varidpl,varidap,varidbp
1271	! Variables auxiliaires NetCDF:
1272	INTEGER, DIMENSION(4) :: start,count
1273	INTEGER :: status,rcode
1274	! Variables for 3D extension:
1275	REAL, DIMENSION (jjp1,llm) :: zu
1276	REAL, DIMENSION (jjm,llm) :: zv
1277	INTEGER :: i
1278
1279	CHARACTER (len = 80) :: abort_message
1280	CHARACTER (len = 20) :: modname = 'guide_read2D'
1281	! -----------------------------------------------------------------
1282	! Premier appel: initialisation de la lecture des fichiers
1283	! -----------------------------------------------------------------
1284	if (first) then
1285	ncidpl=-99
1286	print*,'Guide: ouverture des fichiers guidage '
1287	! Niveaux de pression si non constants
1288	if (guide_modele) then
1289	print *,'Lecture du guidage sur niveaux modele'
1290	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1291	IF (rcode.NE.NF_NOERR) THEN
1292	print *,'Guide: probleme -> pas de fichier apbp.nc'
1293	CALL abort_gcm(modname,abort_message,1)
1294	ENDIF
1295	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1296	IF (rcode.NE.NF_NOERR) THEN
1297	print *,'Guide: probleme -> pas de variable AP, fichier apbp.nc'
1298	CALL abort_gcm(modname,abort_message,1)
1299	ENDIF
1300	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1301	IF (rcode.NE.NF_NOERR) THEN
1302	print *,'Guide: probleme -> pas de variable BP, fichier apbp.nc'
1303	CALL abort_gcm(modname,abort_message,1)
1304	ENDIF
1305	print*,'ncidpl,varidap',ncidpl,varidap
1306	endif
1307	! Vent zonal
1308	if (guide_u) then
1309	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1310	IF (rcode.NE.NF_NOERR) THEN
1311	print *,'Guide: probleme -> pas de fichier u.nc'
1312	CALL abort_gcm(modname,abort_message,1)
1313	ENDIF
1314	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1315	IF (rcode.NE.NF_NOERR) THEN
1316	print *,'Guide: probleme -> pas de variable UWND, fichier u.nc'
1317	CALL abort_gcm(modname,abort_message,1)
1318	ENDIF
1319	print*,'ncidu,varidu',ncidu,varidu
1320	if (ncidpl.eq.-99) ncidpl=ncidu
1321	endif
1322	! Vent meridien
1323	if (guide_v) then
1324	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1325	IF (rcode.NE.NF_NOERR) THEN
1326	print *,'Guide: probleme -> pas de fichier v.nc'
1327	CALL abort_gcm(modname,abort_message,1)
1328	ENDIF
1329	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1330	IF (rcode.NE.NF_NOERR) THEN
1331	print *,'Guide: probleme -> pas de variable VWND, fichier v.nc'
1332	CALL abort_gcm(modname,abort_message,1)
1333	ENDIF
1334	print*,'ncidv,varidv',ncidv,varidv
1335	if (ncidpl.eq.-99) ncidpl=ncidv
1336	endif
1337	! Temperature
1338	if (guide_T) then
1339	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1340	IF (rcode.NE.NF_NOERR) THEN
1341	print *,'Guide: probleme -> pas de fichier T.nc'
1342	CALL abort_gcm(modname,abort_message,1)
1343	ENDIF
1344	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1345	IF (rcode.NE.NF_NOERR) THEN
1346	print *,'Guide: probleme -> pas de variable AIR, fichier T.nc'
1347	CALL abort_gcm(modname,abort_message,1)
1348	ENDIF
1349	print*,'ncidT,varidT',ncidt,varidt
1350	if (ncidpl.eq.-99) ncidpl=ncidt
1351	endif
1352	! Humidite
1353	if (guide_Q) then
1354	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1355	IF (rcode.NE.NF_NOERR) THEN
1356	print *,'Guide: probleme -> pas de fichier hur.nc'
1357	CALL abort_gcm(modname,abort_message,1)
1358	ENDIF
1359	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1360	IF (rcode.NE.NF_NOERR) THEN
1361	print *,'Guide: probleme -> pas de variable RH, fichier hur.nc'
1362	CALL abort_gcm(modname,abort_message,1)
1363	ENDIF
1364	print*,'ncidQ,varidQ',ncidQ,varidQ
1365	if (ncidpl.eq.-99) ncidpl=ncidQ
1366	endif
1367	! Pression de surface
1368	if ((guide_P).OR.(guide_modele)) then
1369	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1370	IF (rcode.NE.NF_NOERR) THEN
1371	print *,'Guide: probleme -> pas de fichier ps.nc'
1372	CALL abort_gcm(modname,abort_message,1)
1373	ENDIF
1374	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1375	IF (rcode.NE.NF_NOERR) THEN
1376	print *,'Guide: probleme -> pas de variable SP, fichier ps.nc'
1377	CALL abort_gcm(modname,abort_message,1)
1378	ENDIF
1379	print*,'ncidps,varidps',ncidps,varidps
1380	endif
1381	! Coordonnee verticale
1382	if (.not.guide_modele) then
1383	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1384	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1385	print*,'ncidpl,varidpl',ncidpl,varidpl
1386	endif
1387	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1388	if (guide_modele) then
1389	#ifdef NC_DOUBLE
1390	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1391	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1392	#else
1393	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1394	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1395	#endif
1396	else
1397	#ifdef NC_DOUBLE
1398	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1399	#else
1400	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1401	#endif
1402	apnc=apnc*100.! conversion en Pascals
1403	bpnc(:)=0.
1404	endif
1405	first=.FALSE.
1406	endif ! (first)
1407
1408	! -----------------------------------------------------------------
1409	! lecture des champs u, v, T, Q, ps
1410	! -----------------------------------------------------------------
1411
1412	! dimensions pour les champs scalaires et le vent zonal
1413	start(1)=1
1414	start(2)=1
1415	start(3)=1
1416	start(4)=timestep
1417
1418	count(1)=1
1419	count(2)=jjp1
1420	count(3)=nlevnc
1421	count(4)=1
1422
1423	! Vent zonal
1424	if (guide_u) then
1425	#ifdef NC_DOUBLE
1426	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
1427	#else
1428	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
1429	#endif
1430	DO i=1,iip1
1431	unat2(i,:,:)=zu(:,:)
1432	ENDDO
1433
1434	IF (invert_y) THEN
1435	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1436	ENDIF
1437
1438	endif
1439
1440	! Temperature
1441	if (guide_T) then
1442	#ifdef NC_DOUBLE
1443	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
1444	#else
1445	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
1446	#endif
1447	DO i=1,iip1
1448	tnat2(i,:,:)=zu(:,:)
1449	ENDDO
1450
1451	IF (invert_y) THEN
1452	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1453	ENDIF
1454
1455	endif
1456
1457	! Humidite
1458	if (guide_Q) then
1459	#ifdef NC_DOUBLE
1460	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
1461	#else
1462	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
1463	#endif
1464	DO i=1,iip1
1465	qnat2(i,:,:)=zu(:,:)
1466	ENDDO
1467
1468	IF (invert_y) THEN
1469	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1470	ENDIF
1471
1472	endif
1473
1474	! Vent meridien
1475	if (guide_v) then
1476	count(2)=jjm
1477	#ifdef NC_DOUBLE
1478	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
1479	#else
1480	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
1481	#endif
1482	DO i=1,iip1
1483	vnat2(i,:,:)=zv(:,:)
1484	ENDDO
1485
1486	IF (invert_y) THEN
1487	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1488	ENDIF
1489
1490	endif
1491
1492	! Pression de surface
1493	if ((guide_P).OR.(guide_modele)) then
1494	start(3)=timestep
1495	start(4)=0
1496	count(2)=jjp1
1497	count(3)=1
1498	count(4)=0
1499	#ifdef NC_DOUBLE
1500	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
1501	#else
1502	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
1503	#endif
1504	DO i=1,iip1
1505	psnat2(i,:)=zu(:,1)
1506	ENDDO
1507
1508	IF (invert_y) THEN
1509	CALL invert_lat(iip1,jjp1,1,psnat2)
1510	ENDIF
1511
1512	endif
1513
1514	END SUBROUTINE guide_read2D
1515
1516	!=======================================================================
1517	SUBROUTINE guide_out(varname,hsize,vsize,field)
1518
1519	USE comvert_mod, ONLY: presnivs
1520	USE comconst_mod, ONLY: pi
1521
1522	IMPLICIT NONE
1523
1524	INCLUDE "dimensions.h"
1525	INCLUDE "paramet.h"
1526	INCLUDE "netcdf.inc"
1527	INCLUDE "comgeom2.h"
1528
1529	! Variables entree
1530	CHARACTER(), INTENT(IN) :: varname
1531	INTEGER, INTENT (IN) :: hsize,vsize
1532	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1533
1534	! Variables locales
1535	INTEGER, SAVE :: timestep=0
1536	! Identites fichier netcdf
1537	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1538	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1539	INTEGER :: vid_au,vid_av
1540	INTEGER, DIMENSION (3) :: dim3
1541	INTEGER, DIMENSION (4) :: dim4,count,start
1542	INTEGER :: ierr, varid,l
1543	REAL, DIMENSION (iip1,hsize,vsize) :: field2
1544
1545	print *,'Guide: output timestep',timestep,'var ',varname
1546	IF (timestep.EQ.0) THEN
1547	! ----------------------------------------------
1548	! initialisation fichier de sortie
1549	! ----------------------------------------------
1550	! Ouverture du fichier
1551	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
1552	! Definition des dimensions
1553	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1554	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1555	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1556	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1557	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1558	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1559
1560	! Creation des variables dimensions
1561	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1562	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1563	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1564	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1565	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1566	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1567	ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
1568	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1569	ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
1570
1571	ierr=NF_ENDDEF(nid)
1572
1573	! Enregistrement des variables dimensions
1574	#ifdef NC_DOUBLE
1575	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1576	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1577	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1578	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1579	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1580	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1581	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1582	ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
1583	ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
1584	#else
1585	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1586	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1587	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1588	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1589	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1590	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1591	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1592	ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
1593	ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
1594	#endif
1595	! --------------------------------------------------------------------
1596	! Cr�ation des variables sauvegard�es
1597	! --------------------------------------------------------------------
1598	ierr = NF_REDEF(nid)
1599	! Surface pressure (GCM)
1600	dim3=(/id_lonv,id_latu,id_tim/)
1601	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
1602	! Surface pressure (guidage)
1603	IF (guide_P) THEN
1604	dim3=(/id_lonv,id_latu,id_tim/)
1605	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1606	ENDIF
1607	! Zonal wind
1608	IF (guide_u) THEN
1609	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1610	ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
1611	ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
1612	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1613	ENDIF
1614	! Merid. wind
1615	IF (guide_v) THEN
1616	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1617	ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
1618	ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
1619	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1620	ENDIF
1621	! Pot. Temperature
1622	IF (guide_T) THEN
1623	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1624	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1625	ENDIF
1626	! Specific Humidity
1627	IF (guide_Q) THEN
1628	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1629	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1630	ENDIF
1631
1632	ierr = NF_ENDDEF(nid)
1633	ierr = NF_CLOSE(nid)
1634	ENDIF ! timestep=0
1635
1636	! --------------------------------------------------------------------
1637	! Enregistrement du champ
1638	! --------------------------------------------------------------------
1639	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1640
1641	IF (varname=="SP") timestep=timestep+1
1642
1643	ierr = NF_INQ_VARID(nid,varname,varid)
1644	SELECT CASE (varname)
1645	CASE ("SP","ps")
1646	start=(/1,1,timestep,0/)
1647	count=(/iip1,jjp1,1,0/)
1648	CASE ("v","va","vcov")
1649	start=(/1,1,1,timestep/)
1650	count=(/iip1,jjm,llm,1/)
1651	CASE DEFAULT
1652	start=(/1,1,1,timestep/)
1653	count=(/iip1,jjp1,llm,1/)
1654	END SELECT
1655
1656	SELECT CASE (varname)
1657	CASE("u","ua")
1658	DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
1659	field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
1660	CASE("v","va")
1661	DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
1662	CASE DEFAULT
1663	field2=field
1664	END SELECT
1665
1666
1667	#ifdef NC_DOUBLE
1668	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
1669	#else
1670	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
1671	#endif
1672
1673	ierr = NF_CLOSE(nid)
1674
1675	END SUBROUTINE guide_out
1676
1677
1678	!===========================================================================
1679	subroutine correctbid(iim,nl,x)
1680	integer iim,nl
1681	real x(iim+1,nl)
1682	integer i,l
1683	real zz
1684
1685	do l=1,nl
1686	do i=2,iim-1
1687	if(abs(x(i,l)).gt.1.e10) then
1688	zz=0.5*(x(i-1,l)+x(i+1,l))
1689	print*,'correction ',i,l,x(i,l),zz
1690	x(i,l)=zz
1691	endif
1692	enddo
1693	enddo
1694	return
1695	end subroutine correctbid
1696
1697	!===========================================================================
1698	END MODULE guide_mod

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

Download in other formats:

Original Format