Context Navigation

guide_p_mod.F90 @ 1200

Last change on this file since 1200 was 1188, checked in by lguez, 15 years ago
Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6 Fortran 90 interface. Corrected a few comments in "output.def".
Property svn:executable set to ``*
File size: 56.0 KB

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

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

Download in other formats:

Original Format