Context Navigation

guide_loc_mod.F90 @ 1720

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

Fin du phasage de la dynamique parallele localisee (petite memoire) avec le tronc LMDZ5 r1671
Il reste quelques routines a verifier (en particulier ce qui touche a l'etude des cas academiques)
et la validation a effectuer

End of the phasing of the localised (low memory) parallel dynamics package with the
LMDZ5 trunk (r1671)
Some routines still need some checking (in particular the academic cases) and some
validation is still required

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 (pressure_exner) 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