srtm_cldprop.F90

Go to the documentation of this file.

SUBROUTINE srtm_cldprop &
 & ( klev, k_icldatm, k_inflag, k_iceflag, k_liqflag, k_nstr,&
 & p_cldfrac, p_clddat1, p_clddat2, p_clddat3, p_clddat4, p_clddatmom,&
 & p_taucldorig, p_taucloud, p_ssacloud, p_xmom &
 & )  

!  path:      $Source: /storm/rc1/cvsroot/rc/rrtm_sw/src/cldprop_sw.f,v $
!  author:    $Author: jdelamer $
!  revision:  $Revision: 2.6 $
!  created:   $Date: 2002/04/04 18:29:47 $

!  PURPOSE:  COMPUTE THE CLOUD OPTICAL DEPTH(S) FOR EACH CLOUDY
!  LAYER.  NOTE:  ONLY INFLAG = 0 AND INFLAG=2/LIQFLAG=1,
!  ICEFLAG=3
!  (HU & STAMNES, Q. FU) ARE IMPLEMENTED.

USE parkind1  ,ONLY : jpim     ,jprb
USE yomhook   ,ONLY : lhook,   dr_hook

USE parsrtm ,  ONLY : jplay, jpband, jpb1, jpb2
USE yoesrtop,  ONLY : extliq1, ssaliq1, asyliq1 &
 & , extice3, ssaice3, asyice3, fdlice3 &
 & , fdelta , extcoice, ssacoice, gice, forwice &
 & , extcoliq, ssacoliq, gliq, forwliq  

!     ------------------------------------------------------------------

IMPLICIT NONE

!-- real arguments

