*IF DEF,A03_3A                                                             ASJ4F401.5      
C ******************************COPYRIGHT******************************    GTS2F400.4429   
C (c) CROWN COPYRIGHT 1995, METEOROLOGICAL OFFICE, All Rights Reserved.    GTS2F400.4430   
C                                                                          GTS2F400.4431   
C Use, duplication or disclosure of this code is subject to the            GTS2F400.4432   
C restrictions as set forth in the contract.                               GTS2F400.4433   
C                                                                          GTS2F400.4434   
C                Meteorological Office                                     GTS2F400.4435   
C                London Road                                               GTS2F400.4436   
C                BRACKNELL                                                 GTS2F400.4437   
C                Berkshire UK                                              GTS2F400.4438   
C                RG12 2SZ                                                  GTS2F400.4439   
C                                                                          GTS2F400.4440   
C If no contract has been raised with this copy of the code, the use,      GTS2F400.4441   
C duplication or disclosure of it is strictly prohibited.  Permission      GTS2F400.4442   
C to do so must first be obtained in writing from the Head of Numerical    GTS2F400.4443   
C Modelling at the above address.                                          GTS2F400.4444   
C ******************************COPYRIGHT******************************    GTS2F400.4445   
C                                                                          GTS2F400.4446   
C*LL  SUBROUTINE IMPL_CAL ----------------------------------------------   IMPLCA2C.3      
CLL                                                                        IMPLCA2C.4      
CLL  Purpose: Calculate increments for surface temperature and             IMPLCA2C.5      
CLL           atmospheric variables in the boundary layer, using an        IMPLCA2C.6      
CLL           implicit numerical scheme.  The tridiagonal matrices are     IMPLCA2C.7      
CLL           inverted using simple Gaussian elimination.  See written     IMPLCA2C.8      
CLL           documentation for a very detailed explanation.               IMPLCA2C.9      
CLL                                                                        IMPLCA2C.10     
CLL  Suitable for single column use; activate *IF definition IBM.          IMPLCA2C.11     
CLL                                                                        IMPLCA2C.12     
CLL  Model           Modification history from model version 3.0           IMPLCA2C.13     
CLL version  Date                                                          IMPLCA2C.14     
CLL  3.1    12/01/93 Alternative, more complete implicit numerical         IMPLCA2C.15     
CLL                  scheme made available under *IF DEF,A03_2C.           IMPLCA2C.16     
CLL  3.3  07/07/93  Only versions 1B & 2C carried forward to UM3.3         RB070793.2      
CLL  3.4  06/06/94  DEF TIMER replaced by LOGICAL LTIMER                   ASJ1F304.360    
CLL                  Argument LTIMER added                                 ASJ1F304.361    
CLL                                                 S.J.Swarbrick          ASJ1F304.362    
!     3.5    9/5/95   MPP code: Change updateable area  P.Burton           APB1F305.329    
CLL   4.1   08/05/96  decks A03_2C and A03_3B removed                      ASJ4F401.6      
CLL                                     S D Jackson                        ASJ4F401.7      
CLL   4.2   Oct. 96   T3E migration - *DEF CRAY removed                    GSS1F402.62     
CLL                                     S J Swarbrick                      GSS1F402.63     
!LL   4.3  14/01/97   MPP code : Corrected setting of polar rows           GPB1F403.89     
!LL                                                     P.Burton           GPB1F403.90     
CLL  4.5    Jul. 98  Kill the IBM specific lines (JCThil)                  AJC1F405.163    
CLL                                                                        IMPLCA2C.17     
CLL                                                                        ASJ1F304.363    
CLL  FH, RNBS  <- Programmers of some or all of previous code or changes   IMPLCA2C.18     
CLL                                                                        IMPLCA2C.19     
CLL  Programming standard: UM Documentation Paper No 4, Version 2,         IMPLCA2C.20     
CLL                        dated 18/1/90                                   IMPLCA2C.21     
CLL                                                                        IMPLCA2C.22     
CLL  System component covered: P244                                        IMPLCA2C.23     
CLL                                                                        IMPLCA2C.24     
CLL  Project task: P24                                                     IMPLCA2C.25     
CLL                                                                        IMPLCA2C.26     
CLL  Documentation: UM Documentation Paper No 24.                          IMPLCA2C.27     
CLL                                                                        IMPLCA2C.28     
C*----------------------------------------------------------------------   IMPLCA2C.29     
C*L  Arguments :-                                                          IMPLCA2C.30     

      SUBROUTINE IMPL_CAL (                                                 2,8IMPLCA2C.31     
     + P_FIELD,U_FIELD,P1,U1,                                              IMPLCA2C.32     
     + P_POINTS,U_POINTS,BL_LEVELS,ROW_LENGTH,P_ROWS,U_ROWS                IMPLCA2C.33     
     +,CDR10M,U0,V0,SU10,SV10                                              IMPLCA2C.34     
     +,ASOIL_1,DELTA_AK,DELTA_BK,FQW_LEAD,FTL_LEAD                         IMPLCA2C.35     
     +,GAMMA_RHOKE,GAMMA_RHOKEA,GAMMA_RHOKES,GAMMA_RHOKESL                 IMPLCA2C.36     
     +,GAMMA_RHOKH_RDZ,GAMMA_RHOKM_RDZUV                                   IMPLCA2C.37     
     +,ICE_FRACT,LAND_MASK,LYING_SNOW,PSTAR,RADNET                         IMPLCA2C.38     
     +,SEA_ICE_HTF,SOIL_HT_FLUX,TIMESTEP,TSTAR_NL                          IMPLCA2C.39     
     +,EA,ES,ESL,FQW,FTL,QW,GAMMA_RHOKH_1,GAMMA_RHOKM_1,TSTAR,TL,U,V       IMPLCA2C.40     
     +,DQW_1,DQW_RML,DTRDZ,DTRDZ_RML,DTSTAR,E_SEA,H_SEA,TAUX,TAUY          IMPLCA2C.41     
     +,U10M,V10M                                                           IMPLCA2C.42     
     +,NRML                                                                IMPLCA2C.43     
     +,ERROR,LTIMER                                                        ASJ1F304.364    
     +)                                                                    IMPLCA2C.45     
      IMPLICIT NONE                                                        IMPLCA2C.46     
      LOGICAL LTIMER                                                       ASJ1F304.365    
      INTEGER                                                              IMPLCA2C.47     
     + P_FIELD                     ! IN No. of points in P-grid.           IMPLCA2C.48     
     +,U_FIELD                     ! IN No. of points in UV-grid.          IMPLCA2C.49     
     +,P1                          ! IN First point to be processed in     IMPLCA2C.50     
C                                  !    P-grid.                            IMPLCA2C.51     
     +,U1                          ! IN First point to be processed in     IMPLCA2C.52     
C                                  !    UV-grid.                           IMPLCA2C.53     
     +,P_POINTS                    ! IN Number of P-grid points to be      IMPLCA2C.54     
C                                  !    processed.                         IMPLCA2C.55     
     +,U_POINTS                    ! IN Number of UV-grid points.          IMPLCA2C.59     
     +,BL_LEVELS                   ! IN No. of atmospheric levels for      IMPLCA2C.63     
C                                  !    which boundary layer fluxes are    IMPLCA2C.64     
C                                  !    calculated.                        IMPLCA2C.65     
     +,ROW_LENGTH                  ! IN No. of points in latitude row.     IMPLCA2C.66     
     +,P_ROWS                      ! IN No. of P-rows of data to be        IMPLCA2C.70     
C                                  !    processed.                         IMPLCA2C.71     
     +,U_ROWS                      ! IN No. of UV-rows of data to be       IMPLCA2C.75     
C                                  !    processed.                         IMPLCA2C.76     
      REAL                                                                 IMPLCA2C.80     
     + CDR10M(U_FIELD)             ! IN Used in calc of 10m wind - from    IMPLCA2C.81     
C                                  !    P243 (routine SF_EXCH).  First     IMPLCA2C.82     
C                                  !    and last rows are "missing data"   IMPLCA2C.83     
C                                  !    and not used.  UV-grid.            IMPLCA2C.84     
     +,U0(U_FIELD)                 ! IN Westerly component of surface      IMPLCA2C.85     
C                                  !    current (m/s; 0 over land). UVG.   IMPLCA2C.86     
     +,V0(U_FIELD)                 ! IN Southerly component of surface     IMPLCA2C.87     
C                                  !    current (m/s; 0 over land). UVG.   IMPLCA2C.88     
      REAL                                                                 IMPLCA2C.89     
     + ASOIL_1(P_FIELD)            ! IN Soil thermodynamic coefficient     IMPLCA2C.90     
C                                  !    from P242 (SOIL_HTF).  Sq m K      IMPLCA2C.91     
C                                  !    per Joule * timestep.              IMPLCA2C.92     
     +,DELTA_AK(BL_LEVELS)         ! IN Difference of hybrid 'A' across    IMPLCA2C.93     
C                                  !    boundary layers (K-1/2 to K+1/2)   IMPLCA2C.94     
C                                  !    (upper minus lower).               IMPLCA2C.95     
     +,DELTA_BK(BL_LEVELS)         ! IN Difference of hybrid 'B' across    IMPLCA2C.96     
C                                  !    boundary layers (K-1/2 to K+1/2)   IMPLCA2C.97     
C                                  !    (upper minus lower).               IMPLCA2C.98     
     +,FQW_LEAD(P_FIELD)           ! IN "Explicit" surface flux of         IMPLCA2C.99     
C                                  !     QW (i.e. evaporation), on         IMPLCA2C.100    
C                                  !     P-grid for leads fraction         IMPLCA2C.101    
C                                  !     of gridsquare (kg/m2/s).          IMPLCA2C.102    
C                                  !     Missing data at non-sea-ice pts   IMPLCA2C.103    
     +,FTL_LEAD(P_FIELD)           ! IN "Explicit surface flux of TL       IMPLCA2C.104    
C                                  !     = H/CP (sensible heat/Cp) for     IMPLCA2C.105    
C                                  !     leads fraction of gridsquare.     IMPLCA2C.106    
C                                  !     Missing data at non-sea-ice pts   IMPLCA2C.107    
     +,GAMMA_RHOKE(P_FIELD)        ! IN Surface exchange coeff. for FQW,   IMPLCA2C.108    
C                                  !    =GAMMA(1)*RHOKE from SF_EXCH.      IMPLCA2C.109    
     +,GAMMA_RHOKEA(P_FIELD)       ! IN Surface exchange coeff. for EA,    IMPLCA2C.110    
C                                  !    =GAMMA(1)*RHOKEA from SF_EXCH.     IMPLCA2C.111    
     +,GAMMA_RHOKES(P_FIELD)       ! IN Surface exchange coeff. for ES,    IMPLCA2C.112    
C                                  !    =GAMMA(1)*RHOKES from SF_EXCH.     IMPLCA2C.113    
     +,GAMMA_RHOKESL(P_FIELD)      ! IN Surface exchange coeff. for ESL,   IMPLCA2C.114    
C                                  !    =GAMMA(1)*RHOKESL from SF_EXCH.    IMPLCA2C.115    
     +,GAMMA_RHOKH_RDZ(P_FIELD,2:BL_LEVELS)                                IMPLCA2C.116    
C                                  ! IN Exchange coeff for FTL above       IMPLCA2C.117    
C                                  !    surface, =GAMMA(K)*RHOKH(,K)       IMPLCA2C.118    
C                                  !    *RDZ(K) for K>=2 (from KMKH).      IMPLCA2C.119    
     +,GAMMA_RHOKM_RDZUV(U_FIELD,2:BL_LEVELS)                              IMPLCA2C.120    
C                                  ! IN Exchange coefficients for          IMPLCA2C.121    
C                                  !    momentum, on UV-grid with          IMPLCA2C.122    
C                                  !    first and last rows ignored.       IMPLCA2C.123    
C                                  !     =GAMMA(K)*RHOKM(,K)*RDZUV(,K)     IMPLCA2C.124    
C                                  !    for K>=2 (from KMKH).              IMPLCA2C.125    
      REAL                         ! Split to avoid > 19 continuations.    IMPLCA2C.126    
     + ICE_FRACT(P_FIELD)          ! IN Fraction of grid-box which is      IMPLCA2C.127    
C                                  !    sea-ice (decimal fraction).        IMPLCA2C.128    
     +,LYING_SNOW(P_FIELD)         ! IN Lying snow (kg per sq m ie "mm")   IMPLCA2C.129    
     +,PSTAR(P_FIELD)              ! IN Surface pressure (Pa).             IMPLCA2C.130    
     +,RADNET(P_FIELD)             ! IN Area weighted ice component        IMPLCA2C.131    
C                                  !    of surface net radiation.          IMPLCA2C.132    
C                                  !    (+ve downwards, W per sq m)        IMPLCA2C.133    
     +,SEA_ICE_HTF(P_FIELD)        ! IN Heat flux through sea-ice (W per   IMPLCA2C.134    
C                                  !    sq m, +ve downwards).  From P241   IMPLCA2C.135    
C                                  !    (routine SICE_HTF).                IMPLCA2C.136    
     +,SOIL_HT_FLUX(P_FIELD)       ! IN Heat flux from soil layer 1 to     IMPLCA2C.137    
C                                  !    soil layer 2 (deep soil layer 1)   IMPLCA2C.138    
C                                  !    i.e. +ve is downwards (W per sq    IMPLCA2C.139    
C                                  !    m). From P242 (SOIL_HTF).          IMPLCA2C.140    
     +,TIMESTEP                    ! IN Timestep in seconds.               IMPLCA2C.141    
     +,TSTAR_NL(P_FIELD)           ! IN TSTAR No Leads: surface temp-      IMPLCA2C.142    
C                                  !    erature of the ice at sea-ice      IMPLCA2C.143    
C                                  !    points and =TSTAR elsewhere.(K)    IMPLCA2C.144    
      LOGICAL                                                              IMPLCA2C.145    
     + LAND_MASK(P_FIELD)          ! IN T for land, F elsewhere.           IMPLCA2C.146    
     +,SU10                        ! IN STASH flag for 10-metre W wind.    IMPLCA2C.147    
     +,SV10                        ! IN STASH flag for 10-metre S wind.    IMPLCA2C.148    
C                                                                          IMPLCA2C.149    
C  Next 5 arrays are all IN as "explicit" fluxes from P243 (SF_EXCH and    IMPLCA2C.150    
C  possibly KMKH), and OUT as "implicit" fluxes.                           IMPLCA2C.151    
C                                                                          IMPLCA2C.152    
      REAL                                                                 IMPLCA2C.153    
     + EA(P_FIELD)                 ! INOUT Surface evaporation with only   IMPLCA2C.154    
