Context Navigation

guide_p_mod.F90 @ 5360

Last change on this file since 5360 was 2025, checked in by fhourdin, 11 years ago
Enrichissement des sorties pour le guidage Enriched outputs for nudging
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 ``*
File size: 71.6 KB

Line
1	!
2	! $Id$
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