example.txt 100644 000371 000325 00000016022 06541271015 014135 0 ustar 00russd stennis 000000 000000 Format read from file sounding.txt is (14X,F7.1,F7.3) Original sounding plus ozone: lev pres temp wvmr ozmr 1 .1 239.40 .003 .6164 2 .2 251.20 .003 1.0800 3 .5 264.50 .003 2.0939 4 1.0 267.80 .003 3.5141 5 1.5 263.90 .003 4.7302 6 2.0 261.20 .003 5.7203 7 3.0 254.20 .003 6.8769 8 4.0 248.60 .003 7.3300 9 5.0 244.20 .003 7.4627 10 7.0 238.90 .003 7.2790 11 10.0 233.80 .003 6.6234 12 15.0 229.20 .003 5.8500 13 20.0 226.20 .003 5.3962 14 25.0 224.20 .003 4.9888 15 30.0 222.80 .003 4.5007 16 50.0 220.30 .003 3.0250 17 60.0 219.70 .003 2.3849 18 70.0 219.40 .003 1.8667 19 85.0 219.20 .003 1.3210 20 100.0 219.00 .003 1.0185 21 115.0 218.90 .005 .8524 22 135.0 218.90 .008 .7173 23 150.0 218.40 .011 .6449 24 200.0 216.40 .025 .3857 25 250.0 220.00 .050 .2368 26 300.0 227.70 .086 .1484 27 350.0 236.00 .164 .1004 28 400.0 243.30 .214 .0789 29 430.0 247.30 .300 .0691 30 475.0 252.60 .479 .0577 31 500.0 255.30 .619 .0537 32 570.0 261.60 1.063 .0441 33 620.0 265.40 1.531 .0389 34 670.0 268.60 2.124 .0361 35 700.0 270.40 2.529 .0345 36 780.0 274.20 3.614 .0310 37 850.0 277.00 4.549 .0294 38 920.0 279.30 5.046 .0282 39 950.0 279.50 5.293 .0278 40 1000.0 277.40 5.144 .0270 Transmittances chan 1 2 3 4 5 6 7 8 9 10 11 12 1 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 2 .9940 .9955 .9975 .9973 .9997 .9998 1.0000 1.0000 .9999 .9997 .9992 1.0000 3 .9870 .9901 .9946 .9941 .9989 .9995 .9999 1.0000 .9998 .9990 .9977 1.0000 4 .9767 .9823 .9911 .9907 .9978 .9990 .9998 1.0000 .9995 .9982 .9959 1.0000 5 .9671 .9751 .9880 .9876 .9969 .9985 .9997 1.0000 .9992 .9977 .9945 1.0000 6 .9576 .9682 .9849 .9845 .9960 .9979 .9996 .9999 .9990 .9974 .9932 1.0000 7 .9390 .9550 .9790 .9785 .9943 .9968 .9994 .9999 .9986 .9967 .9909 1.0000 8 .9207 .9422 .9732 .9726 .9928 .9955 .9991 .9998 .9983 .9962 .9885 1.0000 9 .9026 .9298 .9677 .9670 .9913 .9942 .9989 .9998 .9980 .9957 .9862 .9999 10 .8674 .9059 .9570 .9562 .9885 .9916 .9985 .9997 .9974 .9947 .9817 .9999 11 .8165 .8718 .9418 .9412 .9847 .9879 .9979 .9997 .9967 .9932 .9753 .9998 12 .7369 .8183 .9180 .9187 .9791 .9818 .9971 .9996 .9954 .9906 .9650 .9996 13 .6636 .7682 .8954 .8984 .9738 .9757 .9964 .9995 .9942 .9880 .9548 .9993 14 .5964 .7208 .8737 .8799 .9689 .9697 .9958 .9994 .9931 .9853 .9448 .9991 15 .5351 .6761 .8529 .8628 .9643 .9636 .9952 .9993 .9919 .9826 .9347 .9988 16 .3413 .5195 .7751 .8041 .9477 .9372 .9933 .9991 .9871 .9714 .8929 .9975 17 .2693 .4541 .7395 .7792 .9404 .9229 .9926 .9990 .9846 .9658 .8712 .9968 18 .2102 .3962 .7056 .7563 .9336 .9074 .9919 .9990 .9821 .9602 .8486 .9959 19 .1416 .3221 .6575 .7248 .9241 .8829 .9911 .9989 .9783 .9516 .8140 .9945 20 .0926 .2614 .6122 .6958 .9152 .8576 .9903 .9988 .9744 .9431 .7792 .9929 21 .0590 .2121 .5696 .6685 .9066 .8323 .9896 .9987 .9700 .9319 .7452 .9912 22 .0315 .1608 .5165 .6342 .8956 .7993 .9887 .9986 .9626 .9105 .7016 .9888 23 .0197 .1310 .4796 .6100 .8875 .7749 .9880 .9985 .9557 .8901 .6701 .9869 24 .0051 .0684 .3740 .5373 .8624 .6963 .9860 .9983 .9255 .8003 .5722 .9800 25 .0020 .0369 .2881 .4722 .8377 .6223 .9838 .9981 .8820 .6812 .4849 .9723 26 .0009 .0201 .2151 .4093 .8105 .5544 .9809 .9978 .8247 .5476 .4083 .9641 27 .0004 .0110 .1542 .3479 .7789 .4918 .9765 .9972 .7430 .3975 .3415 .9553 28 .0002 .0061 .1063 .2906 .7437 .4346 .9703 .9964 .6513 .2674 .2840 .9460 29 .0001 .0044 .0836 .2590 .7209 .4028 .9654 .9958 .5897 .2001 .2538 .9402 30 .0000 .0027 .0571 .2164 .6849 .3587 .9555 .9943 .4830 .1119 .2141 .9312 31 .0000 .0021 .0457 .1952 .6640 .3358 .9486 .9930 .4183 .0743 .1946 .9260 32 .0000 .0010 .0240 .1450 .6033 .2781 .9261 .9879 .2468 .0178 .1487 .9108 33 .0000 .0006 .0148 .1162 .5574 .2419 .9024 .9813 .1429 .0046 .1224 .8992 34 .0000 .0004 .0091 .0922 .5094 .2096 .8718 .9711 .0692 .0009 .1007 .8870 35 .0000 .0003 .0067 .0796 .4794 .1919 .8495 .9624 .0409 .0003 .0895 .8794 36 .0000 .0001 .0030 .0524 .3971 .1507 .7753 .9287 .0077 .0001 .0651 .8578 37 .0000 .0001 .0015 .0347 .3237 .1211 .6917 .8836 .0013 .0000 .0491 .8375 38 .0000 .0000 .0007 .0221 .2560 .0968 .6008 .8283 .0002 .0000 .0369 .8162 39 .0000 .0000 .0005 .0181 .2294 .0878 .5608 .8018 .0002 .0000 .0327 .8068 40 .0000 .0000 .0003 .0128 .1906 .0747 .4978 .7571 .0001 .0000 .0266 .7911 RAD 51.716 48.676 49.779 58.510 82.023 .443 94.629 81.652 11.689 4.585 .270 .351 mW/m^2-str-cm^-1 TBB 227.89 225.85 228.39 238.05 261.22 253.07 274.20 276.75 251.17 236.30 245.11 275.37 degrees Kelvin Sounding retrieved from VAS TBB: lev pres temp wvmr 1 .1 240.75 .003 2 .2 251.94 .003 3 .5 264.50 .003 4 1.0 267.47 .003 5 1.5 263.62 .003 6 2.0 260.84 .003 7 3.0 253.75 .003 8 4.0 248.05 .003 9 5.0 243.11 .003 10 7.0 237.33 .003 11 10.0 232.42 .003 12 15.0 227.75 .003 13 20.0 225.12 .003 14 25.0 223.74 .003 15 30.0 222.77 .003 16 50.0 220.23 .003 17 60.0 219.53 .003 18 70.0 219.08 .003 19 85.0 218.75 .003 20 100.0 219.08 .003 21 115.0 219.63 .005 22 135.0 220.46 .007 23 150.0 220.72 .009 24 200.0 218.17 .022 25 250.0 219.87 .043 26 300.0 226.54 .075 27 350.0 234.31 .175 28 400.0 241.41 .289 29 430.0 245.31 .399 30 475.0 250.66 .586 31 500.0 253.43 .712 32 570.0 260.01 1.110 33 620.0 264.07 1.403 34 670.0 267.64 1.692 35 700.0 269.51 1.867 36 780.0 273.76 2.379 37 850.0 276.60 2.962 38 920.0 278.93 3.516 39 950.0 279.47 3.766 40 1000.0 278.41 4.102 VASFRWRD -- normal termination sounding.txt 100644 000371 000325 00000002723 06533035307 014336 0 ustar 00russd stennis 000000 000000 (14X,F7.1,F7.3) 1 0. 239.4 .003 65670. 2 0. 251.2 .003 60691. 3 0. 264.5 .003 53773. 4 1. 267.8 .003 48371. 5 2. 263.9 .003 45214. 6 2. 261.2 .003 43003. 7 3. 254.2 .003 39943. 8 4. 248.6 .003 37824. 9 5. 244.2 .003 36214. 10 7. 238.9 .003 33834. 11 10. 233.8 .003 31364. 12 15. 229.2 .003 28615. 13 20. 226.2 .003 26697. 14 25. 224.2 .003 25225. 15 30. 222.8 .003 24032. 16 50. 220.3 .003 20718. 17 60. 219.7 .003 19543. 18 70. 219.4 .003 18552. 19 85. 219.2 .003 17305. 20 100. 219.0 .003 16263. 21 115. 218.9 .005 15366. 22 135. 218.9 .008 14339. 23 150. 218.4 .011 13664. 24 200. 216.4 .025 11833. 25 250. 220.0 .050 10407. 26 300. 227.7 .086 9211. 27 350. 236.0 .164 8165. 28 400. 243.3 .214 7227. 29 430. 247.3 .300 6708. 30 475. 252.6 .479 5979. 31 500. 255.3 .619 5598. 32 570. 261.6 1.063 4605. 33 620. 265.4 1.531 3956. 34 670. 268.6 2.124 3349. 35 700. 270.4 2.529 3002. 36 780. 274.2 3.614 2138. 37 850. 277.0 4.549 1442. 38 920. 279.3 5.046 796. 39 950. 279.5 5.293 533. 40 1000. 277.4 5.144 113. zen= 48.00 month 5 satellite 7 39 950.0 279.50 5.293 .0278 40 1000vasfrwrd.f 100644 000371 000325 00000133037 06533035307 013757 0 ustar 00russd stennis 000000 000000 program vasfrwrd c * Forward calculation of VAS radiances from a sounding. c Printed output includes the input sounding, transmittances c for the 12 VAS channels, calculated radiances and blackbody c temperatures, and the regression-retrieved sounding. parameter (iui=51,iuo=61,nc=12,nl=40) character*16 cfmt character*12 ifile,ofile dimension taut(nl),temp(nl),wvmr(nl),ozmr(nl),p(nl), + tauout(nc,nl),tvas(nc),rvas(nc) dimension rtemp(nl),rwvmr(nl) data p/0.1,0.2,0.5,1.0,1.5,2.,3.,4.,5.,7.,10.,15.,20.,25.,30., + 50.,60.,70.,85.,100.,115.,135.,150.,200.,250.,300.,350.,400., + 430.,475.,500.,570.,620.,670.,700.,780.,850.,920.,950.,1000./ c * Specified parameters ksat=7 zena=48.0 kmon=3 rlat=45.0 c * Input file ifile='sounding.txt' c * Output file ofile='example.out' c * Get the climatological ozone mixing ratio. call clmozo(rlat,kmon, ozmr) c * Open the sounding file and read the format from the first line. open(iui,file=ifile,form='formatted',status='old',err=900) c * Open the output (text) file open(iuo,file=ofile,form='formatted',status='unknown') read(iui,'(a16)') cfmt write(iuo,'(''Format read from file '',a12,'' is '',a16)') + ifile,cfmt write(iuo,'(''Original sounding plus ozone:''// + ''lev pres temp wvmr ozmr''/)') numb=0 c * Read in the sounding, print temp, wvmr plus ozmr from climatology do lev=1,nl read(iui,cfmt,err=100,end=100) temp(lev),wvmr(lev) write(iuo,'(i3,f8.1,f8.2,f8.3,f8.4)') lev,p(lev),temp(lev), + wvmr(lev),ozmr(lev) numb=numb+1 enddo 100 close(iui) if(numb.lt.nl) go to 901 c * Input the Planck coefficients for the specified satellite. call pfcvas(ksat,*902) ls=nl c .. ls is the surface level. tsk=temp(ls) c .. tsk is the surface temperature. do kan=1,nc c ..... calculate transmittances for each channel, store for printing call tauvas(temp,wvmr,ozmr,zena,ksat,kan,taut,*903) do lev=1,nl tauout(kan,lev)=taut(lev) enddo c ..... obtain blackbody temperature tvas(kan)=tbcvas(taut,temp,tsk,ls,kan) c ..... convert back to radiance rvas(kan)=plavas(tvas(kan),kan) enddo write(iuo,'(/''Transmittances''//''chan'',12(4x,i4)/)') + (k,k=1,nc) do l=1,nl write(iuo,'(i4,12f8.4)') l,(tauout(k,l),k=1,nc) enddo write(iuo,'(/''RAD '',12f8.3)') (rvas(kan),kan=1,nc) write(iuo,'('' mW/m^2-str-cm^-1'')') write(iuo,'(/''TBB '',12f8.2)') (tvas(kan),kan=1,nc) write(iuo,'('' degrees Kelvin'')') c * Use the blackbody temperatures to retrieve c a sounding via regression. call vasrtw(tvas,zena,kmon,ksat,rtemp,rwvmr,*904) write(iuo,'(/''Sounding retrieved from VAS TBB:''// + ''lev pres temp wvmr''/)') do lev=1,nl write(iuo,'(i3,f8.1,f8.2,f8.3)') lev,p(lev), + rtemp(lev),rwvmr(lev) enddo write(iuo,'(/''VASFRWRD -- normal termination'')') write(0,'(''VASFRWRD -- normal termination'')') go to 999 900 write(iuo,'(/''OPEN error, sounding file '',a12)') ifile go to 990 901 write(iuo,'(/''Did not get the complete sounding.'')') go to 990 902 write(iuo,'(/''Trouble (probably I/O) in PFCVAS'')') go to 990 903 write(iuo,'(/''Trouble (probably I/O) in TAUVAS'')') go to 990 904 write(iuo,'(/''Trouble (probably I/O) in VASRTW'')') 990 write(iuo,'(/''VASFRWRD -- abnormal termination'')') write(0,'(''VASFRWRD -- abnormal termination'')') 999 continue close(iuo) end function tbcvas(tau,tem,tsk,ls,kan) c * Integrate radiative-transfer equation for GOES/VAS c .... version of 10.12.97 c tau = transmittance profile c tem = temperature profile (degK) c tsk = surface-skin temperature (degK) c ls = surface model-level (.le.40) c kan = channel number dimension tau(*),tem(*) rad=0. b1=plavas(tem(1),kan) tau1=tau(1) do i=2,ls b2=plavas(tem(i),kan) tau2=tau(i) rad=rad+.5*(b1+b2)*(tau1-tau2) b1=b2 tau1=tau2 enddo bs=plavas(tsk,kan) rad=rad+bs*tau(ls) tbcvas=brivas(rad,kan) return end function plavas(t,k) c * Planck function for VAS c .... version of 17.02.98 c * Call 'PFCVAS' before using this or 'BRIVAS' c t = temperature, deg K c k = channel number: 1,...,12 parameter (nb=12,nt=2) common/vaspfc/wnum(nb),fk1(nb),fk2(nb),tc(nt,nb) tt=tc(1,k)+tc(2,k)*t expn=exp(fk2(k)/tt) - 1. plavas=fk1(k)/expn return end function brivas(r,k) c * Inverse Planck function for VAS c .... version of 17.02.98 c * Call 'PFCVAS' before using this or 'PLAVAS' c r = radiance, standard units c k = channel number: 1,...,12 parameter (nb=12,nt=2) common/vaspfc/wnum(nb),fk1(nb),fk2(nb),tc(nt,nb) expn=fk1(k)/r+1. tt=fk2(k)/alog(expn) brivas=(tt-tc(1,k))/tc(2,k) return end subroutine pfcvas(ngoes,*) c * Input VAS Planck-function and band-correction coefficients. c .... version of 18.02.98 c * Valid for GOES-4, -5, -6, -7 parameter (iuc=33,lenc=100,nb=12,nt=2,lenp=nb*(nt+3)) parameter (lencb=lenc*4) c common/vaspfc/wnum(nb),fk1(nb),fk2(nb),tc(nt,nb) common/vaspfc/pbuf(lenp) dimension cbuf(lenc) character*12 cfile/'vasxtbnd.dat'/ open(iuc,file=cfile,recl=lencb,access='direct',status='old', * err=100) irec=ngoes-3 read(iuc,rec=irec) cbuf close(iuc) do i=1,lenp pbuf(i)=cbuf(i) enddo return 100 return1 end subroutine tauvas(temp,wvmr,ozmr,theta,ngoes,kan,taut,*) c ... Calculate VAS transmittances c ... version of 18.02.98 c ... Valid for GOES-4, -5, -6, -7 c * Inputs: c temp = temperature profile (degK) c wvmr = water-vapor mixing-ratio profile (g/kg) c ozmr = ozone mixing-ratio profile (ppmv) c theta = local zenith angle (deg) c ngoes = satellite number (4...7) c kan = channel number (1...12) c * Outputs: c taut = total transmittance profile c * = error return for I/O problems c * LarrabeeStrow/HalWoolf/PaulVanDelst regression model based on c * 'FSCD3P' line-by-line transmittances. c * Input temperatures, and water-vapor and ozone mixing ratios, must c * be defined at the pressure levels in array 'pstd' c * (see block data 'screfatm'). c * Units: temperature, deg-K; water vapor, g/kg; ozone, ppmv. c * Logical unit 33 is used for coefficient files. c * Component tau's are returned through common, product in 'taut'. parameter (iuc=33,nk=5,nl=40,nm=nl-1,nr=12) parameter (nxc= 4,ncc=nxc+1,lencc=ncc*nm,lenccb=lencc*4) parameter (nxd= 8,ncd=nxd+1,lencd=ncd*nm,lencdb=lencd*4) parameter (nxo= 9,nco=nxo+1,lenco=nco*nm,lencob=lenco*4) parameter (nxl= 2,ncl=nxl+1,lencl=ncl*nm,lenclb=lencl*4) parameter (nxs=11,ncs=nxs+1,lencs=ncs*nm,lencsb=lencs*4) parameter (nxw=nxl+nxs) common/atmstd/pstd(nl),tstd(nl),wstd(nl),ostd(nl) common/taudwo/taud(nl),tauw(nl),tauo(nl) dimension temp(*),wvmr(*),ozmr(*),taut(*) dimension coefd(ncd,nm,nr),coefo(nco,nm,nr),coefl(ncl,nm,nr) dimension coefs(ncs,nm,nr),coefc(ncc,nm,nr) dimension pavg(nm),tref(nm),wref(nm),oref(nm) dimension tavg(nm),wamt(nm),oamt(nm),secz(nm) dimension tauc(nl),tlas(nl),wlas(nl),olas(nl) dimension xdry(nxd,nm),xozo(nxo,nm),xwet(nxw,nm),xcon(nxc,nm) character*12 cfile/'vasxtcom.dat'/ character*3 comp(nk)/'dry','ozo','wts','wtl','wco'/ integer*4 lengcf(nk)/lencdb,lencob,lencsb,lenclb,lenccb/ logical newang,newatm data init/0/,tlas/nl*0./,wlas/nl*0./,olas/nl*0./,zlas/-999./ secant(z)=1./cos(0.01745329*z) if(ngoes.ne.init) then koff=(ngoes-4)*nr do l=1,nk cfile(6:8)=comp(l) lencf=lengcf(l) open(iuc,file=cfile,recl=lencf,access='direct', * status='old',err=200) if(l.eq.1) then do k=1,nr krec=k+koff read(iuc,rec=krec) ((coefd(i,j,k),i=1,ncd),j=1,nm) enddo elseif(l.eq.2) then do k=1,nr krec=k+koff read(iuc,rec=krec) ((coefo(i,j,k),i=1,nco),j=1,nm) enddo elseif(l.eq.3) then do k=1,nr krec=k+koff read(iuc,rec=krec) ((coefs(i,j,k),i=1,ncs),j=1,nm) enddo elseif(l.eq.4) then do k=1,nr krec=k+koff read(iuc,rec=krec) ((coefl(i,j,k),i=1,ncl),j=1,nm) enddo elseif(l.eq.5) then do k=1,nr krec=k+koff read(iuc,rec=krec) ((coefc(i,j,k),i=1,ncc),j=1,nm) enddo endif close(iuc) enddo call conpir(pstd,tstd,wstd,ostd,nl,1,pavg,tref,wref,oref) init=ngoes endif dt=0. dw=0. do=0. do j=1,nl dt=dt+abs(temp(j)-tlas(j)) tlas(j)=temp(j) dw=dw+abs(wvmr(j)-wlas(j)) wlas(j)=wvmr(j) do=do+abs(ozmr(j)-olas(j)) olas(j)=ozmr(j) taud(j)=1.0 tauw(j)=1.0 tauc(j)=1.0 tauo(j)=1.0 taut(j)=1.0 enddo datm=dt+dw+do newatm=datm.ne.0. if(newatm) then call conpir(pstd,temp,wvmr,ozmr,nl,1,pavg,tavg,wamt,oamt) endif newang=theta.ne.zlas if(newang) then zsec=secant(theta) do l=1,nm secz(l)=zsec enddo zlas=theta endif if(newang.or.newatm) then call calpir(tref,wref,oref,tavg,wamt,oamt,pavg,secz, * nm,nxd,nxw,nxo,nxc,xdry,xwet,xozo,xcon) endif k=kan c * dry call taugen(ncd,nxd,nm,coefd(1,1,k),xdry,taud) c * ozo call taugen(nco,nxo,nm,coefo(1,1,k),xozo,tauo) c * wet call tauwet(ncs,ncl,nxs,nxl,nxw,nm,coefs(1,1,k), * coefl(1,1,k),xwet,tauw) call taugen(ncc,nxc,nm,coefc(1,1,k),xcon,tauc) do j=1,nl tauw(j)=tauw(j)*tauc(j) enddo c * total do j=1,nl taut(j)=taud(j)*tauo(j)*tauw(j) enddo return 200 return1 end block data screfatm c * Reference Atmosphere is 1976 U.S. Standard parameter (nl=40) common/atmstd/pstd(nl),tstd(nl),wstd(nl),ostd(nl) data pstd/ .1,.2,.5,1.,1.5,2.,3.,4.,5.,7.,10.,15.,20.,25.,30., + 50.,60.,70.,85.,100.,115.,135.,150.,200.,250.,300.,350.,400., + 430.,475.,500.,570.,620.,670.,700.,780.,850.,920.,950.,1000./ data tstd/ + 231.70, 245.22, 263.35, 270.63, 264.07, 257.93, 249.51, 243.65, + 239.24, 232.64, 228.07, 225.00, 223.13, 221.72, 220.54, 217.28, + 216.70, 216.70, 216.70, 216.70, 216.70, 216.70, 216.70, 216.72, + 220.85, 228.58, 235.38, 241.45, 244.81, 249.48, 251.95, 258.32, + 262.48, 266.40, 268.61, 274.21, 278.74, 282.97, 284.71, 287.50/ data wstd/ + 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, + 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, 0.003, + 0.003, 0.003, 0.003, 0.003, 0.003, 0.004, 0.005, 0.014, + 0.036, 0.089, 0.212, 0.331, 0.427, 0.588, 0.699, 1.059, + 1.368, 1.752, 1.969, 2.741, 3.366, 3.976, 4.255, 4.701/ data ostd/ + 0.65318,1.04797,2.13548,3.82386,5.26768,6.11313,7.35964,7.75004, + 7.82119,7.56126,6.92006,6.10266,5.55513,5.15298,4.59906,2.86792, + 2.29259,1.80627,1.28988,0.93973,0.72277,0.54848,0.46009,0.29116, + 0.16277,0.09861,0.06369,0.05193,0.04718,0.04097,0.03966,0.03614, + 0.03384,0.03342,0.03319,0.03249,0.03070,0.02878,0.02805,0.02689/ end subroutine calpir(t_avg_ref, amt_wet_ref, amt_ozo_ref, + t_avg, amt_wet, amt_ozo, + p_avg, sec_theta, n_layers, + n_dry_pred, n_wet_pred, n_ozo_pred, + n_con_pred, + pred_dry, pred_wet, pred_ozo, + pred_con) c ... version of 19.09.96 c PURPOSE: c Routine to calculate the predictors for the dry (temperature), c wet and ozone components of a fast transmittance model for a c scanning satellite based instrument. c REFERENCES: c AIRS FTC package science notes and software, S. Hannon and L. Strow, c Uni. of Maryland, Baltimore County (UMBC) c CREATED: c 19-Sep-1996 HMW c ARGUMENTS: c Input c ----------- c t_avg_ref - REAL*4 reference layer average temperature array (K) c amt_wet_ref - REAL*4 reference water vapour amount array (k.mol)/cm^2 c amt_ozo_ref - REAL*4 reference ozone amount array (k.mol)/cm^2 c t_avg - REAL*4 layer average temperature array (K) c amt_wet - REAL*4 water vapour amount array (k.mol)/cm^2 c amt_ozo - REAL*4 ozone amount array (k.mol)/cm^2 c p_avg - REAL*4 layer average pressure array (mb) c sec_theta - REAL*4 secant of the zenith angle array c n_layers - INT*4 Number of atmospheric layers c n_dry_pred - INT*4 number of dry (temperature) predictors c n_wet_pred - INT*4 number of water vapour predictors c n_ozo_pred - INT*4 number of ozone predictors c n_con_pred - INT*4 number of water vapour continuum predictors c Output c ----------- c pred_dry - REAL*4 dry gas (temperature) predictor matrix c pred_wet - REAL*4 water vapour predictor matrix c pred_ozo - REAL*4 ozone predictor matrix c pred_con - REAL*4 water vapour continuum predictor matrix c COMMENTS: c Levels or Layers? c ----------------- c Profile data is input at a number of *LAYERS*. c Layer Numbering pt. A c --------------------- c Layer 1 => Atmosphere between LEVELs 1 & 2 c Layer 2 => Atmosphere between LEVELs 2 & 3 c . c . c . c Layer L-1 => Atmosphere between LEVELs L-1 & L c Layer Numbering pt. B c --------------------- c For the HIS instrument, Layer 1 is at the top of the atmosphere c and Layer L-1 is at the surface. c Layer Numbering pt. C c --------------------- c In this routine the number of *LAYERS* is passed in the argument c list, _not_ the number of LEVELS. This was done to improve c the readability of this code, i.e. loop from 1->L(ayers) c rather than from 1->L(evels)-1. c======================================================================= c----------------------------------------------------------------------- c Turn off implicit type declaration c----------------------------------------------------------------------- implicit none c------------------------------------------------------------------------ c Arguments c------------------------------------------------------------------------ c -- Input integer*4 n_layers, + n_dry_pred, n_wet_pred, n_ozo_pred, n_con_pred real*4 t_avg_ref(*), amt_wet_ref(*), amt_ozo_ref(*), + t_avg(*), amt_wet(*), amt_ozo(*), + p_avg(*), sec_theta(*) c -- Output real*4 pred_dry(n_dry_pred, *), + pred_wet(n_wet_pred, *), + pred_ozo(n_ozo_pred, *), + pred_con(n_con_pred, *) c------------------------------------------------------------------------ c Local variables c------------------------------------------------------------------------ c -- Parameters integer*4 max_layers parameter ( max_layers = 41 ) integer*4 max_dry_pred, max_wet_pred, max_ozo_pred, max_con_pred parameter ( max_dry_pred = 8, + max_wet_pred = 13, + max_ozo_pred = 9, + max_con_pred = 4 ) c -- Scalars integer*4 l c -- Arrays c ....Pressure real*4 p_dp(max_layers), + p_norm(max_layers) c ....Temperature real*4 delta_t(max_layers), + t_ratio(max_layers), + pw_t_ratio(max_layers) ! Pressure weighted c ....Water vapour real*4 wet_ratio(max_layers), + pw_wet(max_layers), ! Pressure weighted + pw_wet_ref(max_layers), ! Pressure weighted + pw_wet_ratio(max_layers) ! Pressure weighted c ....Ozone real*4 ozo_ratio(max_layers), + pw_ozo_ratio(max_layers), ! Pressure weighted + pow_t_ratio(max_layers) ! Pressure/ozone weighted c************************************************************************ c ** Executable code ** c************************************************************************ c------------------------------------------------------------------------ c -- Check that n_layers is o.k. -- c------------------------------------------------------------------------ if( n_layers .gt. max_layers )then write(*,'(/10x,''*** calpir : n_layers > max_layers'')') stop end if c------------------------------------------------------------------------ c -- Check that numbers of predictors is consistent -- c------------------------------------------------------------------------ c --------------------------------- c # of dry (temperature) predictors c --------------------------------- if( n_dry_pred .ne. max_dry_pred )then write(*,'(/10x,''*** calpir : invalid n_dry_pred'')') stop end if c ---------------------------- c # of water vapour predictors c ---------------------------- if( n_wet_pred .ne. max_wet_pred )then write(*,'(/10x,''*** calpir : invalid n_wet_pred'')') stop end if c --------------------- c # of ozone predictors c --------------------- if( n_ozo_pred .ne. max_ozo_pred )then write(*,'(/10x,''*** calpir : invalid n_ozo_pred'')') stop end if c -------------------------------------- c # of water vapour continuum predictors c -------------------------------------- if( n_con_pred .ne. max_con_pred )then write(*,'(/10x,''*** calpir : invalid n_con_pred'')') stop end if c------------------------------------------------------------------------ c -- Calculate ratios, offsets, etc, for top layer -- c------------------------------------------------------------------------ c ------------------ c Pressure variables c ------------------ p_dp(1) = p_avg(1) * ( p_avg(2) - p_avg(1) ) p_norm(1) = 0.0 c --------------------- c Temperature variables c --------------------- delta_t(1) = t_avg(1) - t_avg_ref(1) t_ratio(1) = t_avg(1) / t_avg_ref(1) pw_t_ratio(1) = 0.0 c ---------------- c Amount variables c ---------------- c ....Water vapour wet_ratio(1) = amt_wet(1) / amt_wet_ref(1) pw_wet(1) = p_dp(1) * amt_wet(1) pw_wet_ref(1) = p_dp(1) * amt_wet_ref(1) pw_wet_ratio(1) = wet_ratio(1) c ....Ozone ozo_ratio(1) = amt_ozo(1) / amt_ozo_ref(1) pw_ozo_ratio(1) = 0.0 pow_t_ratio(1) = 0.0 c------------------------------------------------------------------------ c -- Calculate ratios, offsets, etc, for all layers -- c------------------------------------------------------------------------ do l = 2, n_layers c ------------------ c Pressure variables c ------------------ p_dp(l) = p_avg(l) * ( p_avg(l) - p_avg(l-1) ) p_norm(l) = p_norm(l-1) + p_dp(l) c --------------------- c Temperature variables c --------------------- delta_t(l) = t_avg(l) - t_avg_ref(l) t_ratio(l) = t_avg(l) / t_avg_ref(l) pw_t_ratio(l) = pw_t_ratio(l-1) + ( p_dp(l) * t_ratio(l-1) ) c ---------------- c Amount variables c ---------------- c ..Water vapour wet_ratio(l) = amt_wet(l) / amt_wet_ref(l) pw_wet(l) = pw_wet(l-1) + ( p_dp(l) * amt_wet(l) ) pw_wet_ref(l) = pw_wet_ref(l-1) + ( p_dp(l) * amt_wet_ref(l) ) c ..Ozone ozo_ratio(l) = amt_ozo(l) / amt_ozo_ref(l) pw_ozo_ratio(l) = pw_ozo_ratio(l-1) + + ( p_dp(l) * ozo_ratio(l-1) ) pow_t_ratio(l) = pow_t_ratio(l-1) + + ( p_dp(l) * ozo_ratio(l-1) * delta_t(l-1) ) end do c------------------------------------------------------------------------ c -- Scale the pressure dependent variables -- c------------------------------------------------------------------------ do l = 2, n_layers pw_t_ratio(l) = pw_t_ratio(l) / p_norm(l) pw_wet_ratio(l) = pw_wet(l) / pw_wet_ref(l) pw_ozo_ratio(l) = pw_ozo_ratio(l) / p_norm(l) pow_t_ratio(l) = pow_t_ratio(l) / p_norm(l) end do c------------------------------------------------------------------------ c -- Load up predictor arrays -- c------------------------------------------------------------------------ do l = 1, n_layers c ---------------------- c Temperature predictors c ---------------------- pred_dry(1,l) = sec_theta(l) pred_dry(2,l) = sec_theta(l) * sec_theta(l) pred_dry(3,l) = sec_theta(l) * t_ratio(l) pred_dry(4,l) = pred_dry(3,l) * t_ratio(l) pred_dry(5,l) = t_ratio(l) pred_dry(6,l) = t_ratio(l) * t_ratio(l) pred_dry(7,l) = sec_theta(l) * pw_t_ratio(l) pred_dry(8,l) = pred_dry(7,l) / t_ratio(l) c ----------------------- c Water vapour predictors c ----------------------- pred_wet(1,l) = sec_theta(l) * wet_ratio(l) pred_wet(2,l) = sqrt( pred_wet(1,l) ) pred_wet(3,l) = pred_wet(1,l) * delta_t(l) pred_wet(4,l) = pred_wet(1,l) * pred_wet(1,l) pred_wet(5,l) = abs( delta_t(l) ) * delta_t(l) * pred_wet(1,l) pred_wet(6,l) = pred_wet(1,l) * pred_wet(4,l) pred_wet(7,l) = sec_theta(l) * pw_wet_ratio(l) pred_wet(8,l) = pred_wet(2,l) * delta_t(l) pred_wet(9,l) = sqrt( pred_wet(2,l) ) pred_wet(10,l) = pred_wet(7,l) * pred_wet(7,l) pred_wet(11,l) = sqrt( pred_wet(7,l) ) pred_wet(12,l) = pred_wet(1,l) pred_wet(13,l) = pred_wet(2,l) c ---------------- c Ozone predictors c ---------------- pred_ozo(1,l) = sec_theta(l) * ozo_ratio(l) pred_ozo(2,l) = sqrt( pred_ozo(1,l) ) pred_ozo(3,l) = pred_ozo(1,l) * delta_t(l) pred_ozo(4,l) = pred_ozo(1,l) * pred_ozo(1,l) pred_ozo(5,l) = pred_ozo(2,l) * delta_t(l) pred_ozo(6,l) = sec_theta(l) * pw_ozo_ratio(l) pred_ozo(7,l) = sqrt( pred_ozo(6,l) ) * pred_ozo(1,l) pred_ozo(8,l) = pred_ozo(1,l) * pred_wet(1,l) pred_ozo(9,l) = sec_theta(l) * pow_t_ratio(l) * pred_ozo(1,l) c --------------------------------- c Water vapour continuum predictors c --------------------------------- pred_con(1,l) = sec_theta(l) * wet_ratio(l) / * ( t_ratio(l) * t_ratio(l) ) pred_con(2,l) = pred_con(1,l) * pred_con(1,l) / sec_theta(l) pred_con(3,l) = sec_theta(l) * wet_ratio(l) / t_ratio(l) pred_con(4,l) = pred_con(3,l) * wet_ratio(l) end do return end subroutine conpir( p, t, w, o, n_levels, i_dir, + p_avg, t_avg, w_amt, o_amt) c ... version of 19.09.96 c PURPOSE: c Function to convert atmospheric water vapour (g/kg) and ozone (ppmv) c profiles specified at n_levels layer BOUNDARIES to n_levels-1 c integrated layer amounts of units (k.moles)/cm^2. The average c LAYER pressure and temperature are also returned. c REFERENCES: c AIRS LAYERS package science notes, S. Hannon and L. Strow, Uni. of c Maryland, Baltimore County (UMBC) c CREATED: c 19-Sep-1996 HMW c ARGUMENTS: c Input c -------- c p - REAL*4 pressure array (mb) c t - REAL*4 temperature profile array (K) c w - REAL*4 water vapour profile array (g/kg) c o - REAL*4 ozone profile array (ppmv) c n_levels - INT*4 number of elements used in passed arrays c i_dir - INT*4 direction of increasing layer number c i_dir = +1, Level(1) == p(top) } satellite/AC c Level(n_levels) == p(sfc) } case c i_dir = -1, Level(1) == p(sfc) } ground-based c Level(n_levels) == p(top) } case c Output c -------- c p_avg - REAL*4 average LAYER pressure array (mb) c t_avg - REAL*4 average LAYER temperature (K) c w_amt - REAL*4 integrated LAYER water vapour amount array (k.moles)/cm^2 c o_amt - REAL*4 integrated LAYER ozone amount array (k.moles)/cm^2 c ROUTINES: c Subroutines: c ------------ c gphite - calculates geopotential height given profile data. c Functions: c ---------- c NONE c COMMENTS: c Levels or Layers? c ----------------- c Profile data is input at a number of *LEVELS*. Number densitites c are calculated for *LAYERS* that are bounded by these levels. c So, for L levels there are L-1 layers. c Layer Numbering c --------------- c Layer 1 => Atmosphere between LEVELs 1 & 2 c Layer 2 => Atmosphere between LEVELs 2 & 3 c . c . c . c Layer L-1 => Atmosphere between LEVELs L-1 & L c======================================================================= c----------------------------------------------------------------------- c -- Prevent implicit typing of variables -- c----------------------------------------------------------------------- implicit none c----------------------------------------------------------------------- c -- Arguments -- c----------------------------------------------------------------------- c -- Arrays real*4 p(*), t(*), w(*), o(*), + p_avg(*), t_avg(*), w_amt(*), o_amt(*) c -- Scalars integer*4 n_levels, i_dir c----------------------------------------------------------------------- c -- Local variables -- c----------------------------------------------------------------------- c -- Parameters integer*4 max_levels parameter ( max_levels = 50 ) ! Maximum number of layers real*4 r_equator, r_polar, r_avg parameter ( r_equator = 6.378388e+06, ! Earth radius at equator + r_polar = 6.356911e+06, ! Earth radius at pole + r_avg = 0.5*(r_equator+r_polar) ) real*4 g_sfc parameter ( g_sfc = 9.80665 ) ! Gravity at surface real*4 rho_ref parameter ( rho_ref = 1.2027e-12 ) ! Reference air "density" real*4 mw_dryair, mw_h2o, mw_o3 parameter ( mw_dryair = 28.97, ! Molec. wgt. of dry air (g/mol) + mw_h2o = 18.0152, ! Molec. wgt. of water + mw_o3 = 47.9982 ) ! Molec. wgt. of ozone real*4 R_gas, R_air parameter ( R_gas = 8.3143, ! Ideal gas constant (J/mole/K) + R_air = 0.9975*R_gas ) ! Gas constant for air (worst case) c -- Scalars integer*4 l, l_start, l_end, l_indx real*4 rho1, rho2, p1, p2, w1, w2, o1, o2, z1, z2, + c_avg, g_avg, z_avg, w_avg, o_avg, + dz, dp, r_hgt, wg, og, A, B c -- Arrays real*4 z(max_levels), ! Pressure heights (m) + g(max_levels), ! Acc. due to gravity (m/s/s) + mw_air(max_levels), ! Molec. wgt. of air (g/mol) + rho_air(max_levels), ! air mass density (kg.mol)/m^3 + c(max_levels), ! (kg.mol.K)/(N.m) + w_ppmv(max_levels) ! h2o LEVEL amount (ppmv) c*********************************************************************** c ** Executable code ** c*********************************************************************** c----------------------------------------------------------------------- c -- Calculate initial values of pressure heights -- c----------------------------------------------------------------------- call gphite( p, t, w, 0.0, n_levels, i_dir, z) c----------------------------------------------------------------------- c -- Set loop bounds for direction sensitive calculations -- c -- so loop iterates from surface to the top -- c----------------------------------------------------------------------- if( i_dir .gt. 0 )then c -------------------- c Data stored top down c -------------------- l_start = n_levels l_end = 1 else c --------------------- c Data stored bottom up c --------------------- l_start = 1 l_end = n_levels end if c----------------------------------------------------------------------- c -- Air molecular mass and density, and gravity -- c -- as a function of LEVEL -- c----------------------------------------------------------------------- c ----------------------- c Loop from bottom to top c ----------------------- do l = l_start, l_end, -1*i_dir c --------------------------------- c Convert water vapour g/kg -> ppmv c --------------------------------- w_ppmv(l) = 1.0e+03 * w(l) * mw_dryair / mw_h2o c ----------------------------------------- c Calculate molecular weight of air (g/mol) c ---------------------------------------- mw_air(l) = ( ( 1.0 - (w_ppmv(l)/1.0e+6) ) * mw_dryair ) + + ( ( w_ppmv(l)/1.0e+06 ) * mw_h2o ) c ---------------- c Air mass density c ---------------- c(l) = 0.001 * mw_air(l) / R_air ! 0.001 factor for g -> kg rho_air(l) = c(l) * p(l) / t(l) c ------- c Gravity c ------- r_hgt = r_avg + z(l) ! m g(l) = g_sfc - ! m/s^2 + g_sfc*( 1.0 - ( (r_avg*r_avg)/(r_hgt*r_hgt) ) ) end do c----------------------------------------------------------------------- c -- LAYER quantities -- c----------------------------------------------------------------------- c ----------------------- c Loop from bottom to top c ----------------------- do l = l_start, l_end+i_dir, -1*i_dir c ------------------------------------------------------- c Determine output array index. This is done so that the c output data is always ordered from 1 -> L-1 regardless c of the orientation of the input data. This is true by c default only for the bottom-up case. For the top down c case no correction would give output layers from 2 -> L c ------------------------------------------------------- if( i_dir .gt. 0 )then l_indx = l - 1 else l_indx = l end if c --------------------------------------- c Assign current layer boundary densities c --------------------------------------- rho1 = rho_air(l) rho2 = rho_air(l-i_dir) c --------- c Average c c --------- c_avg = ( (rho1*c(l)) + (rho2*c(l-i_dir)) ) / ( rho1 + rho2 ) c --------- c Average t c --------- t_avg(l_indx) = + ( (rho1*t(l)) + (rho2*t(l-i_dir)) ) / ( rho1 + rho2 ) c --------- c Average p c --------- p1 = p(l) p2 = p(l-i_dir) z1 = z(l) z2 = z(l-i_dir) dp = p2 - p1 A = log(p2/p1) / (z2-z1) B = p1 / exp(A*z1) p_avg(l_indx) = dp / log(p2/p1) c ------------------------------------------------ c LAYER thickness (rather long-winded as it is not c assumed the layers are thin) in m. Includes c correction for altitude/gravity. c ------------------------------------------------ c ...Initial values g_avg = g(l) dz = -1.0 * dp * t_avg(l_indx) / ( g_avg*c_avg*p_avg(l_indx) ) c ...Calculate z_avg z_avg = z(l) + ( 0.5*dz ) c ...Calculate new g_avg r_hgt = r_avg + z_avg g_avg = g_sfc - g_sfc*( 1.0 - ( (r_avg*r_avg)/(r_hgt*r_hgt) ) ) c ...Calculate new dz dz = -1.0 * dp * t_avg(l_indx) / ( g_avg*c_avg*p_avg(l_indx) ) c ---------------------------------------- c Calculate LAYER amounts for water vapour c ---------------------------------------- w1 = w_ppmv(l) w2 = w_ppmv(l-i_dir) w_avg = ( (rho1*w1) + (rho2*w2) ) / ( rho1+rho2 ) w_amt(l_indx) = + rho_ref * w_avg * dz * p_avg(l_indx) / t_avg(l_indx) c --------------------------------- c Calculate LAYER amounts for ozone c --------------------------------- o1 = o(l) o2 = o(l-i_dir) o_avg = ( (rho1*o1) + (rho2*o2) ) / ( rho1+rho2 ) o_amt(l_indx) = + rho_ref * o_avg * dz * p_avg(l_indx) / t_avg(l_indx) end do return end subroutine gphite( p, t, w, z_sfc, n_levels, i_dir, z) c ... version of 19.09.96 c PURPOSE: c Routine to compute geopotential height given the atmospheric state. c Includes virtual temperature adjustment. c CREATED: c 19-Sep-1996 Received from Hal Woolf c ARGUMENTS: c Input c -------- c p - REAL*4 pressure array (mb) c t - REAL*4 temperature profile array (K) c w - REAL*4 water vapour profile array (g/kg) c z_sfc - REAL*4 surface height (m). 0.0 if not known. c n_levels - INT*4 number of elements used in passed arrays c i_dir - INT*4 direction of increasing layer number c i_dir = +1, Level(1) == p(top) } satellite/AC c Level(n_levels) == p(sfc) } case c i_dir = -1, Level(1) == p(sfc) } ground-based c Level(n_levels) == p(top) } case c Output c -------- c z - REAL*4 pressure level height array (m) c COMMENTS: c Dimension of height array may not not be the same as that of the c input profile data. c======================================================================= c----------------------------------------------------------------------- c -- Prevent implicit typing of variables -- c----------------------------------------------------------------------- implicit none c----------------------------------------------------------------------- c -- Arguments -- c----------------------------------------------------------------------- c -- Arrays real*4 p(*), t(*), w(*), + z(*) c -- Scalars integer*4 n_levels, i_dir real*4 z_sfc c----------------------------------------------------------------------- c -- Local variables -- c----------------------------------------------------------------------- c -- Parameters real*4 rog, fac parameter ( rog = 29.2898, + fac = 0.5 * rog ) c -- Scalars integer*4 i_start, i_end, l real*4 v_lower, v_upper, algp_lower, algp_upper, hgt c*********************************************************************** c ** Executable code ** c*********************************************************************** c----------------------------------------------------------------------- c -- Calculate virtual temperature adjustment and exponential -- c -- pressure height for level above surface. Also set integration -- c -- loop bounds -- c----------------------------------------------------------------------- if( i_dir .gt. 0 )then c -------------------- c Data stored top down c -------------------- v_lower = t(n_levels) * ( 1.0 + ( 0.00061 * w(n_levels) ) ) algp_lower = alog( p(n_levels) ) i_start = n_levels-1 i_end = 1 else c --------------------- c Data stored bottom up c --------------------- v_lower = t(1) * ( 1.0 + ( 0.00061 * w(1) ) ) algp_lower = alog( p(1) ) i_start = 2 i_end = n_levels end if c----------------------------------------------------------------------- c -- Assign surface height -- c----------------------------------------------------------------------- hgt = z_sfc c----------------------------------------------------------------------- c -- Loop over layers always from sfc -> top -- c----------------------------------------------------------------------- do l = i_start, i_end, -1*i_dir c ---------------------------------------------------- c Apply virtual temperature adjustment for upper level c ---------------------------------------------------- v_upper = t(l) if( p(l) .ge. 300.0 ) + v_upper = v_upper * ( 1.0 + ( 0.00061 * w(l) ) ) c ----------------------------------------------------- c Calculate exponential pressure height for upper layer c ----------------------------------------------------- algp_upper = alog( p(l) ) c ---------------- c Calculate height c ---------------- hgt = hgt + ( fac*(v_upper+v_lower)*(algp_lower-algp_upper) ) c ------------------------------- c Overwrite values for next layer c ------------------------------- v_lower = v_upper algp_lower = algp_upper c --------------------------------------------- c Store heights in same direction as other data c --------------------------------------------- z(l) = hgt end do return end subroutine clmozo(rlat,imon, omix) c * Obtain climatological ozone mixing-ratio profile by c interpolating in latitude and month amongst the c FASCODE model atmospheres. c .... version of 17.02.98 c * pressure coordinate is satellite-standard c input: c rlat = latitude (0-90,+N,-S) c imon = month (1,...,12) c output: c omix = ozone mixing ratio profile (ppmv) parameter (nl=40,ns=2,nz=3) dimension tlat(nz),ozmr(nl,nz,ns),omr(nl,ns),omix(nl) data tlat/15.,40.,65./ c * tropical data (ozmr(i,1,1),i=1,nl)/ + 0.55927,0.98223,1.94681,3.13498,4.30596,5.48908,7.41904,8.55498, + 9.22089,9.76594,9.60463,8.45820,6.99682,5.57585,4.30000,1.69985, + 1.23555,0.81780,0.39171,0.20944,0.14056,0.12283,0.10984,0.08405, + 0.06529,0.05392,0.04764,0.04366,0.04232,0.04052,0.03961,0.03735, + 0.03595,0.03534,0.03514,0.03385,0.03252,0.03107,0.03027,0.02901/ c * midlatitude winter data (ozmr(i,2,1),i=1,nl)/ + 0.58382,1.06610,2.23416,3.87594,5.18854,6.20949,7.04493,7.17104, + 7.10972,6.86107,6.29020,5.71788,5.35257,5.03882,4.57679,3.16918, + 2.47482,1.95801,1.43279,1.11336,0.93068,0.81309,0.74765,0.46038, + 0.25869,0.15587,0.10257,0.07960,0.06900,0.05625,0.05211,0.04120, + 0.03513,0.03297,0.03176,0.02883,0.02821,0.02796,0.02790,0.02781/ c * subarctic winter data (ozmr(i,3,1),i=1,nl)/ + 0.75492,1.20217,2.39283,3.75644,4.96742,5.64285,6.17910,6.22151, + 6.15197,5.88338,5.42542,4.91094,4.76030,4.63605,4.52229,3.70528, + 3.01350,2.39073,1.76424,1.38778,1.12046,0.82870,0.66060,0.36047, + 0.28088,0.17023,0.09235,0.06541,0.05169,0.04171,0.03941,0.03409, + 0.03095,0.02819,0.02673,0.02330,0.02159,0.01999,0.01933,0.01828/ c * tropical data (ozmr(i,1,2),i=1,nl)/ + 0.55927,0.98223,1.94681,3.13498,4.30596,5.48908,7.41904,8.55498, + 9.22089,9.76594,9.60463,8.45820,6.99682,5.57585,4.30000,1.69985, + 1.23555,0.81780,0.39171,0.20944,0.14056,0.12283,0.10984,0.08405, + 0.06529,0.05392,0.04764,0.04366,0.04232,0.04052,0.03961,0.03735, + 0.03595,0.03534,0.03514,0.03385,0.03252,0.03107,0.03027,0.02901/ c * midlatitude summer data (ozmr(i,2,2),i=1,nl)/ + 0.61951,1.07099,1.77685,2.91535,3.98547,5.06939,7.01746,8.18549, + 8.74393,8.73998,7.87205,6.65253,5.84694,5.12966,4.37609,2.46410, + 1.97307,1.47696,0.91259,0.66705,0.57759,0.48610,0.44729,0.24487, + 0.16974,0.12153,0.10001,0.08397,0.07669,0.06661,0.06225,0.05355, + 0.04892,0.04505,0.04291,0.03815,0.03517,0.03283,0.03194,0.03053/ c * subarctic summer data (ozmr(i,3,2),i=1,nl)/ + 0.58825,0.98312,1.64271,2.52996,3.56841,4.56575,6.36529,7.39635, + 7.78289,7.56976,6.69058,5.57509,5.21478,4.73043,4.34708,2.75529, + 2.05541,1.64657,1.17452,0.92611,0.78889,0.65317,0.58224,0.40198, + 0.24128,0.15314,0.10010,0.08052,0.07394,0.06544,0.06065,0.04986, + 0.04448,0.04090,0.03887,0.03446,0.03129,0.02823,0.02682,0.02456/ alat=abs(rlat) if(alat.gt.15.) go to 110 jl=1 kk=1 go to 120 110 if(alat.lt.65.) go to 130 jl=3 kk=2 120 do k=1,kk do i=1,nl omr(i,k)=ozmr(i,jl,k) enddo enddo if(kk.ne.1) go to 140 do i=1,nl omix(i)=omr(i,1) enddo go to 150 130 jl1=1 if(alat.gt.40.) jl1=2 jl2=jl1+1 wt1=(tlat(jl2)-alat)/25. wt2=1.-wt1 do k=1,ns do i=1,nl or1=ozmr(i,jl1,k) or2=ozmr(i,jl2,k) omr(i,k)=wt1*or1+wt2*or2 enddo enddo 140 nmon=imon if(rlat.lt.0.) nmon=nmon+6 if(nmon.gt.12) nmon=nmon-12 kmon=iabs(nmon-7) wt1=float(kmon)/6. wt2=1.-wt1 do i=1,nl omix(i)=wt1*omr(i,1)+wt2*omr(i,2) enddo 150 return end function satmix(p,t) c * Saturation mixing ratio c .... version of 18.02.98 if(t.gt.253.) then es=svpwat(t) else es=svpice(t) endif satmix=622.*es/p return end function svpice(temp) c * Saturation vapor pressure over ice c .... version of 18.02.98 c **** temp must be in degrees Kelvin real*8 t,a0,a1,a2,a3,a4,a5 data a0/.7859063157d0/,a1/.3579242320d-1/,a2/-.1292820828d-3/, * a3/.5937519208d-6/,a4/.4482949133d-9/,a5/.2176664827d-10/ t=temp-273.16 e=a0+t*(a1+t*(a2+t*(a3+t*(a4+t*a5)))) svpice=10.**e return end function svpwat(temp) c * Saturation vapor pressure over water c .... version of 18.02.98 c **** temp must be in degrees Kelvin real*8 b,s,t,a0,a1,a2,a3,a4,a5,a6,a7,a8,a9 data a0/.999996876d0/,a1/-.9082695004d-2/,a2/.7873616869d-4/, * a3/-.6111795727d-6/,a4/.4388418740d-8/,a5/-.2988388486d-10/, * a6/.2187442495d-12/,a7/-.1789232111d-14/,a8/.1111201803d-16/, * a9/-.3099457145d-19/,b/.61078d+1/ t=temp-273.16 s=a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+t*(a6+t*(a7+t*(a8+t*a9)))))))) s=b/s**8 svpwat=s return end subroutine taugen(nc,nx,ny,cc,xx,tau) c ... general transmittance calculation ... Strow-Woolf model c ... version of 19.09.96 dimension cc(nc,ny),xx(nx,ny),tau(*) data trap/-999.99/ taul=1. tau(1)=taul do 100 j=1,ny if(taul.eq.0.) go to 100 yy=cc(nc,j) if(yy.eq.trap) then taul=0. go to 100 endif do i=1,nx yy=yy+cc(i,j)*xx(i,j) enddo tauy=taul*exp(-yy) taul=amin1(tauy,taul) 100 tau(j+1)=taul return end subroutine tauwet(ncs,ncl,nxs,nxl,nxw,ny,ccs,ccl,xx,tau) c ... version of 17.02.98 dimension ccs(ncs,ny),ccl(ncl,ny),xx(nxw,ny),tau(*) data trap/-999.99/ odsum=0. taul=1. tau(1)=taul do 100 j=1,ny if(taul.eq.0.) go to 100 if(odsum.lt.5.) then yy=ccs(ncs,j) if(yy.eq.trap) then taul=0. go to 100 endif do i=1,nxs yy=yy+ccs(i,j)*xx(i,j) enddo odsum=odsum+abs(yy) else yy=ccl(ncl,j) if(yy.eq.trap) then taul=0. go to 100 endif do i=1,nxl yy=yy+ccl(i,j)*xx(i+11,j) enddo odsum=odsum+abs(yy) endif tauy=taul*exp(-yy) taul=amin1(tauy,taul) 100 tau(j+1)=taul return end subroutine vasrtw(tbbv,zena,kmon,ksat, t,w,*) c * Regression temperature & water-vapor mixing ratio c from VAS (GOES-4, -5, -6, -7) c .... version of 18.02.98 c * Inputs: c tbbv = VAS brightness temperatures (degK), channels 1-12 c zena = local (satellite) zenith angle (deg) c kmon = numerical month (1...12) c ksat = GOES satellite number (4...7) c * Outputs: c t = temperature profile (degK) at 40 standard levels c w = water-vapor mixing ratio profile (g/kg) at same levels c * = error return for I/O problems parameter (nk=12,nm=12,nt=40,nv=20) parameter (nx=nk+1,nc=nx+1,ny=nt+nv) parameter (isat=4,nsat=7,iuc=77,lenc=nc*ny*4) dimension tbbv(*),t(*),w(*) dimension coef(nc,ny,nm,isat:nsat),p(nt),x(nx),y(ny) character*12 cfile/'vasregtw.dat'/ data init/1/,wmin/0.003/ data p/0.1,0.2,0.5,1.0,1.5,2.,3.,4.,5.,7.,10.,15.,20.,25.,30., + 50.,60.,70.,85.,100.,115.,135.,150.,200.,250.,300.,350.,400., + 430.,475.,500.,570.,620.,670.,700.,780.,850.,920.,950.,1000./ secant(z) = 1./cos(0.01745329*z) if(init.ne.0) then c * store coefficients for all months and satellites open(iuc,file=cfile,recl=lenc,access='direct', + status='old',err=100) krec=0 do l=isat,nsat do k=1,nm krec=krec+1 read(iuc,rec=krec) ((coef(i,j,k,l), * i=1,nc),j=1,ny) enddo enddo close(iuc) init=0 endif c * define predictors do i=1,nk x(i)=tbbv(i) enddo x(nx)=10.*secant(zena)+200. c * obtain predictands k=kmon l=ksat do j=1,ny sum=coef(nc,j,k,l) do i=1,nx sum=sum+x(i)*coef(i,j,k,l) enddo y(j)=sum enddo c * extract temperature profile do j=1,nt t(j)=y(j) enddo c * extract mixing ratio profile do j=1,nv w(j)=wmin w(j+nv)=y(j+nt) enddo c * guard against supersaturation do j=nv+1,nt wmax=satmix(p(j),t(j)) w(j)=amin1(w(j),wmax) enddo return 100 return1 end er = t(n_levels) * ( 1.0 + ( 0.00061 * w(n_levels) ) ) algp_lower = alog( p(n_levels) ) i_start = n_levels-1 i_end = 1 else c --------------------- c Data stored bottom up c --------------------- v_lower = t(1) * ( 1.0 + ( 0.00061 * w(1) ) ) algp_lower = alog( p(1) ) i_start = 2 i_end = n_levels end if c----------------------------------------------------------------------- c vasregtw.dat 100644 000371 000325 00000473000 06533035310 014275 0 ustar 00russd stennis 000000 000000 U`>' ?614<)t=ug=O=̛e>l>#~%S>n?2N= 2==M>>p2> ?/Xj?Pվ>,U>!=]2=K:N u>^!:ۨ?zȩ?+K??[~I>g>1N>Hb<Ǽj>g#y?3b줾͉:?%8?~脴>!>{>=ѽ37="i#u>b?={-?#[?rTEg> >e=ɟ =c[üψ>ls1O?CLÅa