Context Navigation

guide_p.F @ 3947

Last change on this file since 3947 was 709, checked in by Laurent Fairhead, 18 years ago
Nouvelles versions de la dynamique YM LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.8 KB

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

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

Download in other formats:

Original Format