flake_interface.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     ###############################################################################
SUBROUTINE flake_interface (F, KI, &
! Atmospheric forcing
                               dMsnowdt_in, I_atm_in, Q_atm_lw_in, height_u_in, height_tq_in,     &
                               U_a_in, T_a_in, q_a_in, P_a_in,                                    &
! Constant parameters                           
                               del_time,            &
! Parameters that may change                             
                               albedo,              &
! Surface heat, momentum fluxes, and other diags                       
                               Q_sensible, Q_latent ,Q_momentum, z0t, Qsat, Ri, ustar, Cd_a,      &
                               Q_watvap, Q_latenti, Q_sublim, Q_atm_lw_up, pswe,                  &
! Switches to configure FLake runs
                               PPEW_A_COEF, PPEW_B_COEF, rho_a, HIMPLICIT_WIND                )  
!------------------------------------------------------------------------------
!
! Description:
!
!  The FLake interface is
!  a communication routine between "flake_driver"
!  and a prediction system that uses FLake.
!  It assigns the FLake variables at the previous time step 
!  to their input values given by the driving model,
!  calls a number of routines to compute the heat and radiation fluxes,
!  calls "flake_driver",
!  and returns the updated FLake variables to the driving model.
!  The "flake_interface" does not contain any Flake physics. 
!  It only serves as a convenient means to organize calls of "flake_driver"
!  and of external routines that compute heat and radiation fluxes.
!  The interface may (should) be changed so that to provide 
!  the most convenient use of FLake.
!  Within a 3D atmospheric prediction system,
!  "flake_driver" may be called in a DO loop within "flake_interface" 
!  for each grid-point where a lake is present.
!  In this way, the driving atmospheric model should call "flake_interface"
!  only once, passing the FLake variables to "flake_interface" as 2D fields. 
!
!  Lines embraced with "!_tmp" contain temporary parts of the code.
!  These should be removed prior to using FLake in applications.
!  Lines embraced/marked with "!_dev" may be replaced
!  as improved parameterizations are developed and tested.
!  Lines embraced/marked with "!_dm" are DM's comments
!  that may be helpful to a user.
!
!
! Current Code Owner: DWD, Dmitrii Mironov
!  Phone:  +49-69-8062 2705
!  Fax:    +49-69-8062 3721
!  E-mail: dmitrii.mironov@dwd.de
!
! History:
! Version    Date       Name
! ---------- ---------- ----
! 1.00       2005/11/17 Dmitrii Mironov 
!  Initial release
! !VERSION!  !DATE!     <Your name>
!  <Modification comments>
!
! Code Description:
! Language: Fortran 90.
! Software Standards: "European Standards for Writing and
! Documenting Exchangeable Fortran 90 Code".
!==============================================================================
!
! Declarations:
!
USE modd_flake_n, ONLY : flake_t
!
! Modules used:

!USE modd_data_parameters , ONLY : &
!      ireals,                   &! KIND-type parameter for real variables
!      iintegers                  ! KIND-type parameter for "normal" integer variables  

USE modd_flake_derivedtypes         ! Definitions of several derived TYPEs

USE modd_flake_parameters , ONLY :   &
    tpl_kappa_w                 ,  &! Molecular heat conductivity of water [J m^{-1} s^{-1} K^{-1}]
    tpl_t_f                     ,  &! Fresh water freezing point [K]
    tpl_rho_w_r                 ,  &! Maximum density of fresh water [kg m^{-3}]
    h_snow_min_flk              ,  &! Minimum snow thickness [m]
    h_ice_min_flk               ,  &! Minimum ice thickness [m]  
    h_skinlayer_flk                 ! Skin layer thickness [m]  

USE modd_flake_paramoptic_ref       ! Reference values of the optical characteristics
                               ! of the lake water, lake ice and snow 

