Context Navigation

moldiff.F @ 2883

Last change on this file since 2883 was 2615, checked in by romain.vande, 3 years ago
LMDZ_MARS RV : Open_MP; Put all the "save" variables as "!$OMP THREADPRIVATE" in aeronomars
File size: 18.0 KB

Rev	Line
[1047]	1	subroutine moldiff(ngrid,nlayer,nq,
	2	& pplay,pplev,pt,pdt,pq,pdq,ptimestep,
[38]	3	& zzlay,pdteuv,pdtconduc,pdqdiff)
	4
[1047]	5	use tracer_mod, only: igcm_co2, igcm_co, igcm_o, igcm_o1d,
[1036]	6	& igcm_o2, igcm_o3, igcm_h, igcm_h2, igcm_oh,
	7	& igcm_ho2, igcm_h2o2, igcm_n2, igcm_ar,
	8	& igcm_h2o_vap, mmol
[1047]	9	use conc_mod, only: rnew, mmean
[1226]	10	USE comcstfi_h
[38]	11	implicit none
	12
	13	c
	14	c Input/Output
	15	c
[1047]	16	integer,intent(in) :: ngrid ! number of atmospheric columns
	17	integer,intent(in) :: nlayer ! number of atmospheric layers
	18	integer,intent(in) :: nq ! number of advected tracers
[38]	19	real ptimestep
[1047]	20	real pplay(ngrid,nlayer)
	21	real zzlay(ngrid,nlayer)
	22	real pplev(ngrid,nlayer+1)
	23	real pq(ngrid,nlayer,nq)
	24	real pdq(ngrid,nlayer,nq)
	25	real pt(ngrid,nlayer)
	26	real pdt(ngrid,nlayer)
	27	real pdteuv(ngrid,nlayer)
	28	real pdtconduc(ngrid,nlayer)
	29	real pdqdiff(ngrid,nlayer,nq)
[38]	30	c
	31	c Local
	32	c
	33
[414]	34	integer,parameter :: ncompmoldiff = 14
	35	real hco2(ncompmoldiff),ho
	36
[38]	37	integer ig,nz,l,n,nn,iq
	38	real del1,del2, tmean ,dalfinvdz, d
	39	real hh,dcoef,dcoef1,ptfac, ntot, dens, dens2, dens3
[1047]	40	real hp(nlayer)
	41	real tt(nlayer)
	42	real qq(nlayer,ncompmoldiff)
	43	real dmmeandz(nlayer)
	44	real qnew(nlayer,ncompmoldiff)
	45	real zlocal(nlayer)
[414]	46	real alf(ncompmoldiff-1,ncompmoldiff-1)
[1047]	47	real alfinv(nlayer,ncompmoldiff-1,ncompmoldiff-1)
[414]	48	real indx(ncompmoldiff-1)
[1047]	49	real b(nlayer,ncompmoldiff-1)
[414]	50	real y(ncompmoldiff-1,ncompmoldiff-1)
[1047]	51	real aa(nlayer,ncompmoldiff-1,ncompmoldiff-1)
	52	real bb(nlayer,ncompmoldiff-1,ncompmoldiff-1)
	53	real cc(nlayer,ncompmoldiff-1,ncompmoldiff-1)
	54	real atri(nlayer-2)
	55	real btri(nlayer-2)
	56	real ctri(nlayer-2)
	57	real rtri(nlayer-2)
	58	real qtri(nlayer-2)
[414]	59	real alfdiag(ncompmoldiff-1)
	60	real wi(ncompmoldiff), flux(ncompmoldiff), pote
