Context Navigation

guide_mod.F90 @ 4171

Last change on this file since 4171 was 4171, checked in by lguez, 2 years ago
Sync latest trunk changes to branch LMDZ-ECRAD.
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: 66.7 KB

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

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

Download in other formats:

Original Format