Context Navigation

guide_p.F @ 1098

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

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

Download in other formats:

Original Format