Context Navigation

guide_p_mod.F90 @ 3411

Last change on this file since 3411 was 2641, checked in by Laurent Fairhead, 8 years ago
Merged trunk changes r2593:2640 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``* Property svn:keywords set to `Id`
File size: 72.1 KB

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

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

Download in other formats:

Original Format