*IF DEF,A03_3A                                                             ASJ4F401.8      
C ******************************COPYRIGHT******************************    GTS2F400.8803   
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.    GTS2F400.8804   
C                                                                          GTS2F400.8805   
C Use, duplication or disclosure of this code is subject to the            GTS2F400.8806   
C restrictions as set forth in the contract.                               GTS2F400.8807   
C                                                                          GTS2F400.8808   
C                Meteorological Office                                     GTS2F400.8809   
C                London Road                                               GTS2F400.8810   
C                BRACKNELL                                                 GTS2F400.8811   
C                Berkshire UK                                              GTS2F400.8812   
C                RG12 2SZ                                                  GTS2F400.8813   
C                                                                          GTS2F400.8814   
C If no contract has been raised with this copy of the code, the use,      GTS2F400.8815   
C duplication or disclosure of it is strictly prohibited.  Permission      GTS2F400.8816   
C to do so must first be obtained in writing from the Head of Numerical    GTS2F400.8817   
C Modelling at the above address.                                          GTS2F400.8818   
C ******************************COPYRIGHT******************************    GTS2F400.8819   
C                                                                          GTS2F400.8820   
C*LL  SUBROUTINE SF_EVAP------------------------------------------------   SFEVAP2B.3      
CLL                                                                        SFEVAP2B.4      
CLL  Purpose: Calculate surface evaporation and sublimation amounts        SFEVAP2B.5      
CLL           (without applying them to the surface stores).               SFEVAP2B.6      
CLL           Also calculate heat flux due to sea-ice melting.             SFEVAP2B.7      
CLL           Also calculate 1.5 metre T and Q.                            SFEVAP2B.8      
CLL                                                                        SFEVAP2B.9      
CLL         Change to the calculation of soil evaporation when             SFEVAP2B.10     
CLL    lying snow disappears during a timestep. Also two new               SFEVAP2B.11     
CLL    common decks C_HT_M and C_SICEHC                                    SFEVAP2B.12     
CLL                                                                        SFEVAP2B.13     
CLL  Suitable for single column usage.                                     SFEVAP2B.14     
CLL                                                                        SFEVAP2B.15     
CLL  Model            Modification history:                                SFEVAP2B.16     
CLL version  Date                                                          SFEVAP2B.17     
CLL   3.1   11/1/93   Co-ordinate transformation so that screen            SFEVAP2B.18     
CLL                   temperature is calculated at 1.5m above              SFEVAP2B.19     
CLL                   surface + 1.5m            S.Jackson                  SFEVAP2B.20     
CLL   3.4   06/06/94  DEF TIMER replaced by LOGICAL LTIMER                 ASJ1F304.382    
CLL                   Argument LTIMER added                                ASJ1F304.383    
CLL                                                 S.J.Swarbrick          ASJ1F304.384    
CLL   4.1   08/05/96  decks A03_2C and A03_3B removed                      ASJ4F401.9      
CLL                                     S D Jackson                        ASJ4F401.10     
CLL   4.2   Oct. 96   T3E migration - *DEF CRAY removed                    GSS1F402.68     
CLL                                     S J Swarbrick                      GSS1F402.69     
CLL  4.5    Jul. 98  Kill the IBM specific lines (JCThil)                  AJC1F405.127    
CLL                                                                        SFEVAP2B.21     
CLL  Programming standard: Unified Model Documentation Paper No 4,         SFEVAP2B.22     
CLL                        version 2, dated 18/1/90.                       SFEVAP2B.23     
CLL                                                                        SFEVAP2B.24     
CLL  Logical component covered: P245.                                      SFEVAP2B.25     
CLL                                                                        SFEVAP2B.26     
CLL  System task:                                                          SFEVAP2B.27     
CLL                                                                        SFEVAP2B.28     
CLL  Documentation: UMDP 24                                                SFEVAP2B.29     
CLL                                                                        SFEVAP2B.30     
CLL---------------------------------------------------------------------   SFEVAP2B.31     
C*                                                                         SFEVAP2B.32     
C*L Arguments :---------------------------------------------------------   SFEVAP2B.33     

      SUBROUTINE SF_EVAP (                                                  4,14SFEVAP2B.34     
     + P_FIELD,P1,LAND_FIELD,LAND1,                                        SFEVAP2B.35     
     + POINTS,BL_LEVELS,LAND_MASK,LAND_PTS,LAND_INDEX                      SFEVAP2B.39     
     +,ASOIL_1,CANOPY,CATCH,DTRDZ,DTRDZ_RML,EA,ESL,SOIL_HT_FLUX            SFEVAP2B.41     
     +,E_SEA,ICE_FRACT,LYING_SNOW,RADNET,SMC,TIMESTEP                      SFEVAP2B.42     
     +,CER1P5M,CHR1P5M,PSTAR,Z1,Z0M,Z0H,SQ1P5,ST1P5,SMLT                   SFEVAP2B.43     
     +,NRML,RHOKH_1,DQW_1,DQW_RML,DTSTAR,FTL,FQW,ES,QW,TL,TSTAR            SFEVAP2B.44     
     +,ECAN,EI                                                             SFEVAP2B.45     
     +,SICE_MLT_HTF,Q1P5M,T1P5M,LTIMER                                     ASJ1F304.385    
     +)                                                                    SFEVAP2B.47     
      IMPLICIT NONE                                                        SFEVAP2B.48     
      LOGICAL LTIMER                                                       ASJ1F304.386    
      INTEGER                                                              SFEVAP2B.49     
     + P_FIELD              ! IN No. of gridpoints in the whole grid.      SFEVAP2B.50     
     +,P1                   ! IN 1st P-pt in full field to be processed.   SFEVAP2B.51     
     +,LAND_FIELD           ! IN No. of LANDpoints in the whole grid.      SFEVAP2B.52     
     +,LAND1                ! IN 1st L-pt in full field to be processed.   SFEVAP2B.53     
     +,POINTS               ! IN No. of gridpoints to be processed.        SFEVAP2B.54     
     +,BL_LEVELS            ! IN No. of levels treated by b.l. scheme.     SFEVAP2B.55     
     +,LAND_PTS             ! IN No. of land points to be processed.       SFEVAP2B.56     
      LOGICAL                                                              SFEVAP2B.57     
     + LAND_MASK(P_FIELD)   ! IN T for land points, F otherwise.           SFEVAP2B.58     
      INTEGER                                                              SFEVAP2B.60     
     + LAND_INDEX(P_FIELD)  ! IN Index of land points on the P-grid.       SFEVAP2B.61     
