Context Navigation

guide_p_mod.F90 @ 1243

Last change on this file since 1243 was 1019, checked in by emillour, 11 years ago

Common dynamics; keep up with updates (seq and ) in LMDZ5 (up tio rev 1845):

General stuff:
makelmdz_fcm: add options -j # (compile using # threads) and -full, and to keep up

with Earth model, possibility to compile with various versions of orchidee

bld.cfg: adaptations to enable compiling using multiple threads
build_gcm: adaptations to enable compiling using multiple threads
makelmdz: keep up with Earth model: possibility to compile with various versions of orchidee + cosmetic changes + library directory name change

bibio:
wxios.F90 : Added for possible future use of XIOS library

filtrez:
mkl_dft_type.f90 & mkl_dfti.f90 : MKL (for MKL FFT) interface definitions
filtreg_mod : limit use of FFT to parallel mode
mod_filtre_fft.F90 & mod_filtre_fft_lov.F90 : swich to use parallel_lmdz

dyn3d:
abort_gcm.F : add things for xios
advtrac.F90 : minor change in CFL outputs
ce0l.F90 : indicesol.h is now module indice_sol_mod
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
inigeom.F : move two computations outside loop

dyn3dpar:
parallel.F90 => parallel_lmdz.F90 : and change all the "use parallel" into "use parallel_lmdz" in all files in dyn3dpar
comvert.h : cosmetic change on comments
gcm.F : add xios and use module indice_sol_mod (for INCA)
leapfrog_p.F : add xios + correction for times in Newtonian case
ce0l.F90 : indicesol.h is now module indice_sol_mod
inigeom.F : move two computations outside loop

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

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

Download in other formats:

Original Format