Context Navigation

guide_p_mod.F90 @ 5169

Last change on this file since 5169 was 1222, checked in by Ehouarn Millour, 15 years ago

Changes and cleanups to enable compiling without physics
and without ioipsl.

IOIPSL related cleanups:

bibio/writehist.F encapsulate the routine (which needs IOIPSL to function)

with #ifdef IOIPSL flag.

dyn3d/abort_gcm.F, dyn3dpar/abort_gcm.F and dyn3dpar/getparam.F90: use ioipsl_getincom module when not compiling with IOIPSL library, in order to always be able to use getin() routine.
removed unused "use IOIPSL" in dyn3dpar/guide_p_mod.F90
calendar related issue: Initialize day_ref and annee_ref in iniacademic.F (i.e. when they are not read from start.nc file)

Earth-specific programs/routines/modules:
create_etat0.F, fluxstokenc.F, limit_netcdf.F, startvar.F
(versions in dyn3d and dyn3dpar)
These routines and modules, which by design and porpose are made to function with
Earth physics are encapsulated with #CPP_EARTH cpp flag.

Earth-specific instructions:

calls to qminimum (specific treatment of first 2 tracers, i.e. water) in dyn3d/caladvtrac.F, dyn3d/integrd.F, dyn3dpar/caladvtrac_p.F, dyn3dpar/integrd_p.F only if (planet_type == 'earth')

Interaction with parallel physics:

routine dyn3dpar/parallel.F90 uses "surface_data" module (which is in the physics ...) to know value of "type_ocean" . Encapsulated that with #ifdef CPP_EARTH and set to a default type_ocean="dummy" otherwise.
So far, only Earth physics are parallelized, so all the interaction between parallel dynamics and parallel physics are encapsulated with #ifdef CCP_EARTH (this way we can run parallel without any physics). The (dyn3dpar) routines which contains such interaction are: bands.F90, gr_dyn_fi_p.F, gr_fi_dyn_p.F, mod_interface_dyn_phys.F90 This should later (when improving dyn/phys interface) be encapsulated with a more general and appropriate #ifdef CPP_PHYS cpp flag.

I checked that these changes do not alter results (on the simple
32x24x11 bench) on Ciclad (seq & mpi), Brodie (seq, mpi & omp) and
Vargas (seq, mpi & omp).

Property svn:executable set to *

File size: 55.9 KB

