Context Navigation

guide_mod.F90 @ 2023

Last change on this file since 2023 was 2021, checked in by lguez, 11 years ago

Removed unused variables pks, pk, pkf from main program unit gcm.

Encapsulated procedures exner_hyb, exner_hyb_p, exner_hyb_loc,
exner_milieu, exner_milieu_p and exner_milieu_loc into
modules. (Compulsory to allow optional arguments.)

In the procedures exner_hyb, exner_hyb_p, exner_hyb_loc, donwgraded
arguments alpha and beta to local variables. In exner_milieu,
exner_milieu_p and exner_milieu_loc, removed beta altogether. In the
six procedures exner_*, made pkf an optional argument. Made some
cosmetic modifications in order to keep the six procedures exner_* as
close as possible.

In the six procedures exner_*, removed the averaging of pks at the
poles: this is not useful because ps is already the same at all
longitudes at the poles. This modification changes the results of the
program. Motivation: had to do this for exner_hyb because we call it
from test_disvert with a few surface pressure values.

In all the procedures calling exner_*, removed the variables alpha and
beta. Also removed variables alpha and beta from module leapfrog_mod
and from module call_calfis_mod.

Removed actual argument pkf in call to exner_hyb* and exner_milieu*
from guide_interp, guide_main, iniacademic and iniacademic_loc (pkf
was not used in those procedures).

Argument workvar of startget_dyn is used only if varname is tpot or

When varname is tpot or q, the actual argument associated to

workvar in etat0_netcdf is not y. So y in etat0_netcdf is a
place-holder, never used. So we remove optional argument y in the
calls to exner_hyb and exner_milieu from etat0_netcdf.

Created procedure test_disvert, called only by etat0_netcdf. This
procedure tests the order of pressure values at half-levels and full
levels.

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

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

Download in other formats:

Original Format