Context Navigation

guide_loc_mod.F90 @ 2461

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

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

Download in other formats:

Original Format