C                           !    The ith element contains the position     SFEVAP2B.62     
C                           !    in whole grid of the ith land point.      SFEVAP2B.63     
      REAL                                                                 SFEVAP2B.65     
     + ASOIL_1(P_FIELD)     ! IN Soil coefficient from P242 (sq m K per    SFEVAP2B.66     
C                           !    per Joule * timestep).                    SFEVAP2B.67     
     +,CANOPY(LAND_FIELD)   ! IN Gridbox mean canopy / surface water       SFEVAP2B.68     
C                           !    store (kg per sq m).                      SFEVAP2B.69     
     +,CATCH(LAND_FIELD)    ! IN Canopy / surface water store capacity     SFEVAP2B.70     
C                           !    (kg per sq m).                            SFEVAP2B.71     
     +,DTRDZ(P_FIELD,       ! IN -g.dt/dp for each model layer on p-grid   SFEVAP2B.72     
     +       BL_LEVELS)     !    From P244 ((kg/sq m/s)**-1).              SFEVAP2B.73     
     +,DTRDZ_RML(P_FIELD)   ! IN -g.dt/dp for the rapidly mixing layer     SFEVAP2B.74     
C                           !    (if it exists) on the p-grid from P244.   SFEVAP2B.75     
     +,EA(P_FIELD)          ! IN Surface evaporation with only aero-       SFEVAP2B.76     
C                           !    dynamic resistance (+ve), or condens-     SFEVAP2B.77     
C                           !    ation (-ve), averaged over gridbox.       SFEVAP2B.78     
C                           !    From P243 and P244.  Kg per sq m per s.   SFEVAP2B.79     
     +,ESL(P_FIELD)         ! IN ES (q.v.) without fractional weighting    SFEVAP2B.80     
C                           !    factor FRACS ('L' is for 'local').        SFEVAP2B.81     
C                           !    From P243 and P244.  Kg per sq m per s.   SFEVAP2B.82     
     +,SOIL_HT_FLUX(P_FIELD) ! IN Heat flux from surface to deep soil      SFEVAP2B.83     
C                           !    layer 1 (Watts per sq m).  From P242.     SFEVAP2B.84     
     +,E_SEA(P_FIELD)       ! IN Evaporation from sea (weighted with       SFEVAP2B.85     
C                           !    leads fraction at sea-ice points).        SFEVAP2B.86     
C                           !       Diagnostics defined on land and sea.   SFEVAP2B.88     
     +,ICE_FRACT(P_FIELD)   ! IN Fraction of gridbox which is covered by   SFEVAP2B.90     
C                           !    sea-ice (decimal fraction, but most of    SFEVAP2B.91     
C                           !    this sub-component assumes it to be       SFEVAP2B.92     
C                           !    either 1.0 or 0.0 precisely).             SFEVAP2B.93     
     +,LYING_SNOW(P_FIELD)  ! IN Lying snow (kg per sq m).                 SFEVAP2B.94     
C                           !    NB Dimension is PFIELD not LAND_FIELD f   SFEVAP2B.96     
C                           !    snow on sea-ice in coupled model runs.    SFEVAP2B.97     
     +,RADNET(P_FIELD)      ! IN Surface net radiation (Watts per sq m).   SFEVAP2B.99     
     +,SMC(LAND_FIELD)      ! IN Soil moisture content (kg per sq m).      SFEVAP2B.100    
     +,TIMESTEP             ! IN Timestep (sec).                           SFEVAP2B.101    
      LOGICAL                                                              SFEVAP2B.102    
     + SQ1P5                ! IN STASH flag for 1.5-metre sp humidity.     SFEVAP2B.103    
     +,ST1P5                ! IN STASH flag for 1.5-metre temperature.     SFEVAP2B.104    
     +,SMLT                 ! IN STASH flag for sea-ice melting ht flux.   SFEVAP2B.105    
      REAL                                                                 SFEVAP2B.106    
     + CER1P5M(P_FIELD)     ! IN Transfer coefficient ratio, from P243.    SFEVAP2B.107    
     +,CHR1P5M(P_FIELD)     ! IN Transfer coefficient ratio, from P243.    SFEVAP2B.108    
     +,PSTAR(P_FIELD)       ! IN Surface pressure (Pa).                    SFEVAP2B.109    
     +,Z1(P_FIELD)          ! IN Height of lowest atmospheric level        SFEVAP2B.110    
C                           !    (i.e. middle of lowest layer).  Metres.   SFEVAP2B.111    
     +,Z0M(P_FIELD)         ! IN Roughness length for momentum (m)         SFEVAP2B.112    
     +,Z0H(P_FIELD)         ! IN Roughness length for heat and moisture    SFEVAP2B.113    
      INTEGER                                                              SFEVAP2B.114    
     & NRML(P_FIELD)        ! IN  The Number of model layers in the        SFEVAP2B.115    
C                           !     Rapidly Mixing Layer.                    SFEVAP2B.116    
      REAL                                                                 SFEVAP2B.117    
     + RHOKH_1(P_FIELD)     ! IN    Turbulent surface exchange             SFEVAP2B.118    
C                           !       coefficient for sensible heat.         SFEVAP2B.119    
     +,DQW_1(P_FIELD)       ! IN    Increment for lowest-level total       SFEVAP2B.120    
C                           !       water content.  From P244.             SFEVAP2B.121    
     +,DQW_RML(P_FIELD)     ! IN    Increment to QW due to rapid mixing.   SFEVAP2B.122    
     +,DTSTAR(P_FIELD)      ! IN    Increment for surface temperature.     SFEVAP2B.123    
C                           !       From P244.                             SFEVAP2B.124    
     +,FTL(P_FIELD,         ! INOUT Sensible heat flux from layer k-1 to   SFEVAP2B.125    
     +     BL_LEVELS)       !       layer k (W/sq m).  From P243 and       SFEVAP2B.126    
