*IF DEF,A03_7A                                                             SFOROG7A.2      
C *****************************COPYRIGHT******************************     SFOROG7A.3      
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved.    SFOROG7A.4      
C                                                                          SFOROG7A.5      
C Use, duplication or disclosure of this code is subject to the            SFOROG7A.6      
C restrictions as set forth in the contract.                               SFOROG7A.7      
C                                                                          SFOROG7A.8      
C                Meteorological Office                                     SFOROG7A.9      
C                London Road                                               SFOROG7A.10     
C                BRACKNELL                                                 SFOROG7A.11     
C                Berkshire UK                                              SFOROG7A.12     
C                RG12 2SZ                                                  SFOROG7A.13     
C                                                                          SFOROG7A.14     
C If no contract has been raised with this copy of the code, the use,      SFOROG7A.15     
C duplication or disclosure of it is strictly prohibited.  Permission      SFOROG7A.16     
C to do so must first be obtained in writing from the Head of Numerical    SFOROG7A.17     
C Modelling at the above address.                                          SFOROG7A.18     
C ******************************COPYRIGHT******************************    SFOROG7A.19     
!                                                                          SFOROG7A.20     
!!!  SUBROUTINE SF_OROG------------------------------------------------    SFOROG7A.21     
!!!                                                                        SFOROG7A.22     
!!!  Purpose: Calculate roughness lengths, blending height and wind        SFOROG7A.23     
!!!           profile factor                                               SFOROG7A.24     
!!!                                                                        SFOROG7A.25     
!!! SJ, RE        <- programmer of some or all of previous code changes    SFOROG7A.26     
C Modification History:                                                    AJC1F405.63     
C Version Date     Change                                                  AJC1F405.64     
C  4.5    Jul. 98  Kill the IBM specific lines (JCThil)                    AJC1F405.65     
!!!                                                                        SFOROG7A.27     
!!!--------------------------------------------------------------------    SFOROG7A.28     

      SUBROUTINE SF_OROG(                                                   2,2SFOROG7A.29     
     & P_FIELD,LAND_FIELD,TILE_PTS,LAND_INDEX,TILE_INDEX,                  SFOROG7A.30     
     & L_Z0_OROG,LTIMER,                                                   SFOROG7A.31     
     & HO2R2_OROG,RIB,SIL_OROG,Z0M,Z1,                                     SFOROG7A.32     
     & WIND_PROFILE_FACTOR,Z0M_EFF                                         SFOROG7A.33     
     & )                                                                   SFOROG7A.34     
                                                                           SFOROG7A.35     
      IMPLICIT NONE                                                        SFOROG7A.36     
                                                                           SFOROG7A.37     
      INTEGER                                                              SFOROG7A.38     
     & P_FIELD              ! IN Number of points on P-grid.               SFOROG7A.39     
     &,LAND_FIELD           ! IN Number of land points.                    SFOROG7A.40     
     &,TILE_PTS             ! IN Number of tile points.                    SFOROG7A.41     
     &,LAND_INDEX(P_FIELD)  ! IN Index of land points.                     SFOROG7A.42     
     &,TILE_INDEX(LAND_FIELD)! IN Index of tile points.                    SFOROG7A.43     
                                                                           SFOROG7A.44     
      LOGICAL                                                              SFOROG7A.45     
     & LTIMER               ! IN .TRUE. for timer diagnostics              SFOROG7A.46     
     &,L_Z0_OROG            ! IN .TRUE. to use orographic roughness.       SFOROG7A.47     
                                                                           SFOROG7A.48     
      REAL                                                                 SFOROG7A.49     
     & HO2R2_OROG(LAND_FIELD)!IN Peak to trough height of unresolved       SFOROG7A.50     
!                            !   orography divided by 2SQRT(2) (m).        SFOROG7A.51     
     &,RIB(LAND_FIELD)      ! IN Bulk Richardson number for lowest layer   SFOROG7A.52     
     &,SIL_OROG(LAND_FIELD) ! IN Silhouette area of unresolved orography   SFOROG7A.53     
