Context Navigation

guide_p.F @ 1172

Last change on this file since 1172 was 1167, checked in by Laurent Fairhead, 16 years ago
Les options de compilation pour g95 et pgf90 avec makegcm FH une correction dans guide_p.F AJ LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 17.7 KB

Line
1	!
2	! $Id: guide_p.F 1167 2009-05-29 13:37:46Z ymeurdesoif $
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	cAJ-b 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	cAJ-e
266	endif !mpi_rank
267	cAJ-e
268	c Lecture du premier etat des reanalyses.
269	call Gather_Field(ps,ip1jmp1,1,0)
270
271	if (mpi_rank==0) then
272
273	call read_reanalyse(1,ps
274	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
275	qrea2(:,:)=max(qrea2(:,:),0.1)
276
277	endif
278
279	call Broadcast_Field(ucovrea2,ip1jmp1,llm,0)
280	call Broadcast_Field(vcovrea2,ip1jm,llm,0)
281	call Broadcast_Field(tetarea2,ip1jmp1,llm,0)
282	call Broadcast_Field(qrea2,ip1jmp1,llm,0)
283	call Broadcast_Field(masserea2,ip1jmp1,llm,0)
284	call Broadcast_Field(psrea2,ip1jmp1,1,0)
285
286
287
288	c-----------------------------------------------------------------------
289	c Debut de l'integration temporelle:
290	c ----------------------------------
291
292	endif ! first
293	c
294	C-----------------------------------------------------------------------
295	C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
296	C-----------------------------------------------------------------------
297
298	ditau=real(itau)
299	DDAY_step=real(day_step)
300	write(,)'ditau,dday_step'
301	write(,)ditau,dday_step
302	toto=4*ditau/dday_step
303	reste=toto-aint(toto)
304	c write(,)'toto,reste',toto,reste
305
306	if (reste.eq.0.) then
307	if (itau_test.eq.itau) then
308	write(,)'deuxieme passage de advreel a itau=',itau
309	stop
310	else
311	vcovrea1(:,:)=vcovrea2(:,:)
312	ucovrea1(:,:)=ucovrea2(:,:)
313	tetarea1(:,:)=tetarea2(:,:)
314	qrea1(:,:)=qrea2(:,:)
315
316	print*,'LECTURE REANALYSES, pas ',step_rea
317	s ,'apres ',count_no_rea,' non lectures'
318	step_rea=step_rea+1
319	itau_test=itau
320
321	call Gather_Field(ps,ip1jmp1,1,0)
322
323	if (mpi_rank==0) then
324	call read_reanalyse(step_rea,ps
325	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
326	qrea2(:,:)=max(qrea2(:,:),0.1)
327	endif
328
329	call Broadcast_Field(ucovrea2,ip1jmp1,llm,0)
330	call Broadcast_Field(vcovrea2,ip1jm,llm,0)
331	call Broadcast_Field(tetarea2,ip1jmp1,llm,0)
332	call Broadcast_Field(qrea2,ip1jmp1,llm,0)
333	call Broadcast_Field(masserea2,ip1jmp1,llm,0)
334	call Broadcast_Field(psrea2,ip1jmp1,1,0)
335
336	factt=dtvr*iperiod/daysec
337	ztau(:)=factt/max(alpha_T(:),1.e-10)
338
339	if (mpi_rank==0) then
340	if (ok_gradsfile) then
341	call wrgrads(igrads,1,aire ,'aire ','aire ' )
342	call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
343	call wrgrads(igrads,1,alpha_u,'au ','au ' )
344	call wrgrads(igrads,1,alpha_T,'at ','at ' )
345	call wrgrads(igrads,1,ztau,'taut ','taut ' )
346	call wrgrads(igrads,llm,ucov,'u ','u ' )
347	call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
348	call wrgrads(igrads,llm,teta,'T ','T ' )
349	call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
350	call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
351	call wrgrads(igrads,llm,q,'Q ','Q ' )
352	call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
353	endif !(ok_gradsfile) then
354	endif
355
356	endif
357	else
358	count_no_rea=count_no_rea+1
359	endif
360
361	C-----------------------------------------------------------------------
362	c Guidage
363	c x_gcm = a * x_gcm + (1-a) * x_reanalyses
364	C-----------------------------------------------------------------------
365
366	if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
367
368	ditau=real(itau)
369	dday_step=real(day_step)
370
371
372	tau=4*ditau/dday_step
373	tau=tau-aint(tau)
374
375	c ucov
376	ijb=ij_begin
377	ije=ij_end
378
379	if (guide_u) then
380	do l=1,llm
381	do ij=ijb,ije
382	a=(1.-tau)ucovrea1(ij,l)+tauucovrea2(ij,l)
383	ucov(ij,l)=(1.-alpha_u(ij))*ucov(ij,l)+alpha_u(ij)
384	$ alpha_pcor(l)a
385	if (first.and.ini_anal) ucov(ij,l)=a
386	enddo
387	enddo
388	endif
389
390	c teta
391	if (guide_T) then
392	do l=1,llm
393	do ij=ijb,ije
394	a=(1.-tau)tetarea1(ij,l)+tautetarea2(ij,l)
395	teta(ij,l)=(1.-alpha_T(ij))*teta(ij,l)+alpha_T(ij)
396	$ alpha_pcor(l)a
397	if (first.and.ini_anal) teta(ij,l)=a
398	enddo
399	enddo
400	endif
401
402	c P
403	if (guide_P) then
404	do ij=ijb,ije
405	a=(1.-tau)psrea1(ij)+taupsrea2(ij)
406	ps(ij)=(1.-alpha_P(ij))ps(ij)+alpha_P(ij)a
407	if (first.and.ini_anal) ps(ij)=a
408	enddo
409	CALL pression_p(ip1jmp1,ap,bp,ps,p)
410	CALL massdair_p(p,masse)
411	endif
412
413
414	c q
415	if (guide_Q) then
416	do l=1,llm
417	do ij=ijb,ije
418	a=(1.-tau)qrea1(ij,l)+tauqrea2(ij,l)
419	c hum relative en % -> hum specif
420	a=qsat(ij,l)a0.01
421	q(ij,l)=(1.-alpha_Q(ij))*q(ij,l)+alpha_Q(ij)
422	$ alpha_pcor(l)a
423	if (first.and.ini_anal) q(ij,l)=a
424	enddo
425	enddo
426	endif
427
428	c vcov
429	if (pole_sud) ije=ij_end-iip1
430
431	if (guide_v) then
432	do l=1,llm
433	do ij=ijb,ije
434	a=(1.-tau)vcovrea1(ij,l)+tauvcovrea2(ij,l)
435	vcov(ij,l)=(1.-alpha_v(ij))*vcov(ij,l)+alpha_v(ij)
436	$ alpha_pcor(l)a
437	if (first.and.ini_anal) vcov(ij,l)=a
438	enddo
439	enddo
440	endif
441
442	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
443	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
444	c call dump2d(iip1,jjp1,teta,'TETA ')
445
446	first=.false.
447
448	return
449	end
450
451	c=======================================================================
452	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
453	c=======================================================================
454
455	implicit none
456
457	#include "dimensions.h"
458	#include "paramet.h"
459	#include "comconst.h"
460	#include "comgeom2.h"
461	#include "guide.h"
462	#include "serre.h"
463
464	c arguments :
465	integer type
466	integer pim,pjm
467	real factt,taumin,taumax,dxdymin,dxdymax
468	real dxdy_,alpha(pim,pjm)
469	real dxdy_min,dxdy_max
470
471	c local :
472	real alphamin,alphamax,gamma,xi
473	save gamma
474	integer i,j,ilon,ilat
475
476	logical first
477	save first
478	data first/.true./
479
480	real cus(iip1,jjp1),cvs(iip1,jjp1)
481	real cuv(iip1,jjm),cvu(iip1,jjp1)
482	real zdx(iip1,jjp1),zdy(iip1,jjp1)
483
484	real zlat
485	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
486	common/comdxdy/dxdys,dxdyu,dxdyv
487
488	if (first) then
489	do j=2,jjm
490	do i=2,iip1
491	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
492	enddo
493	zdx(1,j)=zdx(iip1,j)
494	enddo
495	do j=2,jjm
496	do i=1,iip1
497	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
498	enddo
499	enddo
500	do i=1,iip1
501	zdx(i,1)=zdx(i,2)
502	zdx(i,jjp1)=zdx(i,jjm)
503	zdy(i,1)=zdy(i,2)
504	zdy(i,jjp1)=zdy(i,jjm)
505	enddo
506	do j=1,jjp1
507	do i=1,iip1
508	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
509	enddo
510	enddo
511	do j=1,jjp1
512	do i=1,iim
513	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
514	enddo
515	dxdyu(iip1,j)=dxdyu(1,j)
516	enddo
517	do j=1,jjm
518	do i=1,iip1
519	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
520	enddo
521	enddo
522
523	c call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
524	c call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
525	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
526
527	c coordonnees du centre du zoom
528	call coordij(clon,clat,ilon,ilat)
529	c aire de la maille au centre du zoom
530	dxdy_min=dxdys(ilon,ilat)
531	c dxdy maximale de la maille
532	dxdy_max=0.
533	do j=1,jjp1
534	do i=1,iip1
535	dxdy_max=max(dxdy_max,dxdys(i,j))
536	enddo
537	enddo
538
539	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
540	print*,'ATTENTION modele peu zoome'
541	print*,'ATTENTION on prend une constante de guidage cste'
542	gamma=0.
543	else
544	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
545	print*,'gamma=',gamma
546	if (gamma.lt.1.e-5) then
547	print*,'gamma =',gamma,'<1e-5'
548	stop
549	endif
550	gamma=log(0.5)/log(gamma)
551	if (gamma4) then
552	gamma=min(gamma,4.)
553	endif
554	print*,'gamma=',gamma
555	endif
556	endif
557
558	alphamin=factt/taumax
559	alphamax=factt/taumin
560
561	do j=1,pjm
562	do i=1,pim
563	if (type.eq.1) then
564	dxdy_=dxdys(i,j)
565	zlat=rlatu(j)*180./pi
566	elseif (type.eq.2) then
567	dxdy_=dxdyu(i,j)
568	zlat=rlatu(j)*180./pi
569	elseif (type.eq.3) then
570	dxdy_=dxdyv(i,j)
571	zlat=rlatv(j)*180./pi
572	endif
573	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
574	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
575	alpha(i,j)=alphamin
576	else
577	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
578	xi=min(xi,1.)
579	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
580	alpha(i,j)=xialphamin+(1.-xi)alphamax
581	else
582	alpha(i,j)=0.
583	endif
584	endif
585	enddo
586	enddo
587
588
589	return
590	end

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

Download in other formats:

Original Format