Context Navigation

guide.F @ 921

Last change on this file since 921 was 762, checked in by Laurent Fairhead, 18 years ago
Merge entre la version V3_conv et le HEAD YM, JG, LF
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
[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	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