Context Navigation

guide_mod.F90 @ 1907

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

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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

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

Download in other formats:

Original Format