Context Navigation

guide.F @ 1108

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

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

Download in other formats:

Original Format