! $Id $

  ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved.

  ! This code is available without charge with the following conditions:

  ! 1. The code is available for scientific purposes and is not for

  ! commercial use.

  ! 2. Any improvements you make to the code should be made available

  ! to the to the authors for incorporation into a future release.

  ! 3. The code should not be used in any way that brings the authors

  ! or their employers into disrepute.

  IMPLICIT NONE

  ! NOTE:   the maximum number of levels and columns is set by

  ! the following parameter statement

  INTEGER ncolprint

  ! -----

  ! Input

  ! -----

  INTEGER npoints !  number of model points in the horizontal

  ! PARAMETER(npoints=6722)

  INTEGER nlev !  number of model levels in column

  INTEGER ncol !  number of subcolumns

  INTEGER sunlit(npoints) !  1 for day points, 0 for night time

  INTEGER seed(npoints) !  seed value for random number generator

  ! !  ( see Numerical Recipes Chapter 7)

  ! !  It is recommended that the seed is set

  ! !  to a different value for each model

  ! !  gridbox it is called on, as it is

  ! !  possible that the choice of the samec

  ! !  seed value every time may introduce some

  ! !  statistical bias in the results, particularly

  ! !  for low values of NCOL.

  REAL pfull(npoints, nlev) !  pressure of full model levels (Pascals)

  ! !  pfull(npoints,1)    is    top level of model

  ! !  pfull(npoints,nlev) is bottom level of model

  REAL phalf(npoints, nlev+1) !  pressure of half model levels (Pascals)

  ! !  phalf(npoints,1)    is    top       of model

  ! !  phalf(npoints,nlev+1) is the surface pressure

  REAL qv(npoints, nlev) !  water vapor specific humidity (kg vapor/ kg air)

  ! !         on full model levels

  REAL cc(npoints, nlev) !  input cloud cover in each model level (fraction)

  ! !  NOTE:  This is the HORIZONTAL area of each

  ! !         grid box covered by clouds

  REAL conv(npoints, nlev) !  input convective cloud cover in each model level (fraction)

  ! !  NOTE:  This is the HORIZONTAL area of each

  ! !         grid box covered by convective clouds

  REAL dtau_s(npoints, nlev) !  mean 0.67 micron optical depth of stratiform

  ! !  clouds in each model level

  ! !  NOTE:  this the cloud optical depth of only the

  ! !         cloudy part of the grid box, it is not weighted

  ! !         with the 0 cloud optical depth of the clear

  ! !         part of the grid box

  REAL dtau_c(npoints, nlev) !  mean 0.67 micron optical depth of convective

  ! !  clouds in each

  ! !  model level.  Same note applies as in dtau_s.

  INTEGER overlap !  overlap type

  ! 1=max

  ! 2=rand

  ! 3=max/rand

  INTEGER top_height !  1 = adjust top height using both a computed

  ! !  infrared brightness temperature and the visible

  ! !  optical depth to adjust cloud top pressure. Note

  ! !  that this calculation is most appropriate to compare

  ! !  to ISCCP data during sunlit hours.

  ! !  2 = do not adjust top height, that is cloud top

  ! !  pressure is the actual cloud top pressure

  ! !  in the model

  ! !  3 = adjust top height using only the computed

  ! !  infrared brightness temperature. Note that this

  ! !  calculation is most appropriate to compare to ISCCP

  ! !  IR only algortihm (i.e. you can compare to nighttime

  ! !  ISCCP data with this option)

  REAL tautab(0:255) !  ISCCP table for converting count value to

  ! !  optical thickness

  INTEGER invtau(-20:45000) !  ISCCP table for converting optical thickness

  ! !  to count value

  ! The following input variables are used only if top_height = 1 or

  ! top_height = 3

  REAL skt(npoints) !  skin Temperature (K)

  REAL emsfc_lw !  10.5 micron emissivity of surface (fraction)

  REAL at(npoints, nlev) !  temperature in each model level (K)

  REAL dem_s(npoints, nlev) !  10.5 micron longwave emissivity of stratiform

  ! !  clouds in each

  ! !  model level.  Same note applies as in dtau_s.

  REAL dem_c(npoints, nlev) !  10.5 micron longwave emissivity of convective

  ! !  clouds in each

  ! !  model level.  Same note applies as in dtau_s.

  ! IM reg.dyn BEG

  REAL t1, t2

  ! REAL w(npoints)                   !vertical wind at 500 hPa

  ! LOGICAL pct_ocean(npoints)        !TRUE if oceanic point, FALSE otherway

  ! INTEGER iw(npoints) , nw

  ! REAL wmin, pas_w

  ! INTEGER k, l, iwmx

  ! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30)

  ! REAL fq_dynreg(7,7,iwmx)

  ! REAL nfq_dynreg(7,7,iwmx)

  ! LOGICAL pctj(7,7,iwmx)

  ! IM reg.dyn END

  ! ------

  ! Output

  ! ------

  REAL fq_isccp(npoints, 7, 7) !  the fraction of the model grid box covered by

  ! !  each of the 49 ISCCP D level cloud types

  REAL totalcldarea(npoints) !  the fraction of model grid box columns

  ! !  with cloud somewhere in them.  This should

  ! !  equal the sum over all entries of fq_isccp

  ! ! The following three means are averages over the cloudy areas only.  If

  ! no

  ! ! clouds are in grid box all three quantities should equal zero.

  REAL meanptop(npoints) !  mean cloud top pressure (mb) - linear averaging

  ! !  in cloud top pressure.

  REAL meantaucld(npoints) !  mean optical thickness

  ! !  linear averaging in albedo performed.

  REAL boxtau(npoints, ncol) !  optical thickness in each column

  REAL boxptop(npoints, ncol) !  cloud top pressure (mb) in each column

  ! ------

  ! Working variables added when program updated to mimic Mark Webb's PV-Wave

  ! code

  ! ------

  ! ! Equivalent of BOX in original version, but

  ! ! indexed by column then row, rather than

  ! ! by row then column

  ! ! with extra layer of zeroes on top

  ! ! in this version this just contains the values input

  ! ! from cc but with an extra level

  ! ! from conv

  ! ! is chosen

  ! !  emissivity in part of

  ! !  gridbox under consideration

  INTEGER acc(nlev, ncol)

  ! !variables needed for water vapor continuum absorption

  CHARACTER *1 cchar(6), cchar_realtops(6)

  INTEGER icycle

  REAL boxarea

  INTEGER debug ! set to non-zero value to print out inputs

  ! ! with step debug

  INTEGER debugcol ! set to non-zero value to print out column

  ! ! decomposition with step debugcol

  INTEGER index1(npoints), num1, jj

  REAL rec2p13, tauchk

  CHARACTER *10 ftn09

  DATA isccp_taumin/0.3/

  DATA cchar/' ', '-', '1', '+', 'I', '+'/

  DATA cchar_realtops/' ', ' ', '1', '1', 'I', 'I'/

  tauchk = -1.*log(0.9999999)

  rec2p13 = 1./2.13

  ncolprint = 0

  ! IM

  ! PRINT*,' isccp_cloud_types npoints=',npoints

  ! if ( debug.ne.0 ) then

  ! j=1

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write(6,'(a10)') 'debug='

  ! write(6,'(8I10)') debug

  ! write(6,'(a10)') 'debugcol='

  ! write(6,'(8I10)') debugcol

  ! write(6,'(a10)') 'npoints='

  ! write(6,'(8I10)') npoints

  ! write(6,'(a10)') 'nlev='

  ! write(6,'(8I10)') nlev

  ! write(6,'(a10)') 'ncol='

  ! write(6,'(8I10)') ncol

  ! write(6,'(a10)') 'top_height='

  ! write(6,'(8I10)') top_height

  ! write(6,'(a10)') 'overlap='

  ! write(6,'(8I10)') overlap

  ! write(6,'(a10)') 'emsfc_lw='

  ! write(6,'(8f10.2)') emsfc_lw

  ! write(6,'(a10)') 'tautab='

  ! write(6,'(8f10.2)') tautab

  ! write(6,'(a10)') 'invtau(1:100)='

  ! write(6,'(8i10)') (invtau(i),i=1,100)

  ! do j=1,npoints,debug

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write(6,'(a10)') 'sunlit='

  ! write(6,'(8I10)') sunlit(j)

  ! write(6,'(a10)') 'seed='

  ! write(6,'(8I10)') seed(j)

  ! write(6,'(a10)') 'pfull='

  ! write(6,'(8f10.2)') (pfull(j,i),i=1,nlev)

  ! write(6,'(a10)') 'phalf='

  ! write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1)

  ! write(6,'(a10)') 'qv='

  ! write(6,'(8f10.3)') (qv(j,i),i=1,nlev)

  ! write(6,'(a10)') 'cc='

  ! write(6,'(8f10.3)') (cc(j,i),i=1,nlev)

  ! write(6,'(a10)') 'conv='

  ! write(6,'(8f10.2)') (conv(j,i),i=1,nlev)

  ! write(6,'(a10)') 'dtau_s='

  ! write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev)

  ! write(6,'(a10)') 'dtau_c='

  ! write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev)

  ! write(6,'(a10)') 'skt='

  ! write(6,'(8f10.2)') skt(j)

  ! write(6,'(a10)') 'at='

  ! write(6,'(8f10.2)') (at(j,i),i=1,nlev)

  ! write(6,'(a10)') 'dem_s='

  ! write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev)

  ! write(6,'(a10)') 'dem_c='

  ! write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev)

  ! enddo

  ! endif

  ! ---------------------------------------------------!

  ! assign 2d tca array using 1d input array cc

  END DO

    END DO

  END DO

  ! assign 2d cca array using 1d input array conv

    ! IM pas besoin        do ibox=1,ncol

    END DO

    ! IM        enddo

  END DO

  ! IM

  ! do j=1, iwmx

  ! do l=1, 7

  ! do k=1, 7

  ! fq_dynreg(k,l,j) =0.

  ! nfq_dynreg(k,l,j) =0.

  ! enddo !k

  ! enddo !l

  ! enddo !j

  ! IM

  ! IM

  ! if (ncolprint.ne.0) then

  ! do j=1,npoints,1000

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write (6,'(a)') 'seed:'

  ! write (6,'(I3.2)') seed(j)

  ! write (6,'(a)') 'tca_pp_rev:'

  ! write (6,'(8f5.2)')

  ! &   ((tca(j,ilev)),

  ! &      ilev=1,nlev)

  ! write (6,'(a)') 'cca_pp_rev:'

  ! write (6,'(8f5.2)')

  ! &   ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev)

  ! enddo

  ! endif

    END DO

            at(j,ilev)<attropmin(j)) THEN

        END IF

      END DO

    END DO

  END IF

  ! -----------------------------------------------------!

  ! ---------------------------------------------------!

  ! IM

  ! do 13 ilev=1,nlev

  ! num1=0

  ! do j=1,npoints

  ! if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)  ' error = cloud fraction less than zero'

  ! write(6,*) ' or '

  ! write(6,*)  ' error = cloud fraction greater than 1'

  ! write(6,*) 'value at point ',j,' is ',cc(j,ilev)

  ! write(6,*) 'level ',ilev

  ! STOP

  ! enddo

  ! num1=0

  ! do j=1,npoints

  ! if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)

  ! &           ' error = convective cloud fraction less than zero'

  ! write(6,*) ' or '

  ! write(6,*)

  ! &           ' error = convective cloud fraction greater than 1'

  ! write(6,*) 'value at point ',j,' is ',conv(j,ilev)

  ! write(6,*) 'level ',ilev

  ! STOP

  ! enddo

  ! num1=0

  ! do j=1,npoints

  ! if (dtau_s(j,ilev) .lt. 0.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)

  ! &           ' error = stratiform cloud opt. depth less than zero'

  ! write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev)

  ! write(6,*) 'level ',ilev

  ! STOP

  ! enddo

  ! num1=0

  ! do j=1,npoints

  ! if (dtau_c(j,ilev) .lt. 0.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)

  ! &           ' error = convective cloud opt. depth less than zero'

  ! write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev)

  ! write(6,*) 'level ',ilev

  ! STOP

  ! enddo

  ! num1=0

  ! do j=1,npoints

  ! if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)

  ! &           ' error = stratiform cloud emissivity less than zero'

  ! write(6,*)'or'

  ! write(6,*)

  ! &           ' error = stratiform cloud emissivity greater than 1'

  ! write(6,*) 'value at point ',j,' is ',dem_s(j,ilev)

  ! write(6,*) 'level ',ilev

  ! STOP

  ! enddo

  ! num1=0

  ! do j=1,npoints

  ! if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then

  ! num1=num1+1

  ! index1(num1)=j

  ! end if

  ! enddo

  ! do jj=1,num1

  ! j=index1(jj)

  ! write(6,*)

  ! &           ' error = convective cloud emissivity less than zero'

  ! write(6,*)'or'

  ! write(6,*)

  ! &           ' error = convective cloud emissivity greater than 1'

  ! write (6,*)

  ! &          'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev)

  ! STOP

  ! enddo

  ! 13    continue

    END DO

  END DO

  ! ---------------------------------------------------!

  ! Initialise working variables

  ! ---------------------------------------------------!

  ! Initialised frac_out to zero

      END DO

    END DO

  END DO

  ! IM

  ! if (ncolprint.ne.0) then

  ! write (6,'(a)') 'frac_out_pp_rev:'

  ! do j=1,npoints,1000

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write (6,'(8f5.2)')

  ! &     ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)

  ! enddo

  ! write (6,'(a)') 'ncol:'

  ! write (6,'(I3)') ncol

  ! endif

  ! if (ncolprint.ne.0) then

  ! write (6,'(a)') 'last_frac_pp:'

  ! do j=1,npoints,1000

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write (6,'(8f5.2)') (tca(j,0))

  ! enddo

  ! endif

  ! ---------------------------------------------------!

  ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS

  ! frac_out is the array that contains the information

  ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a

  ! convective cloud

    !loop over vertical levels

    ! Initialise threshold

        ! If max overlap

          ! select pixels spread evenly

          ! across the gridbox

          END DO

        END DO

      ELSE

            ! select random pixels from the non-convective

            ! part the gridbox ( some will be converted into

            ! convective pixels below )

          END DO

        END DO

      END IF

      ! IM

      ! IF (ncolprint.ne.0) then

      ! write (6,'(a)') 'threshold_nsf2:'

      ! do j=1,npoints,1000

      ! write(6,'(a10)') 'j='

      ! write(6,'(8I10)') j

      ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)

      ! enddo

      ! ENDIF

    END IF

    ! IF (ncolprint.ne.0) then

    ! write (6,'(a)') 'ilev:'

    ! write (6,'(I2)') ilev

    ! ENDIF

        ! All versions

          ! IM REAL           maxocc(j,ibox) = 1

        ELSE

          ! IM REAL           maxocc(j,ibox) = 0

        END IF

      END DO

        ! Max overlap

          ! IM REAL           maxosc(j,ibox)=1

        END DO

      END IF

        ! Random overlap

          ! IM REAL           maxosc(j,ibox)=0

        END DO

      END IF

        ! Max/Random overlap

          threshold_min(j, ibox) = max(cca(j,ilev), min(tca(j,ilev-1),tca(j, &

          IF (threshold(j,ibox)<min(tca(j,ilev-1),tca(j, &

            ! IM REAL                maxosc(j,ibox)= 1

          ELSE

            ! IM REAL                 maxosc(j,ibox)= 0

          END IF

        END DO

      END IF

        ! Reset threshold

        threshold(j, ibox) = &       !if max overlapped conv cloud

          maxocc(j, ibox)*(boxpos(j,ibox)) + &   !else

          (1-maxocc(j,ibox))*( &     !if max overlapped strat cloud

          (maxosc(j,ibox))*( &       !threshold=boxpos

          threshold(j,ibox))+ &      !else

          (1-maxosc(j,ibox))*( &     !threshold_min=random[thrmin,1]

      END DO

    END DO ! ibox

    ! Fill frac_out with 1's where tca is greater than the threshold

          ! IM REAL             frac_out(j,ibox,ilev)=1

        ELSE

          ! IM REAL             frac_out(j,ibox,ilev)=0

        END IF

      END DO

    END DO

    ! Code to partition boxes into startiform and convective parts

    ! goes here

            ! = 2 IF threshold le cca(j)

          ! IM REAL                  frac_out(j,ibox,ilev) = 2

        ELSE

            ! = the same IF NOT threshold le cca(j)

          frac_out(j, ibox, ilev) = frac_out(j, ibox, ilev)

        END IF

      END DO

    END DO

    ! Set last_frac to tca at this level, so as to be tca

    ! from last level next time round

    ! IM

    ! if (ncolprint.ne.0) then

    ! do j=1,npoints ,1000

    ! write(6,'(a10)') 'j='

    ! write(6,'(8I10)') j

    ! write (6,'(a)') 'last_frac:'

    ! write (6,'(8f5.2)') (tca(j,ilev-1))

    ! write (6,'(a)') 'cca:'

    ! write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)

    ! write (6,'(a)') 'max_overlap_cc:'

    ! write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)

    ! write (6,'(a)') 'max_overlap_sc:'

    ! write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)

    ! write (6,'(a)') 'threshold_min_nsf2:'

    ! write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)

    ! write (6,'(a)') 'threshold_nsf2:'

    ! write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)

    ! write (6,'(a)') 'frac_out_pp_rev:'

    ! write (6,'(8f5.2)')

    ! &       ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)

    ! enddo

    ! endif

  END DO

  ! ---------------------------------------------------!

  ! ---------------------------------------------------!

  ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and

  ! put into vector tau

    !initialize tau and albedocld to zero

  ! loop over nlev

    END DO

  END DO

    !compute total cloud optical depth for each column

      !increment tau for each of the boxes

        ! IM REAL              if (frac_out(j,ibox,ilev).eq.1) then

        END IF

        ! IM REAL              if (frac_out(j,ibox,ilev).eq.2) then

        END IF

      END DO

    END DO ! ibox

  END DO ! ilev

  ! IM

  ! if (ncolprint.ne.0) then

  ! do j=1,npoints ,1000

  ! write(6,'(a10)') 'j='

  ! write(6,'(8I10)') j

  ! write(6,'(i2,1X,8(f7.2,1X))')

  ! &          ilev,

  ! &          (tau(j,ibox),ibox=1,ncolprint)

  ! enddo

  ! endif

  ! ---------------------------------------------------!

  ! ---------------------------------------------------!

  ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES

  ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE

  ! again this is only done if top_height = 1 or 3

  ! fluxtop is the 10.5 micron radiance at the top of the

  ! atmosphere

  ! trans_layers_above is the total transmissivity in the layers

  ! above the current layer

  ! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear

  ! sky versions of these quantities.

      !----------------------------------------------------------------------

      !

      !             DO CLEAR SKY RADIANCE CALCULATION FIRST

      !

      !compute water vapor continuum emissivity

      !this treatment follows Schwarkzopf and Ramasamy

      !JGR 1999,vol 104, pages 9467-9499.

      !the emissivity is calculated at a wavenumber of 955 cm-1,

      !or 10.47 microns

    wtmair = 28.9644

    wtmh20 = 18.01534

    navo = 6.023E+23

    grav = 9.806650E+02

    pstd = 1.013250E+06

    t0 = 296.

    ! IM

    ! if (ncolprint .ne. 0)

    ! &         write(6,*)  'ilev   pw (kg/m2)   tauwv(j)      dem_wv'

          !press and dpress are dyne/cm2 = Pascals *10

          !atmden = g/cm2 = kg/m2 / 10

        tauwv(j) = wk(j)*1.E-20*((0.0224697*rh20s(j)*tmpexp(j))+(3.41817E-7* &

      END DO

      ! IM

      ! if (ncolprint .ne. 0) then

      ! do j=1,npoints ,1000

      ! write(6,'(a10)') 'j='

      ! write(6,'(8I10)') j

      ! write(6,'(i2,1X,3(f8.3,3X))') ilev,

      ! &           qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.),

      ! &           tauwv(j),dem_wv(j,ilev)

      ! enddo

      ! endif

    END DO

      !initialize variables

    END DO

          ! Black body emission at temperature of the layer

          !bb(j)= 5.67e-8*at(j,ilev)**4

          ! increase TOA flux by flux emitted from layer

          ! times total transmittance in layers above

        fluxtop_clrsky(j) = fluxtop_clrsky(j) + dem_wv(j, ilev)*bb(j)* &

          ! update trans_layers_above with transmissivity

          ! from this layer for next time around loop

        trans_layers_above_clrsky(j) = trans_layers_above_clrsky(j)* &

      END DO

      ! IM

      ! if (ncolprint.ne.0) then

      ! do j=1,npoints ,1000

      ! write(6,'(a10)') 'j='

      ! write(6,'(8I10)') j

      ! write (6,'(a)') 'ilev:'

      ! write (6,'(I2)') ilev

      ! write (6,'(a)')

      ! &        'emiss_layer,100.*bb(j),100.*f,total_trans:'

      ! write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j),

      ! &             100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j)

      ! enddo

      ! endif

    END DO !loop over level

        !add in surface emission

        !bb(j)=5.67e-8*skt(j)**4

      fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw*bb(j)* &

    END DO

    ! IM

    ! if (ncolprint.ne.0) then

    ! do j=1,npoints ,1000

    ! write(6,'(a10)') 'j='

    ! write(6,'(8I10)') j

    ! write (6,'(a)') 'id:'

    ! write (6,'(a)') 'surface'

    ! write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:'

    ! write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j),

    ! &      100.*fluxtop_clrsky(j),

    ! &       trans_layers_above_clrsky(j)

    ! enddo

    ! endif

      !

      !           END OF CLEAR SKY CALCULATION

      !

      !----------------------------------------------------------------

    ! IM

    ! if (ncolprint.ne.0) then

    ! do j=1,npoints ,1000

    ! write(6,'(a10)') 'j='

    ! write(6,'(8I10)') j

    ! write (6,'(a)') 'ts:'

    ! write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint)

    ! write (6,'(a)') 'ta_rev:'

    ! write (6,'(8f7.2)')

    ! &       ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)

    ! enddo

    ! endif

      !loop over columns

      END DO

    END DO

          ! Black body emission at temperature of the layer

          !bb(j)= 5.67e-8*at(j,ilev)**4

      END DO

            ! emissivity for point in this layer

          ! IM REAL             if (frac_out(j,ibox,ilev).eq.1) then

            ! IM REAL             else if (frac_out(j,ibox,ilev).eq.2) then

          ELSE

          END IF

            ! increase TOA flux by flux emitted from layer

            ! times total transmittance in layers above

          fluxtop(j, ibox) = fluxtop(j, ibox) + dem(j, ibox)*bb(j)* &

            ! update trans_layers_above with transmissivity

            ! from this layer for next time around loop

          trans_layers_above(j, ibox) = trans_layers_above(j, ibox)* &

        END DO ! j

      END DO ! ibox

      ! IM

      ! if (ncolprint.ne.0) then

      ! do j=1,npoints,1000

      ! write (6,'(a)') 'ilev:'

      ! write (6,'(I2)') ilev

      ! write(6,'(a10)') 'j='

      ! write(6,'(8I10)') j

      ! write (6,'(a)') 'emiss_layer:'

      ! write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)

      ! write (6,'(a)') '100.*bb(j):'

      ! write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint)

      ! write (6,'(a)') '100.*f:'

      ! write (6,'(8f7.2)')

      ! &         (100.*fluxtop(j,ibox),ibox=1,ncolprint)

      ! write (6,'(a)') 'total_trans:'

      ! write (6,'(8f7.2)')

      ! &          (trans_layers_above(j,ibox),ibox=1,ncolprint)

      ! enddo

      ! endif

    END DO ! ilev