USE mode_flake           , ONLY :    &
    flake_driver                ,  &! Subroutine, FLake driver
    flake_radflux               ,  &! Subroutine, computes radiation fluxes at various depths
    flake_snowdensity           ,  &! Function, computes snow density
                                    !
    t_snow_p_flk, t_snow_n_flk  ,  &! Temperature at the air-snow interface [K]
    t_ice_p_flk, t_ice_n_flk    ,  &! Temperature at the snow-ice or air-ice interface [K]
    t_mnw_p_flk, t_mnw_n_flk    ,  &! Mean temperature of the water column [K]
    t_wml_p_flk, t_wml_n_flk    ,  &! Mixed-layer temperature [K]
    t_bot_p_flk, t_bot_n_flk    ,  &! Temperature at the water-bottom sediment interface [K]
    t_b1_p_flk, t_b1_n_flk      ,  &! Temperature at the bottom of the upper layer of the sediments [K]
    c_t_p_flk, c_t_n_flk        ,  &! Shape factor (thermocline)
    h_snow_p_flk, h_snow_n_flk  ,  &! Snow thickness [m]
    h_ice_p_flk, h_ice_n_flk    ,  &! Ice thickness [m]
    h_ml_p_flk, h_ml_n_flk      ,  &! Thickness of the mixed-layer [m]
    h_b1_p_flk, h_b1_n_flk      ,  &! Thickness of the upper layer of bottom sediments [m]
                                  !
    q_snow_flk                  ,  &! Heat flux through the air-snow interface [W m^{-2}]
    q_ice_flk                   ,  &! Heat flux through the snow-ice or air-ice interface [W m^{-2}]
    q_w_flk                     ,  &! Heat flux through the ice-water or air-water interface [W m^{-2}]
    q_bot_flk                   ,  &! Heat flux through the water-bottom sediment interface [W m^{-2}]
    i_atm_flk                   ,  &! Radiation flux at the lower boundary of the atmosphere [W m^{-2}],
                                  ! i.e. the incident radiation flux with no regard for the surface albedo
    i_snow_flk                  ,  &! Radiation flux through the air-snow interface [W m^{-2}]
    i_ice_flk                   ,  &! Radiation flux through the snow-ice or air-ice interface [W m^{-2}]
    i_w_flk                     ,  &! Radiation flux through the ice-water or air-water interface [W m^{-2}]
    i_h_flk                     ,  &! Radiation flux through the mixed-layer-thermocline interface [W m^{-2}]
    i_bot_flk                   ,  &! Radiation flux through the water-bottom sediment interface [W m^{-2}]
    i_intm_0_h_flk              ,  &! Mean radiation flux over the mixed layer [W m^{-1}]
    i_intm_h_d_flk              ,  &! Mean radiation flux over the thermocline [W m^{-1}]
    q_star_flk                  ,  &! A generalized heat flux scale [W m^{-2}]
    u_star_w_flk                ,  &! Friction velocity in the surface layer of lake water [m s^{-1}]
    w_star_sfc_flk              ,  &! Convective velocity scale, using a generalized heat flux scale [m s^{-1}]
    dmsnowdt_flk                    ! The rate of snow accumulation [kg m^{-2} s^{-1}]  


USE mode_sfcflx          , ONLY :    &
    sfcflx_lwradwsfc            ,  &! Function, returns the surface long-wave radiation flux
    sfcflx_momsenlat            ,  &! Subroutine, computes fluxes of momentum and of sensible and latent heat
    z0u_sf                      ,  &! Roughness length with respect to wind velocity [m]
    z0t_sf                          ! Roughness length with respect to potential temperature [m]  


USE modd_flake_configure, ONLY : lflk_botsed_use 
!==============================================================================
!
USE modi_wind_threshold
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
!
IMPLICIT NONE

!==============================================================================
!
! Declarations

!
!*      0.1    declarations of arguments
!
!  Input (procedure arguments)
!
TYPE(flake_t), INTENT(INOUT) :: F
!
INTEGER, INTENT(IN)  :: KI              ! number of points
!
INTEGER :: i ! DO loop index
!
REAL, DIMENSION(KI), INTENT(IN) ::   &
    dMsnowdt_in                       ,  &! The rate of snow accumulation [kg m^{-2} s^{-1}]
    I_atm_in                          ,  &! Solar radiation flux at the surface [W m^{-2}]
    Q_atm_lw_in                       ,  &! Long-wave radiation flux from the atmosphere [W m^{-2}]
    height_u_in                       ,  &! Height above the lake surface where the wind speed is measured [m]
    height_tq_in                      ,  &! Height where temperature and humidity are measured [m]
    U_a_in                            ,  &! Wind speed at z=height_u_in [m s^{-1}]
    T_a_in                            ,  &! Air temperature at z=height_tq_in [K]
    q_a_in                            ,  &! Air specific humidity at z=height_tq_in
    P_a_in                                ! Surface air pressure [N m^{-2} = kg m^{-1} s^{-2}]  

