c $Header: /usr/people/flittner/radtran/tomrad/pv2.2/src/RCS/itratepol.f,v 2.22 2000/08/30 16:38:09 flittner Exp $ subroutine itrate(z,ncp1,itmax,jmuz) c c********************************************************************** cccc c version aug 22,1977 c c purpose c c itrate is a fortran iv routine which performs the iterations c of the dave z matrix. c c method c c using the zeroth order approx for the z matrices from firstz c itrate iterates the z matrix for itmax times and saves the c last three iterated z values for each layer. c in addition the last three computed values of exponential c summations needed by evaltit are also saved. c c calling sequence c c call itrate (qr,z,tnstr,fs,ncp1,itmax,jmuz) c c variable type i/o description c -------- ---- --- ----------- c c qr(10) r*8 i factor used in computing ig c z(202) r*8 i zeroth order approx for z matrix c calculated by firstz c tnstr(15) r*8 o c c fs r*8 c ncp1 i*4 i # layers +1 c itmax i*4 i max # iterations c jmuz i*4 i index specifying present value of c solar zenith angle c c external references c evalit c c common areas referenced c c analysis and programming c k.f. klenk c k.f. klenk,p.m. smith sasc, aug 22 1977 c c modifications (date name purpose) c c last modified 11/19/92 by zia ahmad c modified for single iteration c c last modified 10/19/94 by zia ahmad c modified to print results after each itration c c last modified 03/08/95 by dave flittner c modified to call eva1pol after each iteration if switch c write_iter_file is TRUE c c last modified 03/14/95...dave flittner c purpose: set pressure scale height used in gravity correction c to rayleigh scattering od. Create new variable pscaleforg and c pass in common block consts. c c last mod. Sep. 6, 1998...def c purpose: combine several commons into switches.cmn. c c last mod. Mar. 14, 2000...def c purpose: use definition of itmax as the # of multiple scattering c events. So, mod. such that itmax=0 is single scattering. Also c various minor changes. cccc c*********************************************************************** c implicit none c passed integer ncp1,itmax,jmuz real*8 z(202) c commons include 'switches.cmn' include 'kmat.cmn' c local real*8 z1(202),z2(202),z3(202),z4(202), 1 zj1(202),zj2(202),b1(202),b2(202),b3(202),b4(202), 2 bj1(202),bj2(202),ekary(202,5),zma1(4,202), 3 zma2(4,202),zma3(4,202),zma4(4,202),ej1(4,202),ej2(4,202) real*8 tm1(4,202),tm2(4,202),tm3(4,202),tm4(4,202), 1 tm5(4,202),tm6(4,202),azsum,bzsum,czsum,dzsum,ejsum,fjsum integer itmm2,index,i,kk,kkk,it,m,j,itp,itpp c itmm2=itmax-2 index=1 it=0 itp=1 c c construct zeroth z- matrix and reduced source functions for c azimuth dependent terms c do 250 i=1,ncp1 z1(i)=z(i) z2(i)=0.d0 z3(i)=0.d0 z4(i)=z(i) zj1(i)=z(i) zj2(i)=z(i) tm1(1,i)=z1(i) tm2(1,i)=z2(i) tm3(1,i)=z3(i) tm4(1,i)=z4(i) tm5(1,i)=zj1(i) tm6(1,i)=zj2(i) 250 continue c if(write_iter_file.or.(jprint(4).eq.1)) !def & call eva1pol(tm1,tm2,tm3,tm4,tm5,tm6,jmuz,ncp1,itp) c if (itmax.lt.1) then c single scattering solution do 255 i = 1,ncp1 zma1(index,i) = z1(i) zma2(index,i) = z2(i) zma3(index,i) = z3(i) zma4(index,i) = z4(i) ej1(index,i) = z(i) ej2(index,i) = z(i) 255 continue else c c do first and higher order iterations c do 140 it=1,itmax do 145 m=1,ncp1 b1(m)=z1(m) b2(m)=z2(m) b3(m)=z3(m) b4(m)=z4(m) bj1(m)=zj1(m) bj2(m)=zj2(m) 145 continue do 410 i=1,ncp1 do 5100 kkk = 1,5 do 5101 kk = 1,ncp1 ekary(kk,kkk) = nek(kk,kkk,i) 5101 continue 5100 continue azsum=z(i) bzsum=0.0d0 czsum=0.0d0 dzsum=z(i) ejsum=z(i) fjsum=z(i) c do 290 j=1,ncp1 azsum=azsum+ekary(j,1)*b1(j)+ekary(j,2)*b3(j) bzsum=bzsum+ekary(j,1)*b2(j)+ekary(j,2)*b4(j) czsum=czsum+ekary(j,2)*b1(j)+ekary(j,3)*b3(j) dzsum=dzsum+ekary(j,2)*b2(j)+ekary(j,3)*b4(j) ejsum=ejsum+ekary(j,4)*bj1(j) fjsum=fjsum+ekary(j,5)*bj2(j) 290 continue c z1(i)=azsum z2(i)=bzsum z3(i)=czsum z4(i)=dzsum zj1(i)=ejsum zj2(i)=fjsum c c tm1(1,i)=z1(i) tm2(1,i)=z2(i) tm3(1,i)=z3(i) tm4(1,i)=z4(i) tm5(1,i)=zj1(i) tm6(1,i)=zj2(i) if(it.ge.itmm2)then c c save z matrix and zj1,zj2 of last three iterations c zma1(index,i)=z1(i) zma2(index,i)=z2(i) zma3(index,i)=z3(i) zma4(index,i)=z4(i) ej1(index,i)=zj1(i) ej2(index,i)=zj2(i) endif 410 continue if(it.ge.itmm2)index=index+1 itp=it+1 if(write_iter_file.or.(jprint(4).eq.1)) & call eva1pol(tm1,tm2,tm3,tm4,tm5,tm6,jmuz,ncp1,itp) 140 continue index=index-1 endif if(jprint(4).eq.1) write(33,1000) 1 (i,zma1(index,i),zma2(index,i),zma3(index,i),zma4(index,i), 2 ej1(index,i),ej2(index,i),i=1,ncp1) c -- call evalit to calculate the radiance at the top of the c -- atmosphere and the downward diffuse flux at the ground c itpp=itp+1 if(write_iter_file.or.(jprint(4).eq.1)) & call eva1pol(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1,itpp) call evalit (zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1) if(jprint(4).eq.1) write(33,2000) 1 (i,zma1(4,i),zma2(4,i),zma3(4,i),zma4(4,i), 2 ej1(4,i),ej2(4,i),i=1,ncp1) return 1000 format(1h1,t30,'debug print out from itrate ',10x, 1 'zmatrix after final iteration',////, 2 1h ,1x,' i',4x,'zma1',10x,'zma2',10x,'zma3',10x,'zma4', 3 10x,'ej1',11x,'ej2',4x,///,(1x,i3,1x,d12.6,2x,d12.6, 4 2x,d12.6,2x,d12.6,2x,d12.6,2x,d12.6)) 2000 format(///,1x,'source functions after extrapolation',//, 2 1h ,1x,' i',4x,'zma1',10x,'zma2',10x,'zma3',10x,'zma4', 3 10x,'ej1',11x,'ej2',4x,///,(1x,i3,1x,d12.6,2x,d12.6, 4 2x,d12.6,2x,d12.6,2x,d12.6,2x,d12.6)) end