Context Navigation

guide_mod.F90 @ 1457

Last change on this file since 1457 was 1441, checked in by emillour, 10 years ago

Updates in common dynamics (seq and ) to keep up with updates
in LMDZ5 (up to LMDZ5 trunk, rev 2250):

compilation:
added test in grid/dimension/makdim to check that # of longitudes is a multiple of 8

dyn3d_common:

Bug correction concerning zoom (cf LMDZ5 rev 2218)

coefpoly.F becomes coefpoly_m.F90 (in misc)
fxhyp.F => fxhyp_m.F90 , fyhyp.F => fyhyp_m.F90
new routines for zoom: invert_zoom_x_m.F90 and principal_cshift_m.F90
inigeom.F adapted to new zoom definition routines
fluxstokenc.F : got rid of calls to initial0()

dyn3d:
advtrac.F90 : got rid of calls to initial0()
conf_gcm.F90 : cosmetic changes and change in default dzoomx,dzoomy values
guide_mod.F90 : followed updates from Earth Model
gcm.F is now gcm.F90

dyn3dpar:
advtrac_p.F90, covcont_p.F90, mod_hallo.F90 : cosmetic changes
conf_gcm.F90 : cosmetic and changed in default dzoomx,dzoomy values
parallel_lmdz.F90 : updates to keep up with Earth model

misc:
arth.F90 becomes arth_m.F90
wxios.F90 updated wrt Earth model changes
nrtype.F90 and coefpoly_m.F90 added
ran1.F, sort.F, minmax.F, minmax2.F, juldate.F moved over from dyn3d_common

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

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

Download in other formats:

Original Format