Context Navigation

guide_mod.F90 @ 1270

Last change on this file since 1270 was 1188, checked in by lguez, 15 years ago
Translated calls using NetCDF 2.4 interface to calls using NetCDF 3.6 Fortran 90 interface. Corrected a few comments in "output.def".
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