Context Navigation

guide_p_mod.F90 @ 1632

Last change on this file since 1632 was 1632, checked in by Laurent Fairhead, 12 years ago

Import initial du répertoire dyn3dmem

Attention! ceci n'est qu'une version préliminaire du code "basse mémoire":
le code contenu dans ce répertoire est basé sur la r1320 et a donc besoin
d'être mis à jour par rapport à la dynamique parallèle d'aujourd'hui.
Ce code est toutefois mis à disposition pour circonvenir à des problèmes
de mémoire que certaines configurations du modèle pourraient rencontrer.
Dans l'état, il compile et tourne sur vargas et au CCRT

Initial import of dyn3dmem

Warning! this is just a preliminary version of the memory light code:
it is based on r1320 of the code and thus needs to be updated before
it can replace the present dyn3dpar code. It is nevertheless put at your
disposal to circumvent some memory problems some LMDZ configurations may
encounter. In its present state, it will compile and run on vargas and CCRT

File size: 56.7 KB

Line
1	!
2	! $Id$
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	USE control_mod
69	IMPLICIT NONE
70
71	INCLUDE "dimensions.h"
72	INCLUDE "paramet.h"
73	INCLUDE "netcdf.inc"
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	USE control_mod
277
278	IMPLICIT NONE
279
280	INCLUDE "dimensions.h"
281	INCLUDE "paramet.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/ REAL(iguide_read)
383	ELSE
384	tau= REAL(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)/ REAL(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)/ REAL(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	!ym rustine temporaire pour ne pas depasser 2GB pour la stack
635	!ym => crach compilo avec la version intel 12
636	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: p ! pression intercouches
637	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: pls, pext ! var intermediaire
638	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: pbarx
639	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: pbary
640	! Variables pour fonction Exner (P milieu couche)
641	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: pk, pkf
642	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: alpha, beta
643	REAL, DIMENSION (:,:),ALLOCATABLE,SAVE :: pks
644	REAL :: prefkap,unskap
645	! Pression de vapeur saturante
646	REAL, DIMENSION (ip1jmp1,llm) :: qsat
647	!Variables intermediaires interpolation
648
649	!ym rustine temporaire pour ne pas depasser 2GB pour la stack
650	!ym => crach compilo avec la version intel 12
651	! REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
652	! REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
653	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: zu1,zu2
654	REAL, DIMENSION (:,:,:),ALLOCATABLE,SAVE :: zv1,zv2
655
656	INTEGER :: i,j,l,ij
657	TYPE(Request) :: Req
658
659	print *,'Guide: conversion variables guidage'
660	! -----------------------------------------------------------------
661	! Calcul des niveaux de pression champs guidage
662	! -----------------------------------------------------------------
663	if (guide_modele) then
664	do i=1,iip1
665	do j=jjb_u,jje_u
666	do l=1,nlevnc
667	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
668	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
669	enddo
670	enddo
671	enddo
672	else
673	do i=1,iip1
674	do j=jjb_u,jje_u
675	do l=1,nlevnc
676	plnc2(i,j,l)=apnc(l)
677	plnc1(i,j,l)=apnc(l)
678	enddo
679	enddo
680	enddo
681
682	endif
683	if (first) then
684	first=.FALSE.
685	print*,'Guide: verification ordre niveaux verticaux'
686	print*,'LMDZ :'
687	do l=1,llm
688	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
689	+psi(1,jje_u)*(bp(l)+bp(l+1))/2.
690	enddo
691	print*,'Fichiers guidage'
692	do l=1,nlevnc
693	print*,'PL(',l,')=',plnc2(1,jjb_u,l)
694	enddo
695	print *,'inversion de l''ordre: invert_p=',invert_p
696	if (guide_u) then
697	do l=1,nlevnc
698	print*,'U(',l,')=',unat2(1,jjb_u,l)
699	enddo
700	endif
701	if (guide_T) then
702	do l=1,nlevnc
703	print*,'T(',l,')=',tnat2(1,jjb_u,l)
704	enddo
705	endif
706	ALLOCATE( zu1(iip1,jjp1,llm),zu2(iip1,jjp1,llm))
707	ALLOCATE( zv1(iip1,jjm,llm),zv2(iip1,jjm,llm))
708	ALLOCATE( p(iip1,jjp1,llmp1) )
709	ALLOCATE( pls(iip1,jjp1,llm), pext(iip1,jjp1,llm) )
710	ALLOCATE( pbarx(iip1,jjp1,llm) )
711	ALLOCATE( pbary(iip1,jjm,llm) )
712	ALLOCATE( pk(iip1,jjp1,llm),pkf(iip1,jjp1,llm) )
713	ALLOCATE( alpha(iip1,jjp1,llm),beta(iip1,jjp1,llm) )
714	ALLOCATE( pks(iip1,jjp1) )
715	endif
716
717	! -----------------------------------------------------------------
718	! Calcul niveaux pression modele
719	! -----------------------------------------------------------------
720	CALL pression_p( ip1jmp1, ap, bp, psi, p )
721	CALL exner_hyb_p(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
722
723	! .... Calcul de pls , pression au milieu des couches ,en Pascals
724	unskap=1./kappa
725	prefkap = preff ** kappa
726	DO l = 1, llm
727	DO j=jjb_u,jje_u
728	DO i =1, iip1
729	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
730	ENDDO
731	ENDDO
732	ENDDO
733
734	! calcul des pressions pour les grilles u et v
735	do l=1,llm
736	do j=jjb_u,jje_u
737	do i=1,iip1
738	pext(i,j,l)=pls(i,j,l)*aire(i,j)
739	enddo
740	enddo
741	enddo
742
743	CALL Register_SwapFieldHallo(pext,pext,ip1jmp1,llm,jj_Nb_caldyn,1,2,Req)
744	CALL SendRequest(Req)
745	CALL WaitRequest(Req)
746
747	call massbar_p(pext, pbarx, pbary )
748	do l=1,llm
749	do j=jjb_u,jje_u
750	do i=1,iip1
751	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
752	plsnc(i,j,l)=pls(i,j,l)
753	enddo
754	enddo
755	enddo
756	do l=1,llm
757	do j=jjb_v,jje_v
758	do i=1,iip1
759	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
760	enddo
761	enddo
762	enddo
763
764	! -----------------------------------------------------------------
765	! Interpolation champs guidage sur niveaux modele (+inversion N/S)
766	! Conversion en variables gcm (ucov, vcov...)
767	! -----------------------------------------------------------------
768	if (guide_P) then
769	do j=jjb_u,jje_u
770	do i=1,iim
771	ij=(j-1)*iip1+i
772	psgui1(ij)=psnat1(i,j)
773	psgui2(ij)=psnat2(i,j)
774	enddo
775	psgui1(iip1*j)=psnat1(1,j)
776	psgui2(iip1*j)=psnat2(1,j)
777	enddo
778	endif
779
780	IF (guide_u) THEN
781	CALL pres2lev(unat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
782	plnc1(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
783	CALL pres2lev(unat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
784	plnc2(:,jjb_u:jje_u,:),plunc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
785
786	do l=1,llm
787	do j=jjb_u,jje_u
788	do i=1,iim
789	ij=(j-1)*iip1+i
790	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
791	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
792	enddo
793	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
794	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
795	enddo
796	do i=1,iip1
797	ugui1(i,l)=0.
798	ugui1(ip1jm+i,l)=0.
799	ugui2(i,l)=0.
800	ugui2(ip1jm+i,l)=0.
801	enddo
802	enddo
803	ENDIF
804
805	IF (guide_T) THEN
806	CALL pres2lev(tnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
807	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
808	CALL pres2lev(tnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
809	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
810
811	do l=1,llm
812	do j=jjb_u,jje_u
813	IF (guide_teta) THEN
814	do i=1,iim
815	ij=(j-1)*iip1+i
816	tgui1(ij,l)=zu1(i,j,l)
817	tgui2(ij,l)=zu2(i,j,l)
818	enddo
819	ELSE
820	do i=1,iim
821	ij=(j-1)*iip1+i
822	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
823	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
824	enddo
825	ENDIF
826	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
827	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
828	enddo
829	do i=1,iip1
830	tgui1(i,l)=tgui1(1,l)
831	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
832	tgui2(i,l)=tgui2(1,l)
833	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
834	enddo
835	enddo
836	ENDIF
837
838	IF (guide_v) THEN
839
840	CALL pres2lev(vnat1(:,jjb_v:jje_v,:),zv1(:,jjb_v:jje_v,:),nlevnc,llm, &
841	plnc1(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
842	CALL pres2lev(vnat2(:,jjb_v:jje_v,:),zv2(:,jjb_v:jje_v,:),nlevnc,llm, &
843	plnc2(:,jjb_v:jje_v,:),plvnc(:,jjb_v:jje_v,:),iip1,jjn_v,invert_p)
844
845	do l=1,llm
846	do j=jjb_v,jje_v
847	do i=1,iim
848	ij=(j-1)*iip1+i
849	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
850	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
851	enddo
852	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
853	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
854	enddo
855	enddo
856	ENDIF
857
858	IF (guide_Q) THEN
859	! On suppose qu'on a la bonne variable dans le fichier de guidage:
860	! Hum.Rel si guide_hr, Hum.Spec. sinon.
861	CALL pres2lev(qnat1(:,jjb_u:jje_u,:),zu1(:,jjb_u:jje_u,:),nlevnc,llm, &
862	plnc1(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
863	CALL pres2lev(qnat2(:,jjb_u:jje_u,:),zu2(:,jjb_u:jje_u,:),nlevnc,llm, &
864	plnc2(:,jjb_u:jje_u,:),plsnc(:,jjb_u:jje_u,:),iip1,jjn_u,invert_p)
865
866	do l=1,llm
867	do j=jjb_u,jjb_v
868	do i=1,iim
869	ij=(j-1)*iip1+i
870	qgui1(ij,l)=zu1(i,j,l)
871	qgui2(ij,l)=zu2(i,j,l)
872	enddo
873	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
874	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
875	enddo
876	do i=1,iip1
877	qgui1(i,l)=qgui1(1,l)
878	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
879	qgui2(i,l)=qgui2(1,l)
880	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
881	enddo
882	enddo
883	IF (guide_hr) THEN
884	CALL q_sat(iip1jjn_ullm,teta(:,jjb_u:jje_u,:)*pk(:,jjb_u:jje_u,:)/cpp, &
885	plsnc(:,jjb_u:jje_u,:),qsat(ijb_u:ije_u,:))
886	qgui1(ijb_u:ije_u,:)=qgui1(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01 !hum. rel. en %
887	qgui2(ijb_u:ije_u,:)=qgui2(ijb_u:ije_u,:)qsat(ijb_u:ije_u,:)0.01
888	ENDIF
889	ENDIF
890
891	END SUBROUTINE guide_interp
892
893	!=======================================================================
894	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
895
896	! Calcul des constantes de rappel alpha (=1/tau)
897
898	implicit none
899
900	include "dimensions.h"
901	include "paramet.h"
902	include "comconst.h"
903	include "comgeom2.h"
904	include "serre.h"
905
906	! input arguments :
907	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
908	INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
909	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
910	REAL, INTENT(IN) :: taumin,taumax
911	! output arguments:
912	REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
913
914	! local variables:
915	LOGICAL, SAVE :: first=.TRUE.
916	REAL, SAVE :: gamma,dxdy_min,dxdy_max
917	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
918	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
919	REAL, DIMENSION (iip1,jjm) :: dxdyv
920	real dxdy_
921	real zlat,zlon
922	real alphamin,alphamax,xi
923	integer i,j,ilon,ilat
924
925
926	alphamin=factt/taumax
927	alphamax=factt/taumin
928	IF (guide_reg.OR.guide_add) THEN
929	alpha=alphamax
930	!-----------------------------------------------------------------------
931	! guide_reg: alpha=alpha_min dans region, 0. sinon.
932	!-----------------------------------------------------------------------
933	IF (guide_reg) THEN
934	do j=1,pjm
935	do i=1,pim
936	if (typ.eq.2) then
937	zlat=rlatu(j)*180./pi
938	zlon=rlonu(i)*180./pi
939	elseif (typ.eq.1) then
940	zlat=rlatu(j)*180./pi
941	zlon=rlonv(i)*180./pi
942	elseif (typ.eq.3) then
943	zlat=rlatv(j)*180./pi
944	zlon=rlonv(i)*180./pi
945	endif
946	alpha(i,j)=alphamax/16.* &
947	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
948	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
949	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
950	(1.+tanh((lon_max_g-zlon)/tau_lon))
951	enddo
952	enddo
953	ENDIF
954	ELSE
955	!-----------------------------------------------------------------------
956	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
957	!-----------------------------------------------------------------------
958	!Calcul de l'aire des mailles
959	do j=2,jjm
960	do i=2,iip1
961	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
962	enddo
963	zdx(1,j)=zdx(iip1,j)
964	enddo
965	do j=2,jjm
966	do i=1,iip1
967	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
968	enddo
969	enddo
970	do i=1,iip1
971	zdx(i,1)=zdx(i,2)
972	zdx(i,jjp1)=zdx(i,jjm)
973	zdy(i,1)=zdy(i,2)
974	zdy(i,jjp1)=zdy(i,jjm)
975	enddo
976	do j=1,jjp1
977	do i=1,iip1
978	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
979	enddo
980	enddo
981	IF (typ.EQ.2) THEN
982	do j=1,jjp1
983	do i=1,iim
984	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
985	enddo
986	dxdyu(iip1,j)=dxdyu(1,j)
987	enddo
988	ENDIF
989	IF (typ.EQ.3) THEN
990	do j=1,jjm
991	do i=1,iip1
992	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
993	enddo
994	enddo
995	ENDIF
996	! Premier appel: calcul des aires min et max et de gamma.
997	IF (first) THEN
998	first=.FALSE.
999	! coordonnees du centre du zoom
1000	CALL coordij(clon,clat,ilon,ilat)
1001	! aire de la maille au centre du zoom
1002	dxdy_min=dxdys(ilon,ilat)
1003	! dxdy maximale de la maille
1004	dxdy_max=0.
1005	do j=1,jjp1
1006	do i=1,iip1
1007	dxdy_max=max(dxdy_max,dxdys(i,j))
1008	enddo
1009	enddo
1010	! Calcul de gamma
1011	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1012	print*,'ATTENTION modele peu zoome'
1013	print*,'ATTENTION on prend une constante de guidage cste'
1014	gamma=0.
1015	else
1016	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
1017	print*,'gamma=',gamma
1018	if (gamma.lt.1.e-5) then
1019	print*,'gamma =',gamma,'<1e-5'
1020	stop
1021	endif
1022	gamma=log(0.5)/log(gamma)
1023	if (gamma4) then
1024	gamma=min(gamma,4.)
1025	endif
1026	print*,'gamma=',gamma
1027	endif
1028	ENDIF !first
1029
1030	do j=1,pjm
1031	do i=1,pim
1032	if (typ.eq.1) then
1033	dxdy_=dxdys(i,j)
1034	zlat=rlatu(j)*180./pi
1035	elseif (typ.eq.2) then
1036	dxdy_=dxdyu(i,j)
1037	zlat=rlatu(j)*180./pi
1038	elseif (typ.eq.3) then
1039	dxdy_=dxdyv(i,j)
1040	zlat=rlatv(j)*180./pi
1041	endif
1042	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
1043	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
1044	alpha(i,j)=alphamin
1045	else
1046	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
1047	xi=min(xi,1.)
1048	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
1049	alpha(i,j)=xialphamin+(1.-xi)alphamax
1050	else
1051	alpha(i,j)=0.
1052	endif
1053	endif
1054	enddo
1055	enddo
1056	ENDIF ! guide_reg
1057
1058	END SUBROUTINE tau2alpha
1059
1060	!=======================================================================
1061	SUBROUTINE guide_read(timestep)
1062
1063	IMPLICIT NONE
1064
1065	#include "netcdf.inc"
1066	#include "dimensions.h"
1067	#include "paramet.h"
1068
1069	INTEGER, INTENT(IN) :: timestep
1070
1071	LOGICAL, SAVE :: first=.TRUE.
1072	! Identification fichiers et variables NetCDF:
1073	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1074	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1075	INTEGER :: ncidpl,varidpl,varidap,varidbp
1076	! Variables auxiliaires NetCDF:
1077	INTEGER, DIMENSION(4) :: start,count
1078	INTEGER :: status,rcode
1079
1080	! -----------------------------------------------------------------
1081	! Premier appel: initialisation de la lecture des fichiers
1082	! -----------------------------------------------------------------
1083	if (first) then
1084	ncidpl=-99
1085	print*,'Guide: ouverture des fichiers guidage '
1086	! Niveaux de pression si non constants
1087	if (guide_modele) then
1088	print *,'Lecture du guidage sur niveaux mod�le'
1089	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1090	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1091	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1092	print*,'ncidpl,varidap',ncidpl,varidap
1093	endif
1094	! Vent zonal
1095	if (guide_u) then
1096	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1097	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1098	print*,'ncidu,varidu',ncidu,varidu
1099	if (ncidpl.eq.-99) ncidpl=ncidu
1100	endif
1101	! Vent meridien
1102	if (guide_v) then
1103	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1104	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1105	print*,'ncidv,varidv',ncidv,varidv
1106	if (ncidpl.eq.-99) ncidpl=ncidv
1107	endif
1108	! Temperature
1109	if (guide_T) then
1110	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1111	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1112	print*,'ncidT,varidT',ncidt,varidt
1113	if (ncidpl.eq.-99) ncidpl=ncidt
1114	endif
1115	! Humidite
1116	if (guide_Q) then
1117	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1118	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1119	print*,'ncidQ,varidQ',ncidQ,varidQ
1120	if (ncidpl.eq.-99) ncidpl=ncidQ
1121	endif
1122	! Pression de surface
1123	if ((guide_P).OR.(guide_modele)) then
1124	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1125	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1126	print*,'ncidps,varidps',ncidps,varidps
1127	endif
1128	! Coordonnee verticale
1129	if (.not.guide_modele) then
1130	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1131	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1132	print*,'ncidpl,varidpl',ncidpl,varidpl
1133	endif
1134	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1135	if (guide_modele) then
1136	#ifdef NC_DOUBLE
1137	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1138	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1139	#else
1140	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1141	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1142	#endif
1143	else
1144	#ifdef NC_DOUBLE
1145	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1146	#else
1147	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1148	#endif
1149	apnc=apnc*100.! conversion en Pascals
1150	bpnc(:)=0.
1151	endif
1152	first=.FALSE.
1153	endif ! (first)
1154
1155	! -----------------------------------------------------------------
1156	! lecture des champs u, v, T, Q, ps
1157	! -----------------------------------------------------------------
1158
1159	! dimensions pour les champs scalaires et le vent zonal
1160	start(1)=1
1161	start(2)=1
1162	start(3)=1
1163	start(4)=timestep
1164
1165	count(1)=iip1
1166	count(2)=jjp1
1167	count(3)=nlevnc
1168	count(4)=1
1169
1170	! Vent zonal
1171	if (guide_u) then
1172	#ifdef NC_DOUBLE
1173	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
1174	#else
1175	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
1176	#endif
1177	IF (invert_y) THEN
1178	CALL invert_lat(iip1,jjp1,llm,unat2)
1179	ENDIF
1180
1181	endif
1182
1183	! Temperature
1184	if (guide_T) then
1185	#ifdef NC_DOUBLE
1186	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
1187	#else
1188	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
1189	#endif
1190	IF (invert_y) THEN
1191	CALL invert_lat(iip1,jjp1,llm,tnat2)
1192	ENDIF
1193	endif
1194
1195	! Humidite
1196	if (guide_Q) then
1197	#ifdef NC_DOUBLE
1198	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
1199	#else
1200	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
1201	#endif
1202	IF (invert_y) THEN
1203	CALL invert_lat(iip1,jjp1,llm,qnat2)
1204	ENDIF
1205
1206	endif
1207
1208	! Vent meridien
1209	if (guide_v) then
1210	count(2)=jjm
1211	#ifdef NC_DOUBLE
1212	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
1213	#else
1214	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
1215	#endif
1216	IF (invert_y) THEN
1217	CALL invert_lat(iip1,jjm,llm,vnat2)
1218	ENDIF
1219	endif
1220
1221	! Pression de surface
1222	if ((guide_P).OR.(guide_modele)) then
1223	start(3)=timestep
1224	start(4)=0
1225	count(2)=jjp1
1226	count(3)=1
1227	count(4)=0
1228	#ifdef NC_DOUBLE
1229	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
1230	#else
1231	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
1232	#endif
1233	IF (invert_y) THEN
1234	CALL invert_lat(iip1,jjp1,1,psnat2)
1235	ENDIF
1236	endif
1237
1238	END SUBROUTINE guide_read
1239
1240	!=======================================================================
1241	SUBROUTINE guide_read2D(timestep)
1242
1243	IMPLICIT NONE
1244
1245	#include "netcdf.inc"
1246	#include "dimensions.h"
1247	#include "paramet.h"
1248
1249	INTEGER, INTENT(IN) :: timestep
1250
1251	LOGICAL, SAVE :: first=.TRUE.
1252	! Identification fichiers et variables NetCDF:
1253	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1254	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1255	INTEGER :: ncidpl,varidpl,varidap,varidbp
1256	! Variables auxiliaires NetCDF:
1257	INTEGER, DIMENSION(4) :: start,count
1258	INTEGER :: status,rcode
1259	! Variables for 3D extension:
1260	REAL, DIMENSION (jjp1,llm) :: zu
1261	REAL, DIMENSION (jjm,llm) :: zv
1262	INTEGER :: i
1263
1264	! -----------------------------------------------------------------
1265	! Premier appel: initialisation de la lecture des fichiers
1266	! -----------------------------------------------------------------
1267	if (first) then
1268	ncidpl=-99
1269	print*,'Guide: ouverture des fichiers guidage '
1270	! Niveaux de pression si non constants
1271	if (guide_modele) then
1272	print *,'Lecture du guidage sur niveaux mod�le'
1273	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1274	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1275	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1276	print*,'ncidpl,varidap',ncidpl,varidap
1277	endif
1278	! Vent zonal
1279	if (guide_u) then
1280	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1281	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1282	print*,'ncidu,varidu',ncidu,varidu
1283	if (ncidpl.eq.-99) ncidpl=ncidu
1284	endif
1285	! Vent meridien
1286	if (guide_v) then
1287	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1288	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1289	print*,'ncidv,varidv',ncidv,varidv
1290	if (ncidpl.eq.-99) ncidpl=ncidv
1291	endif
1292	! Temperature
1293	if (guide_T) then
1294	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1295	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1296	print*,'ncidT,varidT',ncidt,varidt
1297	if (ncidpl.eq.-99) ncidpl=ncidt
1298	endif
1299	! Humidite
1300	if (guide_Q) then
1301	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1302	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1303	print*,'ncidQ,varidQ',ncidQ,varidQ
1304	if (ncidpl.eq.-99) ncidpl=ncidQ
1305	endif
1306	! Pression de surface
1307	if ((guide_P).OR.(guide_modele)) then
1308	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1309	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1310	print*,'ncidps,varidps',ncidps,varidps
1311	endif
1312	! Coordonnee verticale
1313	if (.not.guide_modele) then
1314	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1315	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1316	print*,'ncidpl,varidpl',ncidpl,varidpl
1317	endif
1318	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1319	if (guide_modele) then
1320	#ifdef NC_DOUBLE
1321	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1322	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1323	#else
1324	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1325	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1326	#endif
1327	else
1328	#ifdef NC_DOUBLE
1329	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1330	#else
1331	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1332	#endif
1333	apnc=apnc*100.! conversion en Pascals
1334	bpnc(:)=0.
1335	endif
1336	first=.FALSE.
1337	endif ! (first)
1338
1339	! -----------------------------------------------------------------
1340	! lecture des champs u, v, T, Q, ps
1341	! -----------------------------------------------------------------
1342
1343	! dimensions pour les champs scalaires et le vent zonal
1344	start(1)=1
1345	start(2)=1
1346	start(3)=1
1347	start(4)=timestep
1348
1349	count(1)=1
1350	count(2)=jjp1
1351	count(3)=nlevnc
1352	count(4)=1
1353
1354	! Vent zonal
1355	if (guide_u) then
1356	#ifdef NC_DOUBLE
1357	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
1358	#else
1359	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
1360	#endif
1361	DO i=1,iip1
1362	unat2(i,:,:)=zu(:,:)
1363	ENDDO
1364
1365	IF (invert_y) THEN
1366	CALL invert_lat(iip1,jjp1,llm,unat2)
1367	ENDIF
1368
1369	endif
1370
1371	! Temperature
1372	if (guide_T) then
1373	#ifdef NC_DOUBLE
1374	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
1375	#else
1376	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
1377	#endif
1378	DO i=1,iip1
1379	tnat2(i,:,:)=zu(:,:)
1380	ENDDO
1381
1382	IF (invert_y) THEN
1383	CALL invert_lat(iip1,jjp1,llm,tnat2)
1384	ENDIF
1385
1386	endif
1387
1388	! Humidite
1389	if (guide_Q) then
1390	#ifdef NC_DOUBLE
1391	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
1392	#else
1393	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
1394	#endif
1395	DO i=1,iip1
1396	qnat2(i,:,:)=zu(:,:)
1397	ENDDO
1398
1399	IF (invert_y) THEN
1400	CALL invert_lat(iip1,jjp1,llm,qnat2)
1401	ENDIF
1402
1403	endif
1404
1405	! Vent meridien
1406	if (guide_v) then
1407	count(2)=jjm
1408	#ifdef NC_DOUBLE
1409	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
1410	#else
1411	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
1412	#endif
1413	DO i=1,iip1
1414	vnat2(i,:,:)=zv(:,:)
1415	ENDDO
1416
1417	IF (invert_y) THEN
1418	CALL invert_lat(iip1,jjm,llm,vnat2)
1419	ENDIF
1420
1421	endif
1422
1423	! Pression de surface
1424	if ((guide_P).OR.(guide_modele)) then
1425	start(3)=timestep
1426	start(4)=0
1427	count(2)=jjp1
1428	count(3)=1
1429	count(4)=0
1430	#ifdef NC_DOUBLE
1431	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
1432	#else
1433	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
1434	#endif
1435	DO i=1,iip1
1436	psnat2(i,:)=zu(:,1)
1437	ENDDO
1438
1439	IF (invert_y) THEN
1440	CALL invert_lat(iip1,jjp1,1,psnat2)
1441	ENDIF
1442
1443	endif
1444
1445	END SUBROUTINE guide_read2D
1446
1447	!=======================================================================
1448	SUBROUTINE guide_out(varname,hsize,vsize,field)
1449	USE parallel
1450	IMPLICIT NONE
1451
1452	INCLUDE "dimensions.h"
1453	INCLUDE "paramet.h"
1454	INCLUDE "netcdf.inc"
1455	INCLUDE "comgeom2.h"
1456	INCLUDE "comconst.h"
1457	INCLUDE "comvert.h"
1458
1459	! Variables entree
1460	CHARACTER, INTENT(IN) :: varname
1461	INTEGER, INTENT (IN) :: hsize,vsize
1462	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1463
1464	! Variables locales
1465	INTEGER, SAVE :: timestep=0
1466	! Identites fichier netcdf
1467	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1468	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1469	INTEGER, DIMENSION (3) :: dim3
1470	INTEGER, DIMENSION (4) :: dim4,count,start
1471	INTEGER :: ierr, varid
1472
1473	CALL gather_field(field,iip1*hsize,vsize,0)
1474
1475	IF (mpi_rank /= 0) RETURN
1476
1477	print *,'Guide: output timestep',timestep,'var ',varname
1478	IF (timestep.EQ.0) THEN
1479	! ----------------------------------------------
1480	! initialisation fichier de sortie
1481	! ----------------------------------------------
1482	! Ouverture du fichier
1483	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
1484	! Definition des dimensions
1485	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1486	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1487	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1488	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1489	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1490	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1491
1492	! Creation des variables dimensions
1493	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1494	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1495	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1496	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1497	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1498	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1499	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1500
1501	ierr=NF_ENDDEF(nid)
1502
1503	! Enregistrement des variables dimensions
1504	#ifdef NC_DOUBLE
1505	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1506	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1507	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1508	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1509	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1510	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1511	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1512	#else
1513	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1514	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1515	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1516	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1517	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1518	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1519	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1520	#endif
1521	! --------------------------------------------------------------------
1522	! Cr�ation des variables sauvegard�es
1523	! --------------------------------------------------------------------
1524	ierr = NF_REDEF(nid)
1525	! Surface pressure (GCM)
1526	dim3=(/id_lonv,id_latu,id_tim/)
1527	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
1528	! Surface pressure (guidage)
1529	IF (guide_P) THEN
1530	dim3=(/id_lonv,id_latu,id_tim/)
1531	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1532	ENDIF
1533	! Zonal wind
1534	IF (guide_u) THEN
1535	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1536	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1537	ENDIF
1538	! Merid. wind
1539	IF (guide_v) THEN
1540	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1541	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1542	ENDIF
1543	! Pot. Temperature
1544	IF (guide_T) THEN
1545	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1546	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1547	ENDIF
1548	! Specific Humidity
1549	IF (guide_Q) THEN
1550	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1551	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1552	ENDIF
1553
1554	ierr = NF_ENDDEF(nid)
1555	ierr = NF_CLOSE(nid)
1556	ENDIF ! timestep=0
1557
1558	! --------------------------------------------------------------------
1559	! Enregistrement du champ
1560	! --------------------------------------------------------------------
1561	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1562
1563	SELECT CASE (varname)
1564	CASE ("S")
1565	timestep=timestep+1
1566	ierr = NF_INQ_VARID(nid,"SP",varid)
1567	start=(/1,1,timestep,0/)
1568	count=(/iip1,jjp1,1,0/)
1569	#ifdef NC_DOUBLE
1570	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1571	#else
1572	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1573	#endif
1574	CASE ("P")
1575	ierr = NF_INQ_VARID(nid,"ps",varid)
1576	start=(/1,1,timestep,0/)
1577	count=(/iip1,jjp1,1,0/)
1578	#ifdef NC_DOUBLE
1579	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1580	#else
1581	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1582	#endif
1583	CASE ("U")
1584	ierr = NF_INQ_VARID(nid,"ucov",varid)
1585	start=(/1,1,1,timestep/)
1586	count=(/iip1,jjp1,llm,1/)
1587	#ifdef NC_DOUBLE
1588	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1589	#else
1590	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1591	#endif
1592	CASE ("V")
1593	ierr = NF_INQ_VARID(nid,"vcov",varid)
1594	start=(/1,1,1,timestep/)
1595	count=(/iip1,jjm,llm,1/)
1596	#ifdef NC_DOUBLE
1597	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1598	#else
1599	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1600	#endif
1601	CASE ("T")
1602	ierr = NF_INQ_VARID(nid,"teta",varid)
1603	start=(/1,1,1,timestep/)
1604	count=(/iip1,jjp1,llm,1/)
1605	#ifdef NC_DOUBLE
1606	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1607	#else
1608	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1609	#endif
1610	CASE ("Q")
1611	ierr = NF_INQ_VARID(nid,"q",varid)
1612	start=(/1,1,1,timestep/)
1613	count=(/iip1,jjp1,llm,1/)
1614	#ifdef NC_DOUBLE
1615	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field)
1616	#else
1617	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field)
1618	#endif
1619	END SELECT
1620
1621	ierr = NF_CLOSE(nid)
1622
1623	END SUBROUTINE guide_out
1624
1625
1626	!===========================================================================
1627	subroutine correctbid(iim,nl,x)
1628	integer iim,nl
1629	real x(iim+1,nl)
1630	integer i,l
1631	real zz
1632
1633	do l=1,nl
1634	do i=2,iim-1
1635	if(abs(x(i,l)).gt.1.e10) then
1636	zz=0.5*(x(i-1,l)+x(i+1,l))
1637	print*,'correction ',i,l,x(i,l),zz
1638	x(i,l)=zz
1639	endif
1640	enddo
1641	enddo
1642	return
1643	end subroutine correctbid
1644
1645	!===========================================================================
1646	END MODULE guide_p_mod

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

Download in other formats:

Original Format