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) integer*4 nwavl, layr 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 p(101),t(101),tl(101),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 alphp(12), betp(12), xsav(11), xdel real*8 psza(10) character file_name*200 C C Functions C logical get_coef C C Common C include "cwrite_iter.inc" include "emm.inc" include "kmat.inc" include "log.inc" include "consts.inc" include "thkns.inc" include "atmos.inc" include "input.inc" include "inchr.inc" include "in.inc" include "contrl.inc" include "out.inc" include "totals.inc" include "prints.inc" include "buff1.inc" include "buff2.inc" include "buff3.inc" include "buff4.inc" include "csphout.inc" include "cgcorrect.inc" include "crefdir.inc" include "czfunc.inc" include "depolt.inc" c c if lsphout = True, then perform the out going beam using c spherical geometry. If False, then use flat atmosphere. c data lsphout/.TRUE./ c c if lgcorrect = True, then perform the gravity correciton to the c rayleigh scattering optical depth. c data lgcorrect/.TRUE./ c 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, wavlen, 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) = wavlen if (jjprint.ne.0) write (33,7902) wavel(num) lambda = wavlen if (jprint(9).eq.1) write (6,6200) wavlen do while (iter_pos .lt. num_iter .and. * wave_iter(iter_pos+1) .lt. wavlen) iter_pos = iter_pos + 1 end do itmax = iter(iter_pos) c new statements if (jjprint .ne. 0) write(33,2121) wavlen 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) wavlen, 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), * ttlp(j), ttrp(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(' wavlen',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