Context Navigation

guide_p.F @ 1034

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

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

Download in other formats:

Original Format