!+ Subroutine to calculate optical properties of ice clouds.
!
! Method:
!	If the optical properties come from an observational 
!	distribution a separte subroutine is called. Otherwise
!	appropriate mean quantities in the layer are calculated
!	as the parametrization requires and these values are
!	substituted into the parametrization to give the optical
!	properties.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE OPT_PROP_ICE_CLOUD(IERR
     &   , N_PROFILE, N_LAYER, N_CLOUD_TOP
     &   , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &   , L_RESCALE, L_LAYER, L_CLOUD_LAYER
     &   , I_PARAMETRIZATION_ICE, ICE_CLOUD_PARAMETER
     &   , ICE_MASS_FRAC, RADIUS_EFFECT_ICE
     &   , P, T, DENSITY
     &   , N_LAYER_ICE_DATA, ICE_PRESSURE
     &   , ICE_ABSORPTION, ICE_SCATTERING
     &   , ICE_ASYMMETRY
     &   , K_EXT_TOT_CLOUD, K_EXT_SCAT_CLOUD
     &   , ASYMMETRY_CLOUD, FORWARD_SCATTER_CLOUD
     &   )
!
!
      IMPLICIT NONE
!
!
!     INCLUSION OF HEADER FILES.
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_field.ucf'
      INCLUDE 'incl/ice_cloud_parametrization.pcf'
      INCLUDE 'incl/layer_distrib.ucf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY VARIABLES.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_PROFILE
!             NUMBER OF PROFILES
     &   , N_LAYER
!             NUMBER OF LAYERS
     &   , N_CLOUD_TOP
!             TOPMOST CLOUDY LAYER
     &   , I_PARAMETRIZATION_ICE
!             TREATMENT OF ICE CRYSTALS
     &   , N_LAYER_ICE_DATA
!             NUMBER OF ICE CLOUD DATA LAYERS
     &   , N_CLOUD_PROFILE(NPD_LAYER)
!             NUMBER OF CLOUDY PROFILES
     &   , I_CLOUD_PROFILE(NPD_PROFILE, NPD_LAYER)
!             PROFILES CONTAINING CLOUDS
      LOGICAL	!, INTENT(IN)
     &     L_LAYER
!             VARIABLES GIVEN IN LAYERS
     &   , L_CLOUD_LAYER
!             CLOUD VARIABLES GIVEN IN LAYERS
     &   , L_RESCALE
!             DELTA-RESCALING REQUIRED
      REAL*8	!, INTENT(IN)
     &     ICE_CLOUD_PARAMETER(NPD_CLOUD_PARAMETER)
!             ICE CLOUD PARAMETERS
     &   , ICE_MASS_FRAC(NPD_PROFILE, 0: NPD_LAYER)
!             ICE MASS FRACTION
     &   , RADIUS_EFFECT_ICE(NPD_PROFILE, 0: NPD_LAYER)
!             EFFECTIVE RADIUS OF ICE CRYSTALS
     &   , T(NPD_PROFILE, 0: NPD_LAYER)
!             TEMPERATURE
     &   , DENSITY(NPD_PROFILE, 0: NPD_LAYER)
!             DENSITY AT LEVELS
      REAL*8	!, INTENT(IN)
     &     P(NPD_PROFILE, 0: NPD_LAYER)
!             PRESSURE
     &   , ICE_PRESSURE(0: NPD_CLOUD_DATA_LAYER)
!             PRESSURE AT ICE CLOUD DATA LEVELS
     &   , ICE_ABSORPTION(0: NPD_CLOUD_DATA_LAYER)
!             ICE CLOUD ABSORPTION
     &   , ICE_SCATTERING(0: NPD_CLOUD_DATA_LAYER)
!             ICE CLOUD SCATTERING
     &   , ICE_ASYMMETRY(0: NPD_CLOUD_DATA_LAYER)
!             ICE CLOUD ASYMMETRY
      REAL*8	!, INTENT(OUT)
     &     K_EXT_SCAT_CLOUD(NPD_PROFILE, NPD_LAYER)
!             SCATTERING EXTINCTION
     &   , K_EXT_TOT_CLOUD(NPD_PROFILE, NPD_LAYER)
