Context Navigation

guide_loc_mod.F90 @ 5444

Last change on this file since 5444 was 2263, checked in by lguez, 10 years ago

In the sequential version, ini_getparam and fin_getparam were not
called so guide_init created a file "fort.99". Added call to
ini_getparam and fin_getparam in guide_init. The created file is now
"nudging_parameters_out.txt".

In the parallel version, ini_getparam was called from gcm, with the
argument "out.def". So "out.def" was created even without nudging, it
remained empty. Moved the call to ini_getparam from gcm to
guide_init. Also changed the name of the created file to
"nudging_parameters_out.txt", as in the sequential version. "out.def"
was too vague and confusing a name.

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.5 KB

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

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

Download in other formats:

Original Format