[524] | 1 | ! |
---|
[1344] | 2 | ! $Id $ |
---|
[524] | 3 | ! |
---|
| 4 | SUBROUTINE ISCCP_CLOUD_TYPES( |
---|
| 5 | & debug, |
---|
| 6 | & debugcol, |
---|
| 7 | & npoints, |
---|
| 8 | & sunlit, |
---|
| 9 | & nlev, |
---|
| 10 | & ncol, |
---|
| 11 | & seed, |
---|
| 12 | & pfull, |
---|
| 13 | & phalf, |
---|
| 14 | & qv, |
---|
| 15 | & cc, |
---|
| 16 | & conv, |
---|
| 17 | & dtau_s, |
---|
| 18 | & dtau_c, |
---|
| 19 | & top_height, |
---|
| 20 | & overlap, |
---|
| 21 | & tautab, |
---|
| 22 | & invtau, |
---|
| 23 | & skt, |
---|
| 24 | & emsfc_lw, |
---|
| 25 | & at,dem_s,dem_c, |
---|
| 26 | & fq_isccp, |
---|
| 27 | & totalcldarea, |
---|
| 28 | & meanptop, |
---|
| 29 | & meantaucld, |
---|
| 30 | & boxtau, |
---|
| 31 | & boxptop |
---|
| 32 | &) |
---|
[1344] | 33 | |
---|
[524] | 34 | |
---|
| 35 | ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved. |
---|
| 36 | ! |
---|
| 37 | ! This code is available without charge with the following conditions: |
---|
| 38 | ! |
---|
| 39 | ! 1. The code is available for scientific purposes and is not for |
---|
| 40 | ! commercial use. |
---|
| 41 | ! 2. Any improvements you make to the code should be made available |
---|
| 42 | ! to the to the authors for incorporation into a future release. |
---|
| 43 | ! 3. The code should not be used in any way that brings the authors |
---|
| 44 | ! or their employers into disrepute. |
---|
| 45 | |
---|
| 46 | implicit none |
---|
| 47 | |
---|
| 48 | ! NOTE: the maximum number of levels and columns is set by |
---|
| 49 | ! the following parameter statement |
---|
| 50 | |
---|
| 51 | INTEGER ncolprint |
---|
| 52 | |
---|
| 53 | ! ----- |
---|
| 54 | ! Input |
---|
| 55 | ! ----- |
---|
| 56 | |
---|
| 57 | INTEGER npoints ! number of model points in the horizontal |
---|
| 58 | c PARAMETER(npoints=6722) |
---|
| 59 | INTEGER nlev ! number of model levels in column |
---|
| 60 | INTEGER ncol ! number of subcolumns |
---|
| 61 | |
---|
| 62 | INTEGER sunlit(npoints) ! 1 for day points, 0 for night time |
---|
| 63 | |
---|
| 64 | INTEGER seed(npoints) ! seed value for random number generator |
---|
| 65 | c ! ( see Numerical Recipes Chapter 7) |
---|
| 66 | c ! It is recommended that the seed is set |
---|
| 67 | c ! to a different value for each model |
---|
| 68 | c ! gridbox it is called on, as it is |
---|
[1344] | 69 | c ! possible that the choice of the samec |
---|
| 70 | c ! seed value every time may introduce some |
---|
| 71 | c ! statistical bias in the results, particularly |
---|
| 72 | c ! for low values of NCOL. |
---|
[524] | 73 | c |
---|
[1344] | 74 | REAL pfull(npoints,nlev) ! pressure of full model levels (Pascals) |
---|
[524] | 75 | c ! pfull(npoints,1) is top level of model |
---|
| 76 | c ! pfull(npoints,nlev) is bottom level of model |
---|
| 77 | |
---|
| 78 | REAL phalf(npoints,nlev+1) ! pressure of half model levels (Pascals) |
---|
| 79 | c ! phalf(npoints,1) is top of model |
---|
| 80 | c ! phalf(npoints,nlev+1) is the surface pressure |
---|
| 81 | |
---|
| 82 | REAL qv(npoints,nlev) ! water vapor specific humidity (kg vapor/ kg air) |
---|
| 83 | c ! on full model levels |
---|
| 84 | |
---|
| 85 | REAL cc(npoints,nlev) ! input cloud cover in each model level (fraction) |
---|
| 86 | c ! NOTE: This is the HORIZONTAL area of each |
---|
| 87 | c ! grid box covered by clouds |
---|
| 88 | |
---|
| 89 | REAL conv(npoints,nlev) ! input convective cloud cover in each model level (fraction) |
---|
| 90 | c ! NOTE: This is the HORIZONTAL area of each |
---|
| 91 | c ! grid box covered by convective clouds |
---|
| 92 | |
---|
| 93 | REAL dtau_s(npoints,nlev) ! mean 0.67 micron optical depth of stratiform |
---|
[1344] | 94 | c ! clouds in each model level |
---|
[524] | 95 | c ! NOTE: this the cloud optical depth of only the |
---|
| 96 | c ! cloudy part of the grid box, it is not weighted |
---|
| 97 | c ! with the 0 cloud optical depth of the clear |
---|
| 98 | c ! part of the grid box |
---|
| 99 | |
---|
| 100 | REAL dtau_c(npoints,nlev) ! mean 0.67 micron optical depth of convective |
---|
[1344] | 101 | c ! clouds in each |
---|
[524] | 102 | c ! model level. Same note applies as in dtau_s. |
---|
| 103 | |
---|
| 104 | INTEGER overlap ! overlap type |
---|
[1344] | 105 | |
---|
[524] | 106 | ! 1=max |
---|
[1344] | 107 | |
---|
[524] | 108 | ! 2=rand |
---|
| 109 | ! 3=max/rand |
---|
| 110 | |
---|
| 111 | INTEGER top_height ! 1 = adjust top height using both a computed |
---|
| 112 | c ! infrared brightness temperature and the visible |
---|
[1344] | 113 | c ! optical depth to adjust cloud top pressure. Note |
---|
| 114 | c ! that this calculation is most appropriate to compare |
---|
| 115 | c ! to ISCCP data during sunlit hours. |
---|
[524] | 116 | c ! 2 = do not adjust top height, that is cloud top |
---|
| 117 | c ! pressure is the actual cloud top pressure |
---|
| 118 | c ! in the model |
---|
[1344] | 119 | c ! 3 = adjust top height using only the computed |
---|
| 120 | c ! infrared brightness temperature. Note that this |
---|
| 121 | c ! calculation is most appropriate to compare to ISCCP |
---|
| 122 | c ! IR only algortihm (i.e. you can compare to nighttime |
---|
| 123 | c ! ISCCP data with this option) |
---|
[524] | 124 | |
---|
| 125 | REAL tautab(0:255) ! ISCCP table for converting count value to |
---|
| 126 | c ! optical thickness |
---|
| 127 | |
---|
| 128 | INTEGER invtau(-20:45000) ! ISCCP table for converting optical thickness |
---|
| 129 | c ! to count value |
---|
| 130 | ! |
---|
| 131 | ! The following input variables are used only if top_height = 1 or top_height = 3 |
---|
| 132 | ! |
---|
| 133 | REAL skt(npoints) ! skin Temperature (K) |
---|
| 134 | REAL emsfc_lw ! 10.5 micron emissivity of surface (fraction) |
---|
| 135 | REAL at(npoints,nlev) ! temperature in each model level (K) |
---|
| 136 | REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform |
---|
[1344] | 137 | c ! clouds in each |
---|
[524] | 138 | c ! model level. Same note applies as in dtau_s. |
---|
| 139 | REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective |
---|
[1344] | 140 | c ! clouds in each |
---|
[524] | 141 | c ! model level. Same note applies as in dtau_s. |
---|
| 142 | cIM reg.dyn BEG |
---|
| 143 | REAL t1, t2 |
---|
| 144 | ! REAL w(npoints) !vertical wind at 500 hPa |
---|
| 145 | ! LOGICAL pct_ocean(npoints) !TRUE if oceanic point, FALSE otherway |
---|
| 146 | ! INTEGER iw(npoints) , nw |
---|
| 147 | ! REAL wmin, pas_w |
---|
| 148 | ! INTEGER k, l, iwmx |
---|
| 149 | ! PARAMETER(wmin=-100.,pas_w=10.,iwmx=30) |
---|
| 150 | ! REAL fq_dynreg(7,7,iwmx) |
---|
| 151 | ! REAL nfq_dynreg(7,7,iwmx) |
---|
| 152 | ! LOGICAL pctj(7,7,iwmx) |
---|
| 153 | cIM reg.dyn END |
---|
| 154 | ! ------ |
---|
| 155 | ! Output |
---|
| 156 | ! ------ |
---|
| 157 | |
---|
| 158 | REAL fq_isccp(npoints,7,7) ! the fraction of the model grid box covered by |
---|
| 159 | c ! each of the 49 ISCCP D level cloud types |
---|
| 160 | |
---|
| 161 | REAL totalcldarea(npoints) ! the fraction of model grid box columns |
---|
| 162 | c ! with cloud somewhere in them. This should |
---|
[1344] | 163 | c ! equal the sum over all entries of fq_isccp |
---|
| 164 | |
---|
| 165 | |
---|
[524] | 166 | c ! The following three means are averages over the cloudy areas only. If no |
---|
[1344] | 167 | c ! clouds are in grid box all three quantities should equal zero. |
---|
| 168 | |
---|
[524] | 169 | REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging |
---|
| 170 | c ! in cloud top pressure. |
---|
[1344] | 171 | |
---|
[524] | 172 | REAL meantaucld(npoints) ! mean optical thickness |
---|
| 173 | c ! linear averaging in albedo performed. |
---|
| 174 | |
---|
| 175 | REAL boxtau(npoints,ncol) ! optical thickness in each column |
---|
| 176 | |
---|
| 177 | REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column |
---|
[1344] | 178 | |
---|
| 179 | |
---|
[524] | 180 | ! |
---|
| 181 | ! ------ |
---|
| 182 | ! Working variables added when program updated to mimic Mark Webb's PV-Wave code |
---|
| 183 | ! ------ |
---|
| 184 | |
---|
| 185 | REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into |
---|
[1344] | 186 | c ! Equivalent of BOX in original version, but |
---|
| 187 | c ! indexed by column then row, rather than |
---|
| 188 | c ! by row then column |
---|
[524] | 189 | |
---|
| 190 | REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction) |
---|
| 191 | c ! with extra layer of zeroes on top |
---|
| 192 | c ! in this version this just contains the values input |
---|
| 193 | c ! from cc but with an extra level |
---|
| 194 | REAL cca(npoints,nlev) ! convective cloud cover in each model level (fraction) |
---|
| 195 | c ! from conv |
---|
| 196 | |
---|
| 197 | REAL threshold(npoints,ncol) ! pointer to position in gridbox |
---|
| 198 | REAL maxocc(npoints,ncol) ! Flag for max overlapped conv cld |
---|
| 199 | REAL maxosc(npoints,ncol) ! Flag for max overlapped strat cld |
---|
| 200 | |
---|
| 201 | REAL boxpos(npoints,ncol) ! ordered pointer to position in gridbox |
---|
| 202 | |
---|
| 203 | REAL threshold_min(npoints,ncol) ! minimum value to define range in with new threshold |
---|
| 204 | c ! is chosen |
---|
| 205 | |
---|
| 206 | REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave |
---|
[1344] | 207 | c ! emissivity in part of |
---|
| 208 | c ! gridbox under consideration |
---|
[524] | 209 | |
---|
[1344] | 210 | REAL ran(npoints) ! vector of random numbers |
---|
[524] | 211 | REAL ptrop(npoints) |
---|
| 212 | REAL attrop(npoints) |
---|
| 213 | REAL attropmin (npoints) |
---|
| 214 | REAL atmax(npoints) |
---|
| 215 | REAL atmin(npoints) |
---|
| 216 | REAL btcmin(npoints) |
---|
| 217 | REAL transmax(npoints) |
---|
| 218 | |
---|
| 219 | INTEGER i,j,ilev,ibox,itrop(npoints) |
---|
| 220 | INTEGER ipres(npoints) |
---|
| 221 | INTEGER itau(npoints),ilev2 |
---|
| 222 | INTEGER acc(nlev,ncol) |
---|
| 223 | INTEGER match(npoints,nlev-1) |
---|
| 224 | INTEGER nmatch(npoints) |
---|
| 225 | INTEGER levmatch(npoints,ncol) |
---|
| 226 | |
---|
| 227 | c !variables needed for water vapor continuum absorption |
---|
| 228 | real fluxtop_clrsky(npoints),trans_layers_above_clrsky(npoints) |
---|
| 229 | real taumin(npoints) |
---|
| 230 | real dem_wv(npoints,nlev), wtmair, wtmh20, Navo, grav, pstd, t0 |
---|
| 231 | real press(npoints), dpress(npoints), atmden(npoints) |
---|
| 232 | real rvh20(npoints), wk(npoints), rhoave(npoints) |
---|
| 233 | real rh20s(npoints), rfrgn(npoints) |
---|
| 234 | real tmpexp(npoints),tauwv(npoints) |
---|
| 235 | |
---|
| 236 | character*1 cchar(6),cchar_realtops(6) |
---|
| 237 | integer icycle |
---|
| 238 | REAL tau(npoints,ncol) |
---|
| 239 | LOGICAL box_cloudy(npoints,ncol) |
---|
| 240 | REAL tb(npoints,ncol) |
---|
| 241 | REAL ptop(npoints,ncol) |
---|
| 242 | REAL emcld(npoints,ncol) |
---|
| 243 | REAL fluxtop(npoints,ncol) |
---|
| 244 | REAL trans_layers_above(npoints,ncol) |
---|
| 245 | real isccp_taumin,fluxtopinit(npoints),tauir(npoints) |
---|
| 246 | real meanalbedocld(npoints) |
---|
| 247 | REAL albedocld(npoints,ncol) |
---|
| 248 | real boxarea |
---|
| 249 | integer debug ! set to non-zero value to print out inputs |
---|
[1344] | 250 | c ! with step debug |
---|
[524] | 251 | integer debugcol ! set to non-zero value to print out column |
---|
[1344] | 252 | c ! decomposition with step debugcol |
---|
[524] | 253 | |
---|
| 254 | integer index1(npoints),num1,jj |
---|
| 255 | real rec2p13,tauchk |
---|
| 256 | |
---|
| 257 | character*10 ftn09 |
---|
| 258 | |
---|
| 259 | DATA isccp_taumin / 0.3 / |
---|
| 260 | DATA cchar / ' ','-','1','+','I','+'/ |
---|
| 261 | DATA cchar_realtops / ' ',' ','1','1','I','I'/ |
---|
| 262 | |
---|
| 263 | tauchk = -1.*log(0.9999999) |
---|
| 264 | rec2p13=1./2.13 |
---|
| 265 | |
---|
| 266 | ncolprint=0 |
---|
| 267 | |
---|
| 268 | cIM |
---|
| 269 | c PRINT*,' isccp_cloud_types npoints=',npoints |
---|
| 270 | c |
---|
| 271 | c if ( debug.ne.0 ) then |
---|
| 272 | c j=1 |
---|
| 273 | c write(6,'(a10)') 'j=' |
---|
| 274 | c write(6,'(8I10)') j |
---|
| 275 | c write(6,'(a10)') 'debug=' |
---|
| 276 | c write(6,'(8I10)') debug |
---|
| 277 | c write(6,'(a10)') 'debugcol=' |
---|
| 278 | c write(6,'(8I10)') debugcol |
---|
| 279 | c write(6,'(a10)') 'npoints=' |
---|
| 280 | c write(6,'(8I10)') npoints |
---|
| 281 | c write(6,'(a10)') 'nlev=' |
---|
| 282 | c write(6,'(8I10)') nlev |
---|
| 283 | c write(6,'(a10)') 'ncol=' |
---|
| 284 | c write(6,'(8I10)') ncol |
---|
| 285 | c write(6,'(a10)') 'top_height=' |
---|
| 286 | c write(6,'(8I10)') top_height |
---|
| 287 | c write(6,'(a10)') 'overlap=' |
---|
| 288 | c write(6,'(8I10)') overlap |
---|
| 289 | c write(6,'(a10)') 'emsfc_lw=' |
---|
| 290 | c write(6,'(8f10.2)') emsfc_lw |
---|
| 291 | c write(6,'(a10)') 'tautab=' |
---|
| 292 | c write(6,'(8f10.2)') tautab |
---|
| 293 | c write(6,'(a10)') 'invtau(1:100)=' |
---|
| 294 | c write(6,'(8i10)') (invtau(i),i=1,100) |
---|
[1344] | 295 | c do j=1,npoints,debug |
---|
[524] | 296 | c write(6,'(a10)') 'j=' |
---|
| 297 | c write(6,'(8I10)') j |
---|
| 298 | c write(6,'(a10)') 'sunlit=' |
---|
| 299 | c write(6,'(8I10)') sunlit(j) |
---|
| 300 | c write(6,'(a10)') 'seed=' |
---|
| 301 | c write(6,'(8I10)') seed(j) |
---|
| 302 | c write(6,'(a10)') 'pfull=' |
---|
| 303 | c write(6,'(8f10.2)') (pfull(j,i),i=1,nlev) |
---|
| 304 | c write(6,'(a10)') 'phalf=' |
---|
| 305 | c write(6,'(8f10.2)') (phalf(j,i),i=1,nlev+1) |
---|
| 306 | c write(6,'(a10)') 'qv=' |
---|
| 307 | c write(6,'(8f10.3)') (qv(j,i),i=1,nlev) |
---|
| 308 | c write(6,'(a10)') 'cc=' |
---|
| 309 | c write(6,'(8f10.3)') (cc(j,i),i=1,nlev) |
---|
| 310 | c write(6,'(a10)') 'conv=' |
---|
| 311 | c write(6,'(8f10.2)') (conv(j,i),i=1,nlev) |
---|
| 312 | c write(6,'(a10)') 'dtau_s=' |
---|
| 313 | c write(6,'(8g12.5)') (dtau_s(j,i),i=1,nlev) |
---|
| 314 | c write(6,'(a10)') 'dtau_c=' |
---|
| 315 | c write(6,'(8f10.2)') (dtau_c(j,i),i=1,nlev) |
---|
| 316 | c write(6,'(a10)') 'skt=' |
---|
| 317 | c write(6,'(8f10.2)') skt(j) |
---|
| 318 | c write(6,'(a10)') 'at=' |
---|
| 319 | c write(6,'(8f10.2)') (at(j,i),i=1,nlev) |
---|
| 320 | c write(6,'(a10)') 'dem_s=' |
---|
| 321 | c write(6,'(8f10.3)') (dem_s(j,i),i=1,nlev) |
---|
| 322 | c write(6,'(a10)') 'dem_c=' |
---|
| 323 | c write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev) |
---|
[1344] | 324 | c enddo |
---|
[524] | 325 | c endif |
---|
| 326 | |
---|
| 327 | ! ---------------------------------------------------! |
---|
| 328 | |
---|
| 329 | ! assign 2d tca array using 1d input array cc |
---|
| 330 | |
---|
| 331 | do j=1,npoints |
---|
| 332 | tca(j,0)=0 |
---|
| 333 | enddo |
---|
| 334 | |
---|
| 335 | do ilev=1,nlev |
---|
| 336 | do j=1,npoints |
---|
| 337 | tca(j,ilev)=cc(j,ilev) |
---|
| 338 | enddo |
---|
| 339 | enddo |
---|
| 340 | |
---|
| 341 | ! assign 2d cca array using 1d input array conv |
---|
| 342 | |
---|
| 343 | do ilev=1,nlev |
---|
| 344 | cIM pas besoin do ibox=1,ncol |
---|
| 345 | do j=1,npoints |
---|
| 346 | cca(j,ilev)=conv(j,ilev) |
---|
| 347 | enddo |
---|
| 348 | cIM enddo |
---|
| 349 | enddo |
---|
| 350 | |
---|
| 351 | cIM |
---|
| 352 | ! do j=1, iwmx |
---|
| 353 | ! do l=1, 7 |
---|
| 354 | ! do k=1, 7 |
---|
| 355 | ! fq_dynreg(k,l,j) =0. |
---|
| 356 | ! nfq_dynreg(k,l,j) =0. |
---|
| 357 | ! enddo !k |
---|
| 358 | ! enddo !l |
---|
| 359 | ! enddo !j |
---|
| 360 | cIM |
---|
| 361 | cIM |
---|
| 362 | c if (ncolprint.ne.0) then |
---|
[1344] | 363 | c do j=1,npoints,1000 |
---|
[524] | 364 | c write(6,'(a10)') 'j=' |
---|
| 365 | c write(6,'(8I10)') j |
---|
| 366 | c write (6,'(a)') 'seed:' |
---|
| 367 | c write (6,'(I3.2)') seed(j) |
---|
| 368 | |
---|
| 369 | c write (6,'(a)') 'tca_pp_rev:' |
---|
| 370 | c write (6,'(8f5.2)') |
---|
| 371 | c & ((tca(j,ilev)), |
---|
| 372 | c & ilev=1,nlev) |
---|
| 373 | |
---|
| 374 | c write (6,'(a)') 'cca_pp_rev:' |
---|
| 375 | c write (6,'(8f5.2)') |
---|
| 376 | c & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev) |
---|
[1344] | 377 | c enddo |
---|
[524] | 378 | c endif |
---|
| 379 | |
---|
| 380 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
| 381 | |
---|
| 382 | do j=1,npoints |
---|
| 383 | ptrop(j)=5000. |
---|
| 384 | atmin(j) = 400. |
---|
| 385 | attropmin(j) = 400. |
---|
| 386 | atmax(j) = 0. |
---|
| 387 | attrop(j) = 120. |
---|
| 388 | itrop(j) = 1 |
---|
| 389 | enddo |
---|
| 390 | |
---|
| 391 | do 12 ilev=1,nlev |
---|
| 392 | do j=1,npoints |
---|
| 393 | if (pfull(j,ilev) .lt. 40000. .and. |
---|
| 394 | & pfull(j,ilev) .gt. 5000. .and. |
---|
| 395 | & at(j,ilev) .lt. attropmin(j)) then |
---|
| 396 | ptrop(j) = pfull(j,ilev) |
---|
| 397 | attropmin(j) = at(j,ilev) |
---|
| 398 | attrop(j) = attropmin(j) |
---|
| 399 | itrop(j)=ilev |
---|
| 400 | end if |
---|
| 401 | if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev) |
---|
| 402 | if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev) |
---|
| 403 | enddo |
---|
| 404 | 12 continue |
---|
| 405 | |
---|
| 406 | end if |
---|
| 407 | |
---|
| 408 | ! -----------------------------------------------------! |
---|
| 409 | |
---|
| 410 | ! ---------------------------------------------------! |
---|
| 411 | |
---|
| 412 | cIM |
---|
| 413 | c do 13 ilev=1,nlev |
---|
| 414 | cnum1=0 |
---|
| 415 | c do j=1,npoints |
---|
| 416 | c if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then |
---|
[1344] | 417 | c num1=num1+1 |
---|
| 418 | c index1(num1)=j |
---|
[524] | 419 | c end if |
---|
| 420 | c enddo |
---|
| 421 | c do jj=1,num1 |
---|
[1344] | 422 | c j=index1(jj) |
---|
[524] | 423 | c write(6,*) ' error = cloud fraction less than zero' |
---|
[1344] | 424 | c write(6,*) ' or ' |
---|
[524] | 425 | c write(6,*) ' error = cloud fraction greater than 1' |
---|
[1344] | 426 | c write(6,*) 'value at point ',j,' is ',cc(j,ilev) |
---|
| 427 | c write(6,*) 'level ',ilev |
---|
[524] | 428 | c STOP |
---|
| 429 | c enddo |
---|
| 430 | cnum1=0 |
---|
| 431 | c do j=1,npoints |
---|
| 432 | c if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then |
---|
[1344] | 433 | c num1=num1+1 |
---|
| 434 | c index1(num1)=j |
---|
[524] | 435 | c end if |
---|
| 436 | c enddo |
---|
| 437 | c do jj=1,num1 |
---|
[1344] | 438 | c j=index1(jj) |
---|
[524] | 439 | c write(6,*) |
---|
| 440 | c & ' error = convective cloud fraction less than zero' |
---|
[1344] | 441 | c write(6,*) ' or ' |
---|
[524] | 442 | c write(6,*) |
---|
| 443 | c & ' error = convective cloud fraction greater than 1' |
---|
[1344] | 444 | c write(6,*) 'value at point ',j,' is ',conv(j,ilev) |
---|
| 445 | c write(6,*) 'level ',ilev |
---|
[524] | 446 | c STOP |
---|
| 447 | c enddo |
---|
| 448 | |
---|
| 449 | cnum1=0 |
---|
| 450 | c do j=1,npoints |
---|
| 451 | c if (dtau_s(j,ilev) .lt. 0.) then |
---|
[1344] | 452 | c num1=num1+1 |
---|
| 453 | c index1(num1)=j |
---|
[524] | 454 | c end if |
---|
| 455 | c enddo |
---|
| 456 | c do jj=1,num1 |
---|
[1344] | 457 | c j=index1(jj) |
---|
[524] | 458 | c write(6,*) |
---|
| 459 | c & ' error = stratiform cloud opt. depth less than zero' |
---|
[1344] | 460 | c write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) |
---|
| 461 | c write(6,*) 'level ',ilev |
---|
[524] | 462 | c STOP |
---|
| 463 | c enddo |
---|
| 464 | cnum1=0 |
---|
| 465 | c do j=1,npoints |
---|
| 466 | c if (dtau_c(j,ilev) .lt. 0.) then |
---|
[1344] | 467 | c num1=num1+1 |
---|
| 468 | c index1(num1)=j |
---|
[524] | 469 | c end if |
---|
| 470 | c enddo |
---|
| 471 | c do jj=1,num1 |
---|
[1344] | 472 | c j=index1(jj) |
---|
[524] | 473 | c write(6,*) |
---|
| 474 | c & ' error = convective cloud opt. depth less than zero' |
---|
[1344] | 475 | c write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) |
---|
| 476 | c write(6,*) 'level ',ilev |
---|
[524] | 477 | c STOP |
---|
| 478 | c enddo |
---|
| 479 | |
---|
| 480 | cnum1=0 |
---|
| 481 | c do j=1,npoints |
---|
| 482 | c if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then |
---|
[1344] | 483 | c num1=num1+1 |
---|
| 484 | c index1(num1)=j |
---|
[524] | 485 | c end if |
---|
| 486 | c enddo |
---|
| 487 | c do jj=1,num1 |
---|
[1344] | 488 | c j=index1(jj) |
---|
[524] | 489 | c write(6,*) |
---|
| 490 | c & ' error = stratiform cloud emissivity less than zero' |
---|
[1344] | 491 | c write(6,*)'or' |
---|
[524] | 492 | c write(6,*) |
---|
| 493 | c & ' error = stratiform cloud emissivity greater than 1' |
---|
[1344] | 494 | c write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) |
---|
| 495 | c write(6,*) 'level ',ilev |
---|
[524] | 496 | c STOP |
---|
| 497 | c enddo |
---|
| 498 | |
---|
| 499 | cnum1=0 |
---|
| 500 | c do j=1,npoints |
---|
| 501 | c if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then |
---|
[1344] | 502 | c num1=num1+1 |
---|
| 503 | c index1(num1)=j |
---|
[524] | 504 | c end if |
---|
| 505 | c enddo |
---|
| 506 | c do jj=1,num1 |
---|
[1344] | 507 | c j=index1(jj) |
---|
[524] | 508 | c write(6,*) |
---|
| 509 | c & ' error = convective cloud emissivity less than zero' |
---|
[1344] | 510 | c write(6,*)'or' |
---|
[524] | 511 | c write(6,*) |
---|
| 512 | c & ' error = convective cloud emissivity greater than 1' |
---|
| 513 | c write (6,*) |
---|
| 514 | c & 'j=',j,'ilev=',ilev,'dem_c(j,ilev) =',dem_c(j,ilev) |
---|
| 515 | c STOP |
---|
| 516 | c enddo |
---|
| 517 | c13 continue |
---|
| 518 | |
---|
| 519 | |
---|
| 520 | do ibox=1,ncol |
---|
| 521 | do j=1,npoints |
---|
[1344] | 522 | boxpos(j,ibox)=(ibox-.5)/ncol |
---|
[524] | 523 | enddo |
---|
| 524 | enddo |
---|
| 525 | |
---|
| 526 | ! ---------------------------------------------------! |
---|
| 527 | ! Initialise working variables |
---|
| 528 | ! ---------------------------------------------------! |
---|
| 529 | |
---|
| 530 | ! Initialised frac_out to zero |
---|
| 531 | |
---|
[1146] | 532 | do ilev=1,nlev |
---|
| 533 | do ibox=1,ncol |
---|
[524] | 534 | do j=1,npoints |
---|
[1344] | 535 | frac_out(j,ibox,ilev)=0.0 |
---|
[524] | 536 | enddo |
---|
| 537 | enddo |
---|
| 538 | enddo |
---|
| 539 | |
---|
| 540 | cIM |
---|
| 541 | c if (ncolprint.ne.0) then |
---|
| 542 | c write (6,'(a)') 'frac_out_pp_rev:' |
---|
| 543 | c do j=1,npoints,1000 |
---|
| 544 | c write(6,'(a10)') 'j=' |
---|
| 545 | c write(6,'(8I10)') j |
---|
| 546 | c write (6,'(8f5.2)') |
---|
| 547 | c & ((frac_out(j,ibox,ilev),ibox=1,ncolprint),ilev=1,nlev) |
---|
| 548 | |
---|
| 549 | c enddo |
---|
| 550 | c write (6,'(a)') 'ncol:' |
---|
| 551 | c write (6,'(I3)') ncol |
---|
| 552 | c endif |
---|
| 553 | c if (ncolprint.ne.0) then |
---|
| 554 | c write (6,'(a)') 'last_frac_pp:' |
---|
| 555 | c do j=1,npoints,1000 |
---|
| 556 | c write(6,'(a10)') 'j=' |
---|
| 557 | c write(6,'(8I10)') j |
---|
| 558 | c write (6,'(8f5.2)') (tca(j,0)) |
---|
| 559 | c enddo |
---|
| 560 | c endif |
---|
| 561 | |
---|
| 562 | ! ---------------------------------------------------! |
---|
| 563 | ! ALLOCATE CLOUD INTO BOXES, FOR NCOLUMNS, NLEVELS |
---|
| 564 | ! frac_out is the array that contains the information |
---|
| 565 | ! where 0 is no cloud, 1 is a stratiform cloud and 2 is a |
---|
| 566 | ! convective cloud |
---|
| 567 | |
---|
| 568 | !loop over vertical levels |
---|
| 569 | DO 200 ilev = 1,nlev |
---|
| 570 | |
---|
| 571 | ! Initialise threshold |
---|
| 572 | |
---|
| 573 | IF (ilev.eq.1) then |
---|
[1344] | 574 | ! If max overlap |
---|
| 575 | IF (overlap.eq.1) then |
---|
| 576 | ! select pixels spread evenly |
---|
| 577 | ! across the gridbox |
---|
[524] | 578 | DO ibox=1,ncol |
---|
| 579 | do j=1,npoints |
---|
| 580 | threshold(j,ibox)=boxpos(j,ibox) |
---|
| 581 | enddo |
---|
| 582 | enddo |
---|
[1344] | 583 | ELSE |
---|
[524] | 584 | DO ibox=1,ncol |
---|
| 585 | call ran0_vec(npoints,seed,ran) |
---|
[1344] | 586 | ! select random pixels from the non-convective |
---|
| 587 | ! part the gridbox ( some will be converted into |
---|
| 588 | ! convective pixels below ) |
---|
[524] | 589 | do j=1,npoints |
---|
| 590 | threshold(j,ibox)= |
---|
| 591 | & cca(j,ilev)+(1-cca(j,ilev))*ran(j) |
---|
| 592 | enddo |
---|
| 593 | enddo |
---|
| 594 | ENDIF |
---|
| 595 | cIM |
---|
| 596 | c IF (ncolprint.ne.0) then |
---|
| 597 | c write (6,'(a)') 'threshold_nsf2:' |
---|
| 598 | c do j=1,npoints,1000 |
---|
| 599 | c write(6,'(a10)') 'j=' |
---|
| 600 | c write(6,'(8I10)') j |
---|
| 601 | c write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) |
---|
| 602 | c enddo |
---|
| 603 | c ENDIF |
---|
| 604 | ENDIF |
---|
| 605 | |
---|
| 606 | c IF (ncolprint.ne.0) then |
---|
| 607 | c write (6,'(a)') 'ilev:' |
---|
| 608 | c write (6,'(I2)') ilev |
---|
| 609 | c ENDIF |
---|
| 610 | |
---|
| 611 | DO ibox=1,ncol |
---|
| 612 | |
---|
| 613 | ! All versions |
---|
| 614 | do j=1,npoints |
---|
| 615 | if (boxpos(j,ibox).le.cca(j,ilev)) then |
---|
| 616 | cIM REAL maxocc(j,ibox) = 1 |
---|
| 617 | maxocc(j,ibox) = 1.0 |
---|
| 618 | else |
---|
| 619 | cIM REAL maxocc(j,ibox) = 0 |
---|
| 620 | maxocc(j,ibox) = 0.0 |
---|
| 621 | end if |
---|
| 622 | enddo |
---|
| 623 | |
---|
| 624 | ! Max overlap |
---|
| 625 | if (overlap.eq.1) then |
---|
| 626 | do j=1,npoints |
---|
| 627 | threshold_min(j,ibox)=cca(j,ilev) |
---|
| 628 | cIM REAL maxosc(j,ibox)=1 |
---|
| 629 | maxosc(j,ibox)=1.0 |
---|
| 630 | enddo |
---|
| 631 | endif |
---|
| 632 | |
---|
| 633 | ! Random overlap |
---|
| 634 | if (overlap.eq.2) then |
---|
| 635 | do j=1,npoints |
---|
| 636 | threshold_min(j,ibox)=cca(j,ilev) |
---|
| 637 | cIM REAL maxosc(j,ibox)=0 |
---|
| 638 | maxosc(j,ibox)=0.0 |
---|
| 639 | enddo |
---|
| 640 | endif |
---|
| 641 | |
---|
| 642 | ! Max/Random overlap |
---|
| 643 | if (overlap.eq.3) then |
---|
| 644 | do j=1,npoints |
---|
| 645 | threshold_min(j,ibox)=max(cca(j,ilev), |
---|
| 646 | & min(tca(j,ilev-1),tca(j,ilev))) |
---|
| 647 | if (threshold(j,ibox) |
---|
| 648 | & .lt.min(tca(j,ilev-1),tca(j,ilev)) |
---|
| 649 | & .and.(threshold(j,ibox).gt.cca(j,ilev))) then |
---|
| 650 | cIM REAL maxosc(j,ibox)= 1 |
---|
| 651 | maxosc(j,ibox)= 1.0 |
---|
| 652 | else |
---|
| 653 | cIM REAL maxosc(j,ibox)= 0 |
---|
| 654 | maxosc(j,ibox)= 0.0 |
---|
| 655 | end if |
---|
| 656 | enddo |
---|
| 657 | endif |
---|
| 658 | |
---|
| 659 | ! Reset threshold |
---|
| 660 | call ran0_vec(npoints,seed,ran) |
---|
[1344] | 661 | |
---|
[524] | 662 | do j=1,npoints |
---|
| 663 | threshold(j,ibox)= |
---|
| 664 | !if max overlapped conv cloud |
---|
| 665 | & maxocc(j,ibox) * ( |
---|
| 666 | & boxpos(j,ibox) |
---|
| 667 | & ) + |
---|
| 668 | !else |
---|
| 669 | & (1-maxocc(j,ibox)) * ( |
---|
| 670 | !if max overlapped strat cloud |
---|
| 671 | & (maxosc(j,ibox)) * ( |
---|
| 672 | !threshold=boxpos |
---|
| 673 | & threshold(j,ibox) |
---|
| 674 | & ) + |
---|
| 675 | !else |
---|
| 676 | & (1-maxosc(j,ibox)) * ( |
---|
| 677 | !threshold_min=random[thrmin,1] |
---|
| 678 | & threshold_min(j,ibox)+ |
---|
| 679 | & (1-threshold_min(j,ibox))*ran(j) |
---|
| 680 | & ) |
---|
| 681 | & ) |
---|
| 682 | enddo |
---|
| 683 | |
---|
| 684 | ENDDO ! ibox |
---|
| 685 | |
---|
| 686 | ! Fill frac_out with 1's where tca is greater than the threshold |
---|
| 687 | |
---|
| 688 | DO ibox=1,ncol |
---|
| 689 | do j=1,npoints |
---|
| 690 | if (tca(j,ilev).gt.threshold(j,ibox)) then |
---|
| 691 | cIM REAL frac_out(j,ibox,ilev)=1 |
---|
| 692 | frac_out(j,ibox,ilev)=1.0 |
---|
| 693 | else |
---|
| 694 | cIM REAL frac_out(j,ibox,ilev)=0 |
---|
| 695 | frac_out(j,ibox,ilev)=0.0 |
---|
| 696 | end if |
---|
| 697 | enddo |
---|
| 698 | ENDDO |
---|
| 699 | |
---|
[1344] | 700 | ! Code to partition boxes into startiform and convective parts |
---|
| 701 | ! goes here |
---|
[524] | 702 | |
---|
| 703 | DO ibox=1,ncol |
---|
| 704 | do j=1,npoints |
---|
| 705 | if (threshold(j,ibox).le.cca(j,ilev)) then |
---|
| 706 | ! = 2 IF threshold le cca(j) |
---|
| 707 | cIM REAL frac_out(j,ibox,ilev) = 2 |
---|
| 708 | frac_out(j,ibox,ilev) = 2.0 |
---|
| 709 | else |
---|
| 710 | ! = the same IF NOT threshold le cca(j) |
---|
| 711 | frac_out(j,ibox,ilev) = frac_out(j,ibox,ilev) |
---|
| 712 | end if |
---|
| 713 | enddo |
---|
| 714 | ENDDO |
---|
| 715 | |
---|
| 716 | ! Set last_frac to tca at this level, so as to be tca |
---|
| 717 | ! from last level next time round |
---|
| 718 | |
---|
| 719 | cIM |
---|
| 720 | c if (ncolprint.ne.0) then |
---|
| 721 | |
---|
| 722 | c do j=1,npoints ,1000 |
---|
| 723 | c write(6,'(a10)') 'j=' |
---|
| 724 | c write(6,'(8I10)') j |
---|
| 725 | c write (6,'(a)') 'last_frac:' |
---|
| 726 | c write (6,'(8f5.2)') (tca(j,ilev-1)) |
---|
| 727 | |
---|
| 728 | c write (6,'(a)') 'cca:' |
---|
| 729 | c write (6,'(8f5.2)') (cca(j,ilev),ibox=1,ncolprint) |
---|
| 730 | |
---|
| 731 | c write (6,'(a)') 'max_overlap_cc:' |
---|
| 732 | c write (6,'(8f5.2)') (maxocc(j,ibox),ibox=1,ncolprint) |
---|
| 733 | |
---|
| 734 | c write (6,'(a)') 'max_overlap_sc:' |
---|
| 735 | c write (6,'(8f5.2)') (maxosc(j,ibox),ibox=1,ncolprint) |
---|
| 736 | |
---|
| 737 | c write (6,'(a)') 'threshold_min_nsf2:' |
---|
| 738 | c write (6,'(8f5.2)') (threshold_min(j,ibox),ibox=1,ncolprint) |
---|
| 739 | |
---|
| 740 | c write (6,'(a)') 'threshold_nsf2:' |
---|
| 741 | c write (6,'(8f5.2)') (threshold(j,ibox),ibox=1,ncolprint) |
---|
| 742 | |
---|
| 743 | c write (6,'(a)') 'frac_out_pp_rev:' |
---|
| 744 | c write (6,'(8f5.2)') |
---|
| 745 | c & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev) |
---|
[1344] | 746 | c enddo |
---|
[524] | 747 | c endif |
---|
| 748 | |
---|
| 749 | 200 CONTINUE !loop over nlev |
---|
| 750 | |
---|
| 751 | ! |
---|
| 752 | ! ---------------------------------------------------! |
---|
| 753 | |
---|
| 754 | |
---|
| 755 | ! |
---|
| 756 | ! ---------------------------------------------------! |
---|
| 757 | ! COMPUTE CLOUD OPTICAL DEPTH FOR EACH COLUMN and |
---|
| 758 | ! put into vector tau |
---|
| 759 | |
---|
| 760 | !initialize tau and albedocld to zero |
---|
| 761 | do 15 ibox=1,ncol |
---|
| 762 | do j=1,npoints |
---|
| 763 | tau(j,ibox)=0. |
---|
[1344] | 764 | albedocld(j,ibox)=0. |
---|
| 765 | boxtau(j,ibox)=0. |
---|
| 766 | boxptop(j,ibox)=0. |
---|
| 767 | box_cloudy(j,ibox)=.false. |
---|
[524] | 768 | enddo |
---|
| 769 | 15 continue |
---|
| 770 | |
---|
| 771 | !compute total cloud optical depth for each column |
---|
| 772 | do ilev=1,nlev |
---|
| 773 | !increment tau for each of the boxes |
---|
| 774 | do ibox=1,ncol |
---|
| 775 | do j=1,npoints |
---|
| 776 | cIM REAL if (frac_out(j,ibox,ilev).eq.1) then |
---|
| 777 | if (frac_out(j,ibox,ilev).eq.1.0) then |
---|
| 778 | tau(j,ibox)=tau(j,ibox) |
---|
| 779 | & + dtau_s(j,ilev) |
---|
| 780 | endif |
---|
| 781 | cIM REAL if (frac_out(j,ibox,ilev).eq.2) then |
---|
| 782 | if (frac_out(j,ibox,ilev).eq.2.0) then |
---|
| 783 | tau(j,ibox)=tau(j,ibox) |
---|
| 784 | & + dtau_c(j,ilev) |
---|
| 785 | end if |
---|
| 786 | enddo |
---|
| 787 | enddo ! ibox |
---|
| 788 | enddo ! ilev |
---|
| 789 | cIM |
---|
| 790 | c if (ncolprint.ne.0) then |
---|
| 791 | |
---|
| 792 | c do j=1,npoints ,1000 |
---|
| 793 | c write(6,'(a10)') 'j=' |
---|
| 794 | c write(6,'(8I10)') j |
---|
| 795 | c write(6,'(i2,1X,8(f7.2,1X))') |
---|
| 796 | c & ilev, |
---|
| 797 | c & (tau(j,ibox),ibox=1,ncolprint) |
---|
| 798 | c enddo |
---|
| 799 | c endif |
---|
| 800 | ! |
---|
| 801 | ! ---------------------------------------------------! |
---|
| 802 | |
---|
| 803 | |
---|
| 804 | |
---|
| 805 | ! |
---|
| 806 | ! ---------------------------------------------------! |
---|
| 807 | ! COMPUTE INFRARED BRIGHTNESS TEMPERUATRES |
---|
| 808 | ! AND CLOUD TOP TEMPERATURE SATELLITE SHOULD SEE |
---|
| 809 | ! |
---|
| 810 | ! again this is only done if top_height = 1 or 3 |
---|
| 811 | ! |
---|
| 812 | ! fluxtop is the 10.5 micron radiance at the top of the |
---|
| 813 | ! atmosphere |
---|
| 814 | ! trans_layers_above is the total transmissivity in the layers |
---|
| 815 | ! above the current layer |
---|
| 816 | ! fluxtop_clrsky(j) and trans_layers_above_clrsky(j) are the clear |
---|
| 817 | ! sky versions of these quantities. |
---|
| 818 | |
---|
| 819 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
| 820 | |
---|
| 821 | |
---|
| 822 | !---------------------------------------------------------------------- |
---|
| 823 | ! |
---|
| 824 | ! DO CLEAR SKY RADIANCE CALCULATION FIRST |
---|
| 825 | ! |
---|
| 826 | !compute water vapor continuum emissivity |
---|
| 827 | !this treatment follows Schwarkzopf and Ramasamy |
---|
| 828 | !JGR 1999,vol 104, pages 9467-9499. |
---|
| 829 | !the emissivity is calculated at a wavenumber of 955 cm-1, |
---|
| 830 | !or 10.47 microns |
---|
| 831 | wtmair = 28.9644 |
---|
| 832 | wtmh20 = 18.01534 |
---|
| 833 | Navo = 6.023E+23 |
---|
| 834 | grav = 9.806650E+02 |
---|
| 835 | pstd = 1.013250E+06 |
---|
| 836 | t0 = 296. |
---|
| 837 | cIM |
---|
| 838 | c if (ncolprint .ne. 0) |
---|
| 839 | c & write(6,*) 'ilev pw (kg/m2) tauwv(j) dem_wv' |
---|
| 840 | do 125 ilev=1,nlev |
---|
| 841 | do j=1,npoints |
---|
| 842 | !press and dpress are dyne/cm2 = Pascals *10 |
---|
| 843 | press(j) = pfull(j,ilev)*10. |
---|
| 844 | dpress(j) = (phalf(j,ilev+1)-phalf(j,ilev))*10 |
---|
| 845 | !atmden = g/cm2 = kg/m2 / 10 |
---|
| 846 | atmden(j) = dpress(j)/grav |
---|
| 847 | rvh20(j) = qv(j,ilev)*wtmair/wtmh20 |
---|
| 848 | wk(j) = rvh20(j)*Navo*atmden(j)/wtmair |
---|
| 849 | rhoave(j) = (press(j)/pstd)*(t0/at(j,ilev)) |
---|
| 850 | rh20s(j) = rvh20(j)*rhoave(j) |
---|
| 851 | rfrgn(j) = rhoave(j)-rh20s(j) |
---|
| 852 | tmpexp(j) = exp(-0.02*(at(j,ilev)-t0)) |
---|
| 853 | tauwv(j) = wk(j)*1.e-20*( |
---|
| 854 | & (0.0224697*rh20s(j)*tmpexp(j)) + |
---|
| 855 | & (3.41817e-7*rfrgn(j)) )*0.98 |
---|
| 856 | dem_wv(j,ilev) = 1. - exp( -1. * tauwv(j)) |
---|
| 857 | enddo |
---|
| 858 | cIM |
---|
| 859 | c if (ncolprint .ne. 0) then |
---|
| 860 | c do j=1,npoints ,1000 |
---|
| 861 | c write(6,'(a10)') 'j=' |
---|
| 862 | c write(6,'(8I10)') j |
---|
| 863 | c write(6,'(i2,1X,3(f8.3,3X))') ilev, |
---|
| 864 | c & qv(j,ilev)*(phalf(j,ilev+1)-phalf(j,ilev))/(grav/100.), |
---|
| 865 | c & tauwv(j),dem_wv(j,ilev) |
---|
| 866 | c enddo |
---|
[1344] | 867 | c endif |
---|
[524] | 868 | 125 continue |
---|
| 869 | |
---|
| 870 | !initialize variables |
---|
| 871 | do j=1,npoints |
---|
| 872 | fluxtop_clrsky(j) = 0. |
---|
| 873 | trans_layers_above_clrsky(j)=1. |
---|
| 874 | enddo |
---|
| 875 | |
---|
| 876 | do ilev=1,nlev |
---|
| 877 | do j=1,npoints |
---|
| 878 | |
---|
| 879 | ! Black body emission at temperature of the layer |
---|
| 880 | |
---|
[1344] | 881 | bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) |
---|
| 882 | !bb(j)= 5.67e-8*at(j,ilev)**4 |
---|
[524] | 883 | |
---|
[1344] | 884 | ! increase TOA flux by flux emitted from layer |
---|
| 885 | ! times total transmittance in layers above |
---|
[524] | 886 | |
---|
| 887 | fluxtop_clrsky(j) = fluxtop_clrsky(j) |
---|
| 888 | & + dem_wv(j,ilev)*bb(j)*trans_layers_above_clrsky(j) |
---|
| 889 | |
---|
| 890 | ! update trans_layers_above with transmissivity |
---|
[1344] | 891 | ! from this layer for next time around loop |
---|
[524] | 892 | |
---|
| 893 | trans_layers_above_clrsky(j)= |
---|
| 894 | & trans_layers_above_clrsky(j)*(1.-dem_wv(j,ilev)) |
---|
| 895 | |
---|
| 896 | |
---|
| 897 | enddo |
---|
| 898 | cIM |
---|
| 899 | c if (ncolprint.ne.0) then |
---|
| 900 | c do j=1,npoints ,1000 |
---|
| 901 | c write(6,'(a10)') 'j=' |
---|
| 902 | c write(6,'(8I10)') j |
---|
| 903 | c write (6,'(a)') 'ilev:' |
---|
| 904 | c write (6,'(I2)') ilev |
---|
| 905 | |
---|
| 906 | c write (6,'(a)') |
---|
| 907 | c & 'emiss_layer,100.*bb(j),100.*f,total_trans:' |
---|
| 908 | c write (6,'(4(f7.2,1X))') dem_wv(j,ilev),100.*bb(j), |
---|
| 909 | c & 100.*fluxtop_clrsky(j),trans_layers_above_clrsky(j) |
---|
| 910 | c enddo |
---|
| 911 | c endif |
---|
| 912 | |
---|
| 913 | enddo !loop over level |
---|
| 914 | |
---|
| 915 | do j=1,npoints |
---|
| 916 | !add in surface emission |
---|
| 917 | bb(j)=1/( exp(1307.27/skt(j)) - 1. ) |
---|
| 918 | !bb(j)=5.67e-8*skt(j)**4 |
---|
| 919 | |
---|
| 920 | fluxtop_clrsky(j) = fluxtop_clrsky(j) + emsfc_lw * bb(j) |
---|
| 921 | & * trans_layers_above_clrsky(j) |
---|
| 922 | enddo |
---|
| 923 | |
---|
| 924 | cIM |
---|
| 925 | c if (ncolprint.ne.0) then |
---|
| 926 | c do j=1,npoints ,1000 |
---|
| 927 | c write(6,'(a10)') 'j=' |
---|
| 928 | c write(6,'(8I10)') j |
---|
| 929 | c write (6,'(a)') 'id:' |
---|
| 930 | c write (6,'(a)') 'surface' |
---|
| 931 | |
---|
| 932 | c write (6,'(a)') 'emsfc,100.*bb(j),100.*f,total_trans:' |
---|
| 933 | c write (6,'(4(f7.2,1X))') emsfc_lw,100.*bb(j), |
---|
| 934 | c & 100.*fluxtop_clrsky(j), |
---|
| 935 | c & trans_layers_above_clrsky(j) |
---|
| 936 | c enddo |
---|
[1344] | 937 | c endif |
---|
[524] | 938 | |
---|
| 939 | |
---|
| 940 | ! |
---|
| 941 | ! END OF CLEAR SKY CALCULATION |
---|
| 942 | ! |
---|
| 943 | !---------------------------------------------------------------- |
---|
| 944 | |
---|
| 945 | |
---|
| 946 | cIM |
---|
| 947 | c if (ncolprint.ne.0) then |
---|
| 948 | |
---|
| 949 | c do j=1,npoints ,1000 |
---|
| 950 | c write(6,'(a10)') 'j=' |
---|
| 951 | c write(6,'(8I10)') j |
---|
| 952 | c write (6,'(a)') 'ts:' |
---|
| 953 | c write (6,'(8f7.2)') (skt(j),ibox=1,ncolprint) |
---|
| 954 | |
---|
| 955 | c write (6,'(a)') 'ta_rev:' |
---|
| 956 | c write (6,'(8f7.2)') |
---|
| 957 | c & ((at(j,ilev2),ibox=1,ncolprint),ilev2=1,nlev) |
---|
| 958 | |
---|
| 959 | c enddo |
---|
| 960 | c endif |
---|
| 961 | !loop over columns |
---|
| 962 | do ibox=1,ncol |
---|
| 963 | do j=1,npoints |
---|
| 964 | fluxtop(j,ibox)=0. |
---|
| 965 | trans_layers_above(j,ibox)=1. |
---|
| 966 | enddo |
---|
| 967 | enddo |
---|
| 968 | |
---|
| 969 | do ilev=1,nlev |
---|
| 970 | do j=1,npoints |
---|
| 971 | ! Black body emission at temperature of the layer |
---|
| 972 | |
---|
[1344] | 973 | bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) |
---|
| 974 | !bb(j)= 5.67e-8*at(j,ilev)**4 |
---|
[524] | 975 | enddo |
---|
| 976 | |
---|
| 977 | do ibox=1,ncol |
---|
| 978 | do j=1,npoints |
---|
| 979 | |
---|
[1344] | 980 | ! emissivity for point in this layer |
---|
[524] | 981 | cIM REAL if (frac_out(j,ibox,ilev).eq.1) then |
---|
| 982 | if (frac_out(j,ibox,ilev).eq.1.0) then |
---|
| 983 | dem(j,ibox)= 1. - |
---|
| 984 | & ( (1. - dem_wv(j,ilev)) * (1. - dem_s(j,ilev)) ) |
---|
| 985 | cIM REAL else if (frac_out(j,ibox,ilev).eq.2) then |
---|
| 986 | else if (frac_out(j,ibox,ilev).eq.2.0) then |
---|
| 987 | dem(j,ibox)= 1. - |
---|
| 988 | & ( (1. - dem_wv(j,ilev)) * (1. - dem_c(j,ilev)) ) |
---|
| 989 | else |
---|
| 990 | dem(j,ibox)= dem_wv(j,ilev) |
---|
| 991 | end if |
---|
| 992 | |
---|
| 993 | |
---|
| 994 | ! increase TOA flux by flux emitted from layer |
---|
[1344] | 995 | ! times total transmittance in layers above |
---|
[524] | 996 | |
---|
| 997 | fluxtop(j,ibox) = fluxtop(j,ibox) |
---|
| 998 | & + dem(j,ibox) * bb(j) |
---|
| 999 | & * trans_layers_above(j,ibox) |
---|
| 1000 | |
---|
| 1001 | ! update trans_layers_above with transmissivity |
---|
[1344] | 1002 | ! from this layer for next time around loop |
---|
[524] | 1003 | |
---|
| 1004 | trans_layers_above(j,ibox)= |
---|
| 1005 | & trans_layers_above(j,ibox)*(1.-dem(j,ibox)) |
---|
| 1006 | |
---|
| 1007 | enddo ! j |
---|
| 1008 | enddo ! ibox |
---|
| 1009 | |
---|
| 1010 | cIM |
---|
| 1011 | c if (ncolprint.ne.0) then |
---|
| 1012 | c do j=1,npoints,1000 |
---|
| 1013 | c write (6,'(a)') 'ilev:' |
---|
| 1014 | c write (6,'(I2)') ilev |
---|
| 1015 | |
---|
| 1016 | c write(6,'(a10)') 'j=' |
---|
| 1017 | c write(6,'(8I10)') j |
---|
| 1018 | c write (6,'(a)') 'emiss_layer:' |
---|
| 1019 | c write (6,'(8f7.2)') (dem(j,ibox),ibox=1,ncolprint) |
---|
| 1020 | |
---|
| 1021 | c write (6,'(a)') '100.*bb(j):' |
---|
| 1022 | c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) |
---|
| 1023 | |
---|
| 1024 | c write (6,'(a)') '100.*f:' |
---|
| 1025 | c write (6,'(8f7.2)') |
---|
| 1026 | c & (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
| 1027 | |
---|
| 1028 | c write (6,'(a)') 'total_trans:' |
---|
| 1029 | c write (6,'(8f7.2)') |
---|
| 1030 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
[1344] | 1031 | c enddo |
---|
[524] | 1032 | c endif |
---|
| 1033 | |
---|
| 1034 | enddo ! ilev |
---|
| 1035 | |
---|
| 1036 | |
---|
| 1037 | do j=1,npoints |
---|
| 1038 | !add in surface emission |
---|
| 1039 | bb(j)=1/( exp(1307.27/skt(j)) - 1. ) |
---|
| 1040 | !bb(j)=5.67e-8*skt(j)**4 |
---|
| 1041 | end do |
---|
| 1042 | |
---|
| 1043 | do ibox=1,ncol |
---|
| 1044 | do j=1,npoints |
---|
| 1045 | |
---|
| 1046 | !add in surface emission |
---|
| 1047 | |
---|
| 1048 | fluxtop(j,ibox) = fluxtop(j,ibox) |
---|
| 1049 | & + emsfc_lw * bb(j) |
---|
| 1050 | & * trans_layers_above(j,ibox) |
---|
| 1051 | |
---|
| 1052 | end do |
---|
| 1053 | end do |
---|
| 1054 | |
---|
| 1055 | cIM |
---|
| 1056 | c if (ncolprint.ne.0) then |
---|
| 1057 | |
---|
| 1058 | c do j=1,npoints ,1000 |
---|
| 1059 | c write(6,'(a10)') 'j=' |
---|
| 1060 | c write(6,'(8I10)') j |
---|
| 1061 | c write (6,'(a)') 'id:' |
---|
| 1062 | c write (6,'(a)') 'surface' |
---|
| 1063 | |
---|
| 1064 | c write (6,'(a)') 'emiss_layer:' |
---|
| 1065 | c write (6,'(8f7.2)') (dem(1,ibox),ibox=1,ncolprint) |
---|
| 1066 | |
---|
| 1067 | c write (6,'(a)') '100.*bb(j):' |
---|
| 1068 | c write (6,'(8f7.2)') (100.*bb(j),ibox=1,ncolprint) |
---|
| 1069 | |
---|
| 1070 | c write (6,'(a)') '100.*f:' |
---|
| 1071 | c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
| 1072 | c end do |
---|
[1344] | 1073 | c endif |
---|
[524] | 1074 | |
---|
| 1075 | !now that you have the top of atmosphere radiance account |
---|
| 1076 | !for ISCCP procedures to determine cloud top temperature |
---|
| 1077 | |
---|
| 1078 | !account for partially transmitting cloud recompute flux |
---|
| 1079 | !ISCCP would see assuming a single layer cloud |
---|
| 1080 | !note choice here of 2.13, as it is primarily ice |
---|
| 1081 | !clouds which have partial emissivity and need the |
---|
| 1082 | !adjustment performed in this section |
---|
| 1083 | ! |
---|
[1344] | 1084 | !If it turns out that the cloud brightness temperature |
---|
| 1085 | !is greater than 260K, then the liquid cloud conversion |
---|
[524] | 1086 | !factor of 2.56 is used. |
---|
[1344] | 1087 | ! |
---|
[524] | 1088 | !Note that this is discussed on pages 85-87 of |
---|
| 1089 | !the ISCCP D level documentation (Rossow et al. 1996) |
---|
| 1090 | |
---|
| 1091 | do j=1,npoints |
---|
| 1092 | !compute minimum brightness temperature and optical depth |
---|
| 1093 | btcmin(j) = 1. / ( exp(1307.27/(attrop(j)-5.)) - 1. ) |
---|
| 1094 | enddo |
---|
| 1095 | do ibox=1,ncol |
---|
| 1096 | do j=1,npoints |
---|
| 1097 | transmax(j) = (fluxtop(j,ibox)-btcmin(j)) |
---|
| 1098 | & /(fluxtop_clrsky(j)-btcmin(j)) |
---|
[1344] | 1099 | !note that the initial setting of tauir(j) is needed so that |
---|
| 1100 | !tauir(j) has a realistic value should the next if block be |
---|
| 1101 | !bypassed |
---|
[524] | 1102 | tauir(j) = tau(j,ibox) * rec2p13 |
---|
| 1103 | taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001)) |
---|
| 1104 | |
---|
| 1105 | enddo |
---|
| 1106 | |
---|
| 1107 | if (top_height .eq. 1) then |
---|
| 1108 | do j=1,npoints |
---|
| 1109 | if (transmax(j) .gt. 0.001 .and. |
---|
| 1110 | & transmax(j) .le. 0.9999999) then |
---|
| 1111 | fluxtopinit(j) = fluxtop(j,ibox) |
---|
[1344] | 1112 | tauir(j) = tau(j,ibox) *rec2p13 |
---|
[524] | 1113 | endif |
---|
| 1114 | enddo |
---|
| 1115 | do icycle=1,2 |
---|
| 1116 | do j=1,npoints |
---|
| 1117 | if (tau(j,ibox) .gt. (tauchk )) then |
---|
| 1118 | if (transmax(j) .gt. 0.001 .and. |
---|
| 1119 | & transmax(j) .le. 0.9999999) then |
---|
| 1120 | emcld(j,ibox) = 1. - exp(-1. * tauir(j) ) |
---|
| 1121 | fluxtop(j,ibox) = fluxtopinit(j) - |
---|
| 1122 | & ((1.-emcld(j,ibox))*fluxtop_clrsky(j)) |
---|
| 1123 | fluxtop(j,ibox)=max(1.E-06, |
---|
| 1124 | & (fluxtop(j,ibox)/emcld(j,ibox))) |
---|
| 1125 | tb(j,ibox)= 1307.27 |
---|
| 1126 | & / (log(1. + (1./fluxtop(j,ibox)))) |
---|
| 1127 | if (tb(j,ibox) .gt. 260.) then |
---|
[1344] | 1128 | tauir(j) = tau(j,ibox) / 2.56 |
---|
| 1129 | end if |
---|
[524] | 1130 | end if |
---|
| 1131 | end if |
---|
| 1132 | enddo |
---|
| 1133 | enddo |
---|
| 1134 | |
---|
| 1135 | endif |
---|
| 1136 | |
---|
| 1137 | do j=1,npoints |
---|
| 1138 | if (tau(j,ibox) .gt. (tauchk )) then |
---|
| 1139 | !cloudy box |
---|
| 1140 | tb(j,ibox)= 1307.27/ (log(1. + (1./fluxtop(j,ibox)))) |
---|
| 1141 | if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then |
---|
| 1142 | tb(j,ibox) = attrop(j) - 5. |
---|
[1344] | 1143 | tau(j,ibox) = 2.13*taumin(j) |
---|
[524] | 1144 | end if |
---|
| 1145 | else |
---|
| 1146 | !clear sky brightness temperature |
---|
| 1147 | tb(j,ibox) = 1307.27/(log(1.+(1./fluxtop_clrsky(j)))) |
---|
| 1148 | end if |
---|
| 1149 | enddo ! j |
---|
| 1150 | enddo ! ibox |
---|
| 1151 | |
---|
| 1152 | cIM |
---|
| 1153 | c if (ncolprint.ne.0) then |
---|
| 1154 | |
---|
| 1155 | c do j=1,npoints,1000 |
---|
| 1156 | c write(6,'(a10)') 'j=' |
---|
| 1157 | c write(6,'(8I10)') j |
---|
| 1158 | |
---|
| 1159 | c write (6,'(a)') 'attrop:' |
---|
| 1160 | c write (6,'(8f7.2)') (attrop(j)) |
---|
| 1161 | |
---|
| 1162 | c write (6,'(a)') 'btcmin:' |
---|
| 1163 | c write (6,'(8f7.2)') (btcmin(j)) |
---|
| 1164 | |
---|
| 1165 | c write (6,'(a)') 'fluxtop_clrsky*100:' |
---|
| 1166 | c write (6,'(8f7.2)') |
---|
| 1167 | c & (100.*fluxtop_clrsky(j)) |
---|
| 1168 | |
---|
| 1169 | c write (6,'(a)') '100.*f_adj:' |
---|
| 1170 | c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) |
---|
| 1171 | |
---|
| 1172 | c write (6,'(a)') 'transmax:' |
---|
| 1173 | c write (6,'(8f7.2)') (transmax(ibox),ibox=1,ncolprint) |
---|
| 1174 | |
---|
| 1175 | c write (6,'(a)') 'tau:' |
---|
| 1176 | c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) |
---|
| 1177 | |
---|
| 1178 | c write (6,'(a)') 'emcld:' |
---|
| 1179 | c write (6,'(8f7.2)') (emcld(j,ibox),ibox=1,ncolprint) |
---|
| 1180 | |
---|
| 1181 | c write (6,'(a)') 'total_trans:' |
---|
| 1182 | c write (6,'(8f7.2)') |
---|
[1344] | 1183 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
[524] | 1184 | |
---|
| 1185 | c write (6,'(a)') 'total_emiss:' |
---|
| 1186 | c write (6,'(8f7.2)') |
---|
[1344] | 1187 | c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
[524] | 1188 | |
---|
| 1189 | c write (6,'(a)') 'total_trans:' |
---|
| 1190 | c write (6,'(8f7.2)') |
---|
[1344] | 1191 | c & (trans_layers_above(j,ibox),ibox=1,ncolprint) |
---|
[524] | 1192 | |
---|
| 1193 | c write (6,'(a)') 'ppout:' |
---|
| 1194 | c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) |
---|
| 1195 | c enddo ! j |
---|
[1344] | 1196 | c endif |
---|
[524] | 1197 | |
---|
| 1198 | end if |
---|
| 1199 | |
---|
| 1200 | ! ---------------------------------------------------! |
---|
| 1201 | |
---|
| 1202 | ! |
---|
| 1203 | ! ---------------------------------------------------! |
---|
| 1204 | ! DETERMINE CLOUD TOP PRESSURE |
---|
| 1205 | ! |
---|
| 1206 | ! again the 2 methods differ according to whether |
---|
| 1207 | ! or not you use the physical cloud top pressure (top_height = 2) |
---|
| 1208 | ! or the radiatively determined cloud top pressure (top_height = 1 or 3) |
---|
| 1209 | ! |
---|
| 1210 | |
---|
| 1211 | !compute cloud top pressure |
---|
| 1212 | do 30 ibox=1,ncol |
---|
| 1213 | !segregate according to optical thickness |
---|
| 1214 | if (top_height .eq. 1 .or. top_height .eq. 3) then |
---|
| 1215 | !find level whose temperature |
---|
| 1216 | !most closely matches brightness temperature |
---|
| 1217 | do j=1,npoints |
---|
| 1218 | nmatch(j)=0 |
---|
| 1219 | enddo |
---|
| 1220 | do 29 ilev=1,nlev-1 |
---|
[1146] | 1221 | !cdir nodep |
---|
[524] | 1222 | do j=1,npoints |
---|
| 1223 | if ((at(j,ilev) .ge. tb(j,ibox) .and. |
---|
| 1224 | & at(j,ilev+1) .lt. tb(j,ibox)) .or. |
---|
| 1225 | & (at(j,ilev) .le. tb(j,ibox) .and. |
---|
| 1226 | & at(j,ilev+1) .gt. tb(j,ibox))) then |
---|
| 1227 | |
---|
| 1228 | nmatch(j)=nmatch(j)+1 |
---|
| 1229 | if(abs(at(j,ilev)-tb(j,ibox)) .lt. |
---|
| 1230 | & abs(at(j,ilev+1)-tb(j,ibox))) then |
---|
| 1231 | match(j,nmatch(j))=ilev |
---|
| 1232 | else |
---|
| 1233 | match(j,nmatch(j))=ilev+1 |
---|
| 1234 | end if |
---|
| 1235 | end if |
---|
| 1236 | enddo |
---|
| 1237 | 29 continue |
---|
| 1238 | |
---|
| 1239 | do j=1,npoints |
---|
| 1240 | if (nmatch(j) .ge. 1) then |
---|
| 1241 | ptop(j,ibox)=pfull(j,match(j,nmatch(j))) |
---|
| 1242 | levmatch(j,ibox)=match(j,nmatch(j)) |
---|
| 1243 | else |
---|
| 1244 | if (tb(j,ibox) .lt. atmin(j)) then |
---|
| 1245 | ptop(j,ibox)=ptrop(j) |
---|
| 1246 | levmatch(j,ibox)=itrop(j) |
---|
| 1247 | end if |
---|
| 1248 | if (tb(j,ibox) .gt. atmax(j)) then |
---|
| 1249 | ptop(j,ibox)=pfull(j,nlev) |
---|
| 1250 | levmatch(j,ibox)=nlev |
---|
| 1251 | end if |
---|
| 1252 | end if |
---|
| 1253 | enddo ! j |
---|
| 1254 | |
---|
| 1255 | else ! if (top_height .eq. 1 .or. top_height .eq. 3) |
---|
| 1256 | |
---|
| 1257 | do j=1,npoints |
---|
| 1258 | ptop(j,ibox)=0. |
---|
| 1259 | enddo |
---|
| 1260 | do ilev=1,nlev |
---|
| 1261 | do j=1,npoints |
---|
| 1262 | if ((ptop(j,ibox) .eq. 0. ) |
---|
| 1263 | cIM & .and.(frac_out(j,ibox,ilev) .ne. 0)) then |
---|
| 1264 | & .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then |
---|
| 1265 | ptop(j,ibox)=pfull(j,ilev) |
---|
[1344] | 1266 | levmatch(j,ibox)=ilev |
---|
[524] | 1267 | end if |
---|
| 1268 | end do |
---|
| 1269 | end do |
---|
| 1270 | end if |
---|
| 1271 | |
---|
| 1272 | do j=1,npoints |
---|
| 1273 | if (tau(j,ibox) .le. (tauchk )) then |
---|
| 1274 | ptop(j,ibox)=0. |
---|
| 1275 | levmatch(j,ibox)=0 |
---|
| 1276 | endif |
---|
| 1277 | enddo |
---|
| 1278 | |
---|
| 1279 | 30 continue |
---|
| 1280 | |
---|
| 1281 | ! |
---|
| 1282 | ! |
---|
| 1283 | ! ---------------------------------------------------! |
---|
| 1284 | |
---|
| 1285 | |
---|
| 1286 | ! |
---|
| 1287 | ! ---------------------------------------------------! |
---|
| 1288 | ! DETERMINE ISCCP CLOUD TYPE FREQUENCIES |
---|
| 1289 | ! |
---|
| 1290 | ! Now that ptop and tau have been determined, |
---|
| 1291 | ! determine amount of each of the 49 ISCCP cloud |
---|
| 1292 | ! types |
---|
| 1293 | ! |
---|
| 1294 | ! Also compute grid box mean cloud top pressure and |
---|
| 1295 | ! optical thickness. The mean cloud top pressure and |
---|
| 1296 | ! optical thickness are averages over the cloudy |
---|
| 1297 | ! area only. The mean cloud top pressure is a linear |
---|
| 1298 | ! average of the cloud top pressures. The mean cloud |
---|
| 1299 | ! optical thickness is computed by converting optical |
---|
| 1300 | ! thickness to an albedo, averaging in albedo units, |
---|
| 1301 | ! then converting the average albedo back to a mean |
---|
| 1302 | ! optical thickness. |
---|
| 1303 | ! |
---|
| 1304 | |
---|
| 1305 | !compute isccp frequencies |
---|
| 1306 | |
---|
| 1307 | !reset frequencies |
---|
| 1308 | do 38 ilev=1,7 |
---|
| 1309 | do 38 ilev2=1,7 |
---|
| 1310 | do j=1,npoints ! |
---|
| 1311 | fq_isccp(j,ilev,ilev2)=0. |
---|
| 1312 | enddo |
---|
| 1313 | 38 continue |
---|
| 1314 | |
---|
| 1315 | !reset variables need for averaging cloud properties |
---|
| 1316 | do j=1,npoints |
---|
| 1317 | totalcldarea(j) = 0. |
---|
| 1318 | meanalbedocld(j) = 0. |
---|
| 1319 | meanptop(j) = 0. |
---|
| 1320 | meantaucld(j) = 0. |
---|
| 1321 | enddo ! j |
---|
| 1322 | |
---|
| 1323 | boxarea = 1./real(ncol) |
---|
| 1324 | |
---|
| 1325 | !determine optical depth category |
---|
| 1326 | cIM do 39 j=1,npoints |
---|
| 1327 | cIM do ibox=1,ncol |
---|
| 1328 | do 39 ibox=1,ncol |
---|
| 1329 | do j=1,npoints |
---|
| 1330 | |
---|
| 1331 | cIM |
---|
| 1332 | c CALL CPU_time(t1) |
---|
| 1333 | cIM |
---|
| 1334 | |
---|
| 1335 | if (tau(j,ibox) .gt. (tauchk ) |
---|
| 1336 | & .and. ptop(j,ibox) .gt. 0.) then |
---|
| 1337 | box_cloudy(j,ibox)=.true. |
---|
| 1338 | endif |
---|
| 1339 | |
---|
| 1340 | cIM |
---|
| 1341 | c CALL CPU_time(t2) |
---|
| 1342 | c print*,'IF tau t2 - t1',t2 - t1 |
---|
| 1343 | |
---|
| 1344 | c CALL CPU_time(t1) |
---|
| 1345 | cIM |
---|
| 1346 | |
---|
| 1347 | if (box_cloudy(j,ibox)) then |
---|
| 1348 | |
---|
| 1349 | ! totalcldarea always diagnosed day or night |
---|
| 1350 | totalcldarea(j) = totalcldarea(j) + boxarea |
---|
| 1351 | |
---|
| 1352 | if (sunlit(j).eq.1) then |
---|
| 1353 | |
---|
| 1354 | ! tau diagnostics only with sunlight |
---|
| 1355 | |
---|
| 1356 | boxtau(j,ibox) = tau(j,ibox) |
---|
| 1357 | |
---|
| 1358 | !convert optical thickness to albedo |
---|
[1344] | 1359 | albedocld(j,ibox) |
---|
[524] | 1360 | & =real(invtau(min(nint(100.*tau(j,ibox)),45000))) |
---|
[1344] | 1361 | |
---|
[524] | 1362 | !contribute to averaging |
---|
[1344] | 1363 | meanalbedocld(j) = meanalbedocld(j) |
---|
[524] | 1364 | & +albedocld(j,ibox)*boxarea |
---|
| 1365 | |
---|
| 1366 | endif |
---|
| 1367 | |
---|
| 1368 | endif |
---|
| 1369 | |
---|
| 1370 | cIM |
---|
| 1371 | c CALL CPU_time(t2) |
---|
| 1372 | c print*,'IF box_cloudy t2 - t1',t2 - t1 |
---|
| 1373 | |
---|
| 1374 | c CALL CPU_time(t1) |
---|
| 1375 | cIM BEG |
---|
| 1376 | cIM !convert ptop to millibars |
---|
| 1377 | ptop(j,ibox)=ptop(j,ibox) / 100. |
---|
| 1378 | |
---|
| 1379 | cIM !save for output cloud top pressure and optical thickness |
---|
| 1380 | boxptop(j,ibox) = ptop(j,ibox) |
---|
| 1381 | cIM END |
---|
| 1382 | |
---|
| 1383 | cIM BEG |
---|
| 1384 | !reset itau(j), ipres(j) |
---|
| 1385 | itau(j) = 0 |
---|
| 1386 | ipres(j) = 0 |
---|
| 1387 | |
---|
| 1388 | if (tau(j,ibox) .lt. isccp_taumin) then |
---|
| 1389 | itau(j)=1 |
---|
| 1390 | else if (tau(j,ibox) .ge. isccp_taumin |
---|
| 1391 | & |
---|
| 1392 | & .and. tau(j,ibox) .lt. 1.3) then |
---|
| 1393 | itau(j)=2 |
---|
| 1394 | else if (tau(j,ibox) .ge. 1.3 |
---|
| 1395 | & .and. tau(j,ibox) .lt. 3.6) then |
---|
| 1396 | itau(j)=3 |
---|
| 1397 | else if (tau(j,ibox) .ge. 3.6 |
---|
| 1398 | & .and. tau(j,ibox) .lt. 9.4) then |
---|
| 1399 | itau(j)=4 |
---|
| 1400 | else if (tau(j,ibox) .ge. 9.4 |
---|
| 1401 | & .and. tau(j,ibox) .lt. 23.) then |
---|
| 1402 | itau(j)=5 |
---|
| 1403 | else if (tau(j,ibox) .ge. 23. |
---|
| 1404 | & .and. tau(j,ibox) .lt. 60.) then |
---|
| 1405 | itau(j)=6 |
---|
| 1406 | else if (tau(j,ibox) .ge. 60.) then |
---|
| 1407 | itau(j)=7 |
---|
| 1408 | end if |
---|
| 1409 | |
---|
| 1410 | !determine cloud top pressure category |
---|
| 1411 | if ( ptop(j,ibox) .gt. 0. |
---|
| 1412 | & .and.ptop(j,ibox) .lt. 180.) then |
---|
| 1413 | ipres(j)=1 |
---|
| 1414 | else if(ptop(j,ibox) .ge. 180. |
---|
| 1415 | & .and.ptop(j,ibox) .lt. 310.) then |
---|
| 1416 | ipres(j)=2 |
---|
| 1417 | else if(ptop(j,ibox) .ge. 310. |
---|
| 1418 | & .and.ptop(j,ibox) .lt. 440.) then |
---|
| 1419 | ipres(j)=3 |
---|
| 1420 | else if(ptop(j,ibox) .ge. 440. |
---|
| 1421 | & .and.ptop(j,ibox) .lt. 560.) then |
---|
| 1422 | ipres(j)=4 |
---|
| 1423 | else if(ptop(j,ibox) .ge. 560. |
---|
| 1424 | & .and.ptop(j,ibox) .lt. 680.) then |
---|
| 1425 | ipres(j)=5 |
---|
| 1426 | else if(ptop(j,ibox) .ge. 680. |
---|
| 1427 | & .and.ptop(j,ibox) .lt. 800.) then |
---|
| 1428 | ipres(j)=6 |
---|
| 1429 | else if(ptop(j,ibox) .ge. 800.) then |
---|
| 1430 | ipres(j)=7 |
---|
| 1431 | end if |
---|
| 1432 | cIM END |
---|
| 1433 | |
---|
| 1434 | if (sunlit(j).eq.1 .or. top_height .eq. 3) then |
---|
| 1435 | |
---|
| 1436 | cIM !convert ptop to millibars |
---|
| 1437 | cIM ptop(j,ibox)=ptop(j,ibox) / 100. |
---|
| 1438 | |
---|
| 1439 | cIM !save for output cloud top pressure and optical thickness |
---|
| 1440 | cIM boxptop(j,ibox) = ptop(j,ibox) |
---|
| 1441 | |
---|
| 1442 | if (box_cloudy(j,ibox)) then |
---|
[1344] | 1443 | |
---|
[524] | 1444 | meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea |
---|
| 1445 | |
---|
| 1446 | cIM !reset itau(j), ipres(j) |
---|
| 1447 | cIM itau(j) = 0 |
---|
| 1448 | cIM ipres(j) = 0 |
---|
| 1449 | |
---|
| 1450 | c if (tau(j,ibox) .lt. isccp_taumin) then |
---|
| 1451 | c itau(j)=1 |
---|
| 1452 | c else if (tau(j,ibox) .ge. isccp_taumin |
---|
| 1453 | c & |
---|
| 1454 | c & .and. tau(j,ibox) .lt. 1.3) then |
---|
| 1455 | c itau(j)=2 |
---|
| 1456 | c else if (tau(j,ibox) .ge. 1.3 |
---|
| 1457 | c & .and. tau(j,ibox) .lt. 3.6) then |
---|
| 1458 | c itau(j)=3 |
---|
| 1459 | c else if (tau(j,ibox) .ge. 3.6 |
---|
| 1460 | c & .and. tau(j,ibox) .lt. 9.4) then |
---|
| 1461 | c itau(j)=4 |
---|
| 1462 | c else if (tau(j,ibox) .ge. 9.4 |
---|
| 1463 | c & .and. tau(j,ibox) .lt. 23.) then |
---|
| 1464 | c itau(j)=5 |
---|
| 1465 | c else if (tau(j,ibox) .ge. 23. |
---|
| 1466 | c & .and. tau(j,ibox) .lt. 60.) then |
---|
| 1467 | c itau(j)=6 |
---|
| 1468 | c else if (tau(j,ibox) .ge. 60.) then |
---|
| 1469 | c itau(j)=7 |
---|
| 1470 | c end if |
---|
| 1471 | |
---|
| 1472 | c !determine cloud top pressure category |
---|
| 1473 | c if ( ptop(j,ibox) .gt. 0. |
---|
| 1474 | c & .and.ptop(j,ibox) .lt. 180.) then |
---|
| 1475 | c ipres(j)=1 |
---|
| 1476 | c else if(ptop(j,ibox) .ge. 180. |
---|
| 1477 | c & .and.ptop(j,ibox) .lt. 310.) then |
---|
| 1478 | c ipres(j)=2 |
---|
| 1479 | c else if(ptop(j,ibox) .ge. 310. |
---|
| 1480 | c & .and.ptop(j,ibox) .lt. 440.) then |
---|
| 1481 | c ipres(j)=3 |
---|
| 1482 | c else if(ptop(j,ibox) .ge. 440. |
---|
| 1483 | c & .and.ptop(j,ibox) .lt. 560.) then |
---|
| 1484 | c ipres(j)=4 |
---|
| 1485 | c else if(ptop(j,ibox) .ge. 560. |
---|
| 1486 | c & .and.ptop(j,ibox) .lt. 680.) then |
---|
| 1487 | c ipres(j)=5 |
---|
| 1488 | c else if(ptop(j,ibox) .ge. 680. |
---|
| 1489 | c & .and.ptop(j,ibox) .lt. 800.) then |
---|
| 1490 | c ipres(j)=6 |
---|
| 1491 | c else if(ptop(j,ibox) .ge. 800.) then |
---|
| 1492 | c ipres(j)=7 |
---|
| 1493 | c end if |
---|
| 1494 | |
---|
| 1495 | !update frequencies |
---|
| 1496 | if(ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
| 1497 | fq_isccp(j,itau(j),ipres(j))= |
---|
| 1498 | & fq_isccp(j,itau(j),ipres(j))+ boxarea |
---|
| 1499 | end if |
---|
| 1500 | |
---|
| 1501 | cIM calcul stats regime dynamique BEG |
---|
| 1502 | ! iw(j) = int((w(j)-wmin)/pas_w) +1 |
---|
| 1503 | ! pctj(itau(j),ipres(j),iw(j))=.FALSE. |
---|
| 1504 | ! !update frequencies W500 |
---|
| 1505 | ! if (pct_ocean(j)) then |
---|
| 1506 | ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
| 1507 | ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then |
---|
| 1508 | c print*,' ISCCP iw=',iw(j),j |
---|
| 1509 | ! fq_dynreg(itau(j),ipres(j),iw(j))= |
---|
| 1510 | ! & fq_dynreg(itau(j),ipres(j),iw(j))+ |
---|
| 1511 | ! & boxarea |
---|
| 1512 | c & fq_isccp(j,itau(j),ipres(j)) |
---|
| 1513 | ! pctj(itau(j),ipres(j),iw(j))=.TRUE. |
---|
| 1514 | c nfq_dynreg(itau(j),ipres(j),iw(j))= |
---|
| 1515 | c & nfq_dynreg(itau(j),ipres(j),iw(j))+1. |
---|
| 1516 | ! end if |
---|
| 1517 | ! end if |
---|
| 1518 | ! end if |
---|
| 1519 | cIM calcul stats regime dynamique END |
---|
| 1520 | end if !IM boxcloudy |
---|
| 1521 | |
---|
| 1522 | end if !IM sunlit |
---|
| 1523 | |
---|
| 1524 | cIM |
---|
| 1525 | c CALL CPU_time(t2) |
---|
| 1526 | c print*,'IF sunlit boxcloudy t2 - t1',t2 - t1 |
---|
| 1527 | cIM |
---|
| 1528 | enddo !IM ibox/j |
---|
| 1529 | |
---|
| 1530 | |
---|
| 1531 | cIM ajout stats s/ W500 BEG |
---|
| 1532 | cIM ajout stats s/ W500 END |
---|
| 1533 | |
---|
| 1534 | c if(j.EQ.6722) then |
---|
| 1535 | c print*,' ISCCP',w(j),iw(j),ipres(j),itau(j) |
---|
| 1536 | c endif |
---|
| 1537 | |
---|
| 1538 | ! if (pct_ocean(j)) then |
---|
| 1539 | ! if (ipres(j) .gt. 0.and.itau(j) .gt. 0) then |
---|
| 1540 | ! if (iw(j) .gt. int(wmin).and.iw(j) .le. iwmx) then |
---|
| 1541 | ! if(pctj(itau(j),ipres(j),iw(j))) THEN |
---|
| 1542 | ! nfq_dynreg(itau(j),ipres(j),iw(j))= |
---|
| 1543 | ! & nfq_dynreg(itau(j),ipres(j),iw(j))+1. |
---|
| 1544 | c if(itau(j).EQ.4.AND.ipres(j).EQ.2.AND. |
---|
| 1545 | c & iw(j).EQ.10) then |
---|
| 1546 | c PRINT*,' isccp AVANT', |
---|
| 1547 | c & nfq_dynreg(itau(j),ipres(j),iw(j)), |
---|
| 1548 | c & fq_dynreg(itau(j),ipres(j),iw(j)) |
---|
| 1549 | c endif |
---|
| 1550 | ! endif |
---|
| 1551 | ! endif |
---|
| 1552 | ! endif |
---|
| 1553 | ! endif |
---|
| 1554 | 39 continue !IM j/ibox |
---|
| 1555 | |
---|
| 1556 | !compute mean cloud properties |
---|
| 1557 | do j=1,npoints |
---|
[1344] | 1558 | if (totalcldarea(j) .gt. 0.) then |
---|
| 1559 | meanptop(j) = meanptop(j) / totalcldarea(j) |
---|
| 1560 | if (sunlit(j).eq.1) then |
---|
| 1561 | meanalbedocld(j)=meanalbedocld(j) / totalcldarea(j) |
---|
| 1562 | meantaucld(j)=tautab(min(255,max(1,nint(meanalbedocld(j))))) |
---|
| 1563 | end if |
---|
| 1564 | end if |
---|
[524] | 1565 | enddo ! j |
---|
| 1566 | ! |
---|
| 1567 | cIM ajout stats s/ W500 BEG |
---|
| 1568 | ! do nw = 1, iwmx |
---|
| 1569 | ! do l = 1, 7 |
---|
| 1570 | ! do k = 1, 7 |
---|
| 1571 | ! if (nfq_dynreg(k,l,nw).GT.0.) then |
---|
| 1572 | ! fq_dynreg(k,l,nw) = fq_dynreg(k,l,nw)/nfq_dynreg(k,l,nw) |
---|
| 1573 | c if(k.EQ.4.AND.l.EQ.2.AND.nw.EQ.10) then |
---|
| 1574 | c print*,' isccp APRES',nfq_dynreg(k,l,nw), |
---|
| 1575 | c & fq_dynreg(k,l,nw) |
---|
| 1576 | c endif |
---|
| 1577 | ! else |
---|
| 1578 | ! if(fq_dynreg(k,l,nw).NE.0.) then |
---|
| 1579 | ! print*,'nfq_dynreg = 0 tau,pc,nw',k,l,nw,fq_dynreg(k,l,nw) |
---|
| 1580 | ! endif |
---|
| 1581 | c fq_dynreg(k,l,nw) = -1.E+20 |
---|
| 1582 | c nfq_dynreg(k,l,nw) = 1.E+20 |
---|
| 1583 | ! end if |
---|
| 1584 | ! enddo !k |
---|
| 1585 | ! enddo !l |
---|
| 1586 | ! enddo !nw |
---|
| 1587 | cIM ajout stats s/ W500 END |
---|
| 1588 | ! ---------------------------------------------------! |
---|
| 1589 | |
---|
| 1590 | ! ---------------------------------------------------! |
---|
| 1591 | ! OPTIONAL PRINTOUT OF DATA TO CHECK PROGRAM |
---|
| 1592 | ! |
---|
| 1593 | !cIM |
---|
| 1594 | c if (debugcol.ne.0) then |
---|
| 1595 | ! |
---|
| 1596 | c do j=1,npoints,debugcol |
---|
| 1597 | |
---|
| 1598 | c !produce character output |
---|
| 1599 | c do ilev=1,nlev |
---|
| 1600 | c do ibox=1,ncol |
---|
| 1601 | c acc(ilev,ibox)=0 |
---|
| 1602 | c enddo |
---|
| 1603 | c enddo |
---|
| 1604 | |
---|
| 1605 | c do ilev=1,nlev |
---|
| 1606 | c do ibox=1,ncol |
---|
| 1607 | c acc(ilev,ibox)=frac_out(j,ibox,ilev)*2 |
---|
| 1608 | c if (levmatch(j,ibox) .eq. ilev) |
---|
| 1609 | c & acc(ilev,ibox)=acc(ilev,ibox)+1 |
---|
| 1610 | c enddo |
---|
| 1611 | c enddo |
---|
| 1612 | |
---|
| 1613 | !print test |
---|
| 1614 | |
---|
| 1615 | c write(ftn09,11) j |
---|
| 1616 | c11 format('ftn09.',i4.4) |
---|
| 1617 | c open(9, FILE=ftn09, FORM='FORMATTED') |
---|
| 1618 | |
---|
| 1619 | c write(9,'(a1)') ' ' |
---|
| 1620 | c write(9,'(10i5)') |
---|
| 1621 | c & (ilev,ilev=5,nlev,5) |
---|
| 1622 | c write(9,'(a1)') ' ' |
---|
| 1623 | |
---|
| 1624 | c do ibox=1,ncol |
---|
| 1625 | c write(9,'(40(a1),1x,40(a1))') |
---|
| 1626 | c & (cchar_realtops(acc(ilev,ibox)+1),ilev=1,nlev) |
---|
| 1627 | c & ,(cchar(acc(ilev,ibox)+1),ilev=1,nlev) |
---|
| 1628 | c end do |
---|
| 1629 | c close(9) |
---|
| 1630 | c |
---|
| 1631 | cIM |
---|
| 1632 | c if (ncolprint.ne.0) then |
---|
| 1633 | c write(6,'(a1)') ' ' |
---|
| 1634 | c write(6,'(a2,1X,5(a7,1X),a50)') |
---|
| 1635 | c & 'ilev', |
---|
| 1636 | c & 'pfull','at', |
---|
| 1637 | c & 'cc*100','dem_s','dtau_s', |
---|
| 1638 | c & 'cchar' |
---|
| 1639 | |
---|
| 1640 | ! do 4012 ilev=1,nlev |
---|
| 1641 | ! write(6,'(60i2)') (box(i,ilev),i=1,ncolprint) |
---|
| 1642 | ! write(6,'(i2,1X,5(f7.2,1X),50(a1))') |
---|
| 1643 | ! & ilev, |
---|
| 1644 | ! & pfull(j,ilev)/100.,at(j,ilev), |
---|
| 1645 | ! & cc(j,ilev)*100.0,dem_s(j,ilev),dtau_s(j,ilev) |
---|
| 1646 | ! & ,(cchar(acc(ilev,ibox)+1),ibox=1,ncolprint) |
---|
| 1647 | !4012 continue |
---|
| 1648 | c write (6,'(a)') 'skt(j):' |
---|
| 1649 | c write (6,'(8f7.2)') skt(j) |
---|
| 1650 | |
---|
| 1651 | c write (6,'(8I7)') (ibox,ibox=1,ncolprint) |
---|
[1344] | 1652 | |
---|
[524] | 1653 | c write (6,'(a)') 'tau:' |
---|
| 1654 | c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint) |
---|
| 1655 | |
---|
| 1656 | c write (6,'(a)') 'tb:' |
---|
| 1657 | c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) |
---|
| 1658 | |
---|
| 1659 | c write (6,'(a)') 'ptop:' |
---|
| 1660 | c write (6,'(8f7.2)') (ptop(j,ibox),ibox=1,ncolprint) |
---|
| 1661 | c endif |
---|
| 1662 | |
---|
| 1663 | c enddo |
---|
| 1664 | |
---|
| 1665 | c end if |
---|
| 1666 | |
---|
| 1667 | return |
---|
| 1668 | end |
---|