[38]	61
	62	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	63	c tracer numbering in the molecular diffusion
	64	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	65	c Atomic oxygen must always be the LAST species of the list as
	66	c it is the dominant species at high altitudes.
	67	integer,parameter :: i_co = 1
	68	integer,parameter :: i_n2 = 2
	69	integer,parameter :: i_o2 = 3
	70	integer,parameter :: i_co2 = 4
	71	integer,parameter :: i_h2 = 5
	72	integer,parameter :: i_h = 6
	73	integer,parameter :: i_oh = 7
	74	integer,parameter :: i_ho2 = 8
	75	integer,parameter :: i_h2o = 9
	76	integer,parameter :: i_h2o2 = 10
	77	integer,parameter :: i_o1d = 11
	78	integer,parameter :: i_o3 = 12
	79	integer,parameter :: i_ar = 13
	80	integer,parameter :: i_o = 14
	81
	82	! Tracer indexes in the GCM:
	83	integer,save :: g_co2=0
	84	integer,save :: g_co=0
	85	integer,save :: g_o=0
	86	integer,save :: g_o1d=0
	87	integer,save :: g_o2=0
	88	integer,save :: g_o3=0
	89	integer,save :: g_h=0
	90	integer,save :: g_h2=0
	91	integer,save :: g_oh=0
	92	integer,save :: g_ho2=0
	93	integer,save :: g_h2o2=0
	94	integer,save :: g_n2=0
	95	integer,save :: g_ar=0
	96	integer,save :: g_h2o=0
	97
[414]	98	integer,save :: gcmind(ncompmoldiff) ! array of GCM indexes
[38]	99	integer ierr
	100
	101	logical,save :: firstcall=.true.
[414]	102	real abfac(ncompmoldiff)
	103	real,save :: dij(ncompmoldiff,ncompmoldiff)
[38]	104
[2615]	105	!$OMP THREADPRIVATE(g_co2,g_co,g_o,g_o1d,g_o2,g_o3,g_h,g_h2)
	106	!$OMP THREADPRIVATE(g_oh,g_ho2,g_h2o2,g_n2,g_ar,g_h2o,gcmind)
	107	!$OMP THREADPRIVATE(firstcall,dij)
	108
[38]	109	! Initializations at first call
	110	if (firstcall) then
	111	call moldiffcoeff(dij)
	112	print*,'MOLDIFF EXO'
	113
	114	! identify the indexes of the tracers we'll need
	115	g_co2=igcm_co2
	116	if (g_co2.eq.0) then
	117	write(,) "moldiff: Error; no CO2 tracer !!!"
	118	stop
	119	endif
	120	g_co=igcm_co
	121	if (g_co.eq.0) then
	122	write(,) "moldiff: Error; no CO tracer !!!"
	123	stop
	124	endif
	125	g_o=igcm_o
	126	if (g_o.eq.0) then
	127	write(,) "moldiff: Error; no O tracer !!!"
	128	stop
	129	endif
	130	g_o1d=igcm_o1d
	131	if (g_o1d.eq.0) then
	132	write(,) "moldiff: Error; no O1D tracer !!!"
	133	stop
	134	endif
	135	g_o2=igcm_o2
	136	if (g_o2.eq.0) then
	137	write(,) "moldiff: Error; no O2 tracer !!!"
	138	stop
	139	endif
	140	g_o3=igcm_o3
	141	if (g_o3.eq.0) then
	142	write(,) "moldiff: Error; no O3 tracer !!!"
	143	stop
	144	endif
	145	g_h=igcm_h
	146	if (g_h.eq.0) then
	147	write(,) "moldiff: Error; no H tracer !!!"
	148	stop
	149	endif
	150	g_h2=igcm_h2
	151	if (g_h2.eq.0) then
	152	write(,) "moldiff: Error; no H2 tracer !!!"
	153	stop
	154	endif
	155	g_oh=igcm_oh
	156	if (g_oh.eq.0) then
	157	write(,) "moldiff: Error; no OH tracer !!!"
	158	stop
	159	endif
	160	g_ho2=igcm_ho2
	161	if (g_ho2.eq.0) then
	162	write(,) "moldiff: Error; no HO2 tracer !!!"
	163	stop
	164	endif
	165	g_h2o2=igcm_h2o2
	166	if (g_h2o2.eq.0) then
	167	write(,) "moldiff: Error; no H2O2 tracer !!!"
	168	stop
	169	endif
	170	g_n2=igcm_n2
	171	if (g_n2.eq.0) then
	172	write(,) "moldiff: Error; no N2 tracer !!!"
	173	stop
	174	endif
	175	g_ar=igcm_ar
	176	if (g_ar.eq.0) then
	177	write(,) "moldiff: Error; no AR tracer !!!"
	178	stop
	179	endif
	180	g_h2o=igcm_h2o_vap
	181	if (g_h2o.eq.0) then
	182	write(,) "moldiff: Error; no water vapor tracer !!!"
	183	stop
	184	endif
	185
	186	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	187	c fill array to relate local indexes to gcm indexes
	188	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	189
	190	gcmind(i_co) = g_co
	191	gcmind(i_n2) = g_n2
	192	gcmind(i_o2) = g_o2
	193	gcmind(i_co2) = g_co2
	194	gcmind(i_h2) = g_h2
	195	gcmind(i_h) = g_h
	196	gcmind(i_oh) = g_oh
	197	gcmind(i_ho2) = g_ho2
	198	gcmind(i_h2o) = g_h2o
	199	gcmind(i_h2o2) = g_h2o2
	200	gcmind(i_o1d) = g_o1d
	201	gcmind(i_o3) = g_o3
	202	gcmind(i_o) = g_o
	203	gcmind(i_ar) = g_ar
	204
	205	firstcall= .false.
	206	endif ! of if (firstcall)
	207
	208
	209
	210	c
	211	cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	212
