Context Navigation

guide.F @ 1159

Last change on this file since 1159 was 1146, checked in by Laurent Fairhead, 16 years ago
Réintegration dans le tronc des modifications issues de la branche LMDZ-dev comprises entre la révision 1074 et 1145 Validation: une simulation de 1 jour en séquentiel sur PC donne les mêmes résultats entre la trunk et la dev LF
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 16.4 KB

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

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

Download in other formats:

Original Format