C                                  !       aerodynamic resistance (+ve),   IMPLCA2C.155    
C                                  !       or condensation (-ve),          IMPLCA2C.156    
C                                  !       averaged over gridbox.          IMPLCA2C.157    
C                                  !       Kg per sq m per sec.            IMPLCA2C.158    
     +,ES(P_FIELD)                 ! INOUT Surface evapotranspiration      IMPLCA2C.159    
C                                  !       (through a not-entirely         IMPLCA2C.160    
C                                  !       aerodynamic resistance).        IMPLCA2C.161    
C                                  !       Always non-negative.  Kg per    IMPLCA2C.162    
C                                  !       sq m per sec.                   IMPLCA2C.163    
     +,ESL(P_FIELD)                ! INOUT ES without fractional           IMPLCA2C.164    
C                                  !       weighting factor FRACS.  "L"    IMPLCA2C.165    
C                                  !       is for "local".  Kg/sq m/sec.   IMPLCA2C.166    
     +,FQW(P_FIELD,BL_LEVELS)      ! INOUT Flux of QW (ie., for surface,   IMPLCA2C.167    
C                                  !       total evaporation). Kg/sq m/s   IMPLCA2C.168    
     +,FTL(P_FIELD,BL_LEVELS)      ! INOUT Flux of TL (ie., for surface,   IMPLCA2C.169    
C                                  !       H/Cp where H is sensible heat   IMPLCA2C.170    
C                                  !       in W per sq m).                 IMPLCA2C.171    
     +,QW(P_FIELD,BL_LEVELS)       ! INOUT Total water content (kg per     IMPLCA2C.172    
C                                  !       kg air).  From P243.            IMPLCA2C.173    
     +,GAMMA_RHOKH_1(P_FIELD)      ! IN    Surface exchange coeffs for     IMPLCA2C.174    
C                                  !       FTL, =GAMMA(1)*RHOKH(,1)        IMPLCA2C.175    
C                                  !       from P243 (SF_EXCH).            IMPLCA2C.176    
C                                  !   OUT =RHOKH_1 to satisfy STASH.      IMPLCA2C.177    
     +,GAMMA_RHOKM_1(U_FIELD)      ! IN    Surface exchange coeffs for     IMPLCA2C.178    
C                                  !       momentum, on UV-grid            IMPLCA2C.179    
C                                  !       with first and last rows        IMPLCA2C.180    
C                                  !       ignored. =GAMMA(1)*RHOKM(,1)    IMPLCA2C.181    
C                                  !       from P243 (SF_EXCH).            IMPLCA2C.182    
C                                  !   OUT =RHOKM_1 to satisfy STASH.      IMPLCA2C.183    
     +,TSTAR(P_FIELD)              ! INOUT Gridbox mean surface            IMPLCA2C.184    
C                                  !       temperature.                    IMPLCA2C.185    
     +,TL(P_FIELD,BL_LEVELS)       ! INOUT Liquid/frozen water             IMPLCA2C.186    
C                                  !       temperature (K).  From P243.    IMPLCA2C.187    
     +,U(U_FIELD,BL_LEVELS)        ! INOUT Westerly wind (m/s).  On UV-    IMPLCA2C.188    
C                                  !       grid, 1st & last rows unused.   IMPLCA2C.189    
     +,V(U_FIELD,BL_LEVELS)        ! INOUT Southerly wind (m/s).  On UV-   IMPLCA2C.190    
C                                  !       grid, 1st & last rows unused.   IMPLCA2C.191    
      REAL                                                                 IMPLCA2C.192    
     + DQW_1(P_FIELD)              ! OUT Increment for QW(,1) (needed      IMPLCA2C.193    
C                                  !     again in P245).                   IMPLCA2C.194    
     +,DQW_RML(P_FIELD)            ! OUT Rapidly mixing layer increment    IMPLCA2C.195    
C                                  !     to QW (needed again in P245).     IMPLCA2C.196    
     +,DTRDZ(P_FIELD,BL_LEVELS)    ! OUT -g.dt/dp for bottom BL_LEVELS     IMPLCA2C.197    
C                                  !     model layers (needed in P245).    IMPLCA2C.198    
     +,DTRDZ_RML(P_FIELD)          ! OUT -g.dt/dp for the rapidly          IMPLCA2C.199    
C                                  !     mixing layer (needed in P245).    IMPLCA2C.200    
     +,DTSTAR(P_FIELD)             ! OUT Surface temperature increment     IMPLCA2C.201    
C                                  !     (needed again in P245).           IMPLCA2C.202    
     +,E_SEA(P_FIELD)              ! OUT Evaporation from sea times        IMPLCA2C.203    
C                                  !     leads fraction. Zero over land    IMPLCA2C.204    
C                                  !     (kg per square metre per sec).    IMPLCA2C.205    
     +,H_SEA(P_FIELD)              ! OUT Surface sensible heat flux        IMPLCA2C.206    
C                                  !     over sea times leads fraction.    IMPLCA2C.207    
C                                  !     Zero over land. ( W/m2)           IMPLCA2C.208    
     +,TAUX(U_FIELD,BL_LEVELS)     ! INOUT x-component of turbulent        IMPLCA2C.209    
C                                  !       stress at levels k-1/2;         IMPLCA2C.210    
C                                  !     eg. TAUX(,1) is surface stress.   IMPLCA2C.211    
C                                  !     UV-grid, 1st and last rows set    IMPLCA2C.212    
C                                  !     to "missing data". (N/sq m)       IMPLCA2C.213    
     +,TAUY(U_FIELD,BL_LEVELS)     ! INOUT y-component of turbulent        IMPLCA2C.214    
C                                  !       stress at levels k-1/2;         IMPLCA2C.215    
C                                  !     eg. TAUY(,1) is surface stress.   IMPLCA2C.216    
C                                  !     UV-grid, 1st and last rows set    IMPLCA2C.217    
C                                  !     to "missing data". (N/sq m)       IMPLCA2C.218    
     +,U10M(U_FIELD)               ! OUT Westerly wind at 10m (m/s).       IMPLCA2C.219    
C                                  !     1st & last rows "missing data".   IMPLCA2C.220    
     +,V10M(U_FIELD)               ! OUT Southerly wind at 10m (m/s).      IMPLCA2C.221    
C                                  !     1st & last rows "missing data".   IMPLCA2C.222    
      INTEGER                                                              IMPLCA2C.223    
     + NRML(P_FIELD)               ! IN The number of model layers         IMPLCA2C.224    
C                                  !    in the unstable rapidly mixing     IMPLCA2C.225    
C                                  !    layer. Zero if surface layer       IMPLCA2C.226    
C                                  !    is stable.                         IMPLCA2C.227    
     +,ERROR                       ! OUT 1 if bad arguments, else 0.       IMPLCA2C.228    
C*                                                                         IMPLCA2C.229    
C*L  External references :-                                                IMPLCA2C.230    
      EXTERNAL QSAT                                                        IMPLCA2C.231    
*IF -DEF,SCMA                                                              AJC1F405.164    
      EXTERNAL P_TO_UV                                                     IMPLCA2C.233    
*ENDIF                                                                     IMPLCA2C.234    
      EXTERNAL TIMER                                                       IMPLCA2C.236    
C*                                                                         IMPLCA2C.238    
C*L  Local and other symbolic constants :-                                 IMPLCA2C.239    
*CALL C_EPSLON                                                             IMPLCA2C.240    
*CALL C_G                                                                  IMPLCA2C.241    
*CALL C_LHEAT                                                              IMPLCA2C.242    
*CALL C_R_CP                                                               IMPLCA2C.243    
*CALL C_GAMMA                                                              IMPLCA2C.244    
*CALL C_SICEHC      ! Holds integer AI                                     IMPLCA2C.245    
      REAL LS                                                              IMPLCA2C.246    
      PARAMETER (                                                          IMPLCA2C.247    
     + LS=LC+LF     ! Latent heat of sublimation (J per kg).               IMPLCA2C.248    
     +)                                                                    IMPLCA2C.249    
C*                                                                         IMPLCA2C.250    
C*L Workspace :-                                                           IMPLCA2C.251    
C   6*BL_LEVELS + 2 blocks of real workspace are required.                 IMPLCA2C.252    
      REAL                                                                 IMPLCA2C.254    
     + ALPHAS(P_FIELD)             ! Partial dQsat/dT* from Clausius-      IMPLCA2C.255    
C                                  ! Clapeyron relation.                   IMPLCA2C.256    
     +,AQ_AM(U_FIELD,BL_LEVELS)    ! As AQ: "Q" elements of matrix in      IMPLCA2C.257    
C                                  ! eqn P244.79 (modified during          IMPLCA2C.258    
C                                  ! Gaussian elimination process).        IMPLCA2C.259    
C                                  ! As AM: elements of matrix in eqn      IMPLCA2C.260    
C                                  ! P244.80 (also get modified).          IMPLCA2C.261    
     +,AQ_RML(U_FIELD)             ! Matrix element for humidity in        IMPLCA2C.262    
C                                  ! rapidly mixing layer. Then briefly    IMPLCA2C.263    
C                                  ! used for DELTAP on the UV grid.       IMPLCA2C.264    
     +,AT_ATQ(P_FIELD,BL_LEVELS)   ! Elements in atmospheric T rows of     IMPLCA2C.265    
C                                  ! matrix in eqn P244.79 (modified       IMPLCA2C.266    
C                                  ! during Gaussian elimination).         IMPLCA2C.267    
     +,AT_RML(P_FIELD)             ! Matrix element for temperature in     IMPLCA2C.268    
C                                  ! rapidly mixing layer.                 IMPLCA2C.269    
     +,ATSTAR(P_FIELD)             ! Element in surface temperature row    IMPLCA2C.270    
C                                  ! of matrix in eqn P244.79 (modified    IMPLCA2C.271    
C                                  ! during Gaussian elimination).         IMPLCA2C.272    
C                                  ! Also used for saturated sp hum, and   IMPLCA2C.273    
C                                  ! in interpolations to UV-grid.         IMPLCA2C.274    
     +,DELTAP(U_FIELD,BL_LEVELS)   ! Vertical pressure difference across   IMPLCA2C.275    
C                                  ! hybrid layers (upper minus lower)     IMPLCA2C.276    
C                                  ! (Pa).                                 IMPLCA2C.277    
     +,DELTAP_RML(U_FIELD)         ! Vertical pressure difference across   IMPLCA2C.278    
C                                  ! the rapidly mixing layer (Pa).        IMPLCA2C.279    
     +,DQW_DU(U_FIELD,BL_LEVELS)   ! As DQW: delta QW elements of vector   IMPLCA2C.280    
C                                  ! on RHS, then LHS, of eqn P244.79.     IMPLCA2C.281    
C                                  ! As DU: delta U elements of vector     IMPLCA2C.282    
C                                  ! on RHS, then LHS, of eqn P244.80.     IMPLCA2C.283    
     +,DTL_DV(U_FIELD,BL_LEVELS)   ! As DTL: delta TL (for atmosphere)     IMPLCA2C.284    
C                                  ! elements of vector on RHS, then       IMPLCA2C.285    
C                                  ! LHS, of eqn P244.79.                  IMPLCA2C.286    
C                                  ! As DV: delta V elements of vector     IMPLCA2C.287    
C                                  ! on RHS, then LHS, of eqn P244.80.     IMPLCA2C.288    
     +,DTL_RML(U_FIELD)            ! Delta TL for rapidly mixing layer.    IMPLCA2C.289    
     +,DTRDZ_UV(U_FIELD,BL_LEVELS) ! -g.dt/dp for model layers             IMPLCA2C.290    
C                                  ! interpolated to the UV-grid.          IMPLCA2C.291    
     +,FQW_ICE(P_FIELD)            ! "Explicit" surface flux of QW for     IMPLCA2C.292    
C                                  ! sea-ice fraction of gridsquare.       IMPLCA2C.293    
     +,FTL_ICE(P_FIELD)            ! "Explicit" surface flux of TL for     IMPLCA2C.294    