C                           !      P244, units changed in P24 top level.   SFEVAP2B.127    
     +,FQW(P_FIELD,         ! INOUT Turbulent moisture flux from level     SFEVAP2B.128    
     +     BL_LEVELS)       !       k-1 to k (kg/sq m/s). From P243/4.     SFEVAP2B.129    
     +,ES(P_FIELD)          ! INOUT Surface evapotranspiration (through    SFEVAP2B.130    
C                           !       a resistance which is not entirely     SFEVAP2B.131    
C                           !       aerodynamic).  Always non-negative.    SFEVAP2B.132    
C                           !       Kg per sq m per sec.  From P243/4.     SFEVAP2B.133    
C                           !       Diagnostics defined on land and sea.   SFEVAP2B.135    
     +,QW(P_FIELD,BL_LEVELS)! INOUT Total water content (kg(water)/        SFEVAP2B.137    
C                           !       kg(air)).  From P243/4.                SFEVAP2B.138    
     +,TL(P_FIELD,BL_LEVELS)! INOUT Liquid/frozen water temperature (K).   SFEVAP2B.139    
     +,TSTAR(P_FIELD)       ! INOUT Surface temperature (K).               SFEVAP2B.140    
      REAL                                                                 SFEVAP2B.141    
     + ECAN(P_FIELD)        ! OUT Gridbox mean evaporation from canopy/    SFEVAP2B.142    
C                           !     surface store (kg per sq m per s).       SFEVAP2B.143    
C                           !     Zero over sea and sea-ice.               SFEVAP2B.144    
C                           !     Diagnostics defined on land and sea.     SFEVAP2B.146    
     +,EI(P_FIELD)          ! OUT Sublimation from lying snow or sea-      SFEVAP2B.148    
C                           !     ice (kg per sq m per s).                 SFEVAP2B.149    
      REAL                                                                 SFEVAP2B.150    
     + SICE_MLT_HTF(P_FIELD)! OUT Heat flux due to melting of sea-ice      SFEVAP2B.151    
C                           !     (Watts per square metre).                SFEVAP2B.152    
     +,Q1P5M(P_FIELD)       ! OUT Specific humidity at screen height of    SFEVAP2B.153    
C                           !     1.5 metres (kg water per kg air).        SFEVAP2B.154    
     +,T1P5M(P_FIELD)       ! OUT Temperature at 1.5 metres above the      SFEVAP2B.155    
C                           !     surface (K).                             SFEVAP2B.156    
C*                                                                         SFEVAP2B.157    
C*L  External subprogram(s) required :-                                    SFEVAP2B.158    
      EXTERNAL QSAT                                                        SFEVAP2B.159    
      EXTERNAL TIMER                                                       SFEVAP2B.161    
C*                                                                         SFEVAP2B.163    
C*L  Local and other symbolic constants used :-                            SFEVAP2B.164    
*CALL C_0_DG_C                                                             SFEVAP2B.165    
*CALL C_LHEAT                                                              SFEVAP2B.166    
*CALL C_G                                                                  SFEVAP2B.167    
*CALL C_HT_M                   ! Contains Z1P5M                            SFEVAP2B.168    
*CALL C_R_CP                                                               SFEVAP2B.169    
*CALL C_SICEHC                 ! Contains AI                               SFEVAP2B.170    
      REAL GRCP                                                            SFEVAP2B.171    
      PARAMETER (                                                          SFEVAP2B.172    
     + GRCP=G/CP   ! Accn due to gravity / standard heat capacity of       SFEVAP2B.173    
C                  ! air at const pressure.  Used in diagnosis of 1.5      SFEVAP2B.174    
C                  ! metre temperature.                                    SFEVAP2B.175    
     +)                                                                    SFEVAP2B.176    
C*                                                                         SFEVAP2B.177    
      REAL                                                                 SFEVAP2B.179    
     + QS(P_FIELD)           ! Used for saturated specific humidity        SFEVAP2B.180    
C                            ! at surface, in Q1P5M calculation.           SFEVAP2B.181    
      REAL                                                                 SFEVAP2B.182    
     + DTL_RML(P_FIELD)      ! Adjustment increment to TL in the           SFEVAP2B.183    
C                            ! rapidly mixing layer.                       SFEVAP2B.184    
     +,DFTL(P_FIELD,         ! Adjustment increment to the sensible        SFEVAP2B.185    
     +      BL_LEVELS)       ! heat flux.                                  SFEVAP2B.186    
     +,DFQW(P_FIELD,         ! Adjustment increment to the flux of         SFEVAP2B.187    
     +      BL_LEVELS)       ! total water.                                SFEVAP2B.188    
C  Local scalars                                                           SFEVAP2B.202    
      REAL                                                                 SFEVAP2B.203    
     + EADT                ! EA (q.v.) integrated over timestep.           SFEVAP2B.204    
     +,ECANDT              ! ECAN (q.v.) integrated over timestep.         SFEVAP2B.205    
     +,EDT                 ! E=FQW(,1) (q.v.) integrated over timestep.    SFEVAP2B.206    
     +,EIDT                ! EI (q.v.) integrated over timestep.           SFEVAP2B.207    
     +,ESDT                ! ES (q.v.) integrated over timestep.           SFEVAP2B.208    
     +,ESLDT               ! ESL (q.v.) integrated over timestep.          SFEVAP2B.209    
     +,EW                  ! Total surface flux of water, excluding        SFEVAP2B.210    
C                          ! sublimation/frost deposition, over land.      SFEVAP2B.211    
     +,FRACA               ! Fraction of gridbox at which moisture flux    SFEVAP2B.212    
C                          ! is going at the potential rate, i.e. with     SFEVAP2B.213    
C                          ! only aerodynamic resistance.                  SFEVAP2B.214    
     +,FRACS               ! Fraction of gridbox at which moisture flux    SFEVAP2B.215    
C                          ! is additionally impeded by a surface and/or   SFEVAP2B.216    
C                          ! stomatal resistance.                          SFEVAP2B.217    
     +,CT_1                ! Parameter in the implicit adjustment of       SFEVAP2B.218    
C                          ! TSTAR and TL(,1).                             SFEVAP2B.219    
     +,CT_RML              ! Parameter in the implicit adjustment of       SFEVAP2B.220    
C                          ! TSTAR and TL in the rapidly mixing layer.     SFEVAP2B.221    
     +,TSTARMAX            ! Maximum gridbox mean surface temperature in   SFEVAP2B.222    
