Context Navigation

guide_loc_mod.F90 @ 2034

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

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

Download in other formats:

Original Format