[1047]	213	nz=nlayer
[38]	214
[1047]	215	do ig=1,ngrid
[38]	216
	217	do l=2,nz-1
	218	tt(l)=pt(ig,l)+pdt(ig,l)*ptimestep
	219	& +pdteuv(ig,l)*ptimestep
	220	& +pdtconduc(ig,l)*ptimestep
	221
[414]	222	do nn=1,ncompmoldiff
[38]	223	qq(l,nn)=pq(ig,l,gcmind(nn))+pdq(ig,l,gcmind(nn))*ptimestep
	224	qq(l,nn)=max(qq(l,nn),1.e-30)
	225	enddo
	226	hp(l)=-log(pplay(ig,l+1)/pplay(ig,l-1))
	227	dmmeandz(l)=(mmean(ig,l+1)-mmean(ig,l-1))/hp(l)
	228	enddo
	229
	230	tt(1)=pt(ig,1) +pdt(ig,1)*ptimestep
	231	& +pdteuv(ig,1)*ptimestep
	232	& +pdtconduc(ig,1)*ptimestep
	233	tt(nz)=pt(ig,nz)+pdt(ig,nz)*ptimestep
	234	& +pdteuv(ig,nz)*ptimestep
	235	& +pdtconduc(ig,nz)*ptimestep
	236
[414]	237	do nn=1,ncompmoldiff
[38]	238	qq(1,nn)=pq(ig,1,gcmind(nn))+pdq(ig,1,gcmind(nn))*ptimestep
	239	qq(nz,nn)=pq(ig,nz,gcmind(nn))+pdq(ig,nz,gcmind(nn))*ptimestep
	240	qq(1,nn)=max(qq(1,nn),1.e-30)
	241	qq(nz,nn)=max(qq(nz,nn),1.e-30)
	242	enddo
	243	hp(1)=-log(pplay(ig,2)/pplay(ig,1))
	244	dmmeandz(1)=(-3.mmean(ig,1)+4.mmean(ig,2)
	245	& -mmean(ig,3))/(2.*hp(1))
	246	hp(nz)=-log(pplay(ig,nz)/pplay(ig,nz-1))
	247	dmmeandz(nz)=(3.mmean(ig,nz)-4.mmean(ig,nz-1)
	248	& +mmean(ig,nz-2))/(2.*hp(nz))
	249	c
	250	c Setting-up matrix of alfa coefficients from Dickinson and Ridley 1972
	251	c
	252	do l=1,nz
	253	if(abs(dmmeandz(l)) .lt. 1.e-5) dmmeandz(l)=0.0
	254	hh=rnew(ig,l)*tt(l)/g
	255	ptfac=(1.e5/pplay(ig,l))(tt(l)/273)*1.75
	256	ntot=pplay(ig,l)/(1.381e-23*tt(l)) ! in #/m3
	257
