      SUBROUTINE FLAKE (LTSNGAT, LTICGAT, LTAVGAT, LTMXGAT, LTWBGAT,
     1                  LSHPGAT, LZSNGAT, LZICGAT, LDMXGAT, GTGAT,
     2                  SNOGAT,  ANGAT,   SALBGAT, CSALGAT, ZNGAT,
     3                  TISLGAT, SMLTGAT, ROFNGAT, ROFGAT,  LUIMGAT,
     4                  PCSNGAT,
     5                  HLAKGAT, PREGAT,  QFSGAT,  GCGAT,   RHSIGAT,
     6                  SNOWGAT, FSVGAT,  FSIGAT,  FDLGAT,  RADJGAT,
     7                  HFSGAT,  HFLGAT,  VMODGAT, CDMGAT,  RHOAGAT,
     8                  RHONGAT, RPREGAT, SPCPGAT,
     9                  DELT, RHOW, TW, DENI, NBS, ILG, NMK        )
C
C     * DRIVER FOR FLAKE.         DEC 07, 2016 - M.LAZARE.
C
      IMPLICIT NONE
C
C     * I/O FIELDS UNIQUE TO FLAKE:
C
      REAL, DIMENSION(ILG) :: LTSNGAT, LTICGAT, LTAVGAT, LTMXGAT, 
     1                        LTWBGAT, LSHPGAT, LZSNGAT, LZICGAT,
     2                        LDMXGAT
C
C     * I/O FIELDS IN COMMON WITH REST OF PHYSICS:
C
      REAL, DIMENSION(ILG,NBS) :: SALBGAT, CSALGAT
      REAL, DIMENSION(ILG) :: GTGAT,   SNOGAT,  ANGAT,   ZNGAT,
     1                        TISLGAT, SMLTGAT, ROFNGAT, ROFGAT,
     2                        LUIMGAT, PCSNGAT
C
C     * INPUT FIELDS:
C
      REAL, DIMENSION(ILG) :: HLAKGAT, PREGAT,  QFSGAT,  GCGAT,
     1                        RHSIGAT, SNOWGAT, FSVGAT,  FSIGAT, 
     2                        FDLGAT,  RADJGAT, HFSGAT,  HFLGAT,
     3                        VMODGAT, CDMGAT,  RHOAGAT, RHONGAT,
     4                        RPREGAT, SPCPGAT
C
C     * LOCAL ARRAYS/VARIABLE CONSISTENT WITH FLAKE CONSTRUCTS.
C     * "ireals" sets the number of desired significant digits
C     * for real variables, corresponding to 8-byte reals.
C
      INTEGER, PARAMETER       ::
     1         ireals    = SELECTED_REAL_KIND (12,200)

      REAL (KIND=ireals)       :: dMsnowdt, I_atm_in, 
     1                            non_solar_heat, ustar,
     2                            depth_w, depth_bs, T_bs, 
     3                            par_Coriolis, del_time,
     4                            T_snow_in,  T_ice_in,  
     5                            T_mnw_in,  T_wML_in,
     6                            T_bot_in,  T_B1_in,
     7                            C_T_in,  h_snow_in,
     8                            h_ice_in,  h_ML_in,
     9                            H_B1_in, T_sfc_p,
     A                            albedo_water, albedo_ice, albedo_snow,
     B                            T_snow_out, T_ice_out,
     C                            T_mnw_out, T_wML_out,
     D                            T_bot_out, T_B1_out,
     E                            C_T_out, h_snow_out, 
     F                            h_ice_out, h_ML_out,
     G                            H_B1_out, T_sfc_n
C
      REAL DELT,RHOW,TW,DENI,RICEMELT
      INTEGER ILG,NMK,K,NBS,L
C==================================================================
C       * FOR NOW, THE SIMPLEST WAY TO TEST IS TO CALL THE INTERFACE
C       * IN A LOOP OVER THE GATHERED LAKE (TILED) POINTS.
C       * WE WILL LIKELY MOVE THIS INSIDE LATER AS SUGGESTED IN THE
C       * DOCUMENTATION.
C
        DO K=1,NMK
c
c         * Also, we have set "lflk_botsed_use" to ".false." in
c         * module "flake_configure", so "depth_bs" and "T_bs" are not
c         * used and are set to zero here for completeness. The prognostic
c         * "T_B1" and "H_B1" are therefore not done and left as scalars.
c
c         * We remove the passing of albedoes to the Flake interface, 
c         * because we now pass FSG instead of FSS.
c
c         * dMsnowdt=0 is used as per suggestion by Flake documentation.
c         * This is due to concerns about stability of snow module.
c
C         * STORE AGCM FIELDS INTO DATA TYPES FOR FLAKE.
C
          depth_bs=0.
          T_bs=0.
          H_B1_in=0.
          T_B1_in=0.
          dMsnowdt=0.
