Context Navigation

guide_mod.F90 @ 2260

Last change on this file since 2260 was 2249, checked in by lguez, 9 years ago
abort_message was not defined.
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.3 KB

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

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

Download in other formats:

Original Format