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