Context Navigation

guide_mod.F90 @ 1186

Last change on this file since 1186 was 1186, checked in by Ehouarn Millour, 15 years ago

Cleanup around IOIPSL, so that LMDZ dynamics may be used without IOIPSL.

moved ersatz IOIPSL routines (ioipsl_* , taken from IOIPSLv2_1_8, so that 'getin' function may be used even if not using the IOIPSL library) from dyn3d/dyn3dpar to bibio.
enclosed 'use ioipsl' instruction with #ifdef CPP_IOIPSL cpp keys.

Property svn:executable set to *

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

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

Download in other formats:

Original Format