Context Navigation

guide_loc_mod.F90 @ 1910

Last change on this file since 1910 was 1910, checked in by Laurent Fairhead, 10 years ago
Merged trunk changes r1860:1909 into testing branch
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory
File size: 79.1 KB

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

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

Download in other formats:

Original Format