| 1 | ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 2 | ! Copyright (c) 2015, Regents of the University of Colorado |
|---|
| 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 | ! Jul 2007 - A. Bodas-Salcedo - Initial version |
|---|
| 31 | ! Jul 2008 - A. Bodas-Salcedo - Added capability of producing outputs in standard grid |
|---|
| 32 | ! Oct 2008 - J.-L. Dufresne - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns |
|---|
| 33 | ! in COSP_STATS |
|---|
| 34 | ! Oct 2008 - H. Chepfer - Added PARASOL reflectance arguments |
|---|
| 35 | ! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations |
|---|
| 36 | ! Jan 2013 - G. Cesana - Added betaperp and temperature arguments |
|---|
| 37 | ! - Added phase 3D/3Dtemperature/Map output variables in diag_lidar |
|---|
| 38 | ! May 2015 - D. Swales - Modified for cosp2.0 |
|---|
| 39 | ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 40 | MODULE MOD_COSP_STATS |
|---|
| 41 | USE COSP_KINDS, ONLY: wp |
|---|
| 42 | USE MOD_COSP_CONFIG, ONLY: R_UNDEF,R_GROUND |
|---|
| 43 | |
|---|
| 44 | IMPLICIT NONE |
|---|
| 45 | CONTAINS |
|---|
| 46 | |
|---|
| 47 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 48 | !---------- SUBROUTINE COSP_CHANGE_VERTICAL_GRID ---------------- |
|---|
| 49 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 50 | SUBROUTINE COSP_CHANGE_VERTICAL_GRID(Npoints,Ncolumns,Nlevels,zfull,zhalf,y,Nglevels,newgrid_bot,newgrid_top,r,log_units) |
|---|
| 51 | implicit none |
|---|
| 52 | ! Input arguments |
|---|
| 53 | integer,intent(in) :: Npoints !# of grid points |
|---|
| 54 | integer,intent(in) :: Nlevels !# of levels |
|---|
| 55 | integer,intent(in) :: Ncolumns !# of columns |
|---|
| 56 | real(wp),dimension(Npoints,Nlevels),intent(in) :: zfull ! Height at model levels [m] (Bottom of model layer) |
|---|
| 57 | real(wp),dimension(Npoints,Nlevels),intent(in) :: zhalf ! Height at half model levels [m] (Bottom of model layer) |
|---|
| 58 | real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: y ! Variable to be changed to a different grid |
|---|
| 59 | integer,intent(in) :: Nglevels !# levels in the new grid |
|---|
| 60 | real(wp),dimension(Nglevels),intent(in) :: newgrid_bot ! Lower boundary of new levels [m] |
|---|
| 61 | real(wp),dimension(Nglevels),intent(in) :: newgrid_top ! Upper boundary of new levels [m] |
|---|
| 62 | logical,optional,intent(in) :: log_units ! log units, need to convert to linear units |
|---|
| 63 | ! Output |
|---|
| 64 | real(wp),dimension(Npoints,Ncolumns,Nglevels),intent(out) :: r ! Variable on new grid |
|---|
| 65 | |
|---|
| 66 | ! Local variables |
|---|
| 67 | integer :: i,j,k |
|---|
| 68 | logical :: lunits |
|---|
| 69 | integer :: l |
|---|
| 70 | real(wp) :: w ! Weight |
|---|
| 71 | real(wp) :: dbb, dtb, dbt, dtt ! Distances between edges of both grids |
|---|
| 72 | integer :: Nw ! Number of weights |
|---|
| 73 | real(wp) :: wt ! Sum of weights |
|---|
| 74 | real(wp),dimension(Nlevels) :: oldgrid_bot,oldgrid_top ! Lower and upper boundaries of model grid |
|---|
| 75 | real(wp) :: yp ! Local copy of y at a particular point. |
|---|
| 76 | ! This allows for change of units. |
|---|
| 77 | |
|---|
| 78 | lunits=.false. |
|---|
| 79 | if (present(log_units)) lunits=log_units |
|---|
| 80 | |
|---|
| 81 | r = 0._wp |
|---|
| 82 | |
|---|
| 83 | do i=1,Npoints |
|---|
| 84 | ! Calculate tops and bottoms of new and old grids |
|---|
| 85 | oldgrid_bot = zhalf(i,:) |
|---|
| 86 | oldgrid_top(1:Nlevels-1) = oldgrid_bot(2:Nlevels) |
|---|
| 87 | oldgrid_top(Nlevels) = zfull(i,Nlevels) + zfull(i,Nlevels) - zhalf(i,Nlevels) ! Top level symmetric |
|---|
| 88 | l = 0 ! Index of level in the old grid |
|---|
| 89 | ! Loop over levels in the new grid |
|---|
| 90 | do k = 1,Nglevels |
|---|
| 91 | Nw = 0 ! Number of weigths |
|---|
| 92 | wt = 0._wp ! Sum of weights |
|---|
| 93 | ! Loop over levels in the old grid and accumulate total for weighted average |
|---|
| 94 | do |
|---|
| 95 | l = l + 1 |
|---|
| 96 | w = 0.0 ! Initialise weight to 0 |
|---|
| 97 | ! Distances between edges of both grids |
|---|
| 98 | dbb = oldgrid_bot(l) - newgrid_bot(k) |
|---|
| 99 | dtb = oldgrid_top(l) - newgrid_bot(k) |
|---|
| 100 | dbt = oldgrid_bot(l) - newgrid_top(k) |
|---|
| 101 | dtt = oldgrid_top(l) - newgrid_top(k) |
|---|
| 102 | if (dbt >= 0.0) exit ! Do next level in the new grid |
|---|
| 103 | if (dtb > 0.0) then |
|---|
| 104 | if (dbb <= 0.0) then |
|---|
| 105 | if (dtt <= 0) then |
|---|
| 106 | w = dtb |
|---|
| 107 | else |
|---|
| 108 | w = newgrid_top(k) - newgrid_bot(k) |
|---|
| 109 | endif |
|---|
| 110 | else |
|---|
| 111 | if (dtt <= 0) then |
|---|
| 112 | w = oldgrid_top(l) - oldgrid_bot(l) |
|---|
| 113 | else |
|---|
| 114 | w = -dbt |
|---|
| 115 | endif |
|---|
| 116 | endif |
|---|
| 117 | ! If layers overlap (w/=0), then accumulate |
|---|
| 118 | if (w /= 0.0) then |
|---|
| 119 | Nw = Nw + 1 |
|---|
| 120 | wt = wt + w |
|---|
| 121 | do j=1,Ncolumns |
|---|
| 122 | if (lunits) then |
|---|
| 123 | if (y(i,j,l) /= R_UNDEF) then |
|---|
| 124 | yp = 10._wp**(y(i,j,l)/10._wp) |
|---|
| 125 | else |
|---|
| 126 | yp = 0._wp |
|---|
| 127 | endif |
|---|
| 128 | else |
|---|
| 129 | yp = y(i,j,l) |
|---|
| 130 | endif |
|---|
| 131 | r(i,j,k) = r(i,j,k) + w*yp |
|---|
| 132 | enddo |
|---|
| 133 | endif |
|---|
| 134 | endif |
|---|
| 135 | enddo |
|---|
| 136 | l = l - 2 |
|---|
| 137 | if (l < 1) l = 0 |
|---|
| 138 | ! Calculate average in new grid |
|---|
| 139 | if (Nw > 0) then |
|---|
| 140 | do j=1,Ncolumns |
|---|
| 141 | r(i,j,k) = r(i,j,k)/wt |
|---|
| 142 | enddo |
|---|
| 143 | endif |
|---|
| 144 | enddo |
|---|
| 145 | enddo |
|---|
| 146 | |
|---|
| 147 | ! Set points under surface to R_UNDEF, and change to dBZ if necessary |
|---|
| 148 | do k=1,Nglevels |
|---|
| 149 | do j=1,Ncolumns |
|---|
| 150 | do i=1,Npoints |
|---|
| 151 | if (newgrid_top(k) > zhalf(i,1)) then ! Level above model bottom level |
|---|
| 152 | if (lunits) then |
|---|
| 153 | if (r(i,j,k) <= 0.0) then |
|---|
| 154 | r(i,j,k) = R_UNDEF |
|---|
| 155 | else |
|---|
| 156 | r(i,j,k) = 10._wp*log10(r(i,j,k)) |
|---|
| 157 | endif |
|---|
| 158 | endif |
|---|
| 159 | else ! Level below surface |
|---|
| 160 | r(i,j,k) = R_GROUND |
|---|
| 161 | endif |
|---|
| 162 | enddo |
|---|
| 163 | enddo |
|---|
| 164 | enddo |
|---|
| 165 | |
|---|
| 166 | END SUBROUTINE COSP_CHANGE_VERTICAL_GRID |
|---|
| 167 | |
|---|
| 168 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 169 | !------------- SUBROUTINE COSP_LIDAR_ONLY_CLOUD ----------------- |
|---|
| 170 | ! (c) 2008, Lawrence Livermore National Security Limited Liability Corporation. |
|---|
| 171 | ! All rights reserved. |
|---|
| 172 | !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 173 | SUBROUTINE COSP_LIDAR_ONLY_CLOUD(Npoints, Ncolumns, Nlevels, beta_tot, beta_mol, & |
|---|
| 174 | Ze_tot, lidar_only_freq_cloud, tcc, radar_tcc, radar_tcc2) |
|---|
| 175 | ! Inputs |
|---|
| 176 | integer,intent(in) :: & |
|---|
| 177 | Npoints, & ! Number of horizontal gridpoints |
|---|
| 178 | Ncolumns, & ! Number of subcolumns |
|---|
| 179 | Nlevels ! Number of vertical layers |
|---|
| 180 | real(wp),dimension(Npoints,Nlevels),intent(in) :: & |
|---|
| 181 | beta_mol ! Molecular backscatter |
|---|
| 182 | real(wp),dimension(Npoints,Ncolumns,Nlevels),intent(in) :: & |
|---|
| 183 | beta_tot, & ! Total backscattered signal |
|---|
| 184 | Ze_tot ! Radar reflectivity |
|---|
| 185 | ! Outputs |
|---|
| 186 | real(wp),dimension(Npoints,Nlevels),intent(out) :: & |
|---|
| 187 | lidar_only_freq_cloud |
|---|
| 188 | real(wp),dimension(Npoints),intent(out) ::& |
|---|
| 189 | tcc, & ! |
|---|
| 190 | radar_tcc, & ! |
|---|
| 191 | radar_tcc2 ! |
|---|
| 192 | |
|---|
| 193 | ! local variables |
|---|
| 194 | real(wp) :: sc_ratio |
|---|
| 195 | real(wp),parameter :: & |
|---|
| 196 | s_cld=5.0, & |
|---|
| 197 | s_att=0.01 |
|---|
| 198 | integer :: flag_sat,flag_cld,pr,i,j,flag_radarcld,flag_radarcld_no1km,j_1km |
|---|
| 199 | |
|---|
| 200 | lidar_only_freq_cloud = 0._wp |
|---|
| 201 | tcc = 0._wp |
|---|
| 202 | radar_tcc = 0._wp |
|---|
| 203 | radar_tcc2 = 0._wp |
|---|
| 204 | do pr=1,Npoints |
|---|
| 205 | do i=1,Ncolumns |
|---|
| 206 | flag_sat = 0 |
|---|
| 207 | flag_cld = 0 |
|---|
| 208 | flag_radarcld = 0 !+JEK |
|---|
| 209 | flag_radarcld_no1km=0 !+JEK |
|---|
| 210 | ! look for j_1km from bottom to top |
|---|
| 211 | j = 1 |
|---|
| 212 | do while (Ze_tot(pr,i,j) .eq. R_GROUND) |
|---|
| 213 | j = j+1 |
|---|
| 214 | enddo |
|---|
| 215 | j_1km = j+1 !this is the vertical index of 1km above surface |
|---|
| 216 | |
|---|
| 217 | do j=1,Nlevels |
|---|
| 218 | sc_ratio = beta_tot(pr,i,j)/beta_mol(pr,j) |
|---|
| 219 | if ((sc_ratio .le. s_att) .and. (flag_sat .eq. 0)) flag_sat = j |
|---|
| 220 | if (Ze_tot(pr,i,j) .lt. -30.) then !radar can't detect cloud |
|---|
| 221 | if ( (sc_ratio .gt. s_cld) .or. (flag_sat .eq. j) ) then !lidar sense cloud |
|---|
| 222 | lidar_only_freq_cloud(pr,j)=lidar_only_freq_cloud(pr,j)+1. !top->surf |
|---|
| 223 | flag_cld=1 |
|---|
| 224 | endif |
|---|
| 225 | else !radar sense cloud (z%Ze_tot(pr,i,j) .ge. -30.) |
|---|
| 226 | flag_cld=1 |
|---|
| 227 | flag_radarcld=1 |
|---|
| 228 | if (j .gt. j_1km) flag_radarcld_no1km=1 |
|---|
| 229 | endif |
|---|
| 230 | enddo !levels |
|---|
| 231 | if (flag_cld .eq. 1) tcc(pr)=tcc(pr)+1._wp |
|---|
| 232 | if (flag_radarcld .eq. 1) radar_tcc(pr)=radar_tcc(pr)+1. |
|---|
| 233 | if (flag_radarcld_no1km .eq. 1) radar_tcc2(pr)=radar_tcc2(pr)+1. |
|---|
| 234 | enddo !columns |
|---|
| 235 | enddo !points |
|---|
| 236 | lidar_only_freq_cloud=lidar_only_freq_cloud/Ncolumns |
|---|
| 237 | tcc=tcc/Ncolumns |
|---|
| 238 | radar_tcc=radar_tcc/Ncolumns |
|---|
| 239 | radar_tcc2=radar_tcc2/Ncolumns |
|---|
| 240 | |
|---|
| 241 | ! Unit conversion |
|---|
| 242 | where(lidar_only_freq_cloud /= R_UNDEF) & |
|---|
| 243 | lidar_only_freq_cloud = lidar_only_freq_cloud*100._wp |
|---|
| 244 | where(tcc /= R_UNDEF) tcc = tcc*100._wp |
|---|
| 245 | where(radar_tcc /= R_UNDEF) radar_tcc = radar_tcc*100._wp |
|---|
| 246 | where(radar_tcc2 /= R_UNDEF) radar_tcc2 = radar_tcc2*100._wp |
|---|
| 247 | |
|---|
| 248 | END SUBROUTINE COSP_LIDAR_ONLY_CLOUD |
|---|
| 249 | |
|---|
| 250 | ! ###################################################################################### |
|---|
| 251 | ! FUNCTION hist1D |
|---|
| 252 | ! ###################################################################################### |
|---|
| 253 | function hist1d(Npoints,var,nbins,bins) |
|---|
| 254 | ! Inputs |
|---|
| 255 | integer,intent(in) :: & |
|---|
| 256 | Npoints, & ! Number of points in input array |
|---|
| 257 | Nbins ! Number of bins for sorting |
|---|
| 258 | real(wp),intent(in),dimension(Npoints) :: & |
|---|
| 259 | var ! Input variable to be sorted |
|---|
| 260 | real(wp),intent(in),dimension(Nbins+1) :: & |
|---|
| 261 | bins ! Histogram bins [lowest,binTops] |
|---|
| 262 | ! Outputs |
|---|
| 263 | real(wp),dimension(Nbins) :: & |
|---|
| 264 | hist1d ! Output histogram |
|---|
| 265 | ! Local variables |
|---|
| 266 | integer :: ij |
|---|
| 267 | |
|---|
| 268 | do ij=2,Nbins+1 |
|---|
| 269 | hist1D(ij-1) = count(var .ge. bins(ij-1) .and. var .lt. bins(ij)) |
|---|
| 270 | if (count(var .eq. R_GROUND) .ge. 1) hist1D(ij-1)=R_UNDEF |
|---|
| 271 | enddo |
|---|
| 272 | |
|---|
| 273 | end function hist1D |
|---|
| 274 | |
|---|
| 275 | ! ###################################################################################### |
|---|
| 276 | ! SUBROUTINE hist2D |
|---|
| 277 | ! ###################################################################################### |
|---|
| 278 | subroutine hist2D(var1,var2,npts,bin1,nbin1,bin2,nbin2,jointHist) |
|---|
| 279 | implicit none |
|---|
| 280 | |
|---|
| 281 | ! INPUTS |
|---|
| 282 | integer, intent(in) :: & |
|---|
| 283 | npts, & ! Number of data points to be sorted |
|---|
| 284 | nbin1, & ! Number of bins in histogram direction 1 |
|---|
| 285 | nbin2 ! Number of bins in histogram direction 2 |
|---|
| 286 | real(wp),intent(in),dimension(npts) :: & |
|---|
| 287 | var1, & ! Variable 1 to be sorted into bins |
|---|
| 288 | var2 ! variable 2 to be sorted into bins |
|---|
| 289 | real(wp),intent(in),dimension(nbin1+1) :: & |
|---|
| 290 | bin1 ! Histogram bin 1 boundaries |
|---|
| 291 | real(wp),intent(in),dimension(nbin2+1) :: & |
|---|
| 292 | bin2 ! Histogram bin 2 boundaries |
|---|
| 293 | ! OUTPUTS |
|---|
| 294 | real(wp),intent(out),dimension(nbin1,nbin2) :: & |
|---|
| 295 | jointHist |
|---|
| 296 | |
|---|
| 297 | ! LOCAL VARIABLES |
|---|
| 298 | integer :: ij,ik |
|---|
| 299 | |
|---|
| 300 | do ij=2,nbin1+1 |
|---|
| 301 | do ik=2,nbin2+1 |
|---|
| 302 | jointHist(ij-1,ik-1)=count(var1 .ge. bin1(ij-1) .and. var1 .lt. bin1(ij) .and. & |
|---|
| 303 | var2 .ge. bin2(ik-1) .and. var2 .lt. bin2(ik)) |
|---|
| 304 | enddo |
|---|
| 305 | enddo |
|---|
| 306 | end subroutine hist2D |
|---|
| 307 | END MODULE MOD_COSP_STATS |
|---|
