program matscn c author: c. seftor c date: 12-mar-1992 c last modified 9/22/93 .... zia ahmad c purpose: to include the effect of molecular anisotropy in c rayleigh phase function c c purpose: used in conjunction with the mateer scan code to c generate radiance tables for a given ozone profile c and for given scan and azimuth angles (modified from c profil). c c last modified 03/07/95... dave flittner c purpose: compute spherical correction to the out going beam c when integrating after iteration process. Set with lsphout c logical variable. c c last modified 03/08/95... dave flittner c purpose: store Z1=I1/(-3/8*muo*sqrt(1-muo^2)*sqrt(1-mu^2)) and c Z2=I2/(3/32*(1-muo^2)*(1-mu^2)/mu) so can be used as input into c interpolation tables. 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. Also use logical switch lgcorrect c to impliment the gravity correction to the rayleigh scattering c optical depth. c implicit integer*4(i-n), real*8(a-h,o-z) integer*4 max_wave, max_scan, max_sza parameter (max_wave = 2500, max_scan = 6, max_sza = 10) parameter (max_num_iter = 12) integer*4 iter(max_num_iter) integer*4 nthet, nscan, nwavl, layr integer*4 jprint(10) real*8 wave_iter(max_num_iter) C C Setting up virtual memory for output parameters C real*8 wavel (max_wave) real*8 i0_out(max_scan, max_sza, max_wave) real*8 Z1_out(max_scan, max_sza, max_wave) real*8 Z2_out(max_scan, max_sza, max_wave) real*8 t_out(max_scan, max_sza, max_wave) real*8 sb_out(max_wave) C C Local Variables C real*8 x(101),p(101),t(101),tl(101),tsl(101),tt(202),tts(202), 1 ttl(202),dtsp(202),dtts(202),z(202) real*8 sbt real*8 tmp101(101), tmpcf(2) real*8 cofx(4,487),h(487),ps(487) real*8 tr, tensit real*8 xprf(11), tmpprf(11), alphp(12), betp(12), xsav(11) real*8 pres, theta(10), xdel real*8 wavelen, alpha0, beta, pnot, tautot real*8 psza(10) real*8 start_wave, end_wave real*8 scan(15) real*8 thta(15),azmth(8),pshold(101),hhold(101), 1 alb(11),thnot(10) real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),totint(15,8) real*8 nek(202,5,202) character prfnam*8 character wavnum*27 character file_name*200 logical lwav(12) logical write_iter_file common /cwrite_iter/write_iter_file !def real*8 ristar(15,11),tensig(11),total(15,8,11) C C C real*8 ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) common /emm/ematx,admatx,saz,caz,saz2,caz2,emu,emuz C C Functions C logical get_coef C C Common C common /kmat/ nek common /log/ lwav common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /thkns/ tsl,tt,tts,ttl,dtsp,dtts common /atmos/ alpha0,beta,code,pnot,tautot,wavelen common /input/ xprf, tmpprf, pres, wave_iter, iter, theta, * scan, start_wave, end_wave, nthet, nscan, num_iter common /inchr/ prfnam, wavnum common /in/ delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 1 lambda,layer common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,totint common /totals/ ristar,tensig,total common /prints/ jprint c real*8 ddd real*8 sbz,tnstrz,tnstrl,tnstrr real*8 e0za,e1za,e2za,ggz real*8 eitl,eitr,eitu real*8 eittl,eittr,eitotz,polz,ttz c common/buff1/sbz(9),tnstrz(9,15),tnstrl(15),tnstrr(15) common/buff2/e0za(9,15),e1za(9,15,8),e2za(9,15,8),ggz(9) common/buff3/eitl(15,8),eitr(15,8),eitu(15,8), ttl_pol(15), * ttr_pol(15) common/buff4/eittl(15,8,11),eittr(15,8,11),eitotz(15,8,11), 1 polz(15,8,11),ttz(9,15) c logical lsphout,lgcorrect !def common/csphout/lsphout !def common/cgcorrect/lgcorrect !def real*8 refdir(15) !def common/crefdir/refdir !def common/czfunc/Z1func(15),Z2func(15) !def c new statement real*8 rhon,gama,q,q1,q2,delp,sdp,q12s common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c if lsphout = True, then perform the out going beam using !def c spherical geometry. If False, then use flat atmosphere. !def data lsphout/.TRUE./ !def c if lgcorrect = True, then perform the gravity correciton to the !def c rayleigh scattering optical depth. !def data lgcorrect/.TRUE./ !def xdel = 0.0 jjprint = 0 c c --read in all input values c call iniclz call readin c num_arg = iargc() if (num_arg .lt. 4) then file_name = 'profil.dat' else call getarg(4, file_name) end if if (jprint(10).eq.1) 1 open (unit=23,file=file_name,recl=160,status='unknown') if (num_arg .lt. 3) then file_name = '/dev/null' else call getarg(3, file_name) end if open (unit=33,file=file_name,recl=132,status='unknown') if (num_arg .lt. 5) then file_name = 'sumry.dat' else call getarg(5, file_name) end if if (jprint(10).eq.1) 1 open (unit=9,file=file_name,status='unknown') write_iter_file = num_arg .gt. 5 if (write_iter_file) then call get_lun(iter_lun) call getarg(6, file_name) open (iter_lun, file=file_name, status='unknown', * form='unformatted') end if c c Make sure storage buffers do not overflow c if (max_scan .lt. imu) then print*,'Too many scan angles to store. Max num of scan angles' print*,'is limited to ',max_scan,' and ',imu,' were requested.' write(0, '(''Too many ssan angles to process'')') end if if (max_sza .lt. nthet) then print*,'Too many solar zenith angles to store. Max num of SZA' print*,'is limited to ',max_sza,' and ',nthet,' were requested.' write(0, '(''Too many scan angles to process'')') end if c jjprint = jprint(1) + jprint(2) + jprint(3) + jprint(4) 1 + jprint(5) + jprint(6) + jprint(7) + jprint(8) len = index(prfnam,' ') + 1 c c open (unit=34,file=prfnam(1:1)//prfnam(len:len+2)//'.OUT', c 1 form='unformatted',status='unknown') c c open (unit=35,file='TOM2.OUT',form='unformatted', c 1 status='unknown') c c if (jprint(10).eq.1) then write (23,6400) prfnam(1:8) write (23,6450) write (23,6475) pres write (23,6490) nthet write (23,6492) nscan write (23,6493) iazmth write (23,6494) nalb write (23,6496) start_wave, end_wave write (23,6497) write (23,6498) (jprint(j),j=1,8) write (23,6500) (j,j=1,11) write (23,6600) xprf write (23,6700) tmpprf write (23,6725) (j,j=1,10) write (23,6750) (theta(j),j=1,10) write (23,6751) (scan(j),j=1,nscan) write (23,6752) (azmth(j),j=1,iazmth) write (23,6753) (alb(j),j=1,nalb) endif call ctol(lwav) tmp0 = 273. pnot = pres ipath = 0 ipsudo = 1 jw = 0 iter_pos = 1 c --loop over all wavelengths num = 0 do while (get_coef(start_wave, end_wave, wavelen, alpha0, * tmpcf(1), tmpcf(2), beta, rhon) .and. num .lt. max_wave) if (alpha0 .eq. 0.0) alpha0 = 1.0d-10 num = num + 1 wavel(num) = wavelen if (jjprint.ne.0) write (33,7902) wavel(num) lambda = wavelen if (jprint(9).eq.1) write (6,6200) wavelen do while (iter_pos .lt. num_iter .and. * wave_iter(iter_pos+1) .lt. wavelen) iter_pos = iter_pos + 1 end do itmax = iter(iter_pos) c new statements if (jjprint .ne. 0) write(33,2121) wavelen 2121 format('matscn wavelength', 1pe12.4) call iniclz1 call emmat(imu,imuz,iazmth,thnot,thta,azmth) c end of new statements call relayr(xprf(2),xprf(1),tmpprf(2),tmpprf(1),x, 1 tmp101,pnot,xpnot,hhold,pshold) call opthik(x,p,t,tl,cofx,tmp0,tmpcf,tmp101) call dtaus(cofx,nc,xpnot,ncp1) c change the order of the subroutines, placing the call to slant after c dtaus instead of after reflex. !def lmax=82 layer=6 call slant(h,ps,cofx,fracin,lmax,layer,pshold,hhold) if(lsphout)then !def call exponesph(h,ps,cofx,lmax,nc,ncp1,imu) !def else !def do j=1,imu !def refdir(j)=0.0d0 !def if(emu(j).ne.0.d0) refdir(j)=dexp(-tautot/emu(j)) !def enddo !def call expone(nc,ncp1,imu) !def endif !def c --compute effective alpha call reflex(itmax,nek,ncp1) js = 0 do 2000 ijk=1,nthet js = js + 1 if (jprint(9).eq.1) write (6,6300) theta(ijk) psza(js) = theta(ijk) thenot=theta(ijk) cos_theta = cosd(theta(ijk)) cos_sin_theta = cos_theta*sind(theta(ijk)) if(ipath.eq.1) then thenot=dsqrt(1.d0-(1.d0/theta(ijk))**2)/cons thenot=dasin(thenot)/cnvrt endif if (ipsudo.eq.1) then call firstz(h,ps,z,thenot,fracin,lmax,ncp1) else call frstz2(z,thenot,fs,f1,f2,ncp1) endif call itrate(z,ncp1,itmax,ijk) C C Storing output parameters C do j = 1, imu i0_out(j, ijk, num) = eizero(j) sin_view = sind(scan(j)) cos_view = cosd(scan(j)) Z1_out(j, ijk, num) = Z1func(j) !def Z2_out(j, ijk, num) = Z2func(j) !def t_out(j, ijk, num) = (fs+gg)*(tnstr(j)+refdir(j)) !def if (write_iter_file) * write (iter_lun) wavelen, theta(ijk), scan(j), * i0_out(j, ijk, num)/pi, Z1_out(j, ijk, num)/pi, * Z2_out(j, ijk, num)/pi, t_out(j, ijk, num)/pi, * sb(4), itmax + 2, * (tnstrz(iter_num, j), iter_num = 1, itmax+2), * (e0za(iter_num, j), iter_num = 1, itmax+2), * (ggz(iter_num), iter_num = 1, itmax+2), * (sbz(iter_num), iter_num = 1, itmax+2), * fs, gg, exp(-tautot), iazmth, nalb, * (eitl(j, iaz), iaz = 1, iazmth), * (eitr(j, iaz), iaz = 1, iazmth), * (eitu(j, iaz), iaz = 1, iazmth), * ttl_pol(j), ttr_pol(j) end do if (jprint(8).eq.1) call tbprnt(thenot) 2000 continue sb_out(num) = sb(4) end do C C Issue a warning message if there was an attempt to do too many wavelengths C if (num .eq. max_wave) then print*,'Too many wavelenghts were being requested.' print*,'Maximun number of wavelengths is current set to ', * max_wave end if c c c nwavl = num if (num_arg .lt. 2) then file_name = 'tomnval.dat' else call getarg(2, file_name) end if if (jprint(10).eq.1) then open(unit=40, file=file_name ,form='unformatted', * status='unknown') write (40) nscan, nthet, nwavl write (40) (wavel(i), i = 1, nwavl) write (40) (((i0_out(i, j, k), k=1,nwavl), j = 1, nthet), * i=1, nscan) write (40) (((Z1_out(i, j, k), k=1,nwavl), j = 1, nthet), * i=1, nscan) write (40) (((Z2_out(i, j, k), k=1,nwavl), j = 1, nthet), * i=1, nscan) write (40) (((t_out(i, j, k), k=1,nwavl), j = 1, nthet), * i=1, nscan) write (40) (sb_out(j), j = 1, nwavl) endif C C Output stuff here C c 901 format(f4.1,1x,f5.1,1x,10(f7.2,1x)) if (write_iter_file) close(iter_lun) close (40) close (23) close (33) c close (34) 6200 format(' now doing calc for wavelength ',f7.2) 6300 format(' now doing calc for sza ',f7.2) 6400 format(' ***** ',a8,' *****') 6450 format(//,' ***** input *****') 6475 format(//,' pres....',f4.2) 6490 format(' nthet...',i4) 6492 format(' nscan...',i4) 6493 format(' nazim...',i4) 6494 format(' nalb....',i4) 6496 format(' wavelength range...',f8.2,' to ',f8.2) 6497 format(' evalit',' evalrf',' expone',' itrate',' opthik', 1 ' relayr',' slant',' tbprnt') 6498 format(i5,7i7) 6500 format(/,' sbuv layr',t10,11i7) 6550 format(//,'dndw tables for layer',i3,/) 6600 format(' ozone amt',t12,11f7.2) 6700 format(' layer tmp',t12,11f7.2) 6725 format(/,' sza #',t10,10i7) 6750 format(' sza',t12,10f7.2) 6751 format(' scan angles',t12,15f7.2) 6752 format(' azm angles',t12,8f7.2) 6753 format(' albedos ',t12,8f7.2) 6775 format(/,' wavl #',t13,i7,11i12) 7000 format(' wavelen',t11,12f12.3) 7100 format(/,' alpha0',t12,12d12.5) 7200 format(' beta',t12,12d12.5) 7300 format(' tmpcof(1)',t12,12d12.5) 7400 format(' tmpcof(2)',t12,12d12.5) 7425 format(' niter',t12,i7,11i12) 7450 format(//,' ***** output *****') 7700 format(f7.3,t12,12d12.5) 7850 format(' sza') 7800 format(//,'n values for alb of ',f6.2) 7801 format(' scn of ',f6.2,' azm of ',f6.2, 1//,' wvl',t11,12f12.3) 7900 format(//,'dndw values for alb of ',f6.2) 7901 format(' scn of ',f6.2,' azm of ',f6.2, 1//,' wvl',t11,12f12.3) 7902 format(//,' ***** Output for wavelength ',f7.2,' *****',//) end subroutine comrho(wavlen,rhon) c c subroutine comrho determins the 'rhon' (molecular depoarization factor) value for c for the given wavelength 'wavlen'. A table of rhons and wavelengths computed from c King's correction factor (see Bates...Planet. Space Sci. vol 32, 785-790, 1984) c are used for this purpose. A linear interpolation method is used for determining c rhon c implicit real*8 (a-h,o-z) real*8 rnb(24) integer*4 la(24) c data rnb/0.0454,0.0438,0.0422,0.0411,0.0400,0.0389,0.0378,0.0367, 1 0.0356,0.0345,0.0339,0.0328,0.0323,0.0317,0.0317,0.0312, 2 0.0306,0.0306,0.0301,0.0301,0.0301,0.0295,0.0295,0.0295/ c data la /2000,2050,2100,2150,2200,2250,2300,2400,2500,2600,2700, 1 2800,2900,3000,3100,3200,3300,3400,3500,3600,3700,3800, 2 3900,4000/ c c c write(33,222) wavlen 222 format('comrho wavelength',1pe12.4) lam=int(wavlen+0.5) do i=1,24 if(lam .gt. la(i) .and. lam .le. la(i+1)) then grad=(rnb(i+1)-rnb(i))/float(la(i+1)-la(i)) rhon=rnb(i)+grad*(wavlen-float(la(i))) endif enddo c c write(33,100) wavlen,rhon 100 format('wavelength rhon', 1p2e12.4) c return end subroutine ctol(lwav) c --author: c seftor c --date: 28-feb-91 c --purpose: convert character string to integer c --and set values of logical array character prfnam*8, wavnum*27 logical lwav(12) integer iwav, isza common /inchr/ prfnam, wavnum if (wavnum(2:2) .eq. '-') then read (wavnum(1:1),'(i2)') istrt if (wavnum(4:4) .eq. ' ') then read (wavnum(3:3),'(i2)') iend else read (wavnum(3:4),'(i2)') iend endif do 550 i = istrt,iend lwav(i) = .true. 550 continue else if (wavnum(3:3).eq.'-') then read (wavnum(1:2),'(i2)') istrt read (wavnum(4:5),'(i2)') iend do 575 i = istrt,iend lwav(i) = .true. 575 continue else j = 1 do 600 i = 2,21 if (wavnum(i:i) .eq. ',') then read (wavnum(j:i-1),'(i2)') iwav j = i + 1 lwav(iwav) = .true. if (wavnum(j:j) .eq. ' ') goto 700 else if (wavnum(i:i) .eq. ' ') then read (wavnum(j:i-1),'(i2)') iwav lwav(iwav) = .true. goto 700 endif 600 continue 700 continue endif return end function dexpi(x) c*********************************************************************** c this subroutine calculates the exponential integral quantities. c the exponential integrals are calculates in double precision, c and then truncated to single precision values. c computations of exponential integral with fifteen significant c figure accuracy. c c for negative values of the argument c range 1 greater than -1.0d-20. gamma + log( qabs of x ) c range 2 -1.d-20 to -1.5, 3 point gaussian quadrature. c range 3 -1.5 to -4.65, ratio of 2 polynomials each with 7 terms. c range 4 -4.65 to -12.0, ratio of 2 polynomials each with 6 terms. c range 5 -12.0 to -170.0, 12 point gauss-laguerre quadrature. c c for positive values of the argument c range 1 less than 1.0d-20, gamma + log(x) c range 2 1.0d-20 to 40.0, 12 point gaussian quadrature. c range 3 40.0 to 173.0, 12 point gauss-laguerre quadrature. c c c*********************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) real * 8 a1(3),b1(3),a2(7),b2(7),a3(6),b3(6),a4(12),b4(12), 1 a5(12),b5(12) data a1 / 0.1193095930415965d+0, 0.3306046932331323d+0, 1 0.4662347571015760d+0/ data b1 / 0.4679139345726910d+0, 0.3607615730461336d+0, 1 0.1713244923791703d+0/ data a2 / 1 0.2823912701457739d-1, 0.3052042817823067d+1, 1 0.2158885931211323d+2, 0.4104611319636983d+2, 2 0.2785527592726121d+2, 0.7133086969436196d+1, 3 0.5758931590224375d+0/ data b2 / 1 0.1000000000000000d+1, 0.1383869728490638d+2, 1 0.4880858183073600d+2, 0.6348804630786363d+2, 2 0.3441289899236299d+2, 0.7708964199043784d+1, 3 0.5758934565014882d+0/ data a3 / 1 0.7630772325814641d-1, 0.2123699219410890d+1, 1 0.4745350785776186d+1, 0.2966421696379266d+1, 2 0.6444800036068992d+0, 0.4295808082119383d-1/ data b3 / 1 0.1000000000000000d+1, 0.5278950049492932d+1, 1 0.7196111390658517d+1, 0.3567945294128107d+1, 2 0.6874380519301884d+0, 0.4295808112146861d-1/ data a4 / 1 0.1157221173580207d+0, 0.6117574845151307d+0, 1 0.1512610269776419d+1, 0.2833751337743507d+1, 2 0.4599227639418348d+1, 0.6844525453115177d+1, 3 0.9621316842456867d+1, 0.1300605499330635d+2, 4 0.1711685518746226d+2, 0.2215109037939701d+2, 5 0.2848796725098400d+2, 0.3709912104446692d+2/ data b4 / 1 0.2647313710554432d+0, 0.3777592758731380d+0, 1 0.2440820113198776d+0, 0.9044922221168093d-1, 2 0.2010238115463410d-1, 0.2663973541865316d-2, 3 0.2032315926629994d-3, 0.8365055856819799d-5, 4 0.1668493876540910d-6, 0.1342391030515004d-8, 5 0.3061601635035021d-11, 0.8148077467426242d-15/ data a5 / 1 0.3202844643130281d-1, 0.9555943373680816d-1, 1 0.1575213398480817d+0, 0.2168967538130226d+0, 2 0.2727107356944198d+0, 0.3240468259684878d+0, 3 0.3700620957892772d+0, 0.4100009929869515d+0, 4 0.4432077635022005d+0, 0.4691372760013664d+0, 5 0.4873642779856547d+0, 0.4975936099985107d+0/ data b5 / 1 0.1279381953467522d+0, 0.1258374563468283d+0, 1 0.1216704729278034d+0, 0.1155056680537256d+0, 2 0.1074442701159656d+0, 0.9761865210411389d-1, 3 0.8619016153195328d-1, 0.7334648141108031d-1, 4 0.5929858491543678d-1, 0.4427743881741981d-1, 5 0.2853138862893366d-1, 0.1234122979998720d-1/ 10 format(10x,'the argument of expi is very close to zero. it is', 1e25.16//) 20 format(10x,'the argument of expi is very large. it is',e25.16//) c c dexpi = 0.0d+00 if (x) 200,100,300 c 100 write (6,10) x 100 continue return 200 ax = dabs(x) if ( x .gt. -1.0d-20 ) go to 201 if ( x .gt. -1.5 ) go to 205 if ( x .gt. -4.65 ) go to 215 if ( x .gt. -12.0 ) go to 225 if ( x .gt. -170.0 ) go to 235 return 201 dexpi = dlog(ax) +0.57721566490153d+00 return 205 yy = dexp(-0.5 * ax) yz = dexp(a1(1)*ax) sumn=(1.0 -yy/yz)/(0.5 +a1(1))+(1.0 -yy*yz)/(0.5 -a1(1)) yz = dexp(a1(2)*ax) sumd=(1.0 -yy/yz)/(0.5 +a1(2))+(1.0 -yy*yz)/(0.5 -a1(2)) yz = dexp(a1(3)*ax) sumt=(1.0 -yy/yz)/(0.5 +a1(3))+(1.0 -yy*yz)/(0.5 -a1(3)) dexpi= -0.5 *(b1(1)*sumn+b1(2)*sumd+b1(3)*sumt) dexpi = dexpi + dlog(ax) + 0.57721566490153d+00 return 215 sumn =(((((a2(7)*ax+a2(6))*ax+a2(5))*ax+a2(4))*ax+a2(3))*ax+a2(2)) 1*ax+a2(1) sumd =(((((b2(7)*ax+b2(6))*ax+b2(5))*ax+b2(4))*ax+b2(3))*ax+b2(2)) 1*ax+b2(1) 218 dexpi= sumn/(sumd*x) dexpi = dexpi * dexp(x) return 225 sumn=((((a3(6)*ax+a3(5))*ax+a3(4))*ax+a3(3))*ax+a3(2))*ax+a3(1) sumd=((((b3(6)*ax+b3(5))*ax+b3(4))*ax+b3(3))*ax+b3(2))*ax+b3(1) go to 218 235 do 238 j = 1,12 238 dexpi=dexpi + b4(j)/(1.0 + a4(j)/ax) dexpi = (dexp(x)/ax)*(-dexpi) return 300 if ( x .le. 1.0d-20 ) go to 301 if ( x .le. 40.0 ) go to 305 if ( x .le. 173.0 ) go to 335 c write (6,20 ) x return 301 dexpi= dlog(x) + 0.57721566490153d+00 return 305 yy = dexp( 0.5 * x) do 310 j = 1,12 yz = dexp(-a5(j)* x ) dexpi =(( 1.0 - yy/yz)/(0.5 + a5(j)) + ( 1.0 - yy*yz)/ 1 ( 0.5 - a5(j)) )*b5(j) + dexpi 310 continue dexpi = -0.5 * dexpi + dlog(x) + 0.57721566490153d+00 return 335 do 338 j = 1,12 338 dexpi = dexpi + b4(j)/(1.0-a4(j)/x) dexpi = (dexp(x)/x) * dexpi return end subroutine dexpk(a,b) c c*********************************************************************** c cccc c subroutine dexpk c c purpose- c c calculate exponential integrals of orders 1 to 6 and the c elements of k1 matrix (dave's eqn. 4-11) c c method- c recursion relation for the exponential integrals is used c see handbook of math functions-abramowitz and stegun (p229) c c calling sequence- call dexpk(a,b) c c variables type i/o description c --------- ---- --- ----------- c a r*8 i abs(a-b) is the argument of the c b r*8 i exponential integrals c c result of the subroutine is returned in common block 'es' c e(6) r*8 o exponential integrals c eek(6) r*8 o elements of the k-matrix c c external references - dexpi c c modifications- by k.klenk 8/19/77 c c functions k1(1) and k1(2) are also calculated and are returned c in eek(4) and eek(5) respectively (dave's eqn. 4.7) c c last modified 9/4/93....zia ahmad c purpose: to account for molecular depolarization 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. cccc c*********************************************************************** c implicit real *8 (a-h,o-z) real *8 odan(6),e(6),eek(6) real *8 dexpi,dt,dum,dtemp,de(6) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /es/odan,e,eek c c new statement 9/4/93 common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c c*****calculate exponential integrals c dt = dabs(a-b) de(1) = - dexpi(-dt) dum = dexp(-dt) do 70 j = 2, 6 70 de(j) = odan(j)*(dum-dt*de(j-1)) do 80 i=1,6 80 e(i) = de(i) c c*****calculate the k matrix and functions k1(1) and k(2) c new statements added here 9/4/93 if(ipol .eq. 0)then eek(1) = 0.375*(e(1) + e(5)) eek(2) = c38sq2 * (de(3) - de(5)) eek(3) = 0.75 * (de(1)-2.d+00 * de(3) + de(5)) eek(4)=0.375*(e(1)+e(3)-2.0*e(5)) eek(5)=0.1875*(e(1)+2.0*e(3)+e(5)) else eek(1)=q12s*((sdp**2+q**2)*e(1)+2.0d0*q*e(3)+e(5)) eek(2)=q12s*delp*(q*(e(1)-e(3))+e(3)-e(5)) eek(3)=q12s*delp**2*(e(1)-2.0d0*e(3)+e(5)) eek(4)=0.375d0*q2*(e(1)+e(3)-2.0d0*e(5)) eek(5)=0.1875d0*q2*(e(1)+2.0d0*e(3)+e(5)) endif return end subroutine dexpk1(a,b) c c*********************************************************************** c cccc c subroutine dexpk1 c c purpose- c calculate the elements of k1 matrix averaged over the interval c a-b . (dave's eq. 4.11) c c method- c the averages are calculated using eq(5.1.26) in abramowitz and c segun. c c calling sequence - call dexpk1(a,b) c c c variables type i/o description c --------- ---- --- ----------- c a r*8 i abs(a-b) is the argument of the c b r*8 i exponential integrals c c result of the subroutine is returned in common block 'es' c e(6) r*8 o exponential integrals c eek(6) r*8 o elements of the k-matrix c c external references - dexpi c c modifications- by k.klenk 8/19/77 c c functions k1(1) and k1(2) are also calculated and are returned c in eek(4) and eek(5) respectively (dave's eqn. 4.7) c c last modified by zia ahmad 9/10/93 c purpose: to include the effect of molecular anisotropy 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*********************************************************************** c cccc c implicit integer*4(i-n),real*8 (a-h,o-z) real *8 odan(6),de(6),e(6),eek(6) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /es/odan,e,eek c c new statemnet common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c c c*****same as dexpk- calculate expon. integrals c dt = dabs(a-b) de(1) = - dexpi(-dt) dum = dexp(-dt) do 70 j = 2, 6 70 de(j) = odan(j)*(dum-dt*de(j-1)) do 80 i=1,6 80 e(i) = de(i) c c*****calculate elements of k1 matrix (azimuth independent term) c by averaging over dt c also calculate k1(1) and k1(2) functions for azimuth c dependent terms c new statements added here if(ipol.eq.0)then eek(1) = 0.375 *(1.2d+00-de(2)-de(6))/dt eek(2) = c38sq2*(2.0d+00/15.d+00-de(4)+de(6))/dt eek(3)= 0.75d+00*(8.d+00/15.d+00-de(2)+2.d+00*de(4)-de(6))/dt eek(4)=(c1415-e(2)-e(4)+2*e(6))*0.375d0/dt eek(5)=(c2815-e(2)-2*e(4)-e(6))*0.1875d0/dt else eek(1)=q12s*((sdp**2+q**2)*(1.0d0-e(2))+ 1 2.0d0*q*(1.0d0/3.0d0-e(4))+ 2 (1.0d0/5.0d0-e(6)) )/dt eek(2)=q12s*delp*(q*(1.0d0-e(2))+ 1 (1.0d0-q)*(1.0d0/3.0d0-e(4))- 2 (1.0d0/5.0d0-e(6)))/dt eek(3)=q12s*delp**2*((1.0d0-e(2))- 1 2.0d0*(1.0d0/3.0d0-e(4))+ 2 (1.0d0/5.0d0-e(6)))/dt eek(4)=q2*0.375d0*(c1415-e(2)-e(4)+2.0d0*e(6))/dt eek(5)=q2*0.1875d0*(c2815-e(2)-2.0d0*e(4)-e(6))/dt endif c end of new statements return end subroutine dtaus(cofx,nc,xpt,ncp1) c c*********************************************************************** c cccc c subroutine dtaus c c purpose- c 1. calculate the avg. optical thickness of each layer and of c each half layer at mid points between the layers c 2. find the no. of layers between the atmosphere top and a c specified bottom (pnot) c c input/output variables - c c variables type i/o description c --------- ---- --- ----------- c formal parameters c cofx(4,487) r*8 i spline interpolation coeff. c nc i*4 o no. of layers between the top and c the bottom(not incl. bottom layer) c ncp1 i*4 o nc+1 c part of common 'thkns' c dtsp(202) r*8 o avg. optical thickness of each layer c dtts(202) r*8 o avg. optical thickness of each half c layer at mid points between layers c c no external references 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. Perform gravity correction when c computing pnotb. c Also use logical switch lgcorrect c to impliment the gravity correction to the rayleigh scattering c optical depth. cccc c*********************************************************************** c implicit integer*4(i-n),real*8 (a-h,o-z) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101), 1 x(101),cofx(4,487),xpt,alb(11) real*8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) common /thkns/tsl,tt,tts,ttl,dtsp,dtts common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /atmos/ alpha0,beta,code,pnot,tautot,wavlen common /hedout/ rarray(824),psaray(82),iarray(3) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip logical lgcorrect !def common/cgcorrect/lgcorrect !def c c*****calculate dtsp and dtts and find nc c c if(lgcorrect)then pnotb = pnot * beta * (1.0d0-pscaleforg*dlog(pnot))**2 !def else pnotb = pnot * beta endif pnlog = dlog(pnotb) dtts(1) = 0.5*tts(2) dtsp(1) = dtts(1) do 25 i = 2, 202 dtsp(i) = 0.5*(tts(i+1) - tts(i-1)) dtts(i) = 0.5*(tts(i) - tts(i-1)) if (tts(i) .lt. pnotb) go to 25 ncp1 = i nc = i - 1 go to 26 25 continue ncp1 = 202 nc = 201 c c*****calculate tt and log(tt) of the bottom layer c*****(located at pnot)- by spline interpolation c 26 do 35 i = 80, 101 if (tsl(i) .lt. pnlog) go to 30 28 k = i - 1 dum1 = tsl(k+1) - pnlog dum2 = pnlog - tsl(k) dum3 = dum1*(cofx(1,k)*dum1**2 + cofx(3,k)) + dum2*(cofx(2,k)* 1dum2**2 + cofx(4,k)) ttl(ncp1) = dum3 tt(ncp1) = dexp(dum3) tautot = tt(ncp1) tts(ncp1) = pnotb go to 40 30 if (i .eq. 101) go to 28 35 continue c c*****correct the dtts and dtsp of the bottom layer c 40 dtts(ncp1) = 0.5*(tts(ncp1) - tts(nc)) dtsp(ncp1) = dtts(ncp1) dtsp(nc) = dtts(ncp1) + dtts(nc) alfaef = (tautot-pnotb)/xpt c --store data for archive tape rarray(1) = -12.0d0 rarray(203) = -12.0d0 do 100 i = 2,ncp1 j = i+202 rarray(i) = ttl(i) rarray(j)=dlog(tts(i)) 100 continue rarray(405) = code rarray(406) = pnot rarray(408) = alfaef rarray(409) = beta do 300 i = 1,82 psaray(i) = pshold(i) 300 continue do 200 i = 1,imuz j = 409 + i rarray(j) = thnot(i) 200 continue rarray(824) = wavlen iarray(1) = ncp1 iarray(2) = imuz iarray(3) = lambda c if (jprint(8).eq.1) c 1 write(33,6100)ncp1,pnot,(i,tt(i),tts(i),dtts(i), c 2 dtsp(i),i=1,ncp1) c6100 format(//,6hincp1=,i3,5x,5hpnot=f5.3/1h0,2(2x,1h1,8x, c 1 2htt,12x,3htts,10x,4hdtts,10x,4hdtsp,7x)// c 2 (1x,2(i3,4d14.4,4x))) return end subroutine emmat (imu,imuz,iazmth,thnot,thta,azmth) c c*********************************************************************** c cccc c subroutine emmat c version aug. 1, 1977 c purpose c c set up matrix m for each polar backscatter angle c and adjoint of m for each solar zenith angle. stores sines c and cosines of azimuth and 2*azimuth angles c c method c c see dave's paper may 13, 1964 c c calling sequence c c call emmat (imu,imuz,iazmth,thnot,thta,azmth,ematx, c admatx,saz,caz,saz2,caz2,emu,emuz) c c variable type i/o description c -------- ---- --- ----------- c c imu i*4 i # polar backscatter angles c imuz i*4 i # solar zenith angles c iazmth i*4 i # azimuth angles c thnot(n) r*8 i solar zenith angles (degs) c thta(n) r*8 i polar back scatter angles (degs) c azmth(n) r*8 i azimuth angles (degs) c c ematx(3,n) r*8 o m matrix c admatx(3,n) r*8 o adjoint matrix c saz(m) r*8 o sine of azimuth angle c caz(m) r*8 o cosine of azimuth angle c saz2(m) r*8 o 2*sine of azimuth angle c caz2(m) r*8 o 2*cosine of azimuth angle c emu(n) r*8 o cosine of polar backscatter angle c emuz(n) r*8 o cosine of solar zenith angles c c external references c none c c author c p. m. smith,sasc,aug 1, 1977 c modifications c (date name purpose) c last modified by zia ahmad 9/10/93 c purpose: to include the effect of molecular anisotropy 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. cccc c*********************************************************************** c implicit integer*4(i-n),real*8(a-h,o-z) common /emm/ematx,admatx,saz,caz,saz2,caz2,emu,emuz real*8 ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) c common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip c c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c real*8 thnot(10),thta(15),azmth(8) c c set up adjoint c c new statements added in do loops 100 and 200 c do 100 l=1,imuz thet=cnvrt*thnot(l) emuz(l)=dcos(thet) if(ipol.eq.0)then admatx(1,l)=emuz(l)*emuz(l) admatx(2,l)=1.0d0 admatx(3,l)=sq2*(1.0d0-admatx(1,l)) else admatx(1,l)=q+emuz(l)**2 admatx(2,l)=1.0d0+q admatx(3,l)=delp*(1.0d0-emuz(l)**2) endif 100 continue c c set up m matrix c do 200 l=1,imu thet=cnvrt*thta(l) emu(l)=dcos(thet) c write(33,555)l,thta(l),emu(l) 555 format('emmat',i4,1p2e12.4) if(ipol.eq.0)then ematx(1,l)=emu(l)*emu(l) ematx(2,l)=sq2*(1.0d0-ematx(1,l)) ematx(3,l)=1.0d0 else ematx(1,l)=q+emu(l)**2 ematx(2,l)=delp*(1.0d0-emu(l)**2) ematx(3,l)=1+q endif 200 continue c c store trig functions c do 300 i=1,iazmth phi=azmth(i)*cnvrt phi2=2.0d0*phi saz(i)=dsin(-phi) saz2(i)=dsin(-phi2) caz(i)=dcos(phi) caz2(i)=dcos(phi2) 300 continue return end subroutine eva1pol(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1,itz) c*********************************************************************** c c subroutine eva1pol c c purpose 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*********************************************************************** c implicit integer*4 (i-n), real*8(a-h,o-z) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101),alb(11) real*8 x(101) integer jprint(10) real*8 ztop(4,4),a(4),b(4),c(4),d(4),ztp1(4),ztp2(4),top(4), & dum(8),gr(4),gl(4), & zma1(4,202),zma2(4,202),zma3(4,202),zma4(4,202), & ej1(4,202),ej2(4,202),hold(4),xhold(4),eiold(15) c common /out/ sb(4),qr(11),tnstr(15),eizero(15),eiaz1(15,8), & eiaz2(15,8),gg,fs,totint(15,8) common /emm/ ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) c common /eks/ extmu(202,15),ek4(202),ek5(202) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip c common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, c 1 x,pshold,hhold,layer,thnot,alfaef, c 2 lambda common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /prints/ jprint common /in/delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of statement c common/buff2/e0za(9,15),e1za(9,15,8),e2za(9,15,8),ggz(9) data c316/0.1875d0/ c extrapolate source functions layer by layer c c itmaxp=itmax+1 do 300 i=1,ncp1 if (itz .le. itmaxp) then zma1(4,i) = zma1(1,i) zma2(4,i) = zma2(1,i) zma3(4,i) = zma3(1,i) zma4(4,i) = zma4(1,i) ej1(4,i) = ej1(1,i) ej2(4,i) = ej2(1,i) else call geopro(zma1(1,i)) call geopro(zma2(1,i)) call geopro(zma3(1,i)) call geopro(zma4(1,i)) call geopro(ej1(1,i)) call geopro(ej2(1,i)) endif 300 continue c c*****put adjoint matrix for current solar zenith angle into b c b(1)=admatx(1,jmuz) b(2)=admatx(2,jmuz) b(3)=admatx(3,jmuz) b(4)=0.0d0 c do 105 im=1,imu it=4 gr(it)=0.0d0 gl(it)=0.0d0 ztp1(it)=0.0d0 ztp2(it)=0.0d0 do 107 jd=1,4 107 ztop(it,jd)=0.0d0 c c*****do integration over optical thickness c do 110 i=1,ncp1 tmu=extmu(i,im) ztop(it,1)=ztop(it,1)+zma1(it,i)*tmu ztop(it,2)=ztop(it,2)+zma2(it,i)*tmu ztop(it,3)=ztop(it,3)+zma3(it,i)*tmu ztop(it,4)=ztop(it,4)+zma4(it,i)*tmu ztp1(it)=ztp1(it)+ej1(it,i)*tmu ztp2(it)=ztp2(it)+ej2(it,i)*tmu gr(it)=gr(it)+ek4(i)*zma1(it,i)+ek5(i)*zma3(it,i) gl(it)=gl(it)+ek4(i)*zma2(it,i)+ek5(i)*zma4(it,i) 110 continue c c write(33,132)ipol, ztp1(it),ztp2(it) 132 format('debug evalit..ipol,ztp1,ztp2',i2,1x,1p2e12.4) c c*****calculate the azimuth -independent intensity eizero at current c*****value of the polar angle. see eqs(3.7) and (4.16) c c*****multiply matrix top(1-4) times the adjoint matrix admatx(=b). cal c*****call the product c c c c*****multiply c by ematx c a(1)=ematx(1,im) a(2)=ematx(2,im) a(3)=ematx(3,im) a(4)=0.0d0 c c c do 120 i=1,4 c c**** store extrapolated z-matrix in hold(4) c hold(i)=ztop(it,i) 120 continue call matmul(hold,b,c) call matmul(a,c,d) xhold(4)=d(1)+d(2)+d(3)+d(4) 130 continue c new statements added here if(ipol.eq.0) then eizero(im)=0.09375d0*xhold(4)/emu(im) else eizero(im)=0.125d0*(q1/sdp)*xhold(4)/emu(im) endif c e0za(itz,im)=eizero(im) c if(im.eq.1 .and. jprint(4) .eq. 1)then write(33,515)itz,im,emu(im),eizero(im),e0za(itz,im) 515 format('sub eva1pol.itz,im,emu,eizero,e0za',2i5,1x,3f8.5) endif c c*****now calculate azimuth dependent intensities eiaz1 and eiaz2 at all c*****azimuth angles (1 to iazmth) for current value of polar angle. c***** see eqs.(3.7) and (4.8) c if(ipol.eq.0)then as=ematx(1,im)-1.0d0 bs=admatx(1,jmuz)-1.0d0 else as=emu(im)*emu(im)-1.0d0 bs=emuz(jmuz)*emuz(jmuz)-1.0d0 endif c atb=as*bs absq=dsqrt(atb)*emuz(jmuz) atmu=as*emuz(jmuz) c do 225 jaz=1,iazmth c new statements added here if(ipol.eq.0)then eiaz1(im,jaz)=-0.375d0*absq*caz(jaz)*ztp1(4) eiaz2(im,jaz)=0.09375d0*(atb*caz2(jaz))*ztp2(4)/emu(im) else eiaz1(im,jaz)=-0.375d0*q2*absq*caz(jaz)*ztp1(4) eiaz2(im,jaz)=0.09375d0*q2*(atb*caz2(jaz))*ztp2(4)/emu(im) endif c e1za(itz,im,jaz)=eiaz1(im,jaz) e2za(itz,im,jaz)=eiaz2(im,jaz) c c*****compute total intensity at top of atmosphere c totint(im,jaz)=eizero(im)+eiaz1(im,jaz)+eiaz2(im,jaz) c if (jprint(4) .eq. 1) * write(33,133)itz,im,jaz,eizero(im),eiaz1(im,jaz), * eiaz2(im,jaz),totint(im,jaz) 133 format('eva1.itz,im,jaz,ez0,ez1,ez2,tot',3i4,4f9.4) 225 continue 105 continue c c*****compute the downward diffuse intensity gg c if(ipol .eq. 0)then c gg=c316*(gr(4)*(admatx(1,jmuz)+1.0d0)+gl(4)*admatx(3,jmuz)) f1=c316*(admatx(1,jmuz)+1.0d0) f2=c316*admatx(3,jmuz) gg=f1*gr(4)+f2*gl(4) else f1=0.25d0*q1*(admatx(1,jmuz)+1.0d0+q) f2=0.25d0*q1*admatx(3,jmuz) gg=(f1*gr(4)+f2*gl(4))/sdp endif c ggz(itz)=gg c for debugging if (jprint(4) .eq. 1) write(33,212)gl(4),gr(4),f1,f2,gg if (jprint(4) .eq. 1) write(33,312)(e0za(itz,ij),ij=1,imu) if (jprint(4) .eq. 1) write(33,313)ggz(itz) 312 format('e0za',1p6e11.3) 313 format('ggz',1pe11.3) 212 format('eva1pol...gl,gr,f1,f2,gg',5f8.4) c return end subroutine eva2pol(zst1,zst2,ncp1,itz) c c********************************************************************** c c subroutine eva2pol c c this routine is adopted from subroutine evalrf to extract sb,tnstr c after each iteration......zia ahmad 10/21/94 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 c implicit integer*4(i-n),real*8(a-h,o-z) real *8 tenstr(4,15),zst1(4,202),zst2(4,202) real*8 thnot(10),thta(15),azmth(8),hhold(101),alb(11) real*8 x(101) integer jprint(10) common /emm/ ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) common /eks/ extmu(202,15),ek4(202),ek5(202) common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /out/ sb(4),qr(11),tnstr(15),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),gg,fs,fsp(15) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /prints/ jprint c c -- common block in added on 11/19/92 common /in/ delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c common block depolt added on 9/22/93 common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c common/buff1/sbz(9),tnstrz(9,15),tnstrl(15),tnstrr(15) data c316/0.1875d0/ c c write(23,500)rhon,gama,q,q1,q2,delp,sdp,q12s,ipol 500 format('rhon,gama,q,q1,q2,delp,sdp,q12s,ipol'/ 1 1p6e12.4/1p2e12.4,i5) c compute gometrical progression of zstar layer by layer c itmaxp=itmax+1 do i=1,ncp1 if (itz .le. itmaxp) then zst1(4,i) = zst1(1,i) zst2(4,i) = zst2(1,i) else call geopro(zst1(1,i)) call geopro(zst2(1,i)) endif enddo c it=4 c c find sb(it) sum2=0.d0 c sum2l=0.0d0 sum2r=0.0d0 do 110 i=1,ncp1 sum2=sum2+ek4(i)*zst1(it,i)+ek5(i)*zst2(it,i) 110 continue c c new statements added here if(ipol.eq.0)then sb(it)=0.375d0*sum2 else sb(it)=q12s*sum2 endif c sbz(itz)=sb(it) if (jprint(4) .eq. 1) write(33,501)itz,itmax,ipol,sb(it) 501 format('eva2pol...itz,itmax,ipol,sb',3i5,1pe12.4) c find tnstr(it) c do 120 j=1,imu suma=0. sumb=0. c do 130 i=1,ncp1 suma=suma+zst1(it,i)*extmu(i,j) sumb=sumb+zst2(it,i)*extmu(i,j) 130 continue c c multiply by m-matrix see eq(6.7) c tenr=ematx(1,j)*suma+ematx(2,j)*sumb tenl=ematx(3,j)*suma c tensum=tenl+tenr c write(33,565)suma,sumb,ematx(1,j),ematx(2,j),ematx(3,j),emu(j), c 1 tenl,tenr,tensum 565 format('evalrf',9f8.5) c c compute istar/ig of eq (6.7) c c new statements added here if(ipol.eq.0)then tnstr(j)=c316*(tenr+tenl)/emu(j) tnstrl(j)=c316*tenl/emu(j) tnstrr(j)=c316*tenr/emu(j) else tnstr(j)=0.25d0*(q1/sdp)*(tenr+tenl)/emu(j) tnstrl(j)=0.25d0*(q1/sdp)*tenl/emu(j) tnstrr(j)=0.25d0*(q1/sdp)*tenr/emu(j) endif tnstrz(itz,j)=tnstr(j) 120 continue c return end subroutine evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1) c*********************************************************************** c c subroutine evalit c c purpose c evalit is used to perform the integration of the reduced c source functions over the total optical thickness to each c level in the model atmosphere. c evalit is also used to transfer solar angle dependent data c to the archive and tables tapes. c c method c the iterated reduced source functions as calculated by itrate c are integrated over the total cumulative optical thickness.the c resultant integrated quantities (ztop,ztp1,ztp2,gg) are extrapol c extrapolated using a geometric progression. finally certain c matrix multiplications are carried out to transform the c integrated quantities into intensity terms. c c calling sequence c call evalit(zma1,zma2,zma3,zma4,ej1,ej2,jmuz,ncp1) c c variable type i/o description c -------- ---- --- ----------- c c zma1(4,202) r*8 i first element of z matrix c for each source function c level for final three iters. c fourth row of this matrix is c reserved for extrapolated values. c zma2,3,4 r*8 i/o second,third,fourth elements c of z matrix. c ej1(4,202) r*8 o source functions for cos(phi) c dependent terms for last three c iterations of itrate. c ej2(4,202) r*8 i/o source functions for cos(2*phi) c dependent terms. c jmuz i*4 i current index counter for zenith c angle. c ncp1 i*4 i number of levels in model atmos. c common /out/ contains intensity terms computed by c evalrf,evalit and intsum routines. c c common /out/ c c variables type i/o description c --------- ---- --- ----------- c c sb(4) r*8 i fractional atmos. back scatter c factor.sb(4) is extrapolated value. c tnstr(15) r*8 i extrapolated values of istar/ig for c each polar look angle. c qr(10) r*8 i combined reflectivity -sb term for c all albedo values. c eizero(15) r*8 o azimuth independent intensity at each c polar look angle. c eiaz1(15,8) r*8 o cos(phi) azimuth dependent term for c each polar look angle and scan plane c angle. c eiaz2(15,8) r*8 o cos(2*phi) dependent term for each c polar look angle and scan plane angle. c gg r*8 o integrated downward diffuse flux at c ground for current zenith angle. c fs r*8 i direct solar flux at ground for curren c current solar zenith angle. c totint(15,8) r*8 o total intensity at top of atmosphere c for each polar and azimuth angle. c c analysis and programming k. f. klenk, p. m. smith sasc aug. 77 c c modifications (date name purpose) c c last modified 11/18/92 by zia ahmad c modified for single iteration c c last modified 9/22/93 ... zia ahmad c purpose: to add the effect of molecular depolarization c c last modified 10/12/94... zia ahmad c purpose: modified to compute and print polarization of the c diffused radiation. c c last modified 03/08/95... dave flittner c purpose: store Z1func=I1/(-3/8*muo*sqrt(1-muo^2)*sqrt(1-mu^2)) and c Z2func=I2/(3/32*(1-muo^2)*(1-mu^2)/mu) so can be used as input into c interpolation tables. 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 implicit integer*4 (i-n), real*8(a-h,o-z) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101),alb(11) real*8 x(101) integer jprint(10) real*8 ztop(4,4),a(4),b(4),c(4),d(4),ztp1(4),ztp2(4),top(4), & dum(8),gr(4),gl(4), & zma1(4,202),zma2(4,202),zma3(4,202),zma4(4,202), & ej1(4,202),ej2(4,202),hold(4),xhold(4),eiold(15) real*4 out48(48) dimension raray(1214) common /tabbuf/out48 c common /hedout/ rarray(824),psaray(82),iarray(3) common /out/ sb(4),qr(11),tnstr(15),eizero(15),eiaz1(15,8), & eiaz2(15,8),gg,fs,totint(15,8) common /emm/ ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) c common /eks/ extmu(202,15),ek4(202),ek5(202) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip c common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, c 1 x,pshold,hhold,layer,thnot,alfaef, c 2 lambda common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /prints/ jprint common /in/delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of statement c common/buff3/eitl(15,8),eitr(15,8),eitu(15,8), ttl(15), ttr(15) c real*8 eiaz1l(15,8),eiaz1r(15,8),eiaz1u(15,8), 1 eiaz2l(15,8),eiaz2r(15,8),eiaz2u(15,8),pol(15,8), 2 eizrol(15),eizror(15) equivalence (rarray(1),raray(1)) c common/czfunc/Z1func(15),Z2func(15) !def c data c316/0.1875d0/ c extrapolate source functions layer by layer c do 300 i=1,ncp1 if (itmax.eq.1) then zma1(4,i) = zma1(1,i) zma2(4,i) = zma2(1,i) zma3(4,i) = zma3(1,i) zma4(4,i) = zma4(1,i) ej1(4,i) = ej1(1,i) ej2(4,i) = ej2(1,i) else call geopro(zma1(1,i)) call geopro(zma2(1,i)) call geopro(zma3(1,i)) call geopro(zma4(1,i)) call geopro(ej1(1,i)) call geopro(ej2(1,i)) endif 300 continue c c*****put adjoint matrix for current solar zenith angle into b c b(1)=admatx(1,jmuz) b(2)=admatx(2,jmuz) b(3)=admatx(3,jmuz) b(4)=0.0d0 c if(jprint(1).eq.1) write(33,1000) do 105 im=1,imu it=4 gr(it)=0.0d0 gl(it)=0.0d0 ztp1(it)=0.0d0 ztp2(it)=0.0d0 do 107 jd=1,4 107 ztop(it,jd)=0.0d0 c c*****do integration over optical thickness c do 110 i=1,ncp1 tmu=extmu(i,im) ztop(it,1)=ztop(it,1)+zma1(it,i)*tmu ztop(it,2)=ztop(it,2)+zma2(it,i)*tmu ztop(it,3)=ztop(it,3)+zma3(it,i)*tmu ztop(it,4)=ztop(it,4)+zma4(it,i)*tmu ztp1(it)=ztp1(it)+ej1(it,i)*tmu ztp2(it)=ztp2(it)+ej2(it,i)*tmu gr(it)=gr(it)+ek4(i)*zma1(it,i)+ek5(i)*zma3(it,i) gl(it)=gl(it)+ek4(i)*zma2(it,i)+ek5(i)*zma4(it,i) 110 continue if(ipol.eq.0)then !def Z1func(im)=ztp1(4) !def Z2func(im)=ztp2(4) !def else !def Z1func(im)=q2*ztp1(4) !def Z2func(im)=q2*ztp2(4) !def endif !def 100 continue c c write(33,132)ipol, ztp1(it),ztp2(it) 132 format('debug evalit..ipol,ztp1,ztp2',i2,1x,1p2e12.4) c c*****calculate the azimuth -independent intensity eizero at current c*****value of the polar angle. see eqs(3.7) and (4.16) c c*****multiply matrix top(1-4) times the adjoint matrix admatx(=b). cal c*****call the product c c c c*****multiply c by ematx c a(1)=ematx(1,im) a(2)=ematx(2,im) a(3)=ematx(3,im) a(4)=0.0d0 c c c do 120 i=1,4 c c**** store extrapolated z-matrix in hold(4) c hold(i)=ztop(it,i) 120 continue call matmul(hold,b,c) call matmul(a,c,d) xhold(4)=d(1)+d(2)+d(3)+d(4) 130 continue c new statements added here if(ipol.eq.0) then eizero(im)=0.09375d0*xhold(4)/emu(im) eizrol(im)=0.09375d0*(d(1)+d(2))/emu(im) eizror(im)=0.09375d0*(d(3)+d(4))/emu(im) else eizero(im)=0.125d0*(q1/sdp)*xhold(4)/emu(im) eizrol(im)=0.125d0*(q1/sdp)*(d(1)+d(2))/emu(im) eizror(im)=0.125d0*(q1/sdp)*(d(3)+d(4))/emu(im) endif c c if(im.eq.1)then if (jprint(1) .eq. 1) write(33,514)im,emu(im), * d(1),d(2),d(3),d(4) 514 format('im,emu,d1,d2,d3,d4',i3,1x,5f8.4) if (jprint(1) .eq. 1) write(33,515)im,emu(im), * eizero(im),eizrol(im),eizror(im) 515 format('evalitpol..im,emu,ei0,eil,eir',1x,i2,1x,4f8.5) c endif c c*****now calculate azimuth dependent intensities eiaz1 and eiaz2 at all c*****azimuth angles (1 to iazmth) for current value of polar angle. c***** see eqs.(3.7) and (4.8) c if(ipol.eq.0)then as=ematx(1,im)-1.0d0 bs=admatx(1,jmuz)-1.0d0 else as=emu(im)*emu(im)-1.0d0 bs=emuz(jmuz)*emuz(jmuz)-1.0d0 endif c cs=emu(im)*emu(im) atb=as*bs absq=dsqrt(atb)*emuz(jmuz) atmu=as*emuz(jmuz) c c write(33,134) ipol,im,jmuz,as,bs,atb,absq,atmu 134 format('ipol,as,bs,atb,absq,atmu'/3i2,1x,1p5e12.4) c do 225 jaz=1,iazmth c new statements added here if(ipol.eq.0)then eiaz1(im,jaz)=-0.375d0*absq*caz(jaz)*ztp1(4) eiaz2(im,jaz)=0.09375d0*(atb*caz2(jaz))*ztp2(4)/emu(im) c new code for polarization eiaz1l(im,jaz)=eiaz1(im,jaz) eiaz1r(im,jaz)=0.0d0 eiaz1u(im,jaz)=0.375d0*absq*dsqrt(1.0d0-caz(jaz)**2)* 1 ztp1(4)/emu(im) eiaz2l(im,jaz)=0.09375d0*cs*bs*caz2(jaz)*ztp2(4)/emu(im) eiaz2r(im,jaz)=-0.09375d0*bs*caz2(jaz)*ztp2(4)/emu(im) eiaz2u(im,jaz)=-0.09375d0*bs* 1 dsqrt(1.0d0-caz2(jaz)**2)*ztp2(4) else eiaz1(im,jaz)=-0.375d0*q2*absq*caz(jaz)*ztp1(4) eiaz2(im,jaz)=0.09375d0*q2*(atb*caz2(jaz))*ztp2(4)/emu(im) c new code for polarization eiaz1l(im,jaz)=eiaz1(im,jaz) eiaz1r(im,jaz)=0.0d0 eiaz1u(im,jaz)=0.375d0*q2*absq*dsqrt(1.0d0-caz(jaz)**2)* 1 ztp1(4)/emu(im) eiaz2l(im,jaz)=0.09375d0*q2*cs*bs*caz2(jaz)*ztp2(4)/emu(im) eiaz2r(im,jaz)=-0.09375d0*q2*bs*caz2(jaz)*ztp2(4)/emu(im) eiaz2u(im,jaz)=-0.09375d0*q2*bs* 1 dsqrt(1.0d0-caz2(jaz)**2)*ztp2(4) endif c c*****compute total intensity at top of atmosphere c totint(im,jaz)=eizero(im)+eiaz1(im,jaz)+eiaz2(im,jaz) c c compute degree of polarization eitl(im,jaz)=eizrol(im)+eiaz1l(im,jaz)+eiaz2l(im,jaz) eitr(im,jaz)=eizror(im)+eiaz1r(im,jaz)+eiaz2r(im,jaz) eitu(im,jaz)=eiaz1u(im,jaz)+eiaz2u(im,jaz) pol(im,jaz)=dsqrt((eitl(im,jaz)-eitr(im,jaz))**2 1 +eitu(im,jaz)**2)/totint(im,jaz) c if (jprint(1) .eq. 1) write(33,133)im,jaz,eizero(im), * eiaz1(im,jaz),eiaz2(im,jaz),totint(im,jaz),pol(im,jaz) 133 format('evalitpol...im,jaz,ei0,ei1,ei2,totint,pol'/ 1 2i4,1p5e11.3) if (jprint(1) .eq. 1) write(33,136)eitl(im,jaz),eitr(im,jaz), * eitu(im,jaz) 136 format(8x,1p3e11.3) 225 continue if(jprint(1).eq.1) write(33,1100)im,thta(im), 1 (ztop(4,j),j=1,4),ztp1(4),ztp2(4) 105 continue if(jprint(1).eq.1) write(33,2000) gl(4),gr(4) c c*****compute the downward diffuse intensity gg c if(ipol .eq. 0)then c gg=c316*(gr(4)*(admatx(1,jmuz)+1.0d0)+gl(4)*admatx(3,jmuz)) f1=c316*(admatx(1,jmuz)+1.0d0) f2=c316*admatx(3,jmuz) gg=f1*gr(4)+f2*gl(4) else f1=0.25d0*q1*(admatx(1,jmuz)+1.0d0+q) f2=0.25d0*q1*admatx(3,jmuz) gg=(f1*gr(4)+f2*gl(4))/sdp endif c c for debugging c write(33,212)gl(4),gr(4),f1,f2,gg,q 212 format('gl,gr,f1,f2,gg,q',6f8.4) c write(33,213)(admatx(iz,jmuz),iz=1,3) 213 format(1x,'admatx 1-3',3f8.4) c c**** store data for archive tape raray(1)=fs raray(2)=gg k1=3 do 400 i=1,ncp1 k5=810+i k6=1012+i raray(k1)=zma1(4,i) raray(k1+1)=zma2(4,i) raray(k1+2)=zma3(4,i) raray(k1+3)=zma4(4,i) raray(k5)=ej1(4,i) raray(k6)=ej2(4,i) k1=k1+4 400 continue C write (34) raray,iarray c out48(1)=thnot(jmuz) c j=2 c do 420 i=1,15 c out48(j)=eizero(i) c out48(j+15)=eiaz1(i,1) c out48(j+30)=eiaz2(i,1) c j=j+1 c 420 continue c out48(47)=fs c out48(48)=gg c write(35)out48 500 continue c**** call intsum to compute total intensity call intsum if (jprint(10).eq.1) call sumry(jmuz) return 1000 format(1h1,t5,'debug print out from evalit',10x, 1 'source functions(extrapolated) integrated over', 2 1x,'optical thickness',///, 2 1x,'imu',2x,'theta',5x,'ztop1',9x,'ztop2',9x,'ztop3',9x, 3 'ztop4',9x,' ztp1',9x,' ztp2',//) 1100 format(1x,i3,2x,f5.2,2x,d12.6,2x,d12.6,2x,d12.6,2x,d12.6, 1 2x,d12.6,2x,d12.6) 2000 format(///,1x,'gl= ',d12.6,/,1x,'gr= ',d12.6) end subroutine evalrf(zst1,zst2,ncp1) c c********************************************************************** c cccc c subroutine evalrf c c version aug 22,1977 c c purpose c c evalrf is a fortran iv routine which is used to perfform c the integration of eq(6.7) to find tnstr (istar/ig) c c method c c using a trapezoidal integration , evalrf computes tenstr which c represents the sum of all intensities upto and including the c intensity due to itmax+1 order scattering of surface reflected c radiation. c the remaining higher order terms are calulated by gepro. c c calling sequence c c call evalrf (zst1,zst2,tnstr,qr,sb) c c variable type i/o description c -------- ---- --- ----------- c c zst1(4,202) r*8 i last three iterations of c zstar matrix and extrapolated value c zst2(4,202) r** i similar to zst1 c tnstr(4,15) r*8 o integrated intensities from last c three iterations of z-star matrix c and extrapolated intensity. c contains results for each polar c look angle c qr(11) r*8 o factor used in computing ig c sb(4) r*8 o fraction of surface recflected c radiation reaching top of atmosphere c c external references c geopro c c common areas referenced c c c analysis and programming 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 last modified 9/22/93 .... zia ahmad c purpose: to add the effect of molecular anisotropy 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. cccc c*********************************************************************** c c implicit integer*4(i-n),real*8(a-h,o-z) real *8 tenstr(4,15),zst1(4,202),zst2(4,202) real*8 thnot(10),thta(15),azmth(8),hhold(101),alb(11) real*8 x(101) integer jprint(10) real*4 out48(48) common /hedout/ rarray(824),psaray(82),iarray(3) common /tabbuf/ out48 common /emm/ ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) common /eks/ extmu(202,15),ek4(202),ek5(202) c c common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, c 1 x,pshold,hhold,layer,thnot,alfaef, c 2 lambda common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /out/ sb(4),qr(11),tnstr(15),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),gg,fs,fsp(15) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /prints/ jprint c c -- common block in added on 11/19/92 common /in/ delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c common block depolt added on 9/22/93 common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c data c316/0.1875d0/ c c write(23,500)rhon,gama,q,q1,q2,delp,sdp,q12s,ipol 500 format('rhon,gama,q,q1,q2,delp,sdp,q12s,ipol'/ 1 1p6e12.4/1p2e12.4,i5) c compute gometrical progression of zstar layer by layer c do 160 i=1,ncp1 if (itmax.eq.1) then zst1(4,i) = zst1(1,i) zst2(4,i) = zst2(1,i) else call geopro(zst1(1,i)) call geopro(zst2(1,i)) endif 160 continue c it=4 c c find sb(it) sum2=0.d0 c sum2l=0.0d0 sum2r=0.0d0 do 110 i=1,ncp1 sum2=sum2+ek4(i)*zst1(it,i)+ek5(i)*zst2(it,i) 110 continue c c new statements added here if(ipol.eq.0)then sb(it)=0.375d0*sum2 else sb(it)=q12s*sum2 endif c if (jprint(10).eq.1) write(23,501)sb,ipol if (jprint(2) .eq. 1) write(33,501)sb,ipol 501 format('evalrf..sb,ipol',1p4e12.4,i5) c find tnstr(it) c do 120 j=1,imu suma=0. sumb=0. c do 130 i=1,ncp1 suma=suma+zst1(it,i)*extmu(i,j) sumb=sumb+zst2(it,i)*extmu(i,j) 130 continue c c multiply by m-matrix see eq(6.7) c tenr=ematx(1,j)*suma+ematx(2,j)*sumb tenl=ematx(3,j)*suma c tensum=tenl+tenr c write(33,565)suma,sumb,ematx(1,j),ematx(2,j),ematx(3,j),emu(j), c 1 tenl,tenr,tensum 565 format('evalrf',9f8.5) c c compute istar/ig of eq (6.7) c c new statements added here if(ipol.eq.0)then tnstr(j)=c316*(tenr+tenl)/emu(j) else tnstr(j)=0.25d0*(q1/sdp)*(tenr+tenl)/emu(j) endif 120 continue c 100 continue c c do 145 j=1,nalb qr(j)=alb(j)/(1.-sb(4)*alb(j)) 145 continue c c c**** store data for archive tape c**** check if tables data set to be created c rarray(407)=sb(4) j=420 do 200 i=1,ncp1 k=j+1 rarray(j)=zst1(4,i) rarray(k)=zst2(4,i) j=j+2 200 continue c**** write data to tape c c write (34) rarray,psaray,iarray c c**** check if tables data set to be created c c**** load r*4 output buffer for later fwrite c out48(1)=float(lambda) c out48(2)=alfaef c out48(3)=beta c out48(4)=code c out48(5)=pnot c out48(6)=sb(4) c j=7 c do 275 iload=1,15 c out48(j)=tnstr(iload) c j=j+1 c 275 continue c out48(22)=float(imuz) c out48(23)=float(imu) c j=24 c do 260 i=1,15 c out48(j)=thta(i) c j=j+1 c 260 continue c write(35) out48 c c*****optional print out of z-star array c if (jprint(2) .eq. 1) then write(33,1000) write(33,2000) (k,(zst1(j,k),j=1,4),(zst2(j,k), j=1,4), 1 k=1,ncp1) write(33,2500) write(33,3000) (k,thta(k),tnstr(k),k=1,imu) write(33,3500) sb(it) endif c return 1000 format(1h1,t50,'debug printout for evalrf',///,25x, 1 'zst1(i=1,4)',47x,'zst2(i=1,4)') 2000 format(2x,i3,4e13.5,5x,4e13.5) 2500 format(///,2x,'imu',20x,'tnstr(j)') 3000 format(5x,'theta.',5x,'tnstr(k)',//,(1x,i2,3x,f6.2,2x,d12.6)) 3500 format(///,' sb=',e15.5) end subroutine expone (nc,ncp1,imu) c c*********************************************************************** cccc c subroutine expone c version aug 18,1977 c c purpose c c expone is a fortran 4 routine to calulate the iteration c integrals of eq 4.11 and k functions of eq 4.6 c c method c kernals of intgrals to be used in reflex and iterate c are calulated and stored on a disk file c c calling sequence c call expone (ncp,ncp1) c c variable type i/o description c -------- ---- --- ------------ c c nc i*4 i # of layers from top of c atmosphere to reflecting surface c ncp1 i*4 o nc+1 c c external references c dexpk1 c c common areas referenced c consts c eks c emm c es c prints c thkns c c analysis and programming k.f. klenk sasc aug 18,1977 c c modifications (date , name, purpose) 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 cccc c*********************************************************************** c implicit integer*4(i-n),real*8(a-h,o-z) real*8 ekary(202,5),extmux(202,15) real*8 nek(202,5,202) integer jprint(10) common /prints/ jprint common /kmat/nek c c common block area c common /emm/ ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) common /eks/ extmu(202,15),ek4(202),ek5(202) real *8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) common /thkns/ tsl,tt,tts,ttl,dtsp,dtts common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /es/ odan(6),e(6),eek(6) c c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c do 75 i=1,ncp1 c initialize ekary do 200 j=1,5 ekary(i,j)=0.d0 200 continue c c calculate w(t)dt*attenuation factor (see eqn. 3.6,6.7) c do 205 j=1,imu extmu(i,j)=dtsp(i)*dexp(-tt(i)/emu(j)) extmux(i,j)=extmu(i,j)/emu(j) 205 continue if (i.eq.1) go to 55 im1=i-1 c do 210 j=1,5 ekary(im1,j)=0. 210 continue c if(i.eq.2) go to 50 im2=i-2 c c calculates the elements of k1 matrix (see eqn 4.11) c evaluate k1(1),k1(2) of eqn 4.7 for current value of c argument (tt(i)-tt(j)) c do 45 j=1,im2 call dexpk(tt(i),tt(j)) c c ekary array contains k matrix*w(t)dt c do 220 jk=1,5 ekary(j,jk)=dtsp(j)*eek(jk) 220 continue if(i.lt.ncp1) go to 45 c new statements if(ipol.eq.0)then ek4(j)=(e(2)+e(4))*dtsp(j) ek5(j)=sq2*(e(2)-e(4))*dtsp(j) else ek4(j)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(j) ek5(j)=delp*(e(2)-e(4))*dtsp(j) endif c end of new statements 45 continue c call dexpk(tt(i),tt(im1)) do 230 k=1,5 ekary(im1,k)=dtts(im1)*eek(k) 230 continue c if(i.lt.ncp1) go to 50 c new statements if(ipol.eq.0)then ek4(nc)=dtsp(nc)*(e(2)+e(4)) ek5(nc)=sq2*dtsp(nc)*(e(2)-e(4)) else ek4(nc)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(nc) ek5(nc)=delp*(e(2)-e(4))*dtsp(nc) endif c end of new statements 50 call dexpk1(tt(i),tt(im1)) c do 240 k=1,5 ekary(im1,k)=ekary(im1,k)+dtts(i)*eek(k) 240 continue c do 250 k=1,5 ekary(i,k)=dtts(i)*eek(k) 250 continue if(i.lt.ncp1) go to 55 c new statemnets if(ipol.eq.0)then ek4(ncp1)=dtsp(ncp1)*(odan(2)+odan(4)) ek5(ncp1)=sq2*dtsp(ncp1)*(odan(2)-odan(4)) else ek4(ncp1)=((1.0d0+2.0d0*q)*odan(2)+odan(4))*dtsp(ncp1) ek5(ncp1)=delp*(odan(2)-odan(4))*dtsp(ncp1) endif c end of new statements go to 65 55 ip1=i+1 call dexpk1(tt(i),tt(ip1)) c do 260 k=1,5 ekary(i,k)=ekary(i,k)+dtts(ip1)*eek(k) 260 continue c do 270 k=1,5 ekary(ip1,k)=dtts(ip1)*eek(k) 270 continue c if(i.eq.nc) go to 65 ip2=i+2 call dexpk(tt(i),tt(ip1)) c do 280 k=1,5 ekary(ip1,k)=ekary(ip1,k)+dtts(ip2)*eek(k) 280 continue c do 60 j=ip2,ncp1 call dexpk(tt(i),tt(j)) c do 290 k=1,5 ekary(j,k)=dtsp(j)*eek(k) 290 continue 60 continue c c 65 continue do 5000 ii = 1,202 do 5001 iii = 1,5 nek(ii,iii,i) = ekary(ii,iii) 5001 continue 5000 continue 75 continue c c optionally print calculated exponential terms c (zstar matrix * w(t)dt see eqn 6.3) c if (jprint(3).ne.0) then write (33,6120) (j,ek4(j),ek5(j) ,j=1,ncp1) do 295 i=1,imu thta=dacos(emu(i))/cnvrt write(33,7120) thta,(j,extmux(j,i),j=1,ncp1,10) write(33,7130) ncp1,extmux(ncp1,i) 295 continue endif return 6120 format(1h1,t50,'debug print out for expone',///, 1 1h ,2(' j',4x,'ek4(j)',6x,'ek5(j)',6x),//,2(1x,i3,2d12.4,3x)) 7120 format(1x,//,1x,'theta=',f6.2,' *** attenuation factors/*** ',//, 1 1x,'level',5x,'extmu/mu',/,(1x,i5,3x,d12.6)) 7130 format(1x,i5,3x,d12.6) end subroutine exponesph(h,ps,cofx,lmax,nc,ncp1,imu) c c*********************************************************************** cccc c subroutine exponesph c c purpose- c 1. calculate the slant path scattering and total optical depths c for a particular line of sight. This is used in the c integration of the output beam. c 2. Then calculate the trapezoidal integration factors using the c calculated slant path optical depths. c c 3. calculate the iteration integration factors for a flat c atmosphere for eq. 4.11 and k functions of eq 4.6. c This is done the same way as in routine expone. c c method- c 1. using the fine layer scattering and total optical depths in c ps and dxs, and the altitude of the layer boundaries, h, c use simple spherical trig. to calc. the ratio of the distance c through a layer in a spherical atmosphere to that in a flat c atmosphere. The result is stored in the arrays ttsf_sph and c ttf_sph. c 2. The integration factors are computed using extmu=layer scattering c optical thickness * transmission to from level to top of atmo. c Integral = 0.5*Source(1)*Transmission(1)* c scattering optical thickness of layer 1 + c Sum from i=2 to nc of 0.5*(Source(i)*Transmission(i) c +Source(i+1)*Transmission(i+1))* c scattering optical thickness of layer i + c 0.5*Source(ncp1)*Transmission(ncp1) c scattering optical thickness of layer ncp1 c 3. If the scan angle is equal to 0.0 degrees, then the quantities c are computed quickly using a flat atmosphere. Again, note c that the integration factors used in the iteration process c are for a flat atmosphere. Only the final integration of the c out-going beam is done in a spherical atmosphere. c c calling sequence- c call exponesph(h,ps,cofx,lmax,nc,ncp1,imu) c c variable type i/o description c -------- ---- --- ----------- c c h(487) r*8 i height from earth center(fraction of r) c ps(487) r*8 i rayleigh opt. thickness of each layer c of the standard atmosphere c cofx(4,487) r*8 i spline interpolation coeff. c lmax i*4 i number of layers in the standard atmosphere c nc i*4 i # of layers from top of c atmosphere to reflecting surface c ncp1 i*4 i nc+1 c imu i*4 i # of scan angles c passed through common blocks c dxs(487) r*8 i total opt. thickness of each layer c of the standard atmosphere c xs(487) r*8 i log of the total opt. depth to a level c of the standard atmosphere c ttl(202) r*8 i log of the total vertical optical depth c emu(15) r*8 i cosine of the scan angle c extmu(202,15) r*8 o trapezoidal integration factor for each level c c internal arrays of note c dtsp_sph(202) r*8 avg. scattering optical thickness of each c spherical layer c dtts_sph(202) r*8 avg. scattering optical thickness of each c spherical half layer at mid points between c layers c transph(202) r*8 transmission from level to the top of the c spherical atmosphere c c external references c dexpk1 (subroutine) c omerf (function) c c common areas referenced c consts c eks c emm c es c prints c thkns c chpmn c csphout c kmat c crefdir c depolt 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*********************************************************************** cccc implicit none c input integer*4 lmax,nc,ncp1,imu real*8 h(487),ps(487),emu(15),cofx(4,487) c input through common blocks real*8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) real*8 chpn,chxn,sqchp,sqchx,xs(487),dxs(487) logical lsphout c passed through common blocks integer jprint(10),numek,kskip,ipol real*8 ekary(202,5),extmux(202,15) real*8 nek(202,5,202) real*8 ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emuz(10) real*8 ek4(202),ek5(202) real*8 pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg real*8 odan(6),e(6),eek(6) real*8 rhon,gama,q,q1,q2,delp,sdp,q12s c output through common block real*8 extmu(202,15),refdir(15) c internal integer*4 i,j,lmaxm1,im,im1,im2,jk,k,ip1,ip2,ii,iii real*8 sintheta,pc,pcpc,dist1,dist2,ddist,ratio,dum1,dum2,dum3 real*8 ttsf_sph(487),ttf_sph(487) real*8 tts_sph(202) real*8 dtts_sph(202),dtsp_sph(202),transph(202) real*8 chpp,chxx,thta,sum,thetap c c common block area c common /thkns/tsl,tt,tts,ttl,dtsp,dtts common /chpmn/chpn,chxn,sqchp,sqchx,xs,dxs common /csphout/lsphout common /prints/ jprint common /kmat/nek common /emm/ ematx,admatx,saz,caz,saz2,caz2,emu,emuz common /eks/ extmu,ek4,ek5 common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /es/ odan,e,eek common /crefdir/refdir c c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c external functions used real*8 omerf external omerf c end of portion lifted from expone ! def c c loop for each scan angle c do 199 im=1,imu c now loop using the simple trig. relations c calc. the tangent height in terms of earth radii using the scan angle c at the surface sintheta=(1.0d0-emu(im)*emu(im)) if(sintheta.le.0.0d0)then c do not need to correct, for sphericity, so will use the flat atmo. case do i=1,ncp1 extmu(i,im)=dtsp(i)*dexp(-tt(i)) enddo refdir(im)=dexp(-tt(ncp1)) goto 199 else if(sintheta.gt.1.0d0)then sintheta=1.0d0 pc=h(lmax) else sintheta=dsqrt(sintheta) pc=h(lmax)*sintheta endif c start at top of atmosphere and work way down along the line of sight c use chapman function to do the first layer which ranges from h(1) to c infinity. Note that the approximation used to evaluate the Chapman c function for dum1=1.0 (scan angle of 0.0 deg.) returns a value less c than that in a flat atmosphere. dum1=pc/h(1) dum1=dsqrt(1.0d0-dum1**2) chpp=omerf(sqchp*dum1,dum1,chpn) chxx=omerf(sqchx*dum1,dum1,chxn) ttsf_sph(1)=ps(1)*chpp ttf_sph(1)=(dxs(1)-ps(1))*chxx+ps(1)*chpp c now start the loop to create the fine arrays pcpc=pc*pc dist1=dsqrt(h(1)*h(1)-pcpc) do i=2,lmax dist2=dsqrt(h(i)*h(i)-pcpc) ddist=dist1-dist2 ratio=ddist/(h(i-1)-h(i)) ttsf_sph(i)=ttsf_sph(i-1)+ps(i)*ratio ttf_sph(i)=ttf_sph(i-1)+dxs(i)*ratio dist1=dist2 enddo sum=0.0d0 c now spline log of slant total od. as a function of the log of veritcal od. do i=1,lmax ttf_sph(i)=log(ttf_sph(i)) enddo call splset(xs,ttf_sph,cofx,lmax) c c*****for each layer of model atmosphere find two layers in the c*****standard atmosphere straddling it and find total slant optical c*****path to reach each layer of model atmosphere- c*****using spline interpolation and then compute the slant tranmission c*****and store in transph c transph(1) = 1.0 j = 2 lmaxm1 = lmax - 1 do 110 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 110 100 dum1 = xs(i+1) - ttl(j) dum2 = ttl(j) - xs(i) dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) transph(j) = dexp(-dum3) j = j + 1 if (j .gt. ncp1) go to 120 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 100 110 continue c c*****find values for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 120 if (ttl(ncp1) .le. xs(lmax)) go to 125 dum1 = xs(lmax) - ttl(ncp1) dum2 = ttl(ncp1) - xs(lmaxm1) i = lmaxm1 dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) transph(ncp1) = dexp(-dum3) 125 continue refdir(im)=transph(ncp1) c c now need to calculate the spherical scattering optical thickness of c each layer. Already have the scattering optical depth along the slant path, c so can use a spline of log(tts_sph) versus xs and interpolate to the c points ttl and use the same method to compute the average scattering c optical thickness of each layer as is done in dtaus do i=1,lmax ttsf_sph(i)=log(ttsf_sph(i)) enddo call splset(xs,ttsf_sph,cofx,lmax) c now interpolate to ttl to get slant scattering od. c c*****for each layer of model atmosphere find two layers in the c*****standard atmosphere straddling it and find slant scattering optical c*****path to reach each layer of model atmosphere- c*****using spline interpolation and store in tts_sph c tts_sph(1) = 0.0 j = 2 lmaxm1 = lmax - 1 do 140 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 140 130 dum1 = xs(i+1) - ttl(j) dum2 = ttl(j) - xs(i) dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) tts_sph(j) = dum3 j = j + 1 if (j .gt. ncp1) go to 150 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 130 140 continue c c*****find values for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 150 if (ttl(ncp1) .le. xs(lmax)) go to 155 dum1 = xs(lmax) - ttl(ncp1) dum2 = ttl(ncp1) - xs(lmaxm1) i = lmaxm1 dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1 dum2**2 + cofx(4,i)) dum3 = dexp(dum3) tts_sph(ncp1) = dum3 155 continue c c now calc. dtts_sph and dtsp_sph c dtts_sph(1) = 0.5*tts_sph(2) dtsp_sph(1) = 0.5*tts_sph(2) do i=2,nc dtsp_sph(i) = 0.5*(tts_sph(i+1) - tts_sph(i-1)) dtts_sph(i) = 0.5*(tts_sph(i) - tts_sph(i-1)) enddo dtts_sph(ncp1) = 0.5*(tts_sph(ncp1)-tts_sph(nc)) dtsp_sph(ncp1) = dtts_sph(ncp1) c c now calc. the integration factors c dum1=emu(im) do i=1,ncp1 extmu(i,im)=dtsp_sph(i)*transph(i)*dum1 enddo 199 continue c c portion lifted from expone !def c do 400 i=1,ncp1 c initialize ekary do 200 j=1,5 ekary(i,j)=0.d0 200 continue c c calculate w(t)dt*attenuation factor (see eqn. 3.6,6.7) c do 205 j=1,imu extmux(i,j)=extmu(i,j)/emu(j) !def 205 continue if (i.eq.1) go to 255 im1=i-1 c do 210 j=1,5 ekary(im1,j)=0. 210 continue c if(i.eq.2) go to 235 im2=i-2 c c calculates the elements of k1 matrix (see eqn 4.11) c evaluate k1(1),k1(2) of eqn 4.7 for current value of c argument (tt(i)-tt(j)) c do 225 j=1,im2 call dexpk(tt(i),tt(j)) c c ekary array contains k matrix*w(t)dt c do 220 jk=1,5 ekary(j,jk)=dtsp(j)*eek(jk) 220 continue if(i.lt.ncp1) go to 225 c new statements if(ipol.eq.0)then ek4(j)=(e(2)+e(4))*dtsp(j) ek5(j)=sq2*(e(2)-e(4))*dtsp(j) else ek4(j)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(j) ek5(j)=delp*(e(2)-e(4))*dtsp(j) endif c end of new statements 225 continue c call dexpk(tt(i),tt(im1)) do 230 k=1,5 ekary(im1,k)=dtts(im1)*eek(k) 230 continue c if(i.lt.ncp1) go to 235 c new statements if(ipol.eq.0)then ek4(nc)=dtsp(nc)*(e(2)+e(4)) ek5(nc)=sq2*dtsp(nc)*(e(2)-e(4)) else ek4(nc)=((1.0d0+2.0d0*q)*e(2)+e(4))*dtsp(nc) ek5(nc)=delp*(e(2)-e(4))*dtsp(nc) endif c end of new statements 235 call dexpk1(tt(i),tt(im1)) c do 240 k=1,5 ekary(im1,k)=ekary(im1,k)+dtts(i)*eek(k) 240 continue c do 250 k=1,5 ekary(i,k)=dtts(i)*eek(k) 250 continue if(i.lt.ncp1) go to 255 c new statemnets if(ipol.eq.0)then ek4(ncp1)=dtsp(ncp1)*(odan(2)+odan(4)) ek5(ncp1)=sq2*dtsp(ncp1)*(odan(2)-odan(4)) else ek4(ncp1)=((1.0d0+2.0d0*q)*odan(2)+odan(4))*dtsp(ncp1) ek5(ncp1)=delp*(odan(2)-odan(4))*dtsp(ncp1) endif c end of new statements go to 300 255 ip1=i+1 call dexpk1(tt(i),tt(ip1)) c do 260 k=1,5 ekary(i,k)=ekary(i,k)+dtts(ip1)*eek(k) 260 continue c do 270 k=1,5 ekary(ip1,k)=dtts(ip1)*eek(k) 270 continue c if(i.eq.nc) go to 300 ip2=i+2 call dexpk(tt(i),tt(ip1)) c do 280 k=1,5 ekary(ip1,k)=ekary(ip1,k)+dtts(ip2)*eek(k) 280 continue c do 295 j=ip2,ncp1 call dexpk(tt(i),tt(j)) c do 290 k=1,5 ekary(j,k)=dtsp(j)*eek(k) 290 continue 295 continue c c 300 continue do 320 ii = 1,202 do 310 iii = 1,5 nek(ii,iii,i) = ekary(ii,iii) 310 continue 320 continue 400 continue c c optionally print calculated exponential terms c (zstar matrix * w(t)dt see eqn 6.3) c if (jprint(3).ne.0) then write (33,6120) (j,ek4(j),ek5(j) ,j=1,ncp1) do 495 i=1,imu thta=dacos(emu(i))/cnvrt write(33,7120) thta,(j,extmux(j,i),j=1,ncp1,10) write(33,7130) ncp1,extmux(ncp1,i) 495 continue endif return 6120 format(1h1,t50,'debug print out for expone',///, 1 1h ,2(' j',4x,'ek4(j)',6x,'ek5(j)',6x),//,2(1x,i3,2d12.4,3x)) 7120 format(1x,//,1x,'theta=',f6.2,' *** attenuation factors/*** ',//, 1 1x,'level',5x,'extmu/mu',/,(1x,i5,3x,d12.6)) 7130 format(1x,i5,3x,d12.6) c end of portion lifted from expone.f !def end subroutine firstz(h,ps,z,thenot,fracin,lmax,ncp1) c c*********************************************************************** c cccc c subroutine firstz c c purpose- c correct the optical depths calculated for the layers of a plane c parallel atmosphere for the sphericity of the actual atmosphere c c method- c 1.using the optical depths of the layers of the standard atmos. c calculate the slant optical path of the solar rays thru them c using the chapman function. c 2.then using spline interpolation calculate the slant optical c paths in the model atmosphere. c c calling sequence- c call firstz(h,ps,z,thenot,fracin,fs,f1,f2,lmax,ncp1) c c variables type i/o description c --------- ---- --- ----------- c h(487) r*8 i height of layers of standard atmosphere c from earth center (units of r) c ps(487) r*8 i rayleigh optical thickness of c each layer of std. atmosphere c z(202) r8 o total optical path of solar ray to c reach each layer of the model atmos. c thenot r*8 solar zenith angle c fracin r*8 i layer*r (radius of earth in units of c layer thickness) c fs r*8 o cos(thenot)*z(202) c f1 r*8 o 3/16(cos(thenot)**2) c f2 r*8 o 3*sqrt(2)/16(sin(thenot)**2) c lmax i*4 i # layers of standard atmos. c ncp1 i*4 i # layers in model atmos. c c external references c omerf c splset 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. cccc c********************************************************************** c c ****************************************************************** c c using the slant path optical thicknesses and chapman function c this subroutine calculates dave's z matrices for primary c scattering. c c ****************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) real *8 h(487),ps(487),xs(487),dxs(487),zs(487),cofx(4,487) real *8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) 1 ,z(202) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,fsp real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),fsp(15) common /thkns/tsl,tt,tts,ttl,dtsp,dtts common /chpmn/chpn,chxn,sqchp,sqchx,xs,dxs c c*****calculate trig functions to be used in do loops c amuo = dcos(thenot*cnvrt) if (thenot .eq. 90.0) amuo = 0. amuosq = amuo**2 sn = dsqrt(1.-amuosq) f1 = 0.1875*(1. + amuosq) f2 = 0.1875*sq2*(1. - amuosq) c c*****calculate slant opt. path of solar rays to reach c*****layer(i) of standard atmosphere c do 200 i = 1, lmax zs(i) = 0. raycon = sn*h(i) if (i .eq. 1) go to 190 duma= dsqrt((h(1) + raycon)*(h(1) - raycon)) c c*****add slant optical paths thru each layer(j) c*****lying between layer(1) and layer(i) c do 180 j = 2, i dumb= dsqrt((h(j) + raycon)*(h(j) - raycon)) path = fracin*(duma - dumb) zs(i) = zs(i) + path*dxs(j) 180 duma = dumb c c*****find the slant opt. path to reach layer 1 from the top of c*****the atmosphere- using chapman function approximation c 190 amu = raycon/h(1) amu = dsqrt(1. - amu**2) chpp = omerf(sqchp*amu, amu, chpn) chxx = omerf(sqchx*amu, amu, chxn) c c*****add the previous results to obtain total slant path zs to reach c*****each layer i. c zs(i) = zs(i) + chpp*ps(1) + chxx*(dxs(1) - ps(1)) 200 continue c c*****fit a spline between xs=log(vert. path) and zs=log(slant path) c*****through the standard atmosphere c do 215 i = 1, lmax 215 zs(i) = dlog(zs(i)) call splset(xs,zs,cofx,lmax) c c*****for each layer of model atmosphere find two layers in the c*****standard atmosphere straddling it and find total slant optical c*****path (z) to reach each layer of model atmosphere- c***** using spline interpolation c z(1) = 1.0 j = 2 lmaxm1 = lmax - 1 do 230 i = 1, lmaxm1 if (ttl(j) .le. xs(i) .or. ttl(j) .gt. xs(i+1)) go to 230 220 dum1 = xs(i+1) - ttl(j) dum2 = ttl(j) - xs(i) dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1dum2**2 + cofx(4,i)) dum3 = dexp(dum3) z(j) = dexp(-dum3) j = j + 1 if (j .gt. ncp1) go to 240 if (ttl(j) .gt. xs(i) .and. ttl(j) .le. xs(i+1)) go to 220 230 continue c c*****find z for the bottom layer by extrapolation if it was not c*****obtained during the above interpolation c 240 if (ttl(ncp1) .le. xs(lmax)) go to 245 dum1 = xs(lmax) - ttl(ncp1) dum2 = ttl(ncp1) - xs(lmaxm1) i = lmaxm1 dum3 = dum1*(cofx(1,i)*dum1**2 + cofx(3,i)) + dum2*(cofx(2,i)* 1dum2**2 + cofx(4,i)) dum3 = dexp(dum3) z(ncp1) = dexp(-dum3) 245 continue fs = amuo*z(ncp1) c write (6,'(i4,e12.4)') (i,z(i),i=1,ncp1) return end subroutine frstz2(z,thenot,fsx,f1,f2,ncp1) c ****************************************************************** c subroutine frstz2 c c purpose c this routine is used to calculate the zero order source c functions for a parallel plane atmos. model. c c method c the zero order source functions are computed by calculating c the intensity of the incident solar radiation at each source c function level using attenuation factors calculated for a c parallel plane atmosphere. c c calling sequence c call frstz2(z,thnot,fsx,f1,f2,ncp1) c c variable type i/o description c -------- ---- --- ----------- c c z(202) r*8 o zero order reduced source functions c thnot r*8i current solar zenith angle c fsx r*8 o direct flux reaching ground c f1 r*8 o constant c f2 r*8 o computational constant c ncp1 i*4 i # levels in model atmos. c c analysis and programming c k. f. klenk , p. m. smith sasc aug 77 c c modifications (date name purpose) c last modified by zia ahmad c purpose: to include the effect of molecular anisotropy 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 c c implicit integer*4(i-n),real*8 (a-h,o-z) c real *8 h(487),ps(487),xs(487),dxs(487),zs(487),cofx(4,487) real *8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) 1 ,z(202) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /thkns/tsl,tt,tts,ttl,dtsp,dtts common /chpmn/chpn,chxn,sqchp,sqchx,xs,dxs common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,totint c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),totint(15,8) amuo = dcos(thenot*cnvrt) if (thenot .eq. 90.0) amuo = 0. amuosq = amuo**2 sn = dsqrt(1.-amuosq) c new statements if(ipol.eq.0)then f1 = 0.1875d0*(1.d0 + amuosq) f2 = 0.1875d0*sq2*(1.d0 - amuosq) else f1=0.25d0*q1*(1.0d0+amuosq+2.0d0*q) f2=0.25d0*q1*delp*(1.0d0-amuosq) endif c*****parallel plane atmosphere do 200 j=1,ncp1 z(j)=dexp(-tt(j)/amuo) 200 continue c fsx=amuo*z(ncp1) fs=fsx return end subroutine geopro(x) c c********************************************************************** cccc c subroutine geopro c c version aug 22,1977 c c purpose c c compute extrapolated value of argument c c method c c using last three iterated values of argument, a geometric c series approx. is used to compute extrapolated terms. c c calling sequence c c call geopro(x) c c variable type i/o description c -------- ---- --- ----------- c c x(4) r*8 i/o argument to be extrapolated c x(1),x(2),x(3) contain last three c iteration values of argument to be c extrapolated. extrapolated term is c stored in x(4) c c external references c none c c common areas referenced c none c c analysis and programming c k.f. klenk p.m. smith sasc, aug 22 1977 c c modifications (date name purpose) c c last modified 11/18/92 bye zia ahmad c fixed the overflow problem c cccc c********************************************************************** c implicit integer*4(i-n),real*8(a-h,o-z) c real*8 x(4) c a=x(3)-x(2) b=x(2)-x(1) c = dabs(b-a) if (c .gt. 1.0d-9) then x(4) = x(2) + a*b/(b-a) else x(4) = x(3) endif return end logical function get_coef(start, ending, lambda, c0, c1, c2, beta, * rhon) C======================================================================C C C C FUNCTION : GET_COEF C C C C PURPOSE: This subroutine reads the ozone absorption and C C Rayleigh scattering coefficients and splines the C C to the correct wavelenght C C C C USE: C C call get_coef(lambda, c0, c1, c2, beta) C C C C PARAMETERS : C C C C Lambda - Real*4 input parameter. The wavelength in C C angstroms. C C c0, c1, c2 - Real*4 output. Ozone absoption coef. C C Alpha = c0 + c1 * T + c2 * T**2, where C C T is the temperature in Celsius. C C beta - Real*4 Output. Rayleigh coefficient at that C C wavelength. C C C C SOURCE: DONALD J. RICHARDSON C C HUGHES STX C C C C LANGUAGE: FORTRAN 77 C C C C DATE STARTED : September 28, 1994 C C C C LATEST REVISION: September 28, 1994 C C C C C C======================================================================C save implicit none ! ! Parameters ! character*80 data_file_name_def parameter(data_file_name_def = 'BASS.DAT') ! ! Calling Parameters ! real*8 start ! Input -- Starting wavelength real*8 ending ! Input -- ending wavelength real*8 lambda ! Output -- Wavelength real*8 c0 ! Output -- Ozone absorption coef (const term) real*8 c1 ! Output -- Ozone Absorption coef (T(c deg) term) real*8 c2 ! Output -- Ozone Absorption coef (T(c deg)^2 term) real*8 beta ! Output -- Rayleigh scattering coef real*8 rhon ! Output -- depolarization factor ! ! Local Variables ! logical first_time /.true./ integer*4 iunit integer*4 iostat integer*4 num integer*4 additional character*80 junk character*80 data_file_name ! ! If this is the first time read the data file and initialize spline. ! if (first_time) then ! ! Get a free logical unit number ! call get_lun(iunit) ! ! Check to see if all LUN are in use ! if (iunit .lt. 0) then write(*, * '(''No Logical Units available to be open in get_coef'')') stop 'No LUN avaiable' end if ! ! Open coefficient data file ! first_time = .false. call getenv('ABSORP_COEF', data_file_name) if (data_file_name .eq. ' ') then c print*,'Absorption coefficient enviromental variable ', c * '(ABSORP_COEF) was not found.' c print*,'Using the default absoption coef. file', c * data_file_name_def data_file_name = data_file_name_def end if open (iunit, file=data_file_name, status = 'old', err = 200) read(iunit,'(a)') junk ! Header line ! ! locate the first data set ! 90 read(iunit,*,end=100) lambda, c0, c1, c2, beta, rhon if (lambda .lt. start) go to 90 get_coef = lambda .le. ending return end if ! ! Getting next data set ! read (iunit,*,end = 100) lambda, c0, c1, c2, beta, rhon get_coef = lambda .le. ending return ! ! End of file reached.... No more data ! 100 get_coef = .false. return 200 print*, 'Error in opening the coefficient database:', * data_file_name stop 'Error in opening coefficient database' end subroutine get_lun (ivalue) C======================================================================C C C C SUBROUTINE : GET_LUN C C C C PURPOSE: To keep logical unit numbers straight, and issues new C C ones when necessary. The entry point FREE_LUN closes C C file and allows the logical unit number to be reused. C C C C USE: C C CALL GET_LUN (IUNIT) C C CALL FREE_LUN (IUNIT) C C C C PARAMETERS : C C C C IUNIT --- In GET_LUN, it is the returned logical unit C C number. In FREE_LUN, it is the input LUN to C C close and to be reused. C C C C SOURCE: DONALD J. RICHARDSON C C HUGHES STX C C C C LANGUAGE: FORTRAN 77 C C C C DATE STARTED : September 9, 1994 C C C C C LATEST REVISION: September 12, 1994 C C C C C C======================================================================C implicit none integer*4 min_lun, max_lun parameter (min_lun = 10, max_lun=100) logical open_stat integer*4 ivalue integer*4 i do i = min_lun, max_lun inquire (i, opened = open_stat) if (.not. open_stat) then ivalue = i return end if end do write (0,*) 'All the Logical Unit Numbers between ',min_lun, * ' and ', max_lun, ' are all in use.' ivalue = -1 return entry free_lun (ivalue) close (ivalue) return end subroutine iniclz c*********************************************************************** c c this subroutine initializes various numerical constants freque- c ntly used in various subroutines. c c last modified 03/10/95...dave flittner c purpose: set radius of earth equal to value used to find earth c position (6378.145 km). 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*********************************************************************** implicit integer*4(i-n), real*8(a-h,o-z) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /es/ odan(6),e(6),eek(6) c c pi = dacos(-1.0d+00) pi = 3.1415926535d+00 cnvrt = pi / 180.d+00 c c***** r is the radius of the earth in km. c r = 6378.145d+00 !def rinv = 1.d+00 / r cons = r / 6393.d+00 sq2 = dsqrt(2.d+00) c215 = 2.d+00 / 15.0d+00 c815 = 4.0d+00 * c215 c38sq2 = 0.375d+00 * sq2 c1415 = 0.9333333333d0 c2815 = 1.8666666667d0 pscaleforg = 6.95d0/r !def kskip=0 numek = 606 c odan(1)=1.d0 do 10 i=2,6 odan(i)=1.d0/float(i-1) 10 continue c return end subroutine iniclz1 c*********************************************************************** c c last modified by zia ahmad 9/10/93 c purpose: to include the effect of molecular anisotropy c c*********************************************************************** implicit integer*4(i-n), real*8(a-h,o-z) c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of new statement c c new statements if(ipol.eq.1)then c rhon=0.035d0 gama=rhon/(2.0d0-rhon) q2=(1.0d0-gama)/(1.0d0+2.0d0*gama) q1=0.75d0*q2 q=(2.0d0*gama)/(1.0d0-gama) sdp=1.0d0+q delp=dsqrt(2.0d0+3.0d0*q) q12s=q1/(2.0d0*sdp) endif c end of new statements c c write(33,500)ipol 500 format('sub iniclz1 ... ipol=',i2) c write(33,501)rhon,gama,q,q1,q2,sdp,delp,q12s 501 format('sub iniclz..rhon etc.'/1p6e12.4/1p2e12.4) return end subroutine intsum c c last modified 03/07/95...dave flittner c purpose: Spherical correction to the integration of the out-going c beam for each view angle. Lines changed denoted by !def c implicit real*8(a-h,o-z),integer*4(i-n) real*8 ristar(15,11),tensig(11),total(15,8,11) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101) real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),totint(15,8),alb(11) real*8 ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) real*8 alpha0, beta, pnot, tautot real*8 x(101) common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef,lambda,layer common /totals/ ristar,tensig,total common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,totint common /emm/ematx,admatx,saz,caz,saz2,caz2,emu,emuz common /atmos/ alpha0,beta,code,pnot,tautot,wavlen c common /in/delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c common /depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol real*8 tt(15) real*8 totn(15,8,11) real*8 temp(15) c common/buff1/sbz(9),tnstrz(9,15),tnstrl(15),tnstrr(15) common/buff2/e0za(9,15),e1za(9,15,8),e2za(9,15,8),ggz(9) common/buff3/eitl(15,8),eitr(15,8),eitu(15,8), ttl(15), ttr(15) common/buff4/eittl(15,8,11),eittr(15,8,11),eitotz(15,8,11), 1 polz(15,8,11),ttz(9,15) c real*8 refdir(15) !def common/crefdir/refdir !def c c c**** compute total reflected intensity at top of atmosphere ristar c**** loop over albedos and polar look angles c itmaxp=itmax+1 maxpp=itmaxp+1 do 100 i=1,nalb c c**** compute ig, ground reflected radiation c tensig(i)=(fs+gg)*qr(i) do 200 j=1,imu c do 124 it=1,maxpp ttz(it,j)=(fs+ggz(it))*(tnstrz(it,j)+refdir(j)) !def 124 continue tt(j)=(fs+gg)*(tnstr(j)+refdir(j)) !def ttl(j)=(fs+gg)*(tnstrl(j)+0.50d0*refdir(j)) !def ttr(j)=(fs+gg)*(tnstrr(j)+0.50d0*refdir(j)) !def temp(j)=ttl(j)+ttr(j) c write(33,123)fs,gg,refdir(j),tnstr(j),tt(j),temp(j) 123 format('intsum...fs,gg,refdir,tnstr(j),tt,temp',6f8.4) c**** compute sum of direct (refdir) and diffuse (tnstr) reflected inten c ristar(j,i)=tensig(i)*(tnstr(j)+refdir(j)) !def 200 continue 100 continue c c**** sum ristar and totint to calculate the total radiation (scattered c**** and reflected (izero+i1+i2+istar) at top of atmosphere. c do 500 i=1,imu do 400 j=1,iazmth do 300 k=1,nalb c total(i,j,k)=totint(i,j)+ristar(i,k) totn(i,j,k)=-100.0*dlog10(total(i,j,k)) c if(i.eq.1 .and. j.eq.1 .and. k.eq.1)then c write(33,666)i,j,k,thnot(imuz),thta(i),azmth(j),alb(k), c 1 total(i,j,k),totn(i,j,k) 666 format('i,j,k,total,totn',3i3,3f7.1,f7.3,1p2e11.3) c endif 300 continue 400 continue 500 continue c c compute degree of polarization c do i=1,imu do j=1,iazmth do k=1,nalb gistl=qr(k)*ttl(i) gistr=qr(k)*ttr(i) eittl(i,j,k)=eitl(i,j)+gistl eittr(i,j,k)=eitr(i,j)+gistr eitotz(i,j,k)=eittl(i,j,k)+eittr(i,j,k) polz(i,j,k)=dsqrt((eittl(i,j,k)-eittr(i,j,k))**2 1 +eitu(i,j)**2)/eitotz(i,j,k) if (eittl(i,j,k) .gt. eittr(i,j,k)) * polz(i,j,k) = -polz(i,j,k) c c write(33,'(''intsum...eittl,eittr,eitu'',3f9.5)') c * eittl(i,j,k), eittr(i,j,k), eitu(i,j) c write(33,135)i,j,k,tnstrl(i),tnstrr(i),qr(k),ttl(i), c 1 ttr(i),eitotz(i,j,k),polz(i,j,k) enddo enddo enddo 135 format('im,jaz,kalb,tnstrl,tnstrr,qr,ttl,ttr,eitotz,polz'/ 1 3i4,7f9.5) pi=dacos(-1.0d0) c print izero, i1,i2 etc as func. of sza and view angle c c write(33,132) 132 format('theta0','theta',t15,'i0',t25,'i1',t35,'i2',t45,'t', 1 t55,'sb',t65,'rho',t73,'ipol'/) do i=1,imu tfunc=tt(i)/pi xizero=eizero(i)/pi xi1=eiaz1(i,1)/pi xi2=eiaz2(i,1)/pi c write(33,133)thnot(imuz),thta(i),xizero,xi1,xi2,tfunc,sb(4), c 1 rhon,ipol 133 format(2f5.0,1x,6f10.5,1x,i3) enddo c write(33,134) 134 format(/ ) c return end 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. cccc c*********************************************************************** c implicit integer*4(i-n),real*8(a-h,o-z) c real*8 z(202),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 nek(202,5,202) c real*8 tm1(4,202),tm2(4,202),tm3(4,202),tm4(4,202), 1 tm5(4,202),tm6(4,202) c integer jprint(10) logical write_iter_file !def common /cwrite_iter/write_iter_file !def common /prints/ jprint common /kmat/nek common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip c itmm2=itmax-2 index=1 itmaxp=itmax+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,1) c if (itmax.eq.1) then 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 goto 101 endif 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 c c do first and higher order iterations c do 410 i=1,ncp1 do 5100 kk = 1,202 do 5101 kkk = 1,5 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) c save z matrix and zj1,zj2 of last three iterations c if(it.lt.itmm2) go to 410 c 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) 410 continue itp=it+1 if((itp .le. itmaxp).and. !def &(write_iter_file.or.(jprint(4).eq.1)))then !def call eva1pol(tm1,tm2,tm3,tm4,tm5,tm6,jmuz,ncp1,itp) endif if(it.ge.itmm2) index=index+1 140 continue index=index-1 101 continue 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)) !def & 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 subroutine matmul(x,y,z) real*8 x(4),y(4),z(4) c c z(1)=x(1)*y(1)+x(2)*y(3) z(2)=x(1)*y(2)+x(2)*y(4) z(3)=x(3)*y(1)+x(4)*y(3) z(4)=x(3)*y(2)+x(4)*y(4) c return end function omerf (x, y, sq) c ****************************************************************** c c this subroutine calculates chapman function as approximated by c john a. fitzmaurice, appl. opt. 3,640(1964). c c method: c omerf(h,mu)= sqrt(h*pi/2) exp(h*(mu**2)/2) c c where- h=layer height from the center of earth c (in units of scale height of attenuating c species) c c c ****************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) real *8 a(5),dp,dy,dx,dt,dsum,dsq,domerf data dp/0.3275911d+00/ data a / 1 1.061405429d+00, -1.453152027d+00, 2 1.421413741d+00, -0.284496736d+00, 3 0.254829592d+00 / dy = y dx = x dsq = sq dt = 1.0d+00 / (1.0d+00 + dp * dx) dsum = dt * a(1) if (y .ge. 0.2) go to 20 do 10 i=2,5 10 dsum = dt * (a(i) + dsum) domerf = dsum * dsq omerf = domerf return c c*****for all but large solar zenith angles use the series expansion c***** of the chapman function c 20 dt = 0.5d+00 / (dx**2) domerf=(1.0d+00-dt*(1.d+00-3.d+00*dt*(1.d+00-5.d+00*dt)))/dy omerf = domerf return end subroutine opthik(x,p,t,tl,cofx,tmp0,tmpcof,tmp101) c c*********************************************************************** cccc c subroutine opthik c c purpose- c divide the atmosphere into layers and calculate their c absorption,scattering and total optical depths measured from c top of the atmosphere. c c method- c uses spline interpolation to calculate optical thickness values c on finer grid points than those obtained from the input data c c calling sequence- call opthik(x,p,t,tl,cofx) c c variable type i/o description c -------- ---- --- ----------- c c x(101) r*8 i cumulative ozone thicknesses c c p(101) r*8 o pressure at each layer (atmospheres) c t(101) r*8 o total optical depth of each layer c tl(101) r*8 o log(t) c cofx(4,487) r*8 o spline interpolation coeff. c c this subroutine also returns(thru' the common 'thkns' ) c c tsl(101) r*8 o log(scattering optical depth) c tts(202) r*8 o interpolated scattering opt. depth c tt(202) r*8 o interpolated total opt. depth c ttl(202) r*8 o log(tt) c c c external references- c splset c last modified......11/4/94 zia ahmad c purpose............to correct the raleigh optical thickness c for 1/r**2 effect of gravity c last modified......03/10/95 dave flittner c purpose............to set radius of earth equal to value in consts 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 to impliment the gravity correction to the rayleigh scattering c optical depth. c cccc c*********************************************************************** c implicit integer*4(i-n),real*8 (a-h,o-z) real *8 p(101),t(101),tl(101),tsl(101),tts(202),ttl(202), 1 tt(202),dtsp(202),dtts(202),x(101),cofx(4,487),tmpcof(2), & tmp0,tmp101(101),alpha0,beta integer jprint(10) common /prints/ jprint c common /atmos/ alpha0,beta,code,pnot,tautot,wavlen common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /thkns/tsl,tt,tts,ttl,dtsp,dtts logical lgcorrect !def common/cgcorrect/lgcorrect !def c c*****divide the atmosphere into 101 layers by a linear interpolation c*****of log(p) (pmax=1, pmin=1.0e-5) c*****store p and calculate log(ts),t and log(t) c elg = dlog(10.d+00) dum2a = 0.05*elg betal = dlog(beta) do 10 i = 1, 101 plog = dum2a*float(i-101) p(i) = dexp(plog) c c determine the height above the ground(p=1.0) using a constant c scale height of H=6.95 and the g(h). Assume the radius of earth c as: r0=r (in commmon block consts) and then compute the ratio: g0/g(h) c the value of H is set in iniclz and is divided by the radius of !def c the earth. !def if(lgcorrect)then !def gr=(1.0d0-pscaleforg*dlog(p(i)))**2 !def grl=dlog(gr) else !def gr=1.0d0 !def grl=0.0d0 !def endif !def c deltmp=tmp101(i) - tmp0 alpha=alpha0+tmpcof(1)*deltmp+tmpcof(2)*deltmp**2 ta = alpha*x(i) C ts = beta*p(i) ts = beta*p(i)*gr t(i) = ta + ts tl(i) = dlog(t(i)) C tsl(i) = betal + plog tsl(i) = betal + plog + grl if (jprint(5) .eq. 1) 1 write(33,101) i,x(i),p(i),t(i),tsl(i),tl(i),tmp101(i) 10 continue c c*****fit a spline between log(ts) and log(t) c call splset (tsl,tl,cofx,101) c c c*****increase the no. of layers to 202 by preferentially adding c*****layers at the atmosphere bottom using a quadratic function. c cnsta = -3.4219e-3 cnstb = 6.7056e-2 tt(1) = 0. tts(1) = 0. jb = 2 do 20 i = 2, 202 adum = float(202-i) ttsl = betal + cnsta*adum*(1. + cnstb*adum) tts(i) = dexp(ttsl) c c*****search for two layers in the old atmos. stradlling a given c*****layer in the new atmos. and find tt by spline interpolation c do 15 j = jb, 101 if (ttsl .gt. tsl(j)) go to 15 dum1 = tsl(j) - ttsl kb = j - 1 dum2 = ttsl - tsl(kb) dum3 = dum1*(cofx(1,kb)*dum1**2 + cofx(3,kb)) + dum2*(cofx(2,kb)* 1dum2**2 + cofx(4,kb)) ttl(i) = dum3 tt(i) = dexp(dum3) go to 20 15 continue 20 continue return 101 format(i4,6d15.4) end subroutine prfmod(layer,xsav,xdel) c author: c seftor c date: 4-mar-91 c purpose: provide interface for profile modification studies implicit integer*4(i-n), real*8(a-h,o-z) parameter (max_num_iter = 12) real*8 xprf(11), xsav(11), tmpprf(11), theta(10), scan(15) integer*4 nthet, nscan, nwavl integer*4 iter(max_num_iter) real*8 wave_iter(max_num_iter) common /input/ xprf, tmpprf, pres, wave_iter, iter, theta, * scan, wave_start, wave_stop, nthet, nscan, num_iter do 10 i = 1,11 xprf(i) = xsav(i) 10 continue xchng = 0.1 * xprf(layer) xprf(layer) = xprf(layer) - xchng xdel = -xchng return end subroutine readin c author: c. seftor c date: 28-jan-1991 c purpose: reads in all necessary input for the mateer code implicit integer*4(i-n), real*8(a-h,o-z) parameter (max_num_iter = 12) integer*4 nthet, nwavl, iter(max_num_iter), jprint(10), nscan real*8 xprf(11), tmpprf(11), wavel(12), alphp(12), betp(12) real*8 wave_iter(max_num_iter) real*8 tcof1(12), tcof2(12), pres, theta(10), scan(15) real*8 thta(15),azmth(8),x(101),pshold(101),hhold(101), 1 alb(11),azm(8),thnot(10) character prfnam*8 character wavnum*27 character file_name*200 logical lwav(12) common /input/ xprf, tmpprf, pres, wave_iter, iter, theta, * scan, wavel_start, wavel_stop, nthet, nscan, num_iter common /inchr/ prfnam, wavnum common /log/ lwav common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /in/ delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath common /prints/ jprint c c new statement common/depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol c end of statement nwavl = 12 if (iargc() .eq. 0) then file_name = 'PROF' else call getarg(1, file_name) end if open (unit=22, file=file_name, status='old') read (22,'(a)') prfnam read (22,*) pres read (22,*) nthet read (22,*) (theta(num),num=1,nthet) read (22,*) nscan read (22,*) (scan(num), num=1, nscan) read (22,*) naz read (22,*) (azmth(num), num=1, naz) read (22,*) nalb read (22,*) (alb(num), num=1, nalb) read (22, *) wavel_start, wavel_stop read (22,*) (xprf(num),num=11,1,-1) read (22,*) (tmpprf(num),num=11,1,-1) read (22,*) (jprint(i),i=1,10) read (22,*) num_iter if (num_iter .gt. max_num_iter) then print*,'Number of iteration ranges requested was ', num_iter print*,'Maximum allowed is ', max_num_iter,' Program stopped.' stop 'Too many iteration ranges requested.' end if read (22,*) (wave_iter(i),i=1,num_iter) read (22,*) (iter(i),i=1,num_iter) read (22,*) ipol close (22) imuz = nthet imu = nscan iazmth =naz do i=1, nthet thnot(i) = theta(i) end do do i = 1, nscan thta(i) = scan(i) end do if (jprint(9).eq.1) then write (*, '(''PRFNAM = '',a)') prfnam write (*, *) 'PRES = ', pres write (*,*) 'NTHET = ', nthet write (*,*) 'THETA = ', (theta(num),num=1,nthet) write (*,*) 'NSCAN = ',nscan write (*,*) 'SCAN = ', (scan(num),num=1,nscan) write (*,*) 'NAZ = ',naz write (*,*) 'azmth = ', (azmth(num),num=1, naz) write (*,*) 'NALB = ',nalb write (*,*) 'ALB = ', (alb(num),num=1,nalb) write (*, *) 'START, stop = ', wavel_start, wavel_stop write (*,*) 'Ozone profile =',xprf write (*,*) 'Temp profile =',tmpprf write (*,*) 'jprint =', (jprint(i),i=1,8) write (*,*) 'NUM_ITER = ', num_iter write (*,*) (wave_iter(i), i=1,num_iter) write (*,*) (iter(i),i=1,num_iter) write (*,*) 'IPOL = ', ipol endif return end subroutine reflex(itmax,nek,ncp1) c c*********************************************************************** cccc c subroutine reflex c c version aug 22,1977 c c purpose c c reflex is a fortran iv routine to calculate the fract. of surf. c reflected radiation which is back scattered by the atmosphere. c it also computes the ratio of the intensity of light emerging c from the top of the atmosphere (which was surface reflected ) c to the intensity of light reflected by the lambertian surface. c this quantity is tenstr=i-star /i-g,see eq 6.7 c c method c c the scattering of radiations scattered by the ground is c computed.the z-star matrix elements are calculated. c by iteration. c c calling sequence c c call reflex (nek,itmax,ncp1,tnstr,qr) c c variable type i/o description c -------- ---- --- ----------- c c nek i*4 i fortran data set unit # for c file containing exponential integrals c created by expone c itmax i*4 i number of iteration steps to be used c ncp1 i*4 i # layers c tnstr(15) r*8 o intensity ratio described above for c each polar back scatter angle. c qr(10) r*8 o factor to be used in computing i-sb-g c see eqn 6.7 c c external references c evalrf c c common areas referenced c c eks c thkns c analysis and programming c k.f. klenk, p.m. smith sasc,aug 21 1977 c c modifications (date name purpose) c c last modified 11/19/92....zia ahmad c modified for single iteration c last modified 10/19/94....zia ahmad c modified to print results after each itration c last modified 03/07/95....dave flittner c modified to initialize itp to zero, so that one iteration c may be done. 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 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. ccc c********************************************************************** c implicit integer*4(i-n),real*8(a-h,o-z) common /eks/ extmu(202,15),ek4(202),ek5(202) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip c common /thkns/ tsl,tt,tts,ttl,dtsp,dtts real *8 tsl(101),tt(202),ttl(202),dtsp(202),dtts(202),tts(202) c real*8 ekary(202,5),b1(202),b2(202),z1(202),z2(202), 1 zst1(4,202),zst2(4,202),tnstr(15) real*8 nek(202,5,202) real*8 temp1(4,202),temp2(4,202) logical write_iter_file !def common /cwrite_iter/write_iter_file !def c c c construct zeroth z-star vector c itp=0 !def itmm2=itmax-2 index=1 itmaxp=itmax+1 c do 100 i=1,ncp1 z1(i)=ek4(i)/dtsp(i) z2(i)=ek5(i)/dtsp(i) c if (itmax.eq.1) then zst1(1,i) = z1(i) zst2(1,i) = z2(i) temp1(1,i)=z1(i) temp2(1,i)=z2(i) c endif 100 continue if(write_iter_file) call eva2pol(temp1,temp2,ncp1,1) !def if (itmax.eq.1) goto 101 c c do first and higher order iterations c do 140 it=1,itmax c c set b1,b2 equal to new z1,z2 c do 150 k=1,ncp1 b1(k)=z1(k) b2(k)=z2(k) 150 continue c c do 280 i=1,ncp1 sum1=ek4(i)/dtsp(i) sum2=ek5(i)/dtsp(i) c do 5100 ii = 1,202 do 5101 iii = 1,5 ekary(ii,iii) = nek(ii,iii,i) 5101 continue 5100 continue c do 190 j=1,ncp1 sum1=sum1+ekary(j,1)*b1(j)+ekary(j,2)*b2(j) sum2=sum2+ekary(j,2)*b1(j)+ekary(j,3)*b2(j) 190 continue c c z1(i)=sum1 z2(i)=sum2 c temp1(1,i)=z1(i) temp2(1,i)=z2(i) c save zstar matrix of last three iterations c if(it.lt.itmm2) go to 280 c zst1(index,i)=z1(i) zst2(index,i)=z2(i) 280 continue c itp=it+1 if((itp .le. itmaxp).or.write_iter_file)then !def call eva2pol(temp1,temp2,ncp1,itp) endif c if(it.ge.itmm2) index=index+1 140 continue 101 continue c c call evalrf to calculate istar/ig of eq (6.7) c for all polar look angles c itpp=itp+1 if(write_iter_file) call eva2pol(zst1,zst2,ncp1,itpp) !def call evalrf(zst1,zst2,ncp1) c return end subroutine relayr(delx,xtop,tmp10,tmptop,x101,tmp101, 1 pnot,xpnot,h,ps) c ---------------------------------------------------------------- c subroutine relayr c c c version 2.0 june 1984 c algorithm designed by dr. p.k. bhartia ; coded by david lee ,sasc c c purpose - c c find the cumulative ozone amount and ozone-weighted average c temperatures at 101 levels from input values for 10 umkehr c layers using a spline fit. c c c method / procedures c c 1) define input and output pressure levels c 2) compute accumulated ozone amounts at input levels c 3) compute ozone weighted average temperatures above input c pressure levels c 4) perform a spline fit to temperatures and accumulated c ozones at input levels c 5) evaluate spline fit at the lowest 61 pressure levels c 6) find slope for the accumulated ozone at the 41st. level c 7) evaluaate accumulated ozone amounts for the top 40 levels c from the slope computed from the 41st level c 8) compute pressures for altitudes at 1 km intervals c c c calling sequence - call relayr(delx,xtop,tmp10,tmptop,x101,tmp101, c h,ps) c c c variable type i/o description c -------- ---- --- ----------- c c delx(10) r*8 i layer ozone amount (d.u) for c 10 umkehr layers. (layers 0-9, layer c 0 being the bottom half of the atmos- c phere). indexing of the layers are c described in the table below c c xtop r*8 i cumulative ozone amount(m-atm-cm) c above umkehr layer 9 c c tmp10(10) r*8 i average temperatures (degrees kelvin) c for 10 umkehr layers c c tmptop r*8 i average temperature for the c atmosphere above umkehr layer 9 c c p11(11) r*8 pressure levels (atm.) defining the 10 c umkehr layers. indexing of the pressure c levels are described in the table c below c c x11(11) r*8 cumulative ozone amount (m-atm-cm) c at pressure levels described above c for p11 c c pnot r*8 i pressure of reflecting surface c c xpnot r*8 o cumulative ozone at pnot c c p101(101) r*8 pressure (atm.) at 101 atmospheric c levels that define the output layer c ozone amounts and temperatures c c x101(101) r*8 o output cumulative ozone amounts c (atm-cm) at levels corresponding c to p101. c c c tmp101(101) r*8 o output ozone weighted temperatures c (deg-kelvin) for atmospheres above c each of the 101 pressure levels. c c z11(11) r*8 altitude (km) corresponding to p11 c defining the umkehr layers c c h(82) r*8 o altitude array at 1 km interval from c ground to 82 km c c ps(82) r*8 o pressures corresponding to the height c at the h array c c c----------------------------------------------------------------------- c c----------------------------------------------------------------------- c c the table below describes the order of indices for variables c in this subroutine c c pressure umkehr indices indices index c (atm.) layer # for p11,x11 for delx for c & t11ave & tmp10 z11 c c 1/1024 -------------------------1--------------------------11 c 9 1 c 1/512---------------------------2--------------------------10 c 8 2 c 1/256---------------------------3---------------------------9 c 7 3 c 1/128---------------------------4---------------------------8 c 6 4 c 1/64 ---------------------------5---------------------------7 c 5 5 c 1/32 ---------------------------6---------------------------6 c 4 6 c 1/16 ---------------------------7---------------------------5 c 3 7 c 1/8 ---------------------------8---------------------------4 c 2 8 c 1/4 ---------------------------9---------------------------3 c 1 9 c 1/2 --------------------------10---------------------------2 c 0 10 c 1 --------------------------11---------------------------1 c c c last modified 03/10/95...dave flittner c purpose: set earth radius equal to value in common block consts 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-------------------------------------------------------------------- implicit integer*4(i-n),real*8(a-h,o-z) integer jprint(10) real*8 p11(11),p11log(11),p101(101),p101lg(101), 1 x11(11),x11log(11),x101(101),x101lg(101),delx(10), 2 tmp10(10),tmp101(101),t11ave(11) dimension z11(11),plginv(11),ps(82),h(82) dimension cx(11),bparx(4),ct(11),bpart(4),cz(11),bparz(4) common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /prints/ jprint c c*****define presures (atm.) for layer ozone amounts and temperatures. c p11(11) is the at the bottom of the atmosphere c do 10 i=1,11 p11(i)=1.0*2.**(i-11.) p11log(i)=dlog(p11(i)) 10 continue c c*****define pressures (atm.) at 101 levels for output ozone and c temperatures p101(101) is at the bottom of the atmosphere c do 20 i=1,101 p101(i)=1.0*10.**((i-101.)/20.) p101lg(i)=dlog(p101(i)) 20 continue c c*****compute cumulative ozone at 11 pressure levels & convert to atm-cm c x11(1)=xtop/1000. x11log(1)=dlog(x11(1)) do 110 i=2,11 x11(i)=x11(i-1)+delx(i-1)/1000. x11log(i)=dlog(x11(i)) 110 continue c c*****compute ozone weighted average temperatures at 11 pressure level c tsum=tmptop*xtop/1000. t11ave(1)=tmptop do 150 i=2,11 tsum=tsum+tmp10(i-1)*delx(i-1)/1000. t11ave(i)=tsum/x11(i) 150 continue c c*****set up boundary value parametes for cumulative ozone spline fit c bparx(1)=0. bparx(2)=0. bparx(3)=1. c c*****f'(nx)=1.707*layer 0 ozone / total ozone c deltax=x11log(11)-x11log(10) deltap=p11log(11)-p11log(10) fnx=1.707*delx(10)/1000./x11(11) bparx(4)=6.0/deltap*(fnx-deltax/deltap) c c*****evaluate cumulative ozone at the lowest 61 output preure level c using a spline fit to point at input pressure levels c nx=11 ic=10 l1=41 l2=101 m=l2-l1+1 call spline(p11log,x11log,nx,1.d30,fnx,cx) do 175 i = l1,l2 call splint(p11log,x11log,cx,nx,p101lg(i),x101lg(i)) 175 continue c c***** evaluate cumulative ozone at pnot c xpnot=x11(11) if (pnot.ge.p11(11)) go to 180 p0log=dlog(pnot) call splint(p11log,x11log,cx,nx,p0log,xp0log) c if(ier.ne.0) go to 900 xpnot=dexp(xp0log) 180 continue if (jprint(6).eq.1) write(33,2000) xpnot c c*****evaluate cumultive ozone at uppermost 40 levels based on c slope computed from the 41st and 42nd levels c slope=(x101lg(l1)-x101lg(l1+1))/(p101lg(l1)-p101lg(l1+1)) c l1m1=l1-1 do 200 l=1,l1m1 x101lg(l)=x101lg(l1)+slope*(p101lg(l)-p101lg(l1)) 200 continue c c***** convert log of cumulative ozone to cumulative ozone c do 250 l=1,101 x101(l)=dexp(x101lg(l)) 250 continue c c*****evaluate ozone-weighted average temperatures for atmospheres c above lowest 61 output pressure levels using a spline fit to c points at input pressure levels c bpart(1)=0. bpart(2)=0. bpart(3)=0. bpart(4)=0. call spline(p11log,t11ave,nx,1.d30,1.d30,ct) do 275 i = l1,l2 call splint(p11log,t11ave,ct,nx,p101lg(i),tmp101(i)) 275 continue c c*****evaluate weighted temperatures in the upper most 40 layers c based on constant lapse rate of -1.5 degrees per layer) 195 do 300 l=1,l1m1 tmp101(l)=tmp101(l1)-1.5*(l1-l) 300 continue c c*****compute heights z11 (km.) corresponding to the 11 input pressure c levels p11 (in inverse order) using the equations : c dz/dlogp=n*r*t/g c g=g0*(r0/(r0+z))**2 c c constant is g0/(n*r); z11(1) is at bottom of atmosphere c const=1.0/(-34.163) r0=r dzhalf=2.5 z11(1)=0. plginv(1)=p11log(11) do 400 i=2,11 j=12-i plginv(i)=p11log(j) dlogp=plginv(i)-plginv(i-1) z=z11(i-1)+dzhalf t=tmp10(j) dz=const*((r0+z)/r0)**2*t*dlogp z11(i)=z11(i-1)+dz dzhalf=dz/2.0 400 continue c c*****perform a spline fit to log(p) and z c bparz(1)=0. bparz(2)=0. bparz(3)=0. bparz(4)=0. c if(ier.ne.0) go to 900 call spline(z11,plginv,nx,1.d30,1.d30,cz) c c*****set up altitude arrays at 1 km constant intervals c do 320 i=1,82 h(i)=float(i-1)*1.0 320 continue c c*****evaluate log of pressure at 1 km intervals from the ground to c the top of umkehr layer 9 based on spline fit c m=z11(11)+1. c if(ier.ne.0) go to 900 do 345 i = 1,m call splint(z11,plginv,cz,nx,h(i),ps(i)) 345 continue c c*****evaluate log of pressues above umkehr layer 9 based on constant c scale height c slope=(ps(m)-ps(m-1))/(h(m)-h(m-1)) mp1=m+1 do 450 i=mp1,82 ps(i)=ps(m)+slope*(h(i)-h(m)) 450 continue c c***** convert log of pressure to pressure c do 460 i=1,82 ps(i)=dexp(ps(i)) 460 continue c c***** dump out variout input & ouput data c if (jprint(6).eq.1) 1 write(33,6200) tmptop,xtop,(i,p11(i),x11(i),t11ave(i),z11(12-i), 2 tmp10(i),delx(i),i=1,10),(i,p11(i),x11(i),t11ave(i),z11(12-i), 3 i=11,11) if (jprint(6) .eq. 1) write(33,6210) (h(i),ps(i),i=1,82) 900 continue return 2000 format(1x,'xpnot=',f10.5//) 6200 format(//' ozone,temperature & other variables in umkehr layers:'/ 1 /t28,'weighted',t51,'layer',t61,'layer'/1x,'index',t13,'p11',t23, 2 'x11',t32,'temp',t43,'z11',t52,'temp',t61,'ozone' 3 /t46,2f10.2/10(i3,2x,e10.4,f10.3,2f10.2/t46,2f10.2/), 4 i3,2x,e10.4,f10.3,2f10.2) 6210 format(/' pressures for first 82 km altitude :'/5(8x,'h',9x, 1 'ps',2x)/17(5(f10.1,e12.4)/)) end subroutine slant(h,ps,cofx,fracin,lmax,layer,pshold,hhold) c c*********************************************************************** cccc c subroutine slant c c purpose- c 1. subdivide the layers of the standard atmosphere and c calculate the scattering and total optical depths of the c atmosphere from top to each layer. c 2. using the above results calculate the chapman constants. c c method- c 1. the layers are subdivided using exponential interpolation c 2. the total o.d. are calculated using spline interpolation c c calling sequence- c call slant(h,ps,cofx,fracin,juzsfr,lmax,layer) c c variable type i/o description c -------- ---- --- ----------- c c h(487) r*8 i height of the layers above surface(km) c o height from earth center(fraction of r) c ps(487) r*8 i pressure at each layer (in atmosphere) c o rayleigh opt. thickness of each layer c of the standard atmosphere c cofx(4,487) r*8 i spline interpolation coeff. c fracin o layer*radius of earth(km) c juzsfr i*4 i print flag c lmax i*4 i no. of layers in standard atmosphere c layer i*4 i no. of subdivisions to be made of each c layer of the standard atmosphere c other variables are returned thru common block 'chpmn' 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. Perform gravity correction when c computing rayleigh optical depth and storing in array ps. c Also use logical switch lgcorrect c to impliment the gravity correction to the rayleigh scattering c optical depth. c cccc implicit integer*4(i-n),real*8 (a-h,o-z) real *8 h(487),ps(487),xs(487),dxs(487),hhold(100),pshold(101), 1 cofx(4,487) real *8 tsl(101),tt(202),tts(202),ttl(202),dtsp(202),dtts(202) integer jprint(10) common /prints/ jprint common /atmos/ alpha0,beta,code,pnot,tautot,wavlen common /consts/pi,cnvrt,r,rinv,cons,sq2,c215,c815,c38sq2, 1 c1415,c2815,pscaleforg,numek,kskip common /thkns/tsl,tt,tts,ttl,dtsp,dtts common /chpmn/chpn,chxn,sqchp,sqchx,xs,dxs logical lgcorrect !def common/cgcorrect/lgcorrect !def c c*****use the ps and h that were saved before c 80 do 82 i=1,lmax h(i)=hhold(i) 82 ps(i)=pshold(i) c c*****renumber the layers to put layer 1 near the top of the atmos. c***** (standard atmos. is read in a reverse order) c 85 continue lmaxd2 = lmax/2 do 100 i = 1, lmax if(lgcorrect)then !def ps(i) = beta*ps(i)*(1.0d0-pscaleforg*dlog(ps(i)))**2 !def else !def ps(i) = beta*ps(i) endif !def 100 continue do 110 i = 1, lmaxd2 k = lmax - i + 1 hold = h(i) h(i) = h(k) h(k) = hold hold = ps(i) ps(i) = ps(k) 110 ps(k) = hold c c*****subdivide the layers of std. atmos. by interpolation c*****assume pressure dependence between the layers of the form c***** p=pnot**(-h/hnot) c lmax = layer*(lmax - 1) + 1 lmaxm1 = lmax - 1 layrm1 = layer - 1 do 112 i = 1, lmax, layer j = lmax - i + 1 k = j/layer + 1 h(j) = h(k) 112 ps(j) = ps(k) frac = 1./float(layer) fracin = r*float(layer) lmaxml = lmax - layer do 114 i = 1, lmaxml, layer k = i + layer dum = (ps(k)/ps(i))**frac k = i do 114 j = 1, layrm1 k = k + 1 h(k) = h(k-1) - frac 114 ps(k) = ps(i)*dum**j c c*****use spline interpol. to obtain total opt. depth of each layer c***** in the standard atmosphere (xs) c j = 1 dum = dlog(ps(j)) do 130 i = 2, 101 if (dum .gt. tsl(i)) go to 130 k = i - 1 120 dum1 = tsl(i) - dum dum2 = dum - tsl(k) dum3 = dum1*(cofx(1,k)*dum1**2 + cofx(3,k)) + dum2*(cofx(2,k)* 1dum2**2 + cofx(4,k)) xs(j) = dexp(dum3) j = j + 1 if (j .gt. lmax) go to 140 dum = dlog(ps(j)) if (dum .le. tsl(i)) go to 120 130 continue c c*****calculate ps(rayleigh opt. thickness) and dxs(total opt. thick) 140 dxs(1) = xs(1) holdc = ps(1) do 150 i = 2, lmax holdd = ps(i) ps(i) = holdd - holdc holdc = holdd 150 dxs(i) = xs(i) - xs(i-1) do 160 i = 1, lmax xs(i) = dlog(xs(i)) 160 h(i) = 1. + h(i)*rinv if(jprint(7).ne.0) 1 write(33,6400)(h(i),ps(i),xs(i),dxs(i),i=1,lmax) dum = r*h(1) mm=20*layer+2 scalp=20./dlog(ps(mm)/ps(2)) scalx = scalp scalx=20./dlog((dxs(mm)-ps(mm))/(dxs(2)-ps(2))) chp = dum/scalp*0.5 chx = dum/scalx*0.5 chpn = dsqrt(chp*pi) chxn = dsqrt(chx*pi) sqchp = dsqrt(chp) sqchx = dsqrt(chx) 6400 format (1h1,2(5x,1hh,14x,2hps,14x,2hxs,13x,3hdxs,9x)/1h,/, 1 (1h ,2(f10.8, 3d16.4, 5x))) return end subroutine spline(x,y,n,yp1,ypn,y2) c author: c seftor c date: 3 feb 1991 c purpose: given arrays x and y of length n containing a tabulated c function, and given values yp1 and ypn for the first derivative of c the interpolating function at points 1 and n, respectively, this c routinereturns an array y2 of length n which contains the second c derivatives of the interpolating function at the tabulated value c points x. if yp1 and/or ypn are equal to 1 x 10**30 or larger, c routine is signalled to set the corresponding boundary condition for c a natural spline, with zero second derivative on the boundary. c (algorithm taken from numerical recipes (fortran), press et. al.) c implicit integer*4(i-n), real*8(a-h,o-z) parameter (nmax=487) real*8 x(n), y(n), y2(n), u(nmax) real*8 yp1, ypn integer*4 n if (yp1 .gt. .99d30) then y2(1) = 0. u(1) = 0. else y2(1) = -0.5 u(1) = (3./(x(2)-x(1))) * ((y(2)-y(1))/(x(2)-x(1)) - yp1) endif do 11 i = 2,n-1 sig = (x(i)-x(i-1))/(x(i+1)-x(i-1)) p = sig * y2(i-1) + 2. y2(i) = (sig-1.)/p u(i) = (6.*((y(i+1)-y(i))/(x(i+1)-x(i)) - (y(i)-y(i-1)) 1 /(x(i)-x(i-1)))/(x(i+1)-x(i-1)) - sig*u(i-1)) / p 11 continue if (ypn .gt. .99d30) then qn = 0. un = 0. else qn = 0.5 un = (3./(x(n)-x(n-1))) * (ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) endif y2(n) = (un - qn*u(n-1))/(qn*y2(n-1)+1.) do 12 k = n-1,1,-1 y2(k) = y2(k) * y2(k+1) + u(k) 12 continue return end subroutine splint(xa,ya,y2a,n,x,y) c author: c. seftor c date: 3 feb 1991 c purpose: given the arrays xa and ya of length n, which tabulate a c function (with the xa's in order), and given the array y2a, c which is the output from the subroutine spline, and given a value c of x, this routine returns a cubic-spline interpolated value y. c (algorithm taken from numerical recipes (fortran) by press et. al.) c implicit integer*4(i-n), real*8(a-h,o-z) real*8 xa(n), ya(n), y2a(n), x, y real*8 a, b integer*4 n klo = 1 khi = n 10 continue if (khi-klo .gt. 1) then k = (khi+klo)/2 if (xa(k) .gt. x) then khi = k else klo = k endif goto 10 endif h = xa(khi) - xa(klo) if (h .eq. 0.) then write (23,*)' bad xa input' stop endif a = (xa(khi)-x)/h b = (x-xa(klo))/h y = a*ya(klo) + b*ya(khi) + 1 ((a**3-a)*y2a(klo) + (b**3-b)*y2a(khi)) * (h**2)/6. return end subroutine splset(x,y,c,m) c*********************************************************************** c c this subroutine was taken from_ c introductory computer methods and numerical analysisby c ralph h. pennington c the mcmillan company -newyork 1965. c this subroutine will accept up to 487 points. c c***** x = independent variable c***** y = dependent variable c***** m = number of data points in the x and y arrays c***** c = coefficents of the cubic that is fit between adjacent data po c***** the following dimensions must be .ge. m c*********************************************************************** implicit integer*4(i-n),real*8 (a-h,o-z) real*8 d(487), p(487), e(487), a(487,3), b(487), z(487) real*8 x(1),y(1),c(4,1) mm=m-1 do 2 k=1,mm d(k)=x(k+1)-x(k) p(k)=d(k)/6. 2 e(k)=(y(k+1)-y(k))/d(k) do 3 k=2,mm 3 b(k)=e(k)-e(k-1) a(1,2)=-1.-d(1)/d(2) a(1,3)=d(1)/d(2) a(2,3)=p(2)-p(1)*a(1,3) a(2,2)=2.*(p(1)+p(2))-p(1)*a(1,2) a(2,3)=a(2,3)/a(2,2) b(2)=b(2)/a(2,2) do 4 k=3,mm a(k,2)=2.*(p(k-1)+p(k))-p(k-1)*a(k-1,3) b(k)=b(k)-p(k-1)*b(k-1) a(k,3)=p(k)/a(k,2) 4 b(k)=b(k)/a(k,2) q=d(m-2)/d(m-1) a(m,1)=1.+q+a(m-2,3) a(m,2)=-q-a(m,1)*a(m-1,3) b(m)=b(m-2)-a(m,1)*b(m-1) z(m)=b(m)/a(m,2) mn=m-2 do 6 i=1,mn k=m-i 6 z(k)=b(k)-a(k,3)*z(k+1) z(1)=-a(1,2)*z(2)-a(1,3)*z(3) do 7 k=1,mm q=1./(6.*d(k)) c(1,k)=z(k)*q c(2,k)=z(k+1)*q c(3,k)=y(k)/d(k)-z(k)*p(k) 7 c(4,k)=y(k+1)/d(k)-z(k+1)*p(k) return end c subroutine sumry(jmuz) c************************************************************************ c subroutine sumry prints out sb, t, i0 pol etc c************************************************************************ c implicit real*8(a-h,o-z),integer*4(i-n) real*8 ristar(15,11),tensig(11),total(15,8,11) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101) real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),totint(15,8),alb(11) real*8 ematx(3,15),admatx(3,15),saz(8),caz(8), 1 saz2(8),caz2(8),emu(15),emuz(10) real*8 alpha0, beta, pnot, tautot real*8 x(101) common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef,lambda,layer common /totals/ ristar,tensig,total common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,totint common /emm/ematx,admatx,saz,caz,saz2,caz2,emu,emuz common /atmos/ alpha0,beta,code,pnot,tautot,wavlen c common /in/delx(10),xtop,tmp10(10),tmptop,tmp0, 1 itmax,ipsudo,ipath c common /depolt/rhon,gama,q,q1,q2,delp,sdp,q12s,ipol real*8 tt(15) real*8 totn(15,8,11) c common/buff1/sbz(9),tnstrz(9,15),tnstrl(15),tnstrr(15) common/buff2/e0za(9,15),e1za(9,15,8),e2za(9,15,8),ggz(9) common/buff3/eitl(15,8),eitr(15,8),eitu(15,8), ttl(15), ttr(15) common/buff4/eittl(15,8,11),eittr(15,8,11),eitotz(15,8,11), 1 polz(15,8,11),ttz(9,15) c c maxpp=itmax+2 c c write(9,50) do i=1,maxpp do j=1,imu do k=1,iazmth c write(9,100)i,thta(j),azmth(k),e0za(i,j),e1za(i,j,k), c 1 e2za(i,j,k),ttz(i,j),sbz(i),ggz(i) enddo enddo enddo c write(9,150)thnot(jmuz),wavlen write(9,200) do i=1,nalb c do j=1,imu c do k=1,iazmth do k=1,iazmth do j=1,imu write(9,225)thta(j),azmth(k),eittl(j,k,i), 1 eittr(j,k,i),eitotz(j,k,i),polz(j,k,i),alb(i) enddo enddo enddo c return c 50 format(t2,'itr',t10,'the',t15,'phi',t26,'ei0',t37,'ei1', 1 t48,'ei2',t59,'t ',t67,'sb',t77,'gg') 100 format(i5,1x,f7.1,1x,f7.1,1p4e11.3,0pf8.5,1pe11.3) 150 format('solar zenith angle=',f6.1,' wavelength=',f7.1) 200 format(t4,'the',t12,'phi',t24,'eil',t35,'eir',t44,'eitot', 1 t55,'pol',t62,'alb') 225 format(f7.1,1x,f7.1,1x,1p3e11.3,0pf11.3,0pf7.3) c end subroutine tbprnt(thenot) c implicit real*8(a-h,o-z),integer*4(i-n) real*8 thnot(10),thta(15),azmth(8),pshold(101),hhold(101),alb(11) real*8 x(101) integer jprint(10) real*8 ristar(15,11),tensig(11),total(15,8,11) real*8 sb(4),tnstr(15),qr(11),eizero(15),eiaz1(15,8), 1 eiaz2(15,8),totint(15,8) real*8 thenot c common /totals/ ristar,tensig,total common /prints/ jprint common /contrl/ nalb,imu,imuz,iazmth,thta,azmth,alb, 1 x,pshold,hhold,thnot,alfaef, 2 lambda,layer common /out/ sb,qr,tnstr,eizero,eiaz1,eiaz2,gg,fs,totint c c if(thenot.ne.thnot(1)) go to 50 c c**** print current solar zenith angle,fractional back c**** scatter factor and qr values for each albedo. c if (jprint(8) .eq. 1) write(33,1000) thenot,sb(4), * (qr(i),i=1,nalb) c c**** print istar/ig at top of atmosphere. c if (jprint(8) .eq. 1) write(33,1100) * (i,thta(i),tnstr(i),i=1,imu) c c 50 continue c c**** print azimuth independent terms for scattered rad. c if (jprint(8) .eq. 1) write(33,1200) thenot,fs, * gg,(i,thta(i),eizero(i),i=1,imu) c c**** write eiaz1 heading c if (jprint(8) .eq. 1) write(33,2000) thenot, * (azmth(i),i=1,iazmth) c c**** print eiaz1 (cos(phi) azimuth dependent term) c do 100 i=1,imu c if (jprint(8) .eq. 1) write(33,2100) thta(i), * (eiaz1(i,j),j=1,iazmth) 100 continue c c**** print eiaz2 heading c if (jprint(8) .eq. 1) write(33,2200) thenot, * (azmth(j),j=1,iazmth) c c**** print eiaz2 (cos(2*phi) dependent term) c do 200 i=1,imu c if (jprint(8) .eq. 1) write(33,2100) thta(i), * (eiaz2(i,j),j=1,iazmth) 200 continue c c**** print istar total ground reflected radiation at top c**** of atmosphere. c c inc=nalb-1 do 400 k=1,nalb if (jprint(8) .eq. 1) write(33,3200) alb(k) if (jprint(8) .eq. 1) write(33,3300) (ristar(i,k),i=1,imu) c c **** print total intensity izero+i1+i2+istar at top c **** of atmosphere for each albedo and polar look angle. c if (jprint(8) .eq. 1) write(33,3400) thenot,alb(k), * (azmth(m),m=1,iazmth) do 500 i=1,imu if (jprint(8) .eq. 1) write(33,2100) thta(i), * (total(i,j,k),j=1,iazmth) 500 continue 400 continue return 1000 format (1h1,t50,'thnot= ',f6.2,///,1x,'sb(4)= ', 1 4x,(d12.6,1x),//,1x,'qr(11)= ',/, 2 2(13x,5(d12.6,1x),d12.6,/)) c 1100 format (5x,'theta ',5x,' tenstr',//,(1x,i2,3x,f6.2,2x,d12.6)) 1200 format(///,1x,'azimuth angle independent data for thenot= ', 1 f6.2,//,1x,'fs= ',d12.6,5x,'gg= ',d12.6,//, 2 5x,'theta',5x,'eizero',/, 3 (1x,i2,2x,f6.2,2x,d12.6)) c 2000 format(//,11x,'theta/phi table for eiaz1 term', 1 10x,'thenot= ',f6.2,//, 2 1x,'theta',/,10x,8(f6.2,7x)) c 2100 format(1x,f6.2,8(1x,d12.6)) c 2200 format(//,11x,'theta/phi table for eiaz2 term', 1 10x,'thenot= ',f6.2,//, 2 1x,'theta',/,10x,8(f6.2,7x)) c 3000 format(/,' albedo= ',11(f4.2,1x)) c c 3200 format(/,' albedo= ',f4.2) c 3300 format(/,' istar-- total reflected intensity for ', 1 'current albedo and polar look angles',/, 2 ' ristar(15,11)= ',/,3(15x,5(d12.6,1x),/)) 3400 format(//,11x,'theta/phi table for total intensity', 1 ' for thenot= ',f6.2,' and albedo= ',f4.2,//, 2 1x,'theta ',/,10x,8(f6.2,7x)) end