Context Navigation

guide_p.F @ 1155

Last change on this file since 1155 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: 17.6 KB

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

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

Download in other formats:

Original Format