Context Navigation

guide.F @ 2629

Last change on this file since 2629 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.1 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	real aire_min,aire_max
84	integer ilon,ilat
85	real factt,ztau(ip1jmp1)
86
87	INTEGER itau,ij,l,i,j
88	integer ncidt,varidpl,nlev,status
89	integer rcod,rid
90	real ditau,tau,a
91	save nlev
92
93	c TEST SUR QSAT
94	real p(ip1jmp1,llmp1),pk(ip1jmp1,llm),pks(ip1jmp1)
95	real pkf(ip1jmp1,llm)
96	real pres(ip1jmp1,llm)
97	REAL alpha(ip1jmp1,llm),beta(ip1jmp1,llm)
98
99	real qsat(ip1jmp1,llm)
100	real unskap
101	real tnat(ip1jmp1,llm)
102	ccccccccccccccccc
103
104
105	LOGICAL first
106	save first
107	data first/.true./
108
109	save ucovrea1,vcovrea1,tetarea1,masserea1,psrea1,qrea1
110	save ucovrea2,vcovrea2,tetarea2,masserea2,psrea2,qrea2
111
112	save alpha_T,alpha_q,alpha_u,alpha_v,alpha_P,itau_test
113	save step_rea,count_no_rea
114
115	character*10 file
116	integer igrads
117	real dtgrads
118	save igrads,dtgrads
119	data igrads,dtgrads/2,100./
120
121	C-----------------------------------------------------------------------
122	c calcul de l'humidite saturante
123	C-----------------------------------------------------------------------
124	print*,'OK0'
125	CALL pression( ip1jmp1, ap, bp, ps, p )
126	call massdair(p,masse)
127	print*,'OK1'
128	CALL exner_hyb(ip1jmp1,ps,p,alpha,beta,pks,pk,pkf)
129	print*,'OK2'
130	tnat(:,:)=pk(:,:)*teta(:,:)/cpp
131	print*,'OK3'
132	unskap = 1./ kappa
133	pres(:,:)=preff(pk(:,:)/cpp)*unskap
134	print*,'OK4'
135	call q_sat(iip1jjp1llm,tnat,pres,qsat)
136
137	C-----------------------------------------------------------------------
138
139	c-----------------------------------------------------------------------
140	c initialisations pour la lecture des reanalyses.
141	c alpha determine la part des injections de donnees a chaque etape
142	c alpha=1 signifie pas d'injection
143	c alpha=0 signifie injection totale
144	c-----------------------------------------------------------------------
145
146	print*,'ONLINE=',online
147	if(online.eq.-1) then
148	return
149	endif
150
151	if (first) then
152
153	print*,'initialisation du guide '
154	call conf_guide
155	print*,'apres conf_guide'
156
157	file='guide'
158	call inigrads(igrads,iip1
159	s ,rlonv,180./pi,-180.,180.,jjp1,rlatu,-90.,90.,180./pi
160	s ,llm,presnivs,1.
161	s ,dtgrads,file,'dyn_zon ')
162
163	print*
164	s ,'1: en-ligne, 0: hors-ligne (x=x_rea), -1: climat (x=x_gcm)'
165
166	if(online.eq.-1) return
167	if (online.eq.1) then
168
169	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
170	c Constantes de temps de rappel en jour
171	c 0.1 c'est en gros 2h30.
172	c 1e10 est une constante infinie donc en gros pas de guidage
173	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
174	c coordonnees du centre du zoom
175	call coordij(clon,clat,ilon,ilat)
176	c aire de la maille au centre du zoom
177	aire_min=aire(ilon+(ilat-1)*iip1)
178	c aire maximale de la maille
179	aire_max=0.
180	do ij=1,ip1jmp1
181	aire_max=max(aire_max,aire(ij))
182	enddo
183	C factt = pas de temps en fraction de jour
184	factt=dtvr*iperiod/daysec
185
186	c subroutine tau2alpha(type,im,jm,factt,taumin,taumax,alpha)
187	call tau2alpha(3,iip1,jjm ,factt,tau_min_v,tau_max_v,alpha_v)
188	call tau2alpha(2,iip1,jjp1,factt,tau_min_u,tau_max_u,alpha_u)
189	call tau2alpha(1,iip1,jjp1,factt,tau_min_T,tau_max_T,alpha_T)
190	call tau2alpha(1,iip1,jjp1,factt,tau_min_P,tau_max_P,alpha_P)
191	call tau2alpha(1,iip1,jjp1,factt,tau_min_q,tau_max_q,alpha_q)
192
193	call dump2d(iip1,jjp1,aire,'AIRE MAILLe ')
194	call dump2d(iip1,jjp1,alpha_u,'COEFF U ')
195	call dump2d(iip1,jjp1,alpha_T,'COEFF T ')
196
197	cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
198	c Cas ou on force exactement par les variables analysees
199	else
200	alpha_T=0.
201	alpha_u=0.
202	alpha_v=0.
203	alpha_P=0.
204	c physic=.false.
205	endif
206
207	itau_test=1001
208	step_rea=1
209	count_no_rea=0
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	ncidt=NCOPN('T.nc',NCNOWRIT,rcod)
214	if (ncep) then
215	status=NF_INQ_DIMID(ncidt,'LEVEL',rid)
216	else
217	status=NF_INQ_DIMID(ncidt,'PRESSURE',rid)
218	endif
219	status=NF_INQ_DIMLEN(ncidt,rid,nlev)
220	print *,'nlev', nlev
221	call ncclos(ncidt,rcod)
222	c Lecture du premier etat des reanalyses.
223	call read_reanalyse(1,ps
224	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
225	qrea2(:,:)=max(qrea2(:,:),0.1)
226
227
228	c-----------------------------------------------------------------------
229	c Debut de l'integration temporelle:
230	c ----------------------------------
231
232	endif ! first
233	c
234	C-----------------------------------------------------------------------
235	C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
236	C-----------------------------------------------------------------------
237
238	ditau=real(itau)
239	DDAY_step=real(day_step)
240	write(,)'ditau,dday_step'
241	write(,)ditau,dday_step
242	toto=4*ditau/dday_step
243	reste=toto-aint(toto)
244	c write(,)'toto,reste',toto,reste
245
246	if (reste.eq.0.) then
247	if (itau_test.eq.itau) then
248	write(,)'deuxieme passage de advreel a itau=',itau
249	stop
250	else
251	vcovrea1(:,:)=vcovrea2(:,:)
252	ucovrea1(:,:)=ucovrea2(:,:)
253	tetarea1(:,:)=tetarea2(:,:)
254	qrea1(:,:)=qrea2(:,:)
255
256	print*,'LECTURE REANALYSES, pas ',step_rea
257	s ,'apres ',count_no_rea,' non lectures'
258	step_rea=step_rea+1
259	itau_test=itau
260	call read_reanalyse(step_rea,ps
261	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
262	qrea2(:,:)=max(qrea2(:,:),0.1)
263	factt=dtvr*iperiod/daysec
264	ztau(:)=factt/max(alpha_T(:),1.e-10)
265	call wrgrads(igrads,1,aire ,'aire ','aire ' )
266	call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
267	call wrgrads(igrads,1,alpha_u,'au ','au ' )
268	call wrgrads(igrads,1,alpha_T,'at ','at ' )
269	call wrgrads(igrads,1,ztau,'taut ','taut ' )
270	call wrgrads(igrads,llm,ucov,'u ','u ' )
271	call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
272	call wrgrads(igrads,llm,teta,'T ','T ' )
273	call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
274	call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
275	call wrgrads(igrads,llm,q,'Q ','Q ' )
276
277	call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
278
279	endif
280	else
281	count_no_rea=count_no_rea+1
282	endif
283
284	C-----------------------------------------------------------------------
285	c Guidage
286	c x_gcm = a * x_gcm + (1-a) * x_reanalyses
287	C-----------------------------------------------------------------------
288
289	if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
290
291	ditau=real(itau)
292	dday_step=real(day_step)
293
294
295	tau=4*ditau/dday_step
296	tau=tau-aint(tau)
297
298	print*,'ATTENTION !!!! ON NE GUIDE QUE JUSQU A 15N'
299
300	c ucov
301	if (guide_u) then
302	do l=1,llm
303	do ij=1,ip1jmp1
304	a=(1.-tau)ucovrea1(ij,l)+tauucovrea2(ij,l)
305	ucov(ij,l)=(1.-alpha_u(ij))ucov(ij,l)+alpha_u(ij)a
306	if (first.and.ini_anal) ucov(ij,l)=a
307	enddo
308	enddo
309	endif
310
311	c teta
312	if (guide_T) then
313	do l=1,llm
314	do ij=1,ip1jmp1
315	a=(1.-tau)tetarea1(ij,l)+tautetarea2(ij,l)
316	teta(ij,l)=(1.-alpha_T(ij))teta(ij,l)+alpha_T(ij)a
317	if (first.and.ini_anal) teta(ij,l)=a
318	enddo
319	enddo
320	endif
321
322	c P
323	if (guide_P) then
324	do ij=1,ip1jmp1
325	a=(1.-tau)psrea1(ij)+taupsrea2(ij)
326	ps(ij)=(1.-alpha_P(ij))ps(ij)+alpha_P(ij)a
327	if (first.and.ini_anal) ps(ij)=a
328	enddo
329	CALL pression(ip1jmp1,ap,bp,ps,p)
330	CALL massdair(p,masse)
331	endif
332
333
334	c q
335	if (guide_Q) then
336	do l=1,llm
337	do ij=1,ip1jmp1
338	a=(1.-tau)qrea1(ij,l)+tauqrea2(ij,l)
339	c hum relative en % -> hum specif
340	a=qsat(ij,l)a0.01
341	q(ij,l)=(1.-alpha_Q(ij))q(ij,l)+alpha_Q(ij)a
342	if (first.and.ini_anal) q(ij,l)=a
343	enddo
344	enddo
345	endif
346
347	c vcov
348	if (guide_v) then
349	do l=1,llm
350	do ij=1,ip1jm
351	a=(1.-tau)vcovrea1(ij,l)+tauvcovrea2(ij,l)
352	vcov(ij,l)=(1.-alpha_v(ij))vcov(ij,l)+alpha_v(ij)a
353	if (first.and.ini_anal) vcov(ij,l)=a
354	enddo
355	if (first.and.ini_anal) vcov(ij,l)=a
356	enddo
357	endif
358
359	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
360	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
361	c call dump2d(iip1,jjp1,teta,'TETA ')
362
363	first=.false.
364
365	return
366	end
367
368	c=======================================================================
369	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
370	c=======================================================================
371
372	implicit none
373
374	#include "dimensions.h"
375	#include "paramet.h"
376	#include "comconst.h"
377	#include "comgeom2.h"
378	#include "guide.h"
379	#include "serre.h"
380
381	c arguments :
382	integer type
383	integer pim,pjm
384	real factt,taumin,taumax,dxdymin,dxdymax
385	real dxdy_,alpha(pim,pjm)
386	real dxdy_min,dxdy_max
387
388	c local :
389	real alphamin,alphamax,gamma,xi
390	save gamma
391	integer i,j,ilon,ilat
392
393	logical first
394	save first
395	data first/.true./
396
397	real cus(iip1,jjp1),cvs(iip1,jjp1)
398	real cuv(iip1,jjm),cvu(iip1,jjp1)
399	real zdx(iip1,jjp1),zdy(iip1,jjp1)
400
401	real zlat
402	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
403	common/comdxdy/dxdys,dxdyu,dxdyv
404
405	if (first) then
406	do j=2,jjm
407	do i=2,iip1
408	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
409	enddo
410	zdx(1,j)=zdx(iip1,j)
411	enddo
412	do j=2,jjm
413	do i=1,iip1
414	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
415	enddo
416	enddo
417	do i=1,iip1
418	zdx(i,1)=zdx(i,2)
419	zdx(i,jjp1)=zdx(i,jjm)
420	zdy(i,1)=zdy(i,2)
421	zdy(i,jjp1)=zdy(i,jjm)
422	enddo
423	do j=1,jjp1
424	do i=1,iip1
425	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
426	enddo
427	enddo
428	do j=1,jjp1
429	do i=1,iim
430	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
431	enddo
432	dxdyu(iip1,j)=dxdyu(1,j)
433	enddo
434	do j=1,jjm
435	do i=1,iip1
436	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
437	enddo
438	enddo
439
440	call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
441	call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
442	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
443
444	c coordonnees du centre du zoom
445	call coordij(clon,clat,ilon,ilat)
446	c aire de la maille au centre du zoom
447	dxdy_min=dxdys(ilon,ilat)
448	c dxdy maximale de la maille
449	dxdy_max=0.
450	do j=1,jjp1
451	do i=1,iip1
452	dxdy_max=max(dxdy_max,dxdys(i,j))
453	enddo
454	enddo
455
456	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
457	print*,'ATTENTION modele peu zoome'
458	print*,'ATTENTION on prend une constante de guidage cste'
459	gamma=0.
460	else
461	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
462	print*,'gamma=',gamma
463	if (gamma.lt.1.e-5) then
464	print*,'gamma =',gamma,'<1e-5'
465	stop
466	endif
467	print*,'gamma=',gamma
468	gamma=log(0.5)/log(gamma)
469	endif
470	endif
471
472	alphamin=factt/taumax
473	alphamax=factt/taumin
474
475	do j=1,pjm
476	do i=1,pim
477	if (type.eq.1) then
478	dxdy_=dxdys(i,j)
479	zlat=rlatu(j)*180./pi
480	elseif (type.eq.2) then
481	dxdy_=dxdyu(i,j)
482	zlat=rlatu(j)*180./pi
483	elseif (type.eq.3) then
484	dxdy_=dxdyv(i,j)
485	zlat=rlatv(j)*180./pi
486	endif
487	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
488	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
489	xi=min(xi,1.)
490	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
491	alpha(i,j)=xialphamin+(1.-xi)alphamax
492	else
493	alpha(i,j)=0.
494	endif
495	enddo
496	enddo
497
498
499	return
500	end

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

Download in other formats:

Original Format