Context Navigation

guide_p_mod.F90 @ 2569

Last change on this file since 2569 was 2160, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes -r2070:2158 into testing branch. Compilation problems introduced by revision r2155 have been corrected by hand
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 71.9 KB

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

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

Download in other formats:

Original Format