Context Navigation

guide_mod.F90 @ 1422

Last change on this file since 1422 was 1422, checked in by milmd, 10 years ago
In GENERIC, MARS and COMMON models replace some include files by modules (usefull for decoupling physics with dynamics).
File size: 58.5 KB

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

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

Download in other formats:

Original Format