Context Navigation

guide_mod.F90 @ 4233

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

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

Download in other formats:

Original Format