Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

guide_mod.F90

Last change on this file was 4470, checked in by Laurent Fairhead, 15 months ago
Replaced STOP instructions by calls to abort_gcm for a cleaner exit
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