C                                  ! sea-ice fraction of gridsquare.       IMPLCA2C.295    
C*                                                                         IMPLCA2C.340    
C  Local scalars :-                                                        IMPLCA2C.341    
      REAL                                                                 IMPLCA2C.342    
     + CTQ      ! Matrix element in P244.??, for local increments to rml   IMPLCA2C.343    
     +,CM       ! Matrix element in eqn P244.80.                           IMPLCA2C.344    
     +,CQ       ! Matrix element in "Q" row in eqn P244.79.                IMPLCA2C.345    
     +,CQ_RML   ! As above but for rapidly mixing layer increment.         IMPLCA2C.346    
     +,CT       ! Matrix element in "T" row in eqn P244.79.                IMPLCA2C.347    
     +,CT_RML   ! As above but for rapidly mixing layer increment.         IMPLCA2C.348    
     +,RBTQ     ! Reciprocal of B P244.??, for local increments to rml     IMPLCA2C.349    
     +,RBM      ! Reciprocal of BM(') (eqns P244.81, 85, 89).              IMPLCA2C.350    
     +,RBQ      ! Reciprocal of BQ(') (eqns P244.98, 101, 104).            IMPLCA2C.351    
     +,RBQ_RML  ! As above but for the rapidly mixing layer increment.     IMPLCA2C.352    
     +,RBT      ! Reciprocal of BT' (eqns P244.107, 110, 113).             IMPLCA2C.353    
     +,RBT_RML  ! As above but for the rapidly mixing layer increment.     IMPLCA2C.354    
     +,DELTDQ   ! P244.126 times GAMMA(1) (factor required in eqns         IMPLCA2C.355    
C               ! P244.122-4).                                             IMPLCA2C.356    
     +,DTRICE   ! = DTSTAR/ICE_FRACT at sea-ice points.                    IMPLCA2C.357    
     +,DTIMEG   ! TIMESTEP * G (used in several places).                   IMPLCA2C.358    
     +,LAT_HEAT ! Latent heat (either LC or LS) in section 2.2.            IMPLCA2C.359    
      INTEGER                                                              IMPLCA2C.360    
     + BLM1     ! BL_LEVELS minus 1.                                       IMPLCA2C.361    
     +,NRMLP1   ! NRML plus 1.                                             IMPLCA2C.362    
     +,I        ! Loop counter (horizontal field index).                   IMPLCA2C.363    
     +,J        ! Offset version of I.                                     IMPLCA2C.364    
     +,K        ! Loop counter (vertical index).                           IMPLCA2C.365    
     +,KM1      ! K minus 1.                                               IMPLCA2C.366    
     +,KP1      ! K plus 1.                                                IMPLCA2C.367    
*IF DEF,MPP                                                                GPB1F403.91     
! MPP Common block                                                         GPB1F403.92     
*CALL PARVARS                                                              GPB1F403.93     
*ENDIF                                                                     GPB1F403.94     
C                                                                          IMPLCA2C.368    
C-----------------------------------------------------------------------   IMPLCA2C.369    
CL  0.  Check that the scalars input to define the grid are consistent.    IMPLCA2C.370    
C       See comments to routine SF_EXCH for details.                       IMPLCA2C.371    
C-----------------------------------------------------------------------   IMPLCA2C.372    
C                                                                          IMPLCA2C.373    
      IF (LTIMER) THEN                                                     ASJ1F304.366    
      CALL TIMER('IMPLCAL ',3)                                             IMPLCA2C.375    
      ENDIF                                                                ASJ1F304.367    
      ERROR=0                                                              IMPLCA2C.377    
*IF DEF,MPP                                                                AJC1F405.165    
      IF (                                                                 AJC1F405.166    
*ELSEIF DEF,SCMA                                                           AJC1F405.167    
      IF ( U_ROWS .NE. (P_ROWS) .OR.                                       AJC1F405.168    
*ELSE                                                                      AJC1F405.169    
      IF ( U_ROWS .NE. (P_ROWS+1) .OR.                                     AJC1F405.170    
*ENDIF                                                                     AJC1F405.171    
     +     U_POINTS .NE. (U_ROWS*ROW_LENGTH) .OR.                          IMPLCA2C.386    
     +     P_POINTS .NE. (P_ROWS*ROW_LENGTH) )  THEN                       IMPLCA2C.387    
        ERROR=1                                                            IMPLCA2C.388    
        GOTO999                                                            IMPLCA2C.389    
      ENDIF                                                                IMPLCA2C.390    
      DTIMEG = TIMESTEP * G                                                IMPLCA2C.392    
      BLM1 = BL_LEVELS-1                                                   IMPLCA2C.393    
C                                                                          IMPLCA2C.394    
CL----------------------------------------------------------------------   IMPLCA2C.395    
CL (A) Calculations on P-grid.                                             IMPLCA2C.396    
CL----------------------------------------------------------------------   IMPLCA2C.397    
C                                                                          IMPLCA2C.398    
C-----------------------------------------------------------------------   IMPLCA2C.399    
CL 1.  Calculate pressure across layer (in hybrid coordinates), DELTAP,    IMPLCA2C.400    
CL     and then -gdt/dP = dt/rho*dz for use throughout section (A)         IMPLCA2C.401    
C-----------------------------------------------------------------------   IMPLCA2C.402    
C                                                                          IMPLCA2C.403    
      DO 1 K=1,BL_LEVELS                                                   IMPLCA2C.404    
        DO 11 I=P1,P1+P_POINTS-1                                           IMPLCA2C.405    
          DELTAP(I,K) = DELTA_AK(K) + PSTAR(I)*DELTA_BK(K)                 IMPLCA2C.406    
          DTRDZ(I,K) = -DTIMEG / DELTAP(I,K)                               IMPLCA2C.407    
   11   CONTINUE                                                           IMPLCA2C.408    
   1  CONTINUE                                                             IMPLCA2C.409    
C                                                                          IMPLCA2C.410    
C-----------------------------------------------------------------------   IMPLCA2C.411    
CL 2.  Calculate implicit T and Q increments due to local mixing within    IMPLCA2C.412    
CL     the rapidly mixing layer (where it exists).                         IMPLCA2C.413    
CL     The surface fluxes FTL(I,1), FQW(I,1) are used for calculating      IMPLCA2C.414    
CL     the rapidly mixing layer (rml) increments but not here.             IMPLCA2C.415    
CL     Therefore the matrix equation we must solve to find the implicit    IMPLCA2C.416    
CL     T and Q increments due to local mixing within the rml does not      IMPLCA2C.417    
CL     have a "surface" row and we can solve for the T and Q increments    IMPLCA2C.418    
CL     for K = 1 to NRML simultaneously.                                   IMPLCA2C.419    
C-----------------------------------------------------------------------   IMPLCA2C.420    
C                                                                          IMPLCA2C.421    
C-----------------------------------------------------------------------   IMPLCA2C.422    
CL 2.1 Start 'upward sweep' with lowest model layer, which will be the     IMPLCA2C.423    
CL     bottom of the rapidly mixing layer (rml) if it exists.              IMPLCA2C.424    
C-----------------------------------------------------------------------   IMPLCA2C.425    
C                                                                          IMPLCA2C.426    
      DO 21 I=P1,P1+P_POINTS-1                                             IMPLCA2C.427    
        IF (NRML(I) .GE. 2) THEN                                           IMPLCA2C.428    
C                                                                          IMPLCA2C.429    
C  "Explicit" increments due to local mixing within the rml.               IMPLCA2C.430    
C  P244.49/31 but surface flux used in rml increment calculations.         IMPLCA2C.431    
C                                                                          IMPLCA2C.432    
          DQW_DU(I,1) = -DTRDZ(I,1) * FQW(I,2)                             IMPLCA2C.433    
          DTL_DV(I,1) = -DTRDZ(I,1) * FTL(I,2)                             IMPLCA2C.434    
C                                                                          IMPLCA2C.435    
C  Define matrix elements (CTQ always zero for this case).                 IMPLCA2C.436    
C                                                                          IMPLCA2C.437    
          AT_ATQ(I,1) = -DTRDZ(I,1) * GAMMA_RHOKH_RDZ(I,2)    ! P244.28    IMPLCA2C.438    
          RBTQ = 1.0 / ( 1.0 - AT_ATQ(I,1) )    ! Reciprocal of P244.110   IMPLCA2C.439    
C                                                                          IMPLCA2C.440    
C  Now start Gaussian elimination                                          IMPLCA2C.441    
C                                                                          IMPLCA2C.442    
          DQW_DU(I,1) = RBTQ * DQW_DU(I,1)                    ! P244.102   IMPLCA2C.443    
          DTL_DV(I,1) = RBTQ * DTL_DV(I,1)                    ! P244.111   IMPLCA2C.444    
          AT_ATQ(I,1) = RBTQ * AT_ATQ(I,1)                    ! P244.112   IMPLCA2C.445    
C                                                                          IMPLCA2C.446    
C  Start calculating DELTAP_RML. Mid-level depths added in 2.2 below.      IMPLCA2C.447    
C                                                                          IMPLCA2C.448    
          DELTAP_RML(I) = DELTAP(I,1)                                      IMPLCA2C.449    
        ELSE ! No rapidly mixing layer calculations                        IMPLCA2C.450    
          DTRDZ_RML(I) = 1.E30                                             IMPLCA2C.451    
          DQW_RML(I) = 1.E30                                               IMPLCA2C.452    
          DTL_RML(I) = 1.E30                                               IMPLCA2C.453    
          AQ_RML(I) = 1.E30                                                IMPLCA2C.454    
          AT_RML(I) = 1.E30                                                IMPLCA2C.455    
          DELTAP_RML(I) = 1.E30                                            IMPLCA2C.456    
        ENDIF                                                              IMPLCA2C.457    
   21 CONTINUE                                                             IMPLCA2C.458    
C                                                                          IMPLCA2C.459    
C-----------------------------------------------------------------------   IMPLCA2C.460    
CL 2.2 Continue upward sweep through middle of the rapidly mixing layer    IMPLCA2C.461    
CL     (if it exists) and to its top. NB NRML is always <= BLM1.           IMPLCA2C.462    
C-----------------------------------------------------------------------   IMPLCA2C.463    
C                                                                          IMPLCA2C.464    
      DO 22 K=2,BLM1                                                       IMPLCA2C.465    
        KP1 = K+1                                                          IMPLCA2C.466    
        KM1 = K-1                                                          IMPLCA2C.467    
        DO 221 I=P1,P1+P_POINTS-1                                          IMPLCA2C.468    
C                                                                          IMPLCA2C.469    
C   If in the top rapidly mixing layer then do not include flux at its     IMPLCA2C.470    
C   top in the calculation ie FQW(I,NRML+1) and FTL(I,NRML+1) are not      IMPLCA2C.471    
C   included here; they are taken account of in the non-local mixing       IMPLCA2C.472    
C   through the "rapidly mixing layer".                                    IMPLCA2C.473    
C                                                                          IMPLCA2C.474    
          IF ( K .EQ. NRML(I) ) THEN                                       IMPLCA2C.475    
C                                                                          IMPLCA2C.476    
C   Add final DELTAP contribution to DELTAP_RML and then calculate         IMPLCA2C.477    
C   DTRDZ_RML.  Lower level contributions added in 2.1 and below.          IMPLCA2C.478    
C                                                                          IMPLCA2C.479    
            DELTAP_RML(I) = DELTAP_RML(I) + DELTAP(I,K)                    IMPLCA2C.480    
            DTRDZ_RML(I) = -DTIMEG / DELTAP_RML(I)                         IMPLCA2C.481    
C                                                                          IMPLCA2C.482    
C  "Explicit" flux divergence across layer giving explicit                 IMPLCA2C.483    
C  increment due to the local mixing at the top of rml.                    IMPLCA2C.484    
C                                                                          IMPLCA2C.485    
            DQW_DU(I,K) = DTRDZ(I,K) * FQW(I,K)                            IMPLCA2C.486    
            DTL_DV(I,K) = DTRDZ(I,K) * FTL(I,K)                            IMPLCA2C.487    
C                                                                          IMPLCA2C.488    
C  Define matrix elements (A always zero for this case).                   IMPLCA2C.489    
C                                                                          IMPLCA2C.490    
            CTQ = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K)           ! P244.36   IMPLCA2C.491    
            RBTQ = 1.0 / ( 1.0  - CTQ*( 1.0 + AT_ATQ(I,KM1) ) )            IMPLCA2C.492    
C                                             ... Reciprocal of P244.113   IMPLCA2C.493    
C  Now start Gaussian elimination                                          IMPLCA2C.494    
C                                                                          IMPLCA2C.495    
            DQW_DU(I,K) = RBTQ * ( DQW_DU(I,K) - CTQ*DQW_DU(I,KM1) )       IMPLCA2C.496    
            DTL_DV(I,K) = RBTQ * ( DTL_DV(I,K) - CTQ*DTL_DV(I,KM1) )       IMPLCA2C.497    
          ELSEIF (K .LT. NRML(I)) THEN                                     IMPLCA2C.498    
C                                                                          IMPLCA2C.499    
C  Add layer depths to form total rml depth.                               IMPLCA2C.500    
C                                                                          IMPLCA2C.501    
            DELTAP_RML(I) = DELTAP_RML(I) + DELTAP(I,K)                    IMPLCA2C.502    
C                                                                          IMPLCA2C.503    
C  "Explicit" flux divergence across layer giving explicit                 IMPLCA2C.504    
C  increment due to the local mixing.P244.54/38                            IMPLCA2C.505    
C                                                                          IMPLCA2C.506    
            DQW_DU(I,K) = -DTRDZ(I,K) * ( FQW(I,KP1) - FQW(I,K) )          IMPLCA2C.507    
            DTL_DV(I,K) = -DTRDZ(I,K) * ( FTL(I,KP1) - FTL(I,K) )          IMPLCA2C.508    
C                                                                          IMPLCA2C.509    
C  Define matrix elements.                                                 IMPLCA2C.510    
C                                                                          IMPLCA2C.511    
            AT_ATQ(I,K) = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,KP1) ! P244.35   IMPLCA2C.512    
            CTQ = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K)           ! P244.36   IMPLCA2C.513    
            RBTQ = 1.0 / ( 1.0  - AT_ATQ(I,K)                              IMPLCA2C.514    
     &                          - CTQ * ( 1.0 + AT_ATQ(I,KM1) ) )          IMPLCA2C.515    
C                                             ... Reciprocal of P244.113   IMPLCA2C.516    
C  Now start Gaussian elimination                                          IMPLCA2C.517    
C                                                                          IMPLCA2C.518    
            DQW_DU(I,K) = RBTQ * ( DQW_DU(I,K) - CTQ*DQW_DU(I,KM1) )       IMPLCA2C.519    
            DTL_DV(I,K) = RBTQ * ( DTL_DV(I,K) - CTQ*DTL_DV(I,KM1) )       IMPLCA2C.520    
            AT_ATQ(I,K) = RBTQ * AT_ATQ(I,K)                 ! P244.115    IMPLCA2C.521    
          ENDIF                                                            IMPLCA2C.522    
  221   CONTINUE                                                           IMPLCA2C.523    
   22 CONTINUE                                                             IMPLCA2C.524    
C                                                                          IMPLCA2C.525    
C-----------------------------------------------------------------------   IMPLCA2C.526    
CL 2.3 Downward sweep through matrix. Add implicit increments due to       IMPLCA2C.527    
CL     local mixing within the rapidly mixing layer.  Update fluxes of     IMPLCA2C.528    
CL     heat and moisture at the surface and the top-of-rml using           IMPLCA2C.529    
CL     local mixing increments for layers 1 and NRML respectively.         IMPLCA2C.530    
C-----------------------------------------------------------------------   IMPLCA2C.531    
C                                                                          IMPLCA2C.532    
      DO 23 K=BLM1,1,-1                                                    IMPLCA2C.533    
        KP1 = K + 1                                                        IMPLCA2C.534    
        DO 231 I=P1,P1+P_POINTS-1                                          IMPLCA2C.535    
          IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRML(I))) THEN                IMPLCA2C.536    
            QW(I,K) = QW(I,K) + DQW_DU(I,K)                   ! P244.128   IMPLCA2C.537    
            TL(I,K) = TL(I,K) + DTL_DV(I,K)                   ! P244.127   IMPLCA2C.538    
            FQW(I,KP1) = FQW(I,KP1)                                        IMPLCA2C.539    
     &                    + GAMMA_RHOKH_RDZ(I,KP1)*DQW_DU(I,K)             IMPLCA2C.540    
            FTL(I,KP1) = FTL(I,KP1)                                        IMPLCA2C.541    
     &                    + GAMMA_RHOKH_RDZ(I,KP1)*DTL_DV(I,K)             IMPLCA2C.542    
          ELSEIF ((NRML(I) .GE. 2) .AND. (K .LT. NRML(I))) THEN            IMPLCA2C.543    
            DQW_DU(I,K) = DQW_DU(I,K)                                      IMPLCA2C.544    
     &                           - AT_ATQ(I,K)*DQW_DU(I,KP1)  ! P244.???   IMPLCA2C.545    
            DTL_DV(I,K) = DTL_DV(I,K)                                      IMPLCA2C.546    
     &                           - AT_ATQ(I,K)*DTL_DV(I,KP1)  ! P244.???   IMPLCA2C.547    
            QW(I,K) = QW(I,K) + DQW_DU(I,K)                   ! P244.128   IMPLCA2C.548    
            TL(I,K) = TL(I,K) + DTL_DV(I,K)                   ! P244.127   IMPLCA2C.549    
          ENDIF                                                            IMPLCA2C.550    
          IF ((NRML(I) .GE. 2) .AND. (K .EQ. 1)) THEN                      IMPLCA2C.551    
            IF (LAND_MASK(I)) THEN                                         IMPLCA2C.552    
              FTL(I,1) = FTL(I,1) - GAMMA_RHOKH_1(I) * DTL_DV(I,1)         IMPLCA2C.553    
              EA(I) = EA(I) - GAMMA_RHOKEA(I)*DQW_DU(I,1)     ! P244.122   IMPLCA2C.554    
              ESL(I) = ESL(I) - GAMMA_RHOKESL(I)*DQW_DU(I,1)  ! P244.124   IMPLCA2C.555    
              ES(I) = ES(I) - GAMMA_RHOKES(I)*DQW_DU(I,1)     ! P244.123   IMPLCA2C.556    
              FQW(I,1) = EA(I) + ES(I)                        ! P244.125   IMPLCA2C.557    
            ELSEIF (ICE_FRACT(I) .GT. 0.0) THEN                            IMPLCA2C.558    
              FQW_ICE(I) = FQW(I,1) - FQW_LEAD(I)                          IMPLCA2C.559    
     &                      - GAMMA_RHOKE(I) * ICE_FRACT(I)*DQW_DU(I,1)    IMPLCA2C.560    
              FTL_ICE(I) = FTL(I,1) - FTL_LEAD(I)                          IMPLCA2C.561    
     &                      - GAMMA_RHOKH_1(I)*ICE_FRACT(I)*DTL_DV(I,1)    IMPLCA2C.562    
              FQW_LEAD(I) = FQW_LEAD(I) - GAMMA_RHOKE(I)                   IMPLCA2C.563    
     &                        * ( 1.0 - ICE_FRACT(I) ) * DQW_DU(I,1)       IMPLCA2C.564    
              FTL_LEAD(I) = FTL_LEAD(I) - GAMMA_RHOKH_1(I)                 IMPLCA2C.565    
     &                        * ( 1.0 - ICE_FRACT(I) ) * DTL_DV(I,1)       IMPLCA2C.566    
              FTL(I,1) = FTL_ICE(I) + FTL_LEAD(I)                          IMPLCA2C.567    
              FQW(I,1) = FQW_ICE(I) + FQW_LEAD(I)                          IMPLCA2C.568    
            ELSE ! ordinary sea point                                      IMPLCA2C.569    
              FTL(I,1) = FTL(I,1) - GAMMA_RHOKH_1(I)*DTL_DV(I,1)           IMPLCA2C.570    
              FQW(I,1) = FQW(I,1) - GAMMA_RHOKE(I)*DQW_DU(I,1)             IMPLCA2C.571    
            ENDIF                                                          IMPLCA2C.572    
          ENDIF                                                            IMPLCA2C.573    
  231   CONTINUE                                                           IMPLCA2C.574    
   23 CONTINUE                                                             IMPLCA2C.575    
