Context Navigation

guide_mod.F90 @ 5136

Last change on this file since 5136 was 5136, checked in by abarral, 8 weeks ago
Put comgeom.h, comgeom2.h into modules
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: 64.1 KB

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

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

Download in other formats:

Original Format