Context Navigation

guide_mod.F90 @ 5116

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

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

Download in other formats:

Original Format