Index: /LMDZ4/trunk/libf/phylmd/isccp_cloud_types.F
===================================================================
--- /LMDZ4/trunk/libf/phylmd/isccp_cloud_types.F	(revision 1343)
+++ /LMDZ4/trunk/libf/phylmd/isccp_cloud_types.F	(revision 1344)
@@ -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)
