Context Navigation

guide_loc_mod.F90 @ 4067

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

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

Download in other formats:

Original Format