*IF DEF,A03_3A,OR,DEF,A03_5A                                               AJS1F401.1480   
C ******************************COPYRIGHT******************************    GTS2F400.8875   
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.    GTS2F400.8876   
C                                                                          GTS2F400.8877   
C Use, duplication or disclosure of this code is subject to the            GTS2F400.8878   
C restrictions as set forth in the contract.                               GTS2F400.8879   
C                                                                          GTS2F400.8880   
C                Meteorological Office                                     GTS2F400.8881   
C                London Road                                               GTS2F400.8882   
C                BRACKNELL                                                 GTS2F400.8883   
C                Berkshire UK                                              GTS2F400.8884   
C                RG12 2SZ                                                  GTS2F400.8885   
C                                                                          GTS2F400.8886   
C If no contract has been raised with this copy of the code, the use,      GTS2F400.8887   
C duplication or disclosure of it is strictly prohibited.  Permission      GTS2F400.8888   
C to do so must first be obtained in writing from the Head of Numerical    GTS2F400.8889   
C Modelling at the above address.                                          GTS2F400.8890   
C ******************************COPYRIGHT******************************    GTS2F400.8891   
C                                                                          GTS2F400.8892   
CLL  SUBROUTINE SFL_INT------------------------------------------------    SFLINT3A.3      
CLL                                                                        SFLINT3A.4      
CLL  Purpose: To calculate interpolation coefficients for 10m winds        SFLINT3A.5      
CLL           and 1.5m temperature/specific humidity diagnostics           SFLINT3A.6      
CLL           using a generalisation of the method of Geleyn (1988).       SFLINT3A.7      
CLL                                                                        SFLINT3A.8      
CLL  Suitable for single column use (via *IF definition IBM).              SFLINT3A.9      
CLL                                                                        SFLINT3A.10     
CLL  Model            Modification history:                                SFLINT3A.11     
CLL version  Date                                                          SFLINT3A.12     
CLL                                                                        SFLINT3A.13     
CLL   3.4  18/10/94   *DECK inserted into UM version 3.4. S Jackson        SFLINT3A.14     
CLL                                                                        SFLINT3A.15     
CLL   4.0  30/12/94   Modified calculation of 10m wind interpolation       ARN1F400.83     
CLL                   factor when effective roughness length used;         ARN1F400.84     
CLL                   10m wind assumed to lie on "local" profile at        ARN1F400.85     
CLL                   height z0m+10 metres above the surface.              ARN1F400.86     
CLL                                                    R.N.B.Smith         ARN1F400.87     
CLL   4.2   Oct. 96   T3E migration - *DEF CRAY removed                    GSS1F402.81     
CLL                                     S J Swarbrick                      GSS1F402.82     
CLL                                                                        ARN1F400.88     
CLL  Programming standard: Unified Model Documentation Paper No 4,         SFLINT3A.16     
CLL                        Version 2, dated 18/1/90.                       SFLINT3A.17     
CLL                                                                        SFLINT3A.18     
CLL  Logical component covered: Part of P243.                              SFLINT3A.19     
CLL                                                                        SFLINT3A.20     
CLL  System Task:                                                          SFLINT3A.21     
CLL                                                                        SFLINT3A.22     
CLL  External Documentation: UMDP No.24                                    SFLINT3A.23     
CLL                                                                        SFLINT3A.24     
CLL---------------------------------------------------------------------   SFLINT3A.25     
C*L  Arguments :-                                                          SFLINT3A.26     

      SUBROUTINE SFL_INT (                                                  3,8SFLINT3A.27     
     & P_POINTS,U_POINTS,RIB,Z1,Z0M,Z0M_EFF,Z0H,Z0F,CD_STD,CD,CH           ARN1F400.89     
     &,RESFT,WIND_PROFILE_FACTOR                                           ARN1F400.90     
     &,CDR10M,CHR1P5M,CER1P5M                                              ARN1F400.91     
     +,SU10,SV10,ST1P5,SQ1P5,LTIMER                                        SFLINT3A.30     
     +)                                                                    SFLINT3A.31     
      IMPLICIT NONE                                                        SFLINT3A.32     
      INTEGER                                                              SFLINT3A.33     
     + P_POINTS    ! IN No. of P-grid points to be processed.              SFLINT3A.34     
     +,U_POINTS    ! IN No. of UV-grid points to be processed.             SFLINT3A.35     
      REAL                                                                 SFLINT3A.36     
     + RIB(P_POINTS)     ! IN Bulk Richardson number for                   SFLINT3A.37     
