Context Navigation

guide_loc_mod.F90 @ 5151

Last change on this file since 5151 was 5136, checked in by abarral, 7 months 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: 83.9 KB

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

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

Download in other formats:

Original Format