Context Navigation

guide_loc_mod.F90 @ 5136

Last change on this file since 5136 was 5084, checked in by Laurent Fairhead, 4 months ago
Reverting to r4065. Updating fortran standard broke too much stuff. Will do it by smaller chunks AB, LF
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.5 KB

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

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

Download in other formats:

Original Format