Changeset 1343 for LMDZ4/branches/LMDZ4V5.0-dev/libf
- Timestamp:
- Apr 9, 2010, 9:54:52 AM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/isccp_cloud_types.F
r1146 r1343 1 1 ! 2 ! $ Header$2 ! $Id $ 3 3 ! 4 4 SUBROUTINE ISCCP_CLOUD_TYPES( … … 31 31 & boxptop 32 32 &) 33 33 34 34 35 35 ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved. … … 67 67 c ! to a different value for each model 68 68 c ! gridbox it is called on, as it is 69 c 70 c 71 c 72 c 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. 73 73 c 74 REAL pfull(npoints,nlev) 74 REAL pfull(npoints,nlev) ! pressure of full model levels (Pascals) 75 75 c ! pfull(npoints,1) is top level of model 76 76 c ! pfull(npoints,nlev) is bottom level of model … … 92 92 93 93 REAL dtau_s(npoints,nlev) ! mean 0.67 micron optical depth of stratiform 94 c 94 c ! clouds in each model level 95 95 c ! NOTE: this the cloud optical depth of only the 96 96 c ! cloudy part of the grid box, it is not weighted … … 99 99 100 100 REAL dtau_c(npoints,nlev) ! mean 0.67 micron optical depth of convective 101 c 101 c ! clouds in each 102 102 c ! model level. Same note applies as in dtau_s. 103 103 104 104 INTEGER overlap ! overlap type 105 105 106 106 ! 1=max 107 107 108 108 ! 2=rand 109 109 ! 3=max/rand … … 111 111 INTEGER top_height ! 1 = adjust top height using both a computed 112 112 c ! infrared brightness temperature and the visible 113 c 114 c 115 c 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. 116 116 c ! 2 = do not adjust top height, that is cloud top 117 117 c ! pressure is the actual cloud top pressure 118 118 c ! in the model 119 c 120 c 121 c 122 c 123 c 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) 124 124 125 125 REAL tautab(0:255) ! ISCCP table for converting count value to … … 135 135 REAL at(npoints,nlev) ! temperature in each model level (K) 136 136 REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform 137 c 137 c ! clouds in each 138 138 c ! model level. Same note applies as in dtau_s. 139 139 REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective 140 c 140 c ! clouds in each 141 141 c ! model level. Same note applies as in dtau_s. 142 142 cIM reg.dyn BEG … … 161 161 REAL totalcldarea(npoints) ! the fraction of model grid box columns 162 162 c ! with cloud somewhere in them. This should 163 c 164 165 163 c ! equal the sum over all entries of fq_isccp 164 165 166 166 c ! The following three means are averages over the cloudy areas only. If no 167 c ! clouds are in grid box all three quantities should equal zero. 168 167 c ! clouds are in grid box all three quantities should equal zero. 168 169 169 REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging 170 170 c ! in cloud top pressure. 171 171 172 172 REAL meantaucld(npoints) ! mean optical thickness 173 173 c ! linear averaging in albedo performed. … … 176 176 177 177 REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column 178 179 178 179 180 180 ! 181 181 ! ------ … … 184 184 185 185 REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into 186 c 187 c 188 c 186 c ! Equivalent of BOX in original version, but 187 c ! indexed by column then row, rather than 188 c ! by row then column 189 189 190 190 REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction) … … 205 205 206 206 REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave 207 c 208 c 209 210 REAL ran(npoints) 207 c ! emissivity in part of 208 c ! gridbox under consideration 209 210 REAL ran(npoints) ! vector of random numbers 211 211 REAL ptrop(npoints) 212 212 REAL attrop(npoints) … … 248 248 real boxarea 249 249 integer debug ! set to non-zero value to print out inputs 250 c 250 c ! with step debug 251 251 integer debugcol ! set to non-zero value to print out column 252 c 252 c ! decomposition with step debugcol 253 253 254 254 integer index1(npoints),num1,jj … … 293 293 c write(6,'(a10)') 'invtau(1:100)=' 294 294 c write(6,'(8i10)') (invtau(i),i=1,100) 295 c 295 c do j=1,npoints,debug 296 296 c write(6,'(a10)') 'j=' 297 297 c write(6,'(8I10)') j … … 322 322 c write(6,'(a10)') 'dem_c=' 323 323 c write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev) 324 c 324 c enddo 325 325 c endif 326 326 … … 361 361 cIM 362 362 c if (ncolprint.ne.0) then 363 c 363 c do j=1,npoints,1000 364 364 c write(6,'(a10)') 'j=' 365 365 c write(6,'(8I10)') j … … 375 375 c write (6,'(8f5.2)') 376 376 c & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev) 377 c 377 c enddo 378 378 c endif 379 379 … … 415 415 c do j=1,npoints 416 416 c if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then 417 c 418 c 417 c num1=num1+1 418 c index1(num1)=j 419 419 c end if 420 420 c enddo 421 421 c do jj=1,num1 422 c 422 c j=index1(jj) 423 423 c write(6,*) ' error = cloud fraction less than zero' 424 c 424 c write(6,*) ' or ' 425 425 c write(6,*) ' error = cloud fraction greater than 1' 426 c 427 c 426 c write(6,*) 'value at point ',j,' is ',cc(j,ilev) 427 c write(6,*) 'level ',ilev 428 428 c STOP 429 429 c enddo … … 431 431 c do j=1,npoints 432 432 c if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then 433 c 434 c 433 c num1=num1+1 434 c index1(num1)=j 435 435 c end if 436 436 c enddo 437 437 c do jj=1,num1 438 c 438 c j=index1(jj) 439 439 c write(6,*) 440 440 c & ' error = convective cloud fraction less than zero' 441 c 441 c write(6,*) ' or ' 442 442 c write(6,*) 443 443 c & ' error = convective cloud fraction greater than 1' 444 c 445 c 444 c write(6,*) 'value at point ',j,' is ',conv(j,ilev) 445 c write(6,*) 'level ',ilev 446 446 c STOP 447 447 c enddo … … 450 450 c do j=1,npoints 451 451 c if (dtau_s(j,ilev) .lt. 0.) then 452 c 453 c 452 c num1=num1+1 453 c index1(num1)=j 454 454 c end if 455 455 c enddo 456 456 c do jj=1,num1 457 c 457 c j=index1(jj) 458 458 c write(6,*) 459 459 c & ' error = stratiform cloud opt. depth less than zero' 460 c 461 c 460 c write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) 461 c write(6,*) 'level ',ilev 462 462 c STOP 463 463 c enddo … … 465 465 c do j=1,npoints 466 466 c if (dtau_c(j,ilev) .lt. 0.) then 467 c 468 c 467 c num1=num1+1 468 c index1(num1)=j 469 469 c end if 470 470 c enddo 471 471 c do jj=1,num1 472 c 472 c j=index1(jj) 473 473 c write(6,*) 474 474 c & ' error = convective cloud opt. depth less than zero' 475 c 476 c 475 c write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) 476 c write(6,*) 'level ',ilev 477 477 c STOP 478 478 c enddo … … 481 481 c do j=1,npoints 482 482 c if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then 483 c 484 c 483 c num1=num1+1 484 c index1(num1)=j 485 485 c end if 486 486 c enddo 487 487 c do jj=1,num1 488 c 488 c j=index1(jj) 489 489 c write(6,*) 490 490 c & ' error = stratiform cloud emissivity less than zero' 491 c 491 c write(6,*)'or' 492 492 c write(6,*) 493 493 c & ' error = stratiform cloud emissivity greater than 1' 494 c 495 c 494 c write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) 495 c write(6,*) 'level ',ilev 496 496 c STOP 497 497 c enddo … … 500 500 c do j=1,npoints 501 501 c if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then 502 c 503 c 502 c num1=num1+1 503 c index1(num1)=j 504 504 c end if 505 505 c enddo 506 506 c do jj=1,num1 507 c 507 c j=index1(jj) 508 508 c write(6,*) 509 509 c & ' error = convective cloud emissivity less than zero' 510 c 510 c write(6,*)'or' 511 511 c write(6,*) 512 512 c & ' error = convective cloud emissivity greater than 1' … … 520 520 do ibox=1,ncol 521 521 do j=1,npoints 522 522 boxpos(j,ibox)=(ibox-.5)/ncol 523 523 enddo 524 524 enddo … … 533 533 do ibox=1,ncol 534 534 do j=1,npoints 535 535 frac_out(j,ibox,ilev)=0.0 536 536 enddo 537 537 enddo … … 572 572 573 573 IF (ilev.eq.1) then 574 575 576 577 574 ! If max overlap 575 IF (overlap.eq.1) then 576 ! select pixels spread evenly 577 ! across the gridbox 578 578 DO ibox=1,ncol 579 579 do j=1,npoints … … 581 581 enddo 582 582 enddo 583 583 ELSE 584 584 DO ibox=1,ncol 585 585 call ran0_vec(npoints,seed,ran) 586 587 588 586 ! select random pixels from the non-convective 587 ! part the gridbox ( some will be converted into 588 ! convective pixels below ) 589 589 do j=1,npoints 590 590 threshold(j,ibox)= … … 659 659 ! Reset threshold 660 660 call ran0_vec(npoints,seed,ran) 661 661 662 662 do j=1,npoints 663 663 threshold(j,ibox)= … … 698 698 ENDDO 699 699 700 ! 701 ! 700 ! Code to partition boxes into startiform and convective parts 701 ! goes here 702 702 703 703 DO ibox=1,ncol … … 744 744 c write (6,'(8f5.2)') 745 745 c & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev) 746 c 746 c enddo 747 747 c endif 748 748 … … 762 762 do j=1,npoints 763 763 tau(j,ibox)=0. 764 765 766 767 764 albedocld(j,ibox)=0. 765 boxtau(j,ibox)=0. 766 boxptop(j,ibox)=0. 767 box_cloudy(j,ibox)=.false. 768 768 enddo 769 769 15 continue … … 865 865 c & tauwv(j),dem_wv(j,ilev) 866 866 c enddo 867 c 867 c endif 868 868 125 continue 869 869 … … 879 879 ! Black body emission at temperature of the layer 880 880 881 882 883 884 885 881 bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) 882 !bb(j)= 5.67e-8*at(j,ilev)**4 883 884 ! increase TOA flux by flux emitted from layer 885 ! times total transmittance in layers above 886 886 887 887 fluxtop_clrsky(j) = fluxtop_clrsky(j) … … 889 889 890 890 ! update trans_layers_above with transmissivity 891 891 ! from this layer for next time around loop 892 892 893 893 trans_layers_above_clrsky(j)= … … 935 935 c & trans_layers_above_clrsky(j) 936 936 c enddo 937 c 937 c endif 938 938 939 939 … … 971 971 ! Black body emission at temperature of the layer 972 972 973 974 973 bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) 974 !bb(j)= 5.67e-8*at(j,ilev)**4 975 975 enddo 976 976 … … 978 978 do j=1,npoints 979 979 980 980 ! emissivity for point in this layer 981 981 cIM REAL if (frac_out(j,ibox,ilev).eq.1) then 982 982 if (frac_out(j,ibox,ilev).eq.1.0) then … … 993 993 994 994 ! increase TOA flux by flux emitted from layer 995 995 ! times total transmittance in layers above 996 996 997 997 fluxtop(j,ibox) = fluxtop(j,ibox) … … 1000 1000 1001 1001 ! update trans_layers_above with transmissivity 1002 1002 ! from this layer for next time around loop 1003 1003 1004 1004 trans_layers_above(j,ibox)= … … 1029 1029 c write (6,'(8f7.2)') 1030 1030 c & (trans_layers_above(j,ibox),ibox=1,ncolprint) 1031 c 1031 c enddo 1032 1032 c endif 1033 1033 … … 1071 1071 c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) 1072 1072 c end do 1073 c 1073 c endif 1074 1074 1075 1075 !now that you have the top of atmosphere radiance account … … 1082 1082 !adjustment performed in this section 1083 1083 ! 1084 1085 1084 !If it turns out that the cloud brightness temperature 1085 !is greater than 260K, then the liquid cloud conversion 1086 1086 !factor of 2.56 is used. 1087 1087 ! 1088 1088 !Note that this is discussed on pages 85-87 of 1089 1089 !the ISCCP D level documentation (Rossow et al. 1996) … … 1097 1097 transmax(j) = (fluxtop(j,ibox)-btcmin(j)) 1098 1098 & /(fluxtop_clrsky(j)-btcmin(j)) 1099 1100 1101 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 1102 1102 tauir(j) = tau(j,ibox) * rec2p13 1103 1103 taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001)) … … 1110 1110 & transmax(j) .le. 0.9999999) then 1111 1111 fluxtopinit(j) = fluxtop(j,ibox) 1112 1112 tauir(j) = tau(j,ibox) *rec2p13 1113 1113 endif 1114 1114 enddo … … 1126 1126 & / (log(1. + (1./fluxtop(j,ibox)))) 1127 1127 if (tb(j,ibox) .gt. 260.) then 1128 1129 end if 1128 tauir(j) = tau(j,ibox) / 2.56 1129 end if 1130 1130 end if 1131 1131 end if … … 1141 1141 if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then 1142 1142 tb(j,ibox) = attrop(j) - 5. 1143 1143 tau(j,ibox) = 2.13*taumin(j) 1144 1144 end if 1145 1145 else … … 1181 1181 c write (6,'(a)') 'total_trans:' 1182 1182 c write (6,'(8f7.2)') 1183 c & 1183 c & (trans_layers_above(j,ibox),ibox=1,ncolprint) 1184 1184 1185 1185 c write (6,'(a)') 'total_emiss:' 1186 1186 c write (6,'(8f7.2)') 1187 c & 1187 c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) 1188 1188 1189 1189 c write (6,'(a)') 'total_trans:' 1190 1190 c write (6,'(8f7.2)') 1191 c & 1191 c & (trans_layers_above(j,ibox),ibox=1,ncolprint) 1192 1192 1193 1193 c write (6,'(a)') 'ppout:' 1194 1194 c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) 1195 1195 c enddo ! j 1196 c 1196 c endif 1197 1197 1198 1198 end if … … 1264 1264 & .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then 1265 1265 ptop(j,ibox)=pfull(j,ilev) 1266 1266 levmatch(j,ibox)=ilev 1267 1267 end if 1268 1268 end do … … 1357 1357 1358 1358 !convert optical thickness to albedo 1359 1359 albedocld(j,ibox) 1360 1360 & =real(invtau(min(nint(100.*tau(j,ibox)),45000))) 1361 1361 1362 1362 !contribute to averaging 1363 1363 meanalbedocld(j) = meanalbedocld(j) 1364 1364 & +albedocld(j,ibox)*boxarea 1365 1365 … … 1441 1441 1442 1442 if (box_cloudy(j,ibox)) then 1443 1443 1444 1444 meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea 1445 1445 … … 1556 1556 !compute mean cloud properties 1557 1557 do j=1,npoints 1558 1559 1560 1561 meanalbedocld(j) =meanalbedocld(j) / totalcldarea(j)1562 meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j))))) 1563 1564 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 1565 1565 enddo ! j 1566 1566 ! … … 1650 1650 1651 1651 c write (6,'(8I7)') (ibox,ibox=1,ncolprint) 1652 1652 1653 1653 c write (6,'(a)') 'tau:' 1654 1654 c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)
Note: See TracChangeset
for help on using the changeset viewer.
