Context Navigation

guide_p_mod.F90 @ 5017

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

In nudging procedures, replaced explicit Euler integration of nudged
fields by exact integration. This does not change anything if
guide_add is true, but it changes the value of alpha if guide_add is
false. We could have taken into account the variation of the nudging
field during a nudging time step. This would be a small correction. We
choose not to take it into account for the time being. Also, we add a
restriction on zonal nudging: we allow it only for a grid which is
regular in longitude. It does not seem to make sense otherwise and the
exact integration would take more programming for an irregular grid.

In the sequential version, copying the parallel versions, set
iguide_int to 1 when the input value is 0.

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.9 KB

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

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

Download in other formats:

Original Format