Context Navigation

guide_mod.f90 @ 5410

Last change on this file since 5410 was 5285, checked in by abarral, 7 weeks ago
As discussed internally, remove generic ONLY: ... for new _mod_h 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: 65.3 KB

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

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

Download in other formats:

Original Format