source: LMDZ6/trunk/libf/phylmd/cosp2/MISR_simulator.F90 @ 3358

Last change on this file since 3358 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"

File size: 13.0 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2009, Roger Marchand, version 1.2
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. Redistributions in binary form must reproduce the above copyright notice, this list
12!    of conditions and the following disclaimer in the documentation and/or other
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History
30! May 2015 - D. Swales - Modified for COSPv2.0
31! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32MODULE MOD_MISR_SIMULATOR
33  use cosp_kinds,      only: wp
34  use MOD_COSP_STATS,  ONLY: hist2D
35  use mod_cosp_config, ONLY: R_UNDEF,numMISRHgtBins,numMISRTauBins,misr_histHgt,  &
36                             misr_histTau
37  implicit none
38
39  ! Parameters
40  real(wp),parameter :: &
41       misr_taumin = 0.3_wp,            & ! Minimum optical depth for joint-histogram
42       tauchk      = -1.*log(0.9999999)   ! Lower limit on optical depth
43
44contains
45
46  ! ######################################################################################
47  ! SUBROUTINE misr_subcolumn
48  ! ######################################################################################
49  SUBROUTINE MISR_SUBCOLUMN(npoints,ncol,nlev,dtau,zfull,at,sunlit,tauOUT,               &
50                        dist_model_layertops,box_MISR_ztop)
51    ! INPUTS
52    INTEGER, intent(in) :: &
53         npoints,    & ! Number of horizontal gridpoints
54         ncol,       & ! Number of subcolumns
55         nlev          ! Number of vertical layers
56    INTEGER, intent(in),dimension(npoints) :: &
57         sunlit        ! 1 for day points, 0 for night time
58    REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: &
59         dtau          ! Optical thickness
60    REAL(WP),intent(in),dimension(npoints,nlev) :: &
61         zfull,      & ! Height of full model levels (i.e. midpoints), [nlev] is bottom
62         at            ! Temperature (K)
63
64    ! OUTPUTS
65    REAL(WP),intent(out),dimension(npoints,ncol) :: &
66         box_MISR_ztop,     & ! Cloud-top height in each column
67         tauOUT               ! Optical depth in each column
68    REAL(WP),intent(out),dimension(npoints,numMISRHgtBins) :: &
69         dist_model_layertops !
70
71    ! INTERNAL VARIABLES
72    INTEGER :: ilev,j,loop,ibox,thres_crossed_MISR
73    INTEGER :: iMISR_ztop
74    REAL(WP) :: cloud_dtau,MISR_penetration_height,ztest
75
76    ! ############################################################################
77    ! Initialize
78    box_MISR_ztop(1:npoints,1:ncol) = 0._wp 
79
80    do j=1,npoints
81
82       ! Estimate distribution of Model layer tops
83       dist_model_layertops(j,:)=0
84       do ilev=1,nlev
85          ! Define location of "layer top"
86          if(ilev.eq.1 .or. ilev.eq.nlev) then
87             ztest=zfull(j,ilev)
88          else
89             ztest=0.5_wp*(zfull(j,ilev)+zfull(j,ilev-1))
90          endif
91
92          ! Find MISR layer that contains this level
93          ! *NOTE* the first MISR level is "no height" level
94          iMISR_ztop=2
95          do loop=2,numMISRHgtBins
96             if ( ztest .gt. 1000*misr_histHgt(loop+1) ) then
97                iMISR_ztop=loop+1
98             endif
99          enddo
100         
101          dist_model_layertops(j,iMISR_ztop) = dist_model_layertops(j,iMISR_ztop)+1
102       enddo
103
104       ! For each GCM cell or horizontal model grid point   
105       do ibox=1,ncol
106          ! Compute optical depth as a cummulative distribution in the vertical (nlev).
107          tauOUT(j,ibox)=sum(dtau(j,ibox,1:nlev))
108
109          thres_crossed_MISR=0
110          do ilev=1,nlev
111             ! If there a cloud, start the counter and store this height
112             if(thres_crossed_MISR .eq. 0 .and. dtau(j,ibox,ilev) .gt. 0.) then
113                ! First encountered a "cloud"
114                thres_crossed_MISR = 1 
115                cloud_dtau         = 0           
116             endif
117
118             if( thres_crossed_MISR .lt. 99 .and. thres_crossed_MISR .gt. 0 ) then
119                if( dtau(j,ibox,ilev) .eq. 0.) then
120                   ! We have come to the end of the current cloud layer without yet
121                   ! selecting a CTH boundary. Restart cloud tau counter
122                   cloud_dtau=0
123                else
124                   ! Add current optical depth to count for the current cloud layer
125                   cloud_dtau=cloud_dtau+dtau(j,ibox,ilev)
126                endif
127               
128                ! If the cloud is continuous but optically thin (< 1) from above the
129                ! current layer cloud top to the current level then MISR will like
130                ! see a top below the top of the current layer.
131                if( dtau(j,ibox,ilev).gt.0 .and. (cloud_dtau-dtau(j,ibox,ilev)) .lt. 1) then
132                   if(dtau(j,ibox,ilev) .lt. 1 .or. ilev.eq.1 .or. ilev.eq.nlev) then
133                      ! MISR will likely penetrate to some point within this layer ... the middle
134                      MISR_penetration_height=zfull(j,ilev)
135                   else
136                      ! Take the OD = 1.0 level into this layer
137                      MISR_penetration_height=0.5_wp*(zfull(j,ilev)+zfull(j,ilev-1)) - &
138                           0.5_wp*(zfull(j,ilev-1)-zfull(j,ilev+1))/dtau(j,ibox,ilev)
139                   endif
140                   box_MISR_ztop(j,ibox)=MISR_penetration_height
141                endif
142               
143                ! Check for a distinctive water layer
144                if(dtau(j,ibox,ilev) .gt. 1 .and. at(j,ilev) .gt. 273 ) then
145                   ! Must be a water cloud, take this as CTH level
146                   thres_crossed_MISR=99
147                endif
148               
149                ! If the total column optical depth is "large" than MISR can't see
150                ! anything else. Set current point as CTH level
151                if(sum(dtau(j,ibox,1:ilev)) .gt. 5) then
152                   thres_crossed_MISR=99           
153                endif
154             endif
155          enddo 
156         
157          ! Check to see if there was a cloud for which we didn't
158          ! set a MISR cloud top boundary
159          if( thres_crossed_MISR .eq. 1) then
160             ! If the cloud has a total optical depth of greater
161             ! than ~ 0.5 MISR will still likely pick up this cloud
162             ! with a height near the true cloud top
163             ! otherwise there should be no CTH
164             if(sum(dtau(j,ibox,1:nlev)) .gt. 0.5) then
165                ! keep MISR detected CTH
166             elseif(sum(dtau(j,ibox,1:nlev)) .gt. 0.2) then
167                ! MISR may detect but wont likley have a good height
168                box_MISR_ztop(j,ibox)=-1
169             else
170                ! MISR not likely to even detect.
171                ! so set as not cloudy
172                box_MISR_ztop(j,ibox)=0
173             endif
174          endif
175       enddo  ! loop of subcolumns
176       
177    enddo    ! loop of gridpoints
178   
179    ! Modify MISR CTH for satellite spatial / pattern matcher effects
180    ! Code in this region added by roj 5/2006 to account
181    ! for spatial effect of the MISR pattern matcher.
182    ! Basically, if a column is found between two neighbors
183    ! at the same CTH, and that column has no hieght or
184    ! a lower CTH, THEN misr will tend to but place the
185    ! odd column at the same height as it neighbors.
186   
187    ! This setup assumes the columns represent a about a 1 to 4 km scale
188    ! it will need to be modified significantly, otherwise
189!       ! DS2015: Add loop over gridpoints and index accordingly.
190!    if(ncol.eq.1) then
191!       ! Adjust based on neightboring points.
192!       do j=2,npoints-1   
193!          if(box_MISR_ztop(j-1,1) .gt. 0                             .and. &
194!             box_MISR_ztop(j+1,1) .gt. 0                             .and. &
195!             abs(box_MISR_ztop(j-1,1)-box_MISR_ztop(j+1,1)) .lt. 500 .and. &
196!             box_MISR_ztop(j,1) .lt. box_MISR_ztop(j+1,1)) then
197!             box_MISR_ztop(j,1) = box_MISR_ztop(j+1,1)   
198!          endif
199!       enddo
200!    else
201!       ! Adjust based on neighboring subcolumns.
202!       do j=1,npoints
203!          do ibox=2,ncol-1 
204!                 if(box_MISR_ztop(j,ibox-1) .gt. 0                                .and. &
205!                 box_MISR_ztop(j,ibox+1) .gt. 0                                .and. &
206!                 abs(box_MISR_ztop(j,ibox-1)-box_MISR_ztop(j,ibox+1)) .lt. 500 .and. &
207!                 box_MISR_ztop(j,ibox) .lt. box_MISR_ztop(j,ibox+1)) then
208!                 box_MISR_ztop(j,ibox) = box_MISR_ztop(j,ibox+1)   
209!               endif
210!          enddo
211!       enddo
212!    endif
213!    ! DS2015 END
214     
215    ! Fill dark scenes
216    do j=1,numMISRHgtBins
217       where(sunlit .ne. 1) dist_model_layertops(1:npoints,j) = R_UNDEF
218    enddo
219
220  end SUBROUTINE MISR_SUBCOLUMN
221
222  ! ######################################################################################
223  ! SUBROUTINE misr_column
224  ! ######################################################################################
225  SUBROUTINE MISR_COLUMN(npoints,ncol,box_MISR_ztop,sunlit,tau,MISR_cldarea,MISR_mean_ztop,fq_MISR_TAU_v_CTH)
226
227    ! INPUTS
228    INTEGER, intent(in) :: &
229         npoints,        & ! Number of horizontal gridpoints
230         ncol              ! Number of subcolumns
231    INTEGER, intent(in),dimension(npoints) :: &
232         sunlit            ! 1 for day points, 0 for night time
233    REAL(WP),intent(in),dimension(npoints,ncol) :: &
234         box_MISR_ztop,  & ! Cloud-top height in each column
235         tau               ! Column optical thickness
236
237    ! OUTPUTS
238    REAL(WP),intent(inout),dimension(npoints) :: &
239         MISR_cldarea,   & ! Fraction area covered by clouds
240         MISR_mean_ztop    ! Mean cloud top height MISR would observe
241    REAL(WP),intent(inout),dimension(npoints,7,numMISRHgtBins) :: &
242         fq_MISR_TAU_v_CTH ! Joint histogram of cloud-cover and tau
243
244    ! INTERNAL VARIABLES
245    INTEGER :: j
246    LOGICAL,dimension(ncol) :: box_cloudy
247    real(wp),dimension(npoints,ncol) :: tauWRK,box_MISR_ztopWRK
248    ! ############################################################################
249
250    ! Compute column quantities and joint-histogram
251    MISR_cldarea(1:npoints)                       = 0._wp
252    MISR_mean_ztop(1:npoints)                     = 0._wp
253    fq_MISR_TAU_v_CTH(1:npoints,1:7,1:numMISRHgtBins) = 0._wp
254    tauWRK(1:npoints,1:ncol)                      = tau(1:npoints,1:ncol)
255    box_MISR_ztopWRK(1:npoints,1:ncol)            = box_MISR_ztop(1:npoints,1:ncol)
256    do j=1,npoints
257
258       ! Subcolumns that are cloudy(true) and not(false)
259       box_cloudy(1:ncol) = merge(.true.,.false.,tau(j,1:ncol) .gt. tauchk)
260
261       ! Fill optically thin clouds with fill value
262       where(.not. box_cloudy(1:ncol)) tauWRK(j,1:ncol)  = -999._wp
263       where(box_MISR_ztopWRK(j,1:ncol) .eq. 0) box_MISR_ztopWRK(j,1:ncol)=-999._wp
264
265       ! Compute joint histogram and column quantities for points that are sunlit and cloudy
266       if (sunlit(j) .eq. 1) then
267          ! Joint histogram
268          call hist2D(tauWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol),ncol,misr_histTau,numMISRTauBins,&
269               1000*misr_histHgt,numMISRHgtBins,fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins))
270          fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins) =                       &
271             100._wp*fq_MISR_TAU_v_CTH(j,1:numMISRTauBins,1:numMISRHgtBins)/ncol
272
273          ! Column cloud area
274          MISR_cldarea(j)=real(count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.))/ncol
275
276          ! Column cloud-top height
277          if ( count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.) .ne. 0 ) then
278             MISR_mean_ztop(j) = sum(box_MISR_ztopWRK(j,1:ncol),box_MISR_ztopWRK(j,1:ncol) .ne. -999.)/ &
279                  count(box_MISR_ztopWRK(j,1:ncol) .ne. -999.)
280          else
281             MISR_mean_ztop(j) = R_UNDEF
282          endif
283
284       else
285          MISR_cldarea(j)         = R_UNDEF
286          MISR_mean_ztop(npoints) = R_UNDEF
287       endif
288    enddo
289
290  end SUBROUTINE MISR_COLUMN
291
292end MODULE MOD_MISR_SIMULATOR
Note: See TracBrowser for help on using the repository browser.