REAL, DIMENSION(KI), INTENT(IN) ::   &
    del_time                          ,  &! The model time step [s]
    PPEW_A_COEF                       ,  &! coefficient A (m^2 s kg^{-1}) and B (m s^{-1})
    PPEW_B_COEF                       ,  &! for wind implicitation :  V+ = - A * rho_a u'w' + B
    rho_a                                 ! Air density (kg m ^{-3}) (from forcing atm. data)  
!
 CHARACTER(LEN=*),     INTENT(IN)  :: HIMPLICIT_WIND   ! wind implicitation option
!                                                     ! 'OLD' = direct
!                                                     ! 'NEW' = Taylor serie, order 1
!
!
!  Input/Output (procedure arguments)

REAL, DIMENSION(KI), INTENT(IN)    :: albedo  ! surface albedo with respect to the solar radiation
                                              ! (free water, ice or snow; e.g. update_rade_flake.f90)
 
!  Output (procedure arguments)

REAL, DIMENSION(KI), INTENT(INOUT)  ::    &
    Q_sensible             ,  &! Sensible heat flux [W m^{-2}]
    Q_latent               ,  &! Total Latent heat flux [W m^{-2}]
    Q_watvap               ,  &! Total Flux of water vapour [kg m^{-2} s^{-1}] 
    Q_latenti              ,  &! Sublimation Latent heat flux [W m^{-2}]
    Q_sublim               ,  &! Flux of sublimation [kg m^{-2} s^{-1}]    
    Q_momentum             ,  &! Momentum flux [N m^{-2}]
    z0t                    ,  &! Roughness length with respect to potential temperature [m]
    Ri                     ,  &! Gradient Richardson number 
    ustar                  ,  &! air friction velocity  
    Cd_a                       ! wind drag coefficient [no unit]
!
REAL, DIMENSION(KI), INTENT(OUT)  ::    &
    Qsat                   ,  &! specific humidity at saturation [kg.kg-1]
    Q_atm_lw_up            ,  &! Upward longwave flux at t [W m^{-2}]
    pswe                       ! snow water equivalent [kg.m-2]
!
!*      0.2    declarations of local variables
!
INTEGER  :: JL                  ! loop counter on horizontal points

REAL :: T_sfc_n ! Surface temperature at the new time step [K]
REAL :: ustar2  ! square of air friction velocity (m2/s2)  
REAL :: zvmod   ! wind at t+1   
REAL, DIMENSION(KI) ::    &
    zwind             ! thresholded wind

type(opticpar_medium), DIMENSION(KI) ::  &
    opticpar_water                       ,  &! Optical characteristics of water
    opticpar_ice                         ,  &! Optical characteristics of ice
    opticpar_snow                            ! Optical characteristics of snow      
REAL(KIND=JPRB) :: ZHOOK_HANDLE
!
REAL, PARAMETER  :: ZTIMEMAX      = 300.  ! s  Maximum timescale without time spliting
!
INTEGER  :: JDT, INDT
REAL :: zaux, zcond, zloc, ZTSTEP
!
!==============================================================================
!  Start calculations
!------------------------------------------------------------------------------
IF (lhook) CALL dr_hook('FLAKE_INTERFACE',0,zhook_handle)
!
!------------------------------------------------------------------------------
!  Set optical characteristics of the lake water, lake ice and snow
!------------------------------------------------------------------------------
!
! Use default values
!opticpar_water = opticpar_water_ref  ! don't use default values

opticpar_ice   = opticpar_ice_opaque   ! Opaque ice
opticpar_snow  = opticpar_snow_opaque  ! Opaque snow
!
lflk_botsed_use = f%LSEDIMENTS
!
zwind = wind_threshold(u_a_in,height_u_in)
!
h_point_loop: DO jl = 1,ki ! begin of loop on horizontal points
!------------------------------------------------------------------------------
!  Set initial values
!------------------------------------------------------------------------------
   
opticpar_water(jl) = opticpar_medium(1,                       &
      (/1., (0.,i=2,nband_optic_max)/),            &
      (/f%XEXTCOEF_WATER(jl), (1.e+10,i=2,nband_optic_max)/))  
   t_snow_p_flk = f%XT_SNOW(jl)
   t_ice_p_flk  = f%XT_ICE(jl)
   t_mnw_p_flk  = f%XT_MNW(jl)
   t_wml_p_flk  = f%XT_WML(jl)
   t_bot_p_flk  = f%XT_BOT(jl)
   t_b1_p_flk   = f%XT_B1(jl)
   c_t_p_flk    = f%XCT(jl)
   h_snow_p_flk = f%XH_SNOW(jl)
   h_ice_p_flk  = f%XH_ICE(jl)
   h_ml_p_flk   = f%XH_ML(jl)
   h_b1_p_flk   = f%XH_B1(jl)
   
