Context Navigation

guide_p.F @ 865

Last change on this file since 865 was 865, checked in by lsce, 17 years ago
ACo : correction bug qui entrainait un depassement de tableau - deja corrige dans dyn3d
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.7 KB

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

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

Download in other formats:

Original Format