C +-======-+ C Copyright (c) 2003-2007 United States Government as represented by C the Admistrator of the National Aeronautics and Space Administration. C All Rights Reserved. C C THIS OPEN SOURCE AGREEMENT ("AGREEMENT") DEFINES THE RIGHTS OF USE, C REPRODUCTION, DISTRIBUTION, MODIFICATION AND REDISTRIBUTION OF CERTAIN C COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS C REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY"). C THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN C INTENDED THIRD-PARTY BENEFICIARY OF ALL SUBSEQUENT DISTRIBUTIONS OR C REDISTRIBUTIONS OF THE SUBJECT SOFTWARE. ANYONE WHO USES, REPRODUCES, C DISTRIBUTES, MODIFIES OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED C HEREIN, OR ANY PART THEREOF, IS, BY THAT ACTION, ACCEPTING IN FULL THE C RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT. C C Government Agency: National Aeronautics and Space Administration C Government Agency Original Software Designation: GSC-15354-1 C Government Agency Original Software Title: GEOS-5 GCM Modeling Software C User Registration Requested. Please Visit http://opensource.gsfc.nasa.gov C Government Agency Point of Contact for Original Software: C Dale Hithon, SRA Assistant, (301) 286-2691 C C +-======-+ parameter ( im = 144 ) ! x-dimension for unified comparison parameter ( jm = 91 ) ! y-dimension for unified comparison parameter ( nr = 10 ) ! r-dimension (number of geographic regions ) parameter ( nl = 10 ) ! l-dimension (number of upper-air levels ) parameter ( nf = 6 ) ! number of different fields to process integer rc, fcst_file real f(im,jm), fvar(im,jm), varf(im,jm), fstar(im,jm), fsprime(im,jm) ! forecast variables real a(im,jm), avar(im,jm), vara(im,jm), astar(im,jm), asprime(im,jm) ! analysis variables real c(im,jm), cvar(im,jm), varc(im,jm), cstar(im,jm) ! climatology variables real corr(nr,nl,nf,100) ! Note: Hardwired for 100 time periods (Max) real rms(nr,nl,nf,100,5) ! Note: Hardwired for 100 time periods (Max) real*4 dum(nr) real pc(im,jm) real uc(im,jm,nl) real vc(im,jm,nl) real tc(im,jm,nl) real qc(im,jm,nl) real hc(im,jm,nl) real pf(im,jm) real uf(im,jm,nl) real vf(im,jm,nl) real tf(im,jm,nl) real qf(im,jm,nl) real hf(im,jm,nl) real pe(im,jm) real ue(im,jm,nl) real ve(im,jm,nl) real te(im,jm,nl) real qe(im,jm,nl) real he(im,jm,nl) real pa(im,jm) real ua(im,jm,nl) real va(im,jm,nl) real ta(im,jm,nl) real qa(im,jm,nl) real ha(im,jm,nl) real zlev(nl) data zlev / 1000.0, 850.0, 700.0, 500.0, 400.0, . 300.0, 250.0, 200.0, 150.0, 100.0 / real pref real zlat1(nr), lat1 real zlat2(nr), lat2 real zlon1(nr), lon1 real zlon2(nr), lon2 data zlat1 / -90., 20., -20., -80., 0., 0., -90.,-90., 20., 30. / data zlat2 / 90., 80., 20., -20., 90., 90., 0., 0., 60., 60. / data zlon1 /-180.,-180.,-180.,-180.,-180., 0.,-180., 0.,-140., -10. / data zlon2 / 180., 180., 180., 180., 0.,180., 0.,180.,- 60., 30. / character*8 qname(6*nf), name, suffix character*80 qdesc(6*nf) character*1 char data qname(01) /'pcor '/ data qname(02) /'ucor '/ data qname(03) /'vcor '/ data qname(04) /'tcor '/ data qname(05) /'qcor '/ data qname(06) /'hcor '/ data qname(07) /'prms '/ data qname(08) /'urms '/ data qname(09) /'vrms '/ data qname(10) /'trms '/ data qname(11) /'qrms '/ data qname(12) /'hrms '/ data qname(13) /'prms_ran'/ data qname(14) /'urms_ran'/ data qname(15) /'vrms_ran'/ data qname(16) /'trms_ran'/ data qname(17) /'qrms_ran'/ data qname(18) /'hrms_ran'/ data qname(19) /'prms_bar'/ data qname(20) /'urms_bar'/ data qname(21) /'vrms_bar'/ data qname(22) /'trms_bar'/ data qname(23) /'qrms_bar'/ data qname(24) /'hrms_bar'/ data qname(25) /'prms_dis'/ data qname(26) /'urms_dis'/ data qname(27) /'vrms_dis'/ data qname(28) /'trms_dis'/ data qname(29) /'qrms_dis'/ data qname(30) /'hrms_dis'/ data qname(31) /'prms_dsp'/ data qname(32) /'urms_dsp'/ data qname(33) /'vrms_dsp'/ data qname(34) /'trms_dsp'/ data qname(35) /'qrms_dsp'/ data qname(36) /'hrms_dsp'/ data qdesc(01) /'Sea Level Pressure Anomaly Correlation'/ data qdesc(02) /'Upper Air U-Wind Anomaly Correlation'/ data qdesc(03) /'Upper Air V-Wind Anomaly Correlation'/ data qdesc(04) /'Temperature Anomaly Correlation'/ data qdesc(05) /'Specific Humidity Anomaly Correlation'/ data qdesc(06) /'Virtual Heights Anomaly Correlation'/ data qdesc(07) /'Sea Level Pressure Root Mean Square Error'/ data qdesc(08) /'Upper Air U-Wind Root Mean Square Error'/ data qdesc(09) /'Upper Air V-Wind Root Mean Square Error'/ data qdesc(10) /'Temperature Root Mean Square Error'/ data qdesc(11) /'Specific Humidity Root Mean Square Error'/ data qdesc(12) /'Virtual Heights Root Mean Square Error'/ data qdesc(13) /'Sea Level Pressure Root Mean Square Error Random'/ data qdesc(14) /'Upper Air U-Wind Root Mean Square Error Random'/ data qdesc(15) /'Upper Air V-Wind Root Mean Square Error Random'/ data qdesc(16) /'Temperature Root Mean Square Error Random'/ data qdesc(17) /'Specific Humidity Root Mean Square Error Random'/ data qdesc(18) /'Virtual Heights Root Mean Square Error Random'/ data qdesc(19) /'Sea Level Pressure Root Mean Square Error Bias'/ data qdesc(20) /'Upper Air U-Wind Root Mean Square Error Bias'/ data qdesc(21) /'Upper Air V-Wind Root Mean Square Error Bias'/ data qdesc(22) /'Temperature Root Mean Square Error Bias'/ data qdesc(23) /'Specific Humidity Root Mean Square Error Bias'/ data qdesc(24) /'Virtual Heights Root Mean Square Error Bias'/ data qdesc(25) /'Sea Level Pressure Root Mean Square Error Dissipation'/ data qdesc(26) /'Upper Air U-Wind Root Mean Square Error Dissipation'/ data qdesc(27) /'Upper Air V-Wind Root Mean Square Error Dissipation'/ data qdesc(28) /'Temperature Root Mean Square Error Dissipation'/ data qdesc(29) /'Specific Humidity Root Mean Square Error Dissipation'/ data qdesc(30) /'Virtual Heights Root Mean Square Error Dissipation'/ data qdesc(31) /'Sea Level Pressure Root Mean Square Error Dispersion'/ data qdesc(32) /'Upper Air U-Wind Root Mean Square Error Dispersion'/ data qdesc(33) /'Upper Air V-Wind Root Mean Square Error Dispersion'/ data qdesc(34) /'Temperature Root Mean Square Error Dispersion'/ data qdesc(35) /'Specific Humidity Root Mean Square Error Dispersion'/ data qdesc(36) /'Virtual Heights Root Mean Square Error Dispersion'/ character*256 tag character*256, allocatable :: arg(:) character*256, allocatable :: fname(:) character*256, allocatable :: aname(:) character*256, allocatable :: cname(:) character*256 ename ! Systematic Error File integer num_fcst_files integer num_clm_files integer num_ana_files real undef real fundef real eundef real aundef real cundef integer fcst_im,fcst_jm,nargs integer n,nt,time,ntimes,nvars integer yymmdd(100) integer hhmmss(100) character*256 vname(500) logical defined logical syserr data syserr /.false./ real rho,d1,d2,d3 integer nymdb, nhmsb, nymdc integer nymde, nhmse, nhmsc integer month, day, year, hour, minute, hmax, timinc, nfreq, ndt0, ndt, num character*7 xdate character*9 bdate, edate character*4 bhour, ehour character*256 statfile character*256 datafile character*256 ctlfile1 character*256 ctlfile2 character*256 ctlfile3 character*3 months(12) data months /'JAN','FEB','MAR','APR','MAY','JUN', . 'JUL','AUG','SEP','OCT','NOV','DEC'/ C ********************************************************************** C **** Initialization **** C ********************************************************************** call timebeg ('main') call timebeg ('init') pref = 500 nfreq = -999 hmax = -999 tag = "" nargs = iargc() if( nargs.eq.0 ) then call usage() else allocate ( arg(nargs) ) do n=1,nargs call getarg(n,arg(n)) enddo do n=1,nargs if( trim(arg(n)).eq.'-fcst' ) then nfiles = 1 read(arg(n+nfiles),fmt='(a1)') char do while (char.ne.'-' .and. n+nfiles.ne.nargs ) nfiles = nfiles+1 read(arg(n+nfiles),fmt='(a1)') char enddo if( char.eq.'-' ) nfiles = nfiles-1 allocate ( fname(nfiles) ) do m=1,nfiles fname(m) = arg(n+m) enddo num_fcst_files = nfiles endif if( trim(arg(n)).eq.'-ana' ) then nfiles = 1 read(arg(n+nfiles),fmt='(a1)') char do while (char.ne.'-' .and. n+nfiles.ne.nargs ) nfiles = nfiles+1 read(arg(n+nfiles),fmt='(a1)') char enddo if( char.eq.'-' ) nfiles = nfiles-1 allocate ( aname(nfiles) ) do m=1,nfiles aname(m) = arg(n+m) enddo num_ana_files = nfiles endif if( trim(arg(n)).eq.'-cli' ) then nfiles = 1 read(arg(n+nfiles),fmt='(a1)') char do while (char.ne.'-' .and. n+nfiles.ne.nargs ) nfiles = nfiles+1 read(arg(n+nfiles),fmt='(a1)') char enddo if( char.eq.'-' ) nfiles = nfiles-1 allocate ( cname(nfiles) ) do m=1,nfiles cname(m) = arg(n+m) enddo num_clm_files = nfiles endif if( trim(arg(n)).eq.'-syserr' ) then ename = trim(arg(n+1)) syserr = .true. endif if( trim(arg(n)).eq.'-tag' ) tag = trim(arg(n+1)) if( trim(arg(n)).eq.'-hmax' ) read(arg(n+1),*) hmax if( trim(arg(n)).eq.'-nfreq' ) read(arg(n+1),*) nfreq if( trim(arg(n)).eq.'-pref' ) read(arg(n+1),*) pref enddo endif if( trim(tag).ne."" ) tag = trim(tag) // "." print * print *, 'Computing Stats for Forecast File(s):' print *, '-------------------------------------' do n=1,num_fcst_files write(6,1001) trim(fname(n)) enddo 1001 format(1x,a) print * print *, 'Verifying Analysis File(s):' print *, '---------------------------' do n=1,num_ana_files write(6,1001) trim(aname(n)) enddo print * print *, 'Climatology File(s):' print *, '-------------------' do n=1,num_clm_files write(6,1001) trim(cname(n)) enddo print * if( syserr ) then print *, 'Systematic Error File:' print *, '----------------------' write(6,1001) trim(ename) print * endif call timeend ('init') C ********************************************************************** pi = 4.0*atan(1.0) dl = 2*pi/im dp = pi/(jm-1) c Loop over Forecast Files c ------------------------ nt = 0 do fcst_file = 1,num_fcst_files call init_fcst ( fname(fcst_file),fcst_id,fundef,ntimes,yymmdd,hhmmss,timinc ) if( nt.eq.0 ) undef = fundef ndt0 = nsecf(timinc) if( nfreq.eq.-999 ) then ndt = nsecf(timinc) else ndt = nsecf(nfreq) endif c Read Forecast c ------------- do time = 1, ntimes nymd = yymmdd(time) nhms = hhmmss(time) if( nt.eq.0 ) then call read_fcst( fcst_id,nymd,nhms,pf,uf,vf,tf,qf,hf,im,jm,nl,zlev ) ! Tick Clock for Next Desired Time Period ! --------------------------------------- nymdc = nymd nhmsc = nhms nt = nt + 1 call tick (nymdc,nhmsc,ndt) nhmsc = 10000*(nhmsc/10000) ! Kludge to strip off minutes/seconds from Model Time else if( (nymd.lt.nymdc).or. . ((nymd.eq.nymdc).and.(nhms.lt.nhmsc)) ) cycle dowhile( nymd.ne.nymdc .or. nhms.ne.nhmsc ) print * print *, 'Skipping Forecast Time: ',nymdc,nhmsc,' (not available)' c stop c Note: Uncomment following lines if missing Forecast Times should be included as undefs c -------------------------------------------------------------------------------------- c nt = nt + 1 c corr(:,:,:,nt) = fundef c rms(:,:,:,nt,1) = fundef c rms(:,:,:,nt,2) = fundef c rms(:,:,:,nt,3) = fundef c rms(:,:,:,nt,4) = fundef c rms(:,:,:,nt,5) = fundef call tick (nymdc,nhmsc,ndt0) enddo nt = nt + 1 call tick (nymdc,nhmsc,ndt) call read_fcst( fcst_id,nymd,nhms,pf,uf,vf,tf,qf,hf,im,jm,nl,zlev ) endif nhms = 10000*(nhms/10000) ! Kludge to strip off minutes/seconds from Model Time c Find Analysis and Climatology for Forecast Time: nymd,nhms c ---------------------------------------------------------- call read_anal ( nymd,nhms,pa,ua,va,ta,qa,ha,im,jm,nl,zlev,aname,num_ana_files,aundef ) call read_clim_hdf ( nymd,nhms,pc,uc,vc,tc,qc,hc,im,jm,nl,zlev,cname,num_clm_files,cundef ) c Read Systematic Error File and Modify Forecast Fields c ----------------------------------------------------- if( syserr ) then call read_syserr( pe,ue,ve,te,qe,he,im,jm,nl,zlev,ename,ndt,eundef ) do j=1,jm do i=1,im if( defined(pf(i,j),fundef) .and. defined(pe(i,j),eundef) ) pf(i,j) = pf(i,j) - pe(i,j) enddo enddo do L=1,nl do j=1,jm do i=1,im if( defined(uf(i,j,L),fundef) .and. defined(ue(i,j,L),eundef) ) uf(i,j,L) = uf(i,j,L) - ue(i,j,L) if( defined(vf(i,j,L),fundef) .and. defined(ve(i,j,L),eundef) ) vf(i,j,L) = vf(i,j,L) - ve(i,j,L) if( defined(tf(i,j,L),fundef) .and. defined(te(i,j,L),eundef) ) tf(i,j,L) = tf(i,j,L) - te(i,j,L) if( defined(qf(i,j,L),fundef) .and. defined(qe(i,j,L),eundef) ) qf(i,j,L) = qf(i,j,L) - qe(i,j,L) if( defined(hf(i,j,L),fundef) .and. defined(he(i,j,L),eundef) ) hf(i,j,L) = hf(i,j,L) - he(i,j,L) enddo enddo enddo endif if( nt.eq.1 ) then nymdb = nymd nhmsb = nhms write(bdate,3001) nymdb write(xdate,3003) nhmsb statfile = trim(tag) // 'globl.' // bdate // "." // xdate // ".data" open (51,file=trim(statfile),form='unformatted',access='sequential') if( nfreq.ne.-999 ) then print *, 'Forecast Frequency: ',ndt/3600,' (hrs)' print * endif endif call writit ( pf,im,jm,1 ,fundef,undef ) call writit ( uf,im,jm,nl,fundef,undef ) call writit ( vf,im,jm,nl,fundef,undef ) call writit ( tf,im,jm,nl,fundef,undef ) call writit ( qf,im,jm,nl,fundef,undef ) call writit ( hf,im,jm,nl,fundef,undef ) call writit ( pa,im,jm,1 ,aundef,undef ) call writit ( ua,im,jm,nl,aundef,undef ) call writit ( va,im,jm,nl,aundef,undef ) call writit ( ta,im,jm,nl,aundef,undef ) call writit ( qa,im,jm,nl,aundef,undef ) call writit ( ha,im,jm,nl,aundef,undef ) call writit ( pc,im,jm,1 ,cundef,undef ) call writit ( uc,im,jm,nl,cundef,undef ) call writit ( vc,im,jm,nl,cundef,undef ) call writit ( tc,im,jm,nl,cundef,undef ) call writit ( qc,im,jm,nl,cundef,undef ) call writit ( hc,im,jm,nl,cundef,undef ) c Loop over Geographical Regions c ------------------------------ do 1000 iregion = 1,nr lat1 = zlat1(iregion) lat2 = zlat2(iregion) lon1 = zlon1(iregion) lon2 = zlon2(iregion) c Determine beginning and ending i&j for Region c --------------------------------------------- call bounds (lat1,lat2,lon1,lon2, . jbeg,jend,ibeg,iend,im,jm) c Loop over Fields and Levels c --------------------------- do 2000 nfield = 1,nf nlev = nl if( nfield.eq.1 ) nlev = 1 do 3000 lev = 1,nlev c Sea-Level Pressure c ------------------ if( nfield.eq.1 ) then do j=1,jm do i=1,im c(i,j) = pc(i,j) f(i,j) = pf(i,j) a(i,j) = pa(i,j) enddo enddo endif c U-Wind c ------ if( nfield.eq.2 ) then do j=1,jm do i=1,im c(i,j) = uc(i,j,lev) f(i,j) = uf(i,j,lev) a(i,j) = ua(i,j,lev) enddo enddo endif c V-Wind c ------ if( nfield.eq.3 ) then do j=1,jm do i=1,im c(i,j) = vc(i,j,lev) f(i,j) = vf(i,j,lev) a(i,j) = va(i,j,lev) enddo enddo endif c Temperature c ----------- if( nfield.eq.4 ) then do j=1,jm do i=1,im c(i,j) = tc(i,j,lev) f(i,j) = tf(i,j,lev) a(i,j) = ta(i,j,lev) enddo enddo endif c Specific Humidity c ----------------- if( nfield.eq.5 ) then do j=1,jm do i=1,im c(i,j) = qc(i,j,lev) f(i,j) = qf(i,j,lev) a(i,j) = qa(i,j,lev) enddo enddo endif c Virtual Heights c --------------- if( nfield.eq.6 ) then do j=1,jm do i=1,im c(i,j) = hc(i,j,lev) f(i,j) = hf(i,j,lev) a(i,j) = ha(i,j,lev) enddo enddo endif c Compute Regional Area Means c --------------------------- fmean = 0 amean = 0 cmean = 0 area = 0 do j=jbeg,jend phi = -pi/2 + (j-1)*dp cosp = cos(phi) do i=ibeg,iend if( defined(f(i,j),fundef) .and. . defined(a(i,j),aundef) .and. . defined(c(i,j),cundef) ) then fmean = fmean + f(i,j)*cosp amean = amean + a(i,j)*cosp cmean = cmean + c(i,j)*cosp area = area + cosp endif enddo enddo if( area.ne.0.0 ) then fmean = fmean / area amean = amean / area cmean = cmean / area else fmean = fundef amean = aundef cmean = cundef endif c Define Deviations from Area Means c --------------------------------- do j=jbeg,jend do i=ibeg,iend if( defined(f(i,j),fundef) .and. . defined(a(i,j),aundef) .and. . defined(c(i,j),cundef) ) then fstar(i,j) = f(i,j)-fmean astar(i,j) = a(i,j)-amean cstar(i,j) = c(i,j)-cmean else fstar(i,j) = fundef astar(i,j) = aundef cstar(i,j) = cundef endif enddo enddo c Subtract Climatology Deviations and Compute Variances c ----------------------------------------------------- do j=jbeg,jend do i=ibeg,iend if( defined(fstar(i,j),fundef) .and. . defined(astar(i,j),aundef) .and. . defined(cstar(i,j),cundef) ) then fsprime(i,j) = fstar(i,j)- cstar(i,j) asprime(i,j) = astar(i,j)- cstar(i,j) fvar(i,j) = fsprime(i,j)*fsprime(i,j) avar(i,j) = asprime(i,j)*asprime(i,j) cvar(i,j) = fsprime(i,j)*asprime(i,j) varf(i,j) = fstar(i,j)* fstar(i,j) vara(i,j) = astar(i,j)* astar(i,j) varc(i,j) = fstar(i,j)* astar(i,j) else fsprime(i,j) = fundef asprime(i,j) = aundef fvar(i,j) = fundef avar(i,j) = aundef cvar(i,j) = cundef varf(i,j) = fundef vara(i,j) = aundef varc(i,j) = cundef endif enddo enddo c Compute Regional Mean Variances and Anomaly Correlation c ------------------------------------------------------- c1 = 0 c2 = 0 c3 = 0 d1 = 0 d2 = 0 d3 = 0 r1 = 0 r2 = 0 r3 = 0 area = 0 do j=jbeg,jend phi = -pi/2 + (j-1)*dp cosp = cos(phi) do i=ibeg,iend if( defined( f(i,j),fundef) .and. . defined( a(i,j),aundef) .and. . defined(fvar(i,j),fundef) .and. . defined(avar(i,j),aundef) .and. . defined(cvar(i,j),cundef) ) then c1 = c1 + cvar(i,j) * cosp c2 = c2 + fvar(i,j) * cosp c3 = c3 + avar(i,j) * cosp d1 = d1 + varc(i,j) * cosp d2 = d2 + varf(i,j) * cosp d3 = d3 + vara(i,j) * cosp r1 = r1 + ( f(i,j)- a(i,j))**2 * cosp r2 = r2 + (fstar(i,j)-astar(i,j))**2 * cosp r3 = r3 + (fmean -amean )**2 * cosp area = area + cosp endif enddo enddo if( area.ne.0 ) then c1 = c1 / area c2 = c2 / area c3 = c3 / area d1 = d1 / area d2 = d2 / area d3 = d3 / area r1 = r1 / area r2 = r2 / area r3 = r3 / area else c1 = undef c2 = undef c3 = undef d1 = undef d2 = undef d3 = undef r1 = undef r2 = undef r3 = undef endif if( c1.ne.undef ) then corr(iregion,lev,nfield,nt) = c1/(sqrt(c2)*sqrt(c3)) rms(iregion,lev,nfield,nt,1) = sqrt(r1) rms(iregion,lev,nfield,nt,2) = sqrt(r2) rms(iregion,lev,nfield,nt,3) = sqrt(r3) rho = d1/(sqrt(d2)*sqrt(d3)) rms(iregion,lev,nfield,nt,4) = (sqrt(d2)-sqrt(d3))**2 + (fmean-amean)**2 rms(iregion,lev,nfield,nt,5) = 2*(1-rho)*sqrt(d2)*sqrt(d3) else corr(iregion,lev,nfield,nt) = undef rms(iregion,lev,nfield,nt,1) = undef rms(iregion,lev,nfield,nt,2) = undef rms(iregion,lev,nfield,nt,3) = undef rms(iregion,lev,nfield,nt,4) = undef rms(iregion,lev,nfield,nt,5) = undef endif if( nfield.eq.6 .and. iregion.eq.2 .and. zlev(lev).eq.pref ) then if( nt.eq.1 ) then hour = 0 else call interp_time ( nymdb,nhmsb, nymd,nhms, nt, num ) hour = num*(nt-1) endif nymde = nymd nhmse = nhms write(6,1003) int(pref),nymd,nhms,hour,corr(iregion,lev,nfield,nt),fmean,amean,cmean 1003 format(1x,i4,'-mb NH Height at ',i8,2x,i6.6,' (',i3,' hrs) corr: ',f9.7, . 3x,'fcst: ',f7.2,2x,'ana: ',f7.2,2x,'cli: ',f7.2,' (m)') endif 3000 continue ! end level loop 2000 continue ! end field loop 1000 continue ! end region loop if( hour.eq.hmax ) exit enddo ! end time loop for current forecast file call gfio_close(fcst_id,rc) if( hour.eq.hmax ) exit enddo ! end loop for all forecast files print * c Write out correlation and rms data c ---------------------------------- year = nymdb/10000 month = mod(nymdb,10000)/100 day = mod(nymdb,100) hour = nhmsb/10000 minute = mod(nhmsb,10000)/100 write(bdate,3001) nymdb write(bhour,3004) nhmsb/10000 write(edate,3002) nymde write(ehour,3004) nhmse/10000 3001 format('b',i8.8) 3002 format('e',i8.8) 3003 format('x',i6.6) 3004 format('_',i2.2,'z') close(51) datafile = trim(tag) // 'globl.' // bdate // bhour // "." // edate // ehour // ".data" call system ("/bin/mv " // trim(statfile) // " " // trim(datafile) ) statfile = trim(tag) // 'stats.' // bdate // bhour // "." // edate // ehour // ".data" open (85,file=trim(statfile),form='unformatted',access='sequential') ctlfile1 = trim(tag) // 'stats.' // bdate // bhour // "." // edate // ehour // ".ctl1" open (86,file=trim(ctlfile1),form='formatted' ,access='sequential') ctlfile2 = trim(tag) // 'stats.' // bdate // bhour // "." // edate // ehour // ".ctl2" open (87,file=trim(ctlfile2),form='formatted' ,access='sequential') ctlfile3 = trim(tag) // 'globl.' // bdate // bhour // "." // edate // ehour // ".ctl" open (88,file=trim(ctlfile3),form='formatted' ,access='sequential') do ntime =1,nt do nfield=1,nf nlev = nl if(nfield.eq.1) nlev = 1 do level=1,nlev do i=1,nr dum(i) = corr(i,level,nfield,ntime) enddo write(85) dum enddo enddo do m=1,5 do nfield=1,nf nlev = nl if(nfield.eq.1) nlev = 1 do level=1,nlev do i=1,nr dum(i) = rms(i,level,nfield,ntime,m) enddo write(85) dum enddo enddo enddo enddo c Write out correlation tabl ctl1 (uniform times) c ----------------------------------------------- call interp_time ( nymdb,nhmsb, nymde,nhmse, nt, num ) write(86,5001) trim(statfile) write(86,5002) write(86,5003) undef write(86,5004) nr write(86,5005) write(86,5006) nl,(zlev(k),k=1,nl) write(86,5007) nt, months(month),year, num write(86,5008) 6*nf ! Note: 6 Types of Statistical Output per Field (cor, rms, rms_bar, rms_ran, rms_dis, rms_dsp) do n=1,nf lmax = nl if(n.eq.1) lmax = 0 write(86,5009) qname(n),lmax,qdesc(n) enddo do n=nf+1,6*nf lmax = nl if(n.eq. nf+1) lmax = 0 if(n.eq.2*nf+1) lmax = 0 if(n.eq.3*nf+1) lmax = 0 if(n.eq.4*nf+1) lmax = 0 if(n.eq.5*nf+1) lmax = 0 write(86,5009) qname(n),lmax,qdesc(n) enddo write(86,5010) c Write out correlation tabl ctl2 (actual times) c ---------------------------------------------- write(87,5001) trim(statfile) write(87,5002) write(87,5003) undef write(87,5004) nr write(87,5005) write(87,5006) nl,(zlev(k),k=1,nl) write(87,6007) nt,hour,minute,day,months(month),year, num write(87,5008) 6*nf ! Note: 6 Types of Statistical Output per Field (cor, rms, rms_bar, rms_ran, rms_dis, rms_dsp) do n=1,nf lmax = nl if(n.eq.1) lmax = 0 write(87,5009) qname(n),lmax,qdesc(n) enddo do n=nf+1,6*nf lmax = nl if(n.eq. nf+1) lmax = 0 if(n.eq.2*nf+1) lmax = 0 if(n.eq.3*nf+1) lmax = 0 if(n.eq.4*nf+1) lmax = 0 if(n.eq.5*nf+1) lmax = 0 write(87,5009) qname(n),lmax,qdesc(n) enddo write(87,5010) c Write out global data tabl ctl3 (actual times) c ---------------------------------------------- write(88,5001) trim(datafile) write(88,5002) write(88,5003) undef write(88,8004) im, 360.0/float(im) write(88,8005) jm, 180.0/float(jm-1) write(88,5006) nl,(zlev(k),k=1,nl) write(88,6007) nt,hour,minute,day,months(month),year, num write(88,5008) 3*nf ! Note: 3 Sets of Fields (forecast, analysis, climatology) do m=1,3 if( m.eq.1 ) suffix = 'f' if( m.eq.2 ) suffix = 'a' if( m.eq.3 ) suffix = 'c' do n=1,nf lmax = nl if( n.eq.1 ) name = 'p' // trim(suffix) if( n.eq.2 ) name = 'u' // trim(suffix) if( n.eq.3 ) name = 'v' // trim(suffix) if( n.eq.4 ) name = 't' // trim(suffix) if( n.eq.5 ) name = 'q' // trim(suffix) if( n.eq.6 ) name = 'h' // trim(suffix) if( n.eq.1 ) lmax = 0 write(88,5009) trim(name),lmax,trim(name) enddo enddo write(88,5010) 5001 format('DSET ^',a) 5002 format('TITLE Stats',/,'FORMAT sequential big_endian') 5003 format('UNDEF ',g12.6) 5004 format('XDEF ',i3,' LINEAR 1 1') 8004 format('XDEF ',i3,' LINEAR -180 ',f8.5) 5005 format('YDEF 1 LINEAR 1 1') 8005 format('YDEF ',i3,' LINEAR -90 ',f8.5) 5006 format('ZDEF ',i3,' LEVELS ',5(f7.2,1x),/, . 50(18x, 5(f7.2,1x),/) ) 5007 format('TDEF ',i3,' LINEAR ', '00:00', 'Z01' ,a3,i4,' ',i2.2,'hr') 6007 format('TDEF ',i3,' LINEAR ',i2.2,':',i2.2,'Z',i2.2,a3,i4,' ',i2.2,'hr') 5008 format('VARS ',i3) 5009 format(a8,2x,i3,' 0 ',a) 5010 format('ENDVARS') c Write Timing Information c ------------------------ call timeend ('main') call timepri (6) stop end SUBROUTINE BOUNDS (LAT1,LAT2,LONG1,LONG2, . JBEG,JEND,IBEG,IEND,IM,JM) real lat1,lat2,long1,long2 IMP1 = IM+1 PI = 180. DL = 2.*PI/IM DP = PI/(JM-1) jbeg = int((lat1 + 90.)/dp)+1 jend = int((lat2 + 90.)/dp)+1 IF( JBEG.GT.JM) JBEG = JM IF( (JEND-1)*DP-90.LT.LAT2) JEND= JEND+1 IF( JEND.GT.JM) JEND = JM ibeg = int((long1+180.)/dl)+1 iend = int((long2+180.)/dl)+1 IF( IBEG.GT.IM) IBEG = IM IF((IEND-1)*DL-180.LT.LONG2) IEND=IEND+1 IF( IEND.GT.IM) IEND = IM c write(6,100) lat1 ,jbeg, lat2,jend, c . long1,ibeg,long2,iend 100 format(/1x,f6.1,' (',i3,')',2x,f6.1,' (',i3,')',/, . 1x,f6.1,' (',i3,')',2x,f6.1,' (',i3,')',/) return end subroutine read_clim_hdf ( nymd,nhms,pc,uc,vc,tc,qc,hc,idim,jdim,nl,zlev,cli_files,num,undef ) implicit none logical match character*256, parameter :: c_slp = 'slp' character*256, parameter :: c_u = 'u' character*256, parameter :: c_v = 'v' character*256, parameter :: c_t = 't' character*256, parameter :: c_q = 'q' character*256, parameter :: c_h = 'h' integer num integer nymd,nhms integer idim,jdim,nl real pc(idim,jdim) real uc(idim,jdim,nl) real vc(idim,jdim,nl) real tc(idim,jdim,nl) real qc(idim,jdim,nl) real hc(idim,jdim,nl) real zlev(nl) real, allocatable, save :: p1( :,: ) real, allocatable, save :: u1( :,:,: ) real, allocatable, save :: v1( :,:,: ) real, allocatable, save :: t1( :,:,: ) real, allocatable, save :: q1( :,:,: ) real, allocatable, save :: h1( :,:,: ) real, allocatable, save :: p2( :,: ) real, allocatable, save :: u2( :,:,: ) real, allocatable, save :: v2( :,:,: ) real, allocatable, save :: t2( :,:,: ) real, allocatable, save :: q2( :,:,: ) real, allocatable, save :: h2( :,:,: ) integer nid,ncl,im,jm,lm,nvars,rc integer ntime,ngatts,timinc real undef integer L,m,n character*256 cli_files(num) character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) real, allocatable :: q(:,:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: kmvar(:) integer, allocatable, save :: id(:) integer, allocatable, save :: yymmdd(:,:) integer, allocatable, save :: hhmmss(:,:) logical first, found, shift, defined integer LL, i,j,k, loc data shift /.false./ data first /.true./ save integer ISCM, ISCP, ISC real FACM, FACP integer MIDMON integer IMNP, IMNM integer IMONP, IMONM DATA IMONP/0/, IMONM/0/ INTEGER YEAR, MONTH, DAY, SEC INTEGER DAYS(12) DATA DAYS /31,28,31,30,31,30,31,31,30,31,30,31/ INTEGER NSECF, NMONF, NDAYF NSECF(N) = N/10000*3600 + MOD(N,10000)/100* 60 + MOD(N,100) NMONF(N) = MOD(N,10000)/100 NDAYF(N) = MOD(N,100) C********************************************************************* C**** Find Proper Month Boundaries from INPUT Date and Time **** C********************************************************************* SEC = NSECF(NHMS) MONTH = NMONF(NYMD) DAY = NDAYF(NYMD) MIDMON = DAYS(MONTH)/2 + 1 IF(DAY.LT.MIDMON) THEN imnm = month - 1 imnp = month ELSE IMNM = MONTH IMNP = MONTH + 1 ENDIF if( imnm.eq.0 ) imnm = 12 if( imnp.eq.13 ) imnp = 1 C********************************************************************* C**** Open Climatology DataSet and Initialization **** C********************************************************************* if( first ) then allocate ( yymmdd(12,num) ) allocate ( hhmmss(12,num) ) allocate ( id(num) ) do n=1,num call gfio_open ( cli_files(n),1,id(n),rc ) if( rc.ne.0 ) then print *, 'Climatology File: ',trim(cli_files(n)),' NOT found!' stop endif call gfio_diminquire ( id(n),im,jm,lm,ntime,nvars,ngatts,rc ) if( ntime.ne.12 ) then print *, 'Climatology data should consist of 12 monthly means' print *, 'Current file: ',trim(cli_files(n)),' contains ',ntime,' time periods!' stop endif if( first ) then allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) endif call gfio_inquire ( id(n),im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd(1,n),hhmmss(1,n),timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) first = .false. enddo if( lon(1).eq.0.0 ) then print *, 'Climatology data begins at lon: ',lon(1) print *, 'Horizontal Shift will be performed' print * shift = .true. endif allocate ( p1(idim,jdim ), p2(idim,jdim ) ) allocate ( u1(idim,jdim,nl), u2(idim,jdim,nl) ) allocate ( v1(idim,jdim,nl), v2(idim,jdim,nl) ) allocate ( t1(idim,jdim,nl), t2(idim,jdim,nl) ) allocate ( q1(idim,jdim,nl), q2(idim,jdim,nl) ) allocate ( h1(idim,jdim,nl), h2(idim,jdim,nl) ) endif pc = undef uc = undef vc = undef tc = undef qc = undef hc = undef c Determine Climatology TOD Dataset c --------------------------------- if( num.eq.1 ) then nid = id(num) ncl = num else nid = -999 do n=1,num if( nhms.eq.hhmmss(1,n) ) then nid = id(n) ncl = n endif enddo if( nid.eq.-999 ) then print *, 'Climatology Datasets do not have desired TOD: ',nhms stop endif endif C********************************************************************* C**** Read for (-) Month and (+) Month **** C********************************************************************* IMONM = IMNM IMONP = IMNP allocate ( q(im,jm) ) do k=1,ntime MONTH = NMONF(yymmdd(k,ncl)) if( month.eq.imnm .or. month.eq.imnp ) then do n=1,nvars if( match( c_slp,vname(n) ) ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,0,1,q,rc ) if( shift ) call hshift ( q,im,jm ) do j=1,jm do i=1,im if( defined(q(i,j),undef) .and. . q(i,j).gt.10000.0 ) then q(i,j) = q(i,j)*0.01 ! Convert Pa=>mb endif enddo enddo call bin ( q,im,jm,pc,idim,jdim,undef,0 ) endif do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( match( c_u,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,uc(1,1,L),idim,jdim,undef,1 ) endif if( match( c_v,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,vc(1,1,L),idim,jdim,undef,1 ) endif if( match( c_t,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,tc(1,1,L),idim,jdim,undef,0 ) endif if( match( c_q,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,qc(1,1,L),idim,jdim,undef,0 ) endif if( match( c_h,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( nid,vname(n),yymmdd(k,ncl),hhmmss(k,ncl),im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,hc(1,1,L),idim,jdim,undef,0 ) endif enddo enddo if( month.eq.imnm ) then p1 = pc u1 = uc v1 = vc t1 = tc q1 = qc h1 = hc endif if( month.eq.imnp ) then p2 = pc u2 = uc v2 = vc t2 = tc q2 = qc h2 = hc endif endif enddo deallocate ( q ) C********************************************************************* C**** INTERPOLATE DATA TO CURRENT TIME **** C********************************************************************* IF(DAY.LT.MIDMON) THEN ISC = (DAY+DAYS(IMONM))*86400 + SEC ELSE ISC = (DAY )*86400 + SEC ENDIF ISCM = (DAYS(IMONM)/2+1 )*86400 ISCP = (DAYS(IMONP)/2+1 + DAYS(IMONM) )*86400 FACP = (ISC-ISCM) / REAL(ISCP-ISCM) FACM = 1-FACP c write(6,100) imonm,facm,imonp,facp c100 format(1x,'Interpolating Climatology for Month ',i2, c . ' (',f5.3,') and Month ',i2,' (',f5.3,')') do J=1,JDIM do I=1,IDIM if( p1(i,j).ne.undef .and. p2(i,j).ne.undef ) then pc(I,J) = p1(I,J)*FACM + p2(I,J)*FACP else if( p1(i,j).ne.undef ) then pc(I,J) = p1(I,J) else if( p2(i,j).ne.undef ) then pc(I,J) = p2(I,J) else pc(I,J) = undef endif endif endif enddo enddo do L=1,NL do J=1,JDIM do I=1,IDIM if( u1(i,j,L).ne.undef .and. u2(i,j,L).ne.undef ) then uc(I,J,L) = u1(I,J,L)*FACM + u2(I,J,L)*FACP else if( u1(i,j,L).ne.undef ) then uc(I,J,L) = u1(I,J,L) else if( u2(i,j,L).ne.undef ) then uc(I,J,L) = u2(I,J,L) else uc(I,J,L) = undef endif endif endif if( v1(i,j,L).ne.undef .and. v2(i,j,L).ne.undef ) then vc(I,J,L) = v1(I,J,L)*FACM + v2(I,J,L)*FACP else if( v1(i,j,L).ne.undef ) then vc(I,J,L) = v1(I,J,L) else if( v2(i,j,L).ne.undef ) then vc(I,J,L) = v2(I,J,L) else vc(I,J,L) = undef endif endif endif if( t1(i,j,L).ne.undef .and. t2(i,j,L).ne.undef ) then tc(I,J,L) = t1(I,J,L)*FACM + t2(I,J,L)*FACP else if( t1(i,j,L).ne.undef ) then tc(I,J,L) = t1(I,J,L) else if( t2(i,j,L).ne.undef ) then tc(I,J,L) = t2(I,J,L) else tc(I,J,L) = undef endif endif endif if( q1(i,j,L).ne.undef .and. q2(i,j,L).ne.undef ) then qc(I,J,L) = q1(I,J,L)*FACM + q2(I,J,L)*FACP else if( q1(i,j,L).ne.undef ) then qc(I,J,L) = q1(I,J,L) else if( q2(i,j,L).ne.undef ) then qc(I,J,L) = q2(I,J,L) else qc(I,J,L) = undef endif endif endif if( h1(i,j,L).ne.undef .and. h2(i,j,L).ne.undef ) then hc(I,J,L) = h1(I,J,L)*FACM + h2(I,J,L)*FACP else if( h1(i,j,L).ne.undef ) then hc(I,J,L) = h1(I,J,L) else if( h2(i,j,L).ne.undef ) then hc(I,J,L) = h2(I,J,L) else hc(I,J,L) = undef endif endif endif enddo enddo enddo return end subroutine read_clim_bin ( nymd,nhms,p,u,v,t,q,h,idim,jdim,ldim,undef ) C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C* Note: Climatology Data is in Grads Format from the files: * C* ncep_1x1_clim.ctl * C* ncep_1x1_clim.data * C* Climatology Data is stored: January through December * C*********************************************************************** PARAMETER ( IM = 360 ) PARAMETER ( JM = 181 ) PARAMETER ( LM = 10 ) real*4 bum(IM,JM) real dum(IM,JM) real p ( IDIM,JDIM ) real u ( IDIM,JDIM,LDIM ) real v ( IDIM,JDIM,LDIM ) real t ( IDIM,JDIM,LDIM ) real q ( IDIM,JDIM,LDIM ) real h ( IDIM,JDIM,LDIM ) real, allocatable, save :: p1( :,: ) real, allocatable, save :: u1( :,:,: ) real, allocatable, save :: v1( :,:,: ) real, allocatable, save :: t1( :,:,: ) real, allocatable, save :: q1( :,:,: ) real, allocatable, save :: h1( :,:,: ) real, allocatable, save :: p2( :,: ) real, allocatable, save :: u2( :,:,: ) real, allocatable, save :: v2( :,:,: ) real, allocatable, save :: t2( :,:,: ) real, allocatable, save :: q2( :,:,: ) real, allocatable, save :: h2( :,:,: ) integer LL real undef logical first data first /.true./ DATA IMONP/0/, IMONM/0/, ku /90/ SAVE imonp, imonm INTEGER DAYS(12) DATA DAYS /31,28,31,30,31,30,31,31,30,31,30,31/ NSECF(N) = N/10000*3600 + MOD(N,10000)/100* 60 + MOD(N,100) NMONF(N) = MOD(N,10000)/100 NDAYF(N) = MOD(N,100) if( first ) then open (90,file='ncep_1x1_clim.data', . form='unformatted',access='direct',recl=im*jm ) c . form='unformatted',access='direct',recl=im*jm*4, c . convert='little_endian') allocate ( p1(idim,jdim) , p2(idim,jdim) ) allocate ( u1(idim,jdim,ldim), u2(idim,jdim,ldim) ) allocate ( v1(idim,jdim,ldim), v2(idim,jdim,ldim) ) allocate ( t1(idim,jdim,ldim), t2(idim,jdim,ldim) ) allocate ( q1(idim,jdim,ldim), q2(idim,jdim,ldim) ) allocate ( h1(idim,jdim,ldim), h2(idim,jdim,ldim) ) first = .false. endif undef = 1e15 C********************************************************************* C**** FIND PROPER MONTH BOUNDARIES **** C********************************************************************* SEC = NSECF(NHMS) MONTH = NMONF(NYMD) DAY = NDAYF(NYMD) MIDMON = DAYS(MONTH)/2 + 1 IF(DAY.LT.MIDMON) THEN imnm = month - 1 imnp = month ELSE IMNM = MONTH IMNP = MONTH + 1 ENDIF if( imnm.eq.0 ) imnm = 12 if( imnp.eq.13 ) imnp = 1 C********************************************************************* C**** READ DATA SET TWICE TO GET BOTH MONTHS **** C********************************************************************* IF(IMONM.NE.IMNM .OR. IMONP.NE.IMNP) THEN IMONM = IMNM IMONP = IMNP c Read First Month c ---------------- month1 = imonm if( month1.lt.1 ) month1 = month1 + 12 n = (month1-1)*(1+LM*5) read(ku,rec=n+1) bum dum = bum call bin ( dum,im,jm,p1,idim,jdim,undef,0 ) do L=1,LM read(ku,rec=n+1+L) bum dum = bum call bin ( dum,im,jm,u1(1,1,L),idim,jdim,undef,1 ) enddo do L=1,LM read(ku,rec=n+1+L+LM) bum dum = bum call bin ( dum,im,jm,v1(1,1,L),idim,jdim,undef,1 ) enddo do L=1,LM read(ku,rec=n+1+L+2*LM) bum dum = bum call bin ( dum,im,jm,t1(1,1,L),idim,jdim,undef,0 ) enddo do L=1,LM read(ku,rec=n+1+L+3*LM) bum dum = bum call bin ( dum,im,jm,q1(1,1,L),idim,jdim,undef,0 ) enddo do L=1,LM read(ku,rec=n+1+L+4*LM) bum dum = bum call bin ( dum,im,jm,h1(1,1,L),idim,jdim,undef,0 ) enddo c Read Second Month c ----------------- month2 = imonp if( month2.lt.1 ) month2 = month2 + 12 n = (month2-1)*(1+LM*5) read(ku,rec=n+1) bum dum = bum call bin ( dum,im,jm,p2,idim,jdim,undef,0 ) do L=1,LM read(ku,rec=n+1+L) bum dum = bum call bin ( dum,im,jm,u2(1,1,L),idim,jdim,undef,1 ) enddo do L=1,LM read(ku,rec=n+1+L+LM) bum dum = bum call bin ( dum,im,jm,v2(1,1,L),idim,jdim,undef,1 ) enddo do L=1,LM read(ku,rec=n+1+L+2*LM) bum dum = bum call bin ( dum,im,jm,t2(1,1,L),idim,jdim,undef,0 ) enddo do L=1,LM read(ku,rec=n+1+L+3*LM) bum dum = bum call bin ( dum,im,jm,q2(1,1,L),idim,jdim,undef,0 ) enddo do L=1,LM read(ku,rec=n+1+L+4*LM) bum dum = bum call bin ( dum,im,jm,h2(1,1,L),idim,jdim,undef,0 ) enddo ENDIF C********************************************************************* C**** INTERPOLATE DATA TO CURRENT TIME **** C********************************************************************* IF(DAY.LT.MIDMON) THEN ISC = (DAY+DAYS(IMONM))*86400 + SEC ELSE ISC = (DAY )*86400 + SEC ENDIF ISCM = (DAYS(IMONM)/2+1 )*86400 ISCP = (DAYS(IMONP)/2+1 + DAYS(IMONM) )*86400 FACP = (ISC-ISCM) / REAL(ISCP-ISCM) FACM = 1-FACP c write(6,100) imonm,facm,imonp,facp c100 format(1x,'Interpolating Climatology for Month ',i2, c . ' (',f5.3,') and Month ',i2,' (',f5.3,')') do J=1,JDIM do I=1,IDIM p(I,J) = p1(I,J)*FACM + p2(I,J)*FACP enddo enddo do L=1,LDIM do J=1,JDIM do I=1,IDIM u(I,J,L) = u1(I,J,L)*FACM + u2(I,J,L)*FACP v(I,J,L) = v1(I,J,L)*FACM + v2(I,J,L)*FACP t(I,J,L) = t1(I,J,L)*FACM + t2(I,J,L)*FACP q(I,J,L) = q1(I,J,L)*FACM + q2(I,J,L)*FACP h(I,J,L) = h1(I,J,L)*FACM + h2(I,J,L)*FACP enddo enddo enddo RETURN END function defined ( q,undef ) implicit none logical defined real q,undef defined = abs(q-undef).gt.0.1*abs(undef) return end subroutine hshift ( q,im,jm ) real q(im,jm), dum(im,jm) dum(1:im/2,:) = q(1:im/2,:) q(1:im/2,:) = q(1+im/2:im,:) q(1+im/2:im,:) = dum(1:im/2,:) return end subroutine minmax (q,im,jm,undef) real q(im,jm) qmin = q(1,1) qmax = q(1,1) do j=1,jm do i=1,im if(q(i,j).ne.undef) qmin = min( qmin,q(i,j) ) if(q(i,j).ne.undef) qmax = max( qmax,q(i,j) ) enddo enddo print *, ' qmin: ',qmin,' qmax: ',qmax return end subroutine bin ( qin,im_in,jm_in,qout,im_out,jm_out,undef,msgn ) implicit none integer im_in ,jm_in ,msgn integer im_out,jm_out real undef real qin(im_in ,jm_in ) real qout(im_out,jm_out) real q10x10(360*6,180*6) c Parse Arbitray Field (im,jm) to 10'x10' Variable c ------------------------------------------------ call timebeg ('bin') call bin_10x10 ( qin,im_in,jm_in,q10x10 ) c Bin 10'x10' Variable to Output Field (im_out,jm_out) c ---------------------------------------------------- call averaged_10x10 ( q10x10,qout,im_out,jm_out,undef,msgn ) call timeend ('bin') return end subroutine averaged_10x10 ( z10x10,z,im,jm,undef,msgn ) C*********************************************************************** C C PURPOSE: C ======== C Average a (10m X 10m) input array to an output array (im,jm) C C INPUT: C ====== C z10x10 ..... Input array(360*6,180*6) C msgn ....... Integer Flag for scalar (0) or vector (1) C C OUTPUT: C ======= C z .......... Output array(im,jm) C im ......... Longitudinal dimension of z C jm ......... Latitudinal dimension of z C C NOTES: C ====== C Input array z10x10 represents values within a 10min X 10min grid-box. C Each box is referenced by the latitude and longitude of C its southwest corner, not its center point. Thus, C the height associated with a coordinate actually C represents the heights centered to the northeast of that point. C C Output array z(im,jm) is assumed to be on an A-grid. C z(i,j) represents the value at the center of the grid-box. C z(1,j) is located at lon=-180. C z(i,1) is located at lat=-90. C z(i,jm) is located at lat=+90. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer im,jm,msgn real z(im,jm) real dlam(im), dphi(jm) real z10x10(360*6,180*6) integer i,j,ibeg,iend,jbeg,jend integer ii,jj,itmp real sum1,sum2 real zlat,zlon real lon1,lon2,wx real lat1,lat2,wy real lonbeg,lonend,lat,coslat real latbeg,latend real undef real pi,dz real lon_cmp(im) real lat_cmp(jm) logical defined call timebeg ('ave_10x10') pi = 4.*atan(1.) dz = pi/(6.*180) dlam = 2*pi/ im dphi = pi/(jm-1) c Compute Computational Lambda's and Phi's c ---------------------------------------- lon_cmp(1) = -pi do i=2,im lon_cmp(i) = lon_cmp(i-1) + dlam(i-1) enddo lat_cmp(1) = -pi*0.5 do j=2,jm-1 lat_cmp(j) = lat_cmp(j-1) + dphi(j-1) enddo lat_cmp(jm) = pi*0.5 c Compute average away from poles c ------------------------------- do j=2,jm-1 do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat-dphi(j-1)/2 latend = zlat+dphi(j) /2 if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1.+(latbeg+pi/2)/dz jend = 1.+(latend+pi/2)/dz sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo enddo c Compute average at South Pole c ----------------------------- j=1 do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat latend = zlat+dphi(j) /2 if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1 jend = 1.+(latend+pi/2)/dz sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo c Compute average at North Pole c ----------------------------- j=jm do i=1,im zlat = lat_cmp(j) zlon = lon_cmp(i) latbeg = zlat-dphi(j-1)/2 latend = zlat if( i.eq.1 ) then lonbeg = zlon-dlam(im) /2 else lonbeg = zlon-dlam(i-1)/2 endif lonend = zlon+dlam(i) /2 ibeg = 1.+(lonbeg+pi) /dz iend = 1.+(lonend+pi) /dz jbeg = 1.+(latbeg+pi/2)/dz jend = 1080 sum1 = 0 sum2 = 0 do jj=jbeg,jend lat = -pi/2+(jj-0.5)*dz coslat = cos(lat) lat1 = -pi/2 + (jj-1)*dz lat2 = -pi/2 + jj *dz wy = 1.0 if( lat1.lt.latbeg ) wy = (lat2-latbeg)/dz if( lat2.gt.latend ) wy = (latend-lat1)/dz if(ibeg.ge.1) then do ii=ibeg,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo else itmp = 1.+(lonbeg+0.1*dz+3*pi)/dz do ii=itmp,360*6 if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg+2*pi ) wx = (lon2-lonbeg-2*pi)/dz if( lon2.gt.lonend+2*pi ) wx = (2*pi+lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo do ii=1,iend if( defined(z10x10(ii,jj),undef) ) then lon1 = -pi + (ii-1)*dz lon2 = -pi + ii *dz wx = 1.0 if( lon1.lt.lonbeg ) wx = (lon2-lonbeg)/dz if( lon2.gt.lonend ) wx = (lonend-lon1)/dz sum1 = sum1 + z10x10(ii,jj)*coslat*wx*wy sum2 = sum2 + coslat*wx*wy endif enddo endif enddo if( sum2.ne.0.0 ) then z(i,j) = sum1/sum2 else z(i,j) = undef endif enddo c Average Pole Values c ------------------- if( msgn.eq.0 ) then sum1 = 0 j = 0 do i=1,im if( defined(z(i,1),undef) ) then sum1 = sum1 + z(i,1) j = j + 1 endif enddo if( j.ne.0 ) then z(:,1) = sum1/j else z(:,1) = undef endif sum2 = 0 j = 0 do i=1,im if( defined(z(i,jm),undef) ) then sum2 = sum2 + z(i,jm) j = j + 1 endif enddo if( j.ne.0 ) then z(:,jm) = sum2/j else z(:,jm) = undef endif endif call timeend ('ave_10x10') return end subroutine bin_10x10 ( z,im,jm,z10x10 ) C*********************************************************************** C C PURPOSE: C ======== C Compute a (10m X 10m) array binned from an input array (im,jm) C C INPUT: C ====== C z .......... Input array(im,jm) C im ......... Longitudinal dimension of z C jm ......... Latitudinal dimension of z C C OUTPUT: C ======= C z10x10 ..... Output array(360*6,180*6) C C NOTES: C ====== C Input array z(im,jm) is assumed to be on an A-grid. C z(i,j) represents the value at the center of the grid-box. C z(1,j) is located at lon=-180. C z(i,1) is located at lat=-90. C z(i,jm) is located at lat=+90. C C Output array z10x10 represents values within a 10min X 10min grid-box. C Each box is referenced by the latitude and longitude of C its southwest corner, not its center point. Thus, C the height associated with a coordinate actually C represents the heights centered to the northeast of that point. C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer im,jm real z(im,jm) real z10x10(360*6,180*6) integer i,j,ii,jj,ibeg,iend,jbeg,jend real zlatc,zlonc real lonbeg,lonend,lat real latbeg,latend real pi,dl,dp,dz call timebeg ('bin_10x10') pi = 4.*atan(1.) dl = 2*pi/im dp = pi/(jm-1) dz = pi/(6.*180) do j=1,180*6 do i=1,360*6 zlatc = -pi/2+(j-0.5)*dz ! Latitude at center of 10x10 box zlonc = -pi +(i-0.5)*dz ! Longitude at center of 10x10 box c Find bounding lat and lon on IMxJM grid c --------------------------------------- iend = nint( 1.+(zlonc+pi)/dl ) lonend = -pi + (iend-1)*dl if( lonend.ge.zlonc ) then lonbeg = -pi + (iend-2)*dl else iend = iend+1 lonbeg = lonend lonend = -pi + (iend-1)*dl endif ibeg = iend-1 jend = nint( 1.+(zlatc+pi/2)/dp ) latend = -pi/2 + (jend-1)*dp if( latend.ge.zlatc ) then latbeg = -pi/2 + (jend-2)*dp else jend = jend+1 latbeg = latend latend = -pi/2 + (jend-1)*dp endif jbeg = jend-1 if(iend.gt.im) iend=iend-im if( zlonc.le.lonbeg+0.5*dl ) then ii = ibeg else ii = iend endif if( zlatc.le.latbeg+0.5*dp ) then jj = jbeg else jj = jend endif z10x10(i,j) = z(ii,jj) enddo enddo call timeend ('bin_10x10') return end subroutine read_anal( nymd,nhms,pa,ua,va,ta,qa,ha,idim,jdim,nl,zlev,ana_files,num_ana_files,undef ) implicit none logical match character*256, parameter :: c_slp = 'slp' character*256, parameter :: c_ps = 'ps' character*256, parameter :: c_dp = 'dp' character*256, parameter :: c_pl = 'pl' character*256, parameter :: c_ple = 'ple' character*256, parameter :: c_u = 'u' character*256, parameter :: c_v = 'v' character*256, parameter :: c_t = 't' character*256, parameter :: c_q = 'q' character*256, parameter :: c_h = 'h' character*256, parameter :: c_th = 'th' character*256, parameter :: c_thv = 'thv' character*256, parameter :: c_phis = 'phis' integer nymd,nhms integer idim,jdim,nl,num_ana_files real pa(idim,jdim) real ua(idim,jdim,nl) real va(idim,jdim,nl) real ta(idim,jdim,nl) real qa(idim,jdim,nl) real ha(idim,jdim,nl) real zlev(nl) integer id,im,jm,lm,nvars,rc integer ntime,ngatts,timinc real undef integer L,m,n character*256 ana_files(num_ana_files) character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) real, allocatable :: q(:,:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: yymmdd(:) integer, allocatable :: hhmmss(:) integer, allocatable :: kmvar(:) integer, allocatable :: dates(:,:) integer, allocatable :: datez(:,:) logical found_slp logical found_uwnd logical found_vwnd logical found_hght logical found_tmpu logical found_sphu logical shift, defined, first integer num, LL, i,j, loc, ndates data id /0/ data num /0/ data shift /.false./ data first /.true./ save c Read All ANA Files to Gather Date and Time Information c ------------------------------------------------------ if( first ) then call timebeg ('exam_anal') print *, 'Examining ANA Files ...' ndates = 0 do num=1,num_ana_files call gfio_open ( ana_files(num),1,id,rc ) call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( yymmdd(ntime) ) allocate ( hhmmss(ntime) ) allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) if( ndates.eq.0 ) then ndates = ndates + ntime allocate ( dates(3,ndates) ) allocate ( datez(3,ndates) ) dates(1,ndates-ntime+1:ndates) = yymmdd(:) dates(2,ndates-ntime+1:ndates) = hhmmss(:) dates(3,ndates-ntime+1:ndates) = num datez = dates else deallocate(dates) ndates = ndates + ntime allocate ( dates(3,ndates) ) dates(1,1:ndates-ntime) = datez(1,1:ndates-ntime) dates(2,1:ndates-ntime) = datez(2,1:ndates-ntime) dates(3,1:ndates-ntime) = datez(3,1:ndates-ntime) dates(1,ndates-ntime+1:ndates) = yymmdd(:) dates(2,ndates-ntime+1:ndates) = hhmmss(:) dates(3,ndates-ntime+1:ndates) = num deallocate(datez) allocate ( datez(3,ndates) ) datez = dates endif call gfio_close(id,rc) deallocate ( lon,lat,lev,yymmdd,hhmmss,vname,vtitle,vunits,kmvar,vrange,prange ) enddo first = .false. call timeend ('exam_anal') endif c Read Appropriate ANA Files to Get Data c -------------------------------------- call timebeg ('read_anal') found_slp = .false. found_uwnd = .false. found_vwnd = .false. found_hght = .false. found_tmpu = .false. found_sphu = .false. do num=1,ndates if( dates(1,num).eq.nymd .and. dates(2,num).eq.nhms ) then call gfio_open ( ana_files(dates(3,num)),1,id,rc ) call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( yymmdd(ntime) ) allocate ( hhmmss(ntime) ) allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) shift = .false. if( lon(1).eq.0.0 ) then print *, 'Analysis data begins at lon: ',lon(1) print *, 'Horizontal Shift will be performed' print * shift = .true. endif allocate ( q(im,jm) ) do n=1,nvars if( .not.found_slp .and. match( c_slp,vname(n) ) ) then found_slp = .true. call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,0,1,q,rc ) if( shift ) call hshift ( q,im,jm ) do j=1,jm do i=1,im if( defined(q(i,j),undef) .and. . q(i,j).gt.10000.0 ) then q(i,j) = q(i,j)*0.01 ! Convert Pa=>mb endif enddo enddo call bin ( q,im,jm,pa,idim,jdim,undef,0 ) endif if( .not.found_uwnd .and. match( c_u,vname(n) ) ) then found_uwnd = .true. do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ua(1,1,L),idim,jdim,undef,1 ) enddo endif if( .not.found_vwnd .and. match( c_v,vname(n) ) ) then found_vwnd = .true. do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,va(1,1,L),idim,jdim,undef,1 ) enddo endif if( .not.found_tmpu .and. match( c_t,vname(n) ) ) then found_tmpu = .true. do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ta(1,1,L),idim,jdim,undef,0 ) enddo endif if( .not.found_sphu .and. match( c_q,vname(n) ) ) then found_sphu = .true. do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,qa(1,1,L),idim,jdim,undef,0 ) enddo endif if( .not.found_hght .and. match( c_h,vname(n) ) ) then found_hght = .true. do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ha(1,1,L),idim,jdim,undef,0 ) enddo endif enddo call gfio_close(id,rc) deallocate ( q,lon,lat,lev,yymmdd,hhmmss,vname,vtitle,vunits,kmvar,vrange,prange ) endif enddo if( .not.found_slp ) pa = undef if( .not.found_uwnd ) ua = undef if( .not.found_vwnd ) va = undef if( .not.found_tmpu ) ta = undef if( .not.found_sphu ) qa = undef if( .not.found_hght ) ha = undef call timeend ('read_anal') return end subroutine init_fcst( fname,id,undef,ntime,yymmdd,hhmmss,timinc ) implicit none integer id,im,jm,lm,nvars,rc integer ntime,ngatts,timinc real undef integer L,m,n character*256 fname integer yymmdd(100) integer hhmmss(100) character*256 vname(500) character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: kmvar(:) call timebeg ('init_fcst') call gfio_open ( fname,1,id,rc ) if( rc.ne.0 ) then print *, 'Forecast File: ',trim(fname),' NOT found!' stop endif call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) call timeend ('init_fcst') return end subroutine read_fcst( id,nymd,nhms,pa,ua,va,ta,qa,ha,idim,jdim,nl,zlev ) implicit none logical match, defined character*256, parameter :: c_slp = 'slp' character*256, parameter :: c_ps = 'ps' character*256, parameter :: c_dp = 'dp' character*256, parameter :: c_pl = 'pl' character*256, parameter :: c_ple = 'ple' character*256, parameter :: c_u = 'u' character*256, parameter :: c_v = 'v' character*256, parameter :: c_t = 't' character*256, parameter :: c_q = 'q' character*256, parameter :: c_h = 'h' character*256, parameter :: c_th = 'th' character*256, parameter :: c_thv = 'thv' character*256, parameter :: c_phis = 'phis' integer nymd,nhms integer idim,jdim,nl real pa(idim,jdim) real ua(idim,jdim,nl) real va(idim,jdim,nl) real ta(idim,jdim,nl) real qa(idim,jdim,nl) real ha(idim,jdim,nl) real zlev(nl) integer id,im,jm,lm,nvars,rc real undef integer i,j,L,m,n integer LL,loc logical shift real, allocatable :: q(:,:) character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) character*256, allocatable :: vname (:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: kmvar(:) integer, allocatable :: yymmdd(:) integer, allocatable :: hhmmss(:) integer ntime,ngatts,timinc call timebeg ('read_fcst') call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( yymmdd(ntime) ) allocate ( hhmmss(ntime) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( vname(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) shift = .false. c if( lon(1).ne.-180.0 ) then if( lon(1).eq.0.0 ) then print *, 'Forecast data begins at lon: ',lon(1) print *, 'Horizontal Shift will be performed' print * shift = .true. endif allocate ( q(im,jm) ) pa = undef ua = undef va = undef ta = undef qa = undef ha = undef do n=1,nvars if( match( c_slp,vname(n) ) ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,0,1,q,rc ) if( shift ) call hshift ( q,im,jm ) do j=1,jm do i=1,im if( defined(q(i,j),undef) .and. . q(i,j).gt.10000.0 ) then q(i,j) = q(i,j)*0.01 ! Convert Pa=>mb endif enddo enddo call bin ( q,im,jm,pa,idim,jdim,undef,0 ) endif do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( match( c_u,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ua(1,1,L),idim,jdim,undef,1 ) endif if( match( c_v,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,va(1,1,L),idim,jdim,undef,1 ) endif if( match( c_t,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ta(1,1,L),idim,jdim,undef,0 ) endif if( match( c_q,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,qa(1,1,L),idim,jdim,undef,0 ) endif if( match( c_h,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,trim(vname(n)),nymd,nhms,im,jm,loc,1,q,rc ) else q = undef endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ha(1,1,L),idim,jdim,undef,0 ) endif enddo enddo deallocate ( q ) deallocate ( lon,lat,lev,yymmdd,hhmmss,vname,vtitle,vunits,kmvar,vrange,prange ) call timeend ('read_fcst') return end subroutine read_syserr( pa,ua,va,ta,qa,ha,idim,jdim,nl,zlev,efile,ndt,undef ) implicit none logical match character*256, parameter :: c_slp = 'slp' character*256, parameter :: c_ps = 'ps' character*256, parameter :: c_dp = 'dp' character*256, parameter :: c_pl = 'pl' character*256, parameter :: c_ple = 'ple' character*256, parameter :: c_u = 'u' character*256, parameter :: c_v = 'v' character*256, parameter :: c_t = 't' character*256, parameter :: c_q = 'q' character*256, parameter :: c_h = 'h' character*256, parameter :: c_th = 'th' character*256, parameter :: c_thv = 'thv' character*256, parameter :: c_phis = 'phis' integer nymd ,nhms integer nymd0,nhms0 integer idim,jdim,nl real pa(idim,jdim) real ua(idim,jdim,nl) real va(idim,jdim,nl) real ta(idim,jdim,nl) real qa(idim,jdim,nl) real ha(idim,jdim,nl) real zlev(nl) integer id,im,jm,lm,nvars,rc integer ntime,ngatts,timinc integer nsecf real undef integer L,m,n integer ndt,num,nt character*256 efile character*256 title character*256 source character*256 contact character*256 levunits character*256, allocatable :: vname(:) character*256, allocatable :: vtitle(:) character*256, allocatable :: vunits(:) real, allocatable :: q(:,:) real, allocatable :: lat(:) real, allocatable :: lon(:) real, allocatable :: lev(:) real, allocatable :: vrange(:,:) real, allocatable :: prange(:,:) integer, allocatable :: yymmdd(:) integer, allocatable :: hhmmss(:) integer, allocatable :: kmvar(:) logical found, shift, defined integer LL, i,j, loc data id /0/ data num /0/ data shift /.false./ save 10 continue if( id.eq.0 ) then call gfio_open ( efile,1,id,rc ) if( rc.ne.0 ) then print *, 'Systematic Error File: ',trim(efile),' NOT found!' stop endif call gfio_diminquire ( id,im,jm,lm,ntime,nvars,ngatts,rc ) allocate ( lon(im) ) allocate ( lat(jm) ) allocate ( lev(lm) ) allocate ( yymmdd(ntime) ) allocate ( hhmmss(ntime) ) allocate ( vname(nvars) ) allocate ( vtitle(nvars) ) allocate ( vunits(nvars) ) allocate ( kmvar(nvars) ) allocate ( vrange(2,nvars) ) allocate ( prange(2,nvars) ) call gfio_inquire ( id,im,jm,lm,ntime,nvars, . title,source,contact,undef, . lon,lat,lev,levunits, . yymmdd,hhmmss,timinc, . vname,vtitle,vunits,kmvar, . vrange,prange,rc ) if( ndt .ne. nsecf(timinc) ) then print *, 'Error!' print *, 'Forecast File Frequency: ',ndt,' (sec)' print *, 'Systematic Error File Frequency: ',nsecf(timinc),' (sec)' stop endif if( lon(1).eq.0.0 ) then print *, 'Systematic Error data begins at lon: ',lon(1) print *, 'Horizontal Shift will be performed' print * shift = .true. endif nymd0 = yymmdd(1) nhms0 = hhmmss(1) endif ! Tick Clock for Next Desired Time Period ! --------------------------------------- nymd = nymd0 nhms = nhms0 nt = num*ndt call tick (nymd,nhms,nt) num = num+1 found = .false. do n=1,ntime if( yymmdd(n).eq.nymd .and. hhmmss(n).eq.nhms ) found = .true. enddo if( .not.found ) then write(6,100) nymd,nhms pa = 0.0 ua = 0.0 va = 0.0 ta = 0.0 qa = 0.0 ha = 0.0 return endif 100 format(/,1x,'Cannot find matching SysError Time for ',i8,2x,i6.6,' ZEROs will be used',/) allocate ( q(im,jm) ) pa = 0.0 ua = 0.0 va = 0.0 ta = 0.0 qa = 0.0 ha = 0.0 do n=1,nvars if( match( c_slp,vname(n) ) ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,0,1,q,rc ) if( rc.eq.0 ) then if( shift ) call hshift ( q,im,jm ) do j=1,jm do i=1,im if( defined(q(i,j),undef) .and. . q(i,j).gt.10000.0 ) then q(i,j) = q(i,j)*0.01 ! Convert Pa=>mb endif enddo enddo call bin ( q,im,jm,pa,idim,jdim,undef,0 ) endif endif do L=1,nl loc = -1 do LL=1,lm if( lev(LL).eq.zlev(L) ) loc = LL enddo if( match( c_u,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = 0.0 endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ua(1,1,L),idim,jdim,undef,1 ) endif if( match( c_v,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = 0.0 endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,va(1,1,L),idim,jdim,undef,1 ) endif if( match( c_t,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = 0.0 endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ta(1,1,L),idim,jdim,undef,0 ) endif if( match( c_q,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = 0.0 endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,qa(1,1,L),idim,jdim,undef,0 ) endif if( match( c_h,vname(n) ) ) then if( loc.ne.-1 ) then call gfio_getvar ( id,vname(n),nymd,nhms,im,jm,loc,1,q,rc ) else q = 0.0 endif if( shift ) call hshift ( q,im,jm ) call bin ( q,im,jm,ha(1,1,L),idim,jdim,undef,0 ) endif enddo enddo deallocate ( q ) return end subroutine interp_time ( nymd1,nhms1, nymd2,nhms2, ntimes, num ) integer nymd1, nhms1, nymd2, nhms2, ntimes, num INTEGER YEAR1, MONTH1, DAY1, SEC1 INTEGER YEAR2, MONTH2, DAY2, SEC2 real time1, time2 INTEGER DAYSCY PARAMETER (DAYSCY = 365*4+1) REAL MNDY(12,4) DATA MNDY /0,31,60,91,121,152,182,213,244,274,305,335,366,397,34*0/ NSECF(N) = N/10000*3600 + MOD(N,10000)/100* 60 + MOD(N,100) DO I=15,48 MNDY(I,1) = MNDY(I-12,1) + 365 ENDDO YEAR1 = NYMD1 / 10000 MONTH1 = MOD(NYMD1,10000) / 100 DAY1 = MOD(NYMD1,100) SEC1 = NSECF(NHMS1) YEAR2 = NYMD2 / 10000 MONTH2 = MOD(NYMD2,10000) / 100 DAY2 = MOD(NYMD2,100) SEC2 = NSECF(NHMS2) time1 = DAY1 + MNDY(MONTH1,MOD(YEAR1,4)+1) + float(sec1)/86400. time2 = DAY2 + MNDY(MONTH2,MOD(YEAR2,4)+1) + float(sec2)/86400. num = 24.0*(time2-time1)/float(ntimes-1) RETURN END subroutine writit ( q,im,jm,lm,qundef,undef ) implicit none integer im,jm,lm,L real q(im,jm,lm) real*4 dum(im,jm) real*4 qundef, undef do L=1,lm where( q(:,:,L).ne.qundef ) dum(:,:) = q(:,:,L) elsewhere dum(:,:) = undef endwhere write(51) dum enddo return end subroutine usage() print *, "Usage: " print * print *, " stats_$ARCH.x -fcst fcst_fname(s)" print *, " -ana ana_fname(s)" print *, " -cli climatology(s)" print *, " <-tag tag>" print *, " <-pref PREF>" print *, " <-nfreq HHMMSS>" print *, " <-syserror SYSERR>" print * print *, "where:" print * print *, " -fcst fcst_fname(s): Filename(s) of pressure-level forecast data (HDF)" print *, " -ana ana_fname(s): Filename(s) in pressure-level analysis data (HDF)" print *, " -cli climatology(s): Filename(s) in pressure-level climatology data (HDF)" print *, " Note:" print *, " A single Mean Climatology, or a" print *, " Set of 4 Climatologies (00,06,12,18z) may be used" print * print *, " -tag tag : Optional Tag for Output Names" print *, " -pref PREF : Optional Reference Pressure for STD Output Printing" print *, " (Default: 500-mb)" print *, " -nfreq HHMMSS : Optional Frequency for Forecast Time Periods" print *, " -syserr SYSERR File : Optional Systematic Error File" print *, " to be Subtracted from Forecasts" print * print *, "creates:" print * print *, "Tag_Name.stats.b{YYYYMMDD}.e{YYYYMMDD}.data" print *, "Tag_Name.stats.b{YYYYMMDD}.e{YYYYMMDD}.ctl1" print *, "Tag_Name.stats.b{YYYYMMDD}.e{YYYYMMDD}.ctl2" print * call exit(7) end subroutine usage subroutine tick (nymd,nhms,ndt) C*********************************************************************** C Purpose C Tick the Date (nymd) and Time (nhms) by NDT (seconds) C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** IF(NDT.NE.0) THEN NSEC = NSECF(NHMS) + NDT IF (NSEC.GT.86400) THEN DO WHILE (NSEC.GT.86400) NSEC = NSEC - 86400 NYMD = INCYMD (NYMD,1) ENDDO ENDIF IF (NSEC.EQ.86400) THEN NSEC = 0 NYMD = INCYMD (NYMD,1) ENDIF IF (NSEC.LT.00000) THEN DO WHILE (NSEC.LT.0) NSEC = 86400 + NSEC NYMD = INCYMD (NYMD,-1) ENDDO ENDIF NHMS = NHMSF (NSEC) ENDIF RETURN end subroutine tick function incymd (NYMD,M) C*********************************************************************** C PURPOSE C INCYMD: NYMD CHANGED BY ONE DAY C MODYMD: NYMD CONVERTED TO JULIAN DATE C DESCRIPTION OF PARAMETERS C NYMD CURRENT DATE IN YYMMDD FORMAT C M +/- 1 (DAY ADJUSTMENT) C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** INTEGER NDPM(12) DATA NDPM /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/ LOGICAL LEAP LEAP(NY) = MOD(NY,4).EQ.0 .AND. (MOD(NY,100).NE.0 .OR. MOD(NY,400).EQ.0) C*********************************************************************** C NY = NYMD / 10000 NM = MOD(NYMD,10000) / 100 ND = MOD(NYMD,100) + M IF (ND.EQ.0) THEN NM = NM - 1 IF (NM.EQ.0) THEN NM = 12 NY = NY - 1 ENDIF ND = NDPM(NM) IF (NM.EQ.2 .AND. LEAP(NY)) ND = 29 ENDIF IF (ND.EQ.29 .AND. NM.EQ.2 .AND. LEAP(NY)) GO TO 20 IF (ND.GT.NDPM(NM)) THEN ND = 1 NM = NM + 1 IF (NM.GT.12) THEN NM = 1 NY = NY + 1 ENDIF ENDIF 20 CONTINUE INCYMD = NY*10000 + NM*100 + ND RETURN C*********************************************************************** C E N T R Y M O D Y M D C*********************************************************************** ENTRY MODYMD (NYMD) NY = NYMD / 10000 NM = MOD(NYMD,10000) / 100 ND = MOD(NYMD,100) 40 CONTINUE IF (NM.LE.1) GO TO 60 NM = NM - 1 ND = ND + NDPM(NM) IF (NM.EQ.2 .AND. LEAP(NY)) ND = ND + 1 GO TO 40 60 CONTINUE MODYMD = ND RETURN end function incymd function nsecf (nhms) C*********************************************************************** C Purpose C Converts NHMS format to Total Seconds C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer nhms, nsecf nsecf = nhms/10000*3600 + mod(nhms,10000)/100*60 + mod(nhms,100) return end function nsecf function nhmsf (nsec) C*********************************************************************** C Purpose C Converts Total Seconds to NHMS format C C*********************************************************************** C* GODDARD LABORATORY FOR ATMOSPHERES * C*********************************************************************** implicit none integer nhmsf, nsec nhmsf = nsec/3600*10000 + mod(nsec,3600)/60*100 + mod(nsec,60) return end function nhmsf