Context Navigation

guide_mod.F90 @ 3617

Last change on this file since 3617 was 1520, checked in by Ehouarn Millour, 13 years ago

Implementation of a different vertical discretization (from/for planets, but
can in principle also be used for Earth).
Choice of vertical discretization is set by flag 'disvert_type';
'disvert_type=1' is Earth standard (default; ie set to 1 if
planet_type=="earth") case.
With 'disvert_type=2', approximate altitude of layers and reference atmospheric
scale height must be given using an input file ("z2sig.def", first line
should give scale height, in km, following lines must specify the altitude,
in km above surface, of mid-layers, one per line; see disvert_noterre.F).

Checked that these changes do not impact on 'bench' results, on Vargas.

EM.

Property svn:executable set to *

File size: 52.1 KB

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

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

Download in other formats:

Original Format