C                          ! sea-ice melting calculation.  See comments    SFEVAP2B.223    
C                          ! to section 1.6.3(b) or 2.5.3(b) below.        SFEVAP2B.224    
      INTEGER                                                              SFEVAP2B.225    
     + I                   ! Loop counter - full horizontal field index.   SFEVAP2B.226    
     +,L                   ! Loop counter - land field index.              SFEVAP2B.227    
     +,K                   ! Loop counter in the vertical.                 SFEVAP2B.228    
     +,KM1                 ! K - 1                                         SFEVAP2B.229    
     +,NRMLP1              ! NRML + 1                                      SFEVAP2B.230    
C                                                                          SFEVAP2B.231    
C-----------------------------------------------------------------------   SFEVAP2B.232    
CL 1. Initialise some output variables to zero.                            SFEVAP2B.233    
C-----------------------------------------------------------------------   SFEVAP2B.234    
C                                                                          SFEVAP2B.235    
      IF (LTIMER) THEN                                                     ASJ1F304.387    
      CALL TIMER('SFEVAP  ',3)                                             SFEVAP2B.237    
      ENDIF                                                                ASJ1F304.388    
      DO 1 I=P1,P1+POINTS-1                                                SFEVAP2B.239    
        ECAN(I) = 0.0                                                      SFEVAP2B.240    
        EI(I) = 0.0                                                        SFEVAP2B.241    
        IF (SMLT)                                                          SFEVAP2B.242    
     +    SICE_MLT_HTF(I) = 0.0                                            SFEVAP2B.243    
    1 CONTINUE                                                             SFEVAP2B.244    
C                                                                          SFEVAP2B.245    
C---------------------------------------------------------------------     SFEVAP2B.246    
CL 2. Do calculations for land points.                                     SFEVAP2B.247    
C---------------------------------------------------------------------     SFEVAP2B.248    
C                                                                          SFEVAP2B.249    
CMIC$ DO ALL VECTOR SHARED(P_FIELD, LAND_FIELD, BL_LEVELS, LAND1,          SFEVAP2B.250    
CMIC$1   LAND_PTS, LAND_INDEX, ESL, TIMESTEP, ES, LYING_SNOW, ECAN,        SFEVAP2B.251    
CMIC$2   EA, CATCH, CANOPY, SMC, EI, TSTAR, FQW) PRIVATE(I, L, ESLDT,      SFEVAP2B.252    
CMIC$3   ESDT, EADT, EDT, ECANDT, FRACA, FRACS, EW, EIDT)                  SFEVAP2B.253    
CDIR$ IVDEP                                                                SFEVAP2B.254    
! Fujitsu vectorization directive                                          GRB0F405.443    
!OCL NOVREC                                                                GRB0F405.444    
        DO 2 L=LAND1,LAND1+LAND_PTS-1                                      SFEVAP2B.260    
          I = LAND_INDEX(L)                                                SFEVAP2B.261    
C                                                                          SFEVAP2B.263    
C-----------------------------------------------------------------------   SFEVAP2B.264    
CL 2.1 Calculate fluxes integrated over timestep.                          SFEVAP2B.265    
C-----------------------------------------------------------------------   SFEVAP2B.266    
C                                                                          SFEVAP2B.267    
            ESLDT = ESL(I) * TIMESTEP                                      SFEVAP2B.268    
            EADT = EA(I) * TIMESTEP                                        SFEVAP2B.269    
            ESDT = ES(I) * TIMESTEP                                        SFEVAP2B.270    
            EDT = EADT + ESDT                                              SFEVAP2B.271    
C                                                                          SFEVAP2B.272    
C-----------------------------------------------------------------------   SFEVAP2B.273    
CL 2.2 Do calculations for snow-free land.  Canopy processes operate.      SFEVAP2B.274    
CL     LYING_SNOW is defined on sea and land points for snow on sea-ice    SFEVAP2B.276    
CL     in coupled model runs.                                              SFEVAP2B.277    
C-----------------------------------------------------------------------   SFEVAP2B.279    
C                                                                          SFEVAP2B.280    
            IF (LYING_SNOW(I).LE.0.0) THEN                                 SFEVAP2B.281    
              IF (EDT.GE.0.0) THEN                                         SFEVAP2B.282    
C                                                                          SFEVAP2B.283    
C-----------------------------------------------------------------------   SFEVAP2B.284    
CL 2.2.1 Non-negative moisture flux over snow-free land.                   SFEVAP2B.285    
C-----------------------------------------------------------------------   SFEVAP2B.286    
C                                                                          SFEVAP2B.287    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.288    
CL   (a) Water in canopy and soil is assumed to be liquid, so all          SFEVAP2B.289    
CL       positive moisture flux over snow-free land is evaporation         SFEVAP2B.290    
CL       rather than sublimation, even if TSTAR is less than or equal      SFEVAP2B.291    
CL       to TM.                                                            SFEVAP2B.292    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.293    
C                                                                          SFEVAP2B.294    
                ECAN(I) = EA(I)                                            SFEVAP2B.295    
                ECANDT = EADT                                              SFEVAP2B.296    
C                                                                          SFEVAP2B.297    
C  If EDT is non-negative, then ECANDT must be non-negative.               SFEVAP2B.298    
C                                                                          SFEVAP2B.299    
                FRACA = 0.0                                                SFEVAP2B.300    
                IF (CATCH(L).GT.0.0)                                       SFEVAP2B.301    
     +            FRACA = CANOPY(L) / CATCH(L)                             SFEVAP2B.302    
                IF (CANOPY(L).LT.ECANDT) THEN                              SFEVAP2B.303    
C                                                                          SFEVAP2B.304    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.305    
CL   (b) It is assumed that any 'canopy' moisture flux in excess of the    SFEVAP2B.306    
CL       current canopy water amount is in fact soil evaporation.          SFEVAP2B.307    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.308    
C                                                                          SFEVAP2B.309    
C        This situation is highly improbable - it will occur at, at        SFEVAP2B.310    
C        most, a few gridpoints in any given timestep.                     SFEVAP2B.311    
C                                                                          SFEVAP2B.312    
                  FRACS = 1.0 - FRACA*( CANOPY(L) / ECANDT )               SFEVAP2B.313    
                  ESDT = ESLDT * FRACS                                     SFEVAP2B.314    
                  ECANDT = CANOPY(L)                                       SFEVAP2B.315    
                  ECAN(I) = ECANDT / TIMESTEP                              SFEVAP2B.316    
                  ES(I) = ESDT / TIMESTEP                                  SFEVAP2B.317    
                ENDIF                                                      SFEVAP2B.318    
