program checksum *_______________________________________________________________________ * created: 2/17/93 * author : Bill Ridgway (ridgway@climate.gsfc.nasa.gov, 301-286-9138) * purpose: check that spectral fluxes sum to quoted integrated values * routines : * rfm93 : read data in fm93 formatted data from unit 5 * chk93 : write brief (or long) summary of comparison to unit 6 *_______________________________________________________________________ * parameters: * mlayer: maximum number of model layers * mlevel: maximum number of model levels * mbands: maximum number of spectral bands *_______________________________________________________________________ * internal variables containing all model data * header(10): 10-line file header containing 80 characters per line * atmlbl : 3-character atmospheric descriptor * 'trp' tropical * 'mls' mid-latitude summer * 'mlw' mid-latitude winter * 'sas' sub-arctic summer * 'saw' sub-arctic winter * 'iso' isothermal * iphase : icrccm phase (see report) * icase : icrccm case number (see report) * nlayer : actual number of model layers * nlevel : actual number of model levels = nlayer+1 * nbands : actual number of spectral bands * psurf : surface pressure (mb) * tsurf : surface temperature (K) * ptrop : tropopause pressure (mb) * flgctm : flag to indicate h2o continuum is/is not (1/0) included * ppmco2 : co2 mixing ratio (ppm) * ppmch4 : ch4 mixing ratio (ppm) * ppmn2o : n2o mixing ratio (ppm) * pthair : total gas column (#air molecules/cm**2) * pthh2o : total h2o column (#h2o molecules/cm**2) * pthco2 : co2 path (#co2 molecules/cm**2) * ptho3 : ozone path (#ozone molecules/cm**2) * pthch4 : ch4 path (#ch4 molecules/cm**2) * pthn2o : n2o path (#n2o molecules/cm**2) * bandv1 : wavenumber lower limit of each band * bandv2 : wavenumber upper limit of each band * plevel : pressure levels starting at top (mb) * tlevel : temperature at pressure levels (K) * player : pressure at mid-point of each layer (mb) * tlayer : temperature characteristic of each layer (K) * alayer : air mass of layer (#air molecules/cm**2) * wlayer : water vapor volume mixing ratio of each layer (ppm) * olayer : ozone mixing ratio of each layer (ppm) * fluxup : radiative flux upward by level and band (w/m**2) * fluxdn : radiative flux downward by level and band (w/m**2) * fluxnt : net radiative flux upward by level and band (w/m**2) * heatrt : layer heating rates indexed by layer and band (K/day) * tropup : tropopause radiative flux upward by band (w/m**2) * tropdn : tropopause radiative flux downward by band (w/m**2) * tropnt : tropopause net radiative flux upward by band (w/m**2) *_______________________________________________________________________ * parameters defined in driver and passed to all subroutines: integer iounit,mlayer,mlevel,mbands parameter (mlayer=150,mlevel=151,mbands=600) * variables passed to or created by subroutines: character*80 header(10) character*3 atmlbl integer iphase,icase,nlayer,nlevel,nbands real psurf,tsurf,ptrop,flgctm real ppmco2,ppmch4,ppmn2o real pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o real bandv1(0:mbands),bandv2(0:mbands) real plevel(mlevel),tlevel(mlevel) real player(mlayer),tlayer(mlayer) real alayer(mlayer),wlayer(mlayer),olayer(mlayer) real fluxup(mlevel,0:mbands) real fluxdn(mlevel,0:mbands) real fluxnt(mlevel,0:mbands) real heatrt(mlayer,0:mbands) real tropup( 0:mbands) real tropdn( 0:mbands) real tropnt( 0:mbands) *_______________________________________________________________________ iounit=5 call rfm93 + (iounit,mlayer,mlevel,mbands + ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands + ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o + ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o + ,bandv1,bandv2,plevel,tlevel + ,player,tlayer,alayer,wlayer,olayer + ,fluxup,fluxdn,fluxnt,heatrt + ,tropup,tropdn,tropnt) iounit=6 call chk93 + (iounit,mlayer,mlevel,mbands + ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands + ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o + ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o + ,bandv1,bandv2,plevel,tlevel + ,player,tlayer,alayer,wlayer,olayer + ,fluxup,fluxdn,fluxnt,heatrt + ,tropup,tropdn,tropnt) stop end subroutine rfm93 + (iounit,mlayer,mlevel,mbands + ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands + ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o + ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o + ,bandv1,bandv2,plevel,tlevel + ,player,tlayer,alayer,wlayer,olayer + ,fluxup,fluxdn,fluxnt,heatrt + ,tropup,tropdn,tropnt) *_______________________________________________________________________ * purpose: rfm93 - read 1993-format icrccm data file, return data arrays *_______________________________________________________________________ * parameters defined in driver and passed to all subroutines: integer iounit,mlayer,mlevel,mbands * variables passed to or created by subroutines: character*80 header(10) character*3 atmlbl integer iphase,icase,nlayer,nlevel,nbands real psurf,tsurf,ptrop,flgctm real ppmco2,ppmch4,ppmn2o real pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o real bandv1(0:mbands),bandv2(0:mbands) real plevel(mlevel),tlevel(mlevel) real player(mlayer),tlayer(mlayer) real alayer(mlayer),wlayer(mlayer),olayer(mlayer) real fluxup(mlevel,0:mbands) real fluxdn(mlevel,0:mbands) real fluxnt(mlevel,0:mbands) real heatrt(mlayer,0:mbands) real tropup( 0:mbands) real tropdn( 0:mbands) real tropnt( 0:mbands) *_______________________________________________________________________ * blank out header initially do k=1,10 write(header(k),'(80a1)') (' ',j=1,80) end do * read 10 lines of header info from fm93 file read (iounit,'(a80)') header * read key parameters which identify this case read (iounit,910) iphase,icase,atmlbl,nlayer,nlevel,nbands, + psurf,tsurf,ptrop,flgctm, + ppmco2,ppmch4,ppmn2o, + pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o 910 format(24x,i16, + /,24x,i16, + /,24x,13x,a3, + 3(/,24x,i16), + 3(/,24x,f16.2), + 1(/,24x,f16.0), + 3(/,24x,f16.2), + 6(/,24x,f16.4)) * ingest spectrally summed results (k=0 only) k=0 read (iounit,'(24x,f16.2)') bandv1(k),bandv2(k) read (iounit,'(24x,e16.4)') fluxup(nlevel,k),fluxdn(nlevel,k), + tropup( k),tropdn( k), + fluxup( 1,k) tropnt(k)=tropup(k)-tropdn(k) read (iounit,'(1x)') read (iounit,'(1x)') do j=1,nlayer read (iounit,'(6e12.4)') + player(j),tlayer(j),alayer(j),wlayer(j),olayer(j),heatrt(j,k) end do read (iounit,'(1x)') read (iounit,'(1x)') do j=1,nlevel read (iounit,'(5e12.4)') + plevel(j),tlevel(j),fluxup(j,k),fluxdn(j,k),fluxnt(j,k) end do * read band-by-band flux summaries, fluxes, heating rates do k=1,nbands read(iounit,'(2f15.2)') bandv1(k),bandv2(k) read(iounit,'(5e15.6)') fluxup(nlevel,k),fluxdn(nlevel,k), + tropup( k),tropdn( k), + fluxup( 1,k) read(iounit,'(5e15.6)') (fluxup(j,k),j=1,nlevel) read(iounit,'(5e15.6)') (fluxdn(j,k),j=1,nlevel) read(iounit,'(5e15.6)') (fluxnt(j,k),j=1,nlevel) read(iounit,'(5e15.6)') (heatrt(j,k),j=1,nlayer) end do return end subroutine chk93 + (iounit,mlayer,mlevel,mbands + ,header,atmlbl,iphase,icase,nlayer,nlevel,nbands + ,psurf,tsurf,ptrop,flgctm,ppmco2,ppmch4,ppmn2o + ,pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o + ,bandv1,bandv2,plevel,tlevel + ,player,tlayer,alayer,wlayer,olayer + ,fluxup,fluxdn,fluxnt,heatrt + ,tropup,tropdn,tropnt) *_______________________________________________________________________ * purpose: chk93 - check that spectral fluxes sum to correct sums *_______________________________________________________________________ * parameters defined in driver and passed to all subroutines: integer iounit,mlayer,mlevel,mbands * variables passed to or created by subroutines: character*80 header(10) character*3 atmlbl integer iphase,icase,nlayer,nlevel,nbands real psurf,tsurf,ptrop,flgctm real ppmco2,ppmch4,ppmn2o real pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o real bandv1(0:mbands),bandv2(0:mbands) real plevel(mlevel),tlevel(mlevel) real player(mlayer),tlayer(mlayer) real alayer(mlayer),wlayer(mlayer),olayer(mlayer) real fluxup(mlevel,0:mbands) real fluxdn(mlevel,0:mbands) real fluxnt(mlevel,0:mbands) real heatrt(mlayer,0:mbands) real tropup( 0:mbands) real tropdn( 0:mbands) real tropnt( 0:mbands) parameter (maxx=300) real sumfup(maxx) real sumfdn(maxx) real sumfnt(maxx) real sumhrt(maxx) real sumtup real sumtdn real sumtnt *_______________________________________________________________________ * initialize sums to zero sumtup=0. sumtdn=0. sumtnt=0. do j=1,nlayer sumhrt(j)=0. end do do j=1,nlevel sumfup(j)=0. sumfdn(j)=0. sumfnt(j)=0. end do * sum key fluxes and heating rates by band do k=1,nbands sumtup=sumtup+tropup(k) sumtdn=sumtdn+tropdn(k) sumtnt=sumtnt+tropnt(k) do j=1,nlayer sumhrt(j)=sumhrt(j)+heatrt(j,k) end do do j=1,nlevel sumfup(j)=sumfup(j)+fluxup(j,k) sumfdn(j)=sumfdn(j)+fluxdn(j,k) sumfnt(j)=sumfnt(j)+fluxnt(j,k) end do end do * compute ratios with stored sums, min and max ratio rmin=1. rmax=1. if(sumtup.ne.0.) sumtup=sumtup/tropup(0) if(sumtdn.ne.0.) sumtdn=sumtdn/tropdn(0) if(sumtnt.ne.0.) sumtnt=sumtnt/tropnt(0) if(sumtup.ne.0.) rmin=min(rmin,sumtup) if(sumtdn.ne.0.) rmin=min(rmin,sumtdn) if(sumtnt.ne.0.) rmin=min(rmin,sumtnt) rmax=max(rmax,sumtup) rmax=max(rmax,sumtdn) rmax=max(rmax,sumtnt) do j=1,nlayer if(sumhrt(j).ne.0.) sumhrt(j)=sumhrt(j)/heatrt(j,0) if(sumhrt(j).ne.0.) rmin=min(rmin,sumhrt(j)) rmax=max(rmax,sumhrt(j)) end do do j=1,nlevel if(sumfup(j).ne.0.) sumfup(j)=sumfup(j)/fluxup(j,0) if(sumfdn(j).ne.0.) sumfdn(j)=sumfdn(j)/fluxdn(j,0) if(sumfnt(j).ne.0.) sumfnt(j)=sumfnt(j)/fluxnt(j,0) if(sumfup(j).ne.0.) rmin=min(rmin,sumfup(j)) if(sumfdn(j).ne.0.) rmin=min(rmin,sumfdn(j)) if(sumfnt(j).ne.0.) rmin=min(rmin,sumfnt(j)) rmax=max(rmax,sumfup(j)) rmax=max(rmax,sumfdn(j)) rmax=max(rmax,sumfnt(j)) end do write(*,'(15x,"phase=",i2, & 5x,"case=",i2, & 5x,"ratio in range [",2f9.6," ]")') & iphase,icase,rmin,rmax if(0.999.lt.rmin.and.rmax.lt.1.001) stop * print file header with 10 lines of info in any format write(iounit,'(a80)') header * print descriptive parameters write(iounit,910) iphase,icase,atmlbl,nlayer,nlevel,nbands, + psurf,tsurf,ptrop,flgctm, + ppmco2,ppmch4,ppmn2o, + pthair,pthh2o,pthco2,ptho3,pthch4,pthn2o 910 format('.....icrccm.phase.number',i16 + ,/,'......icrccm.case.number',i16 + ,/,'.......atmospheric.label',a16 + ,/,'.....number.model.layers',i16 + ,/,'.....number.model.levels',i16 + ,/,'...number.spectral.bands',i16 + ,/,'........surface.pressure',0p,f16.2,' mb' + ,/,'.....surface.temperature',0p,f16.2,' k' + ,/,'.....tropopause.pressure',0p,f16.2,' mb' + ,/,'......h2o.continuum.flag',0p,f16.0,' (0=off,1=included)' + ,/,'........co2.mixing.ratio',0p,f16.2,' ppm' + ,/,'........ch4.mixing.ratio',0p,f16.2,' ppm' + ,/,'........n2o.mixing.ratio',0p,f16.2,' ppm' + ,/,'....total.dry.air.column',1p,e16.4,' molecules/cm**2' + ,/,'........total.h2o.column',1p,e16.4,' molecules/cm**2' + ,/,'........total.co2.column',1p,e16.4,' molecules/cm**2' + ,/,'.........total.o3.column',1p,e16.4,' molecules/cm**2' + ,/,'........total.ch4.column',1p,e16.4,' molecules/cm**2' + ,/,'........total.n2o.column',1p,e16.4,' molecules/cm**2') * print spectrally itegrated summary of icrccm fluxes k=0 write(iounit,920) bandv1(k),bandv2(k) 920 format('....spectral.region.from',0p,f16.2,' cm-1 ' + ,/,'......................to',0p,f16.2,' cm-1 ') write(iounit,930) sumfup(nlevel ),sumfdn(nlevel ), + sumtup ,sumtdn , + sumfup( 1 ) 930 format('.......flux.up.@.surface',1p,e16.4,' watts/m**2' + ,/,'.....flux.down.@.surface',1p,e16.4,' watts/m**2' + ,/,'....flux.up.@.tropopause',1p,e16.4,' watts/m**2' + ,/,'..flux.down.@.tropopause',1p,e16.4,' watts/m**2' + ,/,'..flux.up.@.top.of.atmos',1p,e16.4,' watts/m**2') * print table of atmospheric profile data and layer heating rates write(iounit,'(2a)') + ' layer >> p t air mass' +,' h2o o3 heating' +,' (mb) (K) (#/cm**2)' +,' (ppm) (ppm) (deg/day)' do j=1,nlayer write(iounit,'(1p,e12.4,3p,e12.4,1p,4e12.4)') + player(j),tlayer(j),alayer(j),wlayer(j),olayer(j),sumhrt(j) end do * print table of computed fluxes write(iounit,'(a)') + ' level >> p t flux up flux down net flux' +,' (mb) (k) (w/m**2) (w/m**2) (w/m**2)' do j=1,nlevel write(iounit,'(1p,e12.4,3p,e12.4,1p,3e12.4)') + plevel(j),tlevel(j),sumfup(j),sumfdn(j),sumfnt(j) end do return end