[414]	258	do nn=1,ncompmoldiff-1 ! rows
[38]	259	alfdiag(nn)=0.
	260	dcoef1=dij(nn,i_o)*ptfac
[414]	261	do n=1,ncompmoldiff-1 ! columns
[38]	262	y(nn,n)=0.
	263	dcoef=dij(nn,n)*ptfac
	264	alf(nn,n)=qq(l,nn)/ntot/1.66e-27
	265	& (1./(mmol(gcmind(n))dcoef)-1./(mmol(g_o)*dcoef1))
	266	alfdiag(nn)=alfdiag(nn)
	267	& +(1./(mmol(gcmind(n))dcoef)-1./(mmol(g_o)dcoef1))
	268	& *qq(l,n)
	269	enddo
	270	dcoef=dij(nn,nn)*ptfac
	271	alfdiag(nn)=alfdiag(nn)
	272	& -(1./(mmol(gcmind(nn))dcoef)-1./(mmol(g_o)dcoef1))
	273	& *qq(l,nn)
	274	alf(nn,nn)=-(alfdiag(nn)
	275	& +1./(mmol(g_o)*dcoef1))/ntot/1.66e-27
	276	y(nn,nn)=1.
	277	b(l,nn)=-(dmmeandz(l)/mmean(ig,l)
	278	& +mmol(gcmind(nn))/mmean(ig,l)-1.)
	279	enddo
	280
	281	c
	282	c Inverting the alfa matrix
	283	c
[658]	284	call ludcmp_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,d,ierr)
[38]	285
	286	c TEMPORAIRE *****************************
	287	if (ierr.ne.0) then
[658]	288	write(,)'In moldiff: Problem in LUDCMP_SP with matrix alf'
[38]	289	write(,) 'Singular matrix ?'
	290	c write(,) 'Matrix alf = ', alf
	291	write(,) 'ig, l=',ig, l
	292	write(,) 'No molecular diffusion this time !'
[1047]	293	pdqdiff(1:ngrid,1:nlayer,1:nq)=0
[38]	294	return
	295	c stop
	296	end if
	297	c *******************************************
[414]	298	do n=1,ncompmoldiff-1
[690]	299	call lubksb_sp(alf,ncompmoldiff-1,ncompmoldiff-1,indx,y(1,n))
[414]	300	do nn=1,ncompmoldiff-1
[38]	301	alfinv(l,nn,n)=y(nn,n)/hh
	302	enddo
	303	enddo
	304	enddo
	305
	306	c
	307	c Calculating coefficients of the system
	308	c
	309
	310	c zlocal(1)=-log(pplay(ig,1)/pplev(ig,1))* Rnew(ig,1)*tt(1)/g
	311	zlocal(1)=zzlay(ig,1)
	312
	313	do l=2,nz-1
	314	del1=hp(l)pplay(ig,l)/(gptimestep)
	315	del2=(hp(l)/2)*2pplay(ig,l)/(g*ptimestep)
[414]	316	do nn=1,ncompmoldiff-1
	317	do n=1,ncompmoldiff-1
[38]	318	dalfinvdz=(alfinv(l+1,nn,n)-alfinv(l-1,nn,n))/hp(l)
	319	aa(l,nn,n)=-dalfinvdz/del1+alfinv(l,nn,n)/del2
	320	& +alfinv(l-1,nn,n)*b(l-1,n)/del1
	321	bb(l,nn,n)=-2.*alfinv(l,nn,n)/del2
	322	cc(l,nn,n)=dalfinvdz/del1+alfinv(l,nn,n)/del2
	323	& -alfinv(l+1,nn,n)*b(l+1,n)/del1
	324	enddo
	325	enddo
	326
	327	c tmean=tt(l)
	328	c if(tt(l).ne.tt(l-1))
	329	c & tmean=(tt(l)-tt(l-1))/log(tt(l)/tt(l-1))
	330	c zlocal(l)= zlocal(l-1)
	331	c & -log(pplay(ig,l)/pplay(ig,l-1))rnew(ig,l)tmean/g
	332	zlocal(l)=zzlay(ig,l)
	333	enddo
	334
	335	c zlocal(nz)= zlocal(nz-1)
	336	c & -log(pplay(ig,nz)/pplay(ig,nz-1))rnew(ig,nz)tmean/g
	337	zlocal(nz)=zzlay(ig,nz)
	338
	339	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	340	c Escape velocity from Jeans equation for the flux of H and H2
	341	c (Hunten 1973, eq. 5)
	342
