Context Navigation

guide_p_mod.F90 @ 1523

Last change on this file since 1523 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:executable set to *

File size: 64.1 KB

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

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

Download in other formats:

Original Format