C                                                                          SFEVAP2B.319    
C  (The canopy store is depleted by evaporation in P252, and not here,     SFEVAP2B.320    
C   according to the formula: CANOPY=CANOPY-ECANDT)                        SFEVAP2B.321    
C                                                                          SFEVAP2B.322    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.323    
CL   (c) Adjustments to evaporation from soil as calculated so far :-      SFEVAP2B.324    
C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -   SFEVAP2B.325    
C                                                                          SFEVAP2B.326    
              IF (SMC(L).LE.0.0) THEN                                      SFEVAP2B.327    
C                                                                          SFEVAP2B.328    
CL   (i) If there is currently no soil moisture, there must be no          SFEVAP2B.329    
CL       evaporation of soil moisture, so this flux is set to zero.        SFEVAP2B.330    
C                                                                          SFEVAP2B.331    
                  ESDT = 0.0                                               SFEVAP2B.332    
                  ES(I) = 0.0                                              SFEVAP2B.333    
                ELSEIF (SMC(L).LT.ESDT) THEN                               SFEVAP2B.334    
C                                                                          SFEVAP2B.335    
CL  (ii) Ensure that the soil evaporation is not greater than the          SFEVAP2B.336    
CL       current soil moisture store.                                      SFEVAP2B.337    
C        This situation is extremely unlikely at any given gridpoint       SFEVAP2B.338    
C        at any given timestep.                                            SFEVAP2B.339    
C                                                                          SFEVAP2B.340    
                  ESDT = SMC(L)                                            SFEVAP2B.341    
                  ES(I) = ESDT / TIMESTEP                                  SFEVAP2B.342    
                ENDIF                                                      SFEVAP2B.343    
C                                                                          SFEVAP2B.344    
C  (The soil moisture store is depleted by evaporation in P253, and not    SFEVAP2B.345    
C   here, using the formula:  SMC=SMC-ESDT)                                SFEVAP2B.346    
C                                                                          SFEVAP2B.347    
                EW = ECAN(I) + ES(I)                                       SFEVAP2B.348    
                EI(I) = 0.0                                                SFEVAP2B.349    
C                                                                          SFEVAP2B.350    
C-----------------------------------------------------------------------   SFEVAP2B.351    
CL 2.2.2 Negative moisture flux onto snow-free land above freezing         SFEVAP2B.352    
C-----------------------------------------------------------------------   SFEVAP2B.353    
CL       (i.e. condensation onto snow-free land).  The whole flux is       SFEVAP2B.354    
CL       into the surface/canopy store.                                    SFEVAP2B.355    
C                                                                          SFEVAP2B.356    
              ELSEIF (TSTAR(I).GT.TM) THEN     ! ELSE of evaporation /     SFEVAP2B.357    
C                                              ! condensation block.       SFEVAP2B.358    
C                                                                          SFEVAP2B.359    
C  Condensation implies ES=0, so ECAN=EA=EW=E (=FQW(,1))                   SFEVAP2B.360    
C                                                                          SFEVAP2B.361    
                ECAN(I) = FQW(I,1)                                         SFEVAP2B.362    
                ES(I) = 0.0                                                SFEVAP2B.363    
                EW = ECAN(I)                                               SFEVAP2B.364    
                EI(I) = 0.0                                                SFEVAP2B.365    
C                                                                          SFEVAP2B.366    
C  (The canopy store is augmented by interception of condensation at       SFEVAP2B.367    
C   P252, and not here.)                                                   SFEVAP2B.368    
C                                                                          SFEVAP2B.369    
C-----------------------------------------------------------------------   SFEVAP2B.370    
CL 2.2.3 Negative moisture flux onto snow-free land below freezing         SFEVAP2B.371    
CL       (i.e. deposition of frost).                                       SFEVAP2B.372    
C-----------------------------------------------------------------------   SFEVAP2B.373    
C                                                                          SFEVAP2B.374    
              ELSE      ! ELSE of condensation / frost deposition block.   SFEVAP2B.375    
                EI(I) = FQW(I,1)                                           SFEVAP2B.376    
                ES(I) = 0.0                                                SFEVAP2B.377    
                EW = 0.0                                                   SFEVAP2B.378    
C                                                                          SFEVAP2B.379    
C  (Negative EI is used to increment the snowdepth store - there is        SFEVAP2B.380    
C   no separate "frost" store.  This incrementing is done in P251,         SFEVAP2B.381    
C   according to:  LYING_SNOW = LYING_SNOW - EI*TIMESTEP)                  SFEVAP2B.382    
C                                                                          SFEVAP2B.383    
              ENDIF  ! End of evaporation/condensation/deposition block.   SFEVAP2B.384    
C                                                                          SFEVAP2B.385    
C-----------------------------------------------------------------------   SFEVAP2B.386    
CL 2.3 Do calculations for snow-covered land.                              SFEVAP2B.387    
C-----------------------------------------------------------------------   SFEVAP2B.388    
C                                                                          SFEVAP2B.389    
            ELSEIF (LYING_SNOW(I).LE.EDT) THEN     ! ELSEIF of no-snow.    SFEVAP2B.390    
C                                                                          SFEVAP2B.391    
C-----------------------------------------------------------------------   SFEVAP2B.392    
CL 2.3.1 Shallow snow (lying snow or frost which is being exhausted        SFEVAP2B.393    
CL       by evaporation).  All the snow is sublimated, the remaining       SFEVAP2B.394    
CL       moisture flux being taken from the canopy and soil, with all      SFEVAP2B.395    
CL       the palaver of section 1.2.1 above.                               SFEVAP2B.396    
C-----------------------------------------------------------------------   SFEVAP2B.397    
C                                                                          SFEVAP2B.398    
C        This is extremely unlikely at more than one or two gridpoints     SFEVAP2B.399    
C        at any given timestep, yet the complicated logic probably         SFEVAP2B.400    
C        slows down the routine considerably - this section is a           SFEVAP2B.401    
C        suitable candidate for further consideration as regards           SFEVAP2B.402    
C        making the model optimally efficient.                             SFEVAP2B.403    
C                                                                          SFEVAP2B.404    
              EI(I) = LYING_SNOW(I) / TIMESTEP                             SFEVAP2B.405    
              EIDT = LYING_SNOW(I)                                         SFEVAP2B.406    