!------------------------------------------------------------------------------
!  Set the rate of snow accumulation
!------------------------------------------------------------------------------
   
   dmsnowdt_flk = dmsnowdt_in(jl)  
   
!------------------------------------------------------------------------------
!  Compute solar radiation fluxes (positive downward)
!------------------------------------------------------------------------------
   
   i_atm_flk = i_atm_in(jl)
   CALL flake_radflux ( f%XWATER_DEPTH(jl), albedo(jl), opticpar_water(jl),  &
                          opticpar_ice(jl), opticpar_snow(jl) )  
   
!------------------------------------------------------------------------------
!  Compute long-wave radiation fluxes (positive downward)
!  lwd-lwu = emis*(lwd-sigma*ts**4)
!------------------------------------------------------------------------------
!   
   q_w_flk = f%XEMIS(jl)*q_atm_lw_in(jl) - sfcflx_lwradwsfc(f%XEMIS(jl),f%XTS(jl))
!
   q_atm_lw_up(jl) = q_atm_lw_in(jl) - q_w_flk
!   
!------------------------------------------------------------------------------
!  Compute the surface friction velocity and fluxes of sensible and latent heat 
!------------------------------------------------------------------------------
   
   IF (f%CFLK_FLUX=='FLAKE') THEN
      !
      q_momentum(jl) = - rho_a(jl) * ustar(jl)**2
      !
      CALL sfcflx_momsenlat ( height_u_in(jl), height_tq_in(jl), f%XWATER_FETCH(jl),  &
                             u_a_in(jl), t_a_in(jl), q_a_in(jl), f%XTS(jl), &
                             p_a_in(jl), h_ice_p_flk, q_momentum(jl),       &
                             q_sensible(jl), q_latent(jl), q_watvap(jl),    &
                             ri(jl), f%XZ0(jl), z0t(jl), qsat(jl),             &
                             q_latenti(jl), q_sublim(jl)                    )
      f%XZ0(jl)= z0u_sf
      z0t(jl)=z0t_sf
      ! recomputes the future wind speed and associated momentum flux
      ! in the case wind speed is implicited  (V. Masson, Meteo-France)
      ! 1st step : drag coefficient
      ! It is retrieved assumed a relationship between momentum flux
      !  and previous time-step wind : Q_mom = - rho_a * Cd_a * U_a_in**2
      !
      cd_a(jl) = - q_momentum(jl) / rho_a(jl) / zwind(jl)**2 
      ! 2nd step : friction velocity (for air) computed with future wind speed
      !            (the latter computed using implicit coefficients)
      ustar2 = 0.0
      zvmod  = u_a_in(jl)
      IF(himplicit_wind=='OLD')THEN
      !  old implicitation
         ustar2 = (cd_a(jl)*u_a_in(jl)*ppew_b_coef(jl))/            &
                  (1.0-rho_a(jl)*cd_a(jl)*u_a_in(jl)*ppew_a_coef(jl))
      ELSE
      !  new implicitation
         ustar2 = (cd_a(jl)*u_a_in(jl)*(2.*ppew_b_coef(jl)-u_a_in(jl)))/ &
                  (1.0-2.0*rho_a(jl)*cd_a(jl)*u_a_in(jl)*ppew_a_coef(jl))
         zvmod  = rho_a(jl)*ppew_a_coef(jl)*ustar2 + ppew_b_coef(jl)
         zvmod  = max(0.0,zvmod)
         IF(ppew_a_coef(jl)/= 0.)THEN
           ustar2 = max((zvmod-ppew_b_coef(jl))/(rho_a(jl)*ppew_a_coef(jl)),0.0)
         ENDIF
      ENDIF
      ustar(jl) =sqrt(ustar2)
      ! 3rd step : momentum flux computed with the future wind speed
      q_momentum(jl) = - rho_a(jl) * ustar2

   END IF
   u_star_w_flk = sqrt(-q_momentum(jl)/tpl_rho_w_r)
   
