Context Navigation

guide_loc_mod.F90 @ 4054

Last change on this file since 4054 was 4054, checked in by lguez, 2 years ago

Only compute pk if we are going to need it

In revsion revision 4042, wa added the computation of pk for
guide_plevs == 1. We can restrict a little better the cases where we
need pk.

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: 86.3 KB

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

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

Download in other formats:

Original Format