Context Navigation

guide_p_mod.F90 @ 127

Last change on this file since 127 was 127, checked in by emillour, 14 years ago

Ehouarn: suite de l'implementation de la discretisation verticale,
avec quelques mises a jour pour concorder avec la version terrestre.
-> Finalement, on met un flag "disvert_type" pour fixer la discretisation

disvert_type==1 (defaut si planet_type=="earth") pour cas terrestre
disvert_type==2 (defaut si planet_type!="earth") pour cas planeto (z2sig.def)

-> au passage, pour rester en phase avec modele terrestre on renomme

disvert_terre en disvert (le disvert "alternatif" demeure 'disvert_noterre')

File size: 64.1 KB

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

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

Download in other formats:

Original Format