[414]	343	do n=1,ncompmoldiff
[38]	344	wi(n)=1.
	345	flux(n)=0.
	346	abfac(n)=1.
	347	enddo
	348
	349	dens=pplay(ig,nz)/(rnew(ig,nz)*tt(nz))
	350	c
	351	c For H:
	352	c
	353	pote=(3398000.+zlocal(nz))/
	354	& (1.381e-23tt(nz)/(1.6605e-27mmol(g_h)*g))
	355	wi(i_h)=sqrt(2.1.381e-23tt(nz)/(1.6605e-27*mmol(g_h)))
	356	& /(2.sqrt(3.1415))(1.+pote)*exp(-pote)
	357	flux(i_h)=qq(nz,i_h)dens/(1.6605e-27mmol(g_h))*wi(i_h)
	358	flux(i_h)=flux(i_h)*1.6606e-27
	359	abfac(i_h)=0.
	360	c
	361	c For H2:
	362	c
	363	pote=(3398000.+zlocal(nz))/
	364	& (1.381e-23tt(nz)/(1.6605e-27mmol(g_h2)*g))
	365	wi(i_h2)=sqrt(2.1.381e-23tt(nz)/(1.6605e-27*mmol(g_h2)))
	366	& /(2.sqrt(3.1415))(1.+pote)*exp(-pote)
	367	flux(i_h2)=qq(nz,i_h2)dens/(1.6605e-27mmol(g_h2))*wi(i_h2)
	368	flux(i_h2)=flux(i_h2)*1.6606e-27
	369	abfac(i_h2)=0.
	370
	371	c ********* TEMPORAIRE : no escape for h and h2
	372	c do n=1,ncomptot
	373	c wi(n)=1.
	374	c flux(n)=0.
	375	c abfac(n)=1.
	376	c enddo
	377	c ********************************************
	378
	379
	380	ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	381
	382	c
	383	c Setting coefficients for tridiagonal matrix and solving the system
	384	c
	385
[414]	386	do nn=1,ncompmoldiff-1
[38]	387	do l=2,nz-1
	388	atri(l-1)=aa(l,nn,nn)
	389	btri(l-1)=bb(l,nn,nn)+1.
	390	ctri(l-1)=cc(l,nn,nn)
	391	rtri(l-1)=qq(l,nn)
[414]	392	do n=1,ncompmoldiff-1
[38]	393	rtri(l-1)=rtri(l-1)-(aa(l,nn,n)*qq(l-1,n)
	394	& +bb(l,nn,n)*qq(l,n)
	395	& +cc(l,nn,n)*qq(l+1,n))
	396	enddo
	397	rtri(l-1)=rtri(l-1)+(aa(l,nn,nn)*qq(l-1,nn)
	398	& +bb(l,nn,nn)*qq(l,nn)
	399	& +cc(l,nn,nn)*qq(l+1,nn))
	400	enddo
	401
	402	c
	403	c Boundary conditions:
	404	c Escape flux for H and H2 at top
	405	c Diffusive equilibrium for the other species at top
	406	c Perfect mixing for all at bottom
	407	c
	408
	409	rtri(nz-2)=rtri(nz-2)
	410	& -ctri(nz-2)flux(nn)mmol(gcmind(nn))/(dens*wi(nn))
	411
	412	atri(nz-2)=atri(nz-2)
	413	& -abfac(nn)ctri(nz-2)/(3.-2.hp(nz)*b(nz,nn))
	414	btri(nz-2)=btri(nz-2)
	415	& +abfac(nn)4.ctri(nz-2)/(3.-2.hp(nz)b(nz,nn))
	416
	417	c rtri(1)=rtri(1)-atri(1)*qq(1,nn)
	418	btri(1)=btri(1)+atri(1)
	419
