Context Navigation

guide_loc_mod.F90 @ 2178

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

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

Download in other formats:

Original Format