C                                                                          IMPLCA2C.576    
C-----------------------------------------------------------------------   IMPLCA2C.577    
CL 3.  Calculate those matrix and vector elements on the LHS of eqn        IMPLCA2C.578    
CL     P244.79 which are to do with implicit solution of the moisture      IMPLCA2C.579    
CL     transport problem at the surface, above the rml (if it exists)      IMPLCA2C.580    
CL     and between all levels if it does not.                              IMPLCA2C.581    
CL     Begin with "upward sweep" through lower half of matrix).            IMPLCA2C.582    
C-----------------------------------------------------------------------   IMPLCA2C.583    
CL 3.1 Row of matrix applying to QW transport into top "boundary"          IMPLCA2C.584    
CL     layer of model atmosphere.                                          IMPLCA2C.585    
C-----------------------------------------------------------------------   IMPLCA2C.586    
C                                                                          IMPLCA2C.587    
      DO 31 I=P1,P1+P_POINTS-1                                             IMPLCA2C.588    
        DQW_DU(I,BL_LEVELS) = DTRDZ(I,BL_LEVELS) * FQW(I,BL_LEVELS)        IMPLCA2C.589    
C                                                            ... P244.58   IMPLCA2C.590    
C                                                                          IMPLCA2C.591    
        AQ_AM(I,BL_LEVELS) = -DTRDZ(I,BL_LEVELS)               ! P244.56   IMPLCA2C.592    
     +                        * GAMMA_RHOKH_RDZ(I,BL_LEVELS)               IMPLCA2C.593    
C                                                                          IMPLCA2C.594    
        RBQ = 1.0 / ( 1.0 - AQ_AM(I,BL_LEVELS) )    ! Reciprocal P244.98   IMPLCA2C.595    
        DQW_DU(I,BL_LEVELS) = RBQ * DQW_DU(I,BL_LEVELS)        ! P244.99   IMPLCA2C.596    
        AQ_AM(I,BL_LEVELS) = RBQ * AQ_AM(I,BL_LEVELS)         ! P244.100   IMPLCA2C.597    
   31 CONTINUE                                                             IMPLCA2C.598    
C                                                                          IMPLCA2C.599    
C-----------------------------------------------------------------------   IMPLCA2C.600    
CL 3.2 Rows of matrix applying to "middle of boundary layer" model         IMPLCA2C.601    
CL     layers, i.e. all but the topmost and bottom layers.                 IMPLCA2C.602    
C-----------------------------------------------------------------------   IMPLCA2C.603    
C                                                                          IMPLCA2C.604    
      DO 32 K=BLM1,2,-1                                                    IMPLCA2C.605    
        KP1=K+1                                                            IMPLCA2C.606    
        DO 321 I=P1,P1+P_POINTS-1                                          IMPLCA2C.607    
C                                                                          IMPLCA2C.608    
C  "Explicit" flux divergence across layer giving explicit QW increment.   IMPLCA2C.609    
C                                                                          IMPLCA2C.610    
          IF ( K .GT. NRML(I) ) THEN                                       IMPLCA2C.611    
            DQW_DU(I,K) = -DTRDZ(I,K) * ( FQW(I,KP1) - FQW(I,K) )          IMPLCA2C.612    
C                                                            ... P244.54   IMPLCA2C.613    
C                                                                          IMPLCA2C.614    
            CQ = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,KP1)          ! P244.52   IMPLCA2C.615    
            AQ_AM(I,K) = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K)    ! P244.51   IMPLCA2C.616    
            RBQ = 1.0 / ( 1.0 - AQ_AM(I,K) - CQ*( 1.0 + AQ_AM(I,KP1) ) )   IMPLCA2C.617    
C                       1                       2                    2 1   IMPLCA2C.618    
C                                             ... reciprocal of P244.101   IMPLCA2C.619    
C                                                                          IMPLCA2C.620    
            DQW_DU(I,K) = RBQ * ( DQW_DU(I,K) - CQ*DQW_DU(I,KP1) )         IMPLCA2C.621    
C                                                           ... P244.102   IMPLCA2C.622    
C                                                                          IMPLCA2C.623    
            AQ_AM(I,K) = RBQ * AQ_AM(I,K)                     ! P244.103   IMPLCA2C.624    
          ENDIF                                                            IMPLCA2C.625    
  321   CONTINUE                                                           IMPLCA2C.626    
   32 CONTINUE                                                             IMPLCA2C.627    
C                                                                          IMPLCA2C.628    
C-----------------------------------------------------------------------   IMPLCA2C.629    
CL 3.3 Bottom model layer QW row of matrix equation.                       IMPLCA2C.630    
C-----------------------------------------------------------------------   IMPLCA2C.631    
C                                                                          IMPLCA2C.632    
      DO 33 I=P1,P1+P_POINTS-1                                             IMPLCA2C.633    
        IF ( NRML(I) .GE. 2 ) THEN                                         IMPLCA2C.634    
C                                                                          IMPLCA2C.635    
C-----------------------------------------------------------------------   IMPLCA2C.636    
CL 3.3.1 Start calculating rapidly mixing layer increments.                IMPLCA2C.637    
C-----------------------------------------------------------------------   IMPLCA2C.638    
C                                                                          IMPLCA2C.639    
          NRMLP1 = NRML(I) + 1                                             IMPLCA2C.640    
C                                                                          IMPLCA2C.641    
C  "Explicit" QW increment for the rapidly mixing layer.                   IMPLCA2C.642    
C                                                                          IMPLCA2C.643    
          DQW_RML(I) = -DTRDZ_RML(I) * ( FQW(I,NRMLP1) - FQW(I,1) )        IMPLCA2C.644    
C                                                                          IMPLCA2C.645    
C  Define coefficients A,B,C, for implicit calculations.                   IMPLCA2C.646    
C                                                                          IMPLCA2C.647    
          AQ_RML(I) = -DTRDZ_RML(I) * GAMMA_RHOKE(I)                       IMPLCA2C.648    
          CQ_RML = -DTRDZ_RML(I) * GAMMA_RHOKH_RDZ(I,NRMLP1)               IMPLCA2C.649    
          RBQ_RML = 1.0 / ( 1.0 - AQ_RML(I)                                IMPLCA2C.650    
     &                      - CQ_RML * ( 1.0 + AQ_AM(I,NRMLP1) ) )         IMPLCA2C.651    
          DQW_RML(I) = RBQ_RML * ( DQW_RML(I)                              IMPLCA2C.652    
     &                                 - CQ_RML * DQW_DU(I,NRMLP1) )       IMPLCA2C.653    
          AQ_RML(I) = RBQ_RML * AQ_RML(I)                                  IMPLCA2C.654    
        ELSE                                                               IMPLCA2C.655    
C                                                                          IMPLCA2C.656    
C  "Explicit" increment for QW(1) when there is no rapidly mixing          IMPLCA2C.657    
C  layer or when it is only one model layer in depth.                      IMPLCA2C.658    
C                                                                          IMPLCA2C.659    
          DQW_DU(I,1) = -DTRDZ(I,1) * ( FQW(I,2) - FQW(I,1) )  ! P244.49   IMPLCA2C.660    
          AQ_AM(I,1) = -DTRDZ(I,1) * GAMMA_RHOKE(I)            ! P244.46   IMPLCA2C.661    
          CQ = -DTRDZ(I,1) * GAMMA_RHOKH_RDZ(I,2)              ! P244.47   IMPLCA2C.662    
          RBQ = 1.0 / ( 1.0 - AQ_AM(I,1) - CQ*( 1.0 + AQ_AM(I,2) ) )       IMPLCA2C.663    
C                     1                       2                  2 1       IMPLCA2C.664    
C                                             ... reciprocal of P244.104   IMPLCA2C.665    
C                                                                          IMPLCA2C.666    
          DQW_DU(I,1) = RBQ * ( DQW_DU(I,1) - CQ*DQW_DU(I,2) )! P244.105   IMPLCA2C.667    
          AQ_AM(I,1) = RBQ * AQ_AM(I,1)                       ! P244.106   IMPLCA2C.668    
        ENDIF                                                              IMPLCA2C.669    
C                                                                          IMPLCA2C.670    
C-----------------------------------------------------------------------   IMPLCA2C.671    
CL 4.  Calculate those matrix and vector elements on the LHS of eqn        IMPLCA2C.672    
CL     P244.79 which are to do with implicit solution of the heat          IMPLCA2C.673    
CL     transport problem (i.e. for surface temperature and ice/liquid      IMPLCA2C.674    
CL     water temperatures in atmospheric boundary layer), and begin the    IMPLCA2C.675    
CL     solution algorithm (perform "upward sweep" through upper half of    IMPLCA2C.676    
CL     matrix).                                                            IMPLCA2C.677    
C-----------------------------------------------------------------------   IMPLCA2C.678    
CL 4.1 Surface temperature row of matrix (different treatments for         IMPLCA2C.679    
CL     different surface types).                                           IMPLCA2C.680    
C      ALPHAS used as a temporary store for estimate of TSTAR.             IMPLCA2C.681    
C-----------------------------------------------------------------------   IMPLCA2C.682    
C                                                                          IMPLCA2C.683    
        IF (LAND_MASK(I)) THEN                                             IMPLCA2C.684    
          IF (LYING_SNOW(I).GT.0.0) THEN                                   IMPLCA2C.685    
            LAT_HEAT = LS                                                  IMPLCA2C.686    
          ELSE                                                             IMPLCA2C.687    
            LAT_HEAT = LC                                                  IMPLCA2C.688    
          ENDIF                                                            IMPLCA2C.689    
          DTSTAR(I) = ASOIL_1(I) * ( RADNET(I)                             IMPLCA2C.690    
     &                               - CP*FTL(I,1) - LAT_HEAT*FQW(I,1)     IMPLCA2C.691    
     &                               - SOIL_HT_FLUX(I) )                   IMPLCA2C.692    
C           ... P244.15 (Explicit surface temperature increment;           IMPLCA2C.693    
C                        note that ASOIL_1 contains the TIMESTEP factor)   IMPLCA2C.694    
C                                                                          IMPLCA2C.695    
          ALPHAS(I) = TSTAR_NL(I) + DTSTAR(I)                              IMPLCA2C.696    
C           ...(Estimate of TSTAR at timelevel n+1 based on the            IMPLCA2C.697    
C               explicit surface temperature increment)                    IMPLCA2C.698    
C                                                                          IMPLCA2C.699    
        ELSEIF (ICE_FRACT(I) .GT. 0.0) THEN                                IMPLCA2C.700    
          FTL_ICE(I) = FTL(I,1) - FTL_LEAD(I)                              IMPLCA2C.701    
          FQW_ICE(I) = FQW(I,1) - FQW_LEAD(I)                              IMPLCA2C.702    
          DTSTAR(I) = TIMESTEP * AI * ( RADNET(I)                          IMPLCA2C.703    
     &                              - CP*FTL_ICE(I) - LS*FQW_ICE(I)        IMPLCA2C.704    
     &                              - SEA_ICE_HTF(I) )                     IMPLCA2C.705    
