Context Navigation

guide_mod.F90 @ 2024

Last change on this file since 2024 was 2024, checked in by fhourdin, 10 years ago
Enrichissement des sorties du 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.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
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	IMPLICIT NONE
596
597	include "dimensions.h"
598	include "paramet.h"
599	include "comvert.h"
600	include "comgeom2.h"
601	include "comconst.h"
602
603	REAL, DIMENSION (iip1,jjp1), INTENT(IN) :: psi ! Psol gcm
604	REAL, DIMENSION (iip1,jjp1,llm), INTENT(IN) :: teta ! Temp. Pot. gcm
605
606	LOGICAL, SAVE :: first=.TRUE.
607	! Variables pour niveaux pression:
608	REAL, DIMENSION (iip1,jjp1,nlevnc) :: plnc1,plnc2 !niveaux pression guidage
609	REAL, DIMENSION (iip1,jjp1,llm) :: plunc,plsnc !niveaux pression modele
610	REAL, DIMENSION (iip1,jjm,llm) :: plvnc !niveaux pression modele
611	REAL, DIMENSION (iip1,jjp1,llmp1) :: p ! pression intercouches
612	REAL, DIMENSION (iip1,jjp1,llm) :: pls, pext ! var intermediaire
613	REAL, DIMENSION (iip1,jjp1,llm) :: pbarx
614	REAL, DIMENSION (iip1,jjm,llm) :: pbary
615	! Variables pour fonction Exner (P milieu couche)
616	REAL, DIMENSION (iip1,jjp1,llm) :: pk, pkf
617	REAL, DIMENSION (iip1,jjp1,llm) :: alpha, beta
618	REAL, DIMENSION (iip1,jjp1) :: pks
619	REAL :: prefkap,unskap
620	! Pression de vapeur saturante
621	REAL, DIMENSION (ip1jmp1,llm) :: qsat
622	!Variables intermediaires interpolation
623	REAL, DIMENSION (iip1,jjp1,llm) :: zu1,zu2
624	REAL, DIMENSION (iip1,jjm,llm) :: zv1,zv2
625
626	INTEGER :: i,j,l,ij
627
628	print *,'Guide: conversion variables guidage'
629	! -----------------------------------------------------------------
630	! Calcul des niveaux de pression champs guidage
631	! -----------------------------------------------------------------
632	if (guide_modele) then
633	do i=1,iip1
634	do j=1,jjp1
635	do l=1,nlevnc
636	plnc2(i,j,l)=apnc(l)+bpnc(l)*psnat2(i,j)
637	plnc1(i,j,l)=apnc(l)+bpnc(l)*psnat1(i,j)
638	enddo
639	enddo
640	enddo
641	else
642	do i=1,iip1
643	do j=1,jjp1
644	do l=1,nlevnc
645	plnc2(i,j,l)=apnc(l)
646	plnc1(i,j,l)=apnc(l)
647	enddo
648	enddo
649	enddo
650
651	endif
652	if (first) then
653	first=.FALSE.
654	print*,'Guide: verification ordre niveaux verticaux'
655	print*,'LMDZ :'
656	do l=1,llm
657	print*,'PL(',l,')=',(ap(l)+ap(l+1))/2. &
658	+psi(1,jjp1)*(bp(l)+bp(l+1))/2.
659	enddo
660	print*,'Fichiers guidage'
661	do l=1,nlevnc
662	print*,'PL(',l,')=',plnc2(1,1,l)
663	enddo
664	print *,'inversion de l''ordre: invert_p=',invert_p
665	if (guide_u) then
666	do l=1,nlevnc
667	print*,'U(',l,')=',unat2(1,1,l)
668	enddo
669	endif
670	if (guide_T) then
671	do l=1,nlevnc
672	print*,'T(',l,')=',tnat2(1,1,l)
673	enddo
674	endif
675	endif
676
677	! -----------------------------------------------------------------
678	! Calcul niveaux pression modele
679	! -----------------------------------------------------------------
680	CALL pression( ip1jmp1, ap, bp, psi, p )
681	if (pressure_exner) then
682	CALL exner_hyb(ip1jmp1,psi,p,alpha,beta,pks,pk,pkf)
683	else
684	CALL exner_milieu(ip1jmp1,psi,p,beta,pks,pk,pkf)
685	endif
686	! .... Calcul de pls , pression au milieu des couches ,en Pascals
687	unskap=1./kappa
688	prefkap = preff ** kappa
689	DO l = 1, llm
690	DO j=1,jjp1
691	DO i =1, iip1
692	pls(i,j,l) = preff * ( pk(i,j,l)/cpp) ** unskap
693	ENDDO
694	ENDDO
695	ENDDO
696
697	! calcul des pressions pour les grilles u et v
698	do l=1,llm
699	do j=1,jjp1
700	do i=1,iip1
701	pext(i,j,l)=pls(i,j,l)*aire(i,j)
702	enddo
703	enddo
704	enddo
705	call massbar(pext, pbarx, pbary )
706	do l=1,llm
707	do j=1,jjp1
708	do i=1,iip1
709	plunc(i,j,l)=pbarx(i,j,l)/aireu(i,j)
710	plsnc(i,j,l)=pls(i,j,l)
711	enddo
712	enddo
713	enddo
714	do l=1,llm
715	do j=1,jjm
716	do i=1,iip1
717	plvnc(i,j,l)=pbary(i,j,l)/airev(i,j)
718	enddo
719	enddo
720	enddo
721
722	! -----------------------------------------------------------------
723	! Interpolation champs guidage sur niveaux modele (+inversion N/S)
724	! Conversion en variables gcm (ucov, vcov...)
725	! -----------------------------------------------------------------
726	if (guide_P) then
727	do j=1,jjp1
728	do i=1,iim
729	ij=(j-1)*iip1+i
730	psgui1(ij)=psnat1(i,j)
731	psgui2(ij)=psnat2(i,j)
732	enddo
733	psgui1(iip1*j)=psnat1(1,j)
734	psgui2(iip1*j)=psnat2(1,j)
735	enddo
736	endif
737
738	IF (guide_u) THEN
739	CALL pres2lev(unat1,zu1,nlevnc,llm,plnc1,plunc,iip1,jjp1,invert_p)
740	CALL pres2lev(unat2,zu2,nlevnc,llm,plnc2,plunc,iip1,jjp1,invert_p)
741	do l=1,llm
742	do j=1,jjp1
743	do i=1,iim
744	ij=(j-1)*iip1+i
745	ugui1(ij,l)=zu1(i,j,l)*cu(i,j)
746	ugui2(ij,l)=zu2(i,j,l)*cu(i,j)
747	enddo
748	ugui1(jiip1,l)=ugui1((j-1)iip1+1,l)
749	ugui2(jiip1,l)=ugui2((j-1)iip1+1,l)
750	enddo
751	do i=1,iip1
752	ugui1(i,l)=0.
753	ugui1(ip1jm+i,l)=0.
754	ugui2(i,l)=0.
755	ugui2(ip1jm+i,l)=0.
756	enddo
757	enddo
758	ENDIF
759
760	IF (guide_T) THEN
761	CALL pres2lev(tnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
762	CALL pres2lev(tnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
763	do l=1,llm
764	do j=1,jjp1
765	IF (guide_teta) THEN
766	do i=1,iim
767	ij=(j-1)*iip1+i
768	tgui1(ij,l)=zu1(i,j,l)
769	tgui2(ij,l)=zu2(i,j,l)
770	enddo
771	ELSE
772	do i=1,iim
773	ij=(j-1)*iip1+i
774	tgui1(ij,l)=zu1(i,j,l)*cpp/pk(i,j,l)
775	tgui2(ij,l)=zu2(i,j,l)*cpp/pk(i,j,l)
776	enddo
777	ENDIF
778	tgui1(jiip1,l)=tgui1((j-1)iip1+1,l)
779	tgui2(jiip1,l)=tgui2((j-1)iip1+1,l)
780	enddo
781	do i=1,iip1
782	tgui1(i,l)=tgui1(1,l)
783	tgui1(ip1jm+i,l)=tgui1(ip1jm+1,l)
784	tgui2(i,l)=tgui2(1,l)
785	tgui2(ip1jm+i,l)=tgui2(ip1jm+1,l)
786	enddo
787	enddo
788	ENDIF
789
790	IF (guide_v) THEN
791
792	CALL pres2lev(vnat1,zv1,nlevnc,llm,plnc1(:,:jjm,:),plvnc,iip1,jjm,invert_p)
793	CALL pres2lev(vnat2,zv2,nlevnc,llm,plnc2(:,:jjm,:),plvnc,iip1,jjm,invert_p)
794
795	do l=1,llm
796	do j=1,jjm
797	do i=1,iim
798	ij=(j-1)*iip1+i
799	vgui1(ij,l)=zv1(i,j,l)*cv(i,j)
800	vgui2(ij,l)=zv2(i,j,l)*cv(i,j)
801	enddo
802	vgui1(jiip1,l)=vgui1((j-1)iip1+1,l)
803	vgui2(jiip1,l)=vgui2((j-1)iip1+1,l)
804	enddo
805	enddo
806	ENDIF
807
808	IF (guide_Q) THEN
809	! On suppose qu'on a la bonne variable dans le fichier de guidage:
810	! Hum.Rel si guide_hr, Hum.Spec. sinon.
811	CALL pres2lev(qnat1,zu1,nlevnc,llm,plnc1,plsnc,iip1,jjp1,invert_p)
812	CALL pres2lev(qnat2,zu2,nlevnc,llm,plnc2,plsnc,iip1,jjp1,invert_p)
813	do l=1,llm
814	do j=1,jjp1
815	do i=1,iim
816	ij=(j-1)*iip1+i
817	qgui1(ij,l)=zu1(i,j,l)
818	qgui2(ij,l)=zu2(i,j,l)
819	enddo
820	qgui1(jiip1,l)=qgui1((j-1)iip1+1,l)
821	qgui2(jiip1,l)=qgui2((j-1)iip1+1,l)
822	enddo
823	do i=1,iip1
824	qgui1(i,l)=qgui1(1,l)
825	qgui1(ip1jm+i,l)=qgui1(ip1jm+1,l)
826	qgui2(i,l)=qgui2(1,l)
827	qgui2(ip1jm+i,l)=qgui2(ip1jm+1,l)
828	enddo
829	enddo
830	IF (guide_hr) THEN
831	CALL q_sat(iip1jjp1llm,teta*pk/cpp,plsnc,qsat)
832	qgui1=qgui1qsat0.01 !hum. rel. en %
833	qgui2=qgui2qsat0.01
834	ENDIF
835	ENDIF
836
837	END SUBROUTINE guide_interp
838
839	!=======================================================================
840	SUBROUTINE tau2alpha(typ,pim,pjm,factt,taumin,taumax,alpha)
841
842	! Calcul des constantes de rappel alpha (=1/tau)
843
844	implicit none
845
846	include "dimensions.h"
847	include "paramet.h"
848	include "comconst.h"
849	include "comgeom2.h"
850	include "serre.h"
851
852	! input arguments :
853	INTEGER, INTENT(IN) :: typ ! u(2),v(3), ou scalaire(1)
854	INTEGER, INTENT(IN) :: pim,pjm ! dimensions en lat, lon
855	REAL, INTENT(IN) :: factt ! pas de temps en fraction de jour
856	REAL, INTENT(IN) :: taumin,taumax
857	! output arguments:
858	REAL, DIMENSION(pim,pjm), INTENT(OUT) :: alpha
859
860	! local variables:
861	LOGICAL, SAVE :: first=.TRUE.
862	REAL, SAVE :: gamma,dxdy_min,dxdy_max
863	REAL, DIMENSION (iip1,jjp1) :: zdx,zdy
864	REAL, DIMENSION (iip1,jjp1) :: dxdys,dxdyu
865	REAL, DIMENSION (iip1,jjm) :: dxdyv
866	real dxdy_
867	real zlat,zlon
868	real alphamin,alphamax,xi
869	integer i,j,ilon,ilat
870
871
872	alphamin=factt/taumax
873	alphamax=factt/taumin
874	IF (guide_reg.OR.guide_add) THEN
875	alpha=alphamax
876	!-----------------------------------------------------------------------
877	! guide_reg: alpha=alpha_min dans region, 0. sinon.
878	!-----------------------------------------------------------------------
879	IF (guide_reg) THEN
880	do j=1,pjm
881	do i=1,pim
882	if (typ.eq.2) then
883	zlat=rlatu(j)*180./pi
884	zlon=rlonu(i)*180./pi
885	elseif (typ.eq.1) then
886	zlat=rlatu(j)*180./pi
887	zlon=rlonv(i)*180./pi
888	elseif (typ.eq.3) then
889	zlat=rlatv(j)*180./pi
890	zlon=rlonv(i)*180./pi
891	endif
892	alpha(i,j)=alphamax/16.* &
893	(1.+tanh((zlat-lat_min_g)/tau_lat))* &
894	(1.+tanh((lat_max_g-zlat)/tau_lat))* &
895	(1.+tanh((zlon-lon_min_g)/tau_lon))* &
896	(1.+tanh((lon_max_g-zlon)/tau_lon))
897	enddo
898	enddo
899	ENDIF
900	ELSE
901	!-----------------------------------------------------------------------
902	! Sinon, alpha varie entre alpha_min et alpha_max suivant le zoom.
903	!-----------------------------------------------------------------------
904	!Calcul de l'aire des mailles
905	do j=2,jjm
906	do i=2,iip1
907	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
908	enddo
909	zdx(1,j)=zdx(iip1,j)
910	enddo
911	do j=2,jjm
912	do i=1,iip1
913	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
914	enddo
915	enddo
916	do i=1,iip1
917	zdx(i,1)=zdx(i,2)
918	zdx(i,jjp1)=zdx(i,jjm)
919	zdy(i,1)=zdy(i,2)
920	zdy(i,jjp1)=zdy(i,jjm)
921	enddo
922	do j=1,jjp1
923	do i=1,iip1
924	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
925	enddo
926	enddo
927	IF (typ.EQ.2) THEN
928	do j=1,jjp1
929	do i=1,iim
930	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
931	enddo
932	dxdyu(iip1,j)=dxdyu(1,j)
933	enddo
934	ENDIF
935	IF (typ.EQ.3) THEN
936	do j=1,jjm
937	do i=1,iip1
938	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i,j+1))
939	enddo
940	enddo
941	ENDIF
942	! Premier appel: calcul des aires min et max et de gamma.
943	IF (first) THEN
944	first=.FALSE.
945	! coordonnees du centre du zoom
946	CALL coordij(clon,clat,ilon,ilat)
947	! aire de la maille au centre du zoom
948	dxdy_min=dxdys(ilon,ilat)
949	! dxdy maximale de la maille
950	dxdy_max=0.
951	do j=1,jjp1
952	do i=1,iip1
953	dxdy_max=max(dxdy_max,dxdys(i,j))
954	enddo
955	enddo
956	! Calcul de gamma
957	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
958	print*,'ATTENTION modele peu zoome'
959	print*,'ATTENTION on prend une constante de guidage cste'
960	gamma=0.
961	else
962	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
963	print*,'gamma=',gamma
964	if (gamma.lt.1.e-5) then
965	print*,'gamma =',gamma,'<1e-5'
966	stop
967	endif
968	gamma=log(0.5)/log(gamma)
969	if (gamma4) then
970	gamma=min(gamma,4.)
971	endif
972	print*,'gamma=',gamma
973	endif
974	ENDIF !first
975
976	do j=1,pjm
977	do i=1,pim
978	if (typ.eq.1) then
979	dxdy_=dxdys(i,j)
980	zlat=rlatu(j)*180./pi
981	elseif (typ.eq.2) then
982	dxdy_=dxdyu(i,j)
983	zlat=rlatu(j)*180./pi
984	elseif (typ.eq.3) then
985	dxdy_=dxdyv(i,j)
986	zlat=rlatv(j)*180./pi
987	endif
988	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
989	! pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
990	alpha(i,j)=alphamin
991	else
992	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
993	xi=min(xi,1.)
994	if(lat_min_g.le.zlat .and. zlat.le.lat_max_g) then
995	alpha(i,j)=xialphamin+(1.-xi)alphamax
996	else
997	alpha(i,j)=0.
998	endif
999	endif
1000	enddo
1001	enddo
1002	ENDIF ! guide_reg
1003
1004	END SUBROUTINE tau2alpha
1005
1006	!=======================================================================
1007	SUBROUTINE guide_read(timestep)
1008
1009	IMPLICIT NONE
1010
1011	#include "netcdf.inc"
1012	#include "dimensions.h"
1013	#include "paramet.h"
1014
1015	INTEGER, INTENT(IN) :: timestep
1016
1017	LOGICAL, SAVE :: first=.TRUE.
1018	! Identification fichiers et variables NetCDF:
1019	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1020	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1021	INTEGER :: ncidpl,varidpl,varidap,varidbp
1022	! Variables auxiliaires NetCDF:
1023	INTEGER, DIMENSION(4) :: start,count
1024	INTEGER :: status,rcode
1025
1026	CHARACTER (len = 80) :: abort_message
1027	CHARACTER (len = 20) :: modname = 'guide_read'
1028	! -----------------------------------------------------------------
1029	! Premier appel: initialisation de la lecture des fichiers
1030	! -----------------------------------------------------------------
1031	if (first) then
1032	ncidpl=-99
1033	print*,'Guide: ouverture des fichiers guidage '
1034	! Niveaux de pression si non constants
1035	if (guide_modele) then
1036	print *,'Lecture du guidage sur niveaux modele'
1037	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1038	IF (rcode.NE.NF_NOERR) THEN
1039	print *,'Guide: probleme -> pas de fichier apbp.nc'
1040	CALL abort_gcm(modname,abort_message,1)
1041	ENDIF
1042	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1043	IF (rcode.NE.NF_NOERR) THEN
1044	print *,'Guide: probleme -> pas de variable AP, fichier apbp.nc'
1045	CALL abort_gcm(modname,abort_message,1)
1046	ENDIF
1047	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1048	IF (rcode.NE.NF_NOERR) THEN
1049	print *,'Guide: probleme -> pas de variable BP, fichier apbp.nc'
1050	CALL abort_gcm(modname,abort_message,1)
1051	ENDIF
1052	print*,'ncidpl,varidap',ncidpl,varidap
1053	endif
1054	! Vent zonal
1055	if (guide_u) then
1056	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1057	IF (rcode.NE.NF_NOERR) THEN
1058	print *,'Guide: probleme -> pas de fichier u.nc'
1059	CALL abort_gcm(modname,abort_message,1)
1060	ENDIF
1061	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1062	IF (rcode.NE.NF_NOERR) THEN
1063	print *,'Guide: probleme -> pas de variable UWND, fichier u.nc'
1064	CALL abort_gcm(modname,abort_message,1)
1065	ENDIF
1066	print*,'ncidu,varidu',ncidu,varidu
1067	if (ncidpl.eq.-99) ncidpl=ncidu
1068	endif
1069	! Vent meridien
1070	if (guide_v) then
1071	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1072	IF (rcode.NE.NF_NOERR) THEN
1073	print *,'Guide: probleme -> pas de fichier v.nc'
1074	CALL abort_gcm(modname,abort_message,1)
1075	ENDIF
1076	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1077	IF (rcode.NE.NF_NOERR) THEN
1078	print *,'Guide: probleme -> pas de variable VWND, fichier v.nc'
1079	CALL abort_gcm(modname,abort_message,1)
1080	ENDIF
1081	print*,'ncidv,varidv',ncidv,varidv
1082	if (ncidpl.eq.-99) ncidpl=ncidv
1083	endif
1084	! Temperature
1085	if (guide_T) then
1086	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1087	IF (rcode.NE.NF_NOERR) THEN
1088	print *,'Guide: probleme -> pas de fichier T.nc'
1089	CALL abort_gcm(modname,abort_message,1)
1090	ENDIF
1091	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1092	IF (rcode.NE.NF_NOERR) THEN
1093	print *,'Guide: probleme -> pas de variable AIR, fichier T.nc'
1094	CALL abort_gcm(modname,abort_message,1)
1095	ENDIF
1096	print*,'ncidT,varidT',ncidt,varidt
1097	if (ncidpl.eq.-99) ncidpl=ncidt
1098	endif
1099	! Humidite
1100	if (guide_Q) then
1101	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1102	IF (rcode.NE.NF_NOERR) THEN
1103	print *,'Guide: probleme -> pas de fichier hur.nc'
1104	CALL abort_gcm(modname,abort_message,1)
1105	ENDIF
1106	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1107	IF (rcode.NE.NF_NOERR) THEN
1108	print *,'Guide: probleme -> pas de variable RH, fichier hur.nc'
1109	CALL abort_gcm(modname,abort_message,1)
1110	ENDIF
1111	print*,'ncidQ,varidQ',ncidQ,varidQ
1112	if (ncidpl.eq.-99) ncidpl=ncidQ
1113	endif
1114	! Pression de surface
1115	if ((guide_P).OR.(guide_modele)) then
1116	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1117	IF (rcode.NE.NF_NOERR) THEN
1118	print *,'Guide: probleme -> pas de fichier ps.nc'
1119	CALL abort_gcm(modname,abort_message,1)
1120	ENDIF
1121	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1122	IF (rcode.NE.NF_NOERR) THEN
1123	print *,'Guide: probleme -> pas de variable SP, fichier ps.nc'
1124	CALL abort_gcm(modname,abort_message,1)
1125	ENDIF
1126	print*,'ncidps,varidps',ncidps,varidps
1127	endif
1128	! Coordonnee verticale
1129	if (.not.guide_modele) then
1130	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1131	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1132	print*,'ncidpl,varidpl',ncidpl,varidpl
1133	endif
1134	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1135	if (guide_modele) then
1136	#ifdef NC_DOUBLE
1137	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1138	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1139	#else
1140	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1141	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1142	#endif
1143	else
1144	#ifdef NC_DOUBLE
1145	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1146	#else
1147	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1148	#endif
1149	apnc=apnc*100.! conversion en Pascals
1150	bpnc(:)=0.
1151	endif
1152	first=.FALSE.
1153	endif ! (first)
1154
1155	! -----------------------------------------------------------------
1156	! lecture des champs u, v, T, Q, ps
1157	! -----------------------------------------------------------------
1158
1159	! dimensions pour les champs scalaires et le vent zonal
1160	start(1)=1
1161	start(2)=1
1162	start(3)=1
1163	start(4)=timestep
1164
1165	count(1)=iip1
1166	count(2)=jjp1
1167	count(3)=nlevnc
1168	count(4)=1
1169
1170	! Vent zonal
1171	if (guide_u) then
1172	#ifdef NC_DOUBLE
1173	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,unat2)
1174	#else
1175	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,unat2)
1176	#endif
1177	IF (invert_y) THEN
1178	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1179	ENDIF
1180	endif
1181
1182	! Temperature
1183	if (guide_T) then
1184	#ifdef NC_DOUBLE
1185	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,tnat2)
1186	#else
1187	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,tnat2)
1188	#endif
1189	IF (invert_y) THEN
1190	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1191	ENDIF
1192	endif
1193
1194	! Humidite
1195	if (guide_Q) then
1196	#ifdef NC_DOUBLE
1197	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,qnat2)
1198	#else
1199	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,qnat2)
1200	#endif
1201	IF (invert_y) THEN
1202	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1203	ENDIF
1204
1205	endif
1206
1207	! Vent meridien
1208	if (guide_v) then
1209	count(2)=jjm
1210	#ifdef NC_DOUBLE
1211	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,vnat2)
1212	#else
1213	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,vnat2)
1214	#endif
1215	IF (invert_y) THEN
1216	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1217	ENDIF
1218	endif
1219
1220	! Pression de surface
1221	if ((guide_P).OR.(guide_modele)) then
1222	start(3)=timestep
1223	start(4)=0
1224	count(2)=jjp1
1225	count(3)=1
1226	count(4)=0
1227	#ifdef NC_DOUBLE
1228	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,psnat2)
1229	#else
1230	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,psnat2)
1231	#endif
1232	IF (invert_y) THEN
1233	CALL invert_lat(iip1,jjp1,1,psnat2)
1234	ENDIF
1235	endif
1236
1237	END SUBROUTINE guide_read
1238
1239	!=======================================================================
1240	SUBROUTINE guide_read2D(timestep)
1241
1242	IMPLICIT NONE
1243
1244	#include "netcdf.inc"
1245	#include "dimensions.h"
1246	#include "paramet.h"
1247
1248	INTEGER, INTENT(IN) :: timestep
1249
1250	LOGICAL, SAVE :: first=.TRUE.
1251	! Identification fichiers et variables NetCDF:
1252	INTEGER, SAVE :: ncidu,varidu,ncidv,varidv,ncidQ
1253	INTEGER, SAVE :: varidQ,ncidt,varidt,ncidps,varidps
1254	INTEGER :: ncidpl,varidpl,varidap,varidbp
1255	! Variables auxiliaires NetCDF:
1256	INTEGER, DIMENSION(4) :: start,count
1257	INTEGER :: status,rcode
1258	! Variables for 3D extension:
1259	REAL, DIMENSION (jjp1,llm) :: zu
1260	REAL, DIMENSION (jjm,llm) :: zv
1261	INTEGER :: i
1262
1263	CHARACTER (len = 80) :: abort_message
1264	CHARACTER (len = 20) :: modname = 'guide_read2D'
1265	! -----------------------------------------------------------------
1266	! Premier appel: initialisation de la lecture des fichiers
1267	! -----------------------------------------------------------------
1268	if (first) then
1269	ncidpl=-99
1270	print*,'Guide: ouverture des fichiers guidage '
1271	! Niveaux de pression si non constants
1272	if (guide_modele) then
1273	print *,'Lecture du guidage sur niveaux modele'
1274	rcode = nf90_open('apbp.nc', nf90_nowrite, ncidpl)
1275	IF (rcode.NE.NF_NOERR) THEN
1276	print *,'Guide: probleme -> pas de fichier apbp.nc'
1277	CALL abort_gcm(modname,abort_message,1)
1278	ENDIF
1279	rcode = nf90_inq_varid(ncidpl, 'AP', varidap)
1280	IF (rcode.NE.NF_NOERR) THEN
1281	print *,'Guide: probleme -> pas de variable AP, fichier apbp.nc'
1282	CALL abort_gcm(modname,abort_message,1)
1283	ENDIF
1284	rcode = nf90_inq_varid(ncidpl, 'BP', varidbp)
1285	IF (rcode.NE.NF_NOERR) THEN
1286	print *,'Guide: probleme -> pas de variable BP, fichier apbp.nc'
1287	CALL abort_gcm(modname,abort_message,1)
1288	ENDIF
1289	print*,'ncidpl,varidap',ncidpl,varidap
1290	endif
1291	! Vent zonal
1292	if (guide_u) then
1293	rcode = nf90_open('u.nc', nf90_nowrite, ncidu)
1294	IF (rcode.NE.NF_NOERR) THEN
1295	print *,'Guide: probleme -> pas de fichier u.nc'
1296	CALL abort_gcm(modname,abort_message,1)
1297	ENDIF
1298	rcode = nf90_inq_varid(ncidu, 'UWND', varidu)
1299	IF (rcode.NE.NF_NOERR) THEN
1300	print *,'Guide: probleme -> pas de variable UWND, fichier u.nc'
1301	CALL abort_gcm(modname,abort_message,1)
1302	ENDIF
1303	print*,'ncidu,varidu',ncidu,varidu
1304	if (ncidpl.eq.-99) ncidpl=ncidu
1305	endif
1306	! Vent meridien
1307	if (guide_v) then
1308	rcode = nf90_open('v.nc', nf90_nowrite, ncidv)
1309	IF (rcode.NE.NF_NOERR) THEN
1310	print *,'Guide: probleme -> pas de fichier v.nc'
1311	CALL abort_gcm(modname,abort_message,1)
1312	ENDIF
1313	rcode = nf90_inq_varid(ncidv, 'VWND', varidv)
1314	IF (rcode.NE.NF_NOERR) THEN
1315	print *,'Guide: probleme -> pas de variable VWND, fichier v.nc'
1316	CALL abort_gcm(modname,abort_message,1)
1317	ENDIF
1318	print*,'ncidv,varidv',ncidv,varidv
1319	if (ncidpl.eq.-99) ncidpl=ncidv
1320	endif
1321	! Temperature
1322	if (guide_T) then
1323	rcode = nf90_open('T.nc', nf90_nowrite, ncidt)
1324	IF (rcode.NE.NF_NOERR) THEN
1325	print *,'Guide: probleme -> pas de fichier T.nc'
1326	CALL abort_gcm(modname,abort_message,1)
1327	ENDIF
1328	rcode = nf90_inq_varid(ncidt, 'AIR', varidt)
1329	IF (rcode.NE.NF_NOERR) THEN
1330	print *,'Guide: probleme -> pas de variable AIR, fichier T.nc'
1331	CALL abort_gcm(modname,abort_message,1)
1332	ENDIF
1333	print*,'ncidT,varidT',ncidt,varidt
1334	if (ncidpl.eq.-99) ncidpl=ncidt
1335	endif
1336	! Humidite
1337	if (guide_Q) then
1338	rcode = nf90_open('hur.nc', nf90_nowrite, ncidQ)
1339	IF (rcode.NE.NF_NOERR) THEN
1340	print *,'Guide: probleme -> pas de fichier hur.nc'
1341	CALL abort_gcm(modname,abort_message,1)
1342	ENDIF
1343	rcode = nf90_inq_varid(ncidQ, 'RH', varidQ)
1344	IF (rcode.NE.NF_NOERR) THEN
1345	print *,'Guide: probleme -> pas de variable RH, fichier hur.nc'
1346	CALL abort_gcm(modname,abort_message,1)
1347	ENDIF
1348	print*,'ncidQ,varidQ',ncidQ,varidQ
1349	if (ncidpl.eq.-99) ncidpl=ncidQ
1350	endif
1351	! Pression de surface
1352	if ((guide_P).OR.(guide_modele)) then
1353	rcode = nf90_open('ps.nc', nf90_nowrite, ncidps)
1354	IF (rcode.NE.NF_NOERR) THEN
1355	print *,'Guide: probleme -> pas de fichier ps.nc'
1356	CALL abort_gcm(modname,abort_message,1)
1357	ENDIF
1358	rcode = nf90_inq_varid(ncidps, 'SP', varidps)
1359	IF (rcode.NE.NF_NOERR) THEN
1360	print *,'Guide: probleme -> pas de variable SP, fichier ps.nc'
1361	CALL abort_gcm(modname,abort_message,1)
1362	ENDIF
1363	print*,'ncidps,varidps',ncidps,varidps
1364	endif
1365	! Coordonnee verticale
1366	if (.not.guide_modele) then
1367	rcode = nf90_inq_varid(ncidpl, 'LEVEL', varidpl)
1368	IF (rcode.NE.0) rcode = nf90_inq_varid(ncidpl, 'PRESSURE', varidpl)
1369	print*,'ncidpl,varidpl',ncidpl,varidpl
1370	endif
1371	! Coefs ap, bp pour calcul de la pression aux differents niveaux
1372	if (guide_modele) then
1373	#ifdef NC_DOUBLE
1374	status=NF_GET_VARA_DOUBLE(ncidpl,varidap,1,nlevnc,apnc)
1375	status=NF_GET_VARA_DOUBLE(ncidpl,varidbp,1,nlevnc,bpnc)
1376	#else
1377	status=NF_GET_VARA_REAL(ncidpl,varidap,1,nlevnc,apnc)
1378	status=NF_GET_VARA_REAL(ncidpl,varidbp,1,nlevnc,bpnc)
1379	#endif
1380	else
1381	#ifdef NC_DOUBLE
1382	status=NF_GET_VARA_DOUBLE(ncidpl,varidpl,1,nlevnc,apnc)
1383	#else
1384	status=NF_GET_VARA_REAL(ncidpl,varidpl,1,nlevnc,apnc)
1385	#endif
1386	apnc=apnc*100.! conversion en Pascals
1387	bpnc(:)=0.
1388	endif
1389	first=.FALSE.
1390	endif ! (first)
1391
1392	! -----------------------------------------------------------------
1393	! lecture des champs u, v, T, Q, ps
1394	! -----------------------------------------------------------------
1395
1396	! dimensions pour les champs scalaires et le vent zonal
1397	start(1)=1
1398	start(2)=1
1399	start(3)=1
1400	start(4)=timestep
1401
1402	count(1)=1
1403	count(2)=jjp1
1404	count(3)=nlevnc
1405	count(4)=1
1406
1407	! Vent zonal
1408	if (guide_u) then
1409	#ifdef NC_DOUBLE
1410	status=NF_GET_VARA_DOUBLE(ncidu,varidu,start,count,zu)
1411	#else
1412	status=NF_GET_VARA_REAL(ncidu,varidu,start,count,zu)
1413	#endif
1414	DO i=1,iip1
1415	unat2(i,:,:)=zu(:,:)
1416	ENDDO
1417
1418	IF (invert_y) THEN
1419	CALL invert_lat(iip1,jjp1,nlevnc,unat2)
1420	ENDIF
1421
1422	endif
1423
1424	! Temperature
1425	if (guide_T) then
1426	#ifdef NC_DOUBLE
1427	status=NF_GET_VARA_DOUBLE(ncidt,varidt,start,count,zu)
1428	#else
1429	status=NF_GET_VARA_REAL(ncidt,varidt,start,count,zu)
1430	#endif
1431	DO i=1,iip1
1432	tnat2(i,:,:)=zu(:,:)
1433	ENDDO
1434
1435	IF (invert_y) THEN
1436	CALL invert_lat(iip1,jjp1,nlevnc,tnat2)
1437	ENDIF
1438
1439	endif
1440
1441	! Humidite
1442	if (guide_Q) then
1443	#ifdef NC_DOUBLE
1444	status=NF_GET_VARA_DOUBLE(ncidQ,varidQ,start,count,zu)
1445	#else
1446	status=NF_GET_VARA_REAL(ncidQ,varidQ,start,count,zu)
1447	#endif
1448	DO i=1,iip1
1449	qnat2(i,:,:)=zu(:,:)
1450	ENDDO
1451
1452	IF (invert_y) THEN
1453	CALL invert_lat(iip1,jjp1,nlevnc,qnat2)
1454	ENDIF
1455
1456	endif
1457
1458	! Vent meridien
1459	if (guide_v) then
1460	count(2)=jjm
1461	#ifdef NC_DOUBLE
1462	status=NF_GET_VARA_DOUBLE(ncidv,varidv,start,count,zv)
1463	#else
1464	status=NF_GET_VARA_REAL(ncidv,varidv,start,count,zv)
1465	#endif
1466	DO i=1,iip1
1467	vnat2(i,:,:)=zv(:,:)
1468	ENDDO
1469
1470	IF (invert_y) THEN
1471	CALL invert_lat(iip1,jjm,nlevnc,vnat2)
1472	ENDIF
1473
1474	endif
1475
1476	! Pression de surface
1477	if ((guide_P).OR.(guide_modele)) then
1478	start(3)=timestep
1479	start(4)=0
1480	count(2)=jjp1
1481	count(3)=1
1482	count(4)=0
1483	#ifdef NC_DOUBLE
1484	status=NF_GET_VARA_DOUBLE(ncidps,varidps,start,count,zu(:,1))
1485	#else
1486	status=NF_GET_VARA_REAL(ncidps,varidps,start,count,zu(:,1))
1487	#endif
1488	DO i=1,iip1
1489	psnat2(i,:)=zu(:,1)
1490	ENDDO
1491
1492	IF (invert_y) THEN
1493	CALL invert_lat(iip1,jjp1,1,psnat2)
1494	ENDIF
1495
1496	endif
1497
1498	END SUBROUTINE guide_read2D
1499
1500	!=======================================================================
1501	SUBROUTINE guide_out(varname,hsize,vsize,field)
1502
1503	IMPLICIT NONE
1504
1505	INCLUDE "dimensions.h"
1506	INCLUDE "paramet.h"
1507	INCLUDE "netcdf.inc"
1508	INCLUDE "comgeom2.h"
1509	INCLUDE "comconst.h"
1510	INCLUDE "comvert.h"
1511
1512	! Variables entree
1513	CHARACTER(), INTENT(IN) :: varname
1514	INTEGER, INTENT (IN) :: hsize,vsize
1515	REAL, DIMENSION (iip1,hsize,vsize), INTENT(IN) :: field
1516
1517	! Variables locales
1518	INTEGER, SAVE :: timestep=0
1519	! Identites fichier netcdf
1520	INTEGER :: nid, id_lonu, id_lonv, id_latu, id_latv, id_tim, id_lev
1521	INTEGER :: vid_lonu,vid_lonv,vid_latu,vid_latv,vid_cu,vid_cv,vid_lev
1522	INTEGER :: vid_au,vid_av
1523	INTEGER, DIMENSION (3) :: dim3
1524	INTEGER, DIMENSION (4) :: dim4,count,start
1525	INTEGER :: ierr, varid,l
1526	REAL, DIMENSION (iip1,hsize,vsize) :: field2
1527
1528	print *,'Guide: output timestep',timestep,'var ',varname
1529	IF (timestep.EQ.0) THEN
1530	! ----------------------------------------------
1531	! initialisation fichier de sortie
1532	! ----------------------------------------------
1533	! Ouverture du fichier
1534	ierr=NF_CREATE("guide_ins.nc",NF_CLOBBER,nid)
1535	! Definition des dimensions
1536	ierr=NF_DEF_DIM(nid,"LONU",iip1,id_lonu)
1537	ierr=NF_DEF_DIM(nid,"LONV",iip1,id_lonv)
1538	ierr=NF_DEF_DIM(nid,"LATU",jjp1,id_latu)
1539	ierr=NF_DEF_DIM(nid,"LATV",jjm,id_latv)
1540	ierr=NF_DEF_DIM(nid,"LEVEL",llm,id_lev)
1541	ierr=NF_DEF_DIM(nid,"TIME",NF_UNLIMITED,id_tim)
1542
1543	! Creation des variables dimensions
1544	ierr=NF_DEF_VAR(nid,"LONU",NF_FLOAT,1,id_lonu,vid_lonu)
1545	ierr=NF_DEF_VAR(nid,"LONV",NF_FLOAT,1,id_lonv,vid_lonv)
1546	ierr=NF_DEF_VAR(nid,"LATU",NF_FLOAT,1,id_latu,vid_latu)
1547	ierr=NF_DEF_VAR(nid,"LATV",NF_FLOAT,1,id_latv,vid_latv)
1548	ierr=NF_DEF_VAR(nid,"LEVEL",NF_FLOAT,1,id_lev,vid_lev)
1549	ierr=NF_DEF_VAR(nid,"cu",NF_FLOAT,2,(/id_lonu,id_latu/),vid_cu)
1550	ierr=NF_DEF_VAR(nid,"au",NF_FLOAT,2,(/id_lonu,id_latu/),vid_au)
1551	ierr=NF_DEF_VAR(nid,"cv",NF_FLOAT,2,(/id_lonv,id_latv/),vid_cv)
1552	ierr=NF_DEF_VAR(nid,"av",NF_FLOAT,2,(/id_lonv,id_latv/),vid_av)
1553
1554	ierr=NF_ENDDEF(nid)
1555
1556	! Enregistrement des variables dimensions
1557	#ifdef NC_DOUBLE
1558	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonu,rlonu*180./pi)
1559	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lonv,rlonv*180./pi)
1560	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latu,rlatu*180./pi)
1561	ierr = NF_PUT_VAR_DOUBLE(nid,vid_latv,rlatv*180./pi)
1562	ierr = NF_PUT_VAR_DOUBLE(nid,vid_lev,presnivs)
1563	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cu,cu)
1564	ierr = NF_PUT_VAR_DOUBLE(nid,vid_cv,cv)
1565	ierr = NF_PUT_VAR_DOUBLE(nid,vid_au,alpha_u)
1566	ierr = NF_PUT_VAR_DOUBLE(nid,vid_av,alpha_v)
1567	#else
1568	ierr = NF_PUT_VAR_REAL(nid,vid_lonu,rlonu*180./pi)
1569	ierr = NF_PUT_VAR_REAL(nid,vid_lonv,rlonv*180./pi)
1570	ierr = NF_PUT_VAR_REAL(nid,vid_latu,rlatu*180./pi)
1571	ierr = NF_PUT_VAR_REAL(nid,vid_latv,rlatv*180./pi)
1572	ierr = NF_PUT_VAR_REAL(nid,vid_lev,presnivs)
1573	ierr = NF_PUT_VAR_REAL(nid,vid_cu,cu)
1574	ierr = NF_PUT_VAR_REAL(nid,vid_cv,cv)
1575	ierr = NF_PUT_VAR_REAL(nid,vid_au,alpha_u)
1576	ierr = NF_PUT_VAR_REAL(nid,vid_av,alpha_v)
1577	#endif
1578	! --------------------------------------------------------------------
1579	! Cr�ation des variables sauvegard�es
1580	! --------------------------------------------------------------------
1581	ierr = NF_REDEF(nid)
1582	! Surface pressure (GCM)
1583	dim3=(/id_lonv,id_latu,id_tim/)
1584	ierr = NF_DEF_VAR(nid,"SP",NF_FLOAT,3,dim3,varid)
1585	! Surface pressure (guidage)
1586	IF (guide_P) THEN
1587	dim3=(/id_lonv,id_latu,id_tim/)
1588	ierr = NF_DEF_VAR(nid,"ps",NF_FLOAT,3,dim3,varid)
1589	ENDIF
1590	! Zonal wind
1591	IF (guide_u) THEN
1592	dim4=(/id_lonu,id_latu,id_lev,id_tim/)
1593	ierr = NF_DEF_VAR(nid,"u",NF_FLOAT,4,dim4,varid)
1594	ierr = NF_DEF_VAR(nid,"ua",NF_FLOAT,4,dim4,varid)
1595	ierr = NF_DEF_VAR(nid,"ucov",NF_FLOAT,4,dim4,varid)
1596	ENDIF
1597	! Merid. wind
1598	IF (guide_v) THEN
1599	dim4=(/id_lonv,id_latv,id_lev,id_tim/)
1600	ierr = NF_DEF_VAR(nid,"v",NF_FLOAT,4,dim4,varid)
1601	ierr = NF_DEF_VAR(nid,"va",NF_FLOAT,4,dim4,varid)
1602	ierr = NF_DEF_VAR(nid,"vcov",NF_FLOAT,4,dim4,varid)
1603	ENDIF
1604	! Pot. Temperature
1605	IF (guide_T) THEN
1606	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1607	ierr = NF_DEF_VAR(nid,"teta",NF_FLOAT,4,dim4,varid)
1608	ENDIF
1609	! Specific Humidity
1610	IF (guide_Q) THEN
1611	dim4=(/id_lonv,id_latu,id_lev,id_tim/)
1612	ierr = NF_DEF_VAR(nid,"q",NF_FLOAT,4,dim4,varid)
1613	ENDIF
1614
1615	ierr = NF_ENDDEF(nid)
1616	ierr = NF_CLOSE(nid)
1617	ENDIF ! timestep=0
1618
1619	! --------------------------------------------------------------------
1620	! Enregistrement du champ
1621	! --------------------------------------------------------------------
1622	ierr=NF_OPEN("guide_ins.nc",NF_WRITE,nid)
1623
1624	IF (varname=="SP") timestep=timestep+1
1625
1626	ierr = NF_INQ_VARID(nid,varname,varid)
1627	SELECT CASE (varname)
1628	CASE ("SP","ps")
1629	start=(/1,1,timestep,0/)
1630	count=(/iip1,jjp1,1,0/)
1631	CASE ("v","va","vcov")
1632	start=(/1,1,1,timestep/)
1633	count=(/iip1,jjm,llm,1/)
1634	CASE DEFAULT
1635	start=(/1,1,1,timestep/)
1636	count=(/iip1,jjp1,llm,1/)
1637	END SELECT
1638
1639	SELECT CASE (varname)
1640	CASE("u","ua")
1641	DO l=1,llm ; field2(:,2:jjm,l)=field(:,2:jjm,l)/cu(:,2:jjm) ; ENDDO
1642	field2(:,1,:)=0. ; field2(:,jjp1,:)=0.
1643	CASE("v","va")
1644	DO l=1,llm ; field2(:,:,l)=field(:,:,l)/cv(:,:) ; ENDDO
1645	CASE DEFAULT
1646	field2=field
1647	END SELECT
1648
1649
1650	#ifdef NC_DOUBLE
1651	ierr = NF_PUT_VARA_DOUBLE(nid,varid,start,count,field2)
1652	#else
1653	ierr = NF_PUT_VARA_REAL(nid,varid,start,count,field2)
1654	#endif
1655
1656	ierr = NF_CLOSE(nid)
1657
1658	END SUBROUTINE guide_out
1659
1660
1661	!===========================================================================
1662	subroutine correctbid(iim,nl,x)
1663	integer iim,nl
1664	real x(iim+1,nl)
1665	integer i,l
1666	real zz
1667
1668	do l=1,nl
1669	do i=2,iim-1
1670	if(abs(x(i,l)).gt.1.e10) then
1671	zz=0.5*(x(i-1,l)+x(i+1,l))
1672	print*,'correction ',i,l,x(i,l),zz
1673	x(i,l)=zz
1674	endif
1675	enddo
1676	enddo
1677	return
1678	end subroutine correctbid
1679
1680	!===========================================================================
1681	END MODULE guide_mod

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

Download in other formats:

Original Format