Line
1	!
2	! $Header: /home/cvsroot/LMDZ4/libf/dyn3d/guide.F,v 1.3.4.1 2006/11/06 15:51:16 fairhead Exp $
3	!
4	MODULE guide_p_mod
5
6	!=======================================================================
7	! Auteur: F.Hourdin
8	! F. Codron 01/09
9	!=======================================================================
10
11	USE getparam
12	USE Write_Field_p
13	use netcdf, only: nf90_nowrite, nf90_open, nf90_inq_varid, nf90_close
14
15	IMPLICIT NONE
16
17	! ---------------------------------------------
18	! Declarations des cles logiques et parametres
19	! ---------------------------------------------
20	INTEGER, PRIVATE, SAVE :: iguide_read,iguide_int,iguide_sav
21	INTEGER, PRIVATE, SAVE :: nlevnc
22	LOGICAL, PRIVATE, SAVE :: guide_u,guide_v,guide_T,guide_Q,guide_P
23	LOGICAL, PRIVATE, SAVE :: guide_hr,guide_teta
24	LOGICAL, PRIVATE, SAVE :: guide_BL,guide_reg,guide_add,gamma4,guide_zon
25	LOGICAL, PRIVATE, SAVE :: guide_modele,invert_p,invert_y,ini_anal
26	LOGICAL, PRIVATE, SAVE :: guide_2D,guide_sav
27
28	REAL, PRIVATE, SAVE :: tau_min_u,tau_max_u
29	REAL, PRIVATE, SAVE :: tau_min_v,tau_max_v
30	REAL, PRIVATE, SAVE :: tau_min_T,tau_max_T
31	REAL, PRIVATE, SAVE :: tau_min_Q,tau_max_Q
32	REAL, PRIVATE, SAVE :: tau_min_P,tau_max_P
33
34	REAL, PRIVATE, SAVE :: lat_min_g,lat_max_g
35	REAL, PRIVATE, SAVE :: lon_min_g,lon_max_g
36	REAL, PRIVATE, SAVE :: tau_lon,tau_lat
37
38	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_u,alpha_v
39	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_T,alpha_Q
40	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: alpha_P,alpha_pcor
41
42	! ---------------------------------------------
43	! Variables de guidage
44	! ---------------------------------------------
45	! Variables des fichiers de guidage
46	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: unat1,unat2
47	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: vnat1,vnat2
48	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: tnat1,tnat2
49	REAL, ALLOCATABLE, DIMENSION(:,:,:), PRIVATE, SAVE :: qnat1,qnat2
50	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: psnat1,psnat2
51	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: apnc,bpnc
52	! Variables aux dimensions du modele
53	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: ugui1,ugui2
54	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: vgui1,vgui2
55	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: tgui1,tgui2
56	REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE :: qgui1,qgui2
57	REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: psgui1,psgui2
58
59	INTEGER,SAVE,PRIVATE :: ijb_u,ijb_v,ije_u,ije_v,ijn_u,ijn_v
60	INTEGER,SAVE,PRIVATE :: jjb_u,jjb_v,jje_u,jje_v,jjn_u,jjn_v
61
62
63	CONTAINS
64	!=======================================================================
65
66	SUBROUTINE guide_init
67
68	IMPLICIT NONE
69
70	INCLUDE "dimensions.h"
71	INCLUDE "paramet.h"
72	INCLUDE "netcdf.inc"
73	INCLUDE "control.h"
74
75	INTEGER :: error,ncidpl,rid,rcod
76	CHARACTER (len = 80) :: abort_message
77	CHARACTER (len = 20) :: modname = 'guide_init'
78
79	! ---------------------------------------------
80	! Lecture des parametres:
81	! ---------------------------------------------
82	! Variables guidees
83	CALL getpar('guide_u',.true.,guide_u,'guidage de u')
84	CALL getpar('guide_v',.true.,guide_v,'guidage de v')
85	CALL getpar('guide_T',.true.,guide_T,'guidage de T')
86	CALL getpar('guide_P',.true.,guide_P,'guidage de P')
87	CALL getpar('guide_Q',.true.,guide_Q,'guidage de Q')
88	CALL getpar('guide_hr',.true.,guide_hr,'guidage de Q par H.R')
89	CALL getpar('guide_teta',.false.,guide_teta,'guidage de T par Teta')
90
91	CALL getpar('guide_add',.false.,guide_add,'for�age constant?')
92	CALL getpar('guide_zon',.false.,guide_zon,'guidage moy zonale')
93
94	! Constantes de rappel. Unite : fraction de jour
95	CALL getpar('tau_min_u',0.02,tau_min_u,'Cste de rappel min, u')
96	CALL getpar('tau_max_u', 10.,tau_max_u,'Cste de rappel max, u')
97	CALL getpar('tau_min_v',0.02,tau_min_v,'Cste de rappel min, v')
98	CALL getpar('tau_max_v', 10.,tau_max_v,'Cste de rappel max, v')
99	CALL getpar('tau_min_T',0.02,tau_min_T,'Cste de rappel min, T')
100	CALL getpar('tau_max_T', 10.,tau_max_T,'Cste de rappel max, T')
101	CALL getpar('tau_min_Q',0.02,tau_min_Q,'Cste de rappel min, Q')
102	CALL getpar('tau_max_Q', 10.,tau_max_Q,'Cste de rappel max, Q')
103	CALL getpar('tau_min_P',0.02,tau_min_P,'Cste de rappel min, P')
104	CALL getpar('tau_max_P', 10.,tau_max_P,'Cste de rappel max, P')
105	CALL getpar('gamma4',.false.,gamma4,'Zone sans rappel elargie')
106	CALL getpar('guide_BL',.true.,guide_BL,'guidage dans C.Lim')
107
108	! Sauvegarde du for�age
109	CALL getpar('guide_sav',.false.,guide_sav,'sauvegarde guidage')
110	CALL getpar('iguide_sav',4,iguide_sav,'freq. sauvegarde guidage')
111	! frequences f>0: fx/jour; f<0: tous les f jours; f=0: 1 seule fois.
112	IF (iguide_sav.GT.0) THEN
113	iguide_sav=day_step/iguide_sav
114	ELSE
115	iguide_sav=day_step*iguide_sav
116	ENDIF
117
118	! Guidage regional seulement (sinon constant ou suivant le zoom)
119	CALL getpar('guide_reg',.false.,guide_reg,'guidage regional')
120	CALL getpar('lat_min_g',-90.,lat_min_g,'Latitude mini guidage ')
121	CALL getpar('lat_max_g', 90.,lat_max_g,'Latitude maxi guidage ')
122	CALL getpar('lon_min_g',-180.,lon_min_g,'longitude mini guidage ')
123	CALL getpar('lon_max_g', 180.,lon_max_g,'longitude maxi guidage ')
124	CALL getpar('tau_lat', 5.,tau_lat,'raideur lat guide regional ')
125	CALL getpar('tau_lon', 5.,tau_lon,'raideur lon guide regional ')
126
127	! Parametres pour lecture des fichiers
128	CALL getpar('iguide_read',4,iguide_read,'freq. lecture guidage')
129	CALL getpar('iguide_int',4,iguide_int,'freq. lecture guidage')
130	IF (iguide_int.GT.0) THEN
131	iguide_int=day_step/iguide_int
132	ELSE
133	iguide_int=day_step*iguide_int
134	ENDIF
135	CALL getpar('guide_modele',.false.,guide_modele,'guidage niveaux modele')
136	CALL getpar('ini_anal',.false.,ini_anal,'Etat initial = analyse')
137	CALL getpar('guide_invertp',.true.,invert_p,'niveaux p inverses')
138	CALL getpar('guide_inverty',.true.,invert_y,'inversion N-S')
139	CALL getpar('guide_2D',.false.,guide_2D,'fichier guidage lat-P')
140
141	! ---------------------------------------------
142	! Determination du nombre de niveaux verticaux
143	! des fichiers guidage
144	! ---------------------------------------------
145	ncidpl=-99
146	if (guide_modele) then
147	if (ncidpl.eq.-99) rcod=nf90_open('apbp.nc',Nf90_NOWRITe, ncidpl)
148	else
149	if (guide_u) then
150	if (ncidpl.eq.-99) rcod=nf90_open('u.nc',Nf90_NOWRITe,ncidpl)
151	elseif (guide_v) then
152	if (ncidpl.eq.-99) rcod=nf90_open('v.nc',nf90_nowrite,ncidpl)
153	elseif (guide_T) then
154	if (ncidpl.eq.-99) rcod=nf90_open('T.nc',nf90_nowrite,ncidpl)
155	elseif (guide_Q) then
156	if (ncidpl.eq.-99) rcod=nf90_open('hur.nc',nf90_nowrite, ncidpl)
157	endif
158	endif
159	error=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
160	IF (error.NE.NF_NOERR) error=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
161	IF (error.NE.NF_NOERR) THEN
162	print *,'Guide: probleme lecture niveaux pression'
163	CALL abort_gcm(modname,abort_message,1)
164	ENDIF
165	error=NF_INQ_DIMLEN(ncidpl,rid,nlevnc)
166	print *,'Guide: nombre niveaux vert. nlevnc', nlevnc
167	rcod = nf90_close(ncidpl)
168
169	! ---------------------------------------------
170	! Allocation des variables
171	! ---------------------------------------------
172	abort_message='pb in allocation guide'
173
174	ALLOCATE(apnc(nlevnc), stat = error)
175	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
176	ALLOCATE(bpnc(nlevnc), stat = error)
177	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
178	apnc=0.;bpnc=0.
179
180	ALLOCATE(alpha_pcor(llm), stat = error)
181	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
182	ALLOCATE(alpha_u(ip1jmp1), stat = error)
183	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
184	ALLOCATE(alpha_v(ip1jm), stat = error)
185	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
186	ALLOCATE(alpha_T(ip1jmp1), stat = error)
187	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
188	ALLOCATE(alpha_Q(ip1jmp1), stat = error)
189	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
190	ALLOCATE(alpha_P(ip1jmp1), stat = error)
191	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
192	alpha_u=0.;alpha_v=0;alpha_T=0;alpha_Q=0;alpha_P=0
193
194	IF (guide_u) THEN
195	ALLOCATE(unat1(iip1,jjp1,nlevnc), stat = error)
196	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
197	ALLOCATE(ugui1(ip1jmp1,llm), stat = error)
198	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
199	ALLOCATE(unat2(iip1,jjp1,nlevnc), stat = error)
200	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
201	ALLOCATE(ugui2(ip1jmp1,llm), stat = error)
202	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
203	unat1=0.;unat2=0.;ugui1=0.;ugui2=0.
204	ENDIF
205
206	IF (guide_T) THEN
207	ALLOCATE(tnat1(iip1,jjp1,nlevnc), stat = error)
208	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
209	ALLOCATE(tgui1(ip1jmp1,llm), stat = error)
210	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
211	ALLOCATE(tnat2(iip1,jjp1,nlevnc), stat = error)
212	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
213	ALLOCATE(tgui2(ip1jmp1,llm), stat = error)
214	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
215	tnat1=0.;tnat2=0.;tgui1=0.;tgui2=0.
216	ENDIF
217
218	IF (guide_Q) THEN
219	ALLOCATE(qnat1(iip1,jjp1,nlevnc), stat = error)
220	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
221	ALLOCATE(qgui1(ip1jmp1,llm), stat = error)
222	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
223	ALLOCATE(qnat2(iip1,jjp1,nlevnc), stat = error)
224	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
225	ALLOCATE(qgui2(ip1jmp1,llm), stat = error)
226	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
227	qnat1=0.;qnat2=0.;qgui1=0.;qgui2=0.
228	ENDIF
229
230	IF (guide_v) THEN
231	ALLOCATE(vnat1(iip1,jjm,nlevnc), stat = error)
232	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
233	ALLOCATE(vgui1(ip1jm,llm), stat = error)
234	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
235	ALLOCATE(vnat2(iip1,jjm,nlevnc), stat = error)
236	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
237	ALLOCATE(vgui2(ip1jm,llm), stat = error)
238	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
239	vnat1=0.;vnat2=0.;vgui1=0.;vgui2=0.
240	ENDIF
241
242	IF (guide_P.OR.guide_modele) THEN
243	ALLOCATE(psnat1(iip1,jjp1), stat = error)
244	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
245	ALLOCATE(psnat2(iip1,jjp1), stat = error)
246	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
247	psnat1=0.;psnat2=0.;
248	ENDIF
249	IF (guide_P) THEN
250	ALLOCATE(psgui2(ip1jmp1), stat = error)
251	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
252	ALLOCATE(psgui1(ip1jmp1), stat = error)
253	IF (error /= 0) CALL abort_gcm(modname,abort_message,1)
254	psgui1=0.;psgui2=0.
255	ENDIF
256
257	! ---------------------------------------------
258	! Lecture du premier etat de guidage.
259	! ---------------------------------------------
260	IF (guide_2D) THEN
261	CALL guide_read2D(1)
262	ELSE
263	CALL guide_read(1)
264	ENDIF
265	IF (guide_v) vnat1=vnat2
266	IF (guide_u) unat1=unat2
267	IF (guide_T) tnat1=tnat2
268	IF (guide_Q) qnat1=qnat2
269	IF (guide_P.OR.guide_modele) psnat1=psnat2
270
271	END SUBROUTINE guide_init
272
273	!=======================================================================
274	SUBROUTINE guide_main(itau,ucov,vcov,teta,q,masse,ps)
275	use parallel
276
277	IMPLICIT NONE
278
279	INCLUDE "dimensions.h"
280	INCLUDE "paramet.h"
281	INCLUDE "control.h"
282	INCLUDE "comconst.h"
283	INCLUDE "comvert.h"
284
285	! Variables entree
286	INTEGER, INTENT(IN) :: itau !pas de temps
287	REAL, DIMENSION (ip1jmp1,llm), INTENT(INOUT) :: ucov,teta,q,masse
288	REAL, DIMENSION (ip1jm,llm), INTENT(INOUT) :: vcov
289	REAL, DIMENSION (ip1jmp1), INTENT(INOUT) :: ps
290
291	! Variables locales
292	LOGICAL, SAVE :: first=.TRUE.
293	LOGICAL :: f_out ! sortie guidage
294	REAL, DIMENSION (ip1jmp1,llm) :: f_add ! var aux: champ de guidage
295	REAL, DIMENSION (ip1jmp1,llm) :: p ! besoin si guide_P
296	! Compteurs temps:
297	INTEGER, SAVE :: step_rea,count_no_rea,itau_test ! lecture guidage
298	REAL :: ditau, dday_step
299	REAL :: tau,reste ! position entre 2 etats de guidage
300	REAL, SAVE :: factt ! pas de temps en fraction de jour
301
302	INTEGER :: l
303
304	ijb_u=ij_begin ; ije_u=ij_end ; ijn_u=ije_u-ijb_u+1
305	jjb_u=jj_begin ; jje_u=jj_end ; jjn_u=jje_u-jjb_u+1
306	ijb_v=ij_begin ; ije_v=ij_end ; ijn_v=ije_v-ijb_v+1
307	jjb_v=jj_begin ; jje_v=jj_end ; jjn_v=jje_v-jjb_v+1
308	IF (pole_sud) THEN
309	ije_v=ij_end-iip1
310	jje_v=jj_end-1
311	ijn_v=ije_v-ijb_v+1
312	jjn_v=jje_v-jjb_v+1
313	ENDIF
314
315
316
317	PRINT *,'---> on rentre dans guide_main'
318	! CALL AllGather_Field(ucov,ip1jmp1,llm)
319	! CALL AllGather_Field(vcov,ip1jm,llm)
320	! CALL AllGather_Field(teta,ip1jmp1,llm)
321	! CALL AllGather_Field(ps,ip1jmp1,1)
322	! CALL AllGather_Field(q,ip1jmp1,llm)
323
324	!-----------------------------------------------------------------------
325	! Initialisations au premier passage
326	!-----------------------------------------------------------------------
327
328	IF (first) THEN
329	first=.FALSE.
330	CALL guide_init
331	itau_test=1001
332	step_rea=1
333	count_no_rea=0
334	! Calcul des constantes de rappel
335	factt=dtvr*iperiod/daysec
336	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
337	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
338	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
339	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
340	call tau2alpha(1,iip1,jjp1,factt,tau_min_Q,tau_max_Q,alpha_Q)
341	! correction de rappel dans couche limite
342	if (guide_BL) then
343	alpha_pcor(:)=1.
344	else
345	do l=1,llm
346	alpha_pcor(l)=(1.+tanh((0.85-presnivs(l)/preff)/0.05))/2.
347	enddo
348	endif
349	! ini_anal: etat initial egal au guidage
350	IF (ini_anal) THEN
351	CALL guide_interp(ps,teta)
352	IF (guide_u) ucov(ijb_u:ije_u,:)=ugui2(ijb_u:ije_u,:)
353	IF (guide_v) vcov(ijb_v:ije_v,:)=ugui2(ijb_v:ije_v,:)
354	IF (guide_T) teta(ijb_u:ije_u,:)=tgui2(ijb_u:ije_u,:)
355	IF (guide_Q) q(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)
356	IF (guide_P) THEN
357	ps(ijb_u:ije_u)=psgui2(ijb_u:ije_u)
358	CALL pression_p(ip1jmp1,ap,bp,ps,p)
359	CALL massdair_p(p,masse)
360	ENDIF
361	RETURN
362	ENDIF
363	! Verification structure guidage
364	IF (guide_u) THEN
365	CALL writefield_p('unat',unat1)
366	CALL writefield_p('ucov',RESHAPE(ucov,(/iip1,jjp1,llm/)))
367	ENDIF
368	IF (guide_T) THEN
369	CALL writefield_p('tnat',tnat1)
370	CALL writefield_p('teta',RESHAPE(teta,(/iip1,jjp1,llm/)))
371	ENDIF
372
373	ENDIF !first
374
375	!-----------------------------------------------------------------------
376	! Lecture des fichiers de guidage ?
377	!-----------------------------------------------------------------------
378	IF (iguide_read.NE.0) THEN
379	ditau=real(itau)
380	dday_step=real(day_step)
381	IF (iguide_read.LT.0) THEN
382	tau=ditau/dday_step/FLOAT(iguide_read)
383	ELSE
384	tau=FLOAT(iguide_read)*ditau/dday_step
385	ENDIF
386	reste=tau-AINT(tau)
387	IF (reste.EQ.0.) THEN
388	IF (itau_test.EQ.itau) THEN
389	write(,)'deuxieme passage de advreel a itau=',itau
390	stop
391	ELSE
392	IF (guide_v) vnat1(jjb_v:jje_v,:,:)=vnat2(jjb_v:jje_v,:,:)
393	IF (guide_u) unat1(jjb_u:jje_u,:,:)=unat2(jjb_u:jje_u,:,:)
394	IF (guide_T) tnat1(jjb_u:jje_u,:,:)=tnat2(jjb_u:jje_u,:,:)
395	IF (guide_Q) qnat1(jjb_u:jje_u,:,:)=qnat2(jjb_u:jje_u,:,:)
396	IF (guide_P.OR.guide_modele) psnat1(jjb_u:jje_u,:)=psnat2(jjb_u:jje_u,:)
397	step_rea=step_rea+1
398	itau_test=itau
399	print*,'Lecture fichiers guidage, pas ',step_rea, &
400	'apres ',count_no_rea,' non lectures'
401	IF (guide_2D) THEN
402	CALL guide_read2D(step_rea)
403	ELSE
404	CALL guide_read(step_rea)
405	ENDIF
406	count_no_rea=0
407	ENDIF
408	ELSE
409	count_no_rea=count_no_rea+1
410
411	ENDIF
412	ENDIF !iguide_read=0
413
414	!-----------------------------------------------------------------------
415	! Interpolation et conversion des champs de guidage
416	!-----------------------------------------------------------------------
417	IF (MOD(itau,iguide_int).EQ.0) THEN
418	CALL guide_interp(ps,teta)
419	ENDIF
420	! Repartition entre 2 etats de guidage
421	IF (iguide_read.NE.0) THEN
422	tau=reste
423	ELSE
424	tau=1.
425	ENDIF
426
427	!-----------------------------------------------------------------------
428	! Ajout des champs de guidage
429	!-----------------------------------------------------------------------
430	! Sauvegarde du guidage?
431	f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)
432	IF (f_out) CALL guide_out("S",jjp1,1,ps)
433
434	if (guide_u) then
435	if (guide_add) then
436	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)
437	else
438	f_add(ijb_u:ije_u,:)=(1.-tau)ugui1(ijb_u:ije_u,:)+tauugui2(ijb_u:ije_u,:)-ucov(ijb_u:ije_u,:)
439	endif
440
441	if (guide_zon) CALL guide_zonave(1,jjp1,llm,f_add)
442	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_u)
443	IF (f_out) CALL guide_out("U",jjp1,llm,f_add/factt)
444	ucov(ijb_u:ije_u,:)=ucov(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
445	endif
446
447	if (guide_T) then
448	if (guide_add) then
449	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)
450	else
451	f_add(ijb_u:ije_u,:)=(1.-tau)tgui1(ijb_u:ije_u,:)+tautgui2(ijb_u:ije_u,:)-teta(ijb_u:ije_u,:)
452	endif
453	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
454	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_T)
455	IF (f_out) CALL guide_out("T",jjp1,llm,f_add/factt)
456	teta(ijb_u:ije_u,:)=teta(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
457	endif
458
459	if (guide_P) then
460	if (guide_add) then
461	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)
462	else
463	f_add(ijb_u:ije_u,1)=(1.-tau)psgui1(ijb_u:ije_u)+taupsgui2(ijb_u:ije_u)-ps(ijb_u:ije_u)
464	endif
465	if (guide_zon) CALL guide_zonave(2,jjp1,1,f_add(1:ip1jmp1,1))
466	CALL guide_addfield(ip1jmp1,1,f_add(1:ip1jmp1,1),alpha_P)
467	IF (f_out) CALL guide_out("P",jjp1,1,f_add(1:ip1jmp1,1)/factt)
468	ps(ijb_u:ije_u)=ps(ijb_u:ije_u)+f_add(ijb_u:ije_u,1)
469	CALL pression_p(ip1jmp1,ap,bp,ps,p)
470	CALL massdair_p(p,masse)
471	endif
472
473	if (guide_Q) then
474	if (guide_add) then
475	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)
476	else
477	f_add(ijb_u:ije_u,:)=(1.-tau)qgui1(ijb_u:ije_u,:)+tauqgui2(ijb_u:ije_u,:)-q(ijb_u:ije_u,:)
478	endif
479	if (guide_zon) CALL guide_zonave(2,jjp1,llm,f_add)
480	CALL guide_addfield(ip1jmp1,llm,f_add,alpha_Q)
481	IF (f_out) CALL guide_out("Q",jjp1,llm,f_add/factt)
482	q(ijb_u:ije_u,:)=q(ijb_u:ije_u,:)+f_add(ijb_u:ije_u,:)
483	endif
484
485	if (guide_v) then
486	if (guide_add) then
487	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)
488	else
489	f_add(ijb_v:ije_v,:)=(1.-tau)vgui1(ijb_v:ije_v,:)+tauvgui2(ijb_v:ije_v,:)-vcov(ijb_v:ije_v,:)
490	endif
491
492	if (guide_zon) CALL guide_zonave(2,jjm,llm,f_add(1:ip1jm,:))
493	CALL guide_addfield(ip1jm,llm,f_add(1:ip1jm,:),alpha_v)
494	IF (f_out) CALL guide_out("V",jjm,llm,f_add(1:ip1jm,:)/factt)
495	vcov(ijb_v:ije_v,:)=vcov(ijb_v:ije_v,:)+f_add(ijb_v:ije_v,:)
496	endif
497
498	END SUBROUTINE guide_main
499
500	!=======================================================================
501	SUBROUTINE guide_addfield(hsize,vsize,field,alpha)
502	! field1=afield1+alphafield2
503
504	IMPLICIT NONE
505	INCLUDE "dimensions.h"
506	INCLUDE "paramet.h"
507
508	! input variables
509	INTEGER, INTENT(IN) :: hsize
510	INTEGER, INTENT(IN) :: vsize
511	REAL, DIMENSION(hsize), INTENT(IN) :: alpha
512	REAL, DIMENSION(hsize,vsize), INTENT(INOUT) :: field
513
514	! Local variables
515	INTEGER :: l
516
517	IF (hsize==ip1jm) THEN
518	do l=1,vsize
519	field(ijb_v:ije_v,l)=alpha(ijb_v:ije_v)field(ijb_v:ije_v,l)alpha_pcor(l)
520	enddo
521	ELSE
522	do l=1,vsize
523	field(ijb_u:ije_u,l)=alpha(ijb_u:ije_u)field(ijb_u:ije_u,l)alpha_pcor(l)
524	enddo
525	ENDIF
526
527	END SUBROUTINE guide_addfield
528
529	!=======================================================================
530	SUBROUTINE guide_zonave(typ,hsize,vsize,field)
531
532	IMPLICIT NONE
533
534	INCLUDE "dimensions.h"
535	INCLUDE "paramet.h"
536	INCLUDE "comgeom.h"
537	INCLUDE "comconst.h"
538
539	! input/output variables
540	INTEGER, INTENT(IN) :: typ
541	INTEGER, INTENT(IN) :: vsize
542	INTEGER, INTENT(IN) :: hsize
543	REAL, DIMENSION(hsize*iip1,vsize), INTENT(INOUT) :: field
544
545	! Local variables
546	LOGICAL, SAVE :: first=.TRUE.
547	INTEGER, DIMENSION (2), SAVE :: imin, imax ! averaging domain
548	INTEGER :: i,j,l,ij
549	REAL, DIMENSION (iip1) :: lond ! longitude in Deg.
550	REAL, DIMENSION (hsize,vsize):: fieldm ! zon-averaged field
551
552	IF (first) THEN
553	first=.FALSE.
554	!Compute domain for averaging
555	lond=rlonu*180./pi
556	imin(1)=1;imax(1)=iip1;
557	imin(2)=1;imax(2)=iip1;
558	IF (guide_reg) THEN
559	DO i=1,iim
560	IF (lond(i).LT.lon_min_g) imin(1)=i
561	IF (lond(i).LE.lon_max_g) imax(1)=i
562	ENDDO
563	lond=rlonv*180./pi
564	DO i=1,iim
565	IF (lond(i).LT.lon_min_g) imin(2)=i
566	IF (lond(i).LE.lon_max_g) imax(2)=i
567	ENDDO
568	ENDIF
569	ENDIF
570
571	fieldm=0.
572
573	IF (hsize==jjm) THEN
574	DO l=1,vsize
575	! Compute zonal average
576	DO j=jjb_v,jje_v
577	DO i=imin(typ),imax(typ)
578	ij=(j-1)*iip1+i
579	fieldm(j,l)=fieldm(j,l)+field(ij,l)
580	ENDDO
581	ENDDO
582	fieldm(:,l)=fieldm(:,l)/FLOAT(imax(typ)-imin(typ)+1)
583	! Compute forcing
584	DO j=jjb_v,jje_v
585	DO i=1,iip1
586	ij=(j-1)*iip1+i
587	field(ij,l)=fieldm(j,l)
588	ENDDO
589	ENDDO
590	ENDDO
591	ELSE
592	DO l=1,vsize
593	! Compute zonal average
594	DO j=jjb_v,jje_v
595	DO i=imin(typ),imax(typ)
596	ij=(j-1)*iip1+i
597	fieldm(j,l)=fieldm(j,l)+field(ij,l)
598	ENDDO
599	ENDDO
600	fieldm(:,l)=fieldm(:,l)/FLOAT(imax(typ)-imin(typ)+1)
601	! Compute forcing
602	DO j=jjb_u,jje_u
603	DO i=1,iip1
604	ij=(j-1)*iip1+i
605	field(ij,l)=fieldm(j,l)
606	ENDDO
607	ENDDO
608	ENDDO
609	ENDIF
610
611	END SUBROUTINE guide_zonave
612
613	!=======================================================================
614	SUBROUTINE guide_interp(psi,teta)
615	USE parallel
616	USE mod_hallo
617	USE Bands
618	IMPLICIT NONE
619
620	include "dimensions.h"
621	include "paramet.h"
622	include "comvert.h"
623	include "comgeom2.h"
624	include "comconst.h"
625
626	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
627	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
628
629	LOGICAL, SAVE :: first=.TRUE.
630	! Variables pour niveaux pression:
631	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
632	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
633	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
634	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
635	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
636	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
637	REAL, DIMENSION (iip1,jjm,llm) :: pbary
638	! Variables pour fonction Exner (P milieu couche)
639	REAL, DIMENSION (iip1,jjp1,llm) :: pk, pkf
640	REAL, DIMENSION (iip1,jjp1,llm) :: alpha, beta
641	REAL, DIMENSION (iip1,jjp1) :: pks
642	REAL :: prefkap,unskap
643	! Pression de vapeur saturante
644	REAL, DIMENSION (ip1jmp1,llm) :: qsat
645	!Variables intermediaires interpolation
646	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
647	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
648
649	INTEGER :: i,j,l,ij
650	TYPE(Request) :: Req
651
652	print *,'Guide: conversion variables guidage'
653	! -----------------------------------------------------------------
654	! Calcul des niveaux de pression champs guidage
655	! -----------------------------------------------------------------
656	if (guide_modele) then
657	do i=1,iip1
658	do j=jjb_u,jje_u
659	do l=1,nlevnc
660	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
661	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
662	enddo
663	enddo
664	enddo
665	else
666	do i=1,iip1
667	do j=jjb_u,jje_u
668	do l=1,nlevnc
669	plnc2(i,j,l)=apnc(l)
670	plnc1(i,j,l)=apnc(l)
671	enddo
672	enddo
673	enddo
674
675	endif
676	if (first) then
677	first=.FALSE.
678	print*,'Guide: verification ordre niveaux verticaux'
679	print*,'LMDZ :'
680	do l=1,llm
681	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
682	+psi(1,jje_u)*(bp(l)+bp(l+1))/2.
683	enddo
684	print*,'Fichiers guidage'
685	do l=1,nlevnc
686	print*,'PL(',l,')=',plnc2(1,jjb_u,l)
687	enddo
688	print *,'inversion de l''ordre: invert_p=',invert_p
689	if (guide_u) then
690	do l=1,nlevnc
691	print*,'U(',l,')=',unat2(1,jjb_u,l)
692	enddo
693	endif
694	if (guide_T) then
695	do l=1,nlevnc
696	print*,'T(',l,')=',tnat2(1,jjb_u,l)
697	enddo
698	endif
699	endif
700
701	! -----------------------------------------------------------------
702	! Calcul niveaux pression modele
703	! -----------------------------------------------------------------
704	CALL pression_p( ip1jmp1, ap, bp, psi, p )
705	CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
706
707	! .... Calcul de pls , pression au milieu des couches ,en Pascals
708	unskap=1./kappa
709	prefkap = preff ** kappa
710	DO l = 1, llm
711	DO j=jjb_u,jje_u
712	DO i =1, iip1
713	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
714	ENDDO
715	ENDDO
716	ENDDO
717
718	! calcul des pressions pour les grilles u et v
719	do l=1,llm
720	do j=jjb_u,jje_u
721	do i=1,iip1
722	pext(i,j,l)=pls(i,j,l)*aire(i,j)
723	enddo
724	enddo
725	enddo
726
727	CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
728	CALL SendRequest(Req)
729	CALL WaitRequest(Req)
730
731	call massbar_p(pext, pbarx, pbary )
732	do l=1,llm
733	do j=jjb_u,jje_u
734	do i=1,iip1
735	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
736	plsnc(i,j,l)=pls(i,j,l)
737	enddo
738	enddo
739	enddo
740	do l=1,llm
741	do j=jjb_v,jje_v
742	do i=1,iip1
743	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
744	enddo
745	enddo
746	enddo
747
748	! -----------------------------------------------------------------
749	! Interpolation champs guidage sur niveaux modele (+inversion N/S)
750	! Conversion en variables gcm (ucov, vcov...)
751	! -----------------------------------------------------------------
752	if (guide_P) then
753	do j=jjb_u,jje_u
754	do i=1,iim
755	ij=(j-1)*iip1+i
756	psgui1(ij)=psnat1(i,j)
757	psgui2(ij)=psnat2(i,j)
758	enddo
759	psgui1(iip1*j)=psnat1(1,j)
760	psgui2(iip1*j)=psnat2(1,j)
761	enddo
762	endif
763
764	IF (guide_u) THEN
765	CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
766	plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
767	CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
768	plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
769
770	do l=1,llm
771	do j=jjb_u,jje_u
772	do i=1,iim
773	ij=(j-1)*iip1+i
774	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
775	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
776	enddo
777	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
778	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
779	enddo
780	do i=1,iip1
781	ugui1(i,l)=0.
782	ugui1(ip1jm+i,l)=0.
783	ugui2(i,l)=0.
784	ugui2(ip1jm+i,l)=0.
785	enddo
786	enddo
787	ENDIF
788
789	IF (guide_T) THEN
790	CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
791	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
792	CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
793	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
794
795	do l=1,llm
796	do j=jjb_u,jje_u
797	IF (guide_teta) THEN
798	do i=1,iim
799	ij=(j-1)*iip1+i
800	tgui1(ij,l)=zu1(i,j,l)
801	tgui2(ij,l)=zu2(i,j,l)
802	enddo
803	ELSE
804	do i=1,iim
805	ij=(j-1)*iip1+i
806	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
807	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
808	enddo
809	ENDIF
810	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
811	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
812	enddo
813	do i=1,iip1
814	tgui1(i,l)=tgui1(1,l)
815	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
816	tgui2(i,l)=tgui2(1,l)
817	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
818	enddo
819	enddo
820	ENDIF
821
822	IF (guide_v) THEN
823
824	CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm, &
825	plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
826	CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm, &
827	plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
828
829	do l=1,llm
830	do j=jjb_v,jje_v
831	do i=1,iim
832	ij=(j-1)*iip1+i
833	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
834	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
835	enddo
836	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
837	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
838	enddo
839	enddo
840	ENDIF
841
842	IF (guide_Q) THEN
843	! On suppose qu'on a la bonne variable dans le fichier de guidage:
844	! Hum.Rel si guide_hr, Hum.Spec. sinon.
845	CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
846	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
847	CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
848	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
849
850	do l=1,llm
851	do j=jjb_u,jjb_v
852	do i=1,iim
853	ij=(j-1)*iip1+i
854	qgui1(ij,l)=zu1(i,j,l)
855	qgui2(ij,l)=zu2(i,j,l)
856	enddo
857	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
858	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
859	enddo
860	do i=1,iip1
861	qgui1(i,l)=qgui1(1,l)
862	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
863	qgui2(i,l)=qgui2(1,l)
864	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
865	enddo
866	enddo
867	IF (guide_hr) THEN
868	CALL q_sat(iip1jjn_ullm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp, &
869	plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
870	qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01 !hum. rel. en %
871	qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01
872	ENDIF
873	ENDIF
874
875	END SUBROUTINE guide_interp
876
877	!=======================================================================
878	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
879
880	! Calcul des constantes de rappel alpha (=1/tau)
881
882	implicit none
883
884	include "dimensions.h"
885	include "paramet.h"
886	include "comconst.h"
887	include "comgeom2.h"
888	include "serre.h"
889
890	! input arguments :
891	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
892	INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
893	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
894	REAL, INTENT(IN) :: taumin,taumax
895	! output arguments:
896	REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
897
898	! local variables:
899	LOGICAL, SAVE :: first=.TRUE.
900	REAL, SAVE :: gamma,dxdy_min,dxdy_max
901	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
902	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
903	REAL, DIMENSION (iip1,jjm) :: dxdyv
904	real dxdy_
905	real zlat,zlon
906	real alphamin,alphamax,xi
907	integer i,j,ilon,ilat
908
909
910	alphamin=factt/taumax
911	alphamax=factt/taumin
912	IF (guide_reg.OR.guide_add) THEN
913	alpha=alphamax
914	!-----------------------------------------------------------------------
915	! guide_reg: alpha=alpha_min dans region, 0. sinon.
916	!-----------------------------------------------------------------------
917	IF (guide_reg) THEN
918	do j=1,pjm
919	do i=1,pim
920	if (typ.eq.2) then
921	zlat=rlatu(j)*180./pi
922	zlon=rlonu(i)*180./pi
923	elseif (typ.eq.1) then
924	zlat=rlatu(j)*180./pi
925	zlon=rlonv(i)*180./pi
926	elseif (typ.eq.3) then
927	zlat=rlatv(j)*180./pi
928	zlon=rlonv(i)*180./pi
929	endif
930	alpha(i,j)=alphamax/16.* &
931	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
932	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
933	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
934	(1.+tanh((lon_max_g-zlon)/tau_lon))
935	enddo
936	enddo
937	ENDIF
938	ELSE
939	!-----------------------------------------------------------------------
940	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
941	!-----------------------------------------------------------------------
942	!Calcul de l'aire des mailles
943	do j=2,jjm
944	do i=2,iip1
945	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
946	enddo
947	zdx(1,j)=zdx(iip1,j)
948	enddo
949	do j=2,jjm
950	do i=1,iip1
951	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
952	enddo
953	enddo
954	do i=1,iip1
955	zdx(i,1)=zdx(i,2)
956	zdx(i,jjp1)=zdx(i,jjm)
957	zdy(i,1)=zdy(i,2)
958	zdy(i,jjp1)=zdy(i,jjm)
959	enddo
960	do j=1,jjp1
961	do i=1,iip1
962	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
963	enddo
964	enddo
965	IF (typ.EQ.2) THEN
966	do j=1,jjp1
967	do i=1,iim
968	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
969	enddo
970	dxdyu(iip1,j)=dxdyu(1,j)
971	enddo
972	ENDIF
973	IF (typ.EQ.3) THEN
974	do j=1,jjm
975	do i=1,iip1
976	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
977	enddo
978	enddo
979	ENDIF
980	! Premier appel: calcul des aires min et max et de gamma.
981	IF (first) THEN
982	first=.FALSE.
983	! coordonnees du centre du zoom
984	CALL coordij(clon,clat,ilon,ilat)
985	! aire de la maille au centre du zoom
986	dxdy_min=dxdys(ilon,ilat)
987	! dxdy maximale de la maille
988	dxdy_max=0.
989	do j=1,jjp1
990	do i=1,iip1
991	dxdy_max=max(dxdy_max,dxdys(i,j))
992	enddo
993	enddo
994	! Calcul de gamma
995	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
996	print*,'ATTENTION modele peu zoome'
997	print*,'ATTENTION on prend une constante de guidage cste'
998	gamma=0.
999	else
1000	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
1001	print*,'gamma=',gamma
1002	if (gamma.lt.1.e-5) then
1003	print*,'gamma =',gamma,'<1e-5'
1004	stop
1005	endif
1006	gamma=log(0.5)/log(gamma)
1007	if (gamma4) then
1008	gamma=min(gamma,4.)
1009	endif
1010	print*,'gamma=',gamma
1011	endif
1012	ENDIF !first
1013
1014	do j=1,pjm
1015	do i=1,pim
1016	if (typ.eq.1) then
1017	dxdy_=dxdys(i,j)
1018	zlat=rlatu(j)*180./pi
1019	elseif (typ.eq.2) then
1020	dxdy_=dxdyu(i,j)
1021	zlat=rlatu(j)*180./pi
1022	elseif (typ.eq.3) then
1023	dxdy_=dxdyv(i,j)
1024	zlat=rlatv(j)*180./pi
1025	endif
1026	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1027	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
1028	alpha(i,j)=alphamin
1029	else
1030	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
1031	xi=min(xi,1.)
1032	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
1033	alpha(i,j)=xialphamin+(1.-xi)alphamax
1034	else
1035	alpha(i,j)=0.
1036	endif
1037	endif
1038	enddo
1039	enddo
1040	ENDIF ! guide_reg
1041
1042	END SUBROUTINE tau2alpha
1043
1044	!=======================================================================
1045	SUBROUTINE guide_read(timestep)
1046
1047	IMPLICIT NONE
1048
1049	#include "netcdf.inc"
1050	#include "dimensions.h"
1051	#include "paramet.h"
1052
1053	INTEGER, INTENT(IN) :: timestep
1054
1055	LOGICAL, SAVE :: first=.TRUE.
1056	! Identification fichiers et variables NetCDF:
1057	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1058	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1059	INTEGER :: ncidpl,varidpl,varidap,varidbp
1060	! Variables auxiliaires NetCDF:
1061	INTEGER, DIMENSION(4) :: start,count
1062	INTEGER :: status,rcode
1063
1064	! -----------------------------------------------------------------
1065	! Premier appel: initialisation de la lecture des fichiers
1066	! -----------------------------------------------------------------
1067	if (first) then
1068	ncidpl=-99
1069	print*,'Guide: ouverture des fichiers guidage '
1070	! Niveaux de pression si non constants
1071	if (guide_modele) then
1072	print *,'Lecture du guidage sur niveaux mod�le'
1073	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1074	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1075	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1076	print*,'ncidpl,varidap',ncidpl,varidap
1077	endif
1078	! Vent zonal
1079	if (guide_u) then
1080	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1081	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1082	print*,'ncidu,varidu',ncidu,varidu
1083	if (ncidpl.eq.-99) ncidpl=ncidu
1084	endif
1085	! Vent meridien
1086	if (guide_v) then
1087	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1088	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1089	print*,'ncidv,varidv',ncidv,varidv
1090	if (ncidpl.eq.-99) ncidpl=ncidv
1091	endif
1092	! Temperature
1093	if (guide_T) then
1094	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1095	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1096	print*,'ncidT,varidT',ncidt,varidt
1097	if (ncidpl.eq.-99) ncidpl=ncidt
1098	endif
1099	! Humidite
1100	if (guide_Q) then
1101	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1102	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1103	print*,'ncidQ,varidQ',ncidQ,varidQ
1104	if (ncidpl.eq.-99) ncidpl=ncidQ
1105	endif
1106	! Pression de surface
1107	if ((guide_P).OR.(guide_modele)) then
1108	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1109	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1110	print*,'ncidps,varidps',ncidps,varidps
1111	endif
1112	! Coordonnee verticale
1113	if (.not.guide_modele) then
1114	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1115	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1116	print*,'ncidpl,varidpl',ncidpl,varidpl
1117	endif
1118	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1119	if (guide_modele) then
1120	#ifdef NC_DOUBLE
1121	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1122	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1123	#else
1124	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1125	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1126	#endif
1127	else
1128	#ifdef NC_DOUBLE
1129	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1130	#else
1131	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1132	#endif
1133	apnc=apnc*100.! conversion en Pascals
1134	bpnc(:)=0.
1135	endif
1136	first=.FALSE.
1137	endif ! (first)
1138
1139	! -----------------------------------------------------------------
1140	! lecture des champs u, v, T, Q, ps
1141	! -----------------------------------------------------------------
1142
1143	! dimensions pour les champs scalaires et le vent zonal
1144	start(1)=1
1145	start(2)=1
1146	start(3)=1
1147	start(4)=timestep
1148
1149	count(1)=iip1
1150	count(2)=jjp1
1151	count(3)=nlevnc
1152	count(4)=1
1153
1154	! Vent zonal
1155	if (guide_u) then
1156	#ifdef NC_DOUBLE
1157	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
1158	#else
1159	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
1160	#endif
1161	IF (invert_y) THEN
1162	CALL invert_lat(iip1,jjp1,llm,unat2)
1163	ENDIF
1164
1165	endif
1166
1167	! Temperature
1168	if (guide_T) then
1169	#ifdef NC_DOUBLE
1170	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
1171	#else
1172	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
1173	#endif
1174	IF (invert_y) THEN
1175	CALL invert_lat(iip1,jjp1,llm,tnat2)
1176	ENDIF
1177	endif
1178
1179	! Humidite
1180	if (guide_Q) then
1181	#ifdef NC_DOUBLE
1182	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
1183	#else
1184	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
1185	#endif
1186	IF (invert_y) THEN
1187	CALL invert_lat(iip1,jjp1,llm,qnat2)
1188	ENDIF
1189
1190	endif
1191
1192	! Vent meridien
1193	if (guide_v) then
1194	count(2)=jjm
1195	#ifdef NC_DOUBLE
1196	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
1197	#else
1198	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
1199	#endif
1200	IF (invert_y) THEN
1201	CALL invert_lat(iip1,jjm,llm,vnat2)
1202	ENDIF
1203	endif
1204
1205	! Pression de surface
1206	if ((guide_P).OR.(guide_modele)) then
1207	start(3)=timestep
1208	start(4)=0
1209	count(2)=jjp1
1210	count(3)=1
1211	count(4)=0
1212	#ifdef NC_DOUBLE
1213	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
1214	#else
1215	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
1216	#endif
1217	IF (invert_y) THEN
1218	CALL invert_lat(iip1,jjp1,1,psnat2)
1219	ENDIF
1220	endif
1221
1222	END SUBROUTINE guide_read
1223
1224	!=======================================================================
1225	SUBROUTINE guide_read2D(timestep)
1226
1227	IMPLICIT NONE
1228
1229	#include "netcdf.inc"
1230	#include "dimensions.h"
1231	#include "paramet.h"
1232
1233	INTEGER, INTENT(IN) :: timestep
1234
1235	LOGICAL, SAVE :: first=.TRUE.
1236	! Identification fichiers et variables NetCDF:
1237	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1238	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1239	INTEGER :: ncidpl,varidpl,varidap,varidbp
1240	! Variables auxiliaires NetCDF:
1241	INTEGER, DIMENSION(4) :: start,count
1242	INTEGER :: status,rcode
1243	! Variables for 3D extension:
1244	REAL, DIMENSION (jjp1,llm) :: zu
1245	REAL, DIMENSION (jjm,llm) :: zv
1246	INTEGER :: i
1247
1248	! -----------------------------------------------------------------
1249	! Premier appel: initialisation de la lecture des fichiers
1250	! -----------------------------------------------------------------
1251	if (first) then
1252	ncidpl=-99
1253	print*,'Guide: ouverture des fichiers guidage '
1254	! Niveaux de pression si non constants
1255	if (guide_modele) then
1256	print *,'Lecture du guidage sur niveaux mod�le'
1257	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1258	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1259	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1260	print*,'ncidpl,varidap',ncidpl,varidap
1261	endif
1262	! Vent zonal
1263	if (guide_u) then
1264	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1265	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1266	print*,'ncidu,varidu',ncidu,varidu
1267	if (ncidpl.eq.-99) ncidpl=ncidu
1268	endif
1269	! Vent meridien
1270	if (guide_v) then
1271	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1272	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1273	print*,'ncidv,varidv',ncidv,varidv
1274	if (ncidpl.eq.-99) ncidpl=ncidv
1275	endif
1276	! Temperature
1277	if (guide_T) then
1278	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1279	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1280	print*,'ncidT,varidT',ncidt,varidt
1281	if (ncidpl.eq.-99) ncidpl=ncidt
1282	endif
1283	! Humidite
1284	if (guide_Q) then
1285	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1286	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1287	print*,'ncidQ,varidQ',ncidQ,varidQ
1288	if (ncidpl.eq.-99) ncidpl=ncidQ
1289	endif
1290	! Pression de surface
1291	if ((guide_P).OR.(guide_modele)) then
1292	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1293	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1294	print*,'ncidps,varidps',ncidps,varidps
1295	endif
1296	! Coordonnee verticale
1297	if (.not.guide_modele) then
1298	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1299	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1300	print*,'ncidpl,varidpl',ncidpl,varidpl
1301	endif
1302	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1303	if (guide_modele) then
1304	#ifdef NC_DOUBLE
1305	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1306	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1307	#else
1308	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1309	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1310	#endif
1311	else
1312	#ifdef NC_DOUBLE
1313	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1314	#else
1315	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1316	#endif
1317	apnc=apnc*100.! conversion en Pascals
1318	bpnc(:)=0.
1319	endif
1320	first=.FALSE.
1321	endif ! (first)
1322
1323	! -----------------------------------------------------------------
1324	! lecture des champs u, v, T, Q, ps
1325	! -----------------------------------------------------------------
1326
1327	! dimensions pour les champs scalaires et le vent zonal
1328	start(1)=1
1329	start(2)=1
1330	start(3)=1
1331	start(4)=timestep
1332
1333	count(1)=1
1334	count(2)=jjp1
1335	count(3)=nlevnc
1336	count(4)=1
1337
1338	! Vent zonal
1339	if (guide_u) then
1340	#ifdef NC_DOUBLE
1341	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
1342	#else
1343	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
1344	#endif
1345	DO i=1,iip1
1346	unat2(i,:,:)=zu(:,:)
1347	ENDDO
1348
1349	IF (invert_y) THEN
1350	CALL invert_lat(iip1,jjp1,llm,unat2)
1351	ENDIF
1352
1353	endif
1354
1355	! Temperature
1356	if (guide_T) then
1357	#ifdef NC_DOUBLE
1358	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
1359	#else
1360	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
1361	#endif
1362	DO i=1,iip1
1363	tnat2(i,:,:)=zu(:,:)
1364	ENDDO
1365
1366	IF (invert_y) THEN
1367	CALL invert_lat(iip1,jjp1,llm,tnat2)
1368	ENDIF
1369
1370	endif
1371
1372	! Humidite
1373	if (guide_Q) then
1374	#ifdef NC_DOUBLE
1375	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
1376	#else
1377	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
1378	#endif
1379	DO i=1,iip1
1380	qnat2(i,:,:)=zu(:,:)
1381	ENDDO
1382
1383	IF (invert_y) THEN
1384	CALL invert_lat(iip1,jjp1,llm,qnat2)
1385	ENDIF
1386
1387	endif
1388
1389	! Vent meridien
1390	if (guide_v) then
1391	count(2)=jjm
1392	#ifdef NC_DOUBLE
1393	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
1394	#else
1395	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
1396	#endif
1397	DO i=1,iip1
1398	vnat2(i,:,:)=zv(:,:)
1399	ENDDO
1400
1401	IF (invert_y) THEN
1402	CALL invert_lat(iip1,jjm,llm,vnat2)
1403	ENDIF
1404
1405	endif
1406
1407	! Pression de surface
1408	if ((guide_P).OR.(guide_modele)) then
1409	start(3)=timestep
1410	start(4)=0
1411	count(2)=jjp1
1412	count(3)=1
1413	count(4)=0
1414	#ifdef NC_DOUBLE
1415	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
1416	#else
1417	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
1418	#endif
1419	DO i=1,iip1
1420	psnat2(i,:)=zu(:,1)
1421	ENDDO
1422
1423	IF (invert_y) THEN
1424	CALL invert_lat(iip1,jjp1,1,psnat2)
1425	ENDIF
1426
1427	endif
1428
1429	END SUBROUTINE guide_read2D
1430
1431	!=======================================================================
1432	SUBROUTINE guide_out(varname,hsize,vsize,field)
1433	USE parallel
1434	IMPLICIT NONE
1435
1436	INCLUDE "dimensions.h"
1437	INCLUDE "paramet.h"
1438	INCLUDE "netcdf.inc"
1439	INCLUDE "comgeom2.h"
1440	INCLUDE "comconst.h"
1441	INCLUDE "comvert.h"
1442
1443	! Variables entree
1444	CHARACTER, INTENT(IN) :: varname
1445	INTEGER, INTENT (IN) :: hsize,vsize
1446	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1447
1448	! Variables locales
1449	INTEGER, SAVE :: timestep=0
1450	! Identites fichier netcdf
1451	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1452	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1453	INTEGER, DIMENSION (3) :: dim3
1454	INTEGER, DIMENSION (4) :: dim4,count,start
1455	INTEGER :: ierr, varid
1456
1457	CALL gather_field(field,iip1*hsize,vsize,0)
1458
1459	IF (mpi_rank /= 0) RETURN
1460
1461	print *,'Guide: output timestep',timestep,'var ',varname
1462	IF (timestep.EQ.0) THEN
1463	! ----------------------------------------------
1464	! initialisation fichier de sortie
1465	! ----------------------------------------------
1466	! Ouverture du fichier
1467	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
1468	! Definition des dimensions
1469	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1470	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1471	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1472	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1473	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1474	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1475
1476	! Creation des variables dimensions
1477	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1478	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1479	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1480	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1481	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1482	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1483	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1484
1485	ierr=NF_ENDDEF(nid)
1486
1487	! Enregistrement des variables dimensions
1488	#ifdef NC_DOUBLE
1489	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1490	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1491	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1492	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1493	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1494	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1495	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1496	#else
1497	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1498	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1499	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1500	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1501	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1502	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1503	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1504	#endif
1505	! --------------------------------------------------------------------
1506	! Cr�ation des variables sauvegard�es
1507	! --------------------------------------------------------------------
1508	ierr = NF_REDEF(nid)
1509	! Surface pressure (GCM)
1510	dim3=(/id_lonv,id_latu,id_tim/)
1511	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
1512	! Surface pressure (guidage)
1513	IF (guide_P) THEN
1514	dim3=(/id_lonv,id_latu,id_tim/)
1515	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1516	ENDIF
1517	! Zonal wind
1518	IF (guide_u) THEN
1519	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1520	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1521	ENDIF
1522	! Merid. wind
1523	IF (guide_v) THEN
1524	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1525	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1526	ENDIF
1527	! Pot. Temperature
1528	IF (guide_T) THEN
1529	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1530	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1531	ENDIF
1532	! Specific Humidity
1533	IF (guide_Q) THEN
1534	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1535	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1536	ENDIF
1537
1538	ierr = NF_ENDDEF(nid)
1539	ierr = NF_CLOSE(nid)
1540	ENDIF ! timestep=0
1541
1542	! --------------------------------------------------------------------
1543	! Enregistrement du champ
1544	! --------------------------------------------------------------------
1545	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1546
1547	SELECT CASE (varname)
1548	CASE ("S")
1549	timestep=timestep+1
1550	ierr = NF_INQ_VARID(nid,"SP",varid)
1551	start=(/1,1,timestep,0/)
1552	count=(/iip1,jjp1,1,0/)
1553	#ifdef NC_DOUBLE
1554	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1555	#else
1556	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1557	#endif
1558	CASE ("P")
1559	ierr = NF_INQ_VARID(nid,"ps",varid)
1560	start=(/1,1,timestep,0/)
1561	count=(/iip1,jjp1,1,0/)
1562	#ifdef NC_DOUBLE
1563	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1564	#else
1565	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1566	#endif
1567	CASE ("U")
1568	ierr = NF_INQ_VARID(nid,"ucov",varid)
1569	start=(/1,1,1,timestep/)
1570	count=(/iip1,jjp1,llm,1/)
1571	#ifdef NC_DOUBLE
1572	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1573	#else
1574	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1575	#endif
1576	CASE ("V")
1577	ierr = NF_INQ_VARID(nid,"vcov",varid)
1578	start=(/1,1,1,timestep/)
1579	count=(/iip1,jjm,llm,1/)
1580	#ifdef NC_DOUBLE
1581	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1582	#else
1583	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1584	#endif
1585	CASE ("T")
1586	ierr = NF_INQ_VARID(nid,"teta",varid)
1587	start=(/1,1,1,timestep/)
1588	count=(/iip1,jjp1,llm,1/)
1589	#ifdef NC_DOUBLE
1590	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1591	#else
1592	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1593	#endif
1594	CASE ("Q")
1595	ierr = NF_INQ_VARID(nid,"q",varid)
1596	start=(/1,1,1,timestep/)
1597	count=(/iip1,jjp1,llm,1/)
1598	#ifdef NC_DOUBLE
1599	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1600	#else
1601	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1602	#endif
1603	END SELECT
1604
1605	ierr = NF_CLOSE(nid)
1606
1607	END SUBROUTINE guide_out
1608
1609
1610	!===========================================================================
1611	subroutine correctbid(iim,nl,x)
1612	integer iim,nl
1613	real x(iim+1,nl)
1614	integer i,l
1615	real zz
1616
1617	do l=1,nl
1618	do i=2,iim-1
1619	if(abs(x(i,l)).gt.1.e10) then
1620	zz=0.5*(x(i-1,l)+x(i+1,l))
1621	print*,'correction ',i,l,x(i,l),zz
1622	x(i,l)=zz
1623	endif
1624	enddo
1625	enddo
1626	return
1627	end subroutine correctbid
1628
1629	!===========================================================================
1630	END MODULE guide_p_mod

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

Download in other formats:

Original Format