Context Navigation

guide_loc_mod.F90 @ 5194

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

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

Download in other formats:

Original Format