Context Navigation

guide_loc_mod.F90 @ 3984

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

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

Download in other formats:

Original Format