Context Navigation

guide_mod.F90 @ 1293

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

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

Download in other formats:

Original Format