Context Navigation

guide_mod.F90 @ 2056

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

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

Download in other formats:

Original Format