Context Navigation

guide_mod.F90 @ 4248

Last change on this file since 4248 was 4246, checked in by evignon, 2 years 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