Context Navigation

guide_mod.F90 @ 4248

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

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

Download in other formats:

Original Format