| 1 | ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% |
|---|
| 2 | ! Copyright (c) 2009, Lawrence Livemore National Security Limited Liability |
|---|
| 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_ICARUS |
|---|
| 33 | USE COSP_KINDS, ONLY: wp |
|---|
| 34 | USE COSP_PHYS_CONSTANTS, ONLY: amd,amw,avo,grav |
|---|
| 35 | use MOD_COSP_STATS, ONLY: hist2D |
|---|
| 36 | USE MOD_COSP_CONFIG, ONLY: R_UNDEF,numISCCPTauBins,numISCCPPresBins,isccp_histTau, & |
|---|
| 37 | isccp_histPres |
|---|
| 38 | implicit none |
|---|
| 39 | |
|---|
| 40 | ! Shared Parameters |
|---|
| 41 | integer,parameter :: & |
|---|
| 42 | ncolprint = 0 ! Flag for debug printing (set as parameter to increase performance) |
|---|
| 43 | |
|---|
| 44 | ! Cloud-top height determination |
|---|
| 45 | integer :: & |
|---|
| 46 | isccp_top_height, & ! Top height adjustment method |
|---|
| 47 | isccp_top_height_direction ! Direction for finding atmosphere pressure level |
|---|
| 48 | |
|---|
| 49 | ! Parameters used by icarus |
|---|
| 50 | real(wp),parameter :: & |
|---|
| 51 | tauchk = -1._wp*log(0.9999999_wp), & ! Lower limit on optical depth |
|---|
| 52 | isccp_taumin = 0.3_wp, & ! Minimum optical depth for joint-hostogram |
|---|
| 53 | pstd = 1013250._wp, & ! Mean sea-level pressure (Pa) |
|---|
| 54 | isccp_t0 = 296._wp, & ! Mean surface temperature (K) |
|---|
| 55 | output_missing_value = -1.E+30 ! Missing values |
|---|
| 56 | |
|---|
| 57 | contains |
|---|
| 58 | ! ########################################################################## |
|---|
| 59 | ! ########################################################################## |
|---|
| 60 | SUBROUTINE ICARUS(debug,debugcol,npoints,sunlit,nlev,ncol,pfull, & |
|---|
| 61 | phalf,qv,cc,conv,dtau_s,dtau_c,th,thd,frac_out,skt,emsfc_lw,at,& |
|---|
| 62 | dem_s,dem_c,fq_isccp,totalcldarea, meanptop,meantaucld, & |
|---|
| 63 | meanalbedocld, meantb,meantbclr,boxtau,boxptop,levmatch) |
|---|
| 64 | |
|---|
| 65 | ! INPUTS |
|---|
| 66 | INTEGER,intent(in) :: & ! |
|---|
| 67 | npoints, & ! Number of model points in the horizontal |
|---|
| 68 | nlev, & ! Number of model levels in column |
|---|
| 69 | ncol, & ! Number of subcolumns |
|---|
| 70 | debug, & ! Debug flag |
|---|
| 71 | debugcol ! Debug column flag |
|---|
| 72 | INTEGER,intent(in),dimension(npoints) :: & ! |
|---|
| 73 | sunlit ! 1 for day points, 0 for night time |
|---|
| 74 | REAL(WP),intent(in) :: & ! |
|---|
| 75 | emsfc_lw ! 10.5 micron emissivity of surface (fraction) |
|---|
| 76 | REAL(WP),intent(in),dimension(npoints) :: & ! |
|---|
| 77 | skt ! Skin Temperature (K) |
|---|
| 78 | REAL(WP),intent(in),dimension(npoints,ncol,nlev) :: & ! |
|---|
| 79 | frac_out ! Boxes gridbox divided up into subcolumns |
|---|
| 80 | REAL(WP),intent(in),dimension(npoints,nlev) :: & ! |
|---|
| 81 | pfull, & ! Pressure of full model levels (Pascals) |
|---|
| 82 | qv, & ! Water vapor specific humidity (kg vapor/ kg air) |
|---|
| 83 | cc, & ! Cloud cover in each model level (fraction) |
|---|
| 84 | conv, & ! Convective cloud cover in each model |
|---|
| 85 | at, & ! Temperature in each model level (K) |
|---|
| 86 | dem_c, & ! Emissivity for convective clouds |
|---|
| 87 | dem_s, & ! Emissivity for stratiform clouds |
|---|
| 88 | dtau_c, & ! Optical depth for convective clouds |
|---|
| 89 | dtau_s ! Optical depth for stratiform clouds |
|---|
| 90 | REAL(WP),intent(in),dimension(npoints,nlev+1) :: & ! |
|---|
| 91 | phalf ! Pressure of half model levels (Pascals)! |
|---|
| 92 | integer,intent(in) :: th,thd |
|---|
| 93 | |
|---|
| 94 | ! OUTPUTS |
|---|
| 95 | REAL(WP),intent(out),dimension(npoints,7,7) :: & |
|---|
| 96 | fq_isccp ! The fraction of the model grid box covered by clouds |
|---|
| 97 | REAL(WP),intent(out),dimension(npoints) :: & |
|---|
| 98 | totalcldarea, & ! The fraction of model grid box columns with cloud present |
|---|
| 99 | meanptop, & ! Mean cloud top pressure (mb) - linear averaging |
|---|
| 100 | meantaucld, & ! Mean optical thickness |
|---|
| 101 | meanalbedocld, & ! Mean cloud albedo |
|---|
| 102 | meantb, & ! Mean all-sky 10.5 micron brightness temperature |
|---|
| 103 | meantbclr ! Mean clear-sky 10.5 micron brightness temperature |
|---|
| 104 | REAL(WP),intent(out),dimension(npoints,ncol) :: & |
|---|
| 105 | boxtau, & ! Optical thickness in each column |
|---|
| 106 | boxptop ! Cloud top pressure (mb) in each column |
|---|
| 107 | INTEGER,intent(out),dimension(npoints,ncol) :: & |
|---|
| 108 | levmatch ! Used for icarus unit testing only |
|---|
| 109 | |
|---|
| 110 | |
|---|
| 111 | ! INTERNAL VARIABLES |
|---|
| 112 | CHARACTER(len=10) :: ftn09 |
|---|
| 113 | REAL(WP),dimension(npoints,ncol) :: boxttop |
|---|
| 114 | REAL(WP),dimension(npoints,ncol,nlev) :: dtau,demIN |
|---|
| 115 | INTEGER :: j,ilev,ibox |
|---|
| 116 | INTEGER,dimension(nlev,ncol ) :: acc |
|---|
| 117 | |
|---|
| 118 | ! PARAMETERS |
|---|
| 119 | character ,parameter, dimension(6) :: cchar=(/' ','-','1','+','I','+'/) |
|---|
| 120 | character(len=1),parameter,dimension(6) :: cchar_realtops=(/ ' ',' ','1','1','I','I'/) |
|---|
| 121 | ! ########################################################################## |
|---|
| 122 | |
|---|
| 123 | call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dem_c,dem_s,demIN) |
|---|
| 124 | call cosp_simulator_optics(npoints,ncol,nlev,frac_out,dtau_c,dtau_s,dtau) |
|---|
| 125 | |
|---|
| 126 | call ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demIN,skt,emsfc_lw,qv,at, & |
|---|
| 127 | pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr) |
|---|
| 128 | |
|---|
| 129 | call ICARUS_COLUMN(npoints,ncol,boxtau,boxptop/100._wp,sunlit,boxttop,& |
|---|
| 130 | fq_isccp,meanalbedocld,meanptop,meantaucld,totalcldarea,meantb) |
|---|
| 131 | |
|---|
| 132 | ! ########################################################################## |
|---|
| 133 | ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM |
|---|
| 134 | ! ########################################################################## |
|---|
| 135 | |
|---|
| 136 | if (debugcol.ne.0) then |
|---|
| 137 | do j=1,npoints,debugcol |
|---|
| 138 | |
|---|
| 139 | ! Produce character output |
|---|
| 140 | do ilev=1,nlev |
|---|
| 141 | acc(ilev,1:ncol)=frac_out(j,1:ncol,ilev)*2 |
|---|
| 142 | where(levmatch(j,1:ncol) .eq. ilev) acc(ilev,1:ncol)=acc(ilev,1:ncol)+1 |
|---|
| 143 | enddo |
|---|
| 144 | |
|---|
| 145 | write(ftn09,11) j |
|---|
| 146 | 11 format('ftn09.',i4.4) |
|---|
| 147 | open(9, FILE=ftn09, FORM='FORMATTED') |
|---|
| 148 | |
|---|
| 149 | write(9,'(a1)') ' ' |
|---|
| 150 | write(9,'(10i5)') (ilev,ilev=5,nlev,5) |
|---|
| 151 | write(9,'(a1)') ' ' |
|---|
| 152 | |
|---|
| 153 | do ibox=1,ncol |
|---|
| 154 | write(9,'(40(a1),1x,40(a1))') & |
|---|
| 155 | (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev),& |
|---|
| 156 | (cchar(acc(ilev,ibox)+1),ilev=1,nlev) |
|---|
| 157 | end do |
|---|
| 158 | close(9) |
|---|
| 159 | |
|---|
| 160 | enddo |
|---|
| 161 | end if |
|---|
| 162 | |
|---|
| 163 | return |
|---|
| 164 | end SUBROUTINE ICARUS |
|---|
| 165 | |
|---|
| 166 | ! ############################################################################ |
|---|
| 167 | ! ############################################################################ |
|---|
| 168 | ! ############################################################################ |
|---|
| 169 | SUBROUTINE ICARUS_SUBCOLUMN(npoints,ncol,nlev,sunlit,dtau,demiN,skt,emsfc_lw,qv,at, & |
|---|
| 170 | pfull,phalf,frac_out,levmatch,boxtau,boxptop,boxttop,meantbclr) |
|---|
| 171 | ! Inputs |
|---|
| 172 | INTEGER, intent(in) :: & |
|---|
| 173 | ncol, & ! Number of subcolumns |
|---|
| 174 | npoints, & ! Number of horizontal gridpoints |
|---|
| 175 | nlev ! Number of vertical levels |
|---|
| 176 | INTEGER, intent(in), dimension(npoints) :: & |
|---|
| 177 | sunlit ! 1=day 0=night |
|---|
| 178 | REAL(WP),intent(in) :: & |
|---|
| 179 | emsfc_lw ! 10.5 micron emissivity of surface (fraction) |
|---|
| 180 | REAL(WP),intent(in), dimension(npoints) :: & |
|---|
| 181 | skt ! Skin temperature |
|---|
| 182 | REAL(WP),intent(in), dimension(npoints,nlev) :: & |
|---|
| 183 | at, & ! Temperature |
|---|
| 184 | pfull, & ! Presure |
|---|
| 185 | qv ! Specific humidity |
|---|
| 186 | REAL(WP),intent(in), dimension(npoints,ncol,nlev) :: & |
|---|
| 187 | frac_out, & ! Subcolumn cloud cover |
|---|
| 188 | dtau, & ! Subcolumn optical thickness |
|---|
| 189 | demIN ! Subcolumn emissivity |
|---|
| 190 | REAL(WP),intent(in), dimension(npoints,nlev+1) :: & |
|---|
| 191 | phalf ! Pressure at model half levels |
|---|
| 192 | |
|---|
| 193 | ! Outputs |
|---|
| 194 | REAL(WP),intent(inout),dimension(npoints) :: & |
|---|
| 195 | meantbclr ! Mean clear-sky 10.5 micron brightness temperature |
|---|
| 196 | REAL(WP),intent(inout),dimension(npoints,ncol) :: & |
|---|
| 197 | boxtau, & ! Optical thickness in each column |
|---|
| 198 | boxptop, & ! Cloud top pressure (mb) in each column |
|---|
| 199 | boxttop ! Cloud top temperature in each column |
|---|
| 200 | INTEGER, intent(inout),dimension(npoints,ncol) :: levmatch |
|---|
| 201 | |
|---|
| 202 | ! Local Variables |
|---|
| 203 | INTEGER :: & |
|---|
| 204 | j,ibox,ilev,k1,k2,icycle |
|---|
| 205 | INTEGER,dimension(npoints) :: & |
|---|
| 206 | nmatch,itrop |
|---|
| 207 | INTEGER,dimension(npoints,nlev-1) :: & |
|---|
| 208 | match |
|---|
| 209 | REAL(WP) :: & |
|---|
| 210 | logp,logp1,logp2,atd |
|---|
| 211 | REAL(WP),dimension(npoints) :: & |
|---|
| 212 | bb,attropmin,attrop,ptrop,atmax,btcmin,transmax,tauir,taumin,fluxtopinit,press, & |
|---|
| 213 | dpress,atmden,rvh20,rhoave,rh20s,rfrgn,tmpexp,tauwv,wk,trans_layers_above_clrsky, & |
|---|
| 214 | fluxtop_clrsky |
|---|
| 215 | REAL(WP),dimension(npoints,nlev) :: & |
|---|
| 216 | dem_wv |
|---|
| 217 | REAL(WP),dimension(npoints,ncol) :: & |
|---|
| 218 | trans_layers_above,dem,tb,emcld,fluxtop,tau,ptop |
|---|
| 219 | |
|---|
| 220 | ! #################################################################################### |
|---|
| 221 | ! Compute cloud optical depth for each column by summing up subcolumns |
|---|
| 222 | tau(1:npoints,1:ncol) = 0._wp |
|---|
| 223 | tau(1:npoints,1:ncol) = sum(dtau,dim=3) |
|---|
| 224 | |
|---|
| 225 | ! Set tropopause values |
|---|
| 226 | if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then |
|---|
| 227 | ptrop(1:npoints) = 5000._wp |
|---|
| 228 | attropmin(1:npoints) = 400._wp |
|---|
| 229 | atmax(1:npoints) = 0._wp |
|---|
| 230 | attrop(1:npoints) = 120._wp |
|---|
| 231 | itrop(1:npoints) = 1 |
|---|
| 232 | |
|---|
| 233 | do ilev=1,nlev |
|---|
| 234 | where(pfull(1:npoints,ilev) .lt. 40000. .and. & |
|---|
| 235 | pfull(1:npoints,ilev) .gt. 5000. .and. & |
|---|
| 236 | at(1:npoints,ilev) .lt. attropmin(1:npoints)) |
|---|
| 237 | ptrop(1:npoints) = pfull(1:npoints,ilev) |
|---|
| 238 | attropmin(1:npoints) = at(1:npoints,ilev) |
|---|
| 239 | attrop(1:npoints) = attropmin(1:npoints) |
|---|
| 240 | itrop = ilev |
|---|
| 241 | endwhere |
|---|
| 242 | enddo |
|---|
| 243 | |
|---|
| 244 | do ilev=1,nlev |
|---|
| 245 | atmax(1:npoints) = merge(at(1:npoints,ilev),atmax(1:npoints),& |
|---|
| 246 | at(1:npoints,ilev) .gt. atmax(1:npoints) .and. ilev .ge. itrop(1:npoints)) |
|---|
| 247 | enddo |
|---|
| 248 | end if |
|---|
| 249 | |
|---|
| 250 | if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then |
|---|
| 251 | ! ############################################################################ |
|---|
| 252 | ! Clear-sky radiance calculation |
|---|
| 253 | ! |
|---|
| 254 | ! Compute water vapor continuum emissivity this treatment follows Schwarkzopf |
|---|
| 255 | ! and Ramasamy JGR 1999,vol 104, pages 9467-9499. The emissivity is calculated |
|---|
| 256 | ! at a wavenumber of 955 cm-1, or 10.47 microns |
|---|
| 257 | ! ############################################################################ |
|---|
| 258 | do ilev=1,nlev |
|---|
| 259 | press(1:npoints) = pfull(1:npoints,ilev)*10._wp |
|---|
| 260 | dpress(1:npoints) = (phalf(1:npoints,ilev+1)-phalf(1:npoints,ilev))*10 |
|---|
| 261 | atmden(1:npoints) = dpress(1:npoints)/(grav*100._wp) |
|---|
| 262 | rvh20(1:npoints) = qv(1:npoints,ilev)*amd/amw |
|---|
| 263 | wk(1:npoints) = rvh20(1:npoints)*avo*atmden/amd |
|---|
| 264 | rhoave(1:npoints) = (press(1:npoints)/pstd)*(isccp_t0/at(1:npoints,ilev)) |
|---|
| 265 | rh20s(1:npoints) = rvh20(1:npoints)*rhoave(1:npoints) |
|---|
| 266 | rfrgn(1:npoints) = rhoave(1:npoints)-rh20s(1:npoints) |
|---|
| 267 | tmpexp(1:npoints) = exp(-0.02_wp*(at(1:npoints,ilev)-isccp_t0)) |
|---|
| 268 | tauwv(1:npoints) = wk(1:npoints)*1.e-20*((0.0224697_wp*rh20s(1:npoints)* & |
|---|
| 269 | tmpexp(1:npoints))+(3.41817e-7*rfrgn(1:npoints)))*0.98_wp |
|---|
| 270 | dem_wv(1:npoints,ilev) = 1._wp - exp( -1._wp * tauwv(1:npoints)) |
|---|
| 271 | enddo |
|---|
| 272 | |
|---|
| 273 | fluxtop_clrsky(1:npoints) = 0._wp |
|---|
| 274 | trans_layers_above_clrsky(1:npoints) = 1._wp |
|---|
| 275 | do ilev=1,nlev |
|---|
| 276 | ! Black body emission at temperature of the layer |
|---|
| 277 | bb(1:npoints) = 1._wp / ( exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp ) |
|---|
| 278 | |
|---|
| 279 | ! Increase TOA flux by flux emitted from layer times total transmittance in layers above |
|---|
| 280 | fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + & |
|---|
| 281 | dem_wv(1:npoints,ilev)*bb(1:npoints)*trans_layers_above_clrsky(1:npoints) |
|---|
| 282 | |
|---|
| 283 | ! Update trans_layers_above with transmissivity from this layer for next time around loop |
|---|
| 284 | trans_layers_above_clrsky(1:npoints) = trans_layers_above_clrsky(1:npoints)*& |
|---|
| 285 | (1.-dem_wv(1:npoints,ilev)) |
|---|
| 286 | enddo |
|---|
| 287 | |
|---|
| 288 | ! Add in surface emission |
|---|
| 289 | bb(1:npoints) = 1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp ) |
|---|
| 290 | fluxtop_clrsky(1:npoints) = fluxtop_clrsky(1:npoints) + & |
|---|
| 291 | emsfc_lw * bb(1:npoints)*trans_layers_above_clrsky(1:npoints) |
|---|
| 292 | |
|---|
| 293 | ! Clear Sky brightness temperature |
|---|
| 294 | meantbclr(1:npoints) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop_clrsky(1:npoints)))) |
|---|
| 295 | |
|---|
| 296 | ! ################################################################################# |
|---|
| 297 | ! All-sky radiance calculation |
|---|
| 298 | ! ################################################################################# |
|---|
| 299 | |
|---|
| 300 | fluxtop(1:npoints,1:ncol) = 0._wp |
|---|
| 301 | trans_layers_above(1:npoints,1:ncol) = 1._wp |
|---|
| 302 | do ilev=1,nlev |
|---|
| 303 | ! Black body emission at temperature of the layer |
|---|
| 304 | bb=1._wp/(exp(1307.27_wp/at(1:npoints,ilev)) - 1._wp) |
|---|
| 305 | |
|---|
| 306 | do ibox=1,ncol |
|---|
| 307 | ! Emissivity |
|---|
| 308 | dem(1:npoints,ibox) = merge(dem_wv(1:npoints,ilev), & |
|---|
| 309 | 1._wp-(1._wp-demIN(1:npoints,ibox,ilev))*(1._wp-dem_wv(1:npoints,ilev)), & |
|---|
| 310 | demIN(1:npoints,ibox,ilev) .eq. 0) |
|---|
| 311 | |
|---|
| 312 | ! Increase TOA flux emitted from layer |
|---|
| 313 | fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + dem(1:npoints,ibox)*bb*trans_layers_above(1:npoints,ibox) |
|---|
| 314 | |
|---|
| 315 | ! Update trans_layer by emitted layer from above |
|---|
| 316 | trans_layers_above(1:npoints,ibox) = trans_layers_above(1:npoints,ibox)*(1._wp-dem(1:npoints,ibox)) |
|---|
| 317 | enddo |
|---|
| 318 | enddo |
|---|
| 319 | |
|---|
| 320 | ! Add in surface emission |
|---|
| 321 | bb(1:npoints)=1._wp/( exp(1307.27_wp/skt(1:npoints)) - 1._wp ) |
|---|
| 322 | do ibox=1,ncol |
|---|
| 323 | fluxtop(1:npoints,ibox) = fluxtop(1:npoints,ibox) + emsfc_lw*bb(1:npoints)*trans_layers_above(1:npoints,ibox) |
|---|
| 324 | end do |
|---|
| 325 | |
|---|
| 326 | ! All Sky brightness temperature |
|---|
| 327 | boxttop(1:npoints,1:ncol) = 1307.27_wp/(log(1._wp+(1._wp/fluxtop(1:npoints,1:ncol)))) |
|---|
| 328 | |
|---|
| 329 | ! ################################################################################# |
|---|
| 330 | ! Cloud-Top Temperature |
|---|
| 331 | ! |
|---|
| 332 | ! Now that you have the top of atmosphere radiance, account for ISCCP |
|---|
| 333 | ! procedures to determine cloud top temperature account for partially |
|---|
| 334 | ! transmitting cloud recompute flux ISCCP would see assuming a single layer |
|---|
| 335 | ! cloud. *NOTE* choice here of 2.13, as it is primarily ice clouds which have |
|---|
| 336 | ! partial emissivity and need the adjustment performed in this section. If it |
|---|
| 337 | ! turns out that the cloud brightness temperature is greater than 260K, then |
|---|
| 338 | ! the liquid cloud conversion factor of 2.56 is used. *NOTE* that this is |
|---|
| 339 | ! discussed on pages 85-87 of the ISCCP D level documentation |
|---|
| 340 | ! (Rossow et al. 1996) |
|---|
| 341 | ! ################################################################################# |
|---|
| 342 | |
|---|
| 343 | ! Compute minimum brightness temperature and optical depth |
|---|
| 344 | btcmin(1:npoints) = 1._wp / ( exp(1307.27_wp/(attrop(1:npoints)-5._wp)) - 1._wp ) |
|---|
| 345 | |
|---|
| 346 | do ibox=1,ncol |
|---|
| 347 | transmax(1:npoints) = (fluxtop(1:npoints,ibox)-btcmin) /(fluxtop_clrsky(1:npoints)-btcmin(1:npoints)) |
|---|
| 348 | tauir(1:npoints) = tau(1:npoints,ibox)/2.13_wp |
|---|
| 349 | taumin(1:npoints) = -log(max(min(transmax(1:npoints),0.9999999_wp),0.001_wp)) |
|---|
| 350 | if (isccp_top_height .eq. 1) then |
|---|
| 351 | do j=1,npoints |
|---|
| 352 | if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then |
|---|
| 353 | fluxtopinit(j) = fluxtop(j,ibox) |
|---|
| 354 | tauir(j) = tau(j,ibox)/2.13_wp |
|---|
| 355 | endif |
|---|
| 356 | enddo |
|---|
| 357 | do icycle=1,2 |
|---|
| 358 | do j=1,npoints |
|---|
| 359 | if (tau(j,ibox) .gt. (tauchk)) then |
|---|
| 360 | if (transmax(j) .gt. 0.001 .and. transmax(j) .le. 0.9999999) then |
|---|
| 361 | emcld(j,ibox) = 1._wp - exp(-1._wp * tauir(j) ) |
|---|
| 362 | fluxtop(j,ibox) = fluxtopinit(j) - ((1.-emcld(j,ibox))*fluxtop_clrsky(j)) |
|---|
| 363 | fluxtop(j,ibox)=max(1.E-06_wp,(fluxtop(j,ibox)/emcld(j,ibox))) |
|---|
| 364 | tb(j,ibox)= 1307.27_wp / (log(1._wp + (1._wp/fluxtop(j,ibox)))) |
|---|
| 365 | if (tb(j,ibox) .gt. 260.) then |
|---|
| 366 | tauir(j) = tau(j,ibox) / 2.56_wp |
|---|
| 367 | end if |
|---|
| 368 | end if |
|---|
| 369 | end if |
|---|
| 370 | enddo |
|---|
| 371 | enddo |
|---|
| 372 | endif |
|---|
| 373 | |
|---|
| 374 | ! Cloud-top temperature |
|---|
| 375 | where(tau(1:npoints,ibox) .gt. tauchk) |
|---|
| 376 | tb(1:npoints,ibox)= 1307.27_wp/ (log(1. + (1._wp/fluxtop(1:npoints,ibox)))) |
|---|
| 377 | where (isccp_top_height .eq. 1 .and. tauir(1:npoints) .lt. taumin(1:npoints)) |
|---|
| 378 | tb(1:npoints,ibox) = attrop(1:npoints) - 5._wp |
|---|
| 379 | tau(1:npoints,ibox) = 2.13_wp*taumin(1:npoints) |
|---|
| 380 | endwhere |
|---|
| 381 | endwhere |
|---|
| 382 | |
|---|
| 383 | ! Clear-sky brightness temperature |
|---|
| 384 | where(tau(1:npoints,ibox) .le. tauchk) |
|---|
| 385 | tb(1:npoints,ibox) = meantbclr(1:npoints) |
|---|
| 386 | endwhere |
|---|
| 387 | enddo |
|---|
| 388 | else |
|---|
| 389 | meantbclr(1:npoints) = output_missing_value |
|---|
| 390 | end if |
|---|
| 391 | |
|---|
| 392 | ! #################################################################################### |
|---|
| 393 | ! Cloud-Top Pressure |
|---|
| 394 | ! |
|---|
| 395 | ! The 2 methods differ according to whether or not you use the physical cloud |
|---|
| 396 | ! top pressure (isccp_top_height = 2) or the radiatively determined cloud top |
|---|
| 397 | ! pressure (isccp_top_height = 1 or 3) |
|---|
| 398 | ! #################################################################################### |
|---|
| 399 | do ibox=1,ncol |
|---|
| 400 | !segregate according to optical thickness |
|---|
| 401 | if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then |
|---|
| 402 | |
|---|
| 403 | ! Find level whose temperature most closely matches brightness temperature |
|---|
| 404 | nmatch(1:npoints)=0 |
|---|
| 405 | do k1=1,nlev-1 |
|---|
| 406 | ilev = merge(nlev-k1,k1,isccp_top_height_direction .eq. 2) |
|---|
| 407 | do j=1,npoints |
|---|
| 408 | if (ilev .ge. itrop(j) .and. & |
|---|
| 409 | ((at(j,ilev) .ge. tb(j,ibox) .and. & |
|---|
| 410 | at(j,ilev+1) .le. tb(j,ibox)) .or. & |
|---|
| 411 | (at(j,ilev) .le. tb(j,ibox) .and. & |
|---|
| 412 | at(j,ilev+1) .ge. tb(j,ibox)))) then |
|---|
| 413 | nmatch(j)=nmatch(j)+1 |
|---|
| 414 | match(j,nmatch(j))=ilev |
|---|
| 415 | endif |
|---|
| 416 | enddo |
|---|
| 417 | enddo |
|---|
| 418 | |
|---|
| 419 | do j=1,npoints |
|---|
| 420 | if (nmatch(j) .ge. 1) then |
|---|
| 421 | k1 = match(j,nmatch(j)) |
|---|
| 422 | k2 = k1 + 1 |
|---|
| 423 | logp1 = log(pfull(j,k1)) |
|---|
| 424 | logp2 = log(pfull(j,k2)) |
|---|
| 425 | atd = max(tauchk,abs(at(j,k2) - at(j,k1))) |
|---|
| 426 | logp=logp1+(logp2-logp1)*abs(tb(j,ibox)-at(j,k1))/atd |
|---|
| 427 | ptop(j,ibox) = exp(logp) |
|---|
| 428 | levmatch(j,ibox) = merge(k1,k2,abs(pfull(j,k1)-ptop(j,ibox)) .lt. abs(pfull(j,k2)-ptop(j,ibox))) |
|---|
| 429 | else |
|---|
| 430 | if (tb(j,ibox) .le. attrop(j)) then |
|---|
| 431 | ptop(j,ibox)=ptrop(j) |
|---|
| 432 | levmatch(j,ibox)=itrop(j) |
|---|
| 433 | end if |
|---|
| 434 | if (tb(j,ibox) .ge. atmax(j)) then |
|---|
| 435 | ptop(j,ibox)=pfull(j,nlev) |
|---|
| 436 | levmatch(j,ibox)=nlev |
|---|
| 437 | end if |
|---|
| 438 | end if |
|---|
| 439 | enddo |
|---|
| 440 | else |
|---|
| 441 | ptop(1:npoints,ibox)=0. |
|---|
| 442 | do ilev=1,nlev |
|---|
| 443 | where((ptop(1:npoints,ibox) .eq. 0. ) .and.(frac_out(1:npoints,ibox,ilev) .ne. 0)) |
|---|
| 444 | ptop(1:npoints,ibox)=phalf(1:npoints,ilev) |
|---|
| 445 | levmatch(1:npoints,ibox)=ilev |
|---|
| 446 | endwhere |
|---|
| 447 | end do |
|---|
| 448 | end if |
|---|
| 449 | where(tau(1:npoints,ibox) .le. tauchk) |
|---|
| 450 | ptop(1:npoints,ibox)=0._wp |
|---|
| 451 | levmatch(1:npoints,ibox)=0._wp |
|---|
| 452 | endwhere |
|---|
| 453 | enddo |
|---|
| 454 | |
|---|
| 455 | ! #################################################################################### |
|---|
| 456 | ! Compute subcolumn pressure and optical depth |
|---|
| 457 | ! #################################################################################### |
|---|
| 458 | boxtau(1:npoints,1:ncol) = output_missing_value |
|---|
| 459 | boxptop(1:npoints,1:ncol) = output_missing_value |
|---|
| 460 | do ibox=1,ncol |
|---|
| 461 | do j=1,npoints |
|---|
| 462 | if (tau(j,ibox) .gt. (tauchk) .and. ptop(j,ibox) .gt. 0.) then |
|---|
| 463 | if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then |
|---|
| 464 | boxtau(j,ibox) = tau(j,ibox) |
|---|
| 465 | boxptop(j,ibox) = ptop(j,ibox)!/100._wp |
|---|
| 466 | endif |
|---|
| 467 | endif |
|---|
| 468 | enddo |
|---|
| 469 | enddo |
|---|
| 470 | |
|---|
| 471 | end SUBROUTINE ICARUS_SUBCOLUMN |
|---|
| 472 | |
|---|
| 473 | ! ###################################################################################### |
|---|
| 474 | ! SUBROUTINE icarus_column |
|---|
| 475 | ! ###################################################################################### |
|---|
| 476 | SUBROUTINE ICARUS_column(npoints,ncol,boxtau,boxptop,sunlit,boxttop,fq_isccp, & |
|---|
| 477 | meanalbedocld,meanptop,meantaucld,totalcldarea,meantb) |
|---|
| 478 | ! Inputs |
|---|
| 479 | INTEGER, intent(in) :: & |
|---|
| 480 | ncol, & ! Number of subcolumns |
|---|
| 481 | npoints ! Number of horizontal gridpoints |
|---|
| 482 | INTEGER, intent(in),dimension(npoints) :: & |
|---|
| 483 | sunlit ! day=1 night=0 |
|---|
| 484 | REAL(WP),intent(in),dimension(npoints,ncol) :: & |
|---|
| 485 | boxttop, & ! Subcolumn top temperature |
|---|
| 486 | boxptop, & ! Subcolumn cloud top pressure |
|---|
| 487 | boxtau ! Subcolumn optical depth |
|---|
| 488 | |
|---|
| 489 | ! Outputs |
|---|
| 490 | REAL(WP),intent(inout),dimension(npoints) :: & |
|---|
| 491 | meanalbedocld, & ! Gridmean cloud albedo |
|---|
| 492 | meanptop, & ! Gridmean cloud top pressure (mb) - linear averaging |
|---|
| 493 | meantaucld, & ! Gridmean optical thickness |
|---|
| 494 | totalcldarea, & ! The fraction of model grid box columns with cloud present |
|---|
| 495 | meantb ! Gridmean all-sky 10.5 micron brightness temperature |
|---|
| 496 | REAL(WP),intent(inout),dimension(npoints,7,7) :: & |
|---|
| 497 | fq_isccp ! The fraction of the model grid box covered by clouds |
|---|
| 498 | |
|---|
| 499 | ! Local Variables |
|---|
| 500 | INTEGER :: j,ilev,ilev2 |
|---|
| 501 | REAL(WP),dimension(npoints,ncol) :: albedocld |
|---|
| 502 | LOGICAL, dimension(npoints,ncol) :: box_cloudy |
|---|
| 503 | |
|---|
| 504 | ! Variables for new joint-histogram implementation |
|---|
| 505 | logical,dimension(ncol) :: box_cloudy2 |
|---|
| 506 | |
|---|
| 507 | ! #################################################################################### |
|---|
| 508 | ! Brightness Temperature |
|---|
| 509 | ! #################################################################################### |
|---|
| 510 | if (isccp_top_height .eq. 1 .or. isccp_top_height .eq. 3) then |
|---|
| 511 | meantb(1:npoints)=sum(boxttop,2)/ncol |
|---|
| 512 | else |
|---|
| 513 | meantb(1:npoints) = output_missing_value |
|---|
| 514 | endif |
|---|
| 515 | |
|---|
| 516 | ! #################################################################################### |
|---|
| 517 | ! Determines ISCCP cloud type frequencies |
|---|
| 518 | ! |
|---|
| 519 | ! Now that boxptop and boxtau have been determined, determine amount of each of the |
|---|
| 520 | ! 49 ISCCP cloud types. Also compute grid box mean cloud top pressure and |
|---|
| 521 | ! optical thickness. The mean cloud top pressure and optical thickness are |
|---|
| 522 | ! averages over the cloudy area only. The mean cloud top pressure is a linear |
|---|
| 523 | ! average of the cloud top pressures. The mean cloud optical thickness is |
|---|
| 524 | ! computed by converting optical thickness to an albedo, averaging in albedo |
|---|
| 525 | ! units, then converting the average albedo back to a mean optical thickness. |
|---|
| 526 | ! #################################################################################### |
|---|
| 527 | |
|---|
| 528 | ! Initialize |
|---|
| 529 | albedocld(1:npoints,1:ncol) = 0._wp |
|---|
| 530 | box_cloudy(1:npoints,1:ncol) = .false. |
|---|
| 531 | |
|---|
| 532 | ! Reset frequencies |
|---|
| 533 | !fq_isccp = spread(spread(merge(0._wp,output_missing_value,sunlit .eq. 1 .or. isccp_top_height .eq. 3),2,7),2,7) |
|---|
| 534 | do ilev=1,7 |
|---|
| 535 | do ilev2=1,7 |
|---|
| 536 | do j=1,npoints ! |
|---|
| 537 | if (sunlit(j).eq.1 .or. isccp_top_height .eq. 3) then |
|---|
| 538 | fq_isccp(j,ilev,ilev2)= 0. |
|---|
| 539 | else |
|---|
| 540 | fq_isccp(j,ilev,ilev2)= output_missing_value |
|---|
| 541 | end if |
|---|
| 542 | enddo |
|---|
| 543 | enddo |
|---|
| 544 | enddo |
|---|
| 545 | |
|---|
| 546 | |
|---|
| 547 | ! Reset variables need for averaging cloud properties |
|---|
| 548 | where(sunlit .eq. 1 .or. isccp_top_height .eq. 3) |
|---|
| 549 | totalcldarea(1:npoints) = 0._wp |
|---|
| 550 | meanalbedocld(1:npoints) = 0._wp |
|---|
| 551 | meanptop(1:npoints) = 0._wp |
|---|
| 552 | meantaucld(1:npoints) = 0._wp |
|---|
| 553 | elsewhere |
|---|
| 554 | totalcldarea(1:npoints) = output_missing_value |
|---|
| 555 | meanalbedocld(1:npoints) = output_missing_value |
|---|
| 556 | meanptop(1:npoints) = output_missing_value |
|---|
| 557 | meantaucld(1:npoints) = output_missing_value |
|---|
| 558 | endwhere |
|---|
| 559 | |
|---|
| 560 | ! Compute column quantities and joint-histogram |
|---|
| 561 | do j=1,npoints |
|---|
| 562 | ! Subcolumns that are cloudy(true) and not(false) |
|---|
| 563 | box_cloudy2(1:ncol) = merge(.true.,.false.,boxtau(j,1:ncol) .gt. tauchk .and. boxptop(j,1:ncol) .gt. 0.) |
|---|
| 564 | |
|---|
| 565 | ! Compute joint histogram and column quantities for points that are sunlit and cloudy |
|---|
| 566 | if (sunlit(j) .eq.1 .or. isccp_top_height .eq. 3) then |
|---|
| 567 | ! Joint-histogram |
|---|
| 568 | call hist2D(boxtau(j,1:ncol),boxptop(j,1:ncol),ncol,isccp_histTau,numISCCPTauBins, & |
|---|
| 569 | isccp_histPres,numISCCPPresBins,fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)) |
|---|
| 570 | fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins) = & |
|---|
| 571 | fq_isccp(j,1:numISCCPTauBins,1:numISCCPPresBins)/ncol |
|---|
| 572 | |
|---|
| 573 | ! Column cloud area |
|---|
| 574 | totalcldarea(j) = real(count(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin))/ncol |
|---|
| 575 | |
|---|
| 576 | ! Subcolumn cloud albedo |
|---|
| 577 | !albedocld(j,1:ncol) = merge((boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp),& |
|---|
| 578 | ! 0._wp,box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin) |
|---|
| 579 | where(box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin) |
|---|
| 580 | albedocld(j,1:ncol) = (boxtau(j,1:ncol)**0.895_wp)/((boxtau(j,1:ncol)**0.895_wp)+6.82_wp) |
|---|
| 581 | elsewhere |
|---|
| 582 | albedocld(j,1:ncol) = 0._wp |
|---|
| 583 | endwhere |
|---|
| 584 | |
|---|
| 585 | ! Column cloud albedo |
|---|
| 586 | meanalbedocld(j) = sum(albedocld(j,1:ncol))/ncol |
|---|
| 587 | |
|---|
| 588 | ! Column cloud top pressure |
|---|
| 589 | meanptop(j) = sum(boxptop(j,1:ncol),box_cloudy2(1:ncol) .and. boxtau(j,1:ncol) .gt. isccp_taumin)/ncol |
|---|
| 590 | endif |
|---|
| 591 | enddo |
|---|
| 592 | |
|---|
| 593 | ! Compute mean cloud properties. Set to mssing value in the event that totalcldarea=0 |
|---|
| 594 | where(totalcldarea(1:npoints) .gt. 0) |
|---|
| 595 | meanptop(1:npoints) = 100._wp*meanptop(1:npoints)/totalcldarea(1:npoints) |
|---|
| 596 | meanalbedocld(1:npoints) = meanalbedocld(1:npoints)/totalcldarea(1:npoints) |
|---|
| 597 | meantaucld(1:npoints) = (6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp) |
|---|
| 598 | elsewhere |
|---|
| 599 | meanptop(1:nPoints) = output_missing_value |
|---|
| 600 | meanalbedocld(1:nPoints) = output_missing_value |
|---|
| 601 | meantaucld(1:nPoints) = output_missing_value |
|---|
| 602 | endwhere |
|---|
| 603 | !meanptop(1:npoints) = merge(100._wp*meanptop(1:npoints)/totalcldarea(1:npoints),& |
|---|
| 604 | ! output_missing_value,totalcldarea(1:npoints) .gt. 0) |
|---|
| 605 | !meanalbedocld(1:npoints) = merge(meanalbedocld(1:npoints)/totalcldarea(1:npoints), & |
|---|
| 606 | ! output_missing_value,totalcldarea(1:npoints) .gt. 0) |
|---|
| 607 | !meantaucld(1:npoints) = merge((6.82_wp/((1._wp/meanalbedocld(1:npoints))-1.))**(1._wp/0.895_wp), & |
|---|
| 608 | ! output_missing_value,totalcldarea(1:npoints) .gt. 0) |
|---|
| 609 | |
|---|
| 610 | ! Represent in percent |
|---|
| 611 | where(totalcldarea .ne. output_missing_value) totalcldarea = totalcldarea*100._wp |
|---|
| 612 | where(fq_isccp .ne. output_missing_value) fq_isccp = fq_isccp*100._wp |
|---|
| 613 | |
|---|
| 614 | |
|---|
| 615 | end SUBROUTINE ICARUS_column |
|---|
| 616 | |
|---|
| 617 | subroutine cosp_simulator_optics(dim1,dim2,dim3,flag,varIN1,varIN2,varOUT) |
|---|
| 618 | ! INPUTS |
|---|
| 619 | integer,intent(in) :: & |
|---|
| 620 | dim1, & ! Dimension 1 extent (Horizontal) |
|---|
| 621 | dim2, & ! Dimension 2 extent (Subcolumn) |
|---|
| 622 | dim3 ! Dimension 3 extent (Vertical) |
|---|
| 623 | real(wp),intent(in),dimension(dim1,dim2,dim3) :: & |
|---|
| 624 | flag ! Logical to determine the of merge var1IN and var2IN |
|---|
| 625 | real(wp),intent(in),dimension(dim1, dim3) :: & |
|---|
| 626 | varIN1, & ! Input field 1 |
|---|
| 627 | varIN2 ! Input field 2 |
|---|
| 628 | ! OUTPUTS |
|---|
| 629 | real(wp),intent(out),dimension(dim1,dim2,dim3) :: & |
|---|
| 630 | varOUT ! Merged output field |
|---|
| 631 | ! LOCAL VARIABLES |
|---|
| 632 | integer :: j |
|---|
| 633 | |
|---|
| 634 | varOUT(1:dim1,1:dim2,1:dim3) = 0._wp |
|---|
| 635 | do j=1,dim2 |
|---|
| 636 | where(flag(:,j,:) .eq. 1) |
|---|
| 637 | varOUT(:,j,:) = varIN2 |
|---|
| 638 | endwhere |
|---|
| 639 | where(flag(:,j,:) .eq. 2) |
|---|
| 640 | varOUT(:,j,:) = varIN1 |
|---|
| 641 | endwhere |
|---|
| 642 | enddo |
|---|
| 643 | end subroutine cosp_simulator_optics |
|---|
| 644 | end module MOD_ICARUS |
|---|
| 645 | |
|---|
