Context Navigation

guide_loc_mod.F90 @ 5420

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

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

Download in other formats:

Original Format