Context Navigation

guide_p_mod.F90 @ 2024

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

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

Download in other formats:

Original Format