Index: LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/isccp_cloud_types.F =================================================================== --- LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/isccp_cloud_types.F (revision 1342) +++ LMDZ4/branches/LMDZ4V5.0-dev/libf/phylmd/isccp_cloud_types.F (revision 1343) @@ -1,4 +1,4 @@ ! -! $Header$ +! $Id $ ! SUBROUTINE ISCCP_CLOUD_TYPES( @@ -31,5 +31,5 @@ & boxptop &) - + ! Copyright Steve Klein and Mark Webb 2002 - all rights reserved. @@ -67,10 +67,10 @@ c ! to a different value for each model c ! gridbox it is called on, as it is -c ! possible that the choice of the samec -c ! seed value every time may introduce some -c ! statistical bias in the results, particularly -c ! for low values of NCOL. +c ! possible that the choice of the samec +c ! seed value every time may introduce some +c ! statistical bias in the results, particularly +c ! for low values of NCOL. c - REAL pfull(npoints,nlev) ! pressure of full model levels (Pascals) + REAL pfull(npoints,nlev) ! pressure of full model levels (Pascals) c ! pfull(npoints,1) is top level of model c ! pfull(npoints,nlev) is bottom level of model @@ -92,5 +92,5 @@ REAL dtau_s(npoints,nlev) ! mean 0.67 micron optical depth of stratiform -c ! clouds in each model level +c ! clouds in each model level c ! NOTE: this the cloud optical depth of only the c ! cloudy part of the grid box, it is not weighted @@ -99,11 +99,11 @@ REAL dtau_c(npoints,nlev) ! mean 0.67 micron optical depth of convective -c ! clouds in each +c ! clouds in each c ! model level. Same note applies as in dtau_s. INTEGER overlap ! overlap type - + ! 1=max - + ! 2=rand ! 3=max/rand @@ -111,15 +111,15 @@ INTEGER top_height ! 1 = adjust top height using both a computed c ! infrared brightness temperature and the visible -c ! optical depth to adjust cloud top pressure. Note -c ! that this calculation is most appropriate to compare -c ! to ISCCP data during sunlit hours. +c ! optical depth to adjust cloud top pressure. Note +c ! that this calculation is most appropriate to compare +c ! to ISCCP data during sunlit hours. c ! 2 = do not adjust top height, that is cloud top c ! pressure is the actual cloud top pressure c ! in the model -c ! 3 = adjust top height using only the computed -c ! infrared brightness temperature. Note that this -c ! calculation is most appropriate to compare to ISCCP -c ! IR only algortihm (i.e. you can compare to nighttime -c ! ISCCP data with this option) +c ! 3 = adjust top height using only the computed +c ! infrared brightness temperature. Note that this +c ! calculation is most appropriate to compare to ISCCP +c ! IR only algortihm (i.e. you can compare to nighttime +c ! ISCCP data with this option) REAL tautab(0:255) ! ISCCP table for converting count value to @@ -135,8 +135,8 @@ REAL at(npoints,nlev) ! temperature in each model level (K) REAL dem_s(npoints,nlev) ! 10.5 micron longwave emissivity of stratiform -c ! clouds in each +c ! clouds in each c ! model level. Same note applies as in dtau_s. REAL dem_c(npoints,nlev) ! 10.5 micron longwave emissivity of convective -c ! clouds in each +c ! clouds in each c ! model level. Same note applies as in dtau_s. cIM reg.dyn BEG @@ -161,13 +161,13 @@ REAL totalcldarea(npoints) ! the fraction of model grid box columns c ! with cloud somewhere in them. This should -c ! equal the sum over all entries of fq_isccp - - +c ! equal the sum over all entries of fq_isccp + + c ! The following three means are averages over the cloudy areas only. If no -c ! clouds are in grid box all three quantities should equal zero. - +c ! clouds are in grid box all three quantities should equal zero. + REAL meanptop(npoints) ! mean cloud top pressure (mb) - linear averaging c ! in cloud top pressure. - + REAL meantaucld(npoints) ! mean optical thickness c ! linear averaging in albedo performed. @@ -176,6 +176,6 @@ REAL boxptop(npoints,ncol) ! cloud top pressure (mb) in each column - - + + ! ! ------ @@ -184,7 +184,7 @@ REAL frac_out(npoints,ncol,nlev) ! boxes gridbox divided up into -c ! Equivalent of BOX in original version, but -c ! indexed by column then row, rather than -c ! by row then column +c ! Equivalent of BOX in original version, but +c ! indexed by column then row, rather than +c ! by row then column REAL tca(npoints,0:nlev) ! total cloud cover in each model level (fraction) @@ -205,8 +205,8 @@ REAL dem(npoints,ncol),bb(npoints) ! working variables for 10.5 micron longwave -c ! emissivity in part of -c ! gridbox under consideration - - REAL ran(npoints) ! vector of random numbers +c ! emissivity in part of +c ! gridbox under consideration + + REAL ran(npoints) ! vector of random numbers REAL ptrop(npoints) REAL attrop(npoints) @@ -248,7 +248,7 @@ real boxarea integer debug ! set to non-zero value to print out inputs -c ! with step debug +c ! with step debug integer debugcol ! set to non-zero value to print out column -c ! decomposition with step debugcol +c ! decomposition with step debugcol integer index1(npoints),num1,jj @@ -293,5 +293,5 @@ c write(6,'(a10)') 'invtau(1:100)=' c write(6,'(8i10)') (invtau(i),i=1,100) -c do j=1,npoints,debug +c do j=1,npoints,debug c write(6,'(a10)') 'j=' c write(6,'(8I10)') j @@ -322,5 +322,5 @@ c write(6,'(a10)') 'dem_c=' c write(6,'(8f10.2)') (dem_c(j,i),i=1,nlev) -c enddo +c enddo c endif @@ -361,5 +361,5 @@ cIM c if (ncolprint.ne.0) then -c do j=1,npoints,1000 +c do j=1,npoints,1000 c write(6,'(a10)') 'j=' c write(6,'(8I10)') j @@ -375,5 +375,5 @@ c write (6,'(8f5.2)') c & ((cca(j,ilev),ibox=1,ncolprint),ilev=1,nlev) -c enddo +c enddo c endif @@ -415,15 +415,15 @@ c do j=1,npoints c if (cc(j,ilev) .lt. 0. .or. cc(j,ilev) .gt. 1.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) ' error = cloud fraction less than zero' -c write(6,*) ' or ' +c write(6,*) ' or ' c write(6,*) ' error = cloud fraction greater than 1' -c write(6,*) 'value at point ',j,' is ',cc(j,ilev) -c write(6,*) 'level ',ilev +c write(6,*) 'value at point ',j,' is ',cc(j,ilev) +c write(6,*) 'level ',ilev c STOP c enddo @@ -431,17 +431,17 @@ c do j=1,npoints c if (conv(j,ilev) .lt. 0. .or. conv(j,ilev) .gt. 1.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) c & ' error = convective cloud fraction less than zero' -c write(6,*) ' or ' +c write(6,*) ' or ' c write(6,*) c & ' error = convective cloud fraction greater than 1' -c write(6,*) 'value at point ',j,' is ',conv(j,ilev) -c write(6,*) 'level ',ilev +c write(6,*) 'value at point ',j,' is ',conv(j,ilev) +c write(6,*) 'level ',ilev c STOP c enddo @@ -450,14 +450,14 @@ c do j=1,npoints c if (dtau_s(j,ilev) .lt. 0.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) c & ' error = stratiform cloud opt. depth less than zero' -c write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) -c write(6,*) 'level ',ilev +c write(6,*) 'value at point ',j,' is ',dtau_s(j,ilev) +c write(6,*) 'level ',ilev c STOP c enddo @@ -465,14 +465,14 @@ c do j=1,npoints c if (dtau_c(j,ilev) .lt. 0.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) c & ' error = convective cloud opt. depth less than zero' -c write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) -c write(6,*) 'level ',ilev +c write(6,*) 'value at point ',j,' is ',dtau_c(j,ilev) +c write(6,*) 'level ',ilev c STOP c enddo @@ -481,17 +481,17 @@ c do j=1,npoints c if (dem_s(j,ilev) .lt. 0. .or. dem_s(j,ilev) .gt. 1.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) c & ' error = stratiform cloud emissivity less than zero' -c write(6,*)'or' +c write(6,*)'or' c write(6,*) c & ' error = stratiform cloud emissivity greater than 1' -c write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) -c write(6,*) 'level ',ilev +c write(6,*) 'value at point ',j,' is ',dem_s(j,ilev) +c write(6,*) 'level ',ilev c STOP c enddo @@ -500,13 +500,13 @@ c do j=1,npoints c if (dem_c(j,ilev) .lt. 0. .or. dem_c(j,ilev) .gt. 1.) then -c num1=num1+1 -c index1(num1)=j +c num1=num1+1 +c index1(num1)=j c end if c enddo c do jj=1,num1 -c j=index1(jj) +c j=index1(jj) c write(6,*) c & ' error = convective cloud emissivity less than zero' -c write(6,*)'or' +c write(6,*)'or' c write(6,*) c & ' error = convective cloud emissivity greater than 1' @@ -520,5 +520,5 @@ do ibox=1,ncol do j=1,npoints - boxpos(j,ibox)=(ibox-.5)/ncol + boxpos(j,ibox)=(ibox-.5)/ncol enddo enddo @@ -533,5 +533,5 @@ do ibox=1,ncol do j=1,npoints - frac_out(j,ibox,ilev)=0.0 + frac_out(j,ibox,ilev)=0.0 enddo enddo @@ -572,8 +572,8 @@ IF (ilev.eq.1) then - ! If max overlap - IF (overlap.eq.1) then - ! select pixels spread evenly - ! across the gridbox + ! If max overlap + IF (overlap.eq.1) then + ! select pixels spread evenly + ! across the gridbox DO ibox=1,ncol do j=1,npoints @@ -581,10 +581,10 @@ enddo enddo - ELSE + ELSE DO ibox=1,ncol call ran0_vec(npoints,seed,ran) - ! select random pixels from the non-convective - ! part the gridbox ( some will be converted into - ! convective pixels below ) + ! select random pixels from the non-convective + ! part the gridbox ( some will be converted into + ! convective pixels below ) do j=1,npoints threshold(j,ibox)= @@ -659,5 +659,5 @@ ! Reset threshold call ran0_vec(npoints,seed,ran) - + do j=1,npoints threshold(j,ibox)= @@ -698,6 +698,6 @@ ENDDO -! Code to partition boxes into startiform and convective parts -! goes here +! Code to partition boxes into startiform and convective parts +! goes here DO ibox=1,ncol @@ -744,5 +744,5 @@ c write (6,'(8f5.2)') c & ((frac_out(j,ibox,ilev2),ibox=1,ncolprint),ilev2=1,nlev) -c enddo +c enddo c endif @@ -762,8 +762,8 @@ do j=1,npoints tau(j,ibox)=0. - albedocld(j,ibox)=0. - boxtau(j,ibox)=0. - boxptop(j,ibox)=0. - box_cloudy(j,ibox)=.false. + albedocld(j,ibox)=0. + boxtau(j,ibox)=0. + boxptop(j,ibox)=0. + box_cloudy(j,ibox)=.false. enddo 15 continue @@ -865,5 +865,5 @@ c & tauwv(j),dem_wv(j,ilev) c enddo -c endif +c endif 125 continue @@ -879,9 +879,9 @@ ! Black body emission at temperature of the layer - bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) - !bb(j)= 5.67e-8*at(j,ilev)**4 - - ! increase TOA flux by flux emitted from layer - ! times total transmittance in layers above + bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) + !bb(j)= 5.67e-8*at(j,ilev)**4 + + ! increase TOA flux by flux emitted from layer + ! times total transmittance in layers above fluxtop_clrsky(j) = fluxtop_clrsky(j) @@ -889,5 +889,5 @@ ! update trans_layers_above with transmissivity - ! from this layer for next time around loop + ! from this layer for next time around loop trans_layers_above_clrsky(j)= @@ -935,5 +935,5 @@ c & trans_layers_above_clrsky(j) c enddo -c endif +c endif @@ -971,6 +971,6 @@ ! Black body emission at temperature of the layer - bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) - !bb(j)= 5.67e-8*at(j,ilev)**4 + bb(j)=1 / ( exp(1307.27/at(j,ilev)) - 1. ) + !bb(j)= 5.67e-8*at(j,ilev)**4 enddo @@ -978,5 +978,5 @@ do j=1,npoints - ! emissivity for point in this layer + ! emissivity for point in this layer cIM REAL if (frac_out(j,ibox,ilev).eq.1) then if (frac_out(j,ibox,ilev).eq.1.0) then @@ -993,5 +993,5 @@ ! increase TOA flux by flux emitted from layer - ! times total transmittance in layers above + ! times total transmittance in layers above fluxtop(j,ibox) = fluxtop(j,ibox) @@ -1000,5 +1000,5 @@ ! update trans_layers_above with transmissivity - ! from this layer for next time around loop + ! from this layer for next time around loop trans_layers_above(j,ibox)= @@ -1029,5 +1029,5 @@ c write (6,'(8f7.2)') c & (trans_layers_above(j,ibox),ibox=1,ncolprint) -c enddo +c enddo c endif @@ -1071,5 +1071,5 @@ c write (6,'(8f7.2)') (100.*fluxtop(j,ibox),ibox=1,ncolprint) c end do -c endif +c endif !now that you have the top of atmosphere radiance account @@ -1082,8 +1082,8 @@ !adjustment performed in this section ! - !If it turns out that the cloud brightness temperature - !is greater than 260K, then the liquid cloud conversion + !If it turns out that the cloud brightness temperature + !is greater than 260K, then the liquid cloud conversion !factor of 2.56 is used. - ! + ! !Note that this is discussed on pages 85-87 of !the ISCCP D level documentation (Rossow et al. 1996) @@ -1097,7 +1097,7 @@ transmax(j) = (fluxtop(j,ibox)-btcmin(j)) & /(fluxtop_clrsky(j)-btcmin(j)) - !note that the initial setting of tauir(j) is needed so that - !tauir(j) has a realistic value should the next if block be - !bypassed + !note that the initial setting of tauir(j) is needed so that + !tauir(j) has a realistic value should the next if block be + !bypassed tauir(j) = tau(j,ibox) * rec2p13 taumin(j) = -1. * log(max(min(transmax(j),0.9999999),0.001)) @@ -1110,5 +1110,5 @@ & transmax(j) .le. 0.9999999) then fluxtopinit(j) = fluxtop(j,ibox) - tauir(j) = tau(j,ibox) *rec2p13 + tauir(j) = tau(j,ibox) *rec2p13 endif enddo @@ -1126,6 +1126,6 @@ & / (log(1. + (1./fluxtop(j,ibox)))) if (tb(j,ibox) .gt. 260.) then - tauir(j) = tau(j,ibox) / 2.56 - end if + tauir(j) = tau(j,ibox) / 2.56 + end if end if end if @@ -1141,5 +1141,5 @@ if (top_height.eq.1.and.tauir(j).lt.taumin(j)) then tb(j,ibox) = attrop(j) - 5. - tau(j,ibox) = 2.13*taumin(j) + tau(j,ibox) = 2.13*taumin(j) end if else @@ -1181,18 +1181,18 @@ c write (6,'(a)') 'total_trans:' c write (6,'(8f7.2)') -c & (trans_layers_above(j,ibox),ibox=1,ncolprint) +c & (trans_layers_above(j,ibox),ibox=1,ncolprint) c write (6,'(a)') 'total_emiss:' c write (6,'(8f7.2)') -c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) +c & (1.0-trans_layers_above(j,ibox),ibox=1,ncolprint) c write (6,'(a)') 'total_trans:' c write (6,'(8f7.2)') -c & (trans_layers_above(j,ibox),ibox=1,ncolprint) +c & (trans_layers_above(j,ibox),ibox=1,ncolprint) c write (6,'(a)') 'ppout:' c write (6,'(8f7.2)') (tb(j,ibox),ibox=1,ncolprint) c enddo ! j -c endif +c endif end if @@ -1264,5 +1264,5 @@ & .and.(frac_out(j,ibox,ilev) .ne. 0.0)) then ptop(j,ibox)=pfull(j,ilev) - levmatch(j,ibox)=ilev + levmatch(j,ibox)=ilev end if end do @@ -1357,9 +1357,9 @@ !convert optical thickness to albedo - albedocld(j,ibox) + albedocld(j,ibox) & =real(invtau(min(nint(100.*tau(j,ibox)),45000))) - + !contribute to averaging - meanalbedocld(j) = meanalbedocld(j) + meanalbedocld(j) = meanalbedocld(j) & +albedocld(j,ibox)*boxarea @@ -1441,5 +1441,5 @@ if (box_cloudy(j,ibox)) then - + meanptop(j) = meanptop(j) + ptop(j,ibox)*boxarea @@ -1556,11 +1556,11 @@ !compute mean cloud properties do j=1,npoints - if (totalcldarea(j) .gt. 0.) then - meanptop(j) = meanptop(j) / totalcldarea(j) - if (sunlit(j).eq.1) then - meanalbedocld(j) = meanalbedocld(j) / totalcldarea(j) - meantaucld(j) = tautab(min(255,max(1,nint(meanalbedocld(j))))) - end if - end if + if (totalcldarea(j) .gt. 0.) then + meanptop(j) = meanptop(j) / totalcldarea(j) + if (sunlit(j).eq.1) then + meanalbedocld(j)=meanalbedocld(j) / totalcldarea(j) + meantaucld(j)=tautab(min(255,max(1,nint(meanalbedocld(j))))) + end if + end if enddo ! j ! @@ -1650,5 +1650,5 @@ c write (6,'(8I7)') (ibox,ibox=1,ncolprint) - + c write (6,'(a)') 'tau:' c write (6,'(8f7.2)') (tau(j,ibox),ibox=1,ncolprint)