!             TOTAL EXTINCTION
     &   , ASYMMETRY_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY ASYMMETRIES
     &   , FORWARD_SCATTER_CLOUD(NPD_PROFILE, NPD_LAYER)
!             CLOUDY FORWARD SCATTERING
!
!     LOCAL VARIABLES.
      INTEGER
     &     L
!             LOOP VARIABLE
     &   , LL
!             LOOP VARIABLE
     &   , I
!             LOOP VARIABLE
      REAL*8
     &     ASYMMETRY_PROCESS(NPD_PROFILE)
!             ASYMMETRY FACTOR FOR CURRENT PROC.
     &   , RADIUS_AVE
!             AVERAGE EFFECTIVE RADIUS IN LAYER
     &   , ICE_MASS_FRAC_AVE
!             AVERAGE ICE MASS FRACTION
     &   , DENSITY_AVE
!             AVERAGE DENSITY
     &   , V_M
!             PARAMETER IN CIRRUS ADT SCHEME
     &   , OMEGA_0
!             S. S. ALBEDO
     &   , T_CELSIUS
!             TEMPERATURE IN CELSIUS
     &   , TEMP_CORRECTION
!             TEMPERATURE CORRECTION
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     OBS_CLOUD_PROP
!
!
!
      IF (I_PARAMETRIZATION_ICE.EQ.IP_SLINGO_SCHRECKER_ICE) THEN
!
!        WE CALCULATE AVERAGE PROPERTIES FOR THE LAYER AND PUT THESE
!        INTO THE PARAMETRIZATION, RATHER THAN CALCULATING THE
!        PARAMETRIZATION AT EACH LEVEL: USUALLY THIS IS MORE ACCURATE.
!        IT ALSO FITS MORE NATURALLY WITH CASES WHERE DATA ARE GIVEN
!        IN LAYERS.
!
!        THE TOTAL EXTINCTIONS SHOULD BE INCREMENTED BY THE TOTAL 
!        CONTRIBUTIONS FROM CLOUDS, NOT JUST BY THE ABSORPTIVE 
!        EXTINCTIONS.
!
         DO I=N_CLOUD_TOP, N_LAYER
            DO LL=1, N_CLOUD_PROFILE(I)
               L=I_CLOUD_PROFILE(LL, I)
               IF (L_CLOUD_LAYER) THEN
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
                  RADIUS_AVE=RADIUS_EFFECT_ICE(L, I)
               ELSE
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
     &              +WEIGHT_UPPER_CLOUD
     &              *(ICE_MASS_FRAC(L, I-1)
     &              -ICE_MASS_FRAC(L, I))
                  RADIUS_AVE=RADIUS_EFFECT_ICE(L, I)
     &               +WEIGHT_UPPER_CLOUD
     &               *(RADIUS_EFFECT_ICE(L, I-1)
     &               -RADIUS_EFFECT_ICE(L, I))
               ENDIF
               K_EXT_TOT_CLOUD(L, I)
     &            =ICE_MASS_FRAC_AVE*(ICE_CLOUD_PARAMETER(1)
     &            +ICE_CLOUD_PARAMETER(2)/RADIUS_AVE)
               K_EXT_SCAT_CLOUD(L, I)=K_EXT_TOT_CLOUD(L, I)
     &            *(1.0D+00-ICE_CLOUD_PARAMETER(3)
     &            -ICE_CLOUD_PARAMETER(4)*RADIUS_AVE)
               ASYMMETRY_PROCESS(L)=
     &            ICE_CLOUD_PARAMETER(5)+ICE_CLOUD_PARAMETER(6)
     &            *RADIUS_AVE
               ASYMMETRY_CLOUD(L, I)=
     &            K_EXT_SCAT_CLOUD(L, I)*ASYMMETRY_PROCESS(L)
            ENDDO
            IF (L_RESCALE) THEN
               DO LL=1, N_CLOUD_PROFILE(I)
                  L=I_CLOUD_PROFILE(LL, I)
                  FORWARD_SCATTER_CLOUD(L, I)
     &               =K_EXT_SCAT_CLOUD(L, I)
     &               *ASYMMETRY_PROCESS(L)**2
               ENDDO
            ENDIF
         ENDDO
