Context Navigation

guide_loc_mod.f90 @ 5272

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

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

Download in other formats:

Original Format