module icarus_scops_38 2,5
!-----------------------------------------------------------------------
!Purpose: Create module for:
! Name: ISCCP Simulator ICARUS/SCOPS
! What: Simulate ISCCP cloud products from GCM inputs
! Version: 3.8
! Authors: Steve Klein (klein21@llnl.gov)
! Mark Webb (Mark.Webb@MetOffice.com)
!
!Author: B. Eaton March 2004
!Author: S. Klein September 2008 (update simulator to version 3.8)
!-----------------------------------------------------------------------
use shr_kind_mod
, only: r8 => shr_kind_r8
use abortutils, only: endrun
use marsaglia, only: kissvec
use units
use cam_logfile, only: iulog
implicit none
private
save
! Public interfaces
public :: &
isccp_cloud_types
! Public data
integer, public, parameter :: &
ntau = 7, &! number of optical depth ranges
npres = 7 ! number of pressure ranges
real(r8), public, parameter :: &
prlim(npres+1) = (/0._r8, 180._r8, 310._r8, 440._r8, 560._r8, 680._r8, 800._r8, 1000._r8/), &
taulim(ntau+1) = (/0._r8, 0.3_r8, 1.3_r8, 3.6_r8, 9.4_r8, 23._r8, 60._r8, 379._r8/)
!##############################################################################
CONTAINS
!##############################################################################
!########################################################################################################
!
! Name: ISCCP Simulator ICARUS/SCOPS
! What: Simulate ISCCP cloud products from GCM inputs
! Version: 3.8
! Authors: Steve Klein (klein21@llnl.gov)
! Mark Webb (Mark.Webb@MetOffice.com)
!
! Modifications to downloaded code:
! 1. convert to free format
! 2. convert types of real arguments to real
! 4. add diagnostic meanttop
!
! Steve Klein modified the CAM code for version 3.4. for updates to
! version 3.8
! New diagnostic output variables meantaucld, meanalbedocld, meantb, meantbclr
SUBROUTINE ISCCP_CLOUD_TYPES( & 2,20
debug, &
debugcol, &
npoints, &
sunlit, &
nlev, &
ncol, &
seed1, &
seed2, &
seed3, &
seed4, &
pfull, &
phalf, &
qv, &
cc, &
conv, &
dtau_s, &
dtau_c, &
top_height, &
top_height_direction, &
overlap, &
skt, &
emsfc_lw, &
at, &
dem_s, &
dem_c, &
fq_isccp, &
totalcldarea, &
meanptop, &
meanttop, &
meantaucld, &
meanalbedocld, &
meantb, &
meantbclr, &
boxtau, &
boxptop &
)
! CVSId: isccp_cloud_types.f,v 3.4 2003/05/06 08:45:28 hadmw Exp
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
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 seed1(npoints) ! seed value for random number generator
INTEGER seed2(npoints) ! It is recommended that the seed is set
INTEGER seed3(npoints) ! to a different value for each model
INTEGER seed4(npoints) ! gridbox it is called on, as it is
! possible that the choice of the same
! seed value every time may introduce some
! statistical bias in the results, particularly
! for low values of NCOL.
! four integer seeds are needed by generator
REAL(R8) 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(R8) 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(R8) qv(npoints,nlev) ! water vapor specific humidity (kg vapor/ kg air)
! on full model levels
REAL(R8) 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(R8) 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(R8) 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(R8) 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)
INTEGER top_height_direction ! direction for finding atmosphere pressure level
! with interpolated temperature equal to the radiance
! determined cloud-top temperature
!
! 1 = find the *lowest* altitude (highest pressure) level
! with interpolated temperature equal to the radiance
! determined cloud-top temperature
!
! 2 = find the *highest* altitude (lowest pressure) level
! with interpolated temperature equal to the radiance
! determined cloud-top temperature
!
! ONLY APPLICABLE IF top_height EQUALS 1 or 3
!
! 1 = default setting, and matches all versions of
! ISCCP simulator with versions numbers 3.5.1 and lower
!
! 2 = experimental setting
!
! The following input variables are used only if top_height = 1 or top_height = 3
!
REAL(R8) skt(npoints) ! skin Temperature (K)
REAL(R8) emsfc_lw ! 10.5 micron emissivity of surface (fraction)
REAL(R8) at(npoints,nlev) ! temperature in each model level (K)
REAL(R8) dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform
! clouds in each
! model level. Same note applies as in dtau_s.
REAL(R8) dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective
! clouds in each
! model level. Same note applies as in dtau_s.
! ------
! Output
! ------
REAL(R8) fq_isccp(npoints,7,7) ! the fraction of the model grid box covered by
! each of the 49 ISCCP D level cloud types
REAL(R8) totalcldarea(npoints) ! the fraction of model grid box columns
! with cloud somewhere in them. NOTE: This diagnostic
! does not count model clouds with tau < isccp_taumin
!
! Thus this diagnostic does not equal the sum over all entries of fq_isccp.
! However, this diagnostic does equal the sum over entries of fq_isccp with
! itau = 2:7 (omitting itau = 1)
! The following four means are averages only over the cloudy areas with tau > isccp_taumin.
! If no clouds with tau > isccp_taumin are in grid box all three quantities should equal zero.
REAL(R8) meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging
! in cloud top pressure.
REAL(R8) meanttop(npoints) ! mean cloud top temp (k) - linear averaging
REAL(R8) meantaucld(npoints) ! mean optical thickness
! linear averaging in albedo performed.
REAL(R8) meanalbedocld(npoints) ! mean cloud albedo
! linear averaging in albedo performed.
REAL(R8) meantb(npoints) ! mean all-sky 10.5 micron brightness temperature (K)
REAL(R8) meantbclr(npoints) ! mean clear-sky 10.5 micron brightness temperature (K)
REAL(R8) boxtau(npoints,ncol) ! optical thickness in each column
REAL(R8) boxptop(npoints,ncol) ! cloud top pressure (mb) in each column
!
! ------
! Working variables added when program updated to mimic Mark Webb's PV-Wave code
! ------
REAL(R8) frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into
! Equivalent of BOX in original version, but
! indexed by column then row, rather than
! by row then column
REAL(R8) tca(npoints,0:nlev) ! total cloud cover in each model level (fraction)
! with extra layer of zeroes on top
! in this version this just contains the values input
! from cc but with an extra level
REAL(R8) cca(npoints,nlev) ! convective cloud cover in each model level (fraction)
! from conv
REAL(R8) threshold(npoints,ncol) ! pointer to position in gridbox
REAL(R8) maxocc(npoints,ncol) ! Flag for max overlapped conv cld
REAL(R8) maxosc(npoints,ncol) ! Flag for max overlapped strat cld
REAL(R8) boxpos(npoints,ncol) ! ordered pointer to position in gridbox
REAL(R8) threshold_min(npoints,ncol) ! minimum value to define range in with new threshold
! is chosen
REAL(R8) dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave
! emissivity in part of
! gridbox under consideration
REAL(R8) ran(npoints) ! vector of random numbers
REAL(R8) ptrop(npoints)
REAL(R8) attrop(npoints)
REAL(R8) attropmin (npoints)
REAL(R8) atmax(npoints)
REAL(R8) atmin(npoints)
REAL(R8) btcmin(npoints)
REAL(R8) transmax(npoints)
INTEGER i,j,ilev,ibox,itrop(npoints)
INTEGER ipres(npoints)
INTEGER itau(npoints),ilev2
INTEGER acc(nlev,ncol)
INTEGER match(npoints,nlev-1)
INTEGER nmatch(npoints)
INTEGER levmatch(npoints,ncol)
!variables needed for water vapor continuum absorption
real(r8) fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints)
real(r8) taumin(npoints)
real(r8) dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0
real(r8) press(npoints), dpress(npoints), atmden(npoints)
real(r8) rvh20(npoints), wk(npoints), rhoave(npoints)
real(r8) rh20s(npoints), rfrgn(npoints)
real(r8) tmpexp(npoints),tauwv(npoints)
character*1 cchar(6),cchar_realtops(6)
integer icycle
REAL(R8) tau(npoints,ncol)
LOGICAL box_cloudy(npoints,ncol)
REAL(R8) tb(npoints,ncol)
REAL(R8) ptop(npoints,ncol)
REAL(R8) ttop(npoints,ncol)
REAL(R8) emcld(npoints,ncol)
REAL(R8) fluxtop(npoints,ncol)
REAL(R8) trans_layers_above(npoints,ncol)
real(r8) isccp_taumin,fluxtopinit(npoints),tauir(npoints)
REAL(R8) albedocld(npoints,ncol)
real(r8) 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,k1,k2
real(r8) rec2p13,tauchk,logp,logp1,logp2,atd
real(r8) output_missing_value
integer unitn ! output file unit number
character*10 ftn09
DATA isccp_taumin / 0.3_r8 /
DATA output_missing_value / -1.E+30 /
DATA cchar / ' ','-','1','+','I','+'/
DATA cchar_realtops / ' ',' ','1','1','I','I'/
tauchk = -1._r8*log(0.9999999_r8)
rec2p13=1._r8/2.13_r8
ncolprint=0
if ( debug.ne.0 ) then
j=1
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a10)') 'debug='
write(iulog,'(8I10)') debug
write(iulog,'(a10)') 'debugcol='
write(iulog,'(8I10)') debugcol
write(iulog,'(a10)') 'npoints='
write(iulog,'(8I10)') npoints
write(iulog,'(a10)') 'nlev='
write(iulog,'(8I10)') nlev
write(iulog,'(a10)') 'ncol='
write(iulog,'(8I10)') ncol
write(iulog,'(a11)') 'top_height='
write(iulog,'(8I10)') top_height
write(iulog,'(a21)') 'top_height_direction='
write(iulog,'(8I10)') top_height_direction
write(iulog,'(a10)') 'overlap='
write(iulog,'(8I10)') overlap
write(iulog,'(a10)') 'emsfc_lw='
write(iulog,'(8f10.2)') emsfc_lw
do j=1,npoints,debug
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a10)') 'sunlit='
write(iulog,'(8I10)') sunlit(j)
write(iulog,'(a10)') 'seed1='
write(iulog,'(8I10)') seed1(j)
write(iulog,'(a10)') 'seed2='
write(iulog,'(8I10)') seed2(j)
write(iulog,'(a10)') 'seed3='
write(iulog,'(8I10)') seed3(j)
write(iulog,'(a10)') 'seed4='
write(iulog,'(8I10)') seed4(j)
write(iulog,'(a10)') 'pfull='
write(iulog,'(8f10.2)') (pfull(j,i),i=1,nlev)
write(iulog,'(a10)') 'phalf='
write(iulog,'(8f10.2)') (phalf(j,i),i=1,nlev+1)
write(iulog,'(a10)') 'qv='
write(iulog,'(8f10.3)') (qv(j,i),i=1,nlev)
write(iulog,'(a10)') 'cc='
write(iulog,'(8f10.3)') (cc(j,i),i=1,nlev)
write(iulog,'(a10)') 'conv='
write(iulog,'(8f10.2)') (conv(j,i),i=1,nlev)
write(iulog,'(a10)') 'dtau_s='
write(iulog,'(8g12.5)') (dtau_s(j,i),i=1,nlev)
write(iulog,'(a10)') 'dtau_c='
write(iulog,'(8f10.2)') (dtau_c(j,i),i=1,nlev)
write(iulog,'(a10)') 'skt='
write(iulog,'(8f10.2)') skt(j)
write(iulog,'(a10)') 'at='
write(iulog,'(8f10.2)') (at(j,i),i=1,nlev)
write(iulog,'(a10)') 'dem_s='
write(iulog,'(8f10.3)') (dem_s(j,i),i=1,nlev)
write(iulog,'(a10)') 'dem_c='
write(iulog,'(8f10.2)') (dem_c(j,i),i=1,nlev)
enddo
endif
! ---------------------------------------------------!
! assign 2d tca array using 1d input array cc
do j=1,npoints
tca(j,0)=0
enddo
do ilev=1,nlev
do j=1,npoints
tca(j,ilev)=cc(j,ilev)
enddo
enddo
! assign 2d cca array using 1d input array conv
do ilev=1,nlev
do ibox=1,ncol
do j=1,npoints
cca(j,ilev)=conv(j,ilev)
enddo
enddo
enddo
if (ncolprint.ne.0) then
do j=1,npoints,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'seed1:'
write(iulog,'(I3.2)') seed1(j)
write(iulog,'(a)') 'seed2:'
write(iulog,'(I3.2)') seed2(j)
write(iulog,'(a)') 'seed3:'
write(iulog,'(I3.2)') seed3(j)
write(iulog,'(a)') 'seed4:'
write(iulog,'(I3.2)') seed4(j)
write(iulog,'(a)') 'tca_pp_rev:'
write(iulog,'(8f5.2)') &
((tca(j,ilev)), &
ilev=1,nlev)
write(iulog,'(a)') 'cca_pp_rev:'
write(iulog,'(8f5.2)') &
((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev)
enddo
endif
if (top_height .eq. 1 .or. top_height .eq. 3) then
do j=1,npoints
ptrop(j)=5000._r8
atmin(j) = 400._r8
attropmin(j) = 400._r8
atmax(j) = 0._r8
attrop(j) = 120._r8
itrop(j) = 1
enddo
do 12 ilev=1,nlev
do j=1,npoints
if (pfull(j,ilev) .lt. 40000._r8 .and. &
pfull(j,ilev) .gt. 5000._r8 .and. &
at(j,ilev) .lt. attropmin(j)) then
ptrop(j) = pfull(j,ilev)
attropmin(j) = at(j,ilev)
attrop(j) = attropmin(j)
itrop(j)=ilev
end if
if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev)
if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev)
enddo
12 continue
end if
if (top_height .eq. 1 .or. top_height .eq. 3) then
do j=1,npoints
meantb(j) = 0._r8
meantbclr(j) = 0._r8
end do
else
do j=1,npoints
meantb(j) = output_missing_value
meantbclr(j) = output_missing_value
end do
end if
! -----------------------------------------------------!
! ---------------------------------------------------!
do 13 ilev=1,nlev
num1=0
do j=1,npoints
if (cc(j,ilev) .lt. 0._r8 .or. cc(j,ilev) .gt. 1._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) ' error = cloud fraction less than zero'
write(iulog,*) ' or '
write(iulog,*) ' error = cloud fraction greater than 1'
write(iulog,*) 'value at point ',j,' is ',cc(j,ilev)
write(iulog,*) 'level ',ilev
call endrun('ISCCP_CLOUD_TYPES ERROR: cloud fraction')
enddo
num1=0
do j=1,npoints
if (conv(j,ilev) .lt. 0._r8 .or. conv(j,ilev) .gt. 1._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) &
' error = convective cloud fraction less than zero'
write(iulog,*) ' or '
write(iulog,*) &
' error = convective cloud fraction greater than 1'
write(iulog,*) 'value at point ',j,' is ',conv(j,ilev)
write(iulog,*) 'level ',ilev
call endrun('ISCCP_CLOUD_TYPES ERROR: convective cloud fraction')
enddo
num1=0
do j=1,npoints
if (dtau_s(j,ilev) .lt. 0._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) &
' error = stratiform cloud opt. depth less than zero'
write(iulog,*) 'value at point ',j,' is ',dtau_s(j,ilev)
write(iulog,*) 'level ',ilev
call endrun('ISCCP_CLOUD_TYPES ERROR: stratiform cloud opt. depth')
enddo
num1=0
do j=1,npoints
if (dtau_c(j,ilev) .lt. 0._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) &
' error = convective cloud opt. depth less than zero'
write(iulog,*) 'value at point ',j,' is ',dtau_c(j,ilev)
write(iulog,*) 'level ',ilev
call endrun('ISCCP_CLOUD_TYPES ERROR: convective cloud opt. depth')
enddo
num1=0
do j=1,npoints
if (dem_s(j,ilev) .lt. 0._r8 .or. dem_s(j,ilev) .gt. 1._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) &
' error = stratiform cloud emissivity less than zero'
write(iulog,*)'or'
write(iulog,*) &
' error = stratiform cloud emissivity greater than 1'
write(iulog,*) 'value at point ',j,' is ',dem_s(j,ilev)
write(iulog,*) 'level ',ilev
call endrun('ISCCP_CLOUD_TYPES ERROR: stratiform cloud emissivity')
enddo
num1=0
do j=1,npoints
if (dem_c(j,ilev) .lt. 0._r8 .or. dem_c(j,ilev) .gt. 1._r8) then
num1=num1+1
index1(num1)=j
end if
enddo
do jj=1,num1
j=index1(jj)
write(iulog,*) &
' error = convective cloud emissivity less than zero'
write(iulog,*)'or'
write(iulog,*) &
' error = convective cloud emissivity greater than 1'
write(iulog,*) &
'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev)
call endrun('ISCCP_CLOUD_TYPES ERROR: convective cloud emissivity')
enddo
13 continue
do ibox=1,ncol
do j=1,npoints
boxpos(j,ibox)=(ibox-.5_r8)/ncol
enddo
enddo
! ---------------------------------------------------!
! Initialise working variables
! ---------------------------------------------------!
! Initialised frac_out to zero
do ibox=1,ncol
do ilev=1,nlev
do j=1,npoints
frac_out(j,ibox,ilev)=0.0_r8
enddo
enddo
enddo
if (ncolprint.ne.0) then
write(iulog,'(a)') 'frac_out_pp_rev:'
do j=1,npoints,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(8f5.2)') &
((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev)
enddo
write(iulog,'(a)') 'ncol:'
write(iulog,'(I3)') ncol
endif
if (ncolprint.ne.0) then
write(iulog,'(a)') 'last_frac_pp:'
do j=1,npoints,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(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
DO 200 ilev = 1,nlev
! Initialise threshold
IF (ilev.eq.1) then
! If max overlap
IF (overlap.eq.1) then
! select pixels spread evenly
! across the gridbox
DO ibox=1,ncol
do j=1,npoints
threshold(j,ibox)=boxpos(j,ibox)
enddo
enddo
ELSE
DO ibox=1,ncol
call kissvec(seed1,seed2,seed3,seed4,ran)
! select random pixels from the non-convective
! part the gridbox ( some will be converted into
! convective pixels below )
do j=1,npoints
threshold(j,ibox)= &
cca(j,ilev)+(1-cca(j,ilev))*ran(j)
enddo
enddo
ENDIF
IF (ncolprint.ne.0) then
write(iulog,'(a)') 'threshold_nsf2:'
do j=1,npoints,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
enddo
ENDIF
ENDIF
IF (ncolprint.ne.0) then
write(iulog,'(a)') 'ilev:'
write(iulog,'(I2)') ilev
ENDIF
DO ibox=1,ncol
! All versions
do j=1,npoints
if (boxpos(j,ibox).le.cca(j,ilev)) then
maxocc(j,ibox) = 1
else
maxocc(j,ibox) = 0
end if
enddo
! Max overlap
if (overlap.eq.1) then
do j=1,npoints
threshold_min(j,ibox)=cca(j,ilev)
maxosc(j,ibox)=1
enddo
endif
! Random overlap
if (overlap.eq.2) then
do j=1,npoints
threshold_min(j,ibox)=cca(j,ilev)
maxosc(j,ibox)=0
enddo
endif
! Max/Random overlap
if (overlap.eq.3) then
do j=1,npoints
threshold_min(j,ibox)=max(cca(j,ilev), &
min(tca(j,ilev-1),tca(j,ilev)))
if (threshold(j,ibox) &
.lt.min(tca(j,ilev-1),tca(j,ilev)) &
.and.(threshold(j,ibox).gt.cca(j,ilev))) then
maxosc(j,ibox)= 1
else
maxosc(j,ibox)= 0
end if
enddo
endif
! Reset threshold
call kissvec(seed1,seed2,seed3,seed4,ran)
do j=1,npoints
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]
threshold_min(j,ibox)+ &
(1-threshold_min(j,ibox))*ran(j) &
) &
)
enddo
ENDDO ! ibox
! Fill frac_out with 1's where tca is greater than the threshold
DO ibox=1,ncol
do j=1,npoints
if (tca(j,ilev).gt.threshold(j,ibox)) then
frac_out(j,ibox,ilev)=1
else
frac_out(j,ibox,ilev)=0
end if
enddo
ENDDO
! Code to partition boxes into startiform and convective parts
! goes here
DO ibox=1,ncol
do j=1,npoints
if (threshold(j,ibox).le.cca(j,ilev)) then
! = 2 IF threshold le cca(j)
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
enddo
ENDDO
! Set last_frac to tca at this level, so as to be tca
! from last level next time round
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'last_frac:'
write(iulog,'(8f5.2)') (tca(j,ilev-1))
write(iulog,'(a)') 'cca:'
write(iulog,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint)
write(iulog,'(a)') 'max_overlap_cc:'
write(iulog,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'max_overlap_sc:'
write(iulog,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'threshold_min_nsf2:'
write(iulog,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'threshold_nsf2:'
write(iulog,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'frac_out_pp_rev:'
write(iulog,'(8f5.2)') &
((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
enddo
endif
200 CONTINUE !loop over nlev
!
! ---------------------------------------------------!
!
! ---------------------------------------------------!
! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and
! put into vector tau
!initialize tau and albedocld to zero
do 15 ibox=1,ncol
do j=1,npoints
tau(j,ibox)=0._r8
albedocld(j,ibox)=0._r8
boxtau(j,ibox)=output_missing_value
boxptop(j,ibox)=output_missing_value
box_cloudy(j,ibox)=.false.
enddo
15 continue
!compute total cloud optical depth for each column
do ilev=1,nlev
!increment tau for each of the boxes
do ibox=1,ncol
do j=1,npoints
if (frac_out(j,ibox,ilev).eq.1) then
tau(j,ibox)=tau(j,ibox) &
+ dtau_s(j,ilev)
endif
if (frac_out(j,ibox,ilev).eq.2) then
tau(j,ibox)=tau(j,ibox) &
+ dtau_c(j,ilev)
end if
enddo
enddo ! ibox
enddo ! ilev
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(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.
if (top_height .eq. 1 .or. top_height .eq. 3) then
!----------------------------------------------------------------------
!
! 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_r8
wtmh20 = 18.01534_r8
Navo = 6.023E+23_r8
grav = 9.806650E+02_r8
pstd = 1.013250E+06_r8
t0 = 296._r8
if (ncolprint .ne. 0) &
write(iulog,*) 'ilev pw (kg/m2) tauwv(j) dem_wv'
do 125 ilev=1,nlev
do j=1,npoints
!press and dpress are dyne/cm2 = Pascals *10
press(j) = pfull(j,ilev)*10._r8
dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10
!atmden = g/cm2 = kg/m2 / 10
atmden(j) = dpress(j)/grav
rvh20(j) = qv(j,ilev)*wtmair/wtmh20
wk(j) = rvh20(j)*Navo*atmden(j)/wtmair
rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev))
rh20s(j) = rvh20(j)*rhoave(j)
rfrgn(j) = rhoave(j)-rh20s(j)
tmpexp(j) = exp(-0.02_r8*(at(j,ilev)-t0))
tauwv(j) = wk(j)*1.e-20_r8*( &
(0.0224697_r8*rh20s(j)*tmpexp(j)) + &
(3.41817e-7_r8*rfrgn(j)) )*0.98_r8
dem_wv(j,ilev) = 1._r8 - exp( -1._r8 * tauwv(j))
enddo
if (ncolprint .ne. 0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(i2,1X,3(f8.3,3X))') ilev, &
qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100._r8), &
tauwv(j),dem_wv(j,ilev)
enddo
endif
125 continue
!initialize variables
do j=1,npoints
fluxtop_clrsky(j) = 0._r8
trans_layers_above_clrsky(j)=1._r8
enddo
do ilev=1,nlev
do j=1,npoints
! Black body emission at temperature of the layer
bb(j)=1 / ( exp(1307.27_r8/at(j,ilev)) - 1._r8 )
!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)*trans_layers_above_clrsky(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)*(1._r8-dem_wv(j,ilev))
enddo
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'ilev:'
write(iulog,'(I2)') ilev
write(iulog,'(a)') &
'emiss_layer,100.*bb(j),100.*f,total_trans:'
write(iulog,'(4(f7.2,1X))') dem_wv(j,ilev),100._r8*bb(j), &
100._r8*fluxtop_clrsky(j),trans_layers_above_clrsky(j)
enddo
endif
enddo !loop over level
do j=1,npoints
!add in surface emission
bb(j)=1/( exp(1307.27_r8/skt(j)) - 1._r8 )
!bb(j)=5.67e-8*skt(j)**4
fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j) &
* trans_layers_above_clrsky(j)
!clear sky brightness temperature
meantbclr(j) = 1307.27_r8/(log(1._r8+(1._r8/fluxtop_clrsky(j))))
enddo
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'id:'
write(iulog,'(a)') 'surface'
write(iulog,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:'
write(iulog,'(5(f7.2,1X))') emsfc_lw,100._r8*bb(j), &
100._r8*fluxtop_clrsky(j), &
trans_layers_above_clrsky(j), meantbclr(j)
enddo
endif
!
! END OF CLEAR SKY CALCULATION
!
!----------------------------------------------------------------
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'ts:'
write(iulog,'(8f7.2)') (skt(j),ibox=1,ncolprint)
write(iulog,'(a)') 'ta_rev:'
write(iulog,'(8f7.2)') &
((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev)
enddo
endif
!loop over columns
do ibox=1,ncol
do j=1,npoints
fluxtop(j,ibox)=0._r8
trans_layers_above(j,ibox)=1._r8
enddo
enddo
do ilev=1,nlev
do j=1,npoints
! Black body emission at temperature of the layer
bb(j)=1 / ( exp(1307.27_r8/at(j,ilev)) - 1._r8 )
!bb(j)= 5.67e-8*at(j,ilev)**4
enddo
do ibox=1,ncol
do j=1,npoints
! emissivity for point in this layer
if (frac_out(j,ibox,ilev).eq.1) then
dem(j,ibox)= 1._r8 - &
( (1._r8 - dem_wv(j,ilev)) * (1._r8 - dem_s(j,ilev)) )
else if (frac_out(j,ibox,ilev).eq.2) then
dem(j,ibox)= 1._r8 - &
( (1._r8 - dem_wv(j,ilev)) * (1._r8 - dem_c(j,ilev)) )
else
dem(j,ibox)= dem_wv(j,ilev)
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) &
* trans_layers_above(j,ibox)
! update trans_layers_above with transmissivity
! from this layer for next time around loop
trans_layers_above(j,ibox)= &
trans_layers_above(j,ibox)*(1._r8-dem(j,ibox))
enddo ! j
enddo ! ibox
if (ncolprint.ne.0) then
do j=1,npoints,1000
write(iulog,'(a)') 'ilev:'
write(iulog,'(I2)') ilev
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'emiss_layer:'
write(iulog,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') '100.*bb(j):'
write(iulog,'(8f7.2)') (100._r8*bb(j),ibox=1,ncolprint)
write(iulog,'(a)') '100.*f:'
write(iulog,'(8f7.2)') &
(100._r8*fluxtop(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'total_trans:'
write(iulog,'(8f7.2)') &
(trans_layers_above(j,ibox),ibox=1,ncolprint)
enddo
endif
enddo ! ilev
do j=1,npoints
!add in surface emission
bb(j)=1/( exp(1307.27_r8/skt(j)) - 1._r8 )
!bb(j)=5.67e-8*skt(j)**4
end do
do ibox=1,ncol
do j=1,npoints
!add in surface emission
fluxtop(j,ibox) = fluxtop(j,ibox) &
+ emsfc_lw * bb(j) &
* trans_layers_above(j,ibox)
end do
end do
!calculate mean infrared brightness temperature
do ibox=1,ncol
do j=1,npoints
meantb(j) = meantb(j)+1307.27_r8/(log(1._r8+(1._r8/fluxtop(j,ibox))))
end do
end do
do j=1, npoints
meantb(j) = meantb(j) / real(ncol)
end do
if (ncolprint.ne.0) then
do j=1,npoints ,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'id:'
write(iulog,'(a)') 'surface'
write(iulog,'(a)') 'emiss_layer:'
write(iulog,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint)
write(iulog,'(a)') '100.*bb(j):'
write(iulog,'(8f7.2)') (100._r8*bb(j),ibox=1,ncolprint)
write(iulog,'(a)') '100.*f:'
write(iulog,'(8f7.2)') (100._r8*fluxtop(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'meantb(j):'
write(iulog,'(8f7.2)') (meantb(j),ibox=1,ncolprint)
end do
endif
!now that you have the top of atmosphere radiance account
!for ISCCP procedures to determine cloud top temperature
!account for partially transmitting cloud recompute flux
!ISCCP would see assuming a single layer cloud
!note choice here of 2.13, as it is primarily ice
!clouds which have partial emissivity and need the
!adjustment performed in this section
!
!If it turns out that the cloud brightness temperature
!is greater than 260K, then the liquid cloud conversion
!factor of 2.56 is used.
!
!Note that this is discussed on pages 85-87 of
!the ISCCP D level documentation (Rossow et al. 1996)
do j=1,npoints
!compute minimum brightness temperature and optical depth
btcmin(j) = 1._r8 / ( exp(1307.27_r8/(attrop(j)-5._r8)) - 1._r8 )
enddo
do ibox=1,ncol
do j=1,npoints
transmax(j) = (fluxtop(j,ibox)-btcmin(j)) &
/(fluxtop_clrsky(j)-btcmin(j))
!note that the initial setting of tauir(j) is needed so that
!tauir(j) has a realistic value should the next if block be
!bypassed
tauir(j) = tau(j,ibox) * rec2p13
taumin(j) = -1._r8 * log(max(min(transmax(j),0.9999999_r8),0.001_r8))
enddo
if (top_height .eq. 1) then
do j=1,npoints
if (transmax(j) .gt. 0.001_r8 .and. &
transmax(j) .le. 0.9999999_r8) then
fluxtopinit(j) = fluxtop(j,ibox)
tauir(j) = tau(j,ibox) *rec2p13
endif
enddo
do icycle=1,2
do j=1,npoints
if (tau(j,ibox) .gt. (tauchk )) then
if (transmax(j) .gt. 0.001_r8 .and. &
transmax(j) .le. 0.9999999_r8) then
emcld(j,ibox) = 1._r8 - exp(-1._r8 * tauir(j) )
fluxtop(j,ibox) = fluxtopinit(j) - &
((1._r8-emcld(j,ibox))*fluxtop_clrsky(j))
fluxtop(j,ibox)=max(1.E-06_r8, &
(fluxtop(j,ibox)/emcld(j,ibox)))
tb(j,ibox)= 1307.27_r8 &
/ (log(1._r8 + (1._r8/fluxtop(j,ibox))))
if (tb(j,ibox) .gt. 260._r8) then
tauir(j) = tau(j,ibox) / 2.56_r8
end if
end if
end if
enddo
enddo
endif
do j=1,npoints
if (tau(j,ibox) .gt. (tauchk )) then
!cloudy box
!NOTE: tb is the cloud-top temperature not
!infrared brightness temperature at this point in the code
tb(j,ibox)= 1307.27_r8/ (log(1._r8 + (1._r8/fluxtop(j,ibox))))
if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then
tb(j,ibox) = attrop(j) - 5._r8
tau(j,ibox) = 2.13_r8*taumin(j)
end if
else
!clear sky brightness temperature
tb(j,ibox) = meantbclr(j)
end if
enddo ! j
enddo ! ibox
if (ncolprint.ne.0) then
do j=1,npoints,1000
write(iulog,'(a10)') 'j='
write(iulog,'(8I10)') j
write(iulog,'(a)') 'attrop:'
write(iulog,'(8f7.2)') (attrop(j))
write(iulog,'(a)') 'btcmin:'
write(iulog,'(8f7.2)') (btcmin(j))
write(iulog,'(a)') 'fluxtop_clrsky*100:'
write(iulog,'(8f7.2)') &
(100._r8*fluxtop_clrsky(j))
write(iulog,'(a)') '100.*f_adj:'
write(iulog,'(8f7.2)') (100._r8*fluxtop(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'transmax:'
write(iulog,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'tau:'
write(iulog,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'emcld:'
write(iulog,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'total_trans:'
write(iulog,'(8f7.2)') &
(trans_layers_above(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'total_emiss:'
write(iulog,'(8f7.2)') &
(1.0_r8-trans_layers_above(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'total_trans:'
write(iulog,'(8f7.2)') &
(trans_layers_above(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'ppout:'
write(iulog,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
enddo ! j
endif
end if
! ---------------------------------------------------!
!
! ---------------------------------------------------!
! DETERMINE CLOUD TOP PRESSURE
!
! again the 2 methods differ according to whether
! or not you use the physical cloud top pressure (top_height = 2)
! or the radiatively determined cloud top pressure (top_height = 1 or 3)
!
!compute cloud top pressure
do 30 ibox=1,ncol
!segregate according to optical thickness
if (top_height .eq. 1 .or. top_height .eq. 3) then
!find level whose temperature
!most closely matches brightness temperature
do j=1,npoints
nmatch(j)=0
enddo
do 29 k1=1,nlev-1
if (top_height_direction .eq. 2) then
ilev = nlev - k1
else
ilev = k1
end if
!cdir nodep
do j=1,npoints
if (ilev .ge. itrop(j)) then
if ((at(j,ilev) .ge. tb(j,ibox) .and. &
at(j,ilev+1) .le. tb(j,ibox)) .or. &
(at(j,ilev) .le. tb(j,ibox) .and. &
at(j,ilev+1) .ge. tb(j,ibox))) then
nmatch(j)=nmatch(j)+1
match(j,nmatch(j))=ilev
end if
end if
enddo
29 continue
do j=1,npoints
if (nmatch(j) .ge. 1) then
k1 = match(j,nmatch(j))
k2 = k1 + 1
logp1 = log(pfull(j,k1))
logp2 = log(pfull(j,k2))
atd = max(tauchk,abs(at(j,k2) - at(j,k1)))
logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd
ptop(j,ibox) = exp(logp)
ttop(j,ibox)=tb(j,ibox)
if(abs(pfull(j,k1)-ptop(j,ibox)) .lt. &
abs(pfull(j,k2)-ptop(j,ibox))) then
levmatch(j,ibox)=k1
else
levmatch(j,ibox)=k2
end if
else
if (tb(j,ibox) .le. attrop(j)) then
ptop(j,ibox)=ptrop(j)
ttop(j,ibox)=attrop(j)
levmatch(j,ibox)=itrop(j)
end if
if (tb(j,ibox) .ge. atmax(j)) then
ptop(j,ibox)=pfull(j,nlev)
ttop(j,ibox)=atmax(j)
levmatch(j,ibox)=nlev
end if
!++ bee dbg
! if ((tb(j,ibox) .gt. attrop(j)) .and. (tb(j,ibox) .lt. atmax(j)) ) then
! write(iulog,*)'isccp_cloud_types: warning: nmatch=0 and attrop<tb<atmax'
! end if
!-- bee
end if
enddo ! j
else ! if (top_height .eq. 1 .or. top_height .eq. 3)
do j=1,npoints
ptop(j,ibox)=0._r8
ttop(j,ibox)=0._r8
enddo
do ilev=1,nlev
do j=1,npoints
if ((ptop(j,ibox) .eq. 0._r8 ) &
.and.(frac_out(j,ibox,ilev) .ne. 0)) then
ptop(j,ibox)=phalf(j,ilev)
levmatch(j,ibox)=ilev
end if
end do
end do
end if
do j=1,npoints
if (tau(j,ibox) .le. (tauchk )) then
ptop(j,ibox)=0._r8
levmatch(j,ibox)=0
endif
enddo
30 continue
!
!
! ---------------------------------------------------!
!
! ---------------------------------------------------!
! DETERMINE ISCCP CLOUD TYPE FREQUENCIES
!
! Now that ptop and tau have been determined,
! determine amount of each of the 49 ISCCP cloud
! types
!
! Also compute grid box mean cloud top pressure and
! optical thickness. The mean cloud top pressure and
! optical thickness are averages over the cloudy
! area only. The mean cloud top pressure is a linear
! average of the cloud top pressures. The mean cloud
! optical thickness is computed by converting optical
! thickness to an albedo, averaging in albedo units,
! then converting the average albedo back to a mean
! optical thickness.
!
!compute isccp frequencies
!reset frequencies
do 38 ilev=1,7
do 38 ilev2=1,7
do j=1,npoints !
if (sunlit(j).eq.1 .or. top_height .eq. 3) then
fq_isccp(j,ilev,ilev2)= 0.
else
fq_isccp(j,ilev,ilev2)= output_missing_value
end if
enddo
38 continue
!reset variables need for averaging cloud properties
do j=1,npoints
totalcldarea(j) = 0._r8
meanalbedocld(j) = 0._r8
meanptop(j) = 0._r8
meanttop(j) = 0._r8
meantaucld(j) = 0._r8
enddo ! j
boxarea = 1._r8/real(ncol,r8)
do 39 ibox=1,ncol
do j=1,npoints
if (tau(j,ibox) .gt. (tauchk ) &
.and. ptop(j,ibox) .gt. 0._r8) then
box_cloudy(j,ibox)=.true.
endif
if (box_cloudy(j,ibox)) then
if (sunlit(j).eq.1 .or. top_height .eq. 3) then
boxtau(j,ibox) = tau(j,ibox)
if (tau(j,ibox) .ge. isccp_taumin) then
totalcldarea(j) = totalcldarea(j) + boxarea
!convert optical thickness to albedo
albedocld(j,ibox) &
=(tau(j,ibox)**0.895_r8)/((tau(j,ibox)**0.895_r8)+6.82_r8)
!contribute to averaging
meanalbedocld(j) = meanalbedocld(j) &
+albedocld(j,ibox)*boxarea
end if
endif
endif
if (sunlit(j).eq.1 .or. top_height .eq. 3) then
!convert ptop to millibars
ptop(j,ibox)=ptop(j,ibox) / 100._r8
!save for output cloud top pressure and optical thickness
boxptop(j,ibox) = ptop(j,ibox)
if (box_cloudy(j,ibox)) then
if (tau(j,ibox) .ge. isccp_taumin) then
meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea
meanttop(j) = meanttop(j) + ttop(j,ibox)*boxarea
end if
!reset itau(j), ipres(j)
itau(j) = 0
ipres(j) = 0
!determine optical depth category
if (tau(j,ibox) .lt. isccp_taumin) then
itau(j)=1
else if (tau(j,ibox) .ge. isccp_taumin &
.and. tau(j,ibox) .lt. 1.3_r8) then
itau(j)=2
else if (tau(j,ibox) .ge. 1.3_r8 &
.and. tau(j,ibox) .lt. 3.6_r8) then
itau(j)=3
else if (tau(j,ibox) .ge. 3.6_r8 &
.and. tau(j,ibox) .lt. 9.4_r8) then
itau(j)=4
else if (tau(j,ibox) .ge. 9.4_r8 &
.and. tau(j,ibox) .lt. 23._r8) then
itau(j)=5
else if (tau(j,ibox) .ge. 23._r8 &
.and. tau(j,ibox) .lt. 60._r8) then
itau(j)=6
else if (tau(j,ibox) .ge. 60._r8) then
itau(j)=7
end if
!determine cloud top pressure category
if ( ptop(j,ibox) .gt. 0._r8 &
.and.ptop(j,ibox) .lt. 180._r8) then
ipres(j)=1
else if(ptop(j,ibox) .ge. 180._r8 &
.and.ptop(j,ibox) .lt. 310._r8) then
ipres(j)=2
else if(ptop(j,ibox) .ge. 310._r8 &
.and.ptop(j,ibox) .lt. 440._r8) then
ipres(j)=3
else if(ptop(j,ibox) .ge. 440._r8 &
.and.ptop(j,ibox) .lt. 560._r8) then
ipres(j)=4
else if(ptop(j,ibox) .ge. 560._r8 &
.and.ptop(j,ibox) .lt. 680._r8) then
ipres(j)=5
else if(ptop(j,ibox) .ge. 680._r8 &
.and.ptop(j,ibox) .lt. 800._r8) then
ipres(j)=6
else if(ptop(j,ibox) .ge. 800._r8) then
ipres(j)=7
end if
!update frequencies
if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then
fq_isccp(j,itau(j),ipres(j))= &
fq_isccp(j,itau(j),ipres(j))+ boxarea
end if
end if
end if
enddo ! j
39 continue
!compute mean cloud properties
do j=1,npoints
if (totalcldarea(j) .gt. 0._r8) then
! code above guarantees that totalcldarea > 0
! only if sunlit .eq. 1 .or. top_height = 3
! and applies only to clouds with tau > isccp_taumin
meanptop(j) = meanptop(j) / totalcldarea(j)
meanttop(j) = meanttop(j) / totalcldarea(j)
meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j)
meantaucld(j) = &
(6.82_r8/((1./meanalbedocld(j))-1._r8))**(1._r8/0.895_r8)
else
! this code is necessary so that in the case that totalcldarea = 0.,
! that these variables, which are in-cloud averages, are set to missing
! note that totalcldarea will be 0. if all the clouds in the grid box have
! tau < isccp_taumin
meanptop(j) = output_missing_value
meanttop(j) = output_missing_value
meanalbedocld(j) = output_missing_value
meantaucld(j) = output_missing_value
end if
enddo ! j
!
! ---------------------------------------------------!
! ---------------------------------------------------!
! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM
!
if (debugcol.ne.0) then
!
do j=1,npoints,debugcol
!produce character output
do ilev=1,nlev
do ibox=1,ncol
acc(ilev,ibox)=0
enddo
enddo
do ilev=1,nlev
do ibox=1,ncol
acc(ilev,ibox)=frac_out(j,ibox,ilev)*2
if (levmatch(j,ibox) .eq. ilev) &
acc(ilev,ibox)=acc(ilev,ibox)+1
enddo
enddo
!print test
write(ftn09,11) j
11 format('ftn09.',i4.4)
unitn = getunit()
open(unitn, FILE=ftn09, FORM='FORMATTED')
write(unitn,'(a1)') ' '
write(unitn,'(10i5)') (ilev,ilev=5,nlev,5)
write(unitn,'(a1)') ' '
do ibox=1,ncol
write(unitn,'(40(a1),1x,40(a1))') &
(cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev) &
,(cchar(acc(ilev,ibox)+1),ilev=1,nlev)
end do
close( unitn )
call freeunit( unitn )
if (ncolprint.ne.0) then
write(iulog,'(a1)') ' '
write(iulog,'(a2,1X,5(a7,1X),a50)') &
'ilev', &
'pfull','at', &
'cc*100','dem_s','dtau_s', &
'cchar'
! do 4012 ilev=1,nlev
! write(iulog,'(60i2)') (box(i,ilev),i=1,ncolprint)
! write(iulog,'(i2,1X,5(f7.2,1X),50(a1))') &
! ilev, &
! pfull(j,ilev)/100.,at(j,ilev), &
! cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev) &
! ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint)
!4012 continue
write(iulog,'(a)') 'skt(j):'
write(iulog,'(8f7.2)') skt(j)
write(iulog,'(8I7)') (ibox,ibox=1,ncolprint)
write(iulog,'(a)') 'tau:'
write(iulog,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'tb:'
write(iulog,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint)
write(iulog,'(a)') 'ptop:'
write(iulog,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint)
endif
enddo
end if
return
end subroutine ISCCP_CLOUD_TYPES
end module icarus_scops_38