C         ... P244.19 (Gridbox mean explicit surface temp. increment)      IMPLCA2C.706    
C                                                                          IMPLCA2C.707    
          ALPHAS(I) = TSTAR_NL(I) + DTSTAR(I)/ICE_FRACT(I)                 IMPLCA2C.708    
C          ...(Estimate of T*(I), surface temperature of sea-ice, at       IMPLCA2C.709    
C              timelevel n+1 based on the explicit surface                 IMPLCA2C.710    
C              temperature increment.)                                     IMPLCA2C.711    
C                                                                          IMPLCA2C.712    
        ELSE   ! Sea without sea-ice.                                      IMPLCA2C.713    
          DTSTAR(I) = 0.0   ! Surface temperature not updated.             IMPLCA2C.714    
          ALPHAS(I) = TSTAR_NL(I) + DTSTAR(I)                              IMPLCA2C.715    
C          ...(Estimate of TSTAR at timelevel n+1 based on the             IMPLCA2C.716    
C              explicit surface temperature increment)                     IMPLCA2C.717    
C                                                                          IMPLCA2C.718    
        ENDIF                                                              IMPLCA2C.719    
   33 CONTINUE                                                             IMPLCA2C.720    
      CALL QSAT (ATSTAR(P1),TSTAR_NL(P1),PSTAR(P1),P_POINTS)               IMPLCA2C.721    
C       ...Qsat for T*n (T*n(I) at sea-ice points)  temporarily            IMPLCA2C.722    
C          stored in ATSTAR.                                               IMPLCA2C.723    
C                                                                          IMPLCA2C.724    
      CALL QSAT (AT_ATQ(P1,1),ALPHAS(P1),PSTAR(P1),P_POINTS)               IMPLCA2C.725    
C       ...Qsat for (T*n + explicit deltaT*), or (T*n(I) +                 IMPLCA2C.726    
C          explicit deltaT*(I) at sea-ice points), temporarily             IMPLCA2C.727    
C          storedin AT_ATQ(*,1).                                           IMPLCA2C.728    
C                                                                          IMPLCA2C.729    
      DO 41 I = P1,P1+P_POINTS-1                                           IMPLCA2C.730    
        IF (LAND_MASK(I)) THEN                                             IMPLCA2C.731    
          IF (LYING_SNOW(I).GT.0.0) THEN                                   IMPLCA2C.732    
            LAT_HEAT = LS                                                  IMPLCA2C.733    
          ELSE                                                             IMPLCA2C.734    
            LAT_HEAT = LC                                                  IMPLCA2C.735    
          ENDIF                                                            IMPLCA2C.736    
          IF (DTSTAR(I) .GE. -2.0 .AND. DTSTAR(I) .LE. 2.0) THEN           IMPLCA2C.737    
            ALPHAS(I)=(EPSILON*LAT_HEAT*ATSTAR(I))                         IMPLCA2C.738    
     +                  / (R*TSTAR_NL(I)*TSTAR_NL(I))                      IMPLCA2C.739    
C          ...P244.13a(Clausius-Clapeyron formula for partial dQsat/dT*    IMPLCA2C.740    
C                      used for small explicit T* increments.)             IMPLCA2C.741    
C                                                                          IMPLCA2C.742    
          ELSE                                                             IMPLCA2C.743    
            ALPHAS(I) = ( AT_ATQ(I,1) - ATSTAR(I) ) / DTSTAR(I)            IMPLCA2C.744    
C           ...P244.13b(partial dQsat/dT* calculated as a finite           IMPLCA2C.745    
C               difference for large explicit T* increments.)              IMPLCA2C.746    
C                                                                          IMPLCA2C.747    
          ENDIF                                                            IMPLCA2C.748    
          ATSTAR(I) = -ASOIL_1(I) * CP * GAMMA_RHOKH_1(I)     ! P244.16    IMPLCA2C.749    
          CT = -GAMMA_RHOKE(I) * ASOIL_1(I) * LAT_HEAT                     IMPLCA2C.750    
C                   ... P244.17 (Note: ASOIL_1 includes factor TIMESTEP)   IMPLCA2C.751    
C                                                                          IMPLCA2C.752    
          IF ( NRML(I) .GE. 2 ) THEN                                       IMPLCA2C.753    
            RBT = 1.0 / ( 1.0 - ATSTAR(I) - CT * ALPHAS(I)                 IMPLCA2C.754    
     &                                      * ( 1.0 + AQ_RML(I) ) )        IMPLCA2C.755    
            DTSTAR(I) = RBT * ( DTSTAR(I) - CT * DQW_RML(I) )              IMPLCA2C.756    
          ELSE                                                             IMPLCA2C.757    
            RBT = 1.0 / ( 1.0 - ATSTAR(I) - CT * ALPHAS(I)                 IMPLCA2C.758    
     &                                         * ( 1.0 + AQ_AM(I,1) ) )    IMPLCA2C.759    
C                       1                        2                  2 1    IMPLCA2C.760    
C                                             ... Reciprocal of P244.107   IMPLCA2C.761    
C                                                                          IMPLCA2C.762    
            DTSTAR(I) = RBT * ( DTSTAR(I) - CT * DQW_DU(I,1) )             IMPLCA2C.763    
C                                          ... P244.108 (giving DTSTAR')   IMPLCA2C.764    
C                                                                          IMPLCA2C.765    
          ENDIF                                                            IMPLCA2C.766    
          ATSTAR(I) = RBT * ATSTAR(I)                         ! P244.109   IMPLCA2C.767    
        ELSEIF (ICE_FRACT(I).GT.0.0) THEN                                  IMPLCA2C.768    
          DTRICE = DTSTAR(I)/ICE_FRACT(I)                                  IMPLCA2C.769    
          IF (DTRICE .GE. -2.0 .AND. DTRICE .LE. 2.0) THEN                 IMPLCA2C.770    
            ALPHAS(I) = (EPSILON*LS*ATSTAR(I))                             IMPLCA2C.771    
     +                  / (R*TSTAR_NL(I)*TSTAR_NL(I))                      IMPLCA2C.772    
C          ...P244.13a (Clausius-Clapeyron formula for partial             IMPLCA2C.773    
C             dQsat/dT*(I) used for small explicit T*(I) increments.)      IMPLCA2C.774    
C                                                                          IMPLCA2C.775    
          ELSE                                                             IMPLCA2C.776    
            ALPHAS(I) = ( AT_ATQ(I,1) - ATSTAR(I) ) / DTRICE               IMPLCA2C.777    
C           ...P244.13b (partial dQsat/dT*(I) calculated as a finite       IMPLCA2C.778    
C               difference used for large explicit T*(I) increments.)      IMPLCA2C.779    
C                                                                          IMPLCA2C.780    
          ENDIF                                                            IMPLCA2C.781    
          ATSTAR(I) = -GAMMA_RHOKH_1(I) * TIMESTEP * AI * CP   ! P244.20   IMPLCA2C.782    
          CT = -GAMMA_RHOKE(I) * TIMESTEP * AI * LS            ! P244.21   IMPLCA2C.783    
          IF ( NRML(I) .GE. 2 ) THEN                                       IMPLCA2C.784    
            RBT = 1.0 / ( 1.0 -ATSTAR(I) - CT * ALPHAS(I) *                IMPLCA2C.785    
     &                          ( 1.0 + ICE_FRACT(I) * AQ_RML(I) ) )       IMPLCA2C.786    
            DTSTAR(I) = RBT * ( DTSTAR(I) - ICE_FRACT(I)                   IMPLCA2C.787    
     &                                       * CT * DQW_RML(I) )           IMPLCA2C.788    
          ELSE                                                             IMPLCA2C.789    
            RBT = 1.0 / ( 1.0 - ATSTAR(I) - CT * ALPHAS(I)                 IMPLCA2C.790    
     &                            * ( 1.0 + ICE_FRACT(I)*AQ_AM(I,1) ) )    IMPLCA2C.791    
C                       1           2                               2 1    IMPLCA2C.792    
C                                             ... Reciprocal of P2440.14   IMPLCA2C.793    
C                                                                          IMPLCA2C.794    
            DTSTAR(I) = RBT * ( DTSTAR(I) - ICE_FRACT(I)*CT*DQW_DU(I,1))   IMPLCA2C.795    
C                                           ... P2440.15, giving DTSTAR'   IMPLCA2C.796    
C                                                                          IMPLCA2C.797    
          ENDIF                                                            IMPLCA2C.798    
          ATSTAR(I) = RBT * ATSTAR(I) * ICE_FRACT(I)          ! P2440.16   IMPLCA2C.799    
        ELSE               ! Sea with no sea-ice.                          IMPLCA2C.800    
          ALPHAS(I) = (EPSILON*LC*ATSTAR(I))                               IMPLCA2C.801    
     +                  / (R*TSTAR_NL(I)*TSTAR_NL(I))                      IMPLCA2C.802    
C          ...P244.13a(Clausius-Clapeyron formula for partial dQsat/dT*)   IMPLCA2C.803    
          ATSTAR(I) = 0.0                                                  IMPLCA2C.804    
        ENDIF                                                              IMPLCA2C.805    
C-----------------------------------------------------------------------   IMPLCA2C.806    
CL 4.2 Lowest atmospheric layer TL row of matrix.                          IMPLCA2C.807    
C-----------------------------------------------------------------------   IMPLCA2C.808    
        IF (NRML(I) .GE. 2) THEN                                           IMPLCA2C.809    
          NRMLP1 = NRML(I) + 1                                             IMPLCA2C.810    
C                                                                          IMPLCA2C.811    
C  "Explicit" rapidly mixing layer increment for TL.                       IMPLCA2C.812    
C                                                                          IMPLCA2C.813    
          DTL_RML(I) = -DTRDZ_RML(I) * ( FTL(I,NRMLP1) - FTL(I,1) )        IMPLCA2C.814    
          AT_RML(I) = -DTRDZ_RML(I) * GAMMA_RHOKH_RDZ(I,NRMLP1)            IMPLCA2C.815    
          CT_RML = -DTRDZ_RML(I) * GAMMA_RHOKH_1(I)                        IMPLCA2C.816    
          RBT_RML = 1.0 / ( 1.0 - AT_RML(I) - CT_RML                       IMPLCA2C.817    
     &                                         * ( 1.0 + ATSTAR(I) ) )     IMPLCA2C.818    
          DTL_RML(I) = RBT_RML * ( DTL_RML(I)                              IMPLCA2C.819    
     &                                 - CT_RML * DTSTAR(I) )              IMPLCA2C.820    
          AT_RML(I) = RBT_RML * AT_RML(I)                                  IMPLCA2C.821    
        ELSE                                                               IMPLCA2C.822    
C                                                                          IMPLCA2C.823    
C  "Explicit" increment to TL(1) when there is no rapidly mixing layer     IMPLCA2C.824    
C  or it does not extend beyond the bottom model layer.                    IMPLCA2C.825    
C                                                                          IMPLCA2C.826    
          DTL_DV(I,1) = -DTRDZ(I,1) * ( FTL(I,2) - FTL(I,1) )  ! P244.31   IMPLCA2C.827    
          CT = -DTRDZ(I,1) * GAMMA_RHOKH_1(I)                  ! P244.29   IMPLCA2C.828    
          AT_ATQ(I,1) = -DTRDZ(I,1) * GAMMA_RHOKH_RDZ(I,2)     ! P244.28   IMPLCA2C.829    
          RBT = 1.0 / ( 1.0 - AT_ATQ(I,1) - CT*( 1.0 + ATSTAR(I) ) )       IMPLCA2C.830    
C                     1                        2                 2 1       IMPLCA2C.831    
C                                             ... Reciprocal of P244.110   IMPLCA2C.832    
C                                                                          IMPLCA2C.833    
          DTL_DV(I,1) = RBT * ( DTL_DV(I,1) - CT*DTSTAR(I) )  ! P244.111   IMPLCA2C.834    
          AT_ATQ(I,1) = RBT * AT_ATQ(I,1)                     ! P244.112   IMPLCA2C.835    
        ENDIF                                                              IMPLCA2C.836    
   41 CONTINUE                                                             IMPLCA2C.837    
C-----------------------------------------------------------------------   IMPLCA2C.838    
CL 4.3 Rows of matrix applying to TL transport into model layers in the    IMPLCA2C.839    
CL     "middle" of the "boundary" layer, i.e. all but the bottom layer     IMPLCA2C.840    
CL     and the top "boundary" layer.                                       IMPLCA2C.841    
C-----------------------------------------------------------------------   IMPLCA2C.842    
      DO 43 K=2,BLM1                                                       IMPLCA2C.843    
        KP1 = K+1                                                          IMPLCA2C.844    
        KM1 = K-1                                                          IMPLCA2C.845    
        DO 431 I=P1,P1+P_POINTS-1                                          IMPLCA2C.846    
C                                                                          IMPLCA2C.847    
C   "Explicit" flux divergence across layer giving explicit TL increment   IMPLCA2C.848    
C   due to mixing above rml if it exists or for all levels if it does      IMPLCA2C.849    
C   not.                                                                   IMPLCA2C.850    
C                                                                          IMPLCA2C.851    
          NRMLP1 = NRML(I) + 1                                             IMPLCA2C.852    
          IF (K .GT. NRML(I)) THEN                                         IMPLCA2C.853    
            DTL_DV(I,K) = -DTRDZ(I,K) * ( FTL(I,KP1) - FTL(I,K) )          IMPLCA2C.854    
C                                                            ... P244.38   IMPLCA2C.855    
C                                                                          IMPLCA2C.856    
            AT_ATQ(I,K) = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,KP1) ! P244.35   IMPLCA2C.857    
            CT = -DTRDZ(I,K) * GAMMA_RHOKH_RDZ(I,K)            ! P244.36   IMPLCA2C.858    
            IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRMLP1)) THEN               IMPLCA2C.859    
              RBT = 1.0 / ( 1.0 - AT_ATQ(I,K) - CT*( 1.0 + AT_RML(I) ) )   IMPLCA2C.860    
              DTL_DV(I,K) = RBT * ( DTL_DV(I,K) - CT*DTL_RML(I) )          IMPLCA2C.861    
            ELSE                                                           IMPLCA2C.862    
              RBT = 1.0 / ( 1.0 - AT_ATQ(I,K)                              IMPLCA2C.863    
     &                          - CT*( 1.0 + AT_ATQ(I,KM1) ) )             IMPLCA2C.864    
