Context Navigation

guide_loc_mod.F90 @ 1630

Last change on this file since 1630 was 1630, checked in by Laurent Fairhead, 12 years ago

Importation initiale du rÃ©pertoire dyn3dmem

Initial import of dyn3dmem directory

File size: 68.0 KB

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

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

Download in other formats:

Original Format