Context Navigation

guide_loc_mod.F90 @ 2021

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

Removed unused variables pks, pk, pkf from main program unit gcm.

Encapsulated procedures exner_hyb, exner_hyb_p, exner_hyb_loc,
exner_milieu, exner_milieu_p and exner_milieu_loc into
modules. (Compulsory to allow optional arguments.)

In the procedures exner_hyb, exner_hyb_p, exner_hyb_loc, donwgraded
arguments alpha and beta to local variables. In exner_milieu,
exner_milieu_p and exner_milieu_loc, removed beta altogether. In the
six procedures exner_*, made pkf an optional argument. Made some
cosmetic modifications in order to keep the six procedures exner_* as
close as possible.

In the six procedures exner_*, removed the averaging of pks at the
poles: this is not useful because ps is already the same at all
longitudes at the poles. This modification changes the results of the
program. Motivation: had to do this for exner_hyb because we call it
from test_disvert with a few surface pressure values.

In all the procedures calling exner_*, removed the variables alpha and
beta. Also removed variables alpha and beta from module leapfrog_mod
and from module call_calfis_mod.

Removed actual argument pkf in call to exner_hyb* and exner_milieu*
from guide_interp, guide_main, iniacademic and iniacademic_loc (pkf
was not used in those procedures).

Argument workvar of startget_dyn is used only if varname is tpot or

When varname is tpot or q, the actual argument associated to

workvar in etat0_netcdf is not y. So y in etat0_netcdf is a
place-holder, never used. So we remove optional argument y in the
calls to exner_hyb and exner_milieu from etat0_netcdf.

Created procedure test_disvert, called only by etat0_netcdf. This
procedure tests the order of pressure values at half-levels and full
levels.

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

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

Download in other formats:

Original Format