subroutine cpr (nlon, nlat, nlev, nlevp, numcases, lnorm, verbose)!! $Id: cpr.F,v 1.2.6.1 2002/05/13 17:59:24 erik Exp $!  use precision  use header  use stats  use nldat  implicit none  include 'netcdf.inc'!! Input arguments!  integer, intent(in) :: nlon, nlat, nlev, nlevp, numcases  logical lnorm            ! Whether to print lvl-by-lvl L2 & L-INF norm stats  logical verbose!! Local workspace (static or stack)!  character*(nf_max_name) name  integer nvars            ! number of variables on tape  integer natts  integer n,n2             ! variable index  integer ndims(2)  integer idim  integer xtype(2)  integer dimids(nf_max_dims,2)  integer i, j, k          ! spatial indices  integer numlev,numlev2   ! number of levels for each field  integer isub, nval, ii  integer itime, itime2    ! loop index over time  integer ret              ! return code  integer, dimension(3) :: start2d, count2d    real(r8) wsum  real(r8) pdelbar  real(r8) denom  real(r8) rd  real(r8), parameter :: timeepsilon = 1.e-9  real(r8) :: timediff  logical twocases         ! True => 2 h-tapes are being analyzed  logical found!! Local workspace (dynamic)!  real(r8), dimension(nlon,nlat,nlevp) :: arr, arr2  real(r8) diff(nlon)             ! difference field  real(r8) pdel(nlon,nlev,numcases)  real(r8) pmid(nlon,nlev)  real(r8) rdiff(nlon)            ! relative difference field  real(r8) ps(nlon,nlat,numcases)  real(r8) wbar(nlat)  real(r8) totmass  integer indx(nlon)          ! indices of non-zero diffs!! Externals!  integer ismax, fillarr, lenchr  external ismax, fillarr, lenchr!! Define start, count arrays!  start2d(:) = 1  count2d(1) = nlon  count2d(2) = nlat  count2d(3) = 1  twocases = numcases.eq.2!! Set surface pressure to constant once and for all if requested!  if (psid(1).lt.0 .or. psid(2).lt.0) then    if (twocases .and. (p0(1) /= p0(2))) then      write(6,*)'CPR: p0 values do not match'      stop 99    end if    ps(:,:,:) = p0(1)  end if!! Find number of variables on the 1st tape!  call wrap_inq_nvars (ncid(1), nvars)!! Big loop over number of time slices.  Not written as explicit "do" since! may need to adjust time index in the middle of the loop to match time! levels on tapes.!  itime = 1  itime2 = 110 continue!! Check equality of time steps!  if (matchts .and. twocases) then    timediff = abs (time(itime) - time2(itime2))    do while (timediff .gt. timeepsilon)      write(6,'(1x,a,f12.4,/1x,a,f12.4)')&           'Current times differ:time1=', time(itime), &           '                     time2=', time2(itime2)      if (time(itime).lt.time2(itime2)) then        itime = itime + 1        if (itime.gt.ntime(1)) then          write(6,*)'End of tape 1 reached: stopping'          stop 0        end if      else        itime2 = itime2 + 1        if (itime2.gt.ntime(2)) then          write(6,*)'End of tape 2 reached: stopping'          stop 0        end if      end if      timediff = abs (time(itime) - time2(itime2))    end do  end if  if (twocases) then    write(6,800)trim(case(1)), trim(case(2)),  &                trim(title(1)), trim(title(2)), &#ifndef SUN                date_written(itime,1), date_written(itime2,2), &#else                " ", " ", &#endif                time_written(itime,1), time_written(itime2,2), &                itime, itime2, &                ntime(1), ntime(2), &                nsteph(itime,1), nsteph(itime2,2), &                ncdate(itime,1), ncdate(itime2,2), &                ncsec(itime,1),  ncsec(itime2,2)  else                               ! One tape only.    write(6,801)case(1)(1:lenchr(case(1))),  &                title(1)(1:lenchr(title(1))), &                date_written(itime,1), &                time_written(itime,1), &                itime, &                ntime(1), &                nsteph(itime,1), &                ncdate(itime,1), &                ncsec(itime,1)  end if                    ! twocases!! Read in and check the data fields sequentially!  if (iprs.ne.0) then    write(6,'(/,a,i4,a,i4, a,i3,a,i3, a,i3,a,i3)') &     ' SUMMARY OF FIELD INFORMATION FOR LONGITUDE INDICES ', &     iprs,' through ',ipre,', LATITUDE INDICES ', &     latprs,' through ',latpre,', LEVELS ',kprs,' through ',kpre    write(6,'(a,t10,a,t14,a,t18,a,t23,a,t35,a,t47,a,t70,a)') &     ' FIELD','LON','LAT','LEV','DIFF','RDIFF','CASE1','CASE2'  end if  wsum = 0.  totmass = 0.  if (twocases) then    wbar(:) = 0.5*(gw(:,1) + gw(:,2))    if (psid(1).gt.0 .and. psid(2).gt.0) then      start2d(3) = itime      call wrap_get_vara_realx (ncid(1), psid(1), start2d, count2d, ps(1,1,1))      start2d(3) = itime2      call wrap_get_vara_realx (ncid(2), psid(2), start2d, count2d, ps(1,1,2))    end if    do j=1,nlat      do n=1,2        call plevs0 (nlon, nlev, hyai(1,n), hybi(1,n), hyam(1,n), &                     hybm(1,n), ps(1,j,n), pmid, pdel(1,1,n), p0(n))      end do!! Compute total mass for 2 methods!      do k=1,nlev        do i=1,nlon          pdelbar = 0.5*(pdel(i,k,1) + pdel(i,k,2))          totmass = totmass + pdelbar*wbar(j)        end do      end do      wsum = wsum + nlon*wbar(j)    end do  else    wbar(:) = gw(:,1)  end if!! Loop over all variables on the 1st tape!  do n=1,nvars    dimids(:,1) = -1    call wrap_inq_var (ncid(1), n, name, xtype(1), ndims(1), dimids(1,1), natts)    if (verbose) then      write(6,*) 'processing ', name(1:lenchr(name))    end if!! Determine number of levels for this variable!    found = .false.    numlev = 1    do idim=1,ndims(1)      if (dimids(idim,1).eq.levdimid(1)) then        numlev = nlev        if (found) then          write(6,*) name(1:lenchr(name)),' has 2 level dims: skipping'          cycle        end if        found = .true.      else if (dimids(idim,1).eq.ilevdimid(1)) then        numlev = nlevp        if (found) then          write(6,*) name(1:lenchr(name)),' has 2 level dims: skipping'          cycle        end if        found = .true.      end if    end do    ret = fillarr (ncid(1), n, dimids(1,1), londimid(1), latdimid(1), &                   levdimid(1), ilevdimid(1), nlon, nlat, &                   numlev, arr, itime)!! Skip to the next field if fillarr failed!    if (ret.lt.0) cycle!! Initialize statistics, then modify for this field!    call initstats    xmx(1) = -1.d99    xmn(1) = +1.d99    do k=1,numlev      do j=1,nlat        do i=1,nlon          xbar(1) = xbar(1) + abs(arr(i,j,k))          if (arr(i,j,k).gt.xmx(1)) then            xmx(1) = arr(i,j,k)            imx(1) = i            jmx(1) = j            kmx(1) = k          end if          if (arr(i,j,k).lt.xmn(1)) then            xmn(1) = arr(i,j,k)            imn(1) = i            jmn(1) = j            kmn(1) = k          end if        end do      end do    end do    if (twocases) then      ret = nf_inq_varid (ncid(2), name, n2)      if (ret.ne.NF_NOERR) then        if (verbose) then          write(6,*) name, ' not found on tape 2'        end if        cycle      end if              dimids(:,2) = -1      call wrap_inq_var (ncid(2), n2, name, xtype(2), ndims(2), &                         dimids(1,2), natts)      if (xtype(1).ne.xtype(2)) then        write(6,*)'NOTE: type variables do not match for field ', &                  name(1:lenchr(name))      end if      if (xtype(1).eq.NF_FLOAT .or. xtype(2).eq.NF_FLOAT) then        write(6,*)'NOTE: At least 1 tape has 32-bit values for field ', &                  name(1:lenchr(name))      end if!! Determine number of levels for this variable!      found = .false.      numlev2 = 1      do idim=1,ndims(2)        if (dimids(idim,2).eq.levdimid(2)) then          numlev2 = nlev          if (found) then            write(6,*) name(1:lenchr(name)),' has 2 level dims: skipping'            cycle          end if          found = .true.        else if (dimids(idim,2).eq.ilevdimid(2)) then          numlev2 = nlevp          if (found) then            write(6,*) name(1:lenchr(name)),' has 2 level dims: skipping'            cycle          end if          found = .true.        end if      end do      if (numlev2.ne.numlev) then        write(6,*) name(1:lenchr(name)), ' numlev mismatch: skipping'        cycle      end if      ret = fillarr (ncid(2), n2, dimids(1,2), londimid(2), latdimid(2), &                     levdimid(2), ilevdimid(2), nlon, nlat, &                     numlev, arr2, itime2)!! Skip to the next field if fillarr failed!      if (ret.lt.0) then        write(6,*)'Tape 2 does not have field ',name,': skipping'        cycle      end if      xmx(2) = -1.d99      xmn(2) = +1.d99      do k=1,numlev        do j=1,nlat          npos = npos + nlon          diff(:) = abs (arr(:,j,k) - arr2(:,j,k))!! Only compute mass weighted RMS for variables dimensioned nlev since don't ! know pdel for nlevp!          if (numlev.eq.1) then            do i=1,nlon              diffmw = diffmw + (diff(i)**2)*wbar(j)            end do          else if (numlev.eq.nlev) then            do i=1,nlon              pdelbar = 0.5*(pdel(i,k,1) + pdel(i,k,2))              diffmw = diffmw + (diff(i)**2)*pdelbar*wbar(j)            end do          end if          do i=1,nlon            rms = rms + diff(i)**2            xbar(2) = xbar(2) + abs(arr2(i,j,k))            if (arr2(i,j,k).gt.xmx(2)) then              xmx(2) = arr2(i,j,k)              imx(2) = i              jmx(2) = j              kmx(2) = k            end if                        if (arr2(i,j,k).lt.xmn(2)) then              xmn(2) = arr2(i,j,k)              imn(2) = i              jmn(2) = j              kmn(2) = k            end if          end do          if (iprs.ne.0) then            if (numlev.eq.1 .or. (k.ge.kprs .and. k.le.kpre)) then              if (j.ge.latprs .and. j.le.latpre) then                do i=iprs,ipre                  denom = 0.5*(abs(arr(i,j,k)) + abs(arr2(i,j,k)))                  rd = 0.                  if (denom.ne.0.) rd = diff(i)/denom                  write(6,910) name,i,j,k,diff(i),rd,arr(i,j,k),arr2(i,j,k)                end do              end if            end if          end if!! Test on half of difference field rather than full field due to 0.5 factor! used in computation of "denom" later.!           nval = 0          do i=1,nlon            if (diff(i).ne.0.) then              nval = nval + 1              indx(nval) = i            end if          end do          if (nval.gt.0) then                   ! at least 1 difference            ndif = ndif + nval            isub = ismax(nlon, diff, 1)            if (diff(isub).gt.difmx) then    ! Save max diff info              difmx = diff(isub)              dmxsv(1) = arr(isub,j,k)          ! Save values              dmxsv(2) = arr2(isub,j,k)              idmxsv = isub              jdmxsv = j              kdmxsv = k            end if            rdiff(:) = 0.            do ii=1,nval                        ! Compute relative diffs              i = indx(ii)              denom = max(abs(arr(i,j,k)), abs(arr2(i,j,k)))              rdiff(i) = diff(i)/(2.*denom)              rdbar = rdbar + rdiff(i)              rdlnbar = rdlnbar - log10(rdiff(i))            end do            isub = ismax(nlon,rdiff,1)                        if (rdiff(isub).gt.rdifmx) then  ! Save max relative diff info              rdifmx = rdiff(isub)              rdmxsv(1) = arr(isub,j,k)         ! Save values & indices              rdmxsv(2) = arr2(isub,j,k)              irdmxsv = isub              jrdmxsv = j              krdmxsv = k            end if          end if        end do      end do    end if!! Print stats!    call printstats (twocases, name, numlev, nlev, nlevp, &                     nlon, nlat, wsum, totmass)  end do  ! n=1,nvars!! Proceed to next time slice if more are available!  itime = itime + 1  if (itime.gt.ntime(1)) then    write(6,*)'End of tape 1 reached: stopping'    stop 0  end if  if (twocases) then    itime2 = itime2 + 1    if (itime2.gt.ntime(2)) then      write(6,*)'End of tape 2 reached: stopping'      stop 0    end if  end if  goto 10800 format(' CASE:',2(A,1X),/, &           ' CASE1 TITLE:',a80,/, ' CASE2 TITLE:',a80,/, &           ' date_written: ',2(A8,1X),/, &           ' time_written: ',2(A8,1X),/, &           ' MFILH:  ',2I9,/,' MFILTH: ',2I9,/,' NSTEPH: ',2I9,/, &           ' NCDATE: ',2I9,/, &           ' NCSEC:  ',2I9,/,' MHISF:  ',2I9,/,' MFSTRT: ',2I9,/)801 format(' CASE1:',a,/, &           ' CASE1 TITLE:',a80,/, &           ' date_written: ',1(A8,1X),/, &           ' time_written: ',1(A8,1X),/, &           ' MFILH:  ',1I9,/,' MFILTH: ',1I9,/,' NSTEPH: ',1I9,/, &           ' NCDATE: ',1I9,/, &           ' NCSEC:  ',1I9,/,' MHISF:  ',1I9,/,' MFSTRT: ',1I9,/)910 format(1x,a8,'(',i3,',',i2,',',i2,')=',1p,2e12.4,2e23.15)end subroutine cprinteger function fillarr (ncid, varid, dimids, londimid, latdimid, &                          levdimid, ilevdimid, nlon, nlat, &                          numlev, arr, itime)  use precision  implicit none  include 'netcdf.inc'  integer :: ncid, varid  integer :: dimids(NF_MAX_DIMS)  integer :: londimid, latdimid, levdimid, ilevdimid  integer :: nlon, nlat, numlev  integer :: itime  real(r8) :: arr(nlon,nlat,numlev)!! Local workspace!  integer j, k  integer start(4), count(4)  real(r8), allocatable :: arrxzy(:,:,:)  fillarr = 0       ! Initialize return value to success  if (dimids(1).eq.londimid .and. &      dimids(2).eq.latdimid .and. &      (dimids(3).eq.levdimid .or. dimids(3).eq.ilevdimid)) then    start(1) = 1    start(2) = 1    start(3) = 1    start(4) = itime    count(1) = nlon    count(2) = nlat    count(3) = numlev    count(4) = 1    call wrap_get_vara_realx (ncid, varid, start, count, arr)  else if (dimids(1).eq.londimid .and. &           (dimids(2).eq.levdimid .or. dimids(2).eq.ilevdimid) .and. &           dimids(3).eq.latdimid ) then    start(1) = 1    start(2) = 1    start(3) = 1    start(4) = itime    count(1) = nlon    count(2) = numlev    count(3) = nlat    count(4) = 1    allocate(arrxzy(nlon,numlev,nlat))    call wrap_get_vara_realx (ncid, varid, start, count, arrxzy)    do k=1,numlev      do j=1,nlat        arr(:,j,k) = arrxzy(:,k,j)      end do    end do    deallocate(arrxzy)  else if (dimids(1).eq.londimid .and. dimids(2).eq.latdimid) then    start(1) = 1    start(2) = 1    start(3) = itime    start(4) = 0    count(1) = nlon    count(2) = nlat    count(3) = 1    count(4) = 0    call wrap_get_vara_realx (ncid, varid, start, count, arr)  else    fillarr = -1  end if  returnend function fillarr#if ( ! defined TIMING )subroutine t_startf (xxx)  character*(*) xxx  returnend subroutine t_startfsubroutine t_stopf (xxx)  character*(*) xxx  returnend subroutine t_stopfsubroutine t_prf (xxx)  integer xxx  returnend subroutine t_prf#endif