| 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 | ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 32 | MODULE 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 | |
|---|
| 44 | contains |
|---|
| 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 | |
|---|
| 292 | end MODULE MOD_MISR_SIMULATOR |
|---|
