Context Navigation

guide_p_mod.F90 @ 2021

Last change on this file since 2021 was 2021, checked in by lguez, 10 years ago

Removed unused variables pks, pk, pkf from main program unit gcm.

Encapsulated procedures exner_hyb, exner_hyb_p, exner_hyb_loc,
exner_milieu, exner_milieu_p and exner_milieu_loc into
modules. (Compulsory to allow optional arguments.)

In the procedures exner_hyb, exner_hyb_p, exner_hyb_loc, donwgraded
arguments alpha and beta to local variables. In exner_milieu,
exner_milieu_p and exner_milieu_loc, removed beta altogether. In the
six procedures exner_*, made pkf an optional argument. Made some
cosmetic modifications in order to keep the six procedures exner_* as
close as possible.

In the six procedures exner_*, removed the averaging of pks at the
poles: this is not useful because ps is already the same at all
longitudes at the poles. This modification changes the results of the
program. Motivation: had to do this for exner_hyb because we call it
from test_disvert with a few surface pressure values.

In all the procedures calling exner_*, removed the variables alpha and
beta. Also removed variables alpha and beta from module leapfrog_mod
and from module call_calfis_mod.

Removed actual argument pkf in call to exner_hyb* and exner_milieu*
from guide_interp, guide_main, iniacademic and iniacademic_loc (pkf
was not used in those procedures).

Argument workvar of startget_dyn is used only if varname is tpot or

When varname is tpot or q, the actual argument associated to

workvar in etat0_netcdf is not y. So y in etat0_netcdf is a
place-holder, never used. So we remove optional argument y in the
calls to exner_hyb and exner_milieu from etat0_netcdf.

Created procedure test_disvert, called only by etat0_netcdf. This
procedure tests the order of pressure values at half-levels and full
levels.

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
Property svn:executable set to *

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

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

Download in other formats:

Original Format