!                           !    per unit horizontal area                  SFOROG7A.54     
     &,Z0M(LAND_FIELD)      ! IN Roughness length for momentum (m).        SFOROG7A.55     
     &,Z1(P_FIELD)          ! IN Height of lowest atmospheric level (m).   SFOROG7A.56     
                                                                           SFOROG7A.57     
!  Output variables.                                                       SFOROG7A.58     
                                                                           SFOROG7A.59     
      REAL                                                                 SFOROG7A.60     
     & WIND_PROFILE_FACTOR(LAND_FIELD)                                     SFOROG7A.61     
!                           ! OUT For transforming effective surface       SFOROG7A.62     
!                           !     transfer coefficients to those           SFOROG7A.63     
!                           !     excluding form drag.                     SFOROG7A.64     
     &,Z0M_EFF(LAND_FIELD)  ! OUT Effective roughness length for           SFOROG7A.65     
!                                 momentum (m)                             SFOROG7A.66     
                                                                           SFOROG7A.67     
!  Work Variables                                                          SFOROG7A.68     
                                                                           SFOROG7A.69     
      INTEGER                                                              SFOROG7A.70     
     & I            ! Horizontal field index                               SFOROG7A.71     
     &,J            ! Tile field index                                     SFOROG7A.72     
     &,L            ! Land field index                                     SFOROG7A.73     
                                                                           SFOROG7A.74     
      REAL                                                                 SFOROG7A.75     
     & H_BLEND_OROG ! Blending height                                      SFOROG7A.76     
     &,RIB_FN       ! Interpolation coefficient for 0 < RIB < RI_CRIT      SFOROG7A.77     
     &,ZETA1        ! Work space                                           SFOROG7A.78     
     &,ZETA2        ! More work space                                      SFOROG7A.79     
     &,ZETA3        ! Even more work space                                 SFOROG7A.80     
                                                                           SFOROG7A.81     
!   Common parameters                                                      SFOROG7A.82     
                                                                           SFOROG7A.83     
*CALL C_MDI                                                                SFOROG7A.84     
*CALL C_VKMAN                                                              SFOROG7A.85     
*CALL C_SURF                                                               SFOROG7A.86     
                                                                           SFOROG7A.87     
!   Local parameters                                                       SFOROG7A.88     
                                                                           SFOROG7A.89     
      REAL H_BLEND_MIN,H_BLEND_MAX                                         SFOROG7A.90     
      PARAMETER (                                                          SFOROG7A.91     
     & H_BLEND_MIN=0.0       ! Minimun value of blending height            SFOROG7A.92     
     &,H_BLEND_MAX=1000.0    ! Maximum value of blending height            SFOROG7A.93     
     & )                                                                   SFOROG7A.94     
                                                                           SFOROG7A.95     
      EXTERNAL TIMER                                                       SFOROG7A.96     
                                                                           SFOROG7A.97     
      IF (LTIMER) THEN                                                     SFOROG7A.98     
        CALL TIMER('SF_OROG ',3)                                           SFOROG7A.99     
      ENDIF                                                                SFOROG7A.100    
                                                                           SFOROG7A.101    
! Set blending height, effective roughness length and                      SFOROG7A.102    
! wind profile factor at land points.                                      SFOROG7A.103    
      DO J=1,TILE_PTS                                                      SFOROG7A.104    
        L = TILE_INDEX(J)                                                  SFOROG7A.105    
        I = LAND_INDEX(L)                                                  SFOROG7A.106    
                                                                           SFOROG7A.107    
        WIND_PROFILE_FACTOR(L) = 1.0                                       SFOROG7A.108    
        Z0M_EFF(L) = Z0M(L)                                                SFOROG7A.109    
                                                                           SFOROG7A.110    
        IF (L_Z0_OROG) THEN                                                SFOROG7A.111    
                                                                           SFOROG7A.112    
          ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L)                      SFOROG7A.113    
          ZETA2 = LOG ( 1.0 + Z1(I)/Z0M(L) )                               SFOROG7A.114    
          ZETA3 = 1.0 / SQRT ( ZETA1/(VKMAN*VKMAN) + 1.0/(ZETA2*ZETA2) )   SFOROG7A.115    
          ZETA2 = 1.0 / EXP(ZETA3)                                         SFOROG7A.116    
          H_BLEND_OROG = MAX ( Z1(I) / (1.0 - ZETA2) ,                     SFOROG7A.117    
     &                       HO2R2_OROG(L) * SQRT(2.0) )                   SFOROG7A.118    
          H_BLEND_OROG = MIN ( H_BLEND_MAX, H_BLEND_OROG )                 SFOROG7A.119    
                                                                           SFOROG7A.120    
