Context Navigation

guide_mod.F90 @ 4399

Last change on this file since 4399 was 4256, checked in by lguez, 22 months ago

Remove include netcdf.inc

They have become useless.

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

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

Download in other formats:

Original Format