C                        !    lowest layer.                                SFLINT3A.38     
     +,Z1(P_POINTS)      ! IN Height of lowest model level (m).            SFLINT3A.39     
     +,Z0M(P_POINTS)     ! IN Roughness length for momentum (m).           SFLINT3A.40     
     +,Z0M_EFF(P_POINTS) ! IN Effective roughness length for               SFLINT3A.41     
C                        !    momentum (m).                                SFLINT3A.42     
     +,Z0H(P_POINTS)     ! IN Roughness length for heat and                SFLINT3A.43     
C                        !    moisture (m).                                SFLINT3A.44     
     +,Z0F(P_POINTS)     ! IN Roughness length in the free                 SFLINT3A.45     
C                        !    convective limit (m).                        SFLINT3A.46     
     &,CD_STD(P_POINTS)  ! IN Surface drag coefficient for shear stress    ARN1F400.92     
C                        !    only, i.e. without orographic part of drag   ARN1F400.93     
     &,CD(P_POINTS)      ! IN Effective surface drag coefficient,          ARN1F400.94     
C                        !    including orographic part of drag            ARN1F400.95     
     &,CH(P_POINTS)      ! IN Surface transfer coefficient for heat and    ARN1F400.96     
C                        !    moisture.                                    SFLINT3A.49     
     +,RESFT(P_POINTS)   ! IN Total resistance factor for moisture         SFLINT3A.50     
C                        !    transfer from the surface.                   SFLINT3A.51     
     &,WIND_PROFILE_FACTOR(P_POINTS)                                       ARN1F400.97     
C                        ! IN Factor for converting effective friction     ARN1F400.98     
C                        !    velocity to local one.                       ARN1F400.99     
      LOGICAL                                                              SFLINT3A.52     
     + SU10                      ! IN 10m U-wind diagnostic flag           SFLINT3A.53     
     +,SV10                      ! IN 10m V-wind diagnostic flag           SFLINT3A.54     
     +,ST1P5                     ! IN screen temp diagnostic flag          SFLINT3A.55     
     +,SQ1P5                     ! IN screen specific humidity             SFLINT3A.56     
C                                !    diagnostic flag                      SFLINT3A.57     
     +,LTIMER                    ! IN TIMER diagnostics flag               SFLINT3A.58     
C Output variables                                                         SFLINT3A.59     
C                                                                          SFLINT3A.60     
      REAL                                                                 SFLINT3A.61     
     + CDR10M(U_POINTS)  ! OUT interpolation coefficicent for 10m wind     SFLINT3A.62     
     +,CHR1P5M(P_POINTS) ! OUT Interpolation coefficient for 1.5m          SFLINT3A.63     
C                        !     temperature                                 SFLINT3A.64     
     +,CER1P5M(P_POINTS) ! OUT Interpolation coefficient for 1.5m          SFLINT3A.65     
C                        !     specific humidity                           SFLINT3A.66     
C*                                                                         SFLINT3A.67     
C*L---------------------------------------------------------------------   SFLINT3A.68     
      EXTERNAL TIMER                                                       SFLINT3A.69     
C*                                                                         SFLINT3A.70     
C*L---------------------------------------------------------------------   SFLINT3A.71     
C    Local and other symbolic constants :-                                 SFLINT3A.72     
*CALL C_VKMAN                                                              SFLINT3A.73     
      REAL Z1P5M,Z10M                                                      SFLINT3A.74     
      PARAMETER (                                                          SFLINT3A.75     
     + Z1P5M = 1.5  ! for diagnosis of screen values of temperature        SFLINT3A.76     
C                   ! and humidity (assumed to be at 1.5m).                SFLINT3A.77     
     +,Z10M = 10.0  ! for diagnosis of 10m winds.                          SFLINT3A.78     
     +)                                                                    SFLINT3A.79     
C                                                                          SFLINT3A.80     
C  Define local storage.                                                   SFLINT3A.81     
C                                                                          SFLINT3A.82     
C  (a) Local work arrays.                                                  SFLINT3A.83     
C                                                                          SFLINT3A.84     
      REAL                                                                 SFLINT3A.86     
     + Z1E(P_POINTS)    ! Level 1 height + Z0M_EFF                         SFLINT3A.87     
     +,SQRTCD(P_POINTS) ! Square root of CD                                SFLINT3A.88     
