Context Navigation

guide_mod.F90 @ 2597

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

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

Download in other formats:

Original Format