Context Navigation

guide_p_mod.F90 @ 1957

Last change on this file since 1957 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format