C                                                                          SFLINT3A.95     
C  (b) Scalars.                                                            SFLINT3A.96     
C                                                                          SFLINT3A.97     
      REAL                                                                 SFLINT3A.98     
     + Z1S              ! Level 1 height + Z0M_EFF - Z0H                   SFLINT3A.99     
     +,Z1P5ME           ! Z1P5M + Z0H                                      ARN1F400.100    
     +,Z10ME            ! Z10M + Z0M                                       ARN1F400.101    
     +,SQRTCD_K         ! Temporary storage in calc of 1.5 amd 10m diags   SFLINT3A.102    
     +,Z_OVER_Z1        ! Temporary storage in calc of 1.5 amd 10m diags   SFLINT3A.103    
     +,CDNZ             ! Neutral drag coef. for momentum @ 10m            SFLINT3A.104    
     +,CDNZ1            ! Neutral drag coef. for momentum @ level1         SFLINT3A.105    
     +,CHNZ             ! Neutral drag coef. for heat/moisture @ 1.5m      SFLINT3A.106    
     +,CHNZ1            ! Neutral drag coef. for heat/moisture @ level1    SFLINT3A.107    
     &,CDTEMP1          ! Workspace in calc of interpolation coeffs.       ARN1F400.102    
     &,CDTEMP2          ! Workspace in calc of interpolation coeffs.       ARN1F400.103    
     &,CDTEMP3          ! Workspace in calc of interpolation coeffs.       ARN1F400.104    
                                                                           SFLINT3A.109    
      INTEGER                                                              SFLINT3A.110    
     + I       ! Loop counter (horizontal field index).                    SFLINT3A.111    
C*                                                                         SFLINT3A.112    
      IF (LTIMER) THEN                                                     SFLINT3A.113    
        CALL TIMER('SFL_INT   ',3)                                         SFLINT3A.114    
      ENDIF                                                                SFLINT3A.115    
C                                                                          SFLINT3A.116    
C-----------------------------------------------------------------------   SFLINT3A.117    
CL 1. This routine uses a generalised formulation of Geleyn (1988)         SFLINT3A.118    
CL   to interpolate to screen and 10m height using surface and bottom      SFLINT3A.119    
CL   level values as well as roughness lengths for heat and momentum.      SFLINT3A.120    
CL    Start of main loop. Set up variables needed later (eg height of      SFLINT3A.121    
CL   envelope orography).                                                  SFLINT3A.122    
C-----------------------------------------------------------------------   SFLINT3A.123    
                                                                           SFLINT3A.124    
      DO I=1,P_POINTS                                                      SFLINT3A.125    
        SQRTCD(I) = SQRT(CD(I))                                            SFLINT3A.126    
        Z1E(I) = Z1(I) + Z0M_EFF(I)                                        SFLINT3A.127    
      ENDDO                                                                SFLINT3A.128    
                                                                           SFLINT3A.129    
C-----------------------------------------------------------------------   SFLINT3A.130    
CL 2. If selected calculate interpolation coefficient for 10m winds.       SFLINT3A.131    
C-----------------------------------------------------------------------   SFLINT3A.132    
C                                                                          SFLINT3A.133    
      IF(SU10.OR.SV10) THEN                                                SFLINT3A.134    
        DO I=1,P_POINTS                                                    SFLINT3A.135    
          Z10ME = Z10M + Z0M(I)                                            SFLINT3A.136    
          CDNZ = LOG( Z10ME / Z0M(I) )                                     SFLINT3A.137    
          CDNZ1 = LOG( Z1E(I) / Z0M_EFF(I) )                               ARN1F400.105    
          SQRTCD_K = SQRTCD(I) / VKMAN                                     SFLINT3A.139    
          Z_OVER_Z1 = Z10M  / Z1(I)                                        ARN1F400.106    
                                                                           SFLINT3A.141    
          IF (RIB(I).GE.0.0) THEN                                          SFLINT3A.142    
C                                                                          ARN1F400.107    
C Stable case                                                              ARN1F400.108    
C                                                                          ARN1F400.109    
            CDR10M(I) = Z_OVER_Z1 + SQRTCD_K *                             SFLINT3A.143    
     +                  (CDNZ - Z_OVER_Z1*CDNZ1)                           SFLINT3A.144    
          ELSE                                                             SFLINT3A.145    