C                                                                          SFEVAP2B.407    
C  Set EDT = ( E - SNOSUB ) * TIMESTEP.  This is the moisture in kg per    SFEVAP2B.408    
C  square metre left over to be evaporated from the canopy and soil.       SFEVAP2B.409    
C  N.B.  E=FQW(,1)                                                         SFEVAP2B.410    
C                                                                          SFEVAP2B.411    
              EDT = EDT - EIDT                                             SFEVAP2B.412    
C                                                                          SFEVAP2B.413    
C  (Snowdepth is decreased using EI at P251, and not here.  The formula    SFEVAP2B.414    
C   used is simply:  LYING_SNOW = LYING_SNOW - EI*TIMESTEP.)               SFEVAP2B.415    
C                                                                          SFEVAP2B.416    
C  Now that all the snow has sublimed, canopy processes come into          SFEVAP2B.417    
C  operation (FRACA no longer necessarily equal to 1).                     SFEVAP2B.418    
C                                                                          SFEVAP2B.419    
              FRACA = 0.0                                                  SFEVAP2B.420    
              IF (CATCH(L).GT.0.0)                                         SFEVAP2B.421    
     +          FRACA = CANOPY(L) / CATCH(L)                               SFEVAP2B.422    
              ECANDT = EDT * FRACA                                         SFEVAP2B.423    
              IF (CANOPY(L).LT.ECANDT) THEN                                SFEVAP2B.424    
C                                                                          SFEVAP2B.425    
C  Dry out the canopy completely and assume the remaining moisture flux    SFEVAP2B.426    
C  is soil evaporation.                                                    SFEVAP2B.427    
C                                                                          SFEVAP2B.428    
                FRACS = 1.0 - FRACA*( CANOPY(L) / ECANDT )                 SFEVAP2B.429    
                ESDT = EDT * FRACS                                         SFEVAP2B.430    
                ECANDT = CANOPY(L)                                         SFEVAP2B.431    
              ELSE                                                         SFEVAP2B.432    
C                                                                          SFEVAP2B.433    
C  Calculate soil evaporation.                                             SFEVAP2B.434    
C                                                                          SFEVAP2B.435    
                FRACS = 1.0 - FRACA                                        SFEVAP2B.436    
                ESDT = EDT * FRACS                                         SFEVAP2B.437    
              ENDIF                                                        SFEVAP2B.438    
              ECAN(I) = ECANDT / TIMESTEP                                  SFEVAP2B.439    
              ES(I) = ESDT / TIMESTEP                                      SFEVAP2B.440    
C                                                                          SFEVAP2B.441    
C  (ECAN is used to deplete the canopy store at P252, and not here.  The   SFEVAP2B.442    
C   formula used is simply:  CANOPY = CANOPY - ECAN*TIMESTEP.)             SFEVAP2B.443    
C                                                                          SFEVAP2B.444    
C  Evaporation from soil.                                                  SFEVAP2B.445    
C                                                                          SFEVAP2B.446    
              IF (SMC(L).LE.0.0) THEN                                      SFEVAP2B.447    
C                                                                          SFEVAP2B.448    
C  No evaporation from soil possible when there is no soil moisture.       SFEVAP2B.449    
C                                                                          SFEVAP2B.450    
                ESDT = 0.0                                                 SFEVAP2B.451    
                ES(I) = 0.0                                                SFEVAP2B.452    
              ELSEIF (SMC(L).LT.ESDT) THEN                                 SFEVAP2B.453    
C                                                                          SFEVAP2B.454    
C  Limit evaporation of soil moisture in the extremely unlikely event      SFEVAP2B.455    
C  that soil moisture is exhausted by the evaporation left over from       SFEVAP2B.456    
C  sublimation which exhausted the snow store.                             SFEVAP2B.457    
C                                                                          SFEVAP2B.458    
                ESDT = SMC(L)                                              SFEVAP2B.459    
                ES(I) = ESDT / TIMESTEP                                    SFEVAP2B.460    
              ENDIF                                                        SFEVAP2B.461    
C                                                                          SFEVAP2B.462    
C  (ES is used to deplete the soil moisture store at P253, and not here,   SFEVAP2B.463    
C   according to the formula:  SMC = SMC - ES*TIMESTEP.)                   SFEVAP2B.464    
C                                                                          SFEVAP2B.465    
              EW = ECAN(I) + ES(I)                                         SFEVAP2B.466    
C                                                                          SFEVAP2B.467    
C-----------------------------------------------------------------------   SFEVAP2B.468    
CL 2.3.2 Deep snow (i.e. not being exhausted by evaporation).  This        SFEVAP2B.469    
CL       covers two cases: (a) sublimation from deep snow (if total        SFEVAP2B.470    
CL       moisture flux over the timestep is non-negative but less than     SFEVAP2B.471    
CL       the lying snow amount), and (b) deposition onto an already        SFEVAP2B.472    
CL       snowy surface (if the total moisture flux is negative and         SFEVAP2B.473    
CL       the lying snow amount is positive).                               SFEVAP2B.474    
C-----------------------------------------------------------------------   SFEVAP2B.475    
C                                                                          SFEVAP2B.476    
            ELSE          ! ELSE of shallow snow / deep snow block.        SFEVAP2B.477    
              EI(I) = FQW(I,1)                                             SFEVAP2B.478    
              EW = 0.0                                                     SFEVAP2B.479    
C                                                                          SFEVAP2B.480    
C  (EI is used to increase or decrease the snowdepth at P251, and not      SFEVAP2B.481    
C   here, according to the formula:                                        SFEVAP2B.482    
C            LYING_SNOW = LYING_SNOW - EI*TIMESTEP . )                     SFEVAP2B.483    
C                                                                          SFEVAP2B.484    
            ENDIF         ! End of no snow/shallow snow/deep snow block.   SFEVAP2B.485    
            FQW(I,1) = EW + EI(I)                                          SFEVAP2B.486    
    2   CONTINUE                                                           SFEVAP2B.488    