C                         1          2                     2 1             IMPLCA2C.865    
C                                             ... Reciprocal of P244.113   IMPLCA2C.866    
C                                                                          IMPLCA2C.867    
              DTL_DV(I,K) = RBT * ( DTL_DV(I,K) - CT*DTL_DV(I,KM1) )       IMPLCA2C.868    
C                                                           ... P244.114   IMPLCA2C.869    
C                                                                          IMPLCA2C.870    
            ENDIF                                                          IMPLCA2C.871    
            AT_ATQ(I,K) = RBT * AT_ATQ(I,K)                   ! P244.115   IMPLCA2C.872    
          ENDIF                                                            IMPLCA2C.873    
  431   CONTINUE                                                           IMPLCA2C.874    
   43 CONTINUE                                                             IMPLCA2C.875    
C-----------------------------------------------------------------------   IMPLCA2C.876    
CL 4.4 Top "boundary" layer TL row of matrix.  TL for this layer can       IMPLCA2C.877    
CL     then be, and is, updated.                                           IMPLCA2C.878    
C-----------------------------------------------------------------------   IMPLCA2C.879    
      DO 44 I=P1,P1+P_POINTS-1                                             IMPLCA2C.880    
        DTL_DV(I,BL_LEVELS) = DTRDZ(I,BL_LEVELS) * FTL(I,BL_LEVELS)        IMPLCA2C.881    
C                                                            ... P244.42   IMPLCA2C.882    
C                                                                          IMPLCA2C.883    
        CT = -DTRDZ(I,BL_LEVELS) * GAMMA_RHOKH_RDZ(I,BL_LEVELS)! P244.40   IMPLCA2C.884    
        IF (NRML(I) .EQ. BLM1) THEN                                        IMPLCA2C.885    
          RBT = 1.0 / ( 1.0 - CT*( 1.0 + AT_RML(I) ) )                     IMPLCA2C.886    
          DTL_DV(I,BL_LEVELS) = RBT * ( DTL_DV(I,BL_LEVELS)                IMPLCA2C.887    
     +                                  - CT*DTL_RML(I) )                  IMPLCA2C.888    
        ELSE                                                               IMPLCA2C.889    
         RBT = 1.0 / ( 1.0 - CT*( 1.0 + AT_ATQ(I,BLM1) ) )                 IMPLCA2C.890    
C                    1          2                      2 1                 IMPLCA2C.891    
C                                             ... Reciprocal of P244.116   IMPLCA2C.892    
C                                                                          IMPLCA2C.893    
          DTL_DV(I,BL_LEVELS) = RBT * ( DTL_DV(I,BL_LEVELS)   ! P244.117   IMPLCA2C.894    
     +                                  - CT*DTL_DV(I,BLM1) )              IMPLCA2C.895    
        ENDIF                                                              IMPLCA2C.896    
        TL(I,BL_LEVELS) = TL(I,BL_LEVELS) + DTL_DV(I,BL_LEVELS)            IMPLCA2C.897    
C                                                           ... P244.127   IMPLCA2C.898    
C                                                                          IMPLCA2C.899    
   44 CONTINUE                                                             IMPLCA2C.900    
C                                                                          IMPLCA2C.901    
C-----------------------------------------------------------------------   IMPLCA2C.902    
CL 5.  "Downward sweep" through whole matrix.  TL, QW and TSTAR are        IMPLCA2C.903    
CL     updated when the final implicit increments have been calculated.    IMPLCA2C.904    
C-----------------------------------------------------------------------   IMPLCA2C.905    
CL 5.1 Remaining TL rows of matrix and add implicit increments above       IMPLCA2C.906    
CL     the rml or at all levels if it is less than two layers thick.       IMPLCA2C.907    
C-----------------------------------------------------------------------   IMPLCA2C.908    
C                                                                          IMPLCA2C.909    
      DO 51 K=BLM1,1,-1                                                    IMPLCA2C.910    
        DO 511 I=P1,P1+P_POINTS-1                                          IMPLCA2C.911    
          IF ( (K .GT. NRML(I)) .OR. (NRML(I) .LT. 2) ) THEN               IMPLCA2C.912    
            DTL_DV(I,K) = DTL_DV(I,K) - AT_ATQ(I,K)*DTL_DV(I,K+1)          IMPLCA2C.913    
                                                              ! P244.118   IMPLCA2C.914    
            TL(I,K) = TL(I,K) + DTL_DV(I,K)                   ! P244.127   IMPLCA2C.915    
          ENDIF                                                            IMPLCA2C.916    
  511   CONTINUE                                                           IMPLCA2C.917    
   51 CONTINUE                                                             IMPLCA2C.918    
      DO 52 I=P1,P1+P_POINTS-1                                             IMPLCA2C.919    
C                                                                          IMPLCA2C.920    
C-----------------------------------------------------------------------   IMPLCA2C.921    
CL 5.2 Surface temperature (TSTAR) row of matrix and rapidly mixing        IMPLCA2C.922    
CL     layer increments.                                                   IMPLCA2C.923    
C-----------------------------------------------------------------------   IMPLCA2C.924    
        IF ( NRML(I) .GE. 2 ) THEN                                         IMPLCA2C.925    
          NRMLP1 = NRML(I) + 1                                             IMPLCA2C.926    
          DTL_RML(I) = DTL_RML(I) - AT_RML(I) * DTL_DV(I,NRMLP1)           IMPLCA2C.927    
          DTSTAR(I) = DTSTAR(I) - ATSTAR(I) * DTL_RML(I)                   IMPLCA2C.928    
          DQW_RML(I) = DQW_RML(I)                                          IMPLCA2C.929    
     &                      - AQ_RML(I) * ALPHAS(I) * DTSTAR(I)            IMPLCA2C.930    
          TL(I,1) = TL(I,1) + DTL_RML(I)                                   IMPLCA2C.931    
          QW(I,1) = QW(I,1) + DQW_RML(I)                                   IMPLCA2C.932    
        ELSE                                                               IMPLCA2C.933    
          DTSTAR(I) = DTSTAR(I) - ATSTAR(I)*DTL_DV(I,1)       ! P244.119   IMPLCA2C.934    
C                                                                          IMPLCA2C.935    
C-----------------------------------------------------------------------   IMPLCA2C.936    
CL 5.3 Lowest-level QW row of matrix; local mixing where there is no rml   IMPLCA2C.937    
C-----------------------------------------------------------------------   IMPLCA2C.938    
                                                                           IMPLCA2C.939    
          DQW_DU(I,1) = DQW_DU(I,1) - AQ_AM(I,1)*ALPHAS(I)*DTSTAR(I)       IMPLCA2C.940    
C                                                           ... P244.120   IMPLCA2C.941    
          QW(I,1) = QW(I,1) + DQW_DU(I,1)     ! P244.128                   IMPLCA2C.942    
          DQW_1(I) = DQW_DU(I,1)                                           IMPLCA2C.943    
        ENDIF                                                              IMPLCA2C.944    
        TSTAR(I) = TSTAR(I) + DTSTAR(I)                                    IMPLCA2C.945    
   52 CONTINUE                                                             IMPLCA2C.946    
C                                                                          IMPLCA2C.947    
C-----------------------------------------------------------------------   IMPLCA2C.948    
CL 5.4 Remaining QW rows of matrix + updating of QW's.                     IMPLCA2C.949    
CL     Add implicit increments due to mixing above rml or at all levels    IMPLCA2C.950    
CL     if there it does not exist.                                         IMPLCA2C.951    
C-----------------------------------------------------------------------   IMPLCA2C.952    
      DO 54 K=2,BL_LEVELS                                                  IMPLCA2C.953    
        DO 541 I=P1,P1+P_POINTS-1                                          IMPLCA2C.954    
          IF (K .GT. NRML(I)) THEN                                         IMPLCA2C.955    
            IF ((NRML(I) .GE. 2) .AND. (K-1 .EQ. NRML(I))) THEN            IMPLCA2C.956    
              DQW_DU(I,K) = DQW_DU(I,K) - AQ_AM(I,K)*DQW_RML(I)            IMPLCA2C.957    
            ELSE                                                           IMPLCA2C.958    
              DQW_DU(I,K) = DQW_DU(I,K) - AQ_AM(I,K)*DQW_DU(I,K-1)         IMPLCA2C.959    
C                                                           ... P244.121   IMPLCA2C.960    
            ENDIF                                                          IMPLCA2C.961    
            QW(I,K) = QW(I,K) + DQW_DU(I,K)                  ! P244.128    IMPLCA2C.962    
          ELSE                                                             IMPLCA2C.963    
C                                                                          IMPLCA2C.964    
C  Add the increments due to rapid mixing if in the rapidly mixing layer   IMPLCA2C.965    
C                                                                          IMPLCA2C.966    
            TL(I,K) = TL(I,K) + DTL_RML(I)                                 IMPLCA2C.967    
            QW(I,K) = QW(I,K) + DQW_RML(I)                                 IMPLCA2C.968    
          ENDIF                                                            IMPLCA2C.969    
  541   CONTINUE                                                           IMPLCA2C.970    
   54 CONTINUE                                                             IMPLCA2C.971    
C                                                                          IMPLCA2C.972    
C-----------------------------------------------------------------------   IMPLCA2C.973    
CL 6.  Calculate final implicit fluxes of heat and moisture.               IMPLCA2C.974    
C-----------------------------------------------------------------------   IMPLCA2C.975    
CL 6.1 Surface fluxes for the 3 surface types: land, sea-ice, ordinary     IMPLCA2C.976    
CL     sea. Pass out the value of RHOKH(,1) in GAMMA_RHOKH_1 to satisfy    IMPLCA2C.977    
CL     STASH GAMMA_RHOKH_RDZ will contain precisely that on output.        IMPLCA2C.978    
C-----------------------------------------------------------------------   IMPLCA2C.979    
C                                                                          IMPLCA2C.980    
      DO 61 I=P1,P1+P_POINTS-1                                             IMPLCA2C.981    
        IF (LAND_MASK(I)) THEN                                             IMPLCA2C.982    
          FTL_ICE(I) = 0.0                                                 IMPLCA2C.983    
          H_SEA(I) = 0.0                                                   IMPLCA2C.984    
          FQW_ICE(I) = 0.0                                                 IMPLCA2C.985    
          E_SEA(I) = 0.0                                                   IMPLCA2C.986    
          IF ( NRML(I) .GE. 2 ) THEN                                       IMPLCA2C.987    
            FTL(I,1) = FTL(I,1) - GAMMA_RHOKH_1(I) *                       IMPLCA2C.988    
     &                             ( DTL_RML(I) - DTSTAR(I) )              IMPLCA2C.989    
            DELTDQ = DQW_RML(I) - ALPHAS(I) * DTSTAR(I)                    IMPLCA2C.990    
          ELSE                                                             IMPLCA2C.991    
            FTL(I,1) = FTL(I,1) - GAMMA_RHOKH_1(I) *                       IMPLCA2C.992    
     &                             ( DTL_DV(I,1) - DTSTAR(I) )             IMPLCA2C.993    
C                 ... P244.11 (for FTL = H/Cp, rather than for H itself)   IMPLCA2C.994    
C                                                                          IMPLCA2C.995    
            DELTDQ = DQW_DU(I,1) - ALPHAS(I)*DTSTAR(I)        ! P244.126   IMPLCA2C.996    
          ENDIF                                                            IMPLCA2C.997    
          EA(I) = EA(I) - GAMMA_RHOKEA(I)*DELTDQ              ! P244.122   IMPLCA2C.998    
          ES(I) = ES(I) - GAMMA_RHOKES(I)*DELTDQ              ! P244.123   IMPLCA2C.999    
          ESL(I) = ESL(I) - GAMMA_RHOKESL(I)*DELTDQ           ! P244.124   IMPLCA2C.1000   
          FQW(I,1) = EA(I) + ES(I)                            ! P244.125   IMPLCA2C.1001   
        ELSEIF (ICE_FRACT(I).GT.0.0) THEN                                  IMPLCA2C.1002   
          IF ( NRML(I) .GE. 2 ) THEN                                       IMPLCA2C.1003   
            H_SEA(I) = CP * ( FTL_LEAD(I) - GAMMA_RHOKH_1(I)               IMPLCA2C.1004   
     &                      * ( 1.0 - ICE_FRACT(I) ) * DTL_RML(I) )        IMPLCA2C.1005   
            FTL_ICE(I) = FTL_ICE(I) - GAMMA_RHOKH_1(I)                     IMPLCA2C.1006   
     &                    * ( ICE_FRACT(I)*DTL_RML(I) - DTSTAR(I) )        IMPLCA2C.1007   
            E_SEA(I) = FQW_LEAD(I) - GAMMA_RHOKE(I)                        IMPLCA2C.1008   
     &               * (1.0 - ICE_FRACT(I)) * DQW_RML(I)                   IMPLCA2C.1009   
            DELTDQ = ICE_FRACT(I)*DQW_RML(I) - ALPHAS(I)*DTSTAR(I)         IMPLCA2C.1010   
          ELSE                                                             IMPLCA2C.1011   
            H_SEA(I) = CP * ( FTL_LEAD(I) - GAMMA_RHOKH_1(I)               IMPLCA2C.1012   
     &                     * ( 1.0 - ICE_FRACT(I) ) * DTL_DV(I,1) )        IMPLCA2C.1013   
            FTL_ICE(I) = FTL_ICE(I) - GAMMA_RHOKH_1(I)                     IMPLCA2C.1014   
     &             * ( ICE_FRACT(I)*DTL_DV(I,1) - DTSTAR(I) ) ! P2440.1    IMPLCA2C.1015   
            E_SEA(I) = FQW_LEAD(I) - GAMMA_RHOKE(I)                        IMPLCA2C.1016   
     &               * (1.0 - ICE_FRACT(I)) * DQW_DU(I,1)     ! P2440.2    IMPLCA2C.1017   
            DELTDQ = ICE_FRACT(I)*DQW_DU(I,1) - ALPHAS(I)*DTSTAR(I)        IMPLCA2C.1018   
          ENDIF                                                            IMPLCA2C.1019   
          FQW_ICE(I) = FQW_ICE(I) - GAMMA_RHOKE(I) * DELTDQ                IMPLCA2C.1020   
          FTL(I,1) = FTL_ICE(I) + H_SEA(I)/CP                              IMPLCA2C.1021   
          FQW(I,1) = FQW_ICE(I) + E_SEA(I)                                 IMPLCA2C.1022   
          EA(I) = FQW(I,1)                                                 IMPLCA2C.1023   
          ES(I) = 0.0                                                      IMPLCA2C.1024   
          ESL(I) = 0.0                                                     IMPLCA2C.1025   
        ELSE  ! ordinary sea point                                         IMPLCA2C.1026   
          IF ( NRML(I) .GE. 2) THEN                                        IMPLCA2C.1027   
            H_SEA(I) = CP*( FTL(I,1) - GAMMA_RHOKH_1(I)*DTL_RML(I) )       IMPLCA2C.1028   
            E_SEA(I) = FQW(I,1) - GAMMA_RHOKE(I)*DQW_RML(I)                IMPLCA2C.1029   
          ELSE                                                             IMPLCA2C.1030   
            H_SEA(I) = CP * ( FTL(I,1) - GAMMA_RHOKH_1(I)*DTL_DV(I,1) )    IMPLCA2C.1031   
            E_SEA(I) = FQW(I,1) - GAMMA_RHOKE(I)*DQW_DU(I,1)               IMPLCA2C.1032   
          ENDIF                                                            IMPLCA2C.1033   
          FTL_ICE(I) = 0.0                                                 IMPLCA2C.1034   
          FQW_ICE(I) = 0.0                                                 IMPLCA2C.1035   
          FTL(I,1) = H_SEA(I) / CP                                         IMPLCA2C.1036   
          FQW(I,1) = E_SEA(I)                                              IMPLCA2C.1037   
          EA(I) = FQW(I,1)                                                 IMPLCA2C.1038   
          ES(I) = 0.0                                                      IMPLCA2C.1039   
          ESL(I) = 0.0                                                     IMPLCA2C.1040   
        ENDIF                                                              IMPLCA2C.1041   
        GAMMA_RHOKH_1(I) = GAMMA_RHOKH_1(I) / GAMMA(1)                     IMPLCA2C.1042   
   61 CONTINUE                                                             IMPLCA2C.1043   
