Context Navigation

guide_p.F @ 814

Last change on this file since 814 was 764, checked in by Laurent Fairhead, 17 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.8 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	if (first.and.ini_anal) vcov(ij,l)=a
418	enddo
419	endif
420
421	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
422	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
423	c call dump2d(iip1,jjp1,teta,'TETA ')
424
425	first=.false.
426
427	return
428	end
429
430	c=======================================================================
431	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
432	c=======================================================================
433
434	implicit none
435
436	#include "dimensions.h"
437	#include "paramet.h"
438	#include "comconst.h"
439	#include "comgeom2.h"
440	#include "guide.h"
441	#include "serre.h"
442
443	c arguments :
444	integer type
445	integer pim,pjm
446	real factt,taumin,taumax,dxdymin,dxdymax
447	real dxdy_,alpha(pim,pjm)
448	real dxdy_min,dxdy_max
449
450	c local :
451	real alphamin,alphamax,gamma,xi
452	save gamma
453	integer i,j,ilon,ilat
454
455	logical first
456	save first
457	data first/.true./
458
459	real cus(iip1,jjp1),cvs(iip1,jjp1)
460	real cuv(iip1,jjm),cvu(iip1,jjp1)
461	real zdx(iip1,jjp1),zdy(iip1,jjp1)
462
463	real zlat
464	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
465	common/comdxdy/dxdys,dxdyu,dxdyv
466
467	if (first) then
468	do j=2,jjm
469	do i=2,iip1
470	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
471	enddo
472	zdx(1,j)=zdx(iip1,j)
473	enddo
474	do j=2,jjm
475	do i=1,iip1
476	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
477	enddo
478	enddo
479	do i=1,iip1
480	zdx(i,1)=zdx(i,2)
481	zdx(i,jjp1)=zdx(i,jjm)
482	zdy(i,1)=zdy(i,2)
483	zdy(i,jjp1)=zdy(i,jjm)
484	enddo
485	do j=1,jjp1
486	do i=1,iip1
487	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
488	enddo
489	enddo
490	do j=1,jjp1
491	do i=1,iim
492	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
493	enddo
494	dxdyu(iip1,j)=dxdyu(1,j)
495	enddo
496	do j=1,jjm
497	do i=1,iip1
498	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
499	enddo
500	enddo
501
502	call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
503	call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
504	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
505
506	c coordonnees du centre du zoom
507	call coordij(clon,clat,ilon,ilat)
508	c aire de la maille au centre du zoom
509	dxdy_min=dxdys(ilon,ilat)
510	c dxdy maximale de la maille
511	dxdy_max=0.
512	do j=1,jjp1
513	do i=1,iip1
514	dxdy_max=max(dxdy_max,dxdys(i,j))
515	enddo
516	enddo
517
518	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
519	print*,'ATTENTION modele peu zoome'
520	print*,'ATTENTION on prend une constante de guidage cste'
521	gamma=0.
522	else
523	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
524	print*,'gamma=',gamma
525	if (gamma.lt.1.e-5) then
526	print*,'gamma =',gamma,'<1e-5'
527	stop
528	endif
529	print*,'gamma=',gamma
530	gamma=log(0.5)/log(gamma)
531	endif
532	endif
533
534	alphamin=factt/taumax
535	alphamax=factt/taumin
536
537	do j=1,pjm
538	do i=1,pim
539	if (type.eq.1) then
540	dxdy_=dxdys(i,j)
541	zlat=rlatu(j)*180./pi
542	elseif (type.eq.2) then
543	dxdy_=dxdyu(i,j)
544	zlat=rlatu(j)*180./pi
545	elseif (type.eq.3) then
546	dxdy_=dxdyv(i,j)
547	zlat=rlatv(j)*180./pi
548	endif
549	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
550	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
551	alpha(i,j)=alphamin
552	else
553	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
554	xi=min(xi,1.)
555	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
556	alpha(i,j)=xialphamin+(1.-xi)alphamax
557	else
558	alpha(i,j)=0.
559	endif
560	endif
561	enddo
562	enddo
563
564
565	return
566	end

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

Download in other formats:

Original Format