Context Navigation

guide_mod.F90 @ 2141

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

In nudging procedures, replaced explicit Euler integration of nudged
fields by exact integration. This does not change anything if
guide_add is true, but it changes the value of alpha if guide_add is
false. We could have taken into account the variation of the nudging
field during a nudging time step. This would be a small correction. We
choose not to take it into account for the time being. Also, we add a
restriction on zonal nudging: we allow it only for a grid which is
regular in longitude. It does not seem to make sense otherwise and the
exact integration would take more programming for an irregular grid.

In the sequential version, copying the parallel versions, set
iguide_int to 1 when the input value is 0.

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

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

Download in other formats:

Original Format