Context Navigation

guide.F @ 757

Last change on this file since 757 was 630, checked in by Laurent Fairhead, 20 years ago
Import d'une version parallele de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 15.1 KB

Rev	Line
[630]	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	real aire_min,aire_max
	84	integer ilon,ilat
	85	real factt,ztau(ip1jmp1)
	86
	87	INTEGER itau,ij,l,i,j
	88	integer ncidt,varidpl,nlev,status
	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
	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
	213	ncidt=NCOPN('T.nc',NCNOWRIT,rcod)
	214	if (ncep) then
	215	status=NF_INQ_DIMID(ncidt,'LEVEL',rid)
	216	else
	217	status=NF_INQ_DIMID(ncidt,'PRESSURE',rid)
	218	endif
	219	status=NF_INQ_DIMLEN(ncidt,rid,nlev)
	220	print *,'nlev', nlev
	221	call ncclos(ncidt,rcod)
	222	c Lecture du premier etat des reanalyses.
	223	call read_reanalyse(1,ps
	224	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
	225	qrea2(:,:)=max(qrea2(:,:),0.1)
	226
	227
	228	c-----------------------------------------------------------------------
	229	c Debut de l'integration temporelle:
	230	c ----------------------------------
	231
	232	endif ! first
	233	c
	234	C-----------------------------------------------------------------------
	235	C----- IMPORTATION DES VENTS,PRESSION ET TEMPERATURE REELS:
	236	C-----------------------------------------------------------------------
	237
	238	ditau=real(itau)
	239	DDAY_step=real(day_step)
	240	write(,)'ditau,dday_step'
	241	write(,)ditau,dday_step
	242	toto=4*ditau/dday_step
	243	reste=toto-aint(toto)
	244	c write(,)'toto,reste',toto,reste
	245
	246	if (reste.eq.0.) then
	247	if (itau_test.eq.itau) then
	248	write(,)'deuxieme passage de advreel a itau=',itau
	249	stop
	250	else
	251	vcovrea1(:,:)=vcovrea2(:,:)
	252	ucovrea1(:,:)=ucovrea2(:,:)
	253	tetarea1(:,:)=tetarea2(:,:)
	254	qrea1(:,:)=qrea2(:,:)
	255
	256	print*,'LECTURE REANALYSES, pas ',step_rea
	257	s ,'apres ',count_no_rea,' non lectures'
	258	step_rea=step_rea+1
	259	itau_test=itau
	260	call read_reanalyse(step_rea,ps
	261	s ,ucovrea2,vcovrea2,tetarea2,qrea2,masserea2,psrea2,1,nlev)
	262	qrea2(:,:)=max(qrea2(:,:),0.1)
	263	factt=dtvr*iperiod/daysec
	264	ztau(:)=factt/max(alpha_T(:),1.e-10)
	265	call wrgrads(igrads,1,aire ,'aire ','aire ' )
	266	call wrgrads(igrads,1,dxdys ,'dxdy ','dxdy ' )
	267	call wrgrads(igrads,1,alpha_u,'au ','au ' )
	268	call wrgrads(igrads,1,alpha_T,'at ','at ' )
	269	call wrgrads(igrads,1,ztau,'taut ','taut ' )
	270	call wrgrads(igrads,llm,ucov,'u ','u ' )
	271	call wrgrads(igrads,llm,ucovrea2,'ua ','ua ' )
	272	call wrgrads(igrads,llm,teta,'T ','T ' )
	273	call wrgrads(igrads,llm,tetarea2,'Ta ','Ta ' )
	274	call wrgrads(igrads,llm,qrea2,'Qa ','Qa ' )
	275	call wrgrads(igrads,llm,q,'Q ','Q ' )
	276
	277	call wrgrads(igrads,llm,qsat,'QSAT ','QSAT ' )
	278
	279	endif
	280	else
	281	count_no_rea=count_no_rea+1
	282	endif
	283
	284	C-----------------------------------------------------------------------
	285	c Guidage
	286	c x_gcm = a * x_gcm + (1-a) * x_reanalyses
	287	C-----------------------------------------------------------------------
	288
	289	if(ini_anal) print*,'ATTENTION !!! ON PART DU GUIDAGE'
	290
	291	ditau=real(itau)
	292	dday_step=real(day_step)
	293
	294
	295	tau=4*ditau/dday_step
	296	tau=tau-aint(tau)
	297
	298	print*,'ATTENTION !!!! ON NE GUIDE QUE JUSQU A 15N'
	299
	300	c ucov
	301	if (guide_u) then
	302	do l=1,llm
	303	do ij=1,ip1jmp1
	304	a=(1.-tau)ucovrea1(ij,l)+tauucovrea2(ij,l)
	305	ucov(ij,l)=(1.-alpha_u(ij))ucov(ij,l)+alpha_u(ij)a
	306	if (first.and.ini_anal) ucov(ij,l)=a
	307	enddo
	308	enddo
	309	endif
	310
	311	c teta
	312	if (guide_T) then
	313	do l=1,llm
	314	do ij=1,ip1jmp1
	315	a=(1.-tau)tetarea1(ij,l)+tautetarea2(ij,l)
	316	teta(ij,l)=(1.-alpha_T(ij))teta(ij,l)+alpha_T(ij)a
	317	if (first.and.ini_anal) teta(ij,l)=a
	318	enddo
	319	enddo
	320	endif
	321
	322	c P
	323	if (guide_P) then
	324	do ij=1,ip1jmp1
	325	a=(1.-tau)psrea1(ij)+taupsrea2(ij)
	326	ps(ij)=(1.-alpha_P(ij))ps(ij)+alpha_P(ij)a
	327	if (first.and.ini_anal) ps(ij)=a
	328	enddo
	329	CALL pression(ip1jmp1,ap,bp,ps,p)
	330	CALL massdair(p,masse)
	331	endif
	332
	333
	334	c q
	335	if (guide_Q) then
	336	do l=1,llm
	337	do ij=1,ip1jmp1
	338	a=(1.-tau)qrea1(ij,l)+tauqrea2(ij,l)
	339	c hum relative en % -> hum specif
	340	a=qsat(ij,l)a0.01
	341	q(ij,l)=(1.-alpha_Q(ij))q(ij,l)+alpha_Q(ij)a
	342	if (first.and.ini_anal) q(ij,l)=a
	343	enddo
	344	enddo
	345	endif
	346
	347	c vcov
	348	if (guide_v) then
	349	do l=1,llm
	350	do ij=1,ip1jm
	351	a=(1.-tau)vcovrea1(ij,l)+tauvcovrea2(ij,l)
	352	vcov(ij,l)=(1.-alpha_v(ij))vcov(ij,l)+alpha_v(ij)a
	353	if (first.and.ini_anal) vcov(ij,l)=a
	354	enddo
	355	if (first.and.ini_anal) vcov(ij,l)=a
	356	enddo
	357	endif
	358
	359	c call dump2d(iip1,jjp1,tetarea1,'TETA REA 1 ')
	360	c call dump2d(iip1,jjp1,tetarea2,'TETA REA 2 ')
	361	c call dump2d(iip1,jjp1,teta,'TETA ')
	362
	363	first=.false.
	364
	365	return
	366	end
	367
	368	c=======================================================================
	369	subroutine tau2alpha(type,pim,pjm,factt,taumin,taumax,alpha)
	370	c=======================================================================
	371
	372	implicit none
	373
	374	#include "dimensions.h"
	375	#include "paramet.h"
	376	#include "comconst.h"
	377	#include "comgeom2.h"
	378	#include "guide.h"
	379	#include "serre.h"
	380
	381	c arguments :
	382	integer type
	383	integer pim,pjm
	384	real factt,taumin,taumax,dxdymin,dxdymax
	385	real dxdy_,alpha(pim,pjm)
	386	real dxdy_min,dxdy_max
	387
	388	c local :
	389	real alphamin,alphamax,gamma,xi
	390	save gamma
	391	integer i,j,ilon,ilat
	392
	393	logical first
	394	save first
	395	data first/.true./
	396
	397	real cus(iip1,jjp1),cvs(iip1,jjp1)
	398	real cuv(iip1,jjm),cvu(iip1,jjp1)
	399	real zdx(iip1,jjp1),zdy(iip1,jjp1)
	400
	401	real zlat
	402	real dxdys(iip1,jjp1),dxdyu(iip1,jjp1),dxdyv(iip1,jjm)
	403	common/comdxdy/dxdys,dxdyu,dxdyv
	404
	405	if (first) then
	406	do j=2,jjm
	407	do i=2,iip1
	408	zdx(i,j)=0.5*(cu(i-1,j)+cu(i,j))/cos(rlatu(j))
	409	enddo
	410	zdx(1,j)=zdx(iip1,j)
	411	enddo
	412	do j=2,jjm
	413	do i=1,iip1
	414	zdy(i,j)=0.5*(cv(i,j-1)+cv(i,j))
	415	enddo
	416	enddo
	417	do i=1,iip1
	418	zdx(i,1)=zdx(i,2)
	419	zdx(i,jjp1)=zdx(i,jjm)
	420	zdy(i,1)=zdy(i,2)
	421	zdy(i,jjp1)=zdy(i,jjm)
	422	enddo
	423	do j=1,jjp1
	424	do i=1,iip1
	425	dxdys(i,j)=sqrt(zdx(i,j)zdx(i,j)+zdy(i,j)zdy(i,j))
	426	enddo
	427	enddo
	428	do j=1,jjp1
	429	do i=1,iim
	430	dxdyu(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
	431	enddo
	432	dxdyu(iip1,j)=dxdyu(1,j)
	433	enddo
	434	do j=1,jjm
	435	do i=1,iip1
	436	dxdyv(i,j)=0.5*(dxdys(i,j)+dxdys(i+1,j))
	437	enddo
	438	enddo
	439
	440	call dump2d(iip1,jjp1,dxdys,'DX2DY2 SCAL ')
	441	call dump2d(iip1,jjp1,dxdyu,'DX2DY2 U ')
	442	call dump2d(iip1,jjp1,dxdyv,'DX2DY2 v ')
	443
	444	c coordonnees du centre du zoom
	445	call coordij(clon,clat,ilon,ilat)
	446	c aire de la maille au centre du zoom
	447	dxdy_min=dxdys(ilon,ilat)
	448	c dxdy maximale de la maille
	449	dxdy_max=0.
	450	do j=1,jjp1
	451	do i=1,iip1
	452	dxdy_max=max(dxdy_max,dxdys(i,j))
	453	enddo
	454	enddo
	455
	456	if (abs(grossismx-1.).lt.0.1.or.abs(grossismy-1.).lt.0.1) then
	457	print*,'ATTENTION modele peu zoome'
	458	print*,'ATTENTION on prend une constante de guidage cste'
	459	gamma=0.
	460	else
	461	gamma=(dxdy_max-2.*dxdy_min)/(dxdy_max-dxdy_min)
	462	print*,'gamma=',gamma
	463	if (gamma.lt.1.e-5) then
	464	print*,'gamma =',gamma,'<1e-5'
	465	stop
	466	endif
	467	print*,'gamma=',gamma
	468	gamma=log(0.5)/log(gamma)
	469	endif
	470	endif
	471
	472	alphamin=factt/taumax
	473	alphamax=factt/taumin
	474
	475	do j=1,pjm
	476	do i=1,pim
	477	if (type.eq.1) then
	478	dxdy_=dxdys(i,j)
	479	zlat=rlatu(j)*180./pi
	480	elseif (type.eq.2) then
	481	dxdy_=dxdyu(i,j)
	482	zlat=rlatu(j)*180./pi
	483	elseif (type.eq.3) then
	484	dxdy_=dxdyv(i,j)
	485	zlat=rlatv(j)*180./pi
	486	endif
	487	xi=((dxdy_max-dxdy_)/(dxdy_max-dxdy_min))**gamma
	488	c pour une grille reguliere, xi=xxx**0=1 -> alpha=alphamin
	489	xi=min(xi,1.)
	490	if(lat_min_guide.le.zlat .and. zlat.le.lat_max_guide) then
	491	alpha(i,j)=xialphamin+(1.-xi)alphamax
	492	else
	493	alpha(i,j)=0.
	494	endif
	495	enddo
	496	enddo
	497
	498
	499	return
	500	end

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

Download in other formats:

Original Format