[658]	420	call tridag_sp(atri,btri,ctri,rtri,qtri,nz-2)
[38]	421
	422	do l=2,nz-1
	423	c qnew(l,nn)=qtri(l-1)
	424	qnew(l,nn)=max(qtri(l-1),1.e-30)
	425	enddo
	426
	427	qnew(nz,nn)=flux(nn)mmol(gcmind(nn))/(denswi(nn))
	428	& +abfac(nn)(4.qnew(nz-1,nn)-qnew(nz-2,nn))
	429	& /(3.-2.hp(nz)b(nz,nn))
	430	c qnew(1,nn)=qq(1,nn)
	431	qnew(1,nn)=qnew(2,nn)
	432
	433	qnew(nz,nn)=max(qnew(nz,nn),1.e-30)
	434	qnew(1,nn)=max(qnew(1,nn),1.e-30)
	435
	436	enddo ! loop on species
	437
	438	DO l=1,nz
	439	if(zlocal(l).gt.65000.)then
	440	pdqdiff(ig,l,g_o)=0.
[414]	441	do n=1,ncompmoldiff-1
[38]	442	pdqdiff(ig,l,gcmind(n))=(qnew(l,n)-qq(l,n))
	443	pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)-(qnew(l,n)-qq(l,n))
	444	pdqdiff(ig,l,gcmind(n))=pdqdiff(ig,l,gcmind(n))/ptimestep
	445	enddo
	446	pdqdiff(ig,l,g_o)=pdqdiff(ig,l,g_o)/ptimestep
	447	endif
	448	ENDDO
	449
	450	c do l=2,nz
	451	c do n=1,ncomptot-1
	452	c hco2(n)=qnew(l,n)pplay(ig,l)/(rnew(ig,l)tt(l)) /
	453	c & (qnew(l-1,n)pplay(ig,l-1)/(rnew(ig,l-1)tt(l-1)))
	454	c hco2(n)=-(zlocal(l)-zlocal(l-1))/log(hco2(n))/1000.
	455	c enddo
	456	c write(225,*),l,pt(1,l),(hco2(n),n=1,6)
	457	c write(226,*),l,pt(1,l),(hco2(n),n=7,12)
	458	c enddo
	459
	460	enddo ! ig loop
	461
	462	return
	463	end
	464
	465	c ********************************************************************
	466	c ********************************************************************
	467	c ********************************************************************
	468
[658]	469	subroutine tridag_sp(a,b,c,r,u,n)
	470	c parameter (nmax=100)
[38]	471	c dimension gam(nmax),a(n),b(n),c(n),r(n),u(n)
[658]	472	real gam(n),a(n),b(n),c(n),r(n),u(n)
[38]	473	if(b(1).eq.0.)then
[658]	474	stop 'tridag_sp: error: b(1)=0 !!! '
[38]	475	endif
	476	bet=b(1)
	477	u(1)=r(1)/bet
	478	do 11 j=2,n
	479	gam(j)=c(j-1)/bet
	480	bet=b(j)-a(j)*gam(j)
	481	if(bet.eq.0.) then
[658]	482	stop 'tridag_sp: error: bet=0 !!! '
[38]	483	endif
	484	u(j)=(r(j)-a(j)*u(j-1))/bet
	485	11 continue
	486	do 12 j=n-1,1,-1
	487	u(j)=u(j)-gam(j+1)*u(j+1)
	488	12 continue
	489	return
	490	end
	491
	492	c ********************************************************************
	493	c ********************************************************************
	494	c ********************************************************************
	495
