Context Navigation

guide_p_mod.F90 @ 2597

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

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

Download in other formats:

Original Format