Context Navigation

guide_mod.F90 @ 5087

Last change on this file since 5087 was 2025, checked in by fhourdin, 11 years ago
Enrichissement des sorties pour le guidage Enriched outputs for nudging
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:executable set to ``*
File size: 58.0 KB

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

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

Download in other formats:

Original Format