Context Navigation

guide_p_mod.F90 @ 1443

Last change on this file since 1443 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 72.0 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	USE serre_mod, ONLY: grossismx
72
73	IMPLICIT NONE
74
75	INCLUDE "dimensions.h"
76	INCLUDE "paramet.h"
77	INCLUDE "netcdf.inc"
78
79	INTEGER :: error,ncidpl,rid,rcod
80	CHARACTER (len = 80) :: abort_message
81	CHARACTER (len = 20) :: modname = 'guide_init'
82
83	! ---------------------------------------------
84	! Lecture des parametres:
85	! ---------------------------------------------
86	! Variables guidees
87	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
88	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
89	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
90	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
91	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
92	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
93	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
94
95	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
96	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
97	if (guide_zon .and. abs(grossismx - 1.) > 0.01) &
98	call abort_gcm("guide_init", &
99	"zonal nudging requires grid regular in longitude", 1)
100
101	! Constantes de rappel. Unite : fraction de jour
102	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
103	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
104	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
105	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
106	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
107	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
108	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
109	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
110	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
111	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
112	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
113	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
114
115	! Sauvegarde du for�age
116	CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
117	CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
118	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
119	IF (iguide_sav.GT.0) THEN
120	iguide_sav=day_step/iguide_sav
121	ELSE if (iguide_sav == 0) then
122	iguide_sav = huge(0)
123	ELSE
124	iguide_sav=day_step*iguide_sav
125	ENDIF
126
127	! Guidage regional seulement (sinon constant ou suivant le zoom)
128	CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
129	CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
130	CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
131	CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
132	CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
133	CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
134	CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
135
136	! Parametres pour lecture des fichiers
137	CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
138	CALL getpar('iguide_int',4,iguide_int,'freq. interpolation vert')
139	IF (iguide_int.EQ.0) THEN
140	iguide_int=1
141	ELSEIF (iguide_int.GT.0) THEN
142	iguide_int=day_step/iguide_int
143	ELSE
144	iguide_int=day_step*iguide_int
145	ENDIF
146	CALL getpar('guide_plevs',0,guide_plevs,'niveaux pression fichiers guidage')
147	! Pour compatibilite avec ancienne version avec guide_modele
148	CALL getpar('guide_modele',.false.,guide_modele,'niveaux pression ap+bp*psol')
149	IF (guide_modele) THEN
150	guide_plevs=1
151	ENDIF
152	! Fin raccord
153	CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
154	CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
155	CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
156	CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
157
158	! ---------------------------------------------
159	! Determination du nombre de niveaux verticaux
160	! des fichiers guidage
161	! ---------------------------------------------
162	ncidpl=-99
163	if (guide_plevs.EQ.1) then
164	if (ncidpl.eq.-99) then
165	rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
166	if (rcod.NE.NF_NOERR) THEN
167	print *,'Guide: probleme -> pas de fichier apbp.nc'
168	CALL abort_gcm(modname,abort_message,1)
169	endif
170	endif
171	elseif (guide_plevs.EQ.2) then
172	if (ncidpl.EQ.-99) then
173	rcod=nf90_open('P.nc',Nf90_NOWRITe,ncidpl)
174	if (rcod.NE.NF_NOERR) THEN
175	print *,'Guide: probleme -> pas de fichier P.nc'
176	CALL abort_gcm(modname,abort_message,1)
177	endif
178	endif
179	elseif (guide_u) then
180	if (ncidpl.eq.-99) then
181	rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
182	if (rcod.NE.NF_NOERR) THEN
183	print *,'Guide: probleme -> pas de fichier u.nc'
184	CALL abort_gcm(modname,abort_message,1)
185	endif
186	endif
187	elseif (guide_v) then
188	if (ncidpl.eq.-99) then
189	rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
190	if (rcod.NE.NF_NOERR) THEN
191	print *,'Guide: probleme -> pas de fichier v.nc'
192	CALL abort_gcm(modname,abort_message,1)
193	endif
194	endif
195	elseif (guide_T) then
196	if (ncidpl.eq.-99) then
197	rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
198	if (rcod.NE.NF_NOERR) THEN
199	print *,'Guide: probleme -> pas de fichier T.nc'
200	CALL abort_gcm(modname,abort_message,1)
201	endif
202	endif
203	elseif (guide_Q) then
204	if (ncidpl.eq.-99) then
205	rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
206	if (rcod.NE.NF_NOERR) THEN
207	print *,'Guide: probleme -> pas de fichier hur.nc'
208	CALL abort_gcm(modname,abort_message,1)
209	endif
210	endif
211	endif
212	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
213	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
214	IF (error.NE.NF_NOERR) THEN
215	print *,'Guide: probleme lecture niveaux pression'
216	CALL abort_gcm(modname,abort_message,1)
217	ENDIF
218	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
219	print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
220	rcod = nf90_close(ncidpl)
221
222	! ---------------------------------------------
223	! Allocation des variables
224	! ---------------------------------------------
225	abort_message='pb in allocation guide'
226
227	ALLOCATE(apnc(nlevnc), stat = error)
228	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
229	ALLOCATE(bpnc(nlevnc), stat = error)
230	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
231	apnc=0.;bpnc=0.
232
233	ALLOCATE(alpha_pcor(llm), stat = error)
234	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
235	ALLOCATE(alpha_u(ip1jmp1), stat = error)
236	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237	ALLOCATE(alpha_v(ip1jm), stat = error)
238	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239	ALLOCATE(alpha_T(ip1jmp1), stat = error)
240	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
241	ALLOCATE(alpha_Q(ip1jmp1), stat = error)
242	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
243	ALLOCATE(alpha_P(ip1jmp1), stat = error)
244	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
245	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
246
247	IF (guide_u) THEN
248	ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
249	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
250	ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
251	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
252	ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
253	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
254	ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
255	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
256	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
257	ENDIF
258
259	IF (guide_T) THEN
260	ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
261	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
262	ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
263	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
264	ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
265	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
266	ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
267	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
268	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
269	ENDIF
270
271	IF (guide_Q) THEN
272	ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
273	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
274	ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
275	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
276	ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
277	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
278	ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
279	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
280	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
281	ENDIF
282
283	IF (guide_v) THEN
284	ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
285	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
286	ALLOCATE(vgui1(ip1jm,llm), stat = error)
287	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
288	ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
289	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
290	ALLOCATE(vgui2(ip1jm,llm), stat = error)
291	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
292	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
293	ENDIF
294
295	IF (guide_plevs.EQ.2) THEN
296	ALLOCATE(pnat1(iip1,jjp1,nlevnc), stat = error)
297	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
298	ALLOCATE(pnat2(iip1,jjp1,nlevnc), stat = error)
299	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
300	pnat1=0.;pnat2=0.;
301	ENDIF
302
303	IF (guide_P.OR.guide_plevs.EQ.1) THEN
304	ALLOCATE(psnat1(iip1,jjp1), stat = error)
305	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
306	ALLOCATE(psnat2(iip1,jjp1), stat = error)
307	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
308	psnat1=0.;psnat2=0.;
309	ENDIF
310	IF (guide_P) THEN
311	ALLOCATE(psgui2(ip1jmp1), stat = error)
312	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
313	ALLOCATE(psgui1(ip1jmp1), stat = error)
314	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
315	psgui1=0.;psgui2=0.
316	ENDIF
317
318	! ---------------------------------------------
319	! Lecture du premier etat de guidage.
320	! ---------------------------------------------
321	IF (guide_2D) THEN
322	CALL guide_read2D(1)
323	ELSE
324	CALL guide_read(1)
325	ENDIF
326	IF (guide_v) vnat1=vnat2
327	IF (guide_u) unat1=unat2
328	IF (guide_T) tnat1=tnat2
329	IF (guide_Q) qnat1=qnat2
330	IF (guide_plevs.EQ.2) pnat1=pnat2
331	IF (guide_P.OR.guide_plevs.EQ.1) psnat1=psnat2
332
333	END SUBROUTINE guide_init
334
335	!=======================================================================
336	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
337	use exner_hyb_p_m, only: exner_hyb_p
338	use exner_milieu_p_m, only: exner_milieu_p
339	USE parallel_lmdz
340	USE control_mod
341	USE comvert_mod, ONLY: ap,bp,preff,presnivs,pressure_exner
342	USE comconst_mod, ONLY: daysec,dtvr,kappa,cpp
343
344	IMPLICIT NONE
345
346	INCLUDE "dimensions.h"
347	INCLUDE "paramet.h"
348
349	! Variables entree
350	INTEGER, INTENT(IN) :: itau !pas de temps
351	REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
352	REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
353	REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
354
355	! Variables locales
356	LOGICAL, SAVE :: first=.TRUE.
357	LOGICAL :: f_out ! sortie guidage
358	REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
359	! Variables pour fonction Exner (P milieu couche)
360	REAL, DIMENSION (iip1,jjp1,llm) :: pk
361	REAL, DIMENSION (iip1,jjp1) :: pks
362	REAL :: unskap
363	REAL, DIMENSION (ip1jmp1,llmp1) :: p ! besoin si guide_P
364	! Compteurs temps:
365	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
366	REAL :: ditau, dday_step
367	REAL :: tau,reste ! position entre 2 etats de guidage
368	REAL, SAVE :: factt ! pas de temps en fraction de jour
369
370	INTEGER :: i,j,l
371
372	ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1
373	jjb_u=jj_begin ; jje_u=jj_end ; jjn_u=jje_u-jjb_u+1
374	ijb_v=ij_begin ; ije_v=ij_end ; ijn_v=ije_v-ijb_v+1
375	jjb_v=jj_begin ; jje_v=jj_end ; jjn_v=jje_v-jjb_v+1
376	IF (pole_sud) THEN
377	ije_v=ij_end-iip1
378	jje_v=jj_end-1
379	ijn_v=ije_v-ijb_v+1
380	jjn_v=jje_v-jjb_v+1
381	ENDIF
382
383	! CALL AllGather_Field(ucov,ip1jmp1,llm)
384	! CALL AllGather_Field(vcov,ip1jm,llm)
385	! CALL AllGather_Field(teta,ip1jmp1,llm)
386	! CALL AllGather_Field(ps,ip1jmp1,1)
387	! CALL AllGather_Field(q,ip1jmp1,llm)
388
389	!-----------------------------------------------------------------------
390	! Initialisations au premier passage
391	!-----------------------------------------------------------------------
392
393	IF (first) THEN
394	first=.FALSE.
395	CALL guide_init
396	itau_test=1001
397	step_rea=1
398	count_no_rea=0
399	! Calcul des constantes de rappel
400	factt=dtvr*iperiod/daysec
401	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
402	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
403	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
404	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
405	call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
406	! correction de rappel dans couche limite
407	if (guide_BL) then
408	alpha_pcor(:)=1.
409	else
410	do l=1,llm
411	alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
412	enddo
413	endif
414	! ini_anal: etat initial egal au guidage
415	IF (ini_anal) THEN
416	CALL guide_interp(ps,teta)
417	IF (guide_u) ucov(ijb_u:ije_u,:)=ugui2(ijb_u:ije_u,:)
418	IF (guide_v) vcov(ijb_v:ije_v,:)=ugui2(ijb_v:ije_v,:)
419	IF (guide_T) teta(ijb_u:ije_u,:)=tgui2(ijb_u:ije_u,:)
420	IF (guide_Q) q(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)
421	IF (guide_P) THEN
422	ps(ijb_u:ije_u)=psgui2(ijb_u:ije_u)
423	CALL pression_p(ip1jmp1,ap,bp,ps,p)
424	CALL massdair_p(p,masse)
425	ENDIF
426	RETURN
427	ENDIF
428	! Verification structure guidage
429	IF (guide_u) THEN
430	CALL writefield_p('unat',unat1)
431	CALL writefield_p('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
432	ENDIF
433	IF (guide_T) THEN
434	CALL writefield_p('tnat',tnat1)
435	CALL writefield_p('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
436	ENDIF
437
438	ENDIF !first
439
440	!-----------------------------------------------------------------------
441	! Lecture des fichiers de guidage ?
442	!-----------------------------------------------------------------------
443	IF (iguide_read.NE.0) THEN
444	ditau=real(itau)
445	dday_step=real(day_step)
446	IF (iguide_read.LT.0) THEN
447	tau=ditau/dday_step/REAL(iguide_read)
448	ELSE
449	tau=REAL(iguide_read)*ditau/dday_step
450	ENDIF
451	reste=tau-AINT(tau)
452	IF (reste.EQ.0.) THEN
453	IF (itau_test.EQ.itau) THEN
454	write(,)'deuxieme passage de advreel a itau=',itau
455	stop
456	ELSE
457	IF (guide_v) vnat1(:,jjb_v:jje_v,:)=vnat2(:,jjb_v:jje_v,:)
458	IF (guide_u) unat1(:,jjb_u:jje_u,:)=unat2(:,jjb_u:jje_u,:)
459	IF (guide_T) tnat1(:,jjb_u:jje_u,:)=tnat2(:,jjb_u:jje_u,:)
460	IF (guide_Q) qnat1(:,jjb_u:jje_u,:)=qnat2(:,jjb_u:jje_u,:)
461	IF (guide_plevs.EQ.2) pnat1(:,jjb_u:jje_u,:)=pnat2(:,jjb_u:jje_u,:)
462	IF (guide_P.OR.guide_plevs.EQ.1) psnat1(:,jjb_u:jje_u)=psnat2(:,jjb_u:jje_u)
463	step_rea=step_rea+1
464	itau_test=itau
465	print*,'Lecture fichiers guidage, pas ',step_rea, &
466	'apres ',count_no_rea,' non lectures'
467	IF (guide_2D) THEN
468	CALL guide_read2D(step_rea)
469	ELSE
470	CALL guide_read(step_rea)
471	ENDIF
472	count_no_rea=0
473	ENDIF
474	ELSE
475	count_no_rea=count_no_rea+1
476
477	ENDIF
478	ENDIF !iguide_read=0
479
480	!-----------------------------------------------------------------------
481	! Interpolation et conversion des champs de guidage
482	!-----------------------------------------------------------------------
483	IF (MOD(itau,iguide_int).EQ.0) THEN
484	CALL guide_interp(ps,teta)
485	ENDIF
486	! Repartition entre 2 etats de guidage
487	IF (iguide_read.NE.0) THEN
488	tau=reste
489	ELSE
490	tau=1.
491	ENDIF
492
493	!-----------------------------------------------------------------------
494	! Ajout des champs de guidage
495	!-----------------------------------------------------------------------
496	! Sauvegarde du guidage?
497	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
498	IF (f_out) THEN
499	! Calcul niveaux pression milieu de couches
500	CALL pression_p( ip1jmp1, ap, bp, ps, p )
501	if (pressure_exner) then
502	CALL exner_hyb_p(ip1jmp1,ps,p,pks,pk)
503	else
504	CALL exner_milieu_p(ip1jmp1,ps,p,pks,pk)
505	endif
506	unskap=1./kappa
507	DO l = 1, llm
508	DO j=jjb_u,jje_u
509	DO i =1, iip1
510	p(i+(j-1)iip1,l) = preff ( pk(i,j,l)/cpp) ** unskap
511	ENDDO
512	ENDDO
513	ENDDO
514	CALL guide_out("SP",jjp1,llm,p(1:ip1jmp1,1:llm),1.)
515	ENDIF
516
517	if (guide_u) then
518	if (guide_add) then
519	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)
520	else
521	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)-ucov(ijb_u:ije_u,:)
522	endif
523
524	if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
525	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
526	IF (f_out) CALL guide_out("u",jjp1,llm,ucov,factt)
527	IF (f_out) CALL guide_out("ua",jjp1,llm,(1.-tau)ugui1(:,:)+tauugui2(:,:),factt)
528	IF (f_out) CALL guide_out("ucov",jjp1,llm,f_add(:,:)/factt,factt)
529	ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
530	endif
531
532	if (guide_T) then
533	if (guide_add) then
534	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)
535	else
536	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)-teta(ijb_u:ije_u,:)
537	endif
538	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
539	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
540	IF (f_out) CALL guide_out("teta",jjp1,llm,f_add(:,:)/factt,factt)
541	teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
542	endif
543
544	if (guide_P) then
545	if (guide_add) then
546	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)
547	else
548	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)-ps(ijb_u:ije_u)
549	endif
550	if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
551	CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
552	IF (f_out) CALL guide_out("ps",jjp1,1,f_add(1:ip1jmp1,1)/factt,factt)
553	ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
554	CALL pression_p(ip1jmp1,ap,bp,ps,p)
555	CALL massdair_p(p,masse)
556	endif
557
558	if (guide_Q) then
559	if (guide_add) then
560	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)
561	else
562	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)-q(ijb_u:ije_u,:)
563	endif
564	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
565	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
566	IF (f_out) CALL guide_out("q",jjp1,llm,f_add(:,:)/factt,factt)
567	q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
568	endif
569
570	if (guide_v) then
571	if (guide_add) then
572	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)
573	else
574	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)-vcov(ijb_v:ije_v,:)
575	endif
576
577	if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
578	CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
579	IF (f_out) CALL guide_out("v",jjm,llm,vcov(1:ip1jm,:),factt)
580	IF (f_out) CALL guide_out("va",jjm,llm,(1.-tau)vgui1(:,:)+tauvgui2(:,:),factt)
581	IF (f_out) CALL guide_out("vcov",jjm,llm,f_add(1:ip1jm,:)/factt,factt)
582	vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
583	endif
584
585	END SUBROUTINE guide_main
586
587	!=======================================================================
588	SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
589	! field1=afield1+alphafield2
590
591	IMPLICIT NONE
592	INCLUDE "dimensions.h"
593	INCLUDE "paramet.h"
594
595	! input variables
596	INTEGER, INTENT(IN) :: hsize
597	INTEGER, INTENT(IN) :: vsize
598	REAL, DIMENSION(hsize), INTENT(IN) :: alpha
599	REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
600
601	! Local variables
602	INTEGER :: l
603
604	IF (hsize==ip1jm) THEN
605	do l=1,vsize
606	field(ijb_v:ije_v,l)=alpha(ijb_v:ije_v)field(ijb_v:ije_v,l)alpha_pcor(l)
607	enddo
608	ELSE
609	do l=1,vsize
610	field(ijb_u:ije_u,l)=alpha(ijb_u:ije_u)field(ijb_u:ije_u,l)alpha_pcor(l)
611	enddo
612	ENDIF
613
614	END SUBROUTINE guide_addfield
615
616	!=======================================================================
617	SUBROUTINE guide_zonave(typ,hsize,vsize,field)
618
619	USE comconst_mod, ONLY: pi
620
621	IMPLICIT NONE
622
623	INCLUDE "dimensions.h"
624	INCLUDE "paramet.h"
625	INCLUDE "comgeom.h"
626
627	! input/output variables
628	INTEGER, INTENT(IN) :: typ
629	INTEGER, INTENT(IN) :: vsize
630	INTEGER, INTENT(IN) :: hsize
631	REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
632
633	! Local variables
634	LOGICAL, SAVE :: first=.TRUE.
635	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
636	INTEGER :: i,j,l,ij
637	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
638	REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
639
640	IF (first) THEN
641	first=.FALSE.
642	!Compute domain for averaging
643	lond=rlonu*180./pi
644	imin(1)=1;imax(1)=iip1;
645	imin(2)=1;imax(2)=iip1;
646	IF (guide_reg) THEN
647	DO i=1,iim
648	IF (lond(i).LT.lon_min_g) imin(1)=i
649	IF (lond(i).LE.lon_max_g) imax(1)=i
650	ENDDO
651	lond=rlonv*180./pi
652	DO i=1,iim
653	IF (lond(i).LT.lon_min_g) imin(2)=i
654	IF (lond(i).LE.lon_max_g) imax(2)=i
655	ENDDO
656	ENDIF
657	ENDIF
658
659	fieldm=0.
660
661	IF (hsize==jjm) THEN
662	DO l=1,vsize
663	! Compute zonal average
664	DO j=jjb_v,jje_v
665	DO i=imin(typ),imax(typ)
666	ij=(j-1)*iip1+i
667	fieldm(j,l)=fieldm(j,l)+field(ij,l)
668	ENDDO
669	ENDDO
670	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
671	! Compute forcing
672	DO j=jjb_v,jje_v
673	DO i=1,iip1
674	ij=(j-1)*iip1+i
675	field(ij,l)=fieldm(j,l)
676	ENDDO
677	ENDDO
678	ENDDO
679	ELSE
680	DO l=1,vsize
681	! Compute zonal average
682	DO j=jjb_v,jje_v
683	DO i=imin(typ),imax(typ)
684	ij=(j-1)*iip1+i
685	fieldm(j,l)=fieldm(j,l)+field(ij,l)
686	ENDDO
687	ENDDO
688	fieldm(:,l)=fieldm(:,l)/REAL(imax(typ)-imin(typ)+1)
689	! Compute forcing
690	DO j=jjb_u,jje_u
691	DO i=1,iip1
692	ij=(j-1)*iip1+i
693	field(ij,l)=fieldm(j,l)
694	ENDDO
695	ENDDO
696	ENDDO
697	ENDIF
698
699	END SUBROUTINE guide_zonave
700
701	!=======================================================================
702	SUBROUTINE guide_interp(psi,teta)
703	use exner_hyb_p_m, only: exner_hyb_p
704	use exner_milieu_p_m, only: exner_milieu_p
705	USE parallel_lmdz
706	USE mod_hallo
707	USE Bands
708	USE comvert_mod, ONLY: ap,bp,preff,pressure_exner
709	USE comconst_mod, ONLY: kappa,cpp
710	IMPLICIT NONE
711
712	include "dimensions.h"
713	include "paramet.h"
714	include "comgeom2.h"
715
716	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
717	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
718
719	LOGICAL, SAVE :: first=.TRUE.
720	! Variables pour niveaux pression:
721	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
722	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
723	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
724	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
725	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
726	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
727	REAL, DIMENSION (iip1,jjm,llm) :: pbary
728	! Variables pour fonction Exner (P milieu couche)
729	REAL, DIMENSION (iip1,jjp1,llm) :: pk
730	REAL, DIMENSION (iip1,jjp1) :: pks
731	REAL :: unskap
732	! Pression de vapeur saturante
733	REAL, DIMENSION (ip1jmp1,llm) :: qsat
734	!Variables intermediaires interpolation
735	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
736	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
737
738	INTEGER :: i,j,l,ij
739	TYPE(Request) :: Req
740
741	print *,'Guide: conversion variables guidage'
742	! -----------------------------------------------------------------
743	! Calcul des niveaux de pression champs guidage (pour T et Q)
744	! -----------------------------------------------------------------
745	IF (guide_plevs.EQ.0) THEN
746	DO l=1,nlevnc
747	DO j=jjb_u,jje_u
748	DO i=1,iip1
749	plnc2(i,j,l)=apnc(l)
750	plnc1(i,j,l)=apnc(l)
751	ENDDO
752	ENDDO
753	ENDDO
754	ENDIF
755
756	if (first) then
757	first=.FALSE.
758	print*,'Guide: verification ordre niveaux verticaux'
759	print*,'LMDZ :'
760	do l=1,llm
761	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
762	+psi(1,jje_u)*(bp(l)+bp(l+1))/2.
763	enddo
764	print*,'Fichiers guidage'
765	SELECT CASE (guide_plevs)
766	CASE (0)
767	do l=1,nlevnc
768	print*,'PL(',l,')=',plnc2(1,jjb_u,l)
769	enddo
770	CASE (1)
771	DO l=1,nlevnc
772	print,'PL(',l,')=',apnc(l)+bpnc(l)psnat2(i,jjb_u)
773	ENDDO
774	CASE (2)
775	do l=1,nlevnc
776	print*,'PL(',l,')=',pnat2(1,jjb_u,l)
777	enddo
778	END SELECT
779	print *,'inversion de l''ordre: invert_p=',invert_p
780	if (guide_u) then
781	do l=1,nlevnc
782	print*,'U(',l,')=',unat2(1,jjb_u,l)
783	enddo
784	endif
785	if (guide_T) then
786	do l=1,nlevnc
787	print*,'T(',l,')=',tnat2(1,jjb_u,l)
788	enddo
789	endif
790	endif
791
792	! -----------------------------------------------------------------
793	! Calcul niveaux pression modele
794	! -----------------------------------------------------------------
795
796	! .... Calcul de pls , pression au milieu des couches ,en Pascals
797	IF (guide_plevs.EQ.1) THEN
798	DO l=1,llm
799	DO j=jjb_u,jje_u
800	DO i =1, iip1
801	pls(i,j,l)=(ap(l)+ap(l+1))/2.+psi(i,j)*(bp(l)+bp(l+1))/2.
802	ENDDO
803	ENDDO
804	ENDDO
805	ELSE
806	CALL pression_p( ip1jmp1, ap, bp, psi, p )
807	if (pressure_exner) then
808	CALL exner_hyb_p(ip1jmp1,psi,p,pks,pk)
809	else
810	CALL exner_milieu_p(ip1jmp1,psi,p,pks,pk)
811	endif
812	unskap=1./kappa
813	DO l = 1, llm
814	DO j=jjb_u,jje_u
815	DO i =1, iip1
816	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
817	ENDDO
818	ENDDO
819	ENDDO
820	ENDIF
821
822	! calcul des pressions pour les grilles u et v
823	do l=1,llm
824	do j=jjb_u,jje_u
825	do i=1,iip1
826	pext(i,j,l)=pls(i,j,l)*aire(i,j)
827	enddo
828	enddo
829	enddo
830
831	CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
832	CALL SendRequest(Req)
833	CALL WaitRequest(Req)
834
835	call massbar_p(pext, pbarx, pbary )
836	do l=1,llm
837	do j=jjb_u,jje_u
838	do i=1,iip1
839	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
840	plsnc(i,j,l)=pls(i,j,l)
841	enddo
842	enddo
843	enddo
844	do l=1,llm
845	do j=jjb_v,jje_v
846	do i=1,iip1
847	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
848	enddo
849	enddo
850	enddo
851
852	! -----------------------------------------------------------------
853	! Interpolation verticale champs guidage sur niveaux modele
854	! Conversion en variables gcm (ucov, vcov...)
855	! -----------------------------------------------------------------
856	if (guide_P) then
857	do j=jjb_u,jje_u
858	do i=1,iim
859	ij=(j-1)*iip1+i
860	psgui1(ij)=psnat1(i,j)
861	psgui2(ij)=psnat2(i,j)
862	enddo
863	psgui1(iip1*j)=psnat1(1,j)
864	psgui2(iip1*j)=psnat2(1,j)
865	enddo
866	endif
867
868	IF (guide_T) THEN
869	! Calcul des nouvelles valeurs des niveaux de pression du guidage
870	IF (guide_plevs.EQ.1) THEN
871	DO l=1,nlevnc
872	DO j=jjb_u,jje_u
873	DO i=1,iip1
874	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
875	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
876	ENDDO
877	ENDDO
878	ENDDO
879	ELSE IF (guide_plevs.EQ.2) THEN
880	DO l=1,nlevnc
881	DO j=jjb_u,jje_u
882	DO i=1,iip1
883	plnc2(i,j,l)=pnat2(i,j,l)
884	plnc1(i,j,l)=pnat1(i,j,l)
885	ENDDO
886	ENDDO
887	ENDDO
888	ENDIF
889
890	! Interpolation verticale
891	CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
892	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
893	CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
894	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
895
896	! Conversion en variables GCM
897	do l=1,llm
898	do j=jjb_u,jje_u
899	IF (guide_teta) THEN
900	do i=1,iim
901	ij=(j-1)*iip1+i
902	tgui1(ij,l)=zu1(i,j,l)
903	tgui2(ij,l)=zu2(i,j,l)
904	enddo
905	ELSE
906	do i=1,iim
907	ij=(j-1)*iip1+i
908	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
909	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
910	enddo
911	ENDIF
912	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
913	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
914	enddo
915	do i=1,iip1
916	tgui1(i,l)=tgui1(1,l)
917	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
918	tgui2(i,l)=tgui2(1,l)
919	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
920	enddo
921	enddo
922	ENDIF
923
924	IF (guide_Q) THEN
925	! Calcul des nouvelles valeurs des niveaux de pression du guidage
926	IF (guide_plevs.EQ.1) THEN
927	DO l=1,nlevnc
928	DO j=jjb_u,jje_u
929	DO i=1,iip1
930	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
931	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
932	ENDDO
933	ENDDO
934	ENDDO
935	ELSE IF (guide_plevs.EQ.2) THEN
936	DO l=1,nlevnc
937	DO j=jjb_u,jje_u
938	DO i=1,iip1
939	plnc2(i,j,l)=pnat2(i,j,l)
940	plnc1(i,j,l)=pnat1(i,j,l)
941	ENDDO
942	ENDDO
943	ENDDO
944	ENDIF
945
946	! Interpolation verticale
947	CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
948	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
949	CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
950	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
951
952	! Conversion en variables GCM
953	! On suppose qu'on a la bonne variable dans le fichier de guidage:
954	! Hum.Rel si guide_hr, Hum.Spec. sinon.
955	do l=1,llm
956	do j=jjb_u,jje_u
957	do i=1,iim
958	ij=(j-1)*iip1+i
959	qgui1(ij,l)=zu1(i,j,l)
960	qgui2(ij,l)=zu2(i,j,l)
961	enddo
962	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
963	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
964	enddo
965	do i=1,iip1
966	qgui1(i,l)=qgui1(1,l)
967	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
968	qgui2(i,l)=qgui2(1,l)
969	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
970	enddo
971	enddo
972	IF (guide_hr) THEN
973	CALL q_sat(iip1jjn_ullm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp, &
974	plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
975	qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01 !hum. rel. en %
976	qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01
977	ENDIF
978	ENDIF
979
980	IF (guide_u) THEN
981	! Calcul des nouvelles valeurs des niveaux de pression du guidage
982	IF (guide_plevs.EQ.1) THEN
983	DO l=1,nlevnc
984	DO j=jjb_u,jje_u
985	DO i=1,iim
986	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha1p2(i,j) &
987	& +psnat2(i+1,j)aire(i+1,j)alpha3p4(i+1,j))/aireu(i,j)
988	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha1p2(i,j) &
989	& +psnat1(i+1,j)aire(i+1,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	ELSE IF (guide_plevs.EQ.2) THEN
996	DO l=1,nlevnc
997	DO j=jjb_u,jje_u
998	DO i=1,iim
999	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha1p2(i,j) &
1000	& +pnat2(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
1001	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha1p2(i,j) &
1002	& +pnat1(i+1,j,l)aire(i,j)alpha3p4(i+1,j))/aireu(i,j)
1003	ENDDO
1004	plnc2(iip1,j,l)=plnc2(1,j,l)
1005	plnc1(iip1,j,l)=plnc1(1,j,l)
1006	ENDDO
1007	ENDDO
1008	ENDIF
1009
1010	! Interpolation verticale
1011	CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
1012	plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
1013	CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
1014	plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
1015
1016	! Conversion en variables GCM
1017	do l=1,llm
1018	do j=jjb_u,jje_u
1019	do i=1,iim
1020	ij=(j-1)*iip1+i
1021	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
1022	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
1023	enddo
1024	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
1025	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
1026	enddo
1027	do i=1,iip1
1028	ugui1(i,l)=0.
1029	ugui1(ip1jm+i,l)=0.
1030	ugui2(i,l)=0.
1031	ugui2(ip1jm+i,l)=0.
1032	enddo
1033	enddo
1034	ENDIF
1035
1036	IF (guide_v) THEN
1037	! Calcul des nouvelles valeurs des niveaux de pression du guidage
1038	IF (guide_plevs.EQ.1) THEN
1039	CALL Register_SwapFieldHallo(psnat1,psnat1,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
1040	CALL SendRequest(Req)
1041	CALL WaitRequest(Req)
1042	CALL Register_SwapFieldHallo(psnat2,psnat2,ip1jmp1,1,jj_Nb_caldyn,1,2,Req)
1043	CALL SendRequest(Req)
1044	CALL WaitRequest(Req)
1045	DO l=1,nlevnc
1046	DO j=jjb_v,jje_v
1047	DO i=1,iip1
1048	plnc2(i,j,l)=apnc(l)+bpnc(l)(psnat2(i,j)aire(i,j)*alpha2p3(i,j) &
1049	& +psnat2(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
1050	plnc1(i,j,l)=apnc(l)+bpnc(l)(psnat1(i,j)aire(i,j)*alpha2p3(i,j) &
1051	& +psnat1(i,j+1)aire(i,j+1)alpha1p4(i,j+1))/airev(i,j)
1052	ENDDO
1053	ENDDO
1054	ENDDO
1055	ELSE IF (guide_plevs.EQ.2) THEN
1056	CALL Register_SwapFieldHallo(pnat1,pnat1,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
1057	CALL SendRequest(Req)
1058	CALL WaitRequest(Req)
1059	CALL Register_SwapFieldHallo(pnat2,pnat2,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
1060	CALL SendRequest(Req)
1061	CALL WaitRequest(Req)
1062	DO l=1,nlevnc
1063	DO j=jjb_v,jje_v
1064	DO i=1,iip1
1065	plnc2(i,j,l)=(pnat2(i,j,l)aire(i,j)alpha2p3(i,j) &
1066	& +pnat2(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
1067	plnc1(i,j,l)=(pnat1(i,j,l)aire(i,j)alpha2p3(i,j) &
1068	& +pnat1(i,j+1,l)aire(i,j)alpha1p4(i,j+1))/airev(i,j)
1069	ENDDO
1070	ENDDO
1071	ENDDO
1072	ENDIF
1073	! Interpolation verticale
1074	CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm, &
1075	plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
1076	CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm, &
1077	plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
1078	! Conversion en variables GCM
1079	do l=1,llm
1080	do j=jjb_v,jje_v
1081	do i=1,iim
1082	ij=(j-1)*iip1+i
1083	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
1084	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
1085	enddo
1086	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
1087	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
1088	enddo
1089	enddo
1090	ENDIF
1091
1092
1093	END SUBROUTINE guide_interp
1094
1095	!=======================================================================
1096	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
1097
1098	! Calcul des constantes de rappel alpha (=1/tau)
1099
1100	USE comconst_mod, ONLY: pi
1101	USE serre_mod, ONLY: clon,clat,grossismx,grossismy
1102
1103	implicit none
1104
1105	include "dimensions.h"
1106	include "paramet.h"
1107	include "comgeom2.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	USE comvert_mod, ONLY: presnivs
1816	USE comconst_mod, ONLY: pi
1817	IMPLICIT NONE
1818
1819	INCLUDE "dimensions.h"
1820	INCLUDE "paramet.h"
1821	INCLUDE "netcdf.inc"
1822	INCLUDE "comgeom2.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