Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

guide.F @ 1310

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

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

Download in other formats:

Original Format