Context Navigation

guide_loc_mod.F90 @ 3667

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

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

Download in other formats:

Original Format