Context Navigation

guide.F @ 633

Last change on this file since 633 was 617, checked in by lmdzadmin, 20 years ago

corrections pour le mode guide :

mode non zoome
guidage de certains champs seulement

MAF

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 15.5 KB

Line
1	!
2	! $Header$
3	!
4	subroutine guide(itau,ucov,vcov,teta,q,masse,ps)
5
6	IMPLICIT NONE
7
8	c ...... Version du 10/01/98 ..........
9
10	c avec coordonnees verticales hybrides
11	c avec nouveaux operat. dissipation * ( gradiv2,divgrad2,nxgraro2 )
12
13	c=======================================================================
14	c
15	c Auteur: F.Hourdin
16	c -------
17	c
18	c Objet:
19	c ------
20	c
21	c GCM LMD nouvelle grille
22	c
23	c=======================================================================
24
25	c ... Dans inigeom , nouveaux calculs pour les elongations cu , cv
26	c et possibilite d'appeler une fonction f(y) a derivee tangente
27	c hyperbolique a la place de la fonction a derivee sinusoidale.
28
29	c ... Possibilite de choisir le shema de Van-leer pour l'advection de
30	c q , en faisant iadv = 10 dans traceur (29/04/97) .
31	c
32	c-----------------------------------------------------------------------
33	c Declarations:
34	c -------------
35
36	#include "dimensions.h"
37	#include "paramet.h"
38	#include "comconst.h"
39	#include "comdissnew.h"
40	#include "comvert.h"
41	#include "comgeom.h"
42	#include "logic.h"
43	#include "temps.h"
44	#include "control.h"
45	#include "ener.h"
46	#include "netcdf.inc"
47	#include "description.h"
48	#include "serre.h"
49	#include "tracstoke.h"
50	#include "guide.h"
51
52
53	c variables dynamiques
54	REAL vcov(ip1jm,llm),ucov(ip1jmp1,llm) ! vents covariants
55	REAL teta(ip1jmp1,llm) ! temperature potentielle
56	REAL q(ip1jmp1,llm) ! temperature potentielle
57	REAL ps(ip1jmp1) ! pression au sol
58	REAL masse(ip1jmp1,llm) ! masse d'air
59
60	c common passe pour des sorties
61	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
62	common/comdxdy/dxdys,dxdyu,dxdyv
63
64	c variables dynamiques pour les reanalyses.
65	REAL ucovrea1(ip1jmp1,llm),vcovrea1(ip1jm,llm) !vts cov reas
66	REAL tetarea1(ip1jmp1,llm) ! temp pot reales
67	REAL qrea1(ip1jmp1,llm) ! temp pot reales
68	REAL masserea1(ip1jmp1,llm) ! masse
69	REAL psrea1(ip1jmp1) ! ps
70	REAL ucovrea2(ip1jmp1,llm),vcovrea2(ip1jm,llm) !vts cov reas
71	REAL tetarea2(ip1jmp1,llm) ! temp pot reales
72	REAL qrea2(ip1jmp1,llm) ! temp pot reales
73	REAL masserea2(ip1jmp1,llm) ! masse
74	REAL psrea2(ip1jmp1) ! ps
75	real latmin
76
77	real alpha_q(ip1jmp1)
78	real alpha_T(ip1jmp1),alpha_P(ip1jmp1)
79	real alpha_u(ip1jmp1),alpha_v(ip1jm)
80	real dday_step,toto,reste,itau_test
81	INTEGER step_rea,count_no_rea
82
83	integer ilon,ilat
84	real factt,ztau(ip1jmp1)
85
86	INTEGER itau,ij,l,i,j
87	integer ncidpl,varidpl,nlev,status
88	integer rcod,rid
89	real ditau,tau,a
90	save nlev
91
92	c TEST SUR QSAT
93	real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)
94	real pkf(ip1jmp1,llm)
95	real pres(ip1jmp1,llm)
96	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
97
98	real qsat(ip1jmp1,llm)
99	real unskap
100	real tnat(ip1jmp1,llm)
101	ccccccccccccccccc
102
103
104	LOGICAL first
105	save first
106	data first/.true./
107
108	save ucovrea1,vcovrea1,tetarea1,masserea1,psrea1,qrea1
109	save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2
110
111	save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
112	save step_rea,count_no_rea
113
114	character*10 file
115	integer igrads
116	real dtgrads
117	save igrads,dtgrads
118	data igrads,dtgrads/2,100./
119
120	C-----------------------------------------------------------------------
121	c calcul de l'humidite saturante
122	C-----------------------------------------------------------------------
123	print*,'OK0'
124	CALL pression( ip1jmp1, ap, bp, ps, p )
125	call massdair(p,masse)
126	print*,'OK1'
127	CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
128	print*,'OK2'
129	tnat(:,:)=pk(:,:)*teta(:,:)/cpp
130	print*,'OK3'
131	unskap = 1./ kappa
132	pres(:,:)=preff(pk(:,:)/cpp)*unskap
133	print*,'OK4'
134	call q_sat(iip1jjp1llm,tnat,pres,qsat)
135
136	C-----------------------------------------------------------------------
137
138	c-----------------------------------------------------------------------
139	c initialisations pour la lecture des reanalyses.
140	c alpha determine la part des injections de donnees a chaque etape
141	c alpha=1 signifie pas d'injection
142	c alpha=0 signifie injection totale
143	c-----------------------------------------------------------------------
144
145	print*,'ONLINE=',online
146	if(online.eq.-1) then
147	return
148	endif
149
150	if (first) then
151
152	print*,'initialisation du guide '
153	call conf_guide
154	print*,'apres conf_guide'
155
156	file='guide'
157	call inigrads(igrads,iip1
158	s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
159	s ,llm,presnivs,1.
160	s ,dtgrads,file,'dyn_zon ')
161
162	print*
163	s ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
164
165	if(online.eq.-1) return
166	if (online.eq.1) then
167
168	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
169	c Constantes de temps de rappel en jour
170	c 0.1 c'est en gros 2h30.
171	c 1e10 est une constante infinie donc en gros pas de guidage
172	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
173	c coordonnees du centre du zoom
174	call coordij(clon,clat,ilon,ilat)
175	c aire de la maille au centre du zoom
176	aire_min=aire(ilon+(ilat-1)*iip1)
177	c aire maximale de la maille
178	aire_max=0.
179	do ij=1,ip1jmp1
180	aire_max=max(aire_max,aire(ij))
181	enddo
182	C factt = pas de temps en fraction de jour
183	factt=dtvr*iperiod/daysec
184
185	c subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)
186	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
187	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
188	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
189	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
190	call tau2alpha(1,iip1,jjp1,factt,tau_min_q,tau_max_q,alpha_q)
191
192	call dump2d(iip1,jjp1,aire,'AIRE MAILLe ')
193	call dump2d(iip1,jjp1,alpha_u,'COEFF U ')
194	call dump2d(iip1,jjp1,alpha_T,'COEFF T ')
195
196	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
197	c Cas ou on force exactement par les variables analysees
198	else
199	alpha_T=0.
200	alpha_u=0.
201	alpha_v=0.
202	alpha_P=0.
203	c physic=.false.
204	endif
205
206	itau_test=1001
207	step_rea=1
208	count_no_rea=0
209	ncidpl=-99
210
211	c itau_test montre si l'importation a deja ete faite au rang itau
212	c lecture d'un fichier netcdf pour determiner le nombre de niveaux
213	if (guide_u) then
214	if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
215	endif
216	c
217	if (guide_v) then
218	if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
219	endif
220	c
221	if (guide_T) then
222	if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
223	endif
224	c
225	if (guide_Q) then
226	if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
227	endif
228	c
229	if (ncep) then
230	status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
231	else
232	status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
233	endif
234	status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
235	print *,'nlev', nlev
236	call ncclos(ncidpl,rcod)
237	c Lecture du premier etat des reanalyses.
238	call read_reanalyse(1,ps
239	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
240	qrea2(:,:)=max(qrea2(:,:),0.1)
241
242
243	c-----------------------------------------------------------------------
244	c Debut de l'integration temporelle:
245	c ----------------------------------
246
247	endif ! first
248	c
249	C-----------------------------------------------------------------------
250	C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
251	C-----------------------------------------------------------------------
252
253	ditau=real(itau)
254	DDAY_step=real(day_step)
255	write(,)'ditau,dday_step'
256	write(,)ditau,dday_step
257	toto=4*ditau/dday_step
258	reste=toto-aint(toto)
259	c write(,)'toto,reste',toto,reste
260
261	if (reste.eq.0.) then
262	if (itau_test.eq.itau) then
263	write(,)'deuxieme passage de advreel a itau=',itau
264	stop
265	else
266	vcovrea1(:,:)=vcovrea2(:,:)
267	ucovrea1(:,:)=ucovrea2(:,:)
268	tetarea1(:,:)=tetarea2(:,:)
269	qrea1(:,:)=qrea2(:,:)
270
271	print*,'LECTURE REANALYSES, pas ',step_rea
272	s ,'apres ',count_no_rea,' non lectures'
273	step_rea=step_rea+1
274	itau_test=itau
275	call read_reanalyse(step_rea,ps
276	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
277	qrea2(:,:)=max(qrea2(:,:),0.1)
278	factt=dtvr*iperiod/daysec
279	ztau(:)=factt/max(alpha_T(:),1.e-10)
280	call wrgrads(igrads,1,aire ,'aire ','aire ' )
281	call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
282	call wrgrads(igrads,1,alpha_u,'au ','au ' )
283	call wrgrads(igrads,1,alpha_T,'at ','at ' )
284	call wrgrads(igrads,1,ztau,'taut ','taut ' )
285	call wrgrads(igrads,llm,ucov,'u ','u ' )
286	call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
287	call wrgrads(igrads,llm,teta,'T ','T ' )
288	call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
289	call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
290	call wrgrads(igrads,llm,q,'Q ','Q ' )
291
292	call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
293
294	endif
295	else
296	count_no_rea=count_no_rea+1
297	endif
298
299	C-----------------------------------------------------------------------
300	c Guidage
301	c x_gcm = a * x_gcm + (1-a) * x_reanalyses
302	C-----------------------------------------------------------------------
303
304	if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
305
306	ditau=real(itau)
307	dday_step=real(day_step)
308
309
310	tau=4*ditau/dday_step
311	tau=tau-aint(tau)
312
313	c ucov
314	if (guide_u) then
315	do l=1,llm
316	do ij=1,ip1jmp1
317	a=(1.-tau)ucovrea1(ij,l)+tauucovrea2(ij,l)
318	ucov(ij,l)=(1.-alpha_u(ij))ucov(ij,l)+alpha_u(ij)a
319	if (first.and.ini_anal) ucov(ij,l)=a
320	enddo
321	enddo
322	endif
323
324	c teta
325	if (guide_T) then
326	do l=1,llm
327	do ij=1,ip1jmp1
328	a=(1.-tau)tetarea1(ij,l)+tautetarea2(ij,l)
329	teta(ij,l)=(1.-alpha_T(ij))teta(ij,l)+alpha_T(ij)a
330	if (first.and.ini_anal) teta(ij,l)=a
331	enddo
332	enddo
333	endif
334
335	c P
336	if (guide_P) then
337	do ij=1,ip1jmp1
338	a=(1.-tau)psrea1(ij)+taupsrea2(ij)
339	ps(ij)=(1.-alpha_P(ij))ps(ij)+alpha_P(ij)a
340	if (first.and.ini_anal) ps(ij)=a
341	enddo
342	CALL pression(ip1jmp1,ap,bp,ps,p)
343	CALL massdair(p,masse)
344	endif
345
346
347	c q
348	if (guide_Q) then
349	do l=1,llm
350	do ij=1,ip1jmp1
351	a=(1.-tau)qrea1(ij,l)+tauqrea2(ij,l)
352	c hum relative en % -> hum specif
353	a=qsat(ij,l)a0.01
354	q(ij,l)=(1.-alpha_Q(ij))q(ij,l)+alpha_Q(ij)a
355	if (first.and.ini_anal) q(ij,l)=a
356	enddo
357	enddo
358	endif
359
360	c vcov
361	if (guide_v) then
362	do l=1,llm
363	do ij=1,ip1jm
364	a=(1.-tau)vcovrea1(ij,l)+tauvcovrea2(ij,l)
365	vcov(ij,l)=(1.-alpha_v(ij))vcov(ij,l)+alpha_v(ij)a
366	if (first.and.ini_anal) vcov(ij,l)=a
367	enddo
368	if (first.and.ini_anal) vcov(ij,l)=a
369	enddo
370	endif
371
372	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
373	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
374	c call dump2d(iip1,jjp1,teta,'TETA ')
375
376	first=.false.
377
378	return
379	end
380
381	c=======================================================================
382	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
383	c=======================================================================
384
385	implicit none
386
387	#include "dimensions.h"
388	#include "paramet.h"
389	#include "comconst.h"
390	#include "comgeom2.h"
391	#include "guide.h"
392	#include "serre.h"
393
394	c arguments :
395	integer type
396	integer pim,pjm
397	real factt,taumin,taumax,dxdymin,dxdymax
398	real dxdy_,alpha(pim,pjm)
399	real dxdy_min,dxdy_max
400
401	c local :
402	real alphamin,alphamax,gamma,xi
403	save gamma
404	integer i,j,ilon,ilat
405
406	logical first
407	save first
408	data first/.true./
409
410	real cus(iip1,jjp1),cvs(iip1,jjp1)
411	real cuv(iip1,jjm),cvu(iip1,jjp1)
412	real zdx(iip1,jjp1),zdy(iip1,jjp1)
413
414	real zlat
415	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
416	common/comdxdy/dxdys,dxdyu,dxdyv
417
418	if (first) then
419	do j=2,jjm
420	do i=2,iip1
421	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
422	enddo
423	zdx(1,j)=zdx(iip1,j)
424	enddo
425	do j=2,jjm
426	do i=1,iip1
427	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
428	enddo
429	enddo
430	do i=1,iip1
431	zdx(i,1)=zdx(i,2)
432	zdx(i,jjp1)=zdx(i,jjm)
433	zdy(i,1)=zdy(i,2)
434	zdy(i,jjp1)=zdy(i,jjm)
435	enddo
436	do j=1,jjp1
437	do i=1,iip1
438	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
439	enddo
440	enddo
441	do j=1,jjp1
442	do i=1,iim
443	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
444	enddo
445	dxdyu(iip1,j)=dxdyu(1,j)
446	enddo
447	do j=1,jjm
448	do i=1,iip1
449	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
450	enddo
451	enddo
452
453	call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
454	call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
455	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
456
457	c coordonnees du centre du zoom
458	call coordij(clon,clat,ilon,ilat)
459	c aire de la maille au centre du zoom
460	dxdy_min=dxdys(ilon,ilat)
461	c dxdy maximale de la maille
462	dxdy_max=0.
463	do j=1,jjp1
464	do i=1,iip1
465	dxdy_max=max(dxdy_max,dxdys(i,j))
466	enddo
467	enddo
468
469	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
470	print*,'ATTENTION modele peu zoome'
471	print*,'ATTENTION on prend une constante de guidage cste'
472	gamma=0.
473	else
474	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
475	print*,'gamma=',gamma
476	if (gamma.lt.1.e-5) then
477	print*,'gamma =',gamma,'<1e-5'
478	stop
479	endif
480	print*,'gamma=',gamma
481	gamma=log(0.5)/log(gamma)
482	endif
483	endif
484
485	alphamin=factt/taumax
486	alphamax=factt/taumin
487
488	do j=1,pjm
489	do i=1,pim
490	if (type.eq.1) then
491	dxdy_=dxdys(i,j)
492	zlat=rlatu(j)*180./pi
493	elseif (type.eq.2) then
494	dxdy_=dxdyu(i,j)
495	zlat=rlatu(j)*180./pi
496	elseif (type.eq.3) then
497	dxdy_=dxdyv(i,j)
498	zlat=rlatv(j)*180./pi
499	endif
500	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
501	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
502	alpha(i,j)=alphamin
503	else
504	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
505	xi=min(xi,1.)
506	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
507	alpha(i,j)=xialphamin+(1.-xi)alphamax
508	else
509	alpha(i,j)=0.
510	endif
511	endif
512	enddo
513	enddo
514
515
516	return
517	end

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

Download in other formats:

Original Format