Context Navigation

moldiff.F @ 3567

Last change on this file since 3567 was 3466, checked in by emillour, 2 months ago

Mars PCM:
More tidying in aeronomars:

remove unused "inv.F" and remove "dtridgl.F" which is not used here and is a duplicate of the "dtridgl" routine in phymars/swr_toon.F
turn aeronomars routines to modules, for those which aren't in modules yet.

File size: 17.9 KB

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

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

Download in other formats:

Original Format