Context Navigation

guide_p_mod.F90 @ 1669

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

Version testing basée sur la r1668

http://lmdz.lmd.jussieu.fr/utilisateurs/distribution-du-modele/versions-intermediaires

Testing release based on r1668

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,nlevnc,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,nlevnc,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,nlevnc,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,nlevnc,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,nlevnc,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,nlevnc,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,nlevnc,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,nlevnc,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