C                                                                          SFEVAP2B.489    
C  Split loop 2 here so that it will vectorise.                            SFEVAP2B.490    
C                                                                          SFEVAP2B.491    
CDIR$ IVDEP                                                                SFEVAP2B.492    
! Fujitsu vectorization directive                                          GRB0F405.445    
!OCL NOVREC                                                                GRB0F405.446    
        DO 24 L=LAND1,LAND1+LAND_PTS-1                                     SFEVAP2B.493    
          I = LAND_INDEX(L)                                                SFEVAP2B.494    
C                                                                          SFEVAP2B.496    
C-----------------------------------------------------------------------   SFEVAP2B.497    
CL 2.4 Do calculations which are the same at all land points.  Calculate   SFEVAP2B.498    
CL     and apply increments to (a) total atmospheric water content QW,     SFEVAP2B.499    
CL     (b) surface temperature TSTAR, (c) TL                               SFEVAP2B.500    
C-----------------------------------------------------------------------   SFEVAP2B.501    
C                                                                          SFEVAP2B.502    
            DTSTAR(I) = ASOIL_1(I) * ( RADNET(I) - FTL(I,1)                SFEVAP2B.503    
     &                  - LC*FQW(I,1) - LF*EI(I) - SOIL_HT_FLUX(I) )       SFEVAP2B.504    
     &                  - DTSTAR(I)                                        SFEVAP2B.505    
            IF ( NRML(I) .GE. 2 ) THEN                                     SFEVAP2B.506    
              NRMLP1 = NRML(I) + 1                                         SFEVAP2B.507    
              DQW_RML(I) = -DTRDZ_RML(I) * ( FQW(I,NRMLP1) - FQW(I,1) )    SFEVAP2B.508    
     &                       -DQW_RML(I)                                   SFEVAP2B.509    
              DFQW(I,1) = DQW_RML(I) / DTRDZ_RML(I)                        SFEVAP2B.510    
              CT_RML = -DTRDZ_RML(I) * RHOKH_1(I)                          SFEVAP2B.511    
              DTSTAR(I) = DTSTAR(I) /                                      SFEVAP2B.512    
     &                     ( 1.0 - CT_RML + ASOIL_1(I)*CP*RHOKH_1(I) )     SFEVAP2B.513    
              DFTL(I,1) = CP * RHOKH_1(I) * DTSTAR(I)                      SFEVAP2B.514    
              FTL(I,1) = FTL(I,1) + DFTL(I,1)                              SFEVAP2B.515    
              DTL_RML(I) = -CT_RML * DTSTAR(I)                             SFEVAP2B.516    
              DTSTAR(I) = DTSTAR(I) + DTL_RML(I)                           SFEVAP2B.517    
            ELSE                                                           SFEVAP2B.518    
              DQW_1(I) = -DTRDZ(I,1) * ( FQW(I,2) - FQW(I,1) )             SFEVAP2B.519    
     &                       - DQW_1(I)                                    SFEVAP2B.520    
              QW(I,1) = QW(I,1) + DQW_1(I)                                 SFEVAP2B.521    
              CT_1 = -DTRDZ(I,1) * RHOKH_1(I)                              SFEVAP2B.522    
              DTSTAR(I) = DTSTAR(I) /                                      SFEVAP2B.523    
     &                      ( 1.0 - CT_1 + ASOIL_1(I)*CP*RHOKH_1(I) )      SFEVAP2B.524    
              FTL(I,1) = FTL(I,1) + CP * RHOKH_1(I) * DTSTAR(I)            SFEVAP2B.525    
              TL(I,1) = TL(I,1) - CT_1 * DTSTAR(I)                         SFEVAP2B.526    
              DTSTAR(I) = (1.0 - CT_1) * DTSTAR(I)                         SFEVAP2B.527    
            ENDIF                                                          SFEVAP2B.528    
            TSTAR(I) = TSTAR(I) + DTSTAR(I)                                SFEVAP2B.529    
   24   CONTINUE                                                           SFEVAP2B.533    
C                                                                          SFEVAP2B.535    
C-----------------------------------------------------------------------   SFEVAP2B.536    
CL 2.5 Do calculations for sea points.                                     SFEVAP2B.537    
C-----------------------------------------------------------------------   SFEVAP2B.538    
C                                                                          SFEVAP2B.539    
CMIC$ DO ALL VECTOR SHARED(P_FIELD, BL_LEVELS, P1, POINTS,                 SFEVAP2B.541    
CMIC$1  LAND_MASK, ES, ECAN, EI, ICE_FRACT, FQW, E_SEA,                    SFEVAP2B.542    
CMIC$2  TSTAR, SMLT, SICE_MLT_HTF, TIMESTEP) PRIVATE(I, TSTARMAX)          SFEVAP2B.543    
CDIR$ IVDEP                                                                SFEVAP2B.544    
! Fujitsu vectorization directive                                          GRB0F405.447    
!OCL NOVREC                                                                GRB0F405.448    
        DO 25 I=P1,P1+POINTS-1                                             SFEVAP2B.545    
          IF (.NOT.LAND_MASK(I)) THEN                                      SFEVAP2B.546    
C                                                                          SFEVAP2B.548    
C-----------------------------------------------------------------------   SFEVAP2B.549    
CL 2.5.1 Set soil and canopy evaporation amounts to zero, and set          SFEVAP2B.550    
CL       sublimation to zero for liquid sea points.                        SFEVAP2B.551    
C-----------------------------------------------------------------------   SFEVAP2B.552    
C                                                                          SFEVAP2B.553    
            ES(I) = 0.0                                                    SFEVAP2B.554    
            ECAN(I) = 0.0                                                  SFEVAP2B.555    
            EI(I) = 0.0                                                    SFEVAP2B.556    
C-----------------------------------------------------------------------   SFEVAP2B.557    
CL 2.5.3 For sea-ice points :-                                             SFEVAP2B.558    
C-----------------------------------------------------------------------   SFEVAP2B.559    
C                                                                          SFEVAP2B.560    
            IF (ICE_FRACT(I).GT.0.0) THEN                                  SFEVAP2B.561    
C                                                                          SFEVAP2B.562    
CL   (a) Set sublimation to total moisture flux less evaporation           SFEVAP2B.563    
CL       from sea.                                                         SFEVAP2B.564    
C                                                                          SFEVAP2B.565    
              EI(I) = FQW(I,1) - E_SEA(I)                                  SFEVAP2B.566    
