Context Navigation

guide_loc_mod.F90 @ 2056

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

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

Download in other formats:

Original Format