Context Navigation

guide_p_mod.F90 @ 2100

Last change on this file since 2100 was 2100, checked in by lguez, 10 years ago

Removed "on rentre dans guide_main" from guide_main in dyn3dpar, was
already commented out in the dyn3dmem version.

Keeping length of lines under 80 characters in physiq (for
readability). Removed wrong comments "ajout des tendances de la
diffusion turbulente". Replaced "con" by "convection" as an argument
of add_phys_tend.

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory
Property svn:executable set to *

File size: 71.6 KB

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

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

Download in other formats:

Original Format