!
      ELSE IF (I_PARAMETRIZATION_ICE.EQ.IP_CIRRUS_ADT) THEN
!
         DO I=N_CLOUD_TOP, N_LAYER
            DO LL=1, N_CLOUD_PROFILE(I)
               L=I_CLOUD_PROFILE(LL, I)
               IF (L_CLOUD_LAYER) THEN
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
                  RADIUS_AVE=RADIUS_EFFECT_ICE(L, I)
               ELSE
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
     &              +WEIGHT_UPPER_CLOUD
     &              *(ICE_MASS_FRAC(L, I-1)
     &              -ICE_MASS_FRAC(L, I))
                  RADIUS_AVE=RADIUS_EFFECT_ICE(L, I)
     &               +WEIGHT_UPPER_CLOUD
     &               *(RADIUS_EFFECT_ICE(L, I-1)
     &               -RADIUS_EFFECT_ICE(L, I))
               ENDIF
!              A NOMINAL DENISTY OF 1.0D+03 IS USED FOR ICE. THIS
!              CANCELS OUT IF IT HAS BEEN SPECIFIED IN THE AVERAGING
!              CODE: IT IS ALSO IMPLICIT IN THE WAY IN WHICH IWC
!              AND RE ARE DEFINED.
               K_EXT_TOT_CLOUD(L, I)=1.5D+00*ICE_MASS_FRAC_AVE
     &            /(1.0D+03*RADIUS_AVE)
               V_M=ICE_CLOUD_PARAMETER(1)*RADIUS_AVE
               OMEGA_0=0.5D+00+(1.0D+00/6.0D+00)
     &            *(V_M+3.0D+00)/(V_M+1.0D+00)**3
               K_EXT_SCAT_CLOUD(L, I)=K_EXT_TOT_CLOUD(L, I)*OMEGA_0
               ASYMMETRY_CLOUD(L, I)=
     &            +K_EXT_SCAT_CLOUD(L, I)*ICE_CLOUD_PARAMETER(2)
            ENDDO
            IF (L_RESCALE) THEN
               DO LL=1, N_CLOUD_PROFILE(I)
                  L=I_CLOUD_PROFILE(LL, I)
                  FORWARD_SCATTER_CLOUD(L, I)
     &               =K_EXT_SCAT_CLOUD(L, I)*ICE_CLOUD_PARAMETER(2)**2
               ENDDO
            ENDIF
         ENDDO
!
      ELSE IF (I_PARAMETRIZATION_ICE.EQ.IP_ICE_UNPARAMETRIZED) THEN
!
         CALL OBS_CLOUD_PROP(IERR
     &      , N_PROFILE, N_LAYER, N_LAYER_ICE_DATA
     &      , L_RESCALE
     &      , N_CLOUD_PROFILE, I_CLOUD_PROFILE
     &      , P, DENSITY, L_LAYER
     &      , ICE_PRESSURE, ICE_ABSORPTION, ICE_SCATTERING
     &      , ICE_ASYMMETRY
     &      , K_EXT_TOT_CLOUD, K_EXT_SCAT_CLOUD, ASYMMETRY_CLOUD
     &      , FORWARD_SCATTER_CLOUD
     &      )
!
!
      ELSE IF (I_PARAMETRIZATION_ICE.EQ.IP_SUN_SHINE_VN2_VIS) THEN
