Context Navigation

guide_mod.F90 @ 5435

Last change on this file since 5435 was 4368, checked in by lguez, 2 years ago
Sync latest trunk changes to Ocean_skin
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