Context Navigation

guide_loc_mod.F90 @ 3995

Last change on this file since 3995 was 3995, checked in by Ehouarn Millour, 3 years ago
Nudging: fixed some indexes in parallel about process domain boundaries and updated the serial nudging routine so that it matches the parallel one (they had diverged at some point). Also added an "is_master" logical in the parallel_lmdz module to ease decreasing the number of messages written to standard output. EM
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 86.2 KB

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

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

Download in other formats:

Original Format