[658]	496	SUBROUTINE LUBKSB_SP(A,N,NP,INDX,B)
[38]	497
	498	implicit none
	499
	500	integer i,j,n,np,ii,ll
	501	real sum
	502	real a(np,np),indx(np),b(np)
	503
	504	c DIMENSION A(NP,NP),INDX(N),B(N)
	505	II=0
	506	DO 12 I=1,N
	507	LL=INDX(I)
	508	SUM=B(LL)
	509	B(LL)=B(I)
	510	IF (II.NE.0)THEN
	511	DO 11 J=II,I-1
	512	SUM=SUM-A(I,J)*B(J)
	513	11 CONTINUE
	514	ELSE IF (SUM.NE.0.) THEN
	515	II=I
	516	ENDIF
	517	B(I)=SUM
	518	12 CONTINUE
	519	DO 14 I=N,1,-1
	520	SUM=B(I)
	521	IF(I.LT.N)THEN
	522	DO 13 J=I+1,N
	523	SUM=SUM-A(I,J)*B(J)
	524	13 CONTINUE
	525	ENDIF
	526	B(I)=SUM/A(I,I)
	527	14 CONTINUE
	528	RETURN
	529	END
	530
	531	c ********************************************************************
	532	c ********************************************************************
	533	c ********************************************************************
	534
[658]	535	SUBROUTINE LUDCMP_SP(A,N,NP,INDX,D,ierr)
[38]	536
	537	implicit none
	538
	539	integer n,np,nmax,i,j,k,imax
	540	real d,tiny,aamax
	541	real a(np,np),indx(np)
	542	integer ierr ! error =0 if OK, =1 if problem
	543
	544	PARAMETER (NMAX=100,TINY=1.0E-20)
	545	c DIMENSION A(NP,NP),INDX(N),VV(NMAX)
	546	real sum,vv(nmax),dum
	547
	548	D=1.
	549	DO 12 I=1,N
	550	AAMAX=0.
	551	DO 11 J=1,N
	552	IF (ABS(A(I,J)).GT.AAMAX) AAMAX=ABS(A(I,J))
	553	11 CONTINUE
	554	IF (AAMAX.EQ.0.) then
[658]	555	write(,) 'In moldiff: Problem in LUDCMP_SP with matrix A'
[38]	556	write(,) 'Singular matrix ?'
	557	c write(,) 'Matrix A = ', A
	558	c TO DEBUG :
	559	ierr =1
	560	return
	561	c stop
	562	END IF
	563
	564	VV(I)=1./AAMAX
	565	12 CONTINUE
	566	DO 19 J=1,N
	567	IF (J.GT.1) THEN
	568	DO 14 I=1,J-1
	569	SUM=A(I,J)
	570	IF (I.GT.1)THEN
	571	DO 13 K=1,I-1
	572	SUM=SUM-A(I,K)*A(K,J)
	573	13 CONTINUE
	574	A(I,J)=SUM
	575	ENDIF
	576	14 CONTINUE
	577	ENDIF
	578	AAMAX=0.
	579	DO 16 I=J,N
	580	SUM=A(I,J)
	581	IF (J.GT.1)THEN
	582	DO 15 K=1,J-1
	583	SUM=SUM-A(I,K)*A(K,J)
	584	15 CONTINUE
	585	A(I,J)=SUM
	586	ENDIF
	587	DUM=VV(I)*ABS(SUM)
	588	IF (DUM.GE.AAMAX) THEN
	589	IMAX=I
	590	AAMAX=DUM
	591	ENDIF
	592	16 CONTINUE
	593	IF (J.NE.IMAX)THEN
	594	DO 17 K=1,N
	595	DUM=A(IMAX,K)
	596	A(IMAX,K)=A(J,K)
	597	A(J,K)=DUM
	598	17 CONTINUE
	599	D=-D
	600	VV(IMAX)=VV(J)
	601	ENDIF
	602	INDX(J)=IMAX
	603	IF(J.NE.N)THEN
	604	IF(A(J,J).EQ.0.)A(J,J)=TINY
	605	DUM=1./A(J,J)
	606	DO 18 I=J+1,N
	607	A(I,J)=A(I,J)*DUM
	608	18 CONTINUE
	609	ENDIF
	610	19 CONTINUE
	611	IF(A(N,N).EQ.0.)A(N,N)=TINY
	612	ierr =0
	613	RETURN
	614	END
	615

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

Download in other formats:

Original Format