C                                                                          IMPLCA2C.1044   
C-----------------------------------------------------------------------   IMPLCA2C.1045   
CL 6.2 Fluxes at layer interfaces above the surface.                       IMPLCA2C.1046   
C-----------------------------------------------------------------------   IMPLCA2C.1047   
      DO 62 K=2,BL_LEVELS                                                  IMPLCA2C.1048   
        KM1 = K-1                                                          IMPLCA2C.1049   
        DO 621 I=P1,P1+P_POINTS-1                                          IMPLCA2C.1050   
C                                                                          IMPLCA2C.1051   
C  Calculate and store fluxes due to local mixing.                         IMPLCA2C.1052   
C  FTL(local mixing) stored in array AT,                                   IMPLCA2C.1053   
C  FQW(local mixing) stored in array AQ_AM.                                IMPLCA2C.1054   
C                                                                          IMPLCA2C.1055   
          NRMLP1 = NRML(I) + 1                                             IMPLCA2C.1056   
          IF ((NRML(I) .GE. 2) .AND. (K .EQ. NRMLP1)) THEN                 IMPLCA2C.1057   
            AT_ATQ(I,K) = FTL(I,K) - GAMMA_RHOKH_RDZ(I,K)                  IMPLCA2C.1058   
     &                              * ( DTL_DV(I,K) - DTL_RML(I) )         IMPLCA2C.1059   
            AQ_AM(I,K) = FQW(I,K) - GAMMA_RHOKH_RDZ(I,K)                   IMPLCA2C.1060   
     &                              * ( DQW_DU(I,K) - DQW_RML(I) )         IMPLCA2C.1061   
          ELSE                                                             IMPLCA2C.1062   
            AT_ATQ(I,K) = FTL(I,K) - GAMMA_RHOKH_RDZ(I,K)                  IMPLCA2C.1063   
     &                              * ( DTL_DV(I,K) - DTL_DV(I,KM1) )      IMPLCA2C.1064   
            AQ_AM(I,K) = FQW(I,K) - GAMMA_RHOKH_RDZ(I,K)                   IMPLCA2C.1065   
     &                              * ( DQW_DU(I,K) - DQW_DU(I,KM1) )      IMPLCA2C.1066   
          ENDIF                                                            IMPLCA2C.1067   
C                                                                          IMPLCA2C.1068   
C  Now calculate the implicit fluxes including both local mixing and       IMPLCA2C.1069   
C  if appropriate also the fluxes due to rapid mixing through layers.      IMPLCA2C.1070   
C                                                                          IMPLCA2C.1071   
          IF ( K .EQ. 2 ) THEN                                             IMPLCA2C.1072   
            IF ( NRML(I) .GE. 2 ) THEN                                     IMPLCA2C.1073   
              FTL(I,K) = AT_ATQ(I,K)                                       IMPLCA2C.1074   
     &                    + FTL(I,KM1) - DTL_RML(I) / DTRDZ(I,KM1)         IMPLCA2C.1075   
              FQW(I,K) = AQ_AM(I,K)                                        IMPLCA2C.1076   
     &                    + FQW(I,KM1) - DQW_RML(I) / DTRDZ(I,KM1)         IMPLCA2C.1077   
            ELSE                                                           IMPLCA2C.1078   
              FTL(I,K) = AT_ATQ(I,K)                                       IMPLCA2C.1079   
              FQW(I,K) = AQ_AM(I,K)                                        IMPLCA2C.1080   
            ENDIF                                                          IMPLCA2C.1081   
          ELSEIF ( K .LE. NRML(I) ) THEN                                   IMPLCA2C.1082   
            FTL(I,K) = AT_ATQ(I,K) - AT_ATQ(I,KM1)                         IMPLCA2C.1083   
     &                    + FTL(I,KM1) - DTL_RML(I) / DTRDZ(I,KM1)         IMPLCA2C.1084   
            FQW(I,K) = AQ_AM(I,K) - AQ_AM(I,KM1)                           IMPLCA2C.1085   
     &                    + FQW(I,KM1) - DQW_RML(I) / DTRDZ(I,KM1)         IMPLCA2C.1086   
          ELSE                                                             IMPLCA2C.1087   
            FTL(I,K) = AT_ATQ(I,K)                                         IMPLCA2C.1088   
            FQW(I,K) = AQ_AM(I,K)                                          IMPLCA2C.1089   
          ENDIF                                                            IMPLCA2C.1090   
  621   CONTINUE                                                           IMPLCA2C.1091   
   62 CONTINUE                                                             IMPLCA2C.1092   
C                                                                          IMPLCA2C.1093   
CL----------------------------------------------------------------------   IMPLCA2C.1094   
CL (B) Calculations on UV-grid.                                            IMPLCA2C.1095   
CL----------------------------------------------------------------------   IMPLCA2C.1096   
C                                                                          IMPLCA2C.1097   
C-----------------------------------------------------------------------   IMPLCA2C.1098   
CL 7.  Solve matrix equation P244.80 for implicit increments to U and V.   IMPLCA2C.1099   
C-----------------------------------------------------------------------   IMPLCA2C.1100   
C                                                                          IMPLCA2C.1101   
*IF -DEF,SCMA                                                              AJC1F405.172    
C-----------------------------------------------------------------------   IMPLCA2C.1103   
C  7.0 Interpolate to UV-grid the pressure changes across the model        IMPLCA2C.1104   
C      layers                                                              IMPLCA2C.1105   
C-----------------------------------------------------------------------   IMPLCA2C.1106   
      DO 70 K=1,BL_LEVELS                                                  IMPLCA2C.1107   
        CALL P_TO_UV(DELTAP(P1,K),AQ_RML(U1+ROW_LENGTH),                   IMPLCA2C.1108   
     +     P_POINTS,P_POINTS-ROW_LENGTH,ROW_LENGTH,P_ROWS)                 IMPLCA2C.1109   
        DO 701 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                    IMPLCA2C.1110   
          DELTAP(I,K) = AQ_RML(I)                                          IMPLCA2C.1111   
  701   CONTINUE                                                           IMPLCA2C.1112   
   70 CONTINUE                                                             IMPLCA2C.1113   
C                                                                          IMPLCA2C.1114   
C-----------------------------------------------------------------------   IMPLCA2C.1115   
CL 7.1 Initial calculations and "upward sweep".                            IMPLCA2C.1116   
CL (a) "Surface" fluxes.                                                   IMPLCA2C.1117   
C-----------------------------------------------------------------------   IMPLCA2C.1118   
C                                                                          IMPLCA2C.1119   
      DO 71 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                       IMPLCA2C.1120   
*ELSE                                                                      IMPLCA2C.1121   
C                                                                          IMPLCA2C.1122   
C-----------------------------------------------------------------------   IMPLCA2C.1123   
CL 7.1 Initial calculations and "upward sweep".                            IMPLCA2C.1124   
C-----------------------------------------------------------------------   IMPLCA2C.1125   
CL (a) "Surface" fluxes.                                                   IMPLCA2C.1126   
C-----------------------------------------------------------------------   IMPLCA2C.1127   
      DO 71 I=1,U_POINTS                                                   IMPLCA2C.1128   
*ENDIF                                                                     IMPLCA2C.1129   
        DTRDZ_UV(I,1) = -DTIMEG / DELTAP(I,1)                              IMPLCA2C.1130   
        DELTAP_RML(I) = 0.0                                                IMPLCA2C.1131   
C                                                                          IMPLCA2C.1132   
C  "Explicit" increments to U(1) and V(1) when there is no rapidly         IMPLCA2C.1133   
C  mixing layer or it does not extend beyond the bottom model layer.       IMPLCA2C.1134   
C                                                                          IMPLCA2C.1135   
        DQW_DU(I,1) = DTRDZ_UV(I,1) * ( TAUX(I,2) - TAUX(I,1) )            IMPLCA2C.1136   
C                                                            ... P244.67   IMPLCA2C.1137   
C                                                                          IMPLCA2C.1138   
        DTL_DV(I,1) = DTRDZ_UV(I,1) * ( TAUY(I,2) - TAUY(I,1) )            IMPLCA2C.1139   
C                                                            ... P244.67   IMPLCA2C.1140   
C                                                                          IMPLCA2C.1141   
        CM = -DTRDZ_UV(I,1) * GAMMA_RHOKM_1(I)               ! P244.66     IMPLCA2C.1142   
        AQ_AM(I,1) = -DTRDZ_UV(I,1) * GAMMA_RHOKM_RDZUV(I,2)   ! P244.64   IMPLCA2C.1143   
        RBM = 1.0 / ( 1.0 - AQ_AM(I,1) - CM )    ! Reciprocal of P244.81   IMPLCA2C.1144   
        DQW_DU(I,1) = RBM * DQW_DU(I,1)                        ! P244.82   IMPLCA2C.1145   
        DTL_DV(I,1) = RBM * DTL_DV(I,1)                        ! P244.83   IMPLCA2C.1146   
        AQ_AM(I,1) = RBM * AQ_AM(I,1)                          ! P244.84   IMPLCA2C.1147   
   71 CONTINUE                                                             IMPLCA2C.1148   
C                                                                          IMPLCA2C.1149   
C                                                                          IMPLCA2C.1150   
C-----------------------------------------------------------------------   IMPLCA2C.1151   
CL (b) Fluxes at (or rows representing) layer interfaces above the         IMPLCA2C.1152   
CL     surface but below the top of the boundary layer.                    IMPLCA2C.1153   
C-----------------------------------------------------------------------   IMPLCA2C.1154   
      DO 72 K=2,BLM1                                                       IMPLCA2C.1155   
        KP1 = K+1                                                          IMPLCA2C.1156   
        KM1 = K-1                                                          IMPLCA2C.1157   
*IF -DEF,SCMA                                                              AJC1F405.173    
        DO 721 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                    IMPLCA2C.1159   
*ELSE                                                                      IMPLCA2C.1160   
        DO 721 I=1,U_POINTS                                                IMPLCA2C.1161   
*ENDIF                                                                     IMPLCA2C.1162   
          DTRDZ_UV(I,K) = -DTIMEG / DELTAP(I,K)                            IMPLCA2C.1163   
            DQW_DU(I,K) = DTRDZ_UV(I,K) * ( TAUX(I,KP1) - TAUX(I,K) )      IMPLCA2C.1164   
C                                                            ... P244.74   IMPLCA2C.1165   
            DTL_DV(I,K) = DTRDZ_UV(I,K) * ( TAUY(I,KP1) - TAUY(I,K) )      IMPLCA2C.1166   
C                                                            ... P244.74   IMPLCA2C.1167   
            AQ_AM(I,K) = -DTRDZ_UV(I,K) * GAMMA_RHOKM_RDZUV(I,KP1)         IMPLCA2C.1168   
C                                                            ... P244.71   IMPLCA2C.1169   
          CM = -DTRDZ_UV(I,K) * GAMMA_RHOKM_RDZUV(I,K)         ! P244.72   IMPLCA2C.1170   
            RBM = 1.0 / ( 1.0 - AQ_AM(I,K) -CM*( 1.0 + AQ_AM(I,KM1) ) )    IMPLCA2C.1171   
C                     1                      2                    2 1      IMPLCA2C.1172   
C                                              ... Reciprocal of P244.85   IMPLCA2C.1173   
C                                                                          IMPLCA2C.1174   
            DQW_DU(I,K) = RBM * ( DQW_DU(I,K) - CM*DQW_DU(I,KM1) )         IMPLCA2C.1175   
C                                                            ... P244.86   IMPLCA2C.1176   
C                                                                          IMPLCA2C.1177   
            DTL_DV(I,K) = RBM * ( DTL_DV(I,K) - CM*DTL_DV(I,KM1) )         IMPLCA2C.1178   
C                                                            ... P244.87   IMPLCA2C.1179   
C                                                                          IMPLCA2C.1180   
          AQ_AM(I,K) = RBM * AQ_AM(I,K)                        ! P244.88   IMPLCA2C.1181   
  721   CONTINUE                                                           IMPLCA2C.1182   
   72 CONTINUE                                                             IMPLCA2C.1183   
C-----------------------------------------------------------------------   IMPLCA2C.1184   
CL (c) Top "boundary" layer; also increment U and V here, as implicit      IMPLCA2C.1185   
CL     flux for this layer is got from "upward sweep" alone.               IMPLCA2C.1186   
C-----------------------------------------------------------------------   IMPLCA2C.1187   
*IF -DEF,SCMA                                                              AJC1F405.174    
      DO 73 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                       IMPLCA2C.1189   
*ELSE                                                                      IMPLCA2C.1190   
      DO 73 I=1,U_POINTS                                                   IMPLCA2C.1191   
