Context Navigation

guide_loc_mod.F90 @ 2597

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

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

Download in other formats:

Original Format