!------------------------------------------------------------------------------
!  Compute heat fluxes Q_snow_flk, Q_ice_flk, Q_w_flk
!------------------------------------------------------------------------------
   
   q_w_flk = q_w_flk - q_sensible(jl) - q_latent(jl) ! Add sensible and latent heat fluxes 
                                                     ! (notice the signs)
   IF(h_ice_p_flk.GE.h_ice_min_flk) THEN         ! Ice exists
     IF(h_snow_p_flk.GE.h_snow_min_flk) THEN     ! There is snow above the ice
       q_snow_flk = q_w_flk
       q_ice_flk  = 0.
       q_w_flk    = 0.
     ELSE                                        ! No snow above the ice
       q_snow_flk = 0.
       q_ice_flk  = q_w_flk
       q_w_flk    = 0.
     END IF
   ELSE                                          ! No ice-snow cover
       q_snow_flk = 0.
       q_ice_flk  = 0.
   END IF

!------------------------------------------------------------------------------
!  Advance FLake variables
!  Time splitting parameter for *very large time steps* 
! -----------------------------------------------------------------------------
!
   indt    = max(1,nint(del_time(jl)/ztimemax))
   ztstep  = del_time(jl)/REAL(indt)

   DO jdt=1,indt
!
   t_snow_p_flk = f%XT_SNOW(jl)
   t_ice_p_flk  = f%XT_ICE(jl)
   t_mnw_p_flk  = f%XT_MNW(jl)
   t_wml_p_flk  = f%XT_WML(jl)
   t_bot_p_flk  = f%XT_BOT(jl)
   t_b1_p_flk   = f%XT_B1(jl)
   c_t_p_flk    = f%XCT(jl)
   h_snow_p_flk = f%XH_SNOW(jl)
   h_ice_p_flk  = f%XH_ICE(jl)
   h_ml_p_flk   = f%XH_ML(jl)
   h_b1_p_flk   = f%XH_B1(jl)
!
   CALL flake_driver ( f%XWATER_DEPTH(jl), f%XDEPTH_BS(jl), f%XT_BS(jl), f%XCORIO(jl), &
                       opticpar_water(jl)%extincoef_optic(1),                 &
                       ztstep, f%XTS(jl), t_sfc_n ) 
!                       
!------------------------------------------------------------------------------
!  Set output values
!------------------------------------------------------------------------------
!   
   f%XT_SNOW(jl) = t_snow_n_flk  
   f%XT_ICE(jl)  = t_ice_n_flk      
   f%XT_MNW(jl)  = t_mnw_n_flk     
   f%XT_WML(jl)  = t_wml_n_flk    
   f%XT_BOT(jl)  = t_bot_n_flk   
   f%XT_B1(jl)   = t_b1_n_flk    
   f%XCT(jl)    =  c_t_n_flk     
   f%XH_SNOW(jl) = h_snow_n_flk   
   f%XH_ICE(jl)  = h_ice_n_flk    
   f%XH_ML(jl)   = h_ml_n_flk     
   f%XH_B1(jl)   = h_b1_n_flk  
   f%XTS(jl)  = t_sfc_n
!
   ENDDO                      
!
   pswe(jl) = f%XH_SNOW(jl)*flake_snowdensity(f%XH_SNOW(jl))
!
!------------------------------------------------------------------------------
!  Compute skin temperature in case of no ice nor snow
!------------------------------------------------------------------------------
!
!  Q_w_flk=LWD-LWU-(QH+QE)      accounts for water phase (water, ice, snow)
!  (1-albedo)*I_atm_flk=SWD-SWU accounts for water phase (water, ice, snow)
!
   IF (f%LSKINTEMP.AND.(f%XH_ICE(jl)<h_ice_min_flk).AND.(f%CFLK_FLUX=='FLAKE')) THEN

         zaux  = (1.0-exp(-opticpar_water(jl)%extincoef_optic(1) * h_skinlayer_flk)) &
                        / opticpar_water(jl)%extincoef_optic(1)

         zcond = tpl_kappa_w

         zloc = ((1.0-albedo(jl))*i_atm_flk+q_w_flk)*h_skinlayer_flk-(1.0-albedo(jl))*i_atm_flk*zaux

         f%XTS(jl)  = t_sfc_n + zloc / zcond

   ENDIF
!   
ENDDO h_point_loop
!
IF (lhook) CALL dr_hook('FLAKE_INTERFACE',1,zhook_handle)
!
!------------------------------------------------------------------------------
!  End calculations
!==============================================================================

END SUBROUTINE flake_interface