!>>>>>>>>read J's from cloud-JX code (includes fractional cloud treatments) ver 7.3 (2/2015)<<<<<<<<<<<<

      program standalone

      implicit none

      integer,  parameter         :: L_=57, J_=4, N_=640, LT_=33

      character*10, dimension(4), parameter :: TITJ =  &
        ['J-O3(1D)  ', &
         'J-NO2     ', &
         'J-NO3     ', &
         'J-HNO3    ']
      character*30, dimension(5), parameter :: TITJ8 =  &
        ['fort8_G6_c33_v73c.out         ', &
         'fort8_G6_c33_v73.out          ', &
         'fort8_G6_c99_v73.out          ', &
         'fort8_G3_c00_v73.out          ', &
         'fort8_G0_c00_v73.out          ']

      character*6, dimension(5)            :: TITG
      character*16                         :: FILEJ8
      real*8, dimension(L_,J_,N_,5)  :: VJ, VJLN,VRMS
      real*8, dimension(L_,J_,5)     :: WRMS,WJMX,WJMN, VJMN
      real*8, dimension(L_+1)    :: ETAA,ETAB,P,Z, DP,ZK
      real*8                     :: CC, DPSUM
      integer                    :: I,J,L,N
      integer                    :: ICLD,LNRG,ICA,NICA
      real*8, dimension(J_,5)    :: TRMS,TRMS1,TRMS2,TJMN,TJMN1,TJMN2



! get T,P from atmosphere:
      open (7,file='CTM_GrdCld.dat',status='old',err=91)
      do I=1,14
        read (7,*)
      enddo
      do L = 1,L_+1
        read (7,'(4x,2f11.7,12x,f11.2)') ETAA(L),ETAB(L),Z(L)
      enddo
      close (7)
      do L = 1,L_+1
       P(L) = ETAA(L) + ETAB(L)*1000.d0
      enddo

! do stats on J's only for trop (LT_=33, top of level = 88 hPa)
      do L = 1,LT_
        DP(L) = P(L) - P(L+1)
        ZK(L) = 0.5d-3*(Z(L)+Z(L+1))
      enddo


      I = 1
        FILEJ8 = TITJ8(I)
        TITG(I) = FILEJ8(7:12)

       open (7,file=TITJ8(I),status='old',err=91)
         NICA = 0

       write(6,*) 'opened ', FILEJ8
         VJMN(:,:,:) = 0.d0
       do N=1,N_
         read(7,'(40x,3i8,f8.3)') ICLD,LNRG,ICA,CC
         NICA = NICA + ICA
        do L=1,L_
         read(7,'(4x,4e10.3)') (VJ(L,J,N,I), J=1,4)
         do J = 1,4
          VJMN(L,J,I) = VJMN(L,J,I) + VJ(L,J,N,I)
          VJLN(L,J,N,I) = log(VJ(L,J,N,I))
         enddo
        enddo
       enddo

       do L = 1,L_
        do J = 1,4
          VJMN(L,J,I) = VJMN(L,J,I)/float(N_)
        enddo
       enddo

       write(6,*) TITG(I),NICA
       close (7)

      do I=2,5
         FILEJ8 = TITJ8(I)
         TITG(I) = FILEJ8(7:12)

       open (7,file=TITJ8(I),status='old',err=91)
         NICA = 0
       do N=1,N_
         read(7,'(40x,3i8,f8.3)') ICLD,LNRG,ICA,CC
         NICA = NICA + ICA
        do L=1,L_
         read(7,'(4x,4e10.3)') (VJ(L,J,N,I), J=1,4)
         do J = 1,4
          VJLN(L,J,N,I) = log(VJ(L,J,N,I))
         enddo
        enddo
       enddo
       write(6,*) TITG(I),NICA
       close (7)

      enddo

