Context Navigation

guide_loc_mod.F90 @ 4400

Last change on this file since 4400 was 1848, checked in by yann meurdesoif, 11 years ago

Solve performance problem comming from declarations of a derived type.

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

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

Download in other formats:

Original Format