Context Navigation

guide_mod.F90 @ 1333

Last change on this file since 1333 was 1299, checked in by Laurent Fairhead, 14 years ago

Nettoyage general pour se rapprocher des normes et éviter des erreurs a la
compilation:

tous les FLOAT() sont remplacés par des REAL()
tous les STOP dans phylmd sont remplacés par des appels à abort_gcm
le common control défini dans le fichier control.h est remplacé par le module control_mod pour éviter des messages sur l'alignement des variables dans les déclarations
des $Header$ remplacés par des $Id$ pour svn

Quelques remplacements à faire ont pu m'échapper

General cleanup of the code to try and adhere to norms and to prevent some
compilation errors:

all FLOAT() instructions have been replaced by REAL() instructions
all STOP instructions in phylmd have been replaced by calls to abort_gcm
the common block control defined in the control.h file has been replaced by the control_mod to prevent compilation warnings on the alignement of declared variables
$Header$ replaced by $Id$ for svn

Some changes which should have been made might have escaped me

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

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

Download in other formats:

Original Format