Context Navigation

guide_mod.F90 @ 4106

Last change on this file since 4106 was 1720, checked in by Ehouarn Millour, 12 years ago
Some cleanup around pres2lev: the three dynamical cores now use the module pres2lev_mod.F90 EM
Property svn:executable set to ``*
File size: 52.0 KB

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

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

Download in other formats:

Original Format