! current stats:  rms error in J's across 640 atmospheres
      do I = 2,5
       do L = 1,LT_
        do J = 1,4
           WRMS(L,J,I) = 0.d0
           WJMN(L,J,I) = 0.d0
           WJMX(L,J,I) = 0.d0
          do N = 1,N_
            WRMS(L,J,I) = WRMS(L,J,I) + 1.d4*(VJLN(L,J,N,I)-VJLN(L,J,N,1))**2
            WJMN(L,J,I) = WJMN(L,J,I) + 1.d2*(VJLN(L,J,N,I)-VJLN(L,J,N,1))
            WJMX(L,J,I) = max(WJMX(L,J,I), abs(VJLN(L,J,N,I)-VJLN(L,J,N,1)) )
          enddo
           WRMS(L,J,I) = WRMS(L,J,I)/float(N_)
           WJMN(L,J,I) = WJMN(L,J,I)/float(N_)
        enddo
       enddo
      enddo

      write(6,'(a)') ' Compare LRG:  mean rms error (%)'
      do J = 1,J_
          write(6,'(a10)')   TITJ(J)
          write(6,'(3x,5a8)')   (TITG(I), I=1,5)
       do L = LT_,1,-1
          write(6,'(i3,1p,e9.2,0p,4f8.1)') L, VJMN(L,J,1), (WRMS(L,J,I), I=2,5)
       enddo
      enddo

      write(6,'(a)') ' Compare LRG:  mean error (%)'
      do J = 1,J_
          write(6,'(a10)')   TITJ(J)
          write(6,'(3x,5a8)')   (TITG(I), I=1,5)
       do L = LT_,1,-1
          write(6,'(i3,1p,e9.2,0p,4f8.1)') L, VJMN(L,J,1), (WJMN(L,J,I), I=2,5)
       enddo
      enddo

      write(6,'(a)') ' Compare LRG:  max error (ln)'
      do J = 1,J_
          write(6,'(a10)')   TITJ(J)
          write(6,'(3x,5a8)')   (TITG(I), I=1,5)
       do L = LT_,1,-1
          write(6,'(i3,1p,e9.2,0p,4f8.3)') L, VJMN(L,J,1), (WJMX(L,J,I), I=2,5)
       enddo
      enddo



      do L = LT_,1,-1
       write(6,'(i3,f8.1)') L, ZK(L)
      enddo


      do J=1,4
       do I=1,5
         DPSUM = 0.d0
         TRMS(J,I) = 0.d0
         TJMN(J,I) = 0.d0
        do L=1,LT_
         TRMS(J,I) = TRMS(J,I) + WRMS(L,J,I)*DP(L)
         TJMN(J,I) = TJMN(J,I) + WJMN(L,J,I)*DP(L)
         DPSUM = DPSUM + DP(L)
        enddo
         TRMS(J,I) = TRMS(J,I)/DPSUM
         TJMN(J,I) = TJMN(J,I)/DPSUM

         DPSUM = 0.d0
         TRMS1(J,I) = 0.d0
         TJMN1(J,I) = 0.d0
        do L=1,7
         TRMS1(J,I) = TRMS1(J,I) + WRMS(L,J,I)*DP(L)
         TJMN1(J,I) = TJMN1(J,I) + WJMN(L,J,I)*DP(L)
         DPSUM = DPSUM + DP(L)
        enddo
         TRMS1(J,I) = TRMS1(J,I)/DPSUM
         TJMN1(J,I) = TJMN1(J,I)/DPSUM

         DPSUM = 0.d0
         TRMS2(J,I) = 0.d0
         TJMN2(J,I) = 0.d0
        do L=18,LT_
         TRMS2(J,I) = TRMS2(J,I) + WRMS(L,J,I)*DP(L)
         TJMN2(J,I) = TJMN2(J,I) + WJMN(L,J,I)*DP(L)
         DPSUM = DPSUM + DP(L)
        enddo
         TRMS2(J,I) = TRMS2(J,I)/DPSUM
         TJMN2(J,I) = TJMN2(J,I)/DPSUM

       enddo
      enddo

      write(6,'(a10,5a8)') 'rms 0-16',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TRMS(J,I), I=2,5)
      enddo
      write(6,'(a10,5a8)') 'rms 0-1',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TRMS1(J,I), I=2,5)
      enddo
      write(6,'(a10,5a8)') 'rms 5-16',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TRMS2(J,I), I=2,5)
      enddo

      write(6,'(a10,5a8)') 'jmn 0-16',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TJMN(J,I), I=2,5)
      enddo
      write(6,'(a10,5a8)') 'jmn 0-1',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TJMN1(J,I), I=2,5)
      enddo
      write(6,'(a10,5a8)') 'jmn 5-16',  (TITG(I), I=2,5),TITG(1)
      do J = 1,J_
          write(6,'(a10,4f8.1)')   TITJ(J), (TJMN2(J,I), I=2,5)
      enddo

      stop
   91 write(6,*)'error in opening CTM_GrdCld.dat file ',TITJ8(I)
      stop
      end
