source: LMDZ6/trunk/libf/phylmd/cospv2/MISR_simulator.F90 @ 4744

! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
! Copyright (c) 2009, Roger Marchand, version 1.2
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are 
! permitted provided that the following conditions are met:
!
! 1. Redistributions of source code must retain the above copyright notice, this list of 
!    conditions and the following disclaimer.
!
! 2. 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.
!
! 3. Neither the name of the copyright holder 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 HOLDER 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.
!
! History
! May 2015 - D. Swales - Modified for COSPv2.0
! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MODULE MOD_MISR_SIMULATOR
  use cosp_kinds,      only: wp
  use MOD_COSP_STATS,  ONLY: hist2D
  use mod_cosp_config, ONLY: R_UNDEF,numMISRHgtBins,numMISRTauBins,misr_histHgt,  &
                             misr_histTau
  implicit none 

  ! Parameters
  real(wp),parameter :: &
       misr_taumin = 0.3_wp,            & ! Minimum optical depth for joint-histogram
       tauchk      = -1.*log(0.9999999)   ! Lower limit on optical depth

contains

  ! ######################################################################################
  ! SUBROUTINE misr_subcolumn
  ! ######################################################################################
  SUBROUTINE MISR_SUBCOLUMN(npoints,ncol,nlev,dtau,zfull,at,sunlit,tauOUT,               &
                        dist_model_layertops,box_MISR_ztop)
    ! INPUTS
    INTEGER, intent(in) :: &
         npoints,    & ! Number of horizontal gridpoints
         ncol,       & ! Number of subcolumns
         nlev          ! Number of vertical layers
    INTEGER, intent(in),dimension(npoints) :: &
         sunlit        ! 1 for day points, 0 for night time
    REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: &
         dtau          ! Optical thickness
    REAL(WP),intent(in),dimension(npoints,nlev) :: &
         zfull,      & ! Height of full model levels (i.e. midpoints), [nlev] is bottom
         at            ! Temperature (K)

    ! OUTPUTS
    REAL(WP),intent(out),dimension(npoints,ncol) :: &
         box_MISR_ztop,     & ! Cloud-top height in each column
         tauOUT               ! Optical depth in each column
    REAL(WP),intent(out),dimension(npoints,numMISRHgtBins) :: &
         dist_model_layertops ! 

    ! INTERNAL VARIABLES
    INTEGER :: ilev,j,loop,ibox,thres_crossed_MISR
    INTEGER :: iMISR_ztop
    REAL(WP) :: cloud_dtau,MISR_penetration_height,ztest

    ! ############################################################################
    ! Initialize
    box_MISR_ztop(1:npoints,1:ncol) = 0._wp  

    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_wp*(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,numMISRHgtBins
             if ( ztest .gt. 1000*misr_histHgt(loop+1) ) then
                iMISR_ztop=loop+1
             endif
          enddo
          
          dist_model_layertops(j,iMISR_ztop) = dist_model_layertops(j,iMISR_ztop)+1
       enddo

       ! For each GCM cell or horizontal model grid point   
       do ibox=1,ncol
          ! Compute optical depth as a cummulative distribution in the vertical (nlev).
          tauOUT(j,ibox)=sum(dtau(j,ibox,1:nlev))

          thres_crossed_MISR=0
          do ilev=1,nlev
             ! If there a cloud, start the counter and store this height
             if(thres_crossed_MISR .eq. 0 .and. dtau(j,ibox,ilev) .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(j,ibox,ilev) .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(j,ibox,ilev)
                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(j,ibox,ilev).gt.0 .and. (cloud_dtau-dtau(j,ibox,ilev)) .lt. 1) then
                   if(dtau(j,ibox,ilev) .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_wp*(zfull(j,ilev)+zfull(j,ilev-1)) - &
                           0.5_wp*(zfull(j,ilev-1)-zfull(j,ilev+1))/dtau(j,ibox,ilev) 
                   endif
                   box_MISR_ztop(j,ibox)=MISR_penetration_height
                endif
                
                ! Check for a distinctive water layer
                if(dtau(j,ibox,ilev) .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 see
                ! anything else. Set current point as CTH level
                if(sum(dtau(j,ibox,1:ilev)) .gt. 5) then
                   thres_crossed_MISR=99           
                endif
             endif
          enddo  
          
          ! 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(sum(dtau(j,ibox,1:nlev)) .gt. 0.5) then
                ! keep MISR detected CTH
             elseif(sum(dtau(j,ibox,1:nlev)) .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
!       ! DS2015: Add loop over gridpoints and index accordingly.
!    if(ncol.eq.1) then
!       ! Adjust based on neightboring points.
!       do j=2,npoints-1   
!          if(box_MISR_ztop(j-1,1) .gt. 0                             .and. &
!             box_MISR_ztop(j+1,1) .gt. 0                             .and. &
!             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
!       enddo
!    else
!       ! Adjust based on neighboring subcolumns.
!       do j=1,npoints
!          do ibox=2,ncol-1  
!                 if(box_MISR_ztop(j,ibox-1) .gt. 0                                .and. &
!                 box_MISR_ztop(j,ibox+1) .gt. 0                                .and. &
!                 abs(box_MISR_ztop(j,ibox-1)-box_MISR_ztop(j,ibox+1)) .lt. 500 .and. &
!                 box_MISR_ztop(j,ibox) .lt. box_MISR_ztop(j,ibox+1)) then
!                 box_MISR_ztop(j,ibox) = box_MISR_ztop(j,ibox+1)    
!               endif
!          enddo
!       enddo
!    endif
!    ! DS2015 END
     
    ! Fill dark scenes 
    do j=1,numMISRHgtBins
       where(sunlit .ne. 1) dist_model_layertops(1:npoints,j) = R_UNDEF
    enddo

  end SUBROUTINE MISR_SUBCOLUMN

  ! ######################################################################################
  ! SUBROUTINE misr_column
  ! ######################################################################################
  SUBROUTINE MISR_COLUMN(npoints,ncol,box_MISR_ztop,sunlit,tau,MISR_cldarea,MISR_mean_ztop,fq_MISR_TAU_v_CTH)

    ! INPUTS
    INTEGER, intent(in) :: &
         npoints,        & ! Number of horizontal gridpoints
         ncol              ! Number of subcolumns
    INTEGER, intent(in),dimension(npoints) :: &
         sunlit            ! 1 for day points, 0 for night time
    REAL(WP),intent(in),dimension(npoints,ncol) :: &
         box_MISR_ztop,  & ! Cloud-top height in each column
         tau               ! Column optical thickness

    ! OUTPUTS
    REAL(WP),intent(inout),dimension(npoints) :: &
         MISR_cldarea,   & ! Fraction area covered by clouds
         MISR_mean_ztop    ! Mean cloud top height MISR would observe
    REAL(WP),intent(inout),dimension(npoints,7,numMISRHgtBins) :: &
         fq_MISR_TAU_v_CTH ! Joint histogram of cloud-cover and tau

    ! INTERNAL VARIABLES
    INTEGER :: j
    LOGICAL,dimension(ncol) :: box_cloudy 
    real(wp),dimension(npoints,ncol) :: tauWRK,box_MISR_ztopWRK
    ! ############################################################################

    ! Compute column quantities and joint-histogram
    MISR_cldarea(1:npoints)                       = 0._wp
    MISR_mean_ztop(1:npoints)                     = 0._wp
    fq_MISR_TAU_v_CTH(1:npoints,1:7,1:numMISRHgtBins) = 0._wp
    tauWRK(1:npoints,1:ncol)                      = tau(1:npoints,1:ncol)
    box_MISR_ztopWRK(1:npoints,1:ncol)            = box_MISR_ztop(1:npoints,1:ncol)
    do j=1,npoints

       ! Subcolumns that are cloudy(true) and not(false)
       box_cloudy(1:ncol) = merge(.true.,.false.,tau(j,1:ncol) .gt. tauchk)

       ! Fill optically thin clouds with fill value
       where(.not. box_cloudy(1:ncol)) tauWRK(j,1:ncol)  = -999._wp
       where(box_MISR_ztopWRK(j,1:ncol) .eq. 0) box_MISR_ztopWRK(j,1:ncol)=-999._wp

       ! Compute joint histogram and column quantities for points that are sunlit and cloudy
       if (sunlit(j) .eq. 1) then 
          ! Joint histogram
          call hist2D(tauWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol),ncol,misr_histTau,numMISRTauBins,&
               1000*misr_histHgt,numMISRHgtBins,fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins))
          fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins) =                       &
             100._wp*fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins)/ncol

          ! Column cloud area
          MISR_cldarea(j)=real(count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.))/ncol

          ! Column cloud-top height
          if ( count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.) .ne. 0 ) then
             MISR_mean_ztop(j) = sum(box_MISR_ztopWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol) .ne. -999.)/ &
                  count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.)
          else
             MISR_mean_ztop(j) = R_UNDEF
          endif

       else
          MISR_cldarea(j)         = R_UNDEF
          MISR_mean_ztop(npoints) = R_UNDEF
       endif
    enddo

  end SUBROUTINE MISR_COLUMN

end MODULE MOD_MISR_SIMULATOR