Context Navigation

guide_p.F @ 1146

Last change on this file since 1146 was 1146, checked in by Laurent Fairhead, 15 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.6 KB

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

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

Download in other formats:

Original Format