#include "cosp_defs.H"
#ifdef COSP_GFDL

!---------------------------------------------------------------------
!------------ FMS version number and tagname for this file -----------

! $Id: MISR_simulator.F,v 20.0 2013/12/13 23:15:55 fms Exp $
! $Name: tikal $
! cosp_version = 1.3.2

#endif
! 
! 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:
! 
!     * 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.
!

#ifdef COSP_GFDL
      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,
     &     dtau_col 
     & )
#else
      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
     & )
#endif
    

      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
#ifdef COSP_GFDL
      REAL,optional ::          dtau_col(npoints,ncol,nlev)
                               ! tau values obtained from model
                               ! stochastic columns


#endif
                                 
!     ------
!     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,
     c                    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
     
#ifdef COSP_GFDL
           if (present(dtau_col)) then
                 dtau = dtau_col(j,ibox,ilev)
              else
#endif    
             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 
                 
#ifdef COSP_GFDL
              endif
#endif
             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 ....
#ifdef COSP_GFDL
!RSH ADDS j loop and uses j as index rather than just j=1:
          do j=1,npoints
#endif
         do ibox=2,ncol-1
            
#ifndef COSP_GFDL
            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)    
#else

            if(box_MISR_ztop(j,ibox-1).gt.0 .and. 
     &             box_MISR_ztop(j,ibox+1).gt.0        ) then

                if( 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
                endif

            endif
         enddo
#ifdef COSP_GFDL
         enddo
#endif
      
        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 
