Context Navigation

guide.F @ 644

Last change on this file since 644 was 644, checked in by Laurent Fairhead, 19 years ago
Synchronisation avec tous les diagnostiques de Ionela IM Inclusion du slab ocean IM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.6 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	cIM 180305 real aire_min,aire_max
84	integer ilon,ilat
85	real factt,ztau(ip1jmp1)
86
87	INTEGER itau,ij,l,i,j
88	integer ncidpl,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	ncidpl=-99
211
212	c itau_test montre si l'importation a deja ete faite au rang itau
213	c lecture d'un fichier netcdf pour determiner le nombre de niveaux
214	if (guide_u) then
215	if (ncidpl.eq.-99) ncidpl=NCOPN('u.nc',NCNOWRIT,rcod)
216	endif
217	c
218	if (guide_v) then
219	if (ncidpl.eq.-99) ncidpl=NCOPN('v.nc',NCNOWRIT,rcod)
220	endif
221	c
222	if (guide_T) then
223	if (ncidpl.eq.-99) ncidpl=NCOPN('T.nc',NCNOWRIT,rcod)
224	endif
225	c
226	if (guide_Q) then
227	if (ncidpl.eq.-99) ncidpl=NCOPN('hur.nc',NCNOWRIT,rcod)
228	endif
229	c
230	if (ncep) then
231	status=NF_INQ_DIMID(ncidpl,'LEVEL',rid)
232	else
233	status=NF_INQ_DIMID(ncidpl,'PRESSURE',rid)
234	endif
235	status=NF_INQ_DIMLEN(ncidpl,rid,nlev)
236	print *,'nlev', nlev
237	call ncclos(ncidpl,rcod)
238	c Lecture du premier etat des reanalyses.
239	call read_reanalyse(1,ps
240	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
241	qrea2(:,:)=max(qrea2(:,:),0.1)
242
243
244	c-----------------------------------------------------------------------
245	c Debut de l'integration temporelle:
246	c ----------------------------------
247
248	endif ! first
249	c
250	C-----------------------------------------------------------------------
251	C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
252	C-----------------------------------------------------------------------
253
254	ditau=real(itau)
255	DDAY_step=real(day_step)
256	write(,)'ditau,dday_step'
257	write(,)ditau,dday_step
258	toto=4*ditau/dday_step
259	reste=toto-aint(toto)
260	c write(,)'toto,reste',toto,reste
261
262	if (reste.eq.0.) then
263	if (itau_test.eq.itau) then
264	write(,)'deuxieme passage de advreel a itau=',itau
265	stop
266	else
267	vcovrea1(:,:)=vcovrea2(:,:)
268	ucovrea1(:,:)=ucovrea2(:,:)
269	tetarea1(:,:)=tetarea2(:,:)
270	qrea1(:,:)=qrea2(:,:)
271
272	print*,'LECTURE REANALYSES, pas ',step_rea
273	s ,'apres ',count_no_rea,' non lectures'
274	step_rea=step_rea+1
275	itau_test=itau
276	call read_reanalyse(step_rea,ps
277	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
278	qrea2(:,:)=max(qrea2(:,:),0.1)
279	factt=dtvr*iperiod/daysec
280	ztau(:)=factt/max(alpha_T(:),1.e-10)
281	call wrgrads(igrads,1,aire ,'aire ','aire ' )
282	call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
283	call wrgrads(igrads,1,alpha_u,'au ','au ' )
284	call wrgrads(igrads,1,alpha_T,'at ','at ' )
285	call wrgrads(igrads,1,ztau,'taut ','taut ' )
286	call wrgrads(igrads,llm,ucov,'u ','u ' )
287	call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
288	call wrgrads(igrads,llm,teta,'T ','T ' )
289	call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
290	call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
291	call wrgrads(igrads,llm,q,'Q ','Q ' )
292
293	call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
294
295	endif
296	else
297	count_no_rea=count_no_rea+1
298	endif
299
300	C-----------------------------------------------------------------------
301	c Guidage
302	c x_gcm = a * x_gcm + (1-a) * x_reanalyses
303	C-----------------------------------------------------------------------
304
305	if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
306
307	ditau=real(itau)
308	dday_step=real(day_step)
309
310
311	tau=4*ditau/dday_step
312	tau=tau-aint(tau)
313
314	c ucov
315	if (guide_u) then
316	do l=1,llm
317	do ij=1,ip1jmp1
318	a=(1.-tau)ucovrea1(ij,l)+tauucovrea2(ij,l)
319	ucov(ij,l)=(1.-alpha_u(ij))ucov(ij,l)+alpha_u(ij)a
320	if (first.and.ini_anal) ucov(ij,l)=a
321	enddo
322	enddo
323	endif
324
325	c teta
326	if (guide_T) then
327	do l=1,llm
328	do ij=1,ip1jmp1
329	a=(1.-tau)tetarea1(ij,l)+tautetarea2(ij,l)
330	teta(ij,l)=(1.-alpha_T(ij))teta(ij,l)+alpha_T(ij)a
331	if (first.and.ini_anal) teta(ij,l)=a
332	enddo
333	enddo
334	endif
335
336	c P
337	if (guide_P) then
338	do ij=1,ip1jmp1
339	a=(1.-tau)psrea1(ij)+taupsrea2(ij)
340	ps(ij)=(1.-alpha_P(ij))ps(ij)+alpha_P(ij)a
341	if (first.and.ini_anal) ps(ij)=a
342	enddo
343	CALL pression(ip1jmp1,ap,bp,ps,p)
344	CALL massdair(p,masse)
345	endif
346
347
348	c q
349	if (guide_Q) then
350	do l=1,llm
351	do ij=1,ip1jmp1
352	a=(1.-tau)qrea1(ij,l)+tauqrea2(ij,l)
353	c hum relative en % -> hum specif
354	a=qsat(ij,l)a0.01
355	q(ij,l)=(1.-alpha_Q(ij))q(ij,l)+alpha_Q(ij)a
356	if (first.and.ini_anal) q(ij,l)=a
357	enddo
358	enddo
359	endif
360
361	c vcov
362	if (guide_v) then
363	do l=1,llm
364	do ij=1,ip1jm
365	a=(1.-tau)vcovrea1(ij,l)+tauvcovrea2(ij,l)
366	vcov(ij,l)=(1.-alpha_v(ij))vcov(ij,l)+alpha_v(ij)a
367	if (first.and.ini_anal) vcov(ij,l)=a
368	enddo
369	if (first.and.ini_anal) vcov(ij,l)=a
370	enddo
371	endif
372
373	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
374	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
375	c call dump2d(iip1,jjp1,teta,'TETA ')
376
377	first=.false.
378
379	return
380	end
381
382	c=======================================================================
383	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
384	c=======================================================================
385
386	implicit none
387
388	#include "dimensions.h"
389	#include "paramet.h"
390	#include "comconst.h"
391	#include "comgeom2.h"
392	#include "guide.h"
393	#include "serre.h"
394
395	c arguments :
396	integer type
397	integer pim,pjm
398	real factt,taumin,taumax,dxdymin,dxdymax
399	real dxdy_,alpha(pim,pjm)
400	real dxdy_min,dxdy_max
401
402	c local :
403	real alphamin,alphamax,gamma,xi
404	save gamma
405	integer i,j,ilon,ilat
406
407	logical first
408	save first
409	data first/.true./
410
411	real cus(iip1,jjp1),cvs(iip1,jjp1)
412	real cuv(iip1,jjm),cvu(iip1,jjp1)
413	real zdx(iip1,jjp1),zdy(iip1,jjp1)
414
415	real zlat
416	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
417	common/comdxdy/dxdys,dxdyu,dxdyv
418
419	if (first) then
420	do j=2,jjm
421	do i=2,iip1
422	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
423	enddo
424	zdx(1,j)=zdx(iip1,j)
425	enddo
426	do j=2,jjm
427	do i=1,iip1
428	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
429	enddo
430	enddo
431	do i=1,iip1
432	zdx(i,1)=zdx(i,2)
433	zdx(i,jjp1)=zdx(i,jjm)
434	zdy(i,1)=zdy(i,2)
435	zdy(i,jjp1)=zdy(i,jjm)
436	enddo
437	do j=1,jjp1
438	do i=1,iip1
439	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
440	enddo
441	enddo
442	do j=1,jjp1
443	do i=1,iim
444	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
445	enddo
446	dxdyu(iip1,j)=dxdyu(1,j)
447	enddo
448	do j=1,jjm
449	do i=1,iip1
450	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
451	enddo
452	enddo
453
454	call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
455	call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
456	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
457
458	c coordonnees du centre du zoom
459	call coordij(clon,clat,ilon,ilat)
460	c aire de la maille au centre du zoom
461	dxdy_min=dxdys(ilon,ilat)
462	c dxdy maximale de la maille
463	dxdy_max=0.
464	do j=1,jjp1
465	do i=1,iip1
466	dxdy_max=max(dxdy_max,dxdys(i,j))
467	enddo
468	enddo
469
470	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
471	print*,'ATTENTION modele peu zoome'
472	print*,'ATTENTION on prend une constante de guidage cste'
473	gamma=0.
474	else
475	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
476	print*,'gamma=',gamma
477	if (gamma.lt.1.e-5) then
478	print*,'gamma =',gamma,'<1e-5'
479	stop
480	endif
481	print*,'gamma=',gamma
482	gamma=log(0.5)/log(gamma)
483	endif
484	endif
485
486	alphamin=factt/taumax
487	alphamax=factt/taumin
488
489	do j=1,pjm
490	do i=1,pim
491	if (type.eq.1) then
492	dxdy_=dxdys(i,j)
493	zlat=rlatu(j)*180./pi
494	elseif (type.eq.2) then
495	dxdy_=dxdyu(i,j)
496	zlat=rlatu(j)*180./pi
497	elseif (type.eq.3) then
498	dxdy_=dxdyv(i,j)
499	zlat=rlatv(j)*180./pi
500	endif
501	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
502	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
503	alpha(i,j)=alphamin
504	else
505	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
506	xi=min(xi,1.)
507	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
508	alpha(i,j)=xialphamin+(1.-xi)alphamax
509	else
510	alpha(i,j)=0.
511	endif
512	endif
513	enddo
514	enddo
515
516
517	return
518	end

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

Download in other formats:

Original Format