*ENDIF                                                                     IMPLCA2C.1192   
        DTRDZ_UV(I,BL_LEVELS) = -DTIMEG / DELTAP(I,BL_LEVELS)              IMPLCA2C.1193   
        DQW_DU(I,BL_LEVELS) = -DTRDZ_UV(I,BL_LEVELS) * TAUX(I,BL_LEVELS)   IMPLCA2C.1194   
C                                                            ... P244.78   IMPLCA2C.1195   
C                                                                          IMPLCA2C.1196   
        DTL_DV(I,BL_LEVELS) = -DTRDZ_UV(I,BL_LEVELS) * TAUY(I,BL_LEVELS)   IMPLCA2C.1197   
C                                                            ... P244.78   IMPLCA2C.1198   
C                                                                          IMPLCA2C.1199   
        CM = -DTRDZ_UV(I,BL_LEVELS) * GAMMA_RHOKM_RDZUV(I,BL_LEVELS)       IMPLCA2C.1200   
C                                                            ... P244.76   IMPLCA2C.1201   
        RBM = 1.0 / ( 1.0 - CM*( 1.0 + AQ_AM(I,BLM1) ) )                   IMPLCA2C.1202   
C                     1          2                     2 1                 IMPLCA2C.1203   
C                                              ... Reciprocal of P244.89   IMPLCA2C.1204   
C                                                                          IMPLCA2C.1205   
        DQW_DU(I,BL_LEVELS) = RBM * ( DQW_DU(I,BL_LEVELS)    ! P244.90     IMPLCA2C.1206   
     +                                    - CM*DQW_DU(I,BLM1) )            IMPLCA2C.1207   
        DTL_DV(I,BL_LEVELS) = RBM * ( DTL_DV(I,BL_LEVELS)    ! P244.91     IMPLCA2C.1208   
     +                                    - CM*DTL_DV(I,BLM1) )            IMPLCA2C.1209   
        U(I,BL_LEVELS) = U(I,BL_LEVELS) + DQW_DU(I,BL_LEVELS)  ! P244.94   IMPLCA2C.1210   
        V(I,BL_LEVELS) = V(I,BL_LEVELS) + DTL_DV(I,BL_LEVELS)  ! P244.95   IMPLCA2C.1211   
   73 CONTINUE                                                             IMPLCA2C.1212   
                                                                           IMPLCA2C.1213   
C-----------------------------------------------------------------------   IMPLCA2C.1214   
CL 7.4 Complete solution of matrix equation by performing "downward        IMPLCA2C.1215   
CL     sweep", then update U and V.                                        IMPLCA2C.1216   
C-----------------------------------------------------------------------   IMPLCA2C.1217   
      DO 74 K=BLM1,1,-1                                                    IMPLCA2C.1218   
        KP1 = K+1                                                          IMPLCA2C.1219   
*IF -DEF,SCMA                                                              AJC1F405.175    
        DO 741 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                    IMPLCA2C.1221   
*ELSE                                                                      IMPLCA2C.1222   
        DO 741 I=1,U_POINTS                                                IMPLCA2C.1223   
*ENDIF                                                                     IMPLCA2C.1224   
          DQW_DU(I,K) = DQW_DU(I,K) - AQ_AM(I,K)*DQW_DU(I,KP1) ! P244.92   IMPLCA2C.1225   
          DTL_DV(I,K) = DTL_DV(I,K) - AQ_AM(I,K)*DTL_DV(I,KP1) ! P244.93   IMPLCA2C.1226   
          U(I,K) = U(I,K) + DQW_DU(I,K)                        ! P244.94   IMPLCA2C.1227   
          V(I,K) = V(I,K) + DTL_DV(I,K)                        ! P244.95   IMPLCA2C.1228   
  741   CONTINUE                                                           IMPLCA2C.1229   
   74 CONTINUE                                                             IMPLCA2C.1230   
C                                                                          IMPLCA2C.1231   
C-----------------------------------------------------------------------   IMPLCA2C.1232   
CL 8.  Essentially diagnostic calculations, though some values (i.e. the   IMPLCA2C.1233   
CL     surface wind stresses) are required by the coupled version of the   IMPLCA2C.1234   
CL     model.                                                              IMPLCA2C.1235   
C-----------------------------------------------------------------------   IMPLCA2C.1236   
CL 8.1 Surface wind stress components.                                     IMPLCA2C.1237   
CL     Pass out the value of RHOKM(,1) in GAMMA_RHOKM_1 to satisfy STASH   IMPLCA2C.1238   
CL     GAMMA_RHOKM_RDZ will contain precisely that on output.              IMPLCA2C.1239   
C-----------------------------------------------------------------------   IMPLCA2C.1240   
C                                                                          IMPLCA2C.1241   
*IF -DEF,SCMA                                                              AJC1F405.176    
      DO 81 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                       IMPLCA2C.1243   
*ELSE                                                                      IMPLCA2C.1244   
      DO 81 I=1,U_POINTS                                                   IMPLCA2C.1245   
*ENDIF                                                                     IMPLCA2C.1246   
        TAUX(I,1) = TAUX(I,1) + GAMMA_RHOKM_1(I)*DQW_DU(I,1)               IMPLCA2C.1247   
C                                                            ... P244.61   IMPLCA2C.1248   
C                                                                          IMPLCA2C.1249   
        TAUY(I,1) = TAUY(I,1) + GAMMA_RHOKM_1(I)*DTL_DV(I,1)               IMPLCA2C.1250   
C                                                            ... P244.61   IMPLCA2C.1251   
C                                                                          IMPLCA2C.1252   
        GAMMA_RHOKM_1(I) = GAMMA_RHOKM_1(I) / GAMMA(1)                     IMPLCA2C.1253   
   81 CONTINUE                                                             IMPLCA2C.1254   
*IF -DEF,SCMA                                                              AJC1F405.177    
C                                                                          IMPLCA2C.1256   
C  Set extreme rows to "missing data indicator".                           IMPLCA2C.1257   
C                                                                          IMPLCA2C.1258   
*IF DEF,MPP                                                                GPB1F403.95     
      IF (attop) THEN                                                      GPB1F403.96     
*ENDIF                                                                     GPB1F403.97     
        DO I=U1,U1+ROW_LENGTH-1                                            GPB1F403.98     
          TAUX(I,1)=1.E30                                                  GPB1F403.99     
          TAUY(I,1)=1.E30                                                  GPB1F403.100    
        ENDDO                                                              GPB1F403.101    
*IF DEF,MPP                                                                GPB1F403.102    
      ENDIF                                                                GPB1F403.103    
                                                                           GPB1F403.104    
      IF (atbase) THEN                                                     GPB1F403.105    
*ENDIF                                                                     GPB1F403.106    
        DO I= U1 + (U_ROWS-1)*ROW_LENGTH , U1 + U_ROWS*ROW_LENGTH - 1      GPB1F403.107    
          TAUX(I,1)=1.E30                                                  GPB1F403.108    
          TAUY(I,1)=1.E30                                                  GPB1F403.109    
        ENDDO                                                              GPB1F403.110    
*IF DEF,MPP                                                                GPB1F403.111    
      ENDIF                                                                GPB1F403.112    
*ENDIF                                                                     GPB1F403.113    
*ENDIF                                                                     IMPLCA2C.1283   
C                                                                          IMPLCA2C.1284   
C-----------------------------------------------------------------------   IMPLCA2C.1285   
CL 8.2 Wind stress components at layer interfaces above the surface.       IMPLCA2C.1286   
C-----------------------------------------------------------------------   IMPLCA2C.1287   
      DO 82 K=2,BL_LEVELS                                                  IMPLCA2C.1288   
        KM1 = K-1                                                          IMPLCA2C.1289   
*IF -DEF,SCMA                                                              AJC1F405.178    
        DO 821 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                    IMPLCA2C.1291   
*ELSE                                                                      IMPLCA2C.1292   
        DO 821 I=1,U_POINTS                                                IMPLCA2C.1293   
*ENDIF                                                                     IMPLCA2C.1294   
          AQ_AM(I,K) = TAUX(I,K) + GAMMA_RHOKM_RDZUV(I,K)    ! P244.61     IMPLCA2C.1295   
     +                        *( DQW_DU(I,K) - DQW_DU(I,KM1) )             IMPLCA2C.1296   
          AT_ATQ(I,K) = TAUY(I,K) + GAMMA_RHOKM_RDZUV(I,K)   ! P244.61     IMPLCA2C.1297   
     +                        *( DTL_DV(I,K) - DTL_DV(I,KM1) )             IMPLCA2C.1298   
          TAUX(I,K) = AQ_AM(I,K)                                           IMPLCA2C.1299   
          TAUY(I,K) = AT_ATQ(I,K)                                          IMPLCA2C.1300   
  821   CONTINUE                                                           IMPLCA2C.1301   
*IF -DEF,SCMA                                                              AJC1F405.179    
C                                                                          IMPLCA2C.1303   
C  Set extreme rows to "missing data indicator".                           IMPLCA2C.1304   
C                                                                          IMPLCA2C.1305   
*IF DEF,MPP                                                                GPB1F403.114    
      IF (attop) THEN                                                      GPB1F403.115    
*ENDIF                                                                     GPB1F403.116    
        DO I=U1,U1+ROW_LENGTH-1                                            GPB1F403.117    
          TAUX(I,K)=1.E30                                                  GPB1F403.118    
          TAUY(I,K)=1.E30                                                  GPB1F403.119    
        ENDDO                                                              GPB1F403.120    
*IF DEF,MPP                                                                GPB1F403.121    
      ENDIF                                                                GPB1F403.122    
                                                                           GPB1F403.123    
      IF (atbase) THEN                                                     GPB1F403.124    
*ENDIF                                                                     GPB1F403.125    
        DO I= U1 + (U_ROWS-1)*ROW_LENGTH , U1 + U_ROWS*ROW_LENGTH - 1      GPB1F403.126    
          TAUX(I,K)=1.E30                                                  GPB1F403.127    
          TAUY(I,K)=1.E30                                                  GPB1F403.128    
        ENDDO                                                              GPB1F403.129    
*IF DEF,MPP                                                                GPB1F403.130    
      ENDIF                                                                GPB1F403.131    
*ENDIF                                                                     GPB1F403.132    
*ENDIF                                                                     IMPLCA2C.1313   
   82 CONTINUE                                                             IMPLCA2C.1314   
C                                                                          IMPLCA2C.1315   
C-----------------------------------------------------------------------   IMPLCA2C.1316   
CL 8.3 Wind components at 10 metres above the surface.                     IMPLCA2C.1317   
C-----------------------------------------------------------------------   IMPLCA2C.1318   
      IF (SU10 .OR. SV10) THEN                                             IMPLCA2C.1319   
*IF -DEF,SCMA                                                              AJC1F405.180    
CMIC$ DO ALL VECTOR SHARED(U_POINTS, ROW_LENGTH, U1, SU10, SV10, U_FIELD   IMPLCA2C.1321   
CMIC$1   , U0, CDR10M, BL_LEVELS, U, U10M, V0, V, V10M) PRIVATE(I)         IMPLCA2C.1322   
CDIR$ IVDEP                                                                IMPLCA2C.1323   
! Fujitsu vectorization directive                                          GRB0F405.361    
!OCL NOVREC                                                                GRB0F405.362    
        DO 83 I=U1+ROW_LENGTH,U1+U_POINTS-ROW_LENGTH-1                     IMPLCA2C.1324   
*ELSE                                                                      IMPLCA2C.1325   
        DO 83 I=1,U_POINTS                                                 IMPLCA2C.1326   
*ENDIF                                                                     IMPLCA2C.1327   
          IF (SU10)                                                        IMPLCA2C.1328   
     +      U10M(I) = U0(I) + CDR10M(I)*( U(I,1) - U0(I) )     ! P244.96   IMPLCA2C.1329   
          IF (SV10)                                                        IMPLCA2C.1330   
     +      V10M(I) = V0(I) + CDR10M(I)*( V(I,1) - V0(I) )     ! P244.97   IMPLCA2C.1331   
   83   CONTINUE                                                           IMPLCA2C.1332   
*IF -DEF,SCMA                                                              AJC1F405.181    
C                                                                          IMPLCA2C.1334   
C  Set extreme rows to "missing data indicator".                           IMPLCA2C.1335   
C                                                                          IMPLCA2C.1336   
*IF DEF,MPP                                                                GPB1F403.133    
      IF (attop) THEN                                                      GPB1F403.134    
*ENDIF                                                                     GPB1F403.135    
        DO I=U1,U1+ROW_LENGTH-1                                            GPB1F403.136    
          IF (SU10) THEN                                                   GPB1F403.137    
            U10M(I)=1.E30                                                  GPB1F403.138    
          ENDIF                                                            GPB1F403.139    
          IF (SV10) THEN                                                   GPB1F403.140    
            V10M(I)=1.E30                                                  GPB1F403.141    
          ENDIF                                                            GPB1F403.142    
        ENDDO                                                              GPB1F403.143    
*IF DEF,MPP                                                                GPB1F403.144    
      ENDIF                                                                GPB1F403.145    
                                                                           GPB1F403.146    
      IF (atbase) THEN                                                     GPB1F403.147    
*ENDIF                                                                     GPB1F403.148    
        DO I= U1 + (U_ROWS-1)*ROW_LENGTH , U1 + U_ROWS*ROW_LENGTH - 1      GPB1F403.149    
          IF (SU10) THEN                                                   GPB1F403.150    
            U10M(I)=1.E30                                                  GPB1F403.151    
          ENDIF                                                            GPB1F403.152    
          IF (SV10) THEN                                                   GPB1F403.153    
            V10M(I)=1.E30                                                  GPB1F403.154    
          ENDIF                                                            GPB1F403.155    
        ENDDO                                                              GPB1F403.156    
*IF DEF,MPP                                                                GPB1F403.157    
      ENDIF                                                                GPB1F403.158    
*ENDIF                                                                     GPB1F403.159    
*ENDIF                                                                     IMPLCA2C.1369   
      ENDIF                                                                IMPLCA2C.1370   
  999 CONTINUE   ! Branch for error exit.                                  IMPLCA2C.1371   
      IF (LTIMER) THEN                                                     ASJ1F304.368    
      CALL TIMER('IMPLCAL ',4)                                             IMPLCA2C.1373   
      ENDIF                                                                ASJ1F304.369    
      RETURN                                                               IMPLCA2C.1375   
      END                                                                  IMPLCA2C.1376   
*ENDIF                                                                     IMPLCA2C.1377