Context Navigation

guide_mod.F90 @ 5105

Last change on this file since 5105 was 5105, checked in by abarral, 8 weeks ago
Replace 1DUTILS.h by module lmdz_1dutils.f90 Replace 1DConv.h by module lmdz_old_1dconv.f90 (it's only used by old_* files) Convert .F to .f90 Fix gradsdef.h formatting Remove unnecessary "RETURN" at the end of functions/subroutines
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: 64.8 KB

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

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

Download in other formats:

Original Format