! Apply simple stability correction to form drag if RIB is stable          SFOROG7A.121    
                                                                           SFOROG7A.122    
          IF (SIL_OROG(L) .EQ. RMDI) THEN                                  SFOROG7A.123    
             ZETA1 = 0.0                                                   SFOROG7A.124    
          ELSE                                                             SFOROG7A.125    
             RIB_FN =  ( 1.0 - RIB(L) / RI_CRIT )                          SFOROG7A.126    
             IF (RIB_FN.GT.1.0) RIB_FN = 1.0                               SFOROG7A.127    
             IF (RIB_FN.LT.0.0) RIB_FN = 0.0                               SFOROG7A.128    
             ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) * RIB_FN          SFOROG7A.129    
          ENDIF                                                            SFOROG7A.130    
                                                                           SFOROG7A.131    
          Z0M_EFF(L) = H_BLEND_OROG / EXP ( VKMAN / SQRT ( ZETA1 +         SFOROG7A.132    
     &                 (VKMAN / LOG ( H_BLEND_OROG / Z0M(L) ) ) *          SFOROG7A.133    
     &                 (VKMAN / LOG ( H_BLEND_OROG / Z0M(L) ) ) ) )        SFOROG7A.134    
          IF (RIB(L).GT.RI_CRIT) Z0M_EFF(L)=Z0M(L)                         SFOROG7A.135    
                                                                           SFOROG7A.136    
          WIND_PROFILE_FACTOR(L) = LOG( H_BLEND_OROG / Z0M_EFF(L) ) /      SFOROG7A.137    
     &                             LOG( H_BLEND_OROG / Z0M(L) )            SFOROG7A.138    
                                                                           SFOROG7A.139    
        ENDIF                                                              SFOROG7A.140    
                                                                           SFOROG7A.141    
      ENDDO                                                                SFOROG7A.142    
                                                                           SFOROG7A.143    
      IF (LTIMER) THEN                                                     SFOROG7A.144    
        CALL TIMER('SF_OROG ',4)                                           SFOROG7A.145    
      ENDIF                                                                SFOROG7A.146    
                                                                           SFOROG7A.147    
      RETURN                                                               SFOROG7A.148    
      END                                                                  SFOROG7A.149    
                                                                           SFOROG7A.150    