c         dMsnowdt       = SNOWGAT(K) ! No snowfall!
          I_atm_in       = FSVGAT(K)+FSIGAT(K)
          non_solar_heat = FDLGAT(K)-HFSGAT(K)-HFLGAT(K)
          ustar          = VMODGAT(K)*SQRT(CDMGAT(K)*RHOAGAT(K)/RHOW)
          depth_w        = HLAKGAT(K)
          par_Coriolis   = TW*SIN(RADJGAT(K))
          del_time       = DELT
          T_snow_in      = LTSNGAT(K)
          T_ice_in       = LTICGAT(K)  
          T_mnw_in       = LTAVGAT(K)
          T_wML_in       = LTMXGAT(K)
          T_bot_in       = LTWBGAT(K)
          C_T_in         = LSHPGAT(K)
          h_snow_in      = LZSNGAT(K)
          h_ice_in       = LZICGAT(K)
          h_ML_in        = LDMXGAT(K)
          T_sfc_p        = GTGAT(K)
C
          CALL FLAKE_INTERFACE (dMsnowdt, I_atm_in, 
     1                          non_solar_heat, ustar,
     2                          depth_w, depth_bs, T_bs, 
     3                          par_Coriolis, del_time,
     4                          T_snow_in,  T_ice_in,  
     5                          T_mnw_in,  T_wML_in,
     6                          T_bot_in,  T_B1_in,
     7                          C_T_in,  h_snow_in,
     8                          h_ice_in,  h_ML_in,
     9                          H_B1_in, T_sfc_p,
     A                          albedo_water, albedo_ice, albedo_snow,
     B                          T_snow_out, T_ice_out,
     C                          T_mnw_out, T_wML_out,
     D                          T_bot_out, T_B1_out,
     E                          C_T_out, h_snow_out, 
     F                          h_ice_out, h_ML_out,
     G                          H_B1_out, T_sfc_n )
C
C         * CALCULATE EFFECTIVE ALBEDO FOR DIAGNOSTIC PUPROSES.
C         * THE CODE IS COPIED DIRECTLY FROM THE FLAKE ROUTINES.
C
          IF(h_ice_in.ge.1.e-9) THEN                      ! Ice exists
           IF(h_snow_in.ge.1.e-5) THEN                    ! There is snow above the ice
            ANGAT   (K) = albedo_snow
            DO L = 1,NBS
             SALBGAT(K,L) = albedo_snow
             CSALGAT(K,L) = albedo_snow
            ENDDO
           ELSE                                           ! No snow above the ice
            DO L = 1,NBS
             SALBGAT(K,L) = albedo_ice
             CSALGAT(K,L) = albedo_ice
            ENDDO
           END IF
          ELSE                                            ! No ice-snow cover
           DO L = 1,NBS
            SALBGAT(K,L) = albedo_water
            CSALGAT(K,L) = albedo_water
           ENDDO
          END IF
C
C         * STORE FLAKE OUTPUT INTO USUAL GCM-TYPE FIELDS.
C
          LTSNGAT(K) = T_snow_out
          LTICGAT(K) = T_ice_out
          LTAVGAT(K) = T_mnw_out
          LTMXGAT(K) = T_wML_out    
          LTWBGAT(K) = T_bot_out
          LSHPGAT(K) = C_T_out
          LZSNGAT(K) = h_snow_out
          LZICGAT(K) = h_ice_out
          LDMXGAT(K) = h_ML_out  
          GTGAT  (K) = T_sfc_n
C
          ZNGAT  (K) = LZSNGAT(K)
          TISLGAT(K) = LTICGAT(K)
          SNOGAT (K) = LZSNGAT(K)*RHONGAT(K)   ! but not consistent density!!!
C
C         * CALCULATE SNOW MELT AND RUNOFF.
C
          IF(LZSNGAT(K).LT.h_snow_in) then
            SMLTGAT(K) = ABS((LZSNGAT(K)-h_snow_in)*RHONGAT(K))
            ROFNGAT(K) = SMLTGAT(K)/DELT
          ELSE
            SMLTGAT(K) = 0.
            ROFNGAT(K) = 0.
          ENDIF
          IF(LZICGAT(K).LT.h_ice_in) then
            RICEMELT   = ABS((LZICGAT(K)-h_ice_in)*DENI/DELT)
          ELSE
            RICEMELT   = 0.
          ENDIF
C
C         * CALCULATE RUNOFF FOR FRESH WATER FLUX.
C
          IF(GCGAT(K).EQ.1.) THEN
            ROFGAT(K) = RPREGAT(K) + ROFNGAT(K) + RICEMELT
          ELSE
            ROFGAT(K) = PREGAT(K)-QFSGAT(K)
          ENDIF
C
C         * CALCULATE ADDITIONAL FIELDS.
C
          LUIMGAT (K) = LZICGAT(K)*DENI
          PCSNGAT (K) = SPCPGAT(K)*RHSIGAT(K)
        ENDDO
C
      RETURN
      END
