Context Navigation

guide_loc_mod.F90 @ 2145

Last change on this file since 2145 was 2134, checked in by lguez, 10 years ago

In nudging procedures, replaced explicit Euler integration of nudged
fields by exact integration. This does not change anything if
guide_add is true, but it changes the value of alpha if guide_add is
false. We could have taken into account the variation of the nudging
field during a nudging time step. This would be a small correction. We
choose not to take it into account for the time being. Also, we add a
restriction on zonal nudging: we allow it only for a grid which is
regular in longitude. It does not seem to make sense otherwise and the
exact integration would take more programming for an irregular grid.

In the sequential version, copying the parallel versions, set
iguide_int to 1 when the input value is 0.

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

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

Download in other formats:

Original Format