! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $ ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/icarus-scops-4.1-bsd/icarus.f $ SUBROUTINE ICARUS( & debug, & debugcol, & npoints, & sunlit, & nlev, & ncol, & pfull, & phalf, & qv, & cc, & conv, & dtau_s, & dtau_c, & top_height, & top_height_direction, & overlap, & frac_out, & skt, & emsfc_lw, & at, & dem_s, & dem_c, & fq_isccp, & totalcldarea, & meanptop, & meantaucld, & meanalbedocld, & meantb, & meantbclr, & boxtau, & boxptop & ) !$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $ ! *****************************COPYRIGHT**************************** ! (c) 2009, Lawrence Livermore National Security Limited Liability ! Corporation. ! All rights reserved. ! ! Redistribution and use in source and binary forms, with or without ! modification, are permitted provided that the ! following conditions are met: ! ! * Redistributions of source code must retain the above ! copyright notice, this list of conditions and the following ! disclaimer. ! * Redistributions in binary form must reproduce the above ! copyright notice, this list of conditions and the following ! disclaimer in the documentation and/or other materials ! provided with the distribution. ! * Neither the name of the Lawrence Livermore National Security ! Limited Liability Corporation nor the names of its ! contributors may be used to endorse or promote products ! derived from this software without specific prior written ! permission. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ! *****************************COPYRIGHT******************************* ! *****************************COPYRIGHT******************************* ! *****************************COPYRIGHT******************************* ! *****************************COPYRIGHT******************************* 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 REAL :: pfull(npoints,nlev) ! ! pressure of full model levels (Pascals) ! ! pfull(npoints,1) is top level of model ! ! pfull(npoints,nlev) is bot 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) 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 = old setting: matches all versions of ! ISCCP simulator with versions numbers 3.5.1 and lower ! ! 2 = default setting: for version numbers 4.0 and higher ! ! 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. REAL :: 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 ! ------ ! 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. 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 three 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 :: 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 :: meanalbedocld(npoints) ! mean cloud albedo ! ! linear averaging in albedo performed real :: meantb(npoints) ! mean all-sky 10.5 micron brightness temperature real :: meantbclr(npoints) ! mean clear-sky 10.5 micron brightness temperature 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 ! ------ REAL :: dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave ! ! emissivity in part of ! ! gridbox under consideration REAL :: ptrop(npoints) REAL :: attrop(npoints) REAL :: attropmin (npoints) REAL :: atmax(npoints) REAL :: btcmin(npoints) REAL :: 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 :: fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints) real :: taumin(npoints) real :: dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0 real :: press(npoints), dpress(npoints), atmden(npoints) real :: rvh20(npoints), wk(npoints), rhoave(npoints) real :: rh20s(npoints), rfrgn(npoints) real :: tmpexp(npoints),tauwv(npoints) character(len=1) :: cchar(6),cchar_realtops(6) integer :: icycle REAL :: tau(npoints,ncol) LOGICAL :: box_cloudy(npoints,ncol) REAL :: tb(npoints,ncol) REAL :: ptop(npoints,ncol) REAL :: emcld(npoints,ncol) REAL :: fluxtop(npoints,ncol) REAL :: trans_layers_above(npoints,ncol) real :: isccp_taumin,fluxtopinit(npoints),tauir(npoints) REAL :: albedocld(npoints,ncol) 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 :: rangevec(npoints),rangeerror integer :: index1(npoints),num1,jj,k1,k2 real :: rec2p13,tauchk,logp,logp1,logp2,atd real :: output_missing_value character(len=10) :: ftn09 DATA isccp_taumin / 0.3 / DATA output_missing_value / -1.E+30 / DATA cchar / ' ','-','1','+','I','+'/ DATA cchar_realtops / ' ',' ','1','1','I','I'/ ! ------ End duplicate definitions common to wrapper routine tauchk = -1.*log(0.9999999) rec2p13=1./2.13 ncolprint=0 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,'(a11)') 'top_height=' write(6,'(8I10)') top_height write(6,'(a21)') 'top_height_direction=' write(6,'(8I10)') top_height_direction write(6,'(a10)') 'overlap=' write(6,'(8I10)') overlap write(6,'(a10)') 'emsfc_lw=' write(6,'(8f10.2)') emsfc_lw 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)') '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.3)') (dem_c(j,i),i=1,nlev) enddo endif ! ---------------------------------------------------! if (ncolprint.ne.0) then do j=1,npoints,1000 write(6,'(a10)') 'j=' write(6,'(8I10)') j enddo endif if (top_height .eq. 1 .or. top_height .eq. 3) then do j=1,npoints ptrop(j)=5000. attropmin(j) = 400. atmax(j) = 0. attrop(j) = 120. itrop(j) = 1 enddo do ilev=1,nlev do j=1,npoints if (pfull(j,ilev) .lt. 40000. .and. & pfull(j,ilev) .gt. 5000. .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 enddo end do do ilev=1,nlev do j=1,npoints if (at(j,ilev) .gt. atmax(j) .and. & ilev .ge. itrop(j)) atmax(j)=at(j,ilev) enddo end do end if if (top_height .eq. 1 .or. top_height .eq. 3) then do j=1,npoints meantb(j) = 0. meantbclr(j) = 0. end do else do j=1,npoints meantb(j) = output_missing_value meantbclr(j) = output_missing_value end do end if ! -----------------------------------------------------! ! ---------------------------------------------------! do ilev=1,nlev do j=1,npoints rangevec(j)=0 if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then ! error = cloud fraction less than zero ! error = cloud fraction greater than 1 rangevec(j)=rangevec(j)+1 endif if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then ! ' error = convective cloud fraction less than zero' ! ' error = convective cloud fraction greater than 1' rangevec(j)=rangevec(j)+2 endif if (dtau_s(j,ilev) .lt. 0.) then ! ' error = stratiform cloud opt. depth less than zero' rangevec(j)=rangevec(j)+4 endif if (dtau_c(j,ilev) .lt. 0.) then ! ' error = convective cloud opt. depth less than zero' rangevec(j)=rangevec(j)+8 endif if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then ! ' error = stratiform cloud emissivity less than zero' ! ' error = stratiform cloud emissivity greater than 1' rangevec(j)=rangevec(j)+16 endif if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then ! ' error = convective cloud emissivity less than zero' ! ' error = convective cloud emissivity greater than 1' rangevec(j)=rangevec(j)+32 endif enddo rangeerror=0 do j=1,npoints rangeerror=rangeerror+rangevec(j) enddo if (rangeerror.ne.0) then write (6,*) 'Input variable out of range' write (6,*) 'rangevec:' write (6,*) rangevec STOP endif enddo ! ! ---------------------------------------------------! ! ! ---------------------------------------------------! ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and ! put into vector tau ! !initialize tau and albedocld to zero do ibox=1,ncol do j=1,npoints tau(j,ibox)=0. albedocld(j,ibox)=0. boxtau(j,ibox)=output_missing_value boxptop(j,ibox)=output_missing_value box_cloudy(j,ibox)=.false. enddo end do ! !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(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. 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 wtmh20 = 18.01534 Navo = 6.023E+23 grav = 9.806650E+02 pstd = 1.013250E+06 t0 = 296. if (ncolprint .ne. 0) & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv' do ilev=1,nlev do j=1,npoints ! !press and dpress are dyne/cm2 = Pascals *10 press(j) = pfull(j,ilev)*10. 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*(at(j,ilev)-t0)) tauwv(j) = wk(j)*1.e-20*( & (0.0224697*rh20s(j)*tmpexp(j)) + & (3.41817e-7*rfrgn(j)) )*0.98 dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j)) enddo 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 do j=1,npoints fluxtop_clrsky(j) = 0. trans_layers_above_clrsky(j)=1. enddo do ilev=1,nlev do j=1,npoints ! ! Black body emission at temperature of the layer bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) ! !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.-dem_wv(j,ilev)) enddo 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 enddo !loop over level do j=1,npoints ! !add in surface emission bb(j)=1/( exp(1307.27/skt(j)) - 1. ) ! !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/(log(1.+(1./fluxtop_clrsky(j)))) enddo 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,'(5(f7.2,1X))') emsfc_lw,100.*bb(j), & 100.*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(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 do ibox=1,ncol do j=1,npoints fluxtop(j,ibox)=0. trans_layers_above(j,ibox)=1. enddo enddo do ilev=1,nlev do j=1,npoints ! ! Black body emission at temperature of the layer bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) ! !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. - & ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) ) else if (frac_out(j,ibox,ilev).eq.2) then dem(j,ibox)= 1. - & ( (1. - dem_wv(j,ilev)) * (1. - 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.-dem(j,ibox)) enddo ! j enddo ! ibox 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 enddo ! ilev do j=1,npoints ! !add in surface emission bb(j)=1/( exp(1307.27/skt(j)) - 1. ) ! !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/(log(1.+(1./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(6,'(a10)') 'j=' write(6,'(8I10)') j write (6,'(a)') 'id:' write (6,'(a)') 'surface' write (6,'(a)') 'emiss_layer:' write (6,'(8f7.2)') (dem(1,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)') 'meantb(j):' write (6,'(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. / ( exp(1307.27/(attrop(j)-5.)) - 1. ) 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. * log(max(min(transmax(j),0.9999999),0.001)) enddo if (top_height .eq. 1) then do j=1,npoints if (transmax(j) .gt. 0.001 .and. & transmax(j) .le. 0.9999999) 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 .and. & transmax(j) .le. 0.9999999) then emcld(j,ibox) = 1. - exp(-1. * tauir(j) ) fluxtop(j,ibox) = fluxtopinit(j) - & ((1.-emcld(j,ibox))*fluxtop_clrsky(j)) fluxtop(j,ibox)=max(1.E-06, & (fluxtop(j,ibox)/emcld(j,ibox))) tb(j,ibox)= 1307.27 & / (log(1. + (1./fluxtop(j,ibox)))) if (tb(j,ibox) .gt. 260.) then tauir(j) = tau(j,ibox) / 2.56 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/ (log(1. + (1./fluxtop(j,ibox)))) if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then tb(j,ibox) = attrop(j) - 5. tau(j,ibox) = 2.13*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(6,'(a10)') 'j=' write(6,'(8I10)') j write (6,'(a)') 'attrop:' write (6,'(8f7.2)') (attrop(j)) write (6,'(a)') 'btcmin:' write (6,'(8f7.2)') (btcmin(j)) write (6,'(a)') 'fluxtop_clrsky*100:' write (6,'(8f7.2)') & (100.*fluxtop_clrsky(j)) write (6,'(a)') '100.*f_adj:' write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'transmax:' write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint) write (6,'(a)') 'tau:' write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'emcld:' write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'total_trans:' write (6,'(8f7.2)') & (trans_layers_above(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'total_emiss:' write (6,'(8f7.2)') & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'total_trans:' write (6,'(8f7.2)') & (trans_layers_above(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'ppout:' write (6,'(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 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 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 end do 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) 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) levmatch(j,ibox)=itrop(j) end if if (tb(j,ibox) .ge. atmax(j)) then ptop(j,ibox)=pfull(j,nlev) levmatch(j,ibox)=nlev end if end if enddo ! j else ! if (top_height .eq. 1 .or. top_height .eq. 3) do j=1,npoints ptop(j,ibox)=0. enddo do ilev=1,nlev do j=1,npoints if ((ptop(j,ibox) .eq. 0. ) & .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. levmatch(j,ibox)=0 endif enddo end do ! ! ! ---------------------------------------------------! ! ! ---------------------------------------------------! ! 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 ilev=1,7 do 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 end do end do ! !reset variables need for averaging cloud properties do j=1,npoints if (sunlit(j).eq.1 .or. top_height .eq. 3) then totalcldarea(j) = 0. meanalbedocld(j) = 0. meanptop(j) = 0. meantaucld(j) = 0. else totalcldarea(j) = output_missing_value meanalbedocld(j) = output_missing_value meanptop(j) = output_missing_value meantaucld(j) = output_missing_value end if enddo ! j boxarea = 1./real(ncol) do ibox=1,ncol do j=1,npoints if (tau(j,ibox) .gt. (tauchk ) & .and. ptop(j,ibox) .gt. 0.) 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)/((tau(j,ibox)**0.895)+6.82) ! !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 if (box_cloudy(j,ibox)) then ! !convert ptop to millibars ptop(j,ibox)=ptop(j,ibox) / 100. ! !save for output cloud top pressure and optical thickness boxptop(j,ibox) = ptop(j,ibox) if (tau(j,ibox) .ge. isccp_taumin) then meanptop(j) = meanptop(j) + ptop(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) then itau(j)=2 else if (tau(j,ibox) .ge. 1.3 & .and. tau(j,ibox) .lt. 3.6) then itau(j)=3 else if (tau(j,ibox) .ge. 3.6 & .and. tau(j,ibox) .lt. 9.4) then itau(j)=4 else if (tau(j,ibox) .ge. 9.4 & .and. tau(j,ibox) .lt. 23.) then itau(j)=5 else if (tau(j,ibox) .ge. 23. & .and. tau(j,ibox) .lt. 60.) then itau(j)=6 else if (tau(j,ibox) .ge. 60.) then itau(j)=7 end if ! !determine cloud top pressure category if ( ptop(j,ibox) .gt. 0. & .and.ptop(j,ibox) .lt. 180.) then ipres(j)=1 else if(ptop(j,ibox) .ge. 180. & .and.ptop(j,ibox) .lt. 310.) then ipres(j)=2 else if(ptop(j,ibox) .ge. 310. & .and.ptop(j,ibox) .lt. 440.) then ipres(j)=3 else if(ptop(j,ibox) .ge. 440. & .and.ptop(j,ibox) .lt. 560.) then ipres(j)=4 else if(ptop(j,ibox) .ge. 560. & .and.ptop(j,ibox) .lt. 680.) then ipres(j)=5 else if(ptop(j,ibox) .ge. 680. & .and.ptop(j,ibox) .lt. 800.) then ipres(j)=6 else if(ptop(j,ibox) .ge. 800.) 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 end do ! !compute mean cloud properties do j=1,npoints if (totalcldarea(j) .gt. 0.) 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) meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j) meantaucld(j) = (6.82/((1./meanalbedocld(j))-1.))**(1./0.895) 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 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) open(9, FILE=ftn09, FORM='FORMATTED') write(9,'(a1)') ' ' write(9,'(10i5)') & (ilev,ilev=5,nlev,5) write(9,'(a1)') ' ' do ibox=1,ncol write(9,'(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(9) if (ncolprint.ne.0) then write(6,'(a1)') ' ' write(6,'(a2,1X,5(a7,1X),a50)') & 'ilev', & 'pfull','at', & 'cc*100','dem_s','dtau_s', & 'cchar' ! do 4012 ilev=1,nlev ! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint) ! write(6,'(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 (6,'(a)') 'skt(j):' write (6,'(8f7.2)') skt(j) write (6,'(8I7)') (ibox,ibox=1,ncolprint) write (6,'(a)') 'tau:' write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'tb:' write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) write (6,'(a)') 'ptop:' write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint) endif enddo end if return end subroutine icarus