!
         DO I=N_CLOUD_TOP, N_LAYER
            DO LL=1, N_CLOUD_PROFILE(I)
               L=I_CLOUD_PROFILE(LL, I)
               IF (L_CLOUD_LAYER) THEN
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
                  DENSITY_AVE=DENSITY(L, I)
               ELSE
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
     &              +WEIGHT_UPPER_CLOUD
     &              *(ICE_MASS_FRAC(L, I-1)-ICE_MASS_FRAC(L, I))
                  DENSITY_AVE=DENSITY(L, I)
     &               +WEIGHT_UPPER_CLOUD
     &               *(DENSITY(L, I-1)-DENSITY(L, I))
               ENDIF
               IF (L_LAYER) THEN
                  T_CELSIUS=T(L, I)-273.16D+00
               ELSE
                  T_CELSIUS=T(L, I)
     &               +WEIGHT_UPPER_CLOUD*(T(L, I-1)-T(L, I))
     &               -273.16D+00
               ENDIF
               TEMP_CORRECTION=1.047D+00+T_CELSIUS*(-9.13D-05+T_CELSIUS
     &            *(2.026D-04-1.056D-05*T_CELSIUS))
               K_EXT_TOT_CLOUD(L, I)=TEMP_CORRECTION*ICE_MASS_FRAC_AVE
     &            /(3.05548D-02
     &            +2.54802D+02*DENSITY_AVE*ICE_MASS_FRAC_AVE)
               K_EXT_SCAT_CLOUD(L, I)=K_EXT_TOT_CLOUD(L, I)
     &            *(1.0D+00-ICE_CLOUD_PARAMETER(1)
     &            *EXP(ICE_CLOUD_PARAMETER(2)
     &            *LOG(DENSITY_AVE*ICE_MASS_FRAC_AVE+1.0D-12)))
     &            *(1.0D+00+ICE_CLOUD_PARAMETER(5)
     &            *(TEMP_CORRECTION-1.0D+00)/TEMP_CORRECTION)
               ASYMMETRY_PROCESS(L)=
     &            ICE_CLOUD_PARAMETER(3)*EXP(ICE_CLOUD_PARAMETER(4)
     &            *LOG(DENSITY_AVE*ICE_MASS_FRAC_AVE+1.0D-12))
     &            *(1.0D+00+ICE_CLOUD_PARAMETER(6)
     &            *(TEMP_CORRECTION-1.0D+00)/TEMP_CORRECTION)
               ASYMMETRY_CLOUD(L, I)=
     &            K_EXT_SCAT_CLOUD(L, I)*ASYMMETRY_PROCESS(L)
            ENDDO
            IF (L_RESCALE) THEN
               DO LL=1, N_CLOUD_PROFILE(I)
                  L=I_CLOUD_PROFILE(LL, I)
                  FORWARD_SCATTER_CLOUD(L, I)
     &               =K_EXT_SCAT_CLOUD(L, I)
     &               *ASYMMETRY_PROCESS(L)**2
               ENDDO
            ENDIF
         ENDDO
!
      ELSE IF (I_PARAMETRIZATION_ICE.EQ.IP_SUN_SHINE_VN2_IR) THEN
!
         DO I=N_CLOUD_TOP, N_LAYER
            DO LL=1, N_CLOUD_PROFILE(I)
               L=I_CLOUD_PROFILE(LL, I)
               IF (L_CLOUD_LAYER) THEN
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
                  DENSITY_AVE=DENSITY(L, I)
               ELSE
                  ICE_MASS_FRAC_AVE=ICE_MASS_FRAC(L, I)
     &              +WEIGHT_UPPER_CLOUD
     &              *(ICE_MASS_FRAC(L, I-1)
     &              -ICE_MASS_FRAC(L, I))
                  DENSITY_AVE=DENSITY(L, I)
     &              +WEIGHT_UPPER_CLOUD
     &              *(DENSITY(L, I-1)
     &              -DENSITY(L, I))
               ENDIF
               IF (L_LAYER) THEN
                  T_CELSIUS=T(L, I)-273.16D+00
               ELSE
                  T_CELSIUS=T(L, I)
     &               +WEIGHT_UPPER_CLOUD*(T(L, I-1)-T(L, I))
     &               -273.16D+00
               ENDIF
               TEMP_CORRECTION=1.047D+00+T_CELSIUS*(-9.13D-05+T_CELSIUS
     &            *(2.026D-04-1.056D-05*T_CELSIUS))
               K_EXT_TOT_CLOUD(L, I)=TEMP_CORRECTION*ICE_MASS_FRAC_AVE
     &            /(6.30689D-02
     &            +2.65874D+02*DENSITY_AVE*ICE_MASS_FRAC_AVE)
            ENDDO
         ENDDO
!
      ELSE
         WRITE(IU_ERR, '(/A)') '*** ERROR: AN INVALID PARAMETRIZATION '
     &      //'OF ICE CRYSTALS HAS BEEN USED..'
         IERR=I_ERR_FATAL
         RETURN
      ENDIF
!
!
!
      RETURN
      END
