Context Navigation

guide_loc_mod.F90 @ 5441

Last change on this file since 5441 was 2039, checked in by fhourdin, 11 years ago

Inclusion d'un appel supplementaire a geopot dans leapfrog dans

les versions dyn3d et dyn3dpar pour garantir la convergence
des 3 versions dyn3d/dyn3dpar/dyn3dmem.
La convergence fonctionne avec un compilateur gfortran 4.6.3 / openmpi
distribué sous ubuntu avec la nouvelle physique (NPv3.2) et le guidage.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Une seule exception : si on guide avec dyn3dpar et openMP.

Correction du guidage dans dyn3dmem.

un print supprime dans cv3_toutines.F90

Rm : On ne devrait perdre la convergence numérique avec les précédentes
versions que pour la nouvelle physique puisque le géopotentiel
n'est utilisé dans la physique que par les theermiques et sisvat.

=====================================================================

Added call to geopot in leapfrog in dyn3d and dyn3dpar in order to

insure numerical convergence with dyn3dmem.
The convergence dyn3d/dyn3dpar/dyn3dmem is satistied with
gfortran 4.6.3 / openmpi, new physics (NPv3.2) and nudging.
options : gfortran -fdefault-real-8 -DNC_DOUBLE -O3
Only one exception : with nudging & dynd3dpar & openMP.

Bug fixing for nudging in dyn3dmem.

print supressed in cv3_routines.F90

Rm : Numerical convergence with previous releases should be lost for new
physics only since the geopotential is used only by thermals and
sisvat in the physics.

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

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

Download in other formats:

Original Format