Context Navigation

guide.F @ 1046

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

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

Download in other formats:

Original Format