Context Navigation

guide_mod.F90 @ 2128

Last change on this file since 2128 was 2124, checked in by lguez, 10 years ago

In guide_main, there is the line:

f_out=((MOD(itau,iguide_sav).EQ.0).AND.guide_sav)

The Fortran standard says that the result of mod(a, p) with p = 0 is
processor dependent. With ifort (on Ada at IDRIS), it produces a
run-time error. This is now fixed. The user must still choose
iguide_sav=0 in order to write only once.

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

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

Download in other formats:

Original Format