! ! Copyright (c) 2009, Roger Marchand, version 1.2 ! All rights reserved. ! $Revision: 23 $, $Date: 2011-03-31 15:41:37 +0200 (jeu. 31 mars 2011) $ ! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/MISR_simulator/MISR_simulator.f $ ! ! 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 University of Washington 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. ! SUBROUTINE MISR_simulator( & npoints, & nlev, & ncol, & sunlit, & zfull, & at, & dtau_s, & dtau_c, & frac_out, & missing_value, & fq_MISR_TAU_v_CTH, & dist_model_layertops, & MISR_mean_ztop, & MISR_cldarea & ) implicit none integer :: n_MISR_CTH parameter(n_MISR_CTH=16) ! ----- ! Input ! ----- INTEGER :: npoints ! if ncol ==1, the number of model points in the horizontal grid ! ! else the number of GCM grid points INTEGER :: nlev ! number of model vertical levels INTEGER :: ncol ! number of model sub columns ! ! (must already be generated in via scops and passed to this ! ! routine via the variable frac_out ) INTEGER :: sunlit(npoints) ! 1 for day points, 0 for night time REAL :: zfull(npoints,nlev) ! height (in meters) of full model levels (i.e. midpoints) ! ! zfull(npoints,1) is top level of model ! ! zfull(npoints,nlev) is bottom level of model (closest point to surface) REAL :: at(npoints,nlev) ! temperature in each model level (K) REAL :: dtau_s(npoints,nlev) ! visible wavelength cloud optical depth ... for "stratiform" condensate ! ! NOTE: this the cloud optical depth of only the ! ! the model cell (i,j) REAL :: dtau_c(npoints,nlev) ! visible wavelength cloud optical depth ... for "convective" condensate ! ! NOTE: this the cloud optical depth of only the ! ! the model cell (i,j) REAL :: frac_out(npoints,ncol,nlev) ! NOTE: only need if columns>1 ... subgrid scheme in use. REAL :: missing_value ! ------ ! Outputs ! ------ REAL :: fq_MISR_TAU_v_CTH(npoints,7,n_MISR_CTH) REAL :: dist_model_layertops(npoints,n_MISR_CTH) REAL :: MISR_cldarea(npoints) ! fractional area coverged by clouds REAL :: MISR_mean_ztop(npoints) ! mean cloud top hieght(m) MISR would observe ! ! NOTE: == 0 if area ==0 ! ------ ! Working variables ! ------ REAL :: tau(npoints,ncol) ! total column optical depth ... INTEGER :: j,ilev,ilev2,ibox,k INTEGER :: itau LOGICAL :: box_cloudy(npoints,ncol) real :: isccp_taumin real :: boxarea real :: tauchk REAL :: box_MISR_ztop(npoints,ncol) ! cloud top hieght(m) MISR would observe integer :: thres_crossed_MISR integer :: loop,iMISR_ztop real :: dtau, cloud_dtau, MISR_penetration_height,ztest real :: MISR_CTH_boundaries(n_MISR_CTH+1) DATA MISR_CTH_boundaries / -99, 0, 0.5, 1, 1.5, 2, 2.5, 3, & 4, 5, 7, 9, 11, 13, 15, 17, 99 / DATA isccp_taumin / 0.3 / tauchk = -1.*log(0.9999999) ! ! ! ! For each GCM cell or horizontal model grid point ... ! ! do j=1,npoints ! ! ! ! estimate distribution of Model layer tops ! ! dist_model_layertops(j,:)=0 do ilev=1,nlev ! ! define location of "layer top" if(ilev.eq.1 .or. ilev.eq.nlev) then ztest=zfull(j,ilev) else ztest=0.5*(zfull(j,ilev)+zfull(j,ilev-1)) endif ! ! find MISR layer that contains this level ! ! note, the first MISR level is "no height" level iMISR_ztop=2 do loop=2,n_MISR_CTH if ( ztest .gt. & 1000*MISR_CTH_boundaries(loop+1) ) then iMISR_ztop=loop+1 endif enddo dist_model_layertops(j,iMISR_ztop)= & dist_model_layertops(j,iMISR_ztop)+1 enddo ! ! ! ! compute total cloud optical depth for each column ! ! do ibox=1,ncol ! ! Initialize tau to zero in each subcolum tau(j,ibox)=0. box_cloudy(j,ibox)=.false. box_MISR_ztop(j,ibox)=0 ! ! initialize threshold detection for each sub column thres_crossed_MISR=0; do ilev=1,nlev dtau=0 if (frac_out(j,ibox,ilev).eq.1) then dtau = dtau_s(j,ilev) endif if (frac_out(j,ibox,ilev).eq.2) then dtau = dtau_c(j,ilev) end if tau(j,ibox)=tau(j,ibox)+ dtau ! ! NOW for MISR .. ! ! if there a cloud ... start the counter ... store this height if(thres_crossed_MISR .eq. 0 .and. dtau .gt. 0.) then ! ! first encountered a "cloud" thres_crossed_MISR=1 cloud_dtau=0 endif if( thres_crossed_MISR .lt. 99 .and. & thres_crossed_MISR .gt. 0 ) then if( dtau .eq. 0.) then ! ! we have come to the end of the current cloud ! ! layer without yet selecting a CTH boundary. ! ! ... restart cloud tau counter cloud_dtau=0 else ! ! add current optical depth to count for ! ! the current cloud layer cloud_dtau=cloud_dtau+dtau endif ! ! if the cloud is continuous but optically thin (< 1) ! ! from above the current layer cloud top to the current level ! ! then MISR will like see a top below the top of the current ! ! layer if( dtau.gt.0 .and. (cloud_dtau-dtau) .lt. 1) then if(dtau .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then ! ! MISR will likely penetrate to some point ! ! within this layer ... the middle MISR_penetration_height=zfull(j,ilev) else ! ! take the OD = 1.0 level into this layer MISR_penetration_height= & 0.5*(zfull(j,ilev)+zfull(j,ilev-1)) - & 0.5*(zfull(j,ilev-1)-zfull(j,ilev+1)) & /dtau endif box_MISR_ztop(j,ibox)=MISR_penetration_height endif ! ! check for a distinctive water layer if(dtau .gt. 1 .and. at(j,ilev).gt.273 ) then ! ! must be a water cloud ... ! ! take this as CTH level thres_crossed_MISR=99 endif ! ! if the total column optical depth is "large" than ! ! MISR can't seen anything else ... set current point as CTH level if(tau(j,ibox) .gt. 5) then thres_crossed_MISR=99 endif endif ! MISR CTH booundary not set enddo !ilev - loop over vertical levesl ! ! written by roj 5/2006 ! ! check to see if there was a cloud for which we didn't ! ! set a MISR cloud top boundary if( thres_crossed_MISR .eq. 1) then ! ! if the cloud has a total optical depth of greater ! ! than ~ 0.5 MISR will still likely pick up this cloud ! ! with a height near the true cloud top ! ! otherwise there should be no CTH if( tau(j,ibox) .gt. 0.5) then ! ! keep MISR detected CTH elseif(tau(j,ibox) .gt. 0.2) then ! ! MISR may detect but wont likley have a good height box_MISR_ztop(j,ibox)=-1 else ! ! MISR not likely to even detect. ! ! so set as not cloudy box_MISR_ztop(j,ibox)=0 endif endif enddo ! loop of subcolumns enddo ! loop of gridpoints ! ! ! ! Modify MISR CTH for satellite spatial / pattern matcher effects ! ! ! ! Code in this region added by roj 5/2006 to account ! ! for spatial effect of the MISR pattern matcher. ! ! Basically, if a column is found between two neighbors ! ! at the same CTH, and that column has no hieght or ! ! a lower CTH, THEN misr will tend to but place the ! ! odd column at the same height as it neighbors. ! ! ! ! This setup assumes the columns represent a about a 1 to 4 km scale ! ! it will need to be modified significantly, otherwise if(ncol.eq.1) then ! ! adjust based on neightboring points ... i.e. only 2D grid was input do j=2,npoints-1 if(box_MISR_ztop(j-1,1).gt.0 .and. & box_MISR_ztop(j+1,1).gt.0 ) then if( abs( box_MISR_ztop(j-1,1) - & box_MISR_ztop(j+1,1) ) .lt. 500 & .and. & box_MISR_ztop(j,1) .lt. & box_MISR_ztop(j+1,1) ) then box_MISR_ztop(j,1) = & box_MISR_ztop(j+1,1) endif endif enddo else ! ! adjust based on neighboring subcolumns .... do ibox=2,ncol-1 if(box_MISR_ztop(1,ibox-1).gt.0 .and. & box_MISR_ztop(1,ibox+1).gt.0 ) then if( abs( box_MISR_ztop(1,ibox-1) - & box_MISR_ztop(1,ibox+1) ) .lt. 500 & .and. & box_MISR_ztop(1,ibox) .lt. & box_MISR_ztop(1,ibox+1) ) then box_MISR_ztop(1,ibox) = & box_MISR_ztop(1,ibox+1) endif endif enddo endif ! ! ! ! DETERMINE CLOUD TYPE FREQUENCIES ! ! ! ! Now that ztop and tau have been determined, ! ! determine amount of each cloud type boxarea=1./real(ncol) do j=1,npoints ! ! reset frequencies -- modified loop structure, roj 5/2006 do ilev=1,7 ! "tau loop" do ilev2=1,n_MISR_CTH fq_MISR_TAU_v_CTH(j,ilev,ilev2)=0. enddo enddo MISR_cldarea(j)=0. MISR_mean_ztop(j)=0. do ibox=1,ncol if (tau(j,ibox) .gt. (tauchk)) then box_cloudy(j,ibox)=.true. endif itau = 0 if (box_cloudy(j,ibox)) then ! !determine optical depth category if (tau(j,ibox) .lt. isccp_taumin) then itau=1 else if (tau(j,ibox) .ge. isccp_taumin & .and. tau(j,ibox) .lt. 1.3) then itau=2 else if (tau(j,ibox) .ge. 1.3 & .and. tau(j,ibox) .lt. 3.6) then itau=3 else if (tau(j,ibox) .ge. 3.6 & .and. tau(j,ibox) .lt. 9.4) then itau=4 else if (tau(j,ibox) .ge. 9.4 & .and. tau(j,ibox) .lt. 23.) then itau=5 else if (tau(j,ibox) .ge. 23. & .and. tau(j,ibox) .lt. 60.) then itau=6 else if (tau(j,ibox) .ge. 60.) then itau=7 endif endif ! ! update MISR histograms and summary metrics - roj 5/2005 if (sunlit(j).eq.1) then ! !if cloudy added by roj 5/2005 if( box_MISR_ztop(j,ibox).eq.0) then ! ! no cloud detected iMISR_ztop=0 elseif( box_MISR_ztop(j,ibox).eq.-1) then ! ! cloud can be detected but too thin to get CTH iMISR_ztop=1 fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)= & fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea else ! ! ! ! determine index for MISR bin set ! ! iMISR_ztop=2 do loop=2,n_MISR_CTH if ( box_MISR_ztop(j,ibox) .gt. & 1000*MISR_CTH_boundaries(loop+1) ) then iMISR_ztop=loop+1 endif enddo if(box_cloudy(j,ibox)) then ! ! there is an isccp clouds so itau(j) is defined fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop)= & fq_MISR_TAU_v_CTH(j,itau,iMISR_ztop) + boxarea else ! ! MISR CTH resolution is trying to fill in a ! ! broken cloud scene where there is no condensate. ! ! The MISR CTH-1D-OD product will only put in a cloud ! ! if the MISR cloud mask indicates cloud. ! ! therefore we will not include this column in the histogram ! ! in reality aerosoal and 3D effects or bright surfaces ! ! could fool the MISR cloud mask ! ! the alternative is to count as very thin cloud ?? ! fq_MISR_TAU_v_CTH(1,iMISR_ztop)= ! & fq_MISR_TAU_v_CTH(1,iMISR_ztop) + boxarea endif MISR_mean_ztop(j)=MISR_mean_ztop(j)+ & box_MISR_ztop(j,ibox)*boxarea MISR_cldarea(j)=MISR_cldarea(j) + boxarea endif else ! ! Set to issing data. A. Bodas - 14/05/2010 do loop=1,n_MISR_CTH do k=1,7 fq_MISR_TAU_v_CTH(j,k,loop) = missing_value enddo dist_model_layertops(j,loop) = missing_value enddo MISR_cldarea(j) = missing_value MISR_mean_ztop(npoints) = missing_value endif ! is sunlight ? enddo ! ibox - loop over subcolumns if( MISR_cldarea(j) .gt. 0.) then MISR_mean_ztop(j)= MISR_mean_ztop(j) / MISR_cldarea(j) ! roj 5/2006 endif enddo ! loop over grid points return end subroutine misr_simulator