Context Navigation

guide_loc_mod.F90 @ 5069

Last change on this file since 5069 was 5069, checked in by abarral, 7 weeks ago
Reduce use of #ifdef NC_DOUBLE to single instance in lmdz_netcdf.F90 Add nf_get_vara_rd in lmdz_netcdf.F90 Remove #ifdef NC_DOUBLE in dynredem_mod.F90 & guide_loc_mod.F90 (minor) fix some casting in ncdf calls in guide_loc_mod.F90 (minor) replace netcdf call & reduncate implicit none in dynredem_mod.F90
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 84.4 KB

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

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

Download in other formats:

Original Format