Context Navigation

guide_mod.f90 @ 5273

Last change on this file since 5273 was 5272, checked in by abarral, 2 days ago
Turn paramet.h into a module
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: 66.5 KB

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

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

Download in other formats:

Original Format