Context Navigation

guide_p_mod.F90 @ 2056

Last change on this file since 2056 was 2056, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1997:2055 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 71.6 KB

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

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

Download in other formats:

Original Format