
!>>>>>>>>cloud-JX code (includes fractional cloud treatments) ver 7.2 (12/2013)<<<<<<<<<<<<
      program atm_columns
! USES:
      USE RRSW_FASTJ_CMN, ONLY : W_rrtmg
      USE FJX_CMN_MOD
      USE FJX_INIT_MOD
      USE CLD_SUB_MOD, ONLY : CLOUD_JX
      USE FJX_SUB_MOD

! CALL sequence
!      call rrtmg_sw_inp(0) 
!      call INIT_FJX (TITLJXX,JVN_,NJXX)
!      call read_1DCOL (read met data)
!      call input_for_cloudJ (generate input needed for cloudj; call ACLIM_FJX in here)
!      call SOLAR_JX(PHOTAU,IDAY,YGRD,XGRD, SZA,U0,SOLF)
!      call rrtmg_sw_inp(1) !obtain optical Depth for gas absortion extending to 3.84 microns
!      call CLOUD_JX (U0,SZA,REFLB,SOLF,FG0,LPRTJ,PPP,ZZZ,TTT, &
!             DDD,RRR,OOO,   LWP,IWP,REFFL,REFFI,     CLF,CWC,   &
!             AERSP,NDXAER,L1_,AN_,VALJXX,JVN_,   &
!             CLDFLAG,NRANDO,IRAN,LNRG,NICA,JCOUNT)
      IMPLICIT NONE
      INTEGER  IUNIT
      CHARACTER*27 FIL
      CHARACTER*80 PATH
      real*8  Time1, Time2
      real*8, dimension(L1_) :: ETAA, ETAB, ZOFL,CLDP,AER1,AER2,QW
      real*8, dimension(L1_+1) :: TI, RI ! these are edge values and put one more at top of atmosphere as in ppp
      real*8, dimension(L_) :: WLC,WIC
      real*8  GMTAU,PHOTAU,ALBEDO, XLNG,YLAT,XGRD,YGRD,PSURF, SCALEH
      real*8  CF,PMID,PDEL,ZDEL,ICWC,F1,ZKM
      integer, dimension(L1_):: NAA1,NAA2
      integer MYEAR, MONTH, IDAY, JLAT, ILON, NHRMET
      integer I,J,K,L,N,JP, NN
      integer NRAN, RANSEED, LTOP, NJXX
      integer CLDFLDS 
      integer ICLDP1, ICLDP2, ISKIP, ICLDFLAG, JCLDFLAG
      real*8, dimension(LWEPAR) :: CLDFRW,CLDIWCW,CLDLWCW   ! WCW=Cloud Water Content (g/g)
      character*6, dimension(JVN_)   ::  TITLJXX
      real*8, dimension(L_,8,4)  :: ZPJCLD, ZPJAVG
      real*8, dimension(21,4)    :: ERRJ, ERRJJ, ERRJ2
      integer JP04,JP09,JP11,JP15, ICLD
      integer, dimension(8)      :: JCNT
      character*11, dimension(4) ::  TITJX      
!     setup cldfiles for rrtmg
      integer inflag, iceflag, liqflag
      integer KKK
      Logical LFINISH ! for call cld_errstats subroutine    

!--------------------key params sent to CLOUD_JX-------------------------
      real*8                     :: U0,SZA,REFLB,SOLF, CLDCOR
      real*8                     :: FG0 
      logical                    :: LPRTJ
      real*8,  dimension(L1_+1)  :: PPP,ZZZ
      real*8,  dimension(L1_  )  :: TTT,DDD,RRR,OOO
      real*8,  dimension(L1_)    :: LWP,IWP,REFFL,REFFI
      real*8,  dimension(L1_)    :: CLF,CWC
      integer, dimension(L1_)    :: CLDIW
      real*8,  dimension(L1_,AN_):: AERSP
      integer, dimension(L1_,AN_):: NDXAER
      real*8,  dimension(L_,JVN_) ::VALJXX
      real*8,  dimension(S_+2, L1_)::SKPERD, ODW
      integer                    :: CLDFLAG,NRANDO,IRAN,LNRG
      integer                    :: NICA,JCOUNT