INTEGER(KIND=JPIM),INTENT(IN)    :: KLEV 
INTEGER(KIND=JPIM),INTENT(OUT)   :: K_ICLDATM 
INTEGER(KIND=JPIM),INTENT(IN)    :: K_INFLAG 
INTEGER(KIND=JPIM),INTENT(IN)    :: K_ICEFLAG 
INTEGER(KIND=JPIM),INTENT(IN)    :: K_LIQFLAG 
INTEGER(KIND=JPIM),INTENT(IN)    :: K_NSTR 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDFRAC(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT1(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT2(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT3(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDAT4(jplay) 
REAL(KIND=JPRB)   ,INTENT(IN)    :: P_CLDDATMOM(0:16,jplay) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: P_TAUCLDORIG(jplay,jpband) 
REAL(KIND=JPRB)   ,INTENT(INOUT) :: P_TAUCLOUD(jplay,jpband) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_SSACLOUD(jplay,jpband) 
REAL(KIND=JPRB)   ,INTENT(OUT)   :: P_XMOM(0:16,jplay,jpband) 
!-- integer arguments

!-- real locals
REAL(KIND=JPRB) :: Z_EPS
REAL(KIND=JPRB) :: Z_TAUCLDORIG_A, Z_FFP, Z_FFP1, Z_FFPSSA, Z_SSACLOUD_A, Z_TAUCLOUD_A
REAL(KIND=JPRB) :: Z_CWP, Z_FICE, Z_RADICE, Z_FACTOR, Z_FINT, Z_FLIQ, Z_RADLIQ
REAL(KIND=JPRB) :: Z_TAUICEORIG, Z_SCATICE, Z_SSAICE, Z_TAUICE &
 & , Z_TAULIQORIG, Z_SCATLIQ, Z_SSALIQ, Z_TAULIQ  

!-- integer locals
INTEGER(KIND=JPIM) :: I_NCBANDS, I_NLAYERS
INTEGER(KIND=JPIM) :: IB, IB1, IB2, I_LAY, ISTR , INDEX
INTEGER(KIND=JPIM) :: I_NDBUG
REAL(KIND=JPRB) :: ZHOOK_HANDLE
 
!     INCLUDE 'param.f'

!      COMMON /CONTROL/   IAER, NSTR, IOUT, ISTART, IEND, ICLD 

!      COMMON /PROFILE/   NLAYERS,PAVEL(MXLAY),TAVEL(MXLAY),
!     &                   PZ(0:MXLAY),TZ(0:MXLAY)

!      COMMON /CLOUDIN/   INFLAG,CLDDAT1(MXLAY),CLDDAT2(MXLAY),
!     &                   ICEFLAG,LIQFLAG,CLDDAT3(MXLAY),CLDDAT4(MXLAY),
!     &                   CLDDATMOM(0:16,MXLAY)

!      COMMON /CLOUDDAT/  NCBANDS,CLDFRAC(MXLAY),
!     &     TAUCLOUD(MXLAY,NBANDS),SSACLOUD(MXLAY,NBANDS),
!     &     XMOM(0:16,MXLAY,NBANDS),TAUCLDORIG(MXLAY,NBANDS)

!      COMMON /HVERSN/    HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
!     *                   HVDUM1(4),HVRUTL,HVREXT,
!     *                   HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
!     *                   HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR

!      CHARACTER*15 HVRRTM,HVRRTR,HVRATM,HVRSET,HVRTAU,
!     *            HVDUM1,HVRUTL,HVREXT,
!     *            HVRD1M,HVRR1M,HVREPK,HVRLPK,HVRAER,HVRBKA,
!     *            HVRBKB,HVRCLD,HVRDIS,HVRLAM,HVRPAR

!      DIMENSION ABSICE0(2), ABSICE1(2,5), ABSICE2(40,16)
!      DIMENSION ABSLIQ1(58,16)

!      DIMENSION ABSCOICE(NBANDS), ABSCOLIQ(NBANDS)
!      DIMENSION EXTCOLIQ(NBANDS),SSACOLIQ(NBANDS),GLIQ(NBANDS)
!      DIMENSION EXTCOICE(NBANDS),SSACOICE(NBANDS),GICE(NBANDS)
!      DIMENSION FORWLIQ(NBANDS),FORWICE(NBANDS)

! Added for SW
!      DIMENSION EXTLIQ1(58,NBANDS), SSALIQ1(58,NBANDS), 
!     &     ASYLIQ1(58,NBANDS)
!      DIMENSION EXTICE3(27,NBANDS), SSAICE3(27,NBANDS),
!     &     ASYICE3(27,NBANDS),FDLICE3(28,NBANDS)
!      DIMENSION FDELTA(NBANDS)

!DATA EPS /1.E-6/

! HVRCLD = '$Revision: 2.6 $'

IF (lhook) CALL dr_hook('SRTM_CLDPROP',0,zhook_handle)
z_eps = 1.e-06_jprb
i_ndbug = 3

k_icldatm = 0
i_ncbands = 29
i_nlayers = klev
ib1     = jpb1
ib2     = jpb2

IF (i_ndbug <= 2) THEN
  print *,'cldprop before loop K_INFLAG, K_ICEFLAG, K_LIQFLAG:',k_inflag,k_iceflag,k_liqflag,ib1,ib2
ENDIF

DO i_lay = 1, i_nlayers

  IF (p_cldfrac(i_lay) >= z_eps) THEN
    k_icldatm = 1

    IF (i_ndbug <= 2) THEN
      print 9101,i_lay,k_icldatm,p_cldfrac(i_lay),p_clddat1(i_lay),p_clddat2(i_lay),p_clddat3(i_lay)&
       & ,p_clddat4(i_lay),(p_clddatmom(istr,i_lay),istr=0,k_nstr)  
      9101  format(1x,'Cld :',2i3,f7.4,7e12.5)
    ENDIF

!   Ice clouds and water clouds combined.
    IF (k_inflag == 0) THEN
      z_taucldorig_a = p_clddat1(i_lay)
      z_ffp = p_clddatmom(k_nstr,i_lay)
      z_ffp1 = 1.0 - z_ffp
      z_ffpssa = 1.0 - z_ffp * p_clddat2(i_lay)
      z_ssacloud_a = z_ffp1*p_clddat2(i_lay)/z_ffpssa
      z_taucloud_a = z_ffpssa*z_taucldorig_a

!       DO IB = 16,NBANDS
      DO ib = ib1 , ib2
        p_taucldorig(i_lay,ib) = z_taucldorig_a
        p_ssacloud(i_lay,ib) = z_ssacloud_a
        p_taucloud(i_lay,ib) = z_taucloud_a

        DO istr = 0,k_nstr
          p_xmom(istr,i_lay,ib) = (p_clddatmom(istr,i_lay) - z_ffp)/ &
           & (z_ffp1)  
        ENDDO
      ENDDO

!   Separate treatement of ice clouds and water clouds.
    ELSEIF(k_inflag == 2) THEN       
      z_cwp = p_clddat1(i_lay)
      z_fice = p_clddat2(i_lay)
      z_radice = p_clddat3(i_lay)

      IF (i_ndbug <= 1) THEN
        print 9102,i_lay,z_cwp,z_fice,z_radice
        9102     format(1x,'A',i3,3e13.6)
      ENDIF

!   Calculation of absorption coefficients due to ice clouds.
      IF (z_fice == 0.0) THEN
!         NCBANDS = 29
!         DO IB = 16, NCBANDS
        DO ib = ib1 , ib2
          extcoice(ib) = 0.0_jprb
          ssacoice(ib) = 1.0_jprb
          gice(ib)     = 1.0_jprb
          forwice(ib)  = 0.0_jprb

          IF (i_ndbug <= 1) THEN
            print 9103,i_lay,z_fice,z_cwp,z_radice,ib,extcoice(ib),ssacoice(ib),gice(ib),forwice(ib)
            9103         format(1x,'B',i3,f6.3,2e13.6,i3,4e12.5)
          ENDIF

        ENDDO

      ELSEIF (k_iceflag == 3) THEN
        IF (z_radice < 10.0 .OR. z_radice > 140.0) stop 'ICE EFFECTIVE SIZE OUT OF BOUNDS'
!         NCBANDS = 29
        z_factor = (z_radice - 5._jprb)/5._jprb
        index = int(z_factor)
        IF (index == 27) index = 26
        z_fint = z_factor - REAL(index)

!         DO IB=16,NCBANDS
        DO ib = ib1 , ib2
          extcoice(ib) = z_fice * (extice3(index,ib) + z_fint * &
           & (extice3(index+1,ib) - extice3(index,ib)))  
          ssacoice(ib) = ssaice3(index,ib) + z_fint * &
           & (ssaice3(index+1,ib) - ssaice3(index,ib))  
          gice(ib) = asyice3(index,ib) + z_fint * &
           & (asyice3(index+1,ib) - asyice3(index,ib))  
          fdelta(ib) = fdlice3(index,ib) + z_fint * &
           & (fdlice3(index+1,ib) - fdlice3(index,ib))  
          if (fdelta(ib) < 0.0) stop 'FDELTA LESS THAN 0.0'
          if (fdelta(ib) > 1.0) stop 'FDELTA GT THAN 1.0'                     
          forwice(ib) = fdelta(ib) + 0.5 / ssacoice(ib)
! See Fu 1996 p. 2067 
          IF (forwice(ib) > gice(ib)) forwice(ib) = gice(ib)
! Check to ensure all calculated quantities are within physical limits.
          if (extcoice(ib) < 0.0_jprb) stop 'ICE EXTINCTION LESS THAN 0.0'
          if (ssacoice(ib) > 1.0_jprb) stop 'ICE SSA GRTR THAN 1.0'
          if (ssacoice(ib) < 0.0_jprb) stop 'ICE SSA LESS THAN 0.0'
          if (gice(ib) > 1.0_jprb) stop 'ICE ASYM GRTR THAN 1.0'
          if (gice(ib) < 0.0_jprb) stop 'ICE ASYM LESS THAN 0.0'

          IF (i_ndbug <= 1) THEN
            print 9104,i_lay,z_fice,z_cwp,z_radice,ib,extcoice(ib),ssacoice(ib),gice(ib),forwice(ib),fdelta(ib)
            9104         format(1x,'C',i3,f5.3,2e13.6,i3,5e12.5)
          ENDIF

        ENDDO
      ENDIF
      print *,'end of ice computations for I_LAY=',i_lay
                  
!  Calculation of absorption coefficients due to water clouds.
      z_fliq = 1. - z_fice
      IF (z_fliq == 0.0) THEN
!         NCBANDS = 29
!         DO IB = 16, NCBANDS
        DO ib = ib1 , ib2
          extcoliq(ib) = 0.0
          ssacoliq(ib) = 1.0
          gliq(ib) = 1.0
          forwliq(ib) = 0.0

          IF (i_ndbug <= 1) THEN
            print 9105,i_lay,z_fliq,z_cwp,ib,extcoliq(ib),ssacoliq(ib),gliq(ib),forwliq(ib)
            9105         format(1x,'D',i3,f5.3,1e13.6,i3,4e12.5)
          ENDIF

        ENDDO

      ELSEIF (k_liqflag == 1) THEN
        z_radliq = p_clddat4(i_lay)
        IF (z_radliq < 1.5 .OR. z_radliq > 60.) stop 'LIQUID EFFECTIVE RADIUS OUT OF BOUNDS'
        index = int(z_radliq - 1.5)
        IF (index == 0) index = 1
        IF (index == 58) index = 57
        z_fint = z_radliq - 1.5 - REAL(index)
!         NCBANDS = 29

!         DO IB = 16, NCBANDS
        DO ib = ib1 , ib2
          extcoliq(ib) = z_fliq * (extliq1(index,ib) + z_fint * &
           & (extliq1(index+1,ib) - (extliq1(index,ib))))  
          ssacoliq(ib) = ssaliq1(index,ib) + z_fint * &
           & (ssaliq1(index+1,ib) - ssaliq1(index,ib))  
          gliq(ib) = asyliq1(index,ib) + z_fint * &
           & (asyliq1(index+1,ib) - asyliq1(index,ib))  
          forwliq(ib) = gliq(ib)**k_nstr
! Check to ensure all calculated quantities are within physical limits.
          if (extcoliq(ib) < 0.0_jprb) stop 'LIQUID EXTINCTION LESS THAN 0.0'
          if (ssacoliq(ib) > 1.0_jprb) stop 'LIQUID SSA GRTR THAN 1.0'
          if (ssacoliq(ib) < 0.0_jprb) stop 'LIQUID SSA LESS THAN 0.0'
          if (gliq(ib) > 1.0_jprb) stop 'LIQUID ASYM GRTR THAN 1.0'
          if (gliq(ib) < 0.0_jprb) stop 'LIQUID ASYM LESS THAN 0.0'

          IF (i_ndbug <= 1) THEN
            print 9106,i_lay,z_fliq,z_cwp,z_radliq,ib,extcoliq(ib),ssacoliq(ib),gliq(ib),forwliq(ib)
            9106         format(1x,'E',i3,f5.3,2e13.6,i3,5e12.5)
          ENDIF

        ENDDO
      ENDIF

      IF (i_ndbug <= 1) THEN
        print *,'end of liquid water computations for I_LAY=',i_lay
      ENDIF

!       DO IB = 16, NCBANDS
      DO ib = ib1 , ib2
        z_tauliqorig = z_cwp * extcoliq(ib)
        z_tauiceorig = z_cwp * extcoice(ib)
        p_taucldorig(i_lay,ib) = z_tauliqorig + z_tauiceorig

        IF (i_ndbug <= 1) THEN
          print 9107,ib,z_tauliqorig,z_tauiceorig,p_taucldorig(i_lay,ib),z_cwp &
           & ,extcoliq(ib),extcoice(ib),ssacoliq(ib),ssacoice(ib) &
           & ,forwliq(ib),forwice(ib)  
          9107       format(1x,'F',i3,10e12.5)
        ENDIF

        z_ssaliq = ssacoliq(ib) * (1. - forwliq(ib)) / &
         & (1. - forwliq(ib) * ssacoliq(ib))  
        z_tauliq =  (1. - forwliq(ib) * ssacoliq(ib)) * &
         & z_tauliqorig  
        z_ssaice = ssacoice(ib) * (1. - forwice(ib)) / &
         & (1. - forwice(ib) * ssacoice(ib))  
        z_tauice =  (1. - forwice(ib) * ssacoice(ib)) * &
         & z_tauiceorig  
        z_scatliq = z_ssaliq * z_tauliq
        z_scatice = z_ssaice * z_tauice
        p_taucloud(i_lay,ib) = z_tauliq + z_tauice
        p_ssacloud(i_lay,ib) = (z_scatliq + z_scatice) / &
         & p_taucloud(i_lay,ib)  
        p_xmom(0,i_lay,ib) = 1.0

        IF (i_ndbug <= 1) THEN
          print 9108,ib,z_tauliqorig,z_tauiceorig,z_ssaliq,z_tauliq,z_scatliq,z_ssaice,z_tauice,z_scatice
          9108       format(1x,'G',i3,8e13.6)
        ENDIF

        DO istr = 1, k_nstr
!This commented code is the standard method for delta-m scaling. In accordance
!  with the 1996 Fu paper, equation A.3, the moments for ice were calculated
!  as in the uncommented code.
!                     XMOM(ISTR,LAY,IB) = (SCATLIQ *  &
!     &                    (GLIQ(IB)**ISTR - FORWLIQ(IB)) / &
!     &                    (1. - FORWLIQ(IB)) &
!     &                    + SCATICE * &
!     &                    (GICE(IB)**ISTR - FORWICE(IB)) /  &
!     &                    (1. - FORWICE(IB)))/(SCATLIQ + SCATICE)

          p_xmom(istr,i_lay,ib) = (1.0/(z_scatliq+z_scatice))* &
           & (z_scatliq*(gliq(ib)**istr - forwliq(ib)) / &
           & (1. - forwliq(ib)) &
           & + z_scatice * &
           & ((gice(ib)-forwice(ib))/(1.0-forwice(ib)))**istr)  
        ENDDO
      ENDDO

    ENDIF

  ENDIF

ENDDO

IF (i_ndbug <= 1) THEN
  print *,'about to leave SRTM_CLDPROP'
ENDIF

!-----------------------------------------------------------------------
IF (lhook) CALL dr_hook('SRTM_CLDPROP',1,zhook_handle)
END SUBROUTINE srtm_cldprop