C                                                                          SFEVAP2B.567    
CL   (b) If the surface temperature is above the melting point of          SFEVAP2B.568    
CL       sea-ice, it is reset to be equal to the melting point.            SFEVAP2B.569    
C        (Strictly speaking, the gridbox mean surface temperature is       SFEVAP2B.570    
C        adjusted so that the temperature of the icy fraction is not       SFEVAP2B.571    
C        above the melting point; unlike the rest of this brick, this      SFEVAP2B.572    
C        sub-section is valid for any ice fraction between 1 and 0         SFEVAP2B.573    
C        inclusive, and not just for fractions of precisely 1 or 0.)       SFEVAP2B.574    
CL       This corresponds to a melting of some of the sea-ice.  If         SFEVAP2B.575    
CL       requested, the heat flux associated with this melting is          SFEVAP2B.576    
CL       calculated.  This code is indifferent to whether the melting      SFEVAP2B.577    
CL       reduces the thickness or areal extent of the sea-ice - neither    SFEVAP2B.578    
CL       is altered here.                                                  SFEVAP2B.579    
C        (This too works correctly for a general ice fraction, by          SFEVAP2B.580    
C        weighting the heat flux with the ice fraction, to give the        SFEVAP2B.581    
C        gridbox mean heat flux.)                                          SFEVAP2B.582    
C                                                                          SFEVAP2B.583    
              TSTARMAX = ICE_FRACT(I)*TM + (1.0-ICE_FRACT(I))*TFS          SFEVAP2B.584    
              IF (TSTAR(I).GT.TSTARMAX) THEN                               SFEVAP2B.585    
                IF (SMLT)                                                  SFEVAP2B.586    
     +            SICE_MLT_HTF(I) = (TSTAR(I)-TSTARMAX) / (AI*TIMESTEP)    SFEVAP2B.587    
                TSTAR(I) = TSTARMAX                                        SFEVAP2B.588    
              ENDIF   ! End of TSTAR adjustment block.                     SFEVAP2B.589    
            ENDIF     ! End of liquid sea/sea-ice block.                   SFEVAP2B.590    
          ENDIF       ! End of sea point calculations.                     SFEVAP2B.595    
   25   CONTINUE                                                           SFEVAP2B.596    
C                                                                          SFEVAP2B.598    
C-----------------------------------------------------------------------   SFEVAP2B.599    
CL 3. Update heat and moisture fluxes over land due to a rapidly           SFEVAP2B.600    
CL    layer.                                                               SFEVAP2B.601    
C-----------------------------------------------------------------------   SFEVAP2B.602    
C                                                                          SFEVAP2B.603    
      DO 3 K = 1,BL_LEVELS-1                                               SFEVAP2B.604    
        KM1 = K - 1                                                        SFEVAP2B.605    
CDIR$ IVDEP                                                                SFEVAP2B.611    
! Fujitsu vectorization directive                                          GRB0F405.449    
!OCL NOVREC                                                                GRB0F405.450    
        DO L = LAND1,LAND1+LAND_PTS-1                                      SFEVAP2B.612    
          I = LAND_INDEX(L)                                                SFEVAP2B.613    
          IF ( (K .LE. NRML(I)) .AND. (NRML(I) .GE. 2) ) THEN              SFEVAP2B.615    
            TL(I,K) = TL(I,K) + DTL_RML(I)                                 SFEVAP2B.616    
            QW(I,K) = QW(I,K) + DQW_RML(I)                                 SFEVAP2B.617    
            IF ( K .GE. 2 ) THEN                                           SFEVAP2B.618    
              DFTL(I,K) = DFTL(I,KM1) - CP * DTL_RML(I) / DTRDZ(I,KM1)     SFEVAP2B.619    
              DFQW(I,K) = DFQW(I,KM1) - DQW_RML(I) / DTRDZ(I,KM1)          SFEVAP2B.620    
              FTL(I,K) = FTL(I,K) + DFTL(I,K)                              SFEVAP2B.621    
              FQW(I,K) = FQW(I,K) + DFQW(I,K)                              SFEVAP2B.622    
            ENDIF ! K .GE. 2                                               SFEVAP2B.623    
          ENDIF  ! Rapidly mixing layer exists and is more than one        SFEVAP2B.624    
C                ! model layer deep and current layer is in it.            SFEVAP2B.625    
        ENDDO ! Loop over points                                           SFEVAP2B.629    
    3 CONTINUE! Loop over levels                                           SFEVAP2B.630    
C                                                                          SFEVAP2B.631    
C-----------------------------------------------------------------------   SFEVAP2B.632    
CL 4. Diagnose temperature and/or specific humidity at screen height       SFEVAP2B.633    
CL    (1.5 metres), as requested via the STASH flags.                      SFEVAP2B.634    
C-----------------------------------------------------------------------   SFEVAP2B.635    
C                                                                          SFEVAP2B.636    
      IF (SQ1P5 .OR. ST1P5) THEN                                           SFEVAP2B.637    
        IF (SQ1P5) CALL QSAT(QS(P1),TSTAR(P1),PSTAR(P1),POINTS)            SFEVAP2B.638    
        DO 4 I=P1,P1+POINTS-1                                              SFEVAP2B.639    
          IF (ST1P5) T1P5M(I) = TSTAR(I) - GRCP*Z1P5M + CHR1P5M(I) *       SFEVAP2B.640    
     +               ( TL(I,1) - TSTAR(I) + GRCP*(Z1(I)+Z0M(I)-Z0H(I)) )   SFEVAP2B.641    
          IF (SQ1P5) Q1P5M(I) = QW(I,1) + CER1P5M(I)*( QW(I,1) - QS(I) )   SFEVAP2B.642    
    4   CONTINUE                                                           SFEVAP2B.643    
      ENDIF                                                                SFEVAP2B.644    
      IF (LTIMER) THEN                                                     ASJ1F304.389    
      CALL TIMER('SFEVAP  ',4)                                             SFEVAP2B.646    
      ENDIF                                                                ASJ1F304.390    
      RETURN                                                               SFEVAP2B.648    
      END                                                                  SFEVAP2B.649    
*ENDIF                                                                     SFEVAP2B.650