Context Navigation

guide_loc_mod.F90 @ 4246

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

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

Download in other formats:

Original Format