Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

guide.F @ 5426

Last change on this file since 5426 was 633, checked in by (none), 20 years ago
This commit was manufactured by cvs2svn to create branch 'LMDZ4_par_0'.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.5 KB

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

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

Download in other formats:

Original Format