Context Navigation

guide_loc_mod.f90 @ 5509

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

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

Download in other formats:

Original Format