C                                                                          ARN1F400.110    
C Unstable Case                                                            ARN1F400.111    
C                                                                          ARN1F400.112    
            CDTEMP1 = EXP( 1.0 / SQRTCD_K )                                ARN1F400.113    
            CDTEMP2 = Z1E(I) / Z0M_EFF(I)                                  ARN1F400.114    
            CDTEMP3 = LOG( ( ( Z1E(I) - Z0M(I) ) * CDTEMP1 -               ARN1F400.115    
     &                       ( Z0M_EFF(I) - Z0M(I) ) * CDTEMP2 ) /         ARN1F400.116    
     &                     ( ( Z1E(I) - Z10ME ) * CDTEMP1 -                ARN1F400.117    
     &                       ( Z0M_EFF(I) - Z10ME ) * CDTEMP2 ) )          ARN1F400.118    
            CDR10M(I) = SQRTCD_K * ( CDNZ + CDTEMP3 )                      ARN1F400.119    
          ENDIF                                                            SFLINT3A.150    
          CDR10M(I) = CDR10M(I) * WIND_PROFILE_FACTOR(I)                   ARN1F400.120    
        ENDDO                                                              SFLINT3A.151    
      ENDIF                                                                SFLINT3A.152    
C                                                                          SFLINT3A.153    
C-----------------------------------------------------------------------   SFLINT3A.154    
CL 3. If selected calculate interpolation coefficient for 1.5m screen      SFLINT3A.155    
CL      temperature and specific humidity.                                 SFLINT3A.156    
C-----------------------------------------------------------------------   SFLINT3A.157    
C                                                                          SFLINT3A.158    
      IF (ST1P5.OR.SQ1P5) THEN                                             SFLINT3A.159    
        DO I=1,P_POINTS                                                    SFLINT3A.160    
C                                                                          SFLINT3A.161    
C variables to be used later                                               SFLINT3A.162    
          Z1S = Z1E(I) - Z0H(I)                                            SFLINT3A.163    
          Z1P5ME = Z1P5M + Z0H(I)                                          SFLINT3A.164    
          CHNZ = LOG( Z1P5ME / Z0H(I) )                                    SFLINT3A.165    
          CHNZ1 = LOG( Z1E(I) / Z0H(I) )                                   SFLINT3A.166    
          SQRTCD(I) = SQRT(CD_STD(I))                                      ARN1F400.121    
          SQRTCD_K =0.0                                                    SFLINT3A.167    
          IF (SQRTCD(I) .GT. 0.0) SQRTCD_K = CH(I) / (SQRTCD(I) * VKMAN)   SFLINT3A.168    
          Z_OVER_Z1 =  Z1P5M  / Z1S                                        SFLINT3A.169    
                                                                           SFLINT3A.170    
C                                                                          SFLINT3A.171    
C Stable case                                                              SFLINT3A.172    
C                                                                          ARN1F400.122    
          IF (RIB(I).GE.0.0) THEN                                          SFLINT3A.173    
            CHR1P5M(I) = Z_OVER_Z1 + SQRTCD_K *                            SFLINT3A.174    
     +                  (CHNZ - Z_OVER_Z1 * CHNZ1)                         SFLINT3A.175    
          ELSE                                                             SFLINT3A.176    
C                                                                          SFLINT3A.177    
C Unstable Case                                                            SFLINT3A.178    
C                                                                          ARN1F400.123    
            CDTEMP1 = EXP(1.0 / SQRTCD_K)                                  ARN1F400.124    
            CDTEMP2 = ( Z1P5M * Z1E(I) +                                   ARN1F400.125    
     &                  Z0H(I) * ( Z1S - Z1P5M ) * CDTEMP1 ) /             ARN1F400.126    
     &                ( Z1S * Z1P5ME )                                     ARN1F400.127    
            CHR1P5M(I) = 1.0 - SQRTCD_K * LOG ( CDTEMP2 )                  ARN1F400.128    
          ENDIF                                                            SFLINT3A.183    
                                                                           SFLINT3A.184    
C                                                                          SFLINT3A.185    
            CER1P5M(I) = ( CHR1P5M(I) - 1.0 ) * RESFT(I)      ! P243.123   SFLINT3A.186    
        ENDDO                                                              SFLINT3A.187    
      ENDIF                                                                SFLINT3A.188    
C                                                                          SFLINT3A.189    
      IF (LTIMER) THEN                                                     SFLINT3A.190    
        CALL TIMER('SFL_INT ',4)                                           SFLINT3A.191    
      ENDIF                                                                SFLINT3A.192    
      RETURN                                                               SFLINT3A.193    
      END                                                                  SFLINT3A.194    
*ENDIF                                                                     SFLINT3A.195