!!!  SUBROUTINE SF_OROG_GB --------------------------------------------    SFOROG7A.151    
!!!                                                                        SFOROG7A.152    
!!!  Purpose: Calculate effective roughness length and blending height     SFOROG7A.153    
!!!                                                                        SFOROG7A.154    
!!! SJ, RE        <- programmer of some or all of previous code changes    SFOROG7A.155    
!!!                                                                        SFOROG7A.156    
!!!--------------------------------------------------------------------    SFOROG7A.157    

      SUBROUTINE SF_OROG_GB(                                                1,2SFOROG7A.158    
     & P_FIELD,P1,P_POINTS,LAND_FIELD,LAND1,LAND_PTS,LAND_INDEX,           SFOROG7A.159    
     & LAND_MASK,L_Z0_OROG,HO2R2_OROG,RIB,SIL_OROG,Z0M,Z1,                 SFOROG7A.160    
     & H_BLEND_OROG,Z0M_EFF,LTIMER                                         SFOROG7A.161    
     & )                                                                   SFOROG7A.162    
                                                                           SFOROG7A.163    
      IMPLICIT NONE                                                        SFOROG7A.164    
                                                                           SFOROG7A.165    
      INTEGER                                                              SFOROG7A.166    
     & P_FIELD              ! IN Number of points on P-grid.               SFOROG7A.167    
     &,P1                   ! IN First P-point to be processed.            SFOROG7A.168    
     &,P_POINTS             ! IN Number of P-grid points to be processed   SFOROG7A.169    
     &,LAND_FIELD           ! IN Number of land points.                    SFOROG7A.170    
     &,LAND1                ! IN First land point to be processed.         SFOROG7A.171    
     &,LAND_PTS             ! IN Number of land points to be processed.    SFOROG7A.172    
     &,LAND_INDEX(P_FIELD)  ! IN Index of land points.                     SFOROG7A.173    
                                                                           SFOROG7A.174    
      LOGICAL                                                              SFOROG7A.175    
     & LAND_MASK(P_FIELD)    ! IN .TRUE. for land; .FALSE. elsewhere.      SFOROG7A.176    
     &,LTIMER               ! IN .TRUE. for timer diagnostics              SFOROG7A.177    
     &,L_Z0_OROG            ! IN .TRUE. to use orographic roughness.       SFOROG7A.178    
                                                                           SFOROG7A.179    
                                                                           SFOROG7A.180    
      REAL                                                                 SFOROG7A.181    
     & HO2R2_OROG(LAND_FIELD)!IN Peak to trough height of unresolved       SFOROG7A.182    
!                            !   orography divided by 2SQRT(2) (m).        SFOROG7A.183    
     &,RIB(P_FIELD)         ! IN GBM Bulk Richardson number for lowest     SFOROG7A.184    
!                           !    layer                                     SFOROG7A.185    
     &,SIL_OROG(LAND_FIELD) ! IN Silhouette area of unresolved orography   SFOROG7A.186    
!                           !    per unit horizontal area                  SFOROG7A.187    
     &,Z0M(LAND_FIELD)      ! IN GBM Roughness length for momentum (m).    SFOROG7A.188    
     &,Z1(P_FIELD)          ! IN Height of lowest atmospheric level (m).   SFOROG7A.189    
                                                                           SFOROG7A.190    
!  Output variables.                                                       SFOROG7A.191    
                                                                           SFOROG7A.192    
      REAL                                                                 SFOROG7A.193    
     & H_BLEND_OROG(P_FIELD)! OUT Blending height                          SFOROG7A.194    
     &,Z0M_EFF(P_FIELD)     ! OUT Effective roughness length for           SFOROG7A.195    
!                                 momentum (m)                             SFOROG7A.196    
                                                                           SFOROG7A.197    
!  Work Variables                                                          SFOROG7A.198    
                                                                           SFOROG7A.199    
      INTEGER                                                              SFOROG7A.200    
     & I            ! Horizontal field index                               SFOROG7A.201    
     &,L            ! Land field index                                     SFOROG7A.202    
                                                                           SFOROG7A.203    
      REAL                                                                 SFOROG7A.204    
     & RIB_FN       ! Interpolation coefficient for 0 < RIB < RI_CRIT      SFOROG7A.205    
     &,ZETA1        ! Work space                                           SFOROG7A.206    
     &,ZETA2        ! More work space                                      SFOROG7A.207    
     &,ZETA3        ! Even more work space                                 SFOROG7A.208    
                                                                           SFOROG7A.209    
!   Common parameters                                                      SFOROG7A.210    
                                                                           SFOROG7A.211    
*CALL C_MDI                                                                SFOROG7A.212    
*CALL C_VKMAN                                                              SFOROG7A.213    
*CALL C_SURF                                                               SFOROG7A.214    
                                                                           SFOROG7A.215    
!   Local parameters                                                       SFOROG7A.216    
                                                                           SFOROG7A.217    
      REAL H_BLEND_MIN,H_BLEND_MAX                                         SFOROG7A.218    
      PARAMETER (                                                          SFOROG7A.219    
     & H_BLEND_MIN=0.0       ! Minimun value of blending height            SFOROG7A.220    
     &,H_BLEND_MAX=1000.0    ! Maximum value of blending height            SFOROG7A.221    
     & )                                                                   SFOROG7A.222    
                                                                           SFOROG7A.223    
      EXTERNAL TIMER                                                       SFOROG7A.224    
                                                                           SFOROG7A.225    
      IF (LTIMER) THEN                                                     SFOROG7A.226    
        CALL TIMER('SF_OROG ',3)                                           SFOROG7A.227    
      ENDIF                                                                SFOROG7A.228    
                                                                           SFOROG7A.229    
! Set blending height, effective roughness length and                      SFOROG7A.230    
! wind profile factor at land points.                                      SFOROG7A.231    
                                                                           SFOROG7A.232    
      DO L = LAND1,LAND1+LAND_PTS-1                                        SFOROG7A.238    
        I = LAND_INDEX(L)                                                  SFOROG7A.239    
                                                                           SFOROG7A.241    
        H_BLEND_OROG(I) = H_BLEND_MIN                                      SFOROG7A.242    
        Z0M_EFF(I) = Z0M(L)                                                SFOROG7A.243    
                                                                           SFOROG7A.244    
        IF (L_Z0_OROG) THEN                                                SFOROG7A.245    
                                                                           SFOROG7A.246    
          ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L)                      SFOROG7A.247    
          ZETA2 = LOG ( 1.0 + Z1(I)/Z0M(L) )                               SFOROG7A.248    
          ZETA3 = 1.0 / SQRT ( ZETA1/(VKMAN*VKMAN) + 1.0/(ZETA2*ZETA2) )   SFOROG7A.249    
          ZETA2 = 1.0 / EXP(ZETA3)                                         SFOROG7A.250    
          H_BLEND_OROG(I) = MAX ( Z1(I) / (1.0 - ZETA2) ,                  ABX1F405.927    
     &                       HO2R2_OROG(L) * SQRT(2.0) )                   SFOROG7A.252    
          H_BLEND_OROG(I) = MIN ( H_BLEND_MAX, H_BLEND_OROG(I) )           ABX1F405.928    
                                                                           SFOROG7A.254    
! Apply simple stability correction to form drag if RIB is stable          SFOROG7A.255    
                                                                           SFOROG7A.256    
          IF (SIL_OROG(L) .EQ. RMDI) THEN                                  SFOROG7A.257    
             ZETA1 = 0.0                                                   SFOROG7A.258    
          ELSE                                                             SFOROG7A.259    
             RIB_FN =  ( 1.0 - RIB(I) / RI_CRIT )                          SFOROG7A.260    
             IF (RIB_FN.GT.1.0) RIB_FN = 1.0                               SFOROG7A.261    
             IF (RIB_FN.LT.0.0) RIB_FN = 0.0                               SFOROG7A.262    
             ZETA1 = 0.5 * OROG_DRAG_PARAM * SIL_OROG(L) * RIB_FN          SFOROG7A.263    
          ENDIF                                                            SFOROG7A.264    
                                                                           SFOROG7A.265    
          Z0M_EFF(I) = H_BLEND_OROG(I) / EXP ( VKMAN / SQRT ( ZETA1 +      ABX1F405.929    
     &                 (VKMAN / LOG ( H_BLEND_OROG(I) / Z0M(L) ) ) *       ABX1F405.930    
     &                 (VKMAN / LOG ( H_BLEND_OROG(I) / Z0M(L) ) ) ) )     ABX1F405.931    
          IF (RIB(I).GT.RI_CRIT) Z0M_EFF(I) = Z0M(L)                       SFOROG7A.269    
                                                                           SFOROG7A.270    
        ENDIF                                                              SFOROG7A.271    
                                                                           SFOROG7A.272    
      ENDDO                                                                SFOROG7A.276    
                                                                           SFOROG7A.277    
      IF (LTIMER) THEN                                                     SFOROG7A.278    
        CALL TIMER('SF_OROG ',4)                                           SFOROG7A.279    
      ENDIF                                                                SFOROG7A.280    
                                                                           SFOROG7A.281    
      RETURN                                                               SFOROG7A.282    
      END                                                                  SFOROG7A.283    
*ENDIF                                                                     SFOROG7A.284