Context Navigation

guide_mod.F90 @ 5127

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

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

Download in other formats:

Original Format