!---U0 = cos (SZA), SZA = solar zenith angle
!---REFLB = Lambertian reflectivity at the Lower Boundary
!---SOLF = solar flux factor for sun-earth distance
!---FG0 = scale for asymmetry factor to get equivalent isotropic (CLDFLAG=3 only)
!---LPRTJ = .true. = turn on internal print in both CLOUD_JX & PHOTO_JX
!--- P = edge press (hPa), Z = edge alt (m), T = layer temp (K)
!--- D = layer dens-path (# molec /cm2), O = layer O3 path (# O3 /cm2)
!--- R = layer rel.hum.(fraction)
!---LWP/IWP = Liquid/Ice water path (g/m2)
!---REFFL/REFFI = R-effective(microns) in liquid/ice cloud
!---CLF = cloud fraction (0.0 to 1.0)
!---CWC = cloud water content (g/m3) both liq & ice
!---CLDIW = integer denoting cloud in layer: 1=water, 2=ice, 3=both
!---AERSP = aerosol path (g/m2) & NDXAER = aerosol index type
!---  aerosols are dimensioned with up to AN_ different types in an ICA layer
!---L1_ = parameter, dim of profile variables, L_+1 for top (non CTM) layer
!---AN_ = parameter, dim of number of aerosols being passed
!---VALJXX = J-values from CLOUD_JX & PHOTO_JX
!---JVN_ = dim of max number of J-s reported out (in the order of fast-JX, not CTM)
!---CLDFLAG = integer index for type of cloud overlap
!---CLOUD_JX:   different cloud schemes
!---CLOUD_JX:   different cloud schemes (4:8 require max-ran overlap algorithm)
!       CLDFLAG = 1  :  Clear sky J's
!       CLDFLAG = 2  :  Averaged cloud cover
!       CLDFLAG = 3  :  cloud-fract**3/2, then average cloud cover
!       CLDFLAG = 4  :  Average direct solar beam over all ICAs, invert to get clouds
!       CLDFLAG = 5  :  Random select NRANDO ICA's from all(Independent Column Atmos.)
!       CLDFLAG = 6  :  Use all (up to 4) quadrature cloud cover QCAs (mid-pts of bin)
!       CLDFLAG = 7  :  Use all (up to 4) QCAs (average clouds within each Q-bin)
!       CLDFLAG = 8  :  Calculate J's for ALL ICAs (up to 20,000 per cell!)
!---NRANDO = number of random ICAs to do J's for (CLDFLAG=4)
!---IRAN = starting index for random number selection
!---LNRG = flag for setting max-ran overlap groups:
!---     =0   break max overlap groups at cloud fraction = 0
!---     =3else fixed 3 layer (1:9, 9:last LWcloud, LWclud+1:LTOP)
!---     else(=6) fixed correlated length max-overlap layers
!---NICA = total number of ICAs
!---fast-JX:  INIT_JX is called only once to read in & store all fast-JX data:
!              also sets up random sequence for cloud-JX
!-----------------------------------------------------------------------
! initiation (call only once)
      CALL CPU_TIME(Time1)
      IF(W_RRTMG .NE. 0) call RRTMG_SW_INP(0)
      CALL INIT_FJX (TITLJXX,JVN_,NJXX)
!-----------------------------------------------------------------------
      LPRTJ = .true.
!
!--Set up atmosphere for a single column and time for J-values calculation
!--Nominally taken from CTM, but for standalone here is read in
!read in met data
      PATH="Data/"
      FIL= "CTM_GrdCld.dat"
      iunit=77
! READ column atmosphere meteorological variables 
      CALL READ_1DCOL(PATH, FIL, iunit, MYEAR, IDAY, MONTH, GMTAU, PHOTAU, ALBEDO, ILON, JLAT, XLNG,YLAT,&
           TI ,RI, QW, ETAA, ETAB, PSURF, NAA1, AER1, NAA2, AER2, & 
           ICLDP1, ICLDP2, ICLDFLAG, JCLDFLAG, LNRG, NRANDO, CLDCOR, FG0)
!convert to radian 
      YGRD = YLAT*CPI180
      XGRD = XLNG*CPI180
      REFLB = ALBEDO
      CLDFLDS= ICLDP2- ICLDP1 +1
! Compute PPP, ZZZ, TTT, RRR, DDD, OOO, NDXAER, AERSP 
      CALL INPUT_FOR_CLOUDJ(YLAT, MONTH, TI, RI, QW, ETAA, ETAB, PSURF, NAA1, AER1, NAA2, AER2, &
           PPP, ZZZ, TTT, RRR, DDD, OOO, NDXAER, AERSP)
! Call Solar_JX to output SZA, U0, SOLF needed for write_to_rrtmg
      CALL SOLAR_JX(PHOTAU,IDAY,YGRD,XGRD, SZA,U0,SOLF)
!        write to an output that will be the meteorological standard input for RRTMG column runs
!      CALL WRITE_INPUTMET_RRTMG(L1_, IDAY, ALBEDO, SZA)
!**important note that WRITE_INPUTMET_RRTMG is to be called before RRTMG_SW_INP(1) because colh2o etc is  scaled by 1.e20 in RRTMG_SW_INP(1) 
! Compute optical depths of trace gases for Fast-JX 
      IF ( W_RRTMG .NE. 0)CALL RRTMG_SW_INP(1)
!---following is readin for cloud data, currently has 160 atmospheres
!   from tropical T319 ECMWF atmosphere used in UCI CTM.
      ZPJAVG(:,:,:) = 0.d0
      ERRJJ(:,:)    = 0.d0
      ERRJ2(:,:)    = 0.d0
      JCNT(:)       = 0
!      LFINISH =.FALSE. ! FOR CALL CLD_ERRSTATS 
! read(77,*) continues in  READ_1DCOL after met data being read
      DO ICLD= ICLDP1, ICLDP2 !# of cloud profiles
         IF (ICLDP1 .GT. 1)THEN
            DO ISKIP =1, ICLDP1-1
!               print*,'ISKIP=', ISKIP
               READ (77,*)
               DO L = LWEPAR,1,-1
                  READ (77,'(I3, 1P, E14.5, 28X, 2E14.5)') &
                       J,CLDFRW(L),CLDLWCW(L),CLDIWCW(L)
               ENDDO
            ENDDO
         ENDIF
!         print*,'cloud reading'
         READ (77,*)
         DO L = LWEPAR,1,-1
            READ (77,'(I3,1P,E14.5, 28X, 2E14.5)') &
                 J,CLDFRW(L),CLDLWCW(L),CLDIWCW(L)
!            write (6,'(I3,1P,E14.5, 28X, 2E14.5)') &
!                 J,CLDFRW(L),CLDLWCW(L),CLDIWCW(L)
         ENDDO
!---convert cloud data from our EC met fields into cloud data for cloud-JX
!---init data = cloud fraction, and water content (g/m3) averaged over cell
!---     needs: cloud fraction, ice- and liq-water path (in cloud)
!---     and R-effective of ice and liquid clouds
         LTOP  = LWEPAR
         if (MAXVAL(CLDFRW) .le. 0.005d0) then
            IWP(:) = 0.d0
            REFFI(:) = 0.d0
            LWP(:) = 0.d0
            REFFL(:) = 0.d0
         ENDIF
         DO L = 1,LTOP
            CLDIW(L) = 0
            CF  = CLDFRW(L)
            IF (CF .GT. 0.005D0) THEN
               CLF(L) = CF
               WLC(L) = CLDLWCW(L) / CF  !in-cloud WLC
               WIC(L) = CLDIWCW(L) / CF  !in cloud WIC
               !  CLDIW is an INT flag: 1 = water cld, 2 = ice cloud, 3 = both
               if (WLC(L) .gt. 1.d-11) CLDIW(L) = 1
               if (WIC(L) .gt. 1.d-11) CLDIW(L) = CLDIW(L) + 2
            ELSE
               CLF(L) = 0.d0
               WLC(L) = 0.d0
               WIC(L) = 0.d0
            ENDIF
         ENDDO
!---derive R-effective for clouds:  the current UCI algorithm - use your own
         DO L = 1,LTOP
!---ice clouds
            IF (WIC(L) .GT. 1.D-12) THEN
               PDEL =    PPP(L) - PPP(L+1)
               ZDEL =   (ZZZ(L+1) - ZZZ(L))*0.01D0  ! M 
               IWP(L) =  1000.D0*WIC(L)*PDEL*G100  ! G/M2
               ICWC =   IWP(L) / ZDEL        ! g/m3
               REFFI(L) = 164.d0 * (ICWC**0.23d0)
!  temporary setups  
!               if(TTT(L) .lt. 234)then
!                  REFFI(L)=60
!               else
!                  REFFI(L)=40
!               endif
            else
               IWP(L) = 0.d0
               REFFI(L) = 0.d0
            endif
!---water clouds
            if (WLC(L) .gt. 1.d-12) then
               PMID = 0.5d0*(PPP(L)+PPP(L+1))
               PDEL = PPP(L) - PPP(L+1)
               F1   = 0.005d0 * (PMID - 610.d0)
               F1   = min(1.d0, max(0.d0, F1))
               LWP(L) = 1000.d0*WLC(L)*PDEL*G100     ! g/m2
               REFFL(L) = 9.6d0*F1 + 12.68d0*(1.d0-F1)
!              temporary fix
!               REFFL(L)= 10
!               print*,'LWP=', LWP(L), 'REFFL(L)=', REFFL(L)
            else
               LWP(L) = 0.d0
               REFFL(L) = 0.d0
            endif
         enddo
!write out standard cloud input for rrtmg's column model
!          CALL WRITE_INPUTCLD_RRTMG(ICLD, LTOP, CLDFRW, IWP, LWP, REFFI, REFFL)  !        
!!         
!---cloud input as interpreted by fast_JX
!---extinction K(m2/g) = 3/4 * Q / [Reff(micron) * density(g/cm3)]
!           ODL = LWP(L) * 0.75d0 * 2.1d0 / REFFL(L)
!           ODI = IWP(L) * 0.75d0 * 2.0d0 / (REFFI(L) * 0.917d0)

!---fast-JX:  SOLAR_JX is called  once per grid-cell to set U0, SZA, SOLF
!--- your CTM code may have its own way of calculating and passing these quantities
!-----------------------------------------------------------------------
!-----------------------------------------------------------------------
!        if (LPRTJ) then
!          write(6,'(a,f8.3,3f8.5)')'solar zenith angle, solar-f' &
!                 ,SZA,SOLF,U0,REFLB
!            write(6,'(a,f8.3,f8.3)') 'lat/lng',YLAT,XLNG
!        call JP_ATM0(PPP,TTT,DDD,OOO,ZZZ, L_)
!         endif
!  locate the position of random number sequence based on year/day/hour
         IRAN = 13+ILON+3*JLAT+5*(MYEAR-1900)+7*IDAY + 11*nint(GMTAU)
         ZPJCLD(:,:,:) = 0.d0
         DO KKK=1,1! redundant loop for cpu testing
         DO CLDFLAG = ICLDFLAG, JCLDFLAG ! 1 to 8
            WRITE(RecTle,'(a, I3, a, I3)') &
                 'CLD PROF #', ICLD, 'CLD FLAG=', CLDFLAG 
!            WRITE(6,'(a, 3I5)') '%call CLOUD_JX w/flag=', CLDFLAG,LNRG,NRANDO
!=======================================================================
!            do L=L_ +1, 1, -1
!              write(6, '(F12.3, F12.3, F10.2, 3E14.5)'), PPP(L), ZZZ(L)/1.d5, TTT(L), DDD(L), RRR(L), OOO(L)
!            enddo
!
            CALL CLOUD_JX (U0,SZA,REFLB,SOLF,FG0,LPRTJ,PPP,ZZZ,TTT, &
                 DDD,RRR,OOO, LWP,IWP,REFFL,REFFI, CLF,CLDCOR,CLDIW, &
                 AERSP,NDXAER,L1_,AN_,VALJXX, SKPERD, ODW, JVN_, &
                 CLDFLAG,NRANDO,IRAN,LNRG,NICA,JCOUNT)
            print*,'end of call cloud_jx'
            JCNT(CLDFLAG) = JCNT(CLDFLAG) + JCOUNT
!
!            Below call output_heatingRate Major diagnostics output 
            IF(LPRTJ)THEN
                CALL OUTPUT_HeatingRate(L1_, ICLD,CLDFLAG, PPP, ZZZ, VALJXX, SKPERD) 
!               CALL OUTPUT_OD(ICLD,CLDFLAG,PPP, ODW)
            ENDIF

!            CALL CLD_ERRSTATS(VALJXX,PPP, CLDFLAG, JP04, JP09, JP11, JP15, ZPJAVG, ERRJ2, ERRJJ, JCNT, CLDFLDS, LFINISH, ZZZ)

!=======================================================================
!            if(CLDFLAG .eq. 8) then
!               write(6,'(a,3i8,f8.3)') ' cloud-JX v7.3 ALL ICAs: ICLD/LNRG/NICAs',&
!                    ICLD,LNRG, NICA,CLDCOR
!               write(6,'(a,3i8,f8.3)') ' cloud-JX v7.3 ALL ICAs: ICLD/LNRG/NICAs',&
!                    &ICLD,LNRG,NICA,CLDCOR
!               do L=1,L_
!                  write(6,'(i4,1p,4e10.3)')&
!                       L,VALJXX(L,JP04),VALJXX(L,JP09),VALJXX(L,JP11),VALJXX(L,JP15)
!               enddo
!            endif
         ENDDO ! CLDFLAG
         ENDDO ! KKK redundant loop for testing cpu time
      ENDDO ! ICLD
!
!      LFINISH =.TRUE.
!      CALL CLD_ERRSTATS (VALJXX,PPP, CLDFLAG, JP04, JP09, JP11, JP15, ZPJAVG, ERRJ2, ERRJJ, JCNT,CLDFLDS, LFINISH, ZZZ)
!
!
      call cpu_time(Time2)
      write(6,*)'Elapsed CPU time = ', Time2 - Time1
      stop
 
      end

      SUBROUTINE READ_1DCOL(PATH, FIL, IU, MYEAR, IDAY, MON, GMTAU, PHOTAU, ALBEDO, ILON, JLAT, XLNG, YLAT,&
           TI, RI, QW, ETAA, ETAB, PSURF, NAA1, AER1, NAA2, AER2, & 
           ICLDP1, ICLDP2, ICLDFLAG, JCLDFLAG, LNRG, NRANDO, CLDCOR, FG0)
         USE FJX_CMN_MOD, ONLY: L1_
         IMPLICIT NONE
         CHARACTER*80, INTENT(IN)::PATH
         CHARACTER*27, INTENT(IN):: FIL
         INTEGER,   INTENT(IN):: IU
         INTEGER,   INTENT(OUT):: IDAY, MON, ICLDP1, ICLDP2, ICLDFLAG, JCLDFLAG, LNRG, NRANDO
         REAL*8,    INTENT(OUT):: GMTAU, PHOTAU, ALBEDO, XLNG, YLAT, PSURF, FG0, CLDCOR
         REAL*8, DIMENSION(L1_), INTENT(OUT) :: ETAA,ETAB,AER1,AER2
         INTEGER,DIMENSION(L1_), INTENT(OUT) :: NAA1,NAA2
! these are edge values and put one more at top of atmosphere as in ppp 
         REAL*8, DIMENSION(L1_+1), INTENT(OUT) :: TI, RI !edge values
         REAL*8, DIMENSION(L1_), INTENT(OUT) :: QW ! layer values
 
         INTEGER MYEAR, NHRMET, JLAT, ILON, J, L 
         REAL*8, DIMENSION(L1_) ::  ZOFL

         open (iu, file =adjustl(trim(PATH))//trim(FIL), &
              form='formatted', err=91)
         read (iu,*)
         read (iu,'(i6)') MYEAR
         read (iu,'(i6)') IDAY
         read (iu,'(i6)') MON
         read (iu,'(i6)') NHRMET
! I have reset PHOTAU to be SZA in fjx_sub_mod.f90.  Should change it back PHOTAU normally should be GMT (Universal Time)
         read (iu,'(2f6.2)') GMTAU, PHOTAU
         read (iu,'(f6.1,i4)') YLAT, JLAT
         read (iu,'(f6.1,i4)') XLNG, ILON
         read (iu,'(f6.1)') PSURF
!         write(6, '(f6.1)') psurf
         read (iu,'(f6.2)') ALBEDO
         read (iu,'(2i5,f5.1)')LNRG, NRANDO, CLDCOR      
         read (iu,'(f5.2)') FG0
         read (iu,'(2i5)') ICLDP1, ICLDP2    
         read (iu,'(2i5)') ICLDFLAG, JCLDFLAG
         read (iu,*)
         print*,'psurf=',psurf
         write(6,'(A)')'L    ETAA        ETAB     T     , Specific Q,  Height      Areosols'
         do L = 1, L1_
            read (iu,'(i3,1x,2f11.7,2x,f5.1,f5.2,1x, e8.2, f10.2,2(f7.3,i4))') &
                 J, ETAA(L), ETAB(L), TI(L),   RI(L), QW(L), ZOFL(L) &
                  , AER1(L), NAA1(L), AER2(L), NAA2(L)
 
            write(6, '(i3,1x,2f11.7, 2x, f5.1, 1x, e8.2, f10.2, 2(f7.3,i4))') &
                 J, ETAA(L), ETAB(L), TI(L), QW(L), ZOFL(L) &
                  , AER1(L), NAA1(L), AER2(L), NAA2(L)
            
         enddo
         RI(L1_+1)=  0.d0
         TI(L1_+1)=  TI(L1_)
         write(6, '(A)') 'end of reading from CTM_GrdCld.dat'
         return 
         
91       stop 'error in opening CTM_GrdCld.dat file'
       END subroutine READ_1DCOL


     
      SUBROUTINE  INPUT_FOR_CLOUDJ(YLAT, MONTH, TI, RI, QW, ETAA, ETAB, PSURF, NAA1, AER1, NAA2, AER2, &
       PPP, ZZZ, TTT, RRR, DDD, OOO, NDXAER, AERSP)

        USE FJX_CMN_MOD, ONLY: L_, L1_, AN_, MASFAC, ZZHT
        USE RRSW_FASTJ_CMN, ONLY : PAVEL, TAVEL, PZ, TZ, ALTZ
        USE FJX_SUB_MOD, ONLY: ACLIM_FJX

        IMPLICIT NONE
        REAL*8,  INTENT(IN):: YLAT
        INTEGER, INTENT(IN):: MONTH
        REAL*8,  DIMENSION(L1_+1), INTENT(IN) :: TI, RI
        REAL*8,  DIMENSION(L1_), INTENT(IN):: QW
        REAL*8,  INTENT(IN):: PSURF 
        REAL*8,  DIMENSION(L1_),   INTENT(IN) :: ETAA, ETAB,AER1,AER2
        INTEGER, DIMENSION(L1_),   INTENT(IN) :: NAA1, NAA2
        REAL*8,  DIMENSION(L1_+1),   INTENT(OUT):: PPP,ZZZ
        REAL*8,  DIMENSION(L1_  ),   INTENT(OUT):: TTT,DDD,RRR, OOO
        REAL*8,  DIMENSION(L1_,AN_), INTENT(OUT):: AERSP
        INTEGER, DIMENSION(L1_,AN_), INTENT(OUT):: NDXAER
        INTEGER L
        REAL SCALEH       

        PPP(:)=0.d0 ! PPP(L1_+1) =0.D0 TOP OF ATMOSPHERE
!        print*, 'psurf=', psurf
        do L = 1, L1_ ! nlayer +1
           PPP(L) = ETAA(L) + ETAB(L)*PSURF
!           print*,'PPP=', etaa(L), Etab(L), PPP(L)
        enddo
!ACLIM_FJX sets up climatologies for O3 (because CTM_GrdCld.dat doesn't have ozone as input) 
!PPP is input from CTM; and one more level p=0 is added;
! set up OOO, colh2o, colch4, colco2, coldry and colo2 for rrtmg in ACLIM_FJX
        CALL ACLIM_FJX (YLAT, MONTH, PPP, QW, OOO, L1_)
! TTT,RRR is the averaged layer Temp and Humidity taking TI and RI from CTM input as the edge values 
        DO L = 1, L1_
           TTT(L) = 0.5d0*(TI(L)+TI(L+1))  
           RRR(L) = 0.5d0*(RI(L)+RI(L+1))
        ENDDO
        ZZZ(1)  = 16.d5*log10(1013.25d0/PPP(1))        ! zzz in cm
        DO L = 1, L1_
           DDD(L)  = (PPP(L)-PPP(L+1))*MASFAC
           SCALEH      = 1.3806d-19*MASFAC*TTT(L)
           ZZZ(L+1) = ZZZ(L) -(log(PPP(L+1))-log(PPP(L))) * SCALEH
!           write(6,'(2f8.2)')ppp(l), zzz(l)/1.e5 
        ENDDO
        ZZZ(L1_+1) = ZZZ(L1_) + ZZHT
! assign to the same variables used in RRTMG-Sw; L1_= L_+1 = LPAR +1       
        PZ(0:L1_) = PPP  ! EDGE PRESSURE
        TZ(0:L1_) = TI   ! EDGE TEMPERATURE
        ALTZ(0:L1_)= ZZZ ! EDGE HEIGHT
        TAVEL(1:L1_) = TTT ! LAYER TEMPERATURE
        DO L=1, L1_        !LAYER PRESSURE
           PAVEL(L) = 0.5D0*(PPP(L)+ PPP(L+1))
        END DO
!---load CTM-based data on ozone and aerosols on top of climatology
!---to convert kg (STT) in grid cell w/AREA (m2) to # molecules/cm^2
!      D_O3(I,J,L) = 6.023d26*STT(I,J,L,1)/48.d0  *1.d-4 /AREAXY(I,J)
!---to convert kg (STT) in grid cell w/AREA (m2) to PATH (g/m^2)
!      P_AERSL(I,J,L) = STT(I,J,L,1)*1.d-3/AREAXY(I,J)
!---this data should be available from the CTM somewhere.
!---   ZZZ(1:L_+1) is geometric altitude (cm), approx is OK.
!---AERSP = aerosol path (g/m2) and NDXAER() must come from CTM
!    the index must match one in the std or UMich data sets (or add your own).
        AERSP(:,:)  = 0.d0
        NDXAER(:,:) = 0
        do L = 1, L_
           NDXAER(L,1) = NAA1(L)
           AERSP(L,1)  = AER1(L)
           NDXAER(L,2) = NAA2(L)
           AERSP(L,2)  = AER2(L)
        enddo

        Return
      END Subroutine INPUT_FOR_CLOUDJ
     
      
      SUBROUTINE WRITE_INPUTMET_RRTMG(L1U, JULDAT, ALBEDO, SZA)
!
        USE FJX_CMN_MOD
        USE RRSW_FASTJ_CMN, ONLY : PAVEL, TAVEL, PZ, TZ, ALTZ, &
             COLCH4, COLH2O, COLO3, COLO2, COLCO2, COLDRY, NLAYERS
        USE PARRRSW ,ONLY : MXLAY
        USE FJX_SUB_MOD, only: JP_ATM0
        
!
        IMPLICIT NONE 
        INTEGER, INTENT(IN) ::  L1U, JULDAT
        REAL*8 , INTENT(IN)  :: ALBEDO, SZA 
!
!	  IAER   (0,10)   flag for aerosols
!                  = 0   no layers contain aerosols
!                  = 6   uses ECMWF global mean aerosol properties for one or all of
!                        six aerosol types.  Aerosol optical thickness at 0.55 micron,
!                        ECAER, must be set manually in the main source module, 
!                        rrtmg_sw.1col.f90, to activate the aerosols with this option.
!                  = 10  one or more layers contain aerosols
!                       (requires the presence of file IN_AER_RRTM)
!	  IATM   (0,1)   flag for RRTATM    1 = yes; 0 read from INPUT_RRTM_fastj
!         ISCAT  (0,1) switch for DISORT or simple two-stream scattering 
!                  = 0  DISORT     (unavailable)
!                  = 1  two-stream (default)
!          ISTRM   flag for number of streams used in DISORT  (ISCAT must be equal to 0)
!                  = 0  - 4 streams  (unavailable)
!                  = 1  - 8 streams  (unavailable)
!                  = 2  - 16 streams (unavailable)
!	  IOUT    = -1 if no output is to be printed out.
!	          =  0 if the only output is for 820-50000 cm-1.
!	          =  n (n = 16-29) if the only output is from band n.
!		       For the wavenumbers for each band, see Table I.
!	          = 98 if output is generated for 15 spectral intervals, one
!		       for the full shortwave spectrum (820-50000 cm-1), and one 
!		       for each of the 14 bands.
! 
!	  IMCA   (0,1) flag for McICA (Monte Carlo Independent Column Approximation)
!                       for statistical representation of sub-grid cloud fraction 
!                       and cloud overlap
!                  = 0  standard forward calculation; do not use McICA (valid for
!                       clear or overcast conditions only)
!                  = 1  use McICA (will perform statistical sample of 200 forward
!                       calculations and output average flux and heating rates)
!
!	  ICLD   (0,1,2,3) flag for clouds     
!                  = 0  no cloudy layers in atmosphere
!                  = 1  one or more cloudy layers present in atmosphere.  Cloud layers
!                       are treated as overcast only for IMCA = 0, or they are treated
!                       using a RANDOM overlap assumption for IMCA = 1.
!                       (requires the presence of file IN_CLD_RRTM for column model)
!                       (available for IMCA = 0 or 1)
!                  = 2  one or more cloudy layers present in atmosphere.  Cloud layers
!		       are treated using a MAXIMUM/RANDOM overlap assmption.
!                       (requires the presence of file IN_CLD_RRTM for column model)
!                       (available only for IMCA = 1)
!                  = 3  one or more cloudy layers present in atmosphere.  Cloud layers
!		       are treated using a MAXIMUM overlap assmption.
!                       (requires the presence of file IN_CLD_RRTM for column model)
!                       (available only for IMCA = 1)
!	  
        INTEGER, PARAMETER :: MXMOL= 7, IAER=0, IATM=0, ISCAT=1, ISTRM=0,&
             IOUT=98, IMCA= 0, ICLD=1, IDELM=0, ICOS=0, IB1=16, IB2=29
            
        REAL*8 :: WKL(MXMOL,MXLAY), WBRODL(MXLAY)
        REAL*8 :: SECNTK, CINP, SUMMOL
        INTEGER:: ISOLVAR, IPTHAK,L, M, IEMIS, IREFLECT, IFORM, NMOL
        REAL*8 :: SOLVAR(16:29), SEMISS(IB1:IB2)
        CHARACTER*80, FORM1(0:1), FORM2(0:1), FORM3(0:1)
!
        open(11,file='input_to_rrtmg.txt',form='formatted')
        iform =1 
        ipthak=1
        form1(0) = '(2f10.4,13x,i2,1x,2(f7.2,f8.3,f7.2))'
        form2(0) = '(2f10.4,13x,i2,23x,(f7.2,f8.3,f7.2))'
        form3(0) = '(8e10.3)'
        form1(1) = '(g15.7,g10.4,13x,i2,1x,2(f7.2,f8.3,f7.2))'
        form2(1) = '(g15.7,g10.4,13x,i2,23x,(f7.2,f8.3,f7.2))'
        form3(1) = '(1p 8e15.7)'                      
        write (11,'(A1)') '$'
        write (11,9011) iaer, iatm, iscat, istrm, iout, imca, icld, idelm, icos
!ISOLVAR      Solar variability option
!		     = 0 each band uses standard solar source function, corresponding
!                         to present day conditions.
!		     = 1 scale solar source function, each band will have the same
!                         scale factor applied, (equal to SOLVAR(16))
!		     = 2 scale solar source function, each band has different scale
!	 	         factors (for band IB, equal to SOLVAR(IB))	
!SOLVAR       the solar source function scale factor for each band.  If 
!	      ISOLVAR = 1, only the first value of SOLVAR (SOLVAR(16)) is 
!             considered.  If ISOLVAR = 2, values must be provided for all bands.
!
        solvar(:)=1.0d0
        isolvar=0
!        write (11,9020) juldat, sza, isolvar, solvar(ib1:ib2)
        write (11,9020) juldat, sza, isolvar
!
!IEMIS  = 0 each band has surface emissivity equal to 1.0
!       = 1 each band has the same surface emissivity (equal to SEMISS(16)) 
!       = 2 each band has different surface emissivity (for band IB, equal to SEMISS(IB))
        semiss(:)=0.d0
        iemis=1
        semiss(ib1)= 1.0d0 -albedo
!
        ireflect= 0 !for Lambertian reflection at surface

        write (11,9012)iemis, ireflect, semiss(ib1)
        nlayers= L1U
        nmol=7
        write (11,9013) iform, nlayers, nmol        
        wkl=0.d0
 !       print*,'colh2o(1) finalwrite', colh2o(1)
 !       print*,'colo2(1)=',  colo2(1)
 !       print*,'coldry(1)=', coldry(1)
        do l=1, L1U
           wkl(1,l) =  colh2o(l)/coldry(l)  
           wkl(2,l) =  colco2(l)/coldry(l)  
           wkl(3,l) =  colo3(l) /coldry(l)   
!          wkl(4,l) =  coln2o(l)/coldry(l)  
!          wkl(5,l)=?
           wkl(6,l) =  colch4(l)/coldry(l)  
           wkl(7,l) =  colo2(l) /coldry(l)   
        end do
        print*, 'wkl(1,1)=', wkl(1,1)
        summol=0.d0
        do l=1, L1U
!               summol=   (colco2(l)+ colo3(l)+ colch4(l)+colo2(l))/coldry(l)
           summol=0.d0
           wbrodl(l)= coldry(l)*(1.d0- summol)
        end do
        write (11, form1(1))&
             pavel(1),tavel(1), ipthak, altz(0)/1.d5, pz(0), &
             tz(0), altz(1)/1.d5, pz(1), tz(1)
        write (11,form3(iform)) (wkl(m,1),m=1,7), wbrodl(1)
        do l = 2, nlayers
           write(11,form2(iform)) pavel(l),tavel(l), &
                ipthak,altz(l)/1.d5, pz(l), tz(l)
           write(11,form3(iform)) (wkl(m,l),m=1,7), wbrodl(l)
        enddo
!
9011    format (18x,i2,29x,i1,32x,i1,1x,i1,2x,i3,3x,i1,i1,3x,i1,i1)
!9012    format (11x,i1,2x,i1,14f5.3)
9012    format (11x,i1, 2x, i1, f5.2)
9013    format (1x,i1,i3,i3)  
!9020    format (12x, i3, 3x, f7.4, 4x, i1, 14f7.5)
 9020    format (12x, i3, 3x, f7.3, 4x, i1)
! 
        return          
      end subroutine WRITE_INPUTMET_RRTMG
 
      SUBROUTINE WRITE_INPUTCLD_RRTMG(ICLD, LTOP, CLDFRAC, RIWP, RLWP, RI, RL)
!
        USE FJX_CMN_MOD, ONLY: L1_

        IMPLICIT NONE
        INTEGER, INTENT(IN):: ICLD, LTOP
        REAL*8, DIMENSION(L1_+1), INTENT(IN):: CLDFRAC, RLWP, RIWP, RI, RL
        REAL*8  FRAC, DAT1,DAT2, DAT3, DAT4
        LOGICAL EXIST
        CHARACTER*15  cldfile
        INTEGER INFLAG,ICEFLAG, LIQFLAG, L 
!   
!        print*,'icld=',icld
        if(icld .lt. 10) then
           write(cldfile, '(a14,i1)')'OUT_CLD_RRTMG0', icld
        else
           write(cldfile, '(a13,i2)')'OUT_CLD_RRTMG',  icld
        endif
        inquire(file=cldfile, exist=exist)
        if (exist) then
          print*,'file existed'
          return
        else
           open(13,file=cldfile,status='new',form='formatted', action='write')
        endif
!  inflag=2, iceflag= 1(Ebert and Curry, 97), 2(Key 2001), 3, Fu(1996) liqflag=1
        inflag= 2
        iceflag=3
        liqflag=1
        write(13,'(I5,I5,I5)')inflag, iceflag, liqflag
        do L=1, LTOP
        if(RIWP(L) .gt. 1.d-12 .or. RLWP(L) .gt. 1.d-12)then
           FRAC= CLDFRAC(L)
           DAT1= (RLWP(L)+ RIWP(L))  !g/m2 RLWP and RIWP 
           DAT2=  RIWP(L)/(RLWP(L)+RIWP(L)) !ice fraction
           DAT3=  RI(L) 
           DAT4=  RL(L)
           WRITE (13,9100) L,FRAC,DAT1,DAT2,DAT3,DAT4
        endif
        enddo
        write(13,'(a)')'%'
        close(13)
! in rrtmg
!      cldfrac(lay) = frac !cldfrac; layer cloud fraction
!      clddat1(lay) = dat1 !cwp
!      clddat2(lay) = dat2 !fice
!      clddat3(lay) = dat3 !=rei (cloud ice particle effective size)
!      clddat4(lay) = dat4 !=rel (cloud liquid particle effective radius (microns))

! dat1*dat2 = cwp* fice = ciwp(l)
! dat1*(1-dat2)= clwp(l)
9100    FORMAT (1x,1x,i3, f10.2, f10.6, f10.2, f10.2, f10.2)
        RETURN
      END SUBROUTINE WRITE_INPUTCLD_RRTMG
!------------------------------------------------------------

      SUBROUTINE OUTPUT_HEATINGRATE(L1U, iprof, ncldflag, PPP, ZZZ, FJ, SJ)
        USE FJX_CMN_MOD, ONLY: S_, X_, JVN_, NW1, NW2, NS2,&
             TITLEJX, WL, EXIST, ParaSummary        
        implicit none
        integer, intent(in) :: L1U, iprof, ncldflag 
        real*8 , intent(in), dimension(L1U+1)    :: PPP, ZZZ
        real*8 , intent(in), dimension(L1U-1, JVN_):: FJ
        
        real*8 , intent(in), dimension(S_+2 , L1U):: SJ

        integer I, L, K, NJ

        inquire(file='sj_1.dat', exist= exist)
        if (exist)then
           open(8,file='sj_1.dat', status='old', form='formatted', position='append', action='write')
        else
           open(8,file='sj_1.dat', status='new', form='formatted', action='write')
        endif
        write(8,'(a)')ParaSummary
        write(8,'(a11,I3,a15,I3)')'%profile # ', iprof, ' ;cloud flag=  ', ncldflag
        write(8,'(a)')'%heating rate profiles in K/day v7.4 180-778nm'
        write(8,'(a5, 18f7.1)')'%wvl ',(WL(I),I=NW1,NW2)
        do L = L1U, 1, -1
           write(8,'(i3, 1x, f9.4, 18 f11.5)') L, PPP(L),(SJ(I,L), I=NW1,NW2)
        enddo
        close(8)
!
        inquire(file='sj_2.dat', exist= exist)
        if (exist)then
           open(9,file='sj_2.dat', status='old', form='formatted', position='append', action='write')
        else
           open(9,file='sj_2.dat', status='new', form='formatted', action='write')
        endif
        write(9, '(a)') ParaSummary
        write(9,'(a11,I3,a15,I3)')'%profile # ', iprof, ' ;cloud flag=  ', ncldflag
        write(9, '(a)')'%heating rate profiles in K/day v7.4 778nm-.12um plus 1:18  19:27 & 1:27'
        write(9, '(a5,18f7.1)')'%wvl ',(WL(I),I=NW2+1, NS2)
        do L = L1U, 1, -1
           write(9,'(i4, 1x, f9.4, 18f11.5)') L, PPP(L) ,(SJ(I,L), I=NW2+1,NS2+2),  &
                SJ(NS2+1,L)+SJ(NS2+2,L)
        enddo
        close(9)
!
        inquire(file= 'fj.dat', exist= exist)
        if (exist)then
           open(10,file='fj.dat', status='old', form='formatted', position='append', action='write')
        else
           open(10,file='fj.dat', status='new', form='formatted', action='write')
        endif
        write(10,'(a)') ParaSummary
        write(10,'(a11,I3,a15,I3)')'%profile # ', iprof, ' ;cloud flag=  ', ncldflag
!        write(10,'(a)') '%fast-JX (7.2)----J-values----'
        write(10,'(a)') '%fast-JX (7.2)----JO3& JNO2----'
!       write(10,'(8(9(a6,3x)/))') (TITLEJX(K), K=1,X_)
        do L = L1U-1, 1, -1
        write(10,'(i4, 1x, f9.4, f9.4, 2es13.5)') L, (zzz(L)+zzz(L+1))/2.e5, 16.*log10(1013.25/((ppp(L)+ppp(L+1))/2.0)), FJ(L,2), FJ(L,9)  
!         write(10,'(A4, i3, 1x, f7.2/8(1p 9e9.2/))')'%L= ', L, PPP(L) ,(FJ(L,K),K=1,X_)
        enddo
        close(10)
        return
      end subroutine OUTPUT_HEATINGRATE

      SUBROUTINE OUTPUT_OD(iprof, ncldflag, PPP, OD)
        USE FJX_CMN_MOD, ONLY: L_, S_, NW1, NW2, NS2,&
             WL, EXIST, ParaSummary        
        implicit none
        integer, intent(in) :: iprof, ncldflag 
        real*8 , intent(in), dimension(L_+1)    :: PPP
        real*8 , intent(in), dimension(S_+2 , L_):: OD

        integer I, L, K, NJ

        inquire(file='OD_1.dat', exist= exist)
        if (exist)then
           open(81,file='OD_1.dat', status='old', form='formatted', position='append', action='write')
        else
           open(81,file='OD_1.dat', status='new', form='formatted', action='write')
        endif
        write(81,'(a)')ParaSummary
        write(81,'(a11,I3,a15,I3)')'%profile # ', iprof, ' ;cloud flag=  ', ncldflag
        write(81,'(a)')'%Layer Optical Depth v7.4 180-778nm'
        write(81,'(a5, 18f7.1)')'%wvl ',(WL(I),I=NW1,NW2)
        do L = L_,1,-1
           write(81,'(i3,1x, 18f8.4)') L, (OD(I,L), I=NW1,NW2)
        enddo
        close(81)
!
        inquire(file='OD_2.dat', exist= exist)
        if (exist)then
           open(91,file='OD_2.dat', status='old', form='formatted', position='append', action='write')
        else
           open(91,file='OD_2.dat', status='new', form='formatted', action='write')
        endif
        write(91, '(a)') ParaSummary
        write(91,'(a11,I3,a15,I3)')'%profile # ', iprof, ' ;cloud flag=  ', ncldflag
        write(91, '(a)')'%heating rate profiles in K/day v7.4 778nm-.12um plus 1:18  19:27 & 1:27'
        write(91, '(a5, 18f7.1)')'%wvl ',(WL(I),I=NW2+1,NS2)

        do L = L_,1,-1
           write(91,'(i4, 18f8.4)') L, (OD(I,L), I=NW2+1,NS2)
        enddo
        close(91)
!
        return
      end subroutine OUTPUT_OD


        SUBROUTINE OUTPUT_ICEOD(ODc, SSAc, ODs, ODa, Phaft)

        USE FJX_CMN_MOD, ONLY: L_, S_, NW1, NW2, NS2,&
             WL, EXIST, ParaSummary, RecTle        
        implicit none
        real*8 , intent(in), dimension(S_ , L_):: ODc, SSAc, ODs, ODa, Phaft
!    
        real*8 :: sumt(18:27)
        integer I, L, K, NJ

        inquire(file='CirrusOptics.dat', exist= exist)

        if (exist)then
           open(81,file='CirrusOptics.dat', status='old', form='formatted', position='append', action='write')
        else
           open(81,file='CirrusOptics.dat', status='new', form='formatted', action='write')
        endif
        write(81, '(a)')RecTle !record icld and icldflag
        write(81,'(a)')ParaSummary
        write(81,'(a)')'% Scattering Optical Depth from bins 18-27'
        write(81,'(a5, 10f7.1)')'%wvl ',(WL(I),I=18, 27)
        write(81,'(a)')'cirrus optical depth'
        sumt(18:27)=0.d0
        do L = 34,21,-1
           write(81,'(i3,1x, 10f8.4)') L, (ODc(I,L), I=18,27)
           do i=18,27
           sumt(i)= sumt(i) + ODc(i,L)
           enddo
        enddo
        write(81,'(A4, 10f8.4)')'sum=', (sumt(I), I=18,27)

        write(81,'(a)')'cirrus ssa'
        do L = 34,21,-1
           write(81,'(i3,1x, 10f8.4)') L, (SSAc(I,L), I=18,27)
        enddo
        sumt(18:27)=0.d0
        write(81,'(a)')'scattering OD'
        do L = 34,21,-1
           write(81,'(i3,1x, 10f8.4)') L, (ODs(I,L), I=18,27)
           do i=18,27
              sumt(i)= sumt(i) + ODs(i,L)
           enddo
        enddo
        write(81,'(A4, 10f8.4)')'sum=', (sumt(I), I=18,27)
!
        sumt(18:27)=0.d0
        write(81,'(a)')'absorption OD'
        do L = 34,21,-1
           write(81,'(i3,1x, 10f8.4)') L, (ODa(I,L), I=18,27)
           do i=18,27
              sumt(i)= sumt(i) + ODa(i,L)
           enddo
        enddo
        write(81,'(A4, 10f8.4)')'sum=', (sumt(I), I=18,27)
        write(81,'(a)')'first phase function P(2)/3'
        do L = 34,21,-1
           write(81,'(i4, 10f12.4)') L, (Phaft(I,L), I=18,27)
        enddo

        close(81)

        end subroutine OUTPUT_ICEOD


      Subroutine CLD_ERRSTATS(VALJXX,PPP,CLDFLAG, JP04, JP09, JP11, JP15, ZPJAVG, ERRJ2, ERRJJ, JJCNT, CLDFLDS, LFINISH, ZZZ)
        USE FJX_CMN_MOD, only: L_, L1_, JVN_, JIND
        IMPLICIT NONE
        REAL*8,  DIMENSION(L_,JVN_), INTENT(IN) :: VALJXX
        REAL*8,  DIMENSION(L1_+1), INTENT(IN):: PPP
        INTEGER, INTENT(IN) :: CLDFLAG
        INTEGER, INTENT(OUT):: JP04,JP09, JP11, JP15      
        REAL*8,  DIMENSION(L_,8,4),INTENT(INOUT):: ZPJAVG
        REAL*8,  DIMENSION(21,4), INTENT(INOUT) :: ERRJJ, ERRJ2
        INTEGER, INTENT(IN):: JJCNT(8), CLDFLDS
        LOGICAL, INTENT(IN):: LFINISH
        REAL*8,  INTENT(IN):: ZZZ(L1_+1) 
        REAL*8,  DIMENSION(L_,8,4)  :: ZPJCLD
        REAL*8,  DIMENSION(21,4):: ERRJ
        REAL*8   ZKM
        INTEGER  ICLD, I, J, L, K
        CHARACTER*11, DIMENSION(4) ::  TITJX
        character*8, dimension(8), parameter :: TITERR =  &
             ['clearsky','avgcloud','cldf^3/2','ICA-beam','ICAs-ran','QCAs-mid',&
             'QCAs-avg','all ICAs']
!
!   4O3        PHOTON    O2        O(1D)     1.000 mapped to FJX:   3 O3(1D)
!   9NO2       PHOTON    N2        O         1.000 mapped to FJX:   9 NO2
!  11NO3       PHOTON    NO        O2        0.114 mapped to FJX:  10 NO3
!  15HNO3      PHOTON    NO2       OH        1.000 mapped to FJX:  13 HNO3
            JP04 = JIND(4)
            JP09 = JIND(9)
            JP11 = JIND(11)
            JP15 = JIND(15)
            TITJX(1) = 'J(O3>O1D)  '
            TITJX(2) = 'J(NO2)     '
            TITJX(3) = 'J(NO3>all) '
            TITJX(4) = 'J(HNO3)    '
            do L = 1,L_
               ZPJCLD(L,CLDFLAG,1) = VALJXX(L,JP04)
               ZPJCLD(L,CLDFLAG,2) = VALJXX(L,JP09)
               ZPJCLD(L,CLDFLAG,3) = VALJXX(L,JP11)
               ZPJCLD(L,CLDFLAG,4) = VALJXX(L,JP15)
               ZPJAVG(L,CLDFLAG,1) = ZPJAVG(L,CLDFLAG,1) + VALJXX(L,JP04)
               ZPJAVG(L,CLDFLAG,2) = ZPJAVG(L,CLDFLAG,2) + VALJXX(L,JP09)
               ZPJAVG(L,CLDFLAG,3) = ZPJAVG(L,CLDFLAG,3) + VALJXX(L,JP11)
               ZPJAVG(L,CLDFLAG,4) = ZPJAVG(L,CLDFLAG,4) + VALJXX(L,JP15)
            enddo
            if(CLDFLAG .eq. 8)then
! look at errors at L=1, L=34 (above clouds), and p-wtd mean (1:33)
               do J=1,4
                  do I=1,8
                     ZPJCLD(L_,I,J) = 0.d0
                     do L = 1,33
                        ZPJCLD(L_,I,J) = ZPJCLD(L_,I,J) + &
                                  ZPJCLD(L,I,J)*(PPP(L)-PPP(L+1))
                     enddo
                     ZPJCLD(L_,I,J) = ZPJCLD(L_,I,J)/(PPP(1)-PPP(34))
                  enddo
               enddo
               L=1
               do J=1,4
                  ERRJ(1,J) = ZPJCLD(L,1,J)/ZPJCLD(L,8,J)
                  ERRJ(2,J) = ZPJCLD(L,2,J)/ZPJCLD(L,8,J)
                  ERRJ(3,J) = ZPJCLD(L,3,J)/ZPJCLD(L,8,J)
                  ERRJ(4,J) = ZPJCLD(L,4,J)/ZPJCLD(L,8,J)
                  ERRJ(5,J) = ZPJCLD(L,5,J)/ZPJCLD(L,8,J)
                  ERRJ(6,J) = ZPJCLD(L,6,J)/ZPJCLD(L,8,J)
                  ERRJ(7,J) = ZPJCLD(L,7,J)/ZPJCLD(L,8,J)
               enddo
               L=34
               do J=1,4
                  ERRJ( 8,J) = ZPJCLD(L,1,J)/ZPJCLD(L,8,J)
                  ERRJ( 9,J) = ZPJCLD(L,2,J)/ZPJCLD(L,8,J)
                  ERRJ(10,J) = ZPJCLD(L,3,J)/ZPJCLD(L,8,J)
                  ERRJ(11,J) = ZPJCLD(L,4,J)/ZPJCLD(L,8,J)
                  ERRJ(12,J) = ZPJCLD(L,5,J)/ZPJCLD(L,8,J)
                  ERRJ(13,J) = ZPJCLD(L,6,J)/ZPJCLD(L,8,J)
                  ERRJ(14,J) = ZPJCLD(L,7,J)/ZPJCLD(L,8,J)
               enddo
               L=L_
               do J=1,4
                  ERRJ(15,J) = ZPJCLD(L,1,J)/ZPJCLD(L,8,J)
                  ERRJ(16,J) = ZPJCLD(L,2,J)/ZPJCLD(L,8,J)
                  ERRJ(17,J) = ZPJCLD(L,3,J)/ZPJCLD(L,8,J)
                  ERRJ(18,J) = ZPJCLD(L,4,J)/ZPJCLD(L,8,J)
                  ERRJ(19,J) = ZPJCLD(L,5,J)/ZPJCLD(L,8,J)
                  ERRJ(20,J) = ZPJCLD(L,6,J)/ZPJCLD(L,8,J)
                  ERRJ(21,J) = ZPJCLD(L,7,J)/ZPJCLD(L,8,J)
               enddo
! accumulate errors for each cloud field
               do J=1,4
                  do I=1,21
                     ERRJ(I,J) = log(ERRJ(I,J))
                     ERRJJ(I,J) = ERRJJ(I,J) + ERRJ(I,J)
                     ERRJ2(I,J) = ERRJ2(I,J) + ERRJ(I,J)**2
                  enddo
               enddo
            endif
! print summary over all CLDFLAGS and all atmospheres (ICLDs)
            if(LFINISH)THEN
            do J=1,4
               do L=1,34
                  do K=1,8
                     ZPJAVG(L,K,J) = ZPJAVG(L,K,J)/float(CLDFLDS)
                  enddo
               enddo
            enddo
            do J=1,4
               write(6,'(a11,i7,8i9)') TITJX(J),(K, K=1,8)
               do L=34,1,-1
                  ZKM = 1.d-5*ZZZ(L)
                  write(6,'(i3,f6.1,1p,72e9.2)') L,ZKM,(ZPJAVG(L,K,J),K=1,8)
               enddo
            enddo
            do J=1,4
               do I=1,21
                  ERRJJ(I,J) = ERRJJ(I,J)/float(CLDFLDS)
                  ERRJ2(I,J) = sqrt(ERRJ2(I,J)/float(CLDFLDS))
               enddo
            enddo
            write(6,'(5a)') ' mean error (%)  L=1  &  L=34  & p-avg ',(TITJX(J),J=1,4)
            write(6,'(a11,i8,2p,4f7.0,5x,4f7.0,5x,4f7.0)')  (TITERR(I),JJCNT(I),  &
                 (ERRJJ(I,J),J=1,4),(ERRJJ(I+7,J),J=1,4),(ERRJJ(I+14,J),J=1,4), I=1,7)
            write(6,'(a11,i8)')  TITERR(8),JJCNT(8)
            write(6,'(5a)') ' RMS error (%)   L=1  &  L=34  & p-avg ',(TITJX(J),J=1,4)
            write(6,'(a11,i8,2p,4f7.0,5x,4f7.0,5x,4f7.0)') ,(TITERR(I),JJCNT(I),  &
                 (ERRJ2(I,J),J=1,4),(ERRJ2(I+7,J),J=1,4),(ERRJ2(I+14,J),J=1,4), I=1,7)
            write(6,'(a11,i8)')  TITERR(8),JJCNT(8)
         endif


! accumulate mean J's

         return
         end subroutine CLD_ERRSTATS
