Context Navigation

moldiff.F @ 1242

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

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

Download in other formats:

Original Format