*IF DEF,A03_7A                                                             RAERO7A.2      
C *****************************COPYRIGHT******************************     RAERO7A.3      
C (c) CROWN COPYRIGHT 1997, METEOROLOGICAL OFFICE, All Rights Reserved.    RAERO7A.4      
C                                                                          RAERO7A.5      
C Use, duplication or disclosure of this code is subject to the            RAERO7A.6      
C restrictions as set forth in the contract.                               RAERO7A.7      
C                                                                          RAERO7A.8      
C                Meteorological Office                                     RAERO7A.9      
C                London Road                                               RAERO7A.10     
C                BRACKNELL                                                 RAERO7A.11     
C                Berkshire UK                                              RAERO7A.12     
C                RG12 2SZ                                                  RAERO7A.13     
C                                                                          RAERO7A.14     
C If no contract has been raised with this copy of the code, the use,      RAERO7A.15     
C duplication or disclosure of it is strictly prohibited.  Permission      RAERO7A.16     
C to do so must first be obtained in writing from the Head of Numerical    RAERO7A.17     
C Modelling at the above address.                                          RAERO7A.18     
C ******************************COPYRIGHT******************************    RAERO7A.19     
!**********************************************************************    RAERO7A.20     
! Routine to calculate the aerodynamic resistance                          RAERO7A.21     
!                                                                          RAERO7A.22     
! Written by Peter Cox (June 1997)                                         RAERO7A.23     
!**********************************************************************    RAERO7A.24     

      SUBROUTINE RAERO (LAND_FIELD,LAND_INDEX,P_FIELD                       1RAERO7A.25     
     &,                 VEG_PTS,VEG_INDEX                                  RAERO7A.26     
     &,                 RIB,WIND,Z0H,Z0M,Z1,RA)                            RAERO7A.27     
                                                                           RAERO7A.28     
      IMPLICIT NONE                                                        RAERO7A.29     
                                                                           RAERO7A.30     
      INTEGER                                                              RAERO7A.31     
     & LAND_FIELD                 ! IN Total number of land points.        RAERO7A.32     
     &,LAND_INDEX(LAND_FIELD)     ! IN Index of land points on the         RAERO7A.33     
!                                 !    P-grid.                             RAERO7A.34     
     &,P_FIELD                    ! IN Total number of P-gridpoints.       RAERO7A.35     
     &,VEG_PTS                    ! IN Number of vegetated points.         RAERO7A.36     
     &,VEG_INDEX(LAND_FIELD)      ! IN Index of vegetated points           RAERO7A.37     
!                                 !    on the land grid.                   RAERO7A.38     
     &,I,J,L                      ! WORK Loop counters.                    RAERO7A.39     
                                                                           RAERO7A.40     
      REAL                                                                 RAERO7A.41     
     & RIB(P_FIELD)               ! IN Bulk Richardson number.             RAERO7A.42     
     &,WIND(P_FIELD)              ! IN Windspeed (m/s).                    RAERO7A.43     
     &,Z0H(LAND_FIELD)            ! IN Roughness length for heat (m).      RAERO7A.44     
     &,Z0M(LAND_FIELD)            ! IN Roughness length for momentum (m)   RAERO7A.45     
     &,Z1(P_FIELD)                ! IN Reference level (m).                RAERO7A.46     
     &,RA(LAND_FIELD)             ! OUT Aerodynamic resistance (s/m).      RAERO7A.47     
     &,BH                         ! WORK Stability coefficient.            RAERO7A.48     
     &,CHN(LAND_FIELD)            ! WORK Neutral drag coefficient.         RAERO7A.49     
     &,FH(LAND_FIELD)             ! WORK Stability factor.                 RAERO7A.50     
     &,ZETAH,ZETAM                ! WORK Tempories in calculation of       RAERO7A.51     
!                                 !      CHN.                              RAERO7A.52     
!-----------------------------------------------------------------------   RAERO7A.53     
! Local parameters                                                         RAERO7A.54     
!-----------------------------------------------------------------------   RAERO7A.55     
      REAL                                                                 RAERO7A.56     
     & AH,CZ                      ! Stability parameters.                  RAERO7A.57     
     &,KARMAN                     ! Von Karman's constant.                 RAERO7A.58     
      PARAMETER (AH = 10.0, CZ = 4.0, KARMAN = 0.4)                        RAERO7A.59     
                                                                           RAERO7A.60     
!-----------------------------------------------------------------------   RAERO7A.61     
! Calculate the neutral bulk tranfer coefficient.                          RAERO7A.62     
!-----------------------------------------------------------------------   RAERO7A.63     
      DO J=1,VEG_PTS                                                       RAERO7A.64     
        L = VEG_INDEX(J)                                                   RAERO7A.65     
        I = LAND_INDEX(L)                                                  RAERO7A.66     
        ZETAM = LOG((Z1(I) + Z0M(L)) / Z0M(L))                             RAERO7A.67     
        ZETAH = LOG((Z1(I) + Z0M(L)) / Z0H(L))                             RAERO7A.68     
        CHN(L) = (KARMAN * KARMAN) / (ZETAH * ZETAM)                       RAERO7A.69     
      ENDDO                                                                RAERO7A.70     
                                                                           RAERO7A.71     
!-----------------------------------------------------------------------   RAERO7A.72     
! Calculate the stability factor.                                          RAERO7A.73     
!-----------------------------------------------------------------------   RAERO7A.74     
      DO J=1,VEG_PTS                                                       RAERO7A.75     
        L = VEG_INDEX(J)                                                   RAERO7A.76     
        I = LAND_INDEX(L)                                                  RAERO7A.77     
        BH = AH * CHN(L) * CZ * SQRT (Z1(I) / Z0H(L))                      RAERO7A.78     
        IF (RIB(I) .GE. 0.0) THEN                                          RAERO7A.79     
          FH(L) = 1.0 / (1 + AH * RIB(I))                                  RAERO7A.80     
        ELSE                                                               RAERO7A.81     
          FH(L) = 1 - AH * RIB(I) / (1 + BH * SQRT(-RIB(I)))               RAERO7A.82     
        ENDIF                                                              RAERO7A.83     
      ENDDO                                                                RAERO7A.84     
                                                                           RAERO7A.85     
!-----------------------------------------------------------------------   RAERO7A.86     
! Calculate the aerodynamic resistance.                                    RAERO7A.87     
!-----------------------------------------------------------------------   RAERO7A.88     
      DO J=1,VEG_PTS                                                       RAERO7A.89     
        L = VEG_INDEX(J)                                                   RAERO7A.90     
        I = LAND_INDEX(L)                                                  RAERO7A.91     
        RA(L) = 1.0 / (FH(L) * CHN(L) * WIND(I))                           RAERO7A.92     
      ENDDO                                                                RAERO7A.93     
                                                                           RAERO7A.94     
      RETURN                                                               RAERO7A.95     
      END                                                                  RAERO7A.96     
*ENDIF                                                                     RAERO7A.97