New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 9019 for branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 – NEMO

Ignore:
Timestamp:
2017-12-13T15:58:53+01:00 (6 years ago)
Author:
timgraham
Message:

Merge of dev_CNRS_2017 into branch

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8666 r9019  
    1515   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk 
    1616   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore 
    17    !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
     17   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)  
    1920   !!---------------------------------------------------------------------- 
    2021 
    … …  
    2324   !!   sbc_blk       : bulk formulation as ocean surface boundary condition 
    2425   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean 
    25    !!   blk_ice       : computes momentum, heat and freshwater fluxes over sea ice 
    2626   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP 
    2727   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air) 
    2828   !!   q_sat         : saturation humidity as a function of SLP and temperature 
    2929   !!   L_vap         : latent heat of vaporization of water as a function of temperature 
     30   !!             sea-ice case only :  
     31   !!   blk_ice_tau   : provide the air-ice stress 
     32   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface 
     33   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     34   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     35   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
    3036   !!---------------------------------------------------------------------- 
    3137   USE oce            ! ocean dynamics and tracers 
    … …  
    4046   USE lib_fortran    ! to use key_nosignedzero 
    4147#if defined key_lim3 
    42    USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b 
    43    USE limthd_dh      ! for CALL lim_thd_snwblow 
    44 #elif defined key_lim2 
    45    USE ice_2   , ONLY :   u_ice, v_ice 
    46    USE par_ice_2      ! LIM-2 parameters 
     48   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s 
     49   USE icethd_dh      ! for CALL ice_thd_snwblow 
    4750#endif 
    4851   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    … …  
    5457   USE in_out_manager ! I/O manager 
    5558   USE lib_mpp        ! distribued memory computing library 
    56    USE wrk_nemo       ! work arrays 
    5759   USE timing         ! Timing 
    5860   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    … …  
    6466   PUBLIC   sbc_blk_init  ! called in sbcmod 
    6567   PUBLIC   sbc_blk       ! called in sbcmod 
    66 #if defined key_lim2 || defined key_lim3 
    67    PUBLIC   blk_ice_tau   ! routine called in sbc_ice_lim module 
    68    PUBLIC   blk_ice_flx   ! routine called in sbc_ice_lim module 
    69 #endif 
     68#if defined key_lim3 
     69   PUBLIC   blk_ice_tau   ! routine called in iceforcing 
     70   PUBLIC   blk_ice_flx   ! routine called in iceforcing 
     71   PUBLIC   blk_ice_qcn   ! routine called in iceforcing 
     72#endif  
    7073 
    7174!!Lolo: should ultimately be moved in the module with all physical constants ? 
    … …  
    9699   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation 
    97100   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant 
    98    REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! iovi 1.63e-3     ! transfer coefficient over ice 
     101   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice 
    99102   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant 
    100103   ! 
    … …  
    107110   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data 
    108111   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation 
    109    REAL(wp) ::   rn_efac        ! multiplication factor for evaporation (clem) 
    110    REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress (clem) 
     112   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation 
     113   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
    111114   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements 
    112115   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements 
    113    LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm and oce-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     116!!gm ref namelist initialize it so remove the setting to false below 
     117   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012) 
     118   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015) 
    114119   ! 
    115    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_oce   ! air-ocean drag (clem) 
     120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau) 
     121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens) 
     122   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat) 
     123   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme) 
     124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme) 
     125   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme 
    116126 
    117127   INTEGER  ::   nblk           ! choice of the bulk algorithm 
    … …  
    135145      !!             ***  ROUTINE sbc_blk_alloc *** 
    136146      !!------------------------------------------------------------------- 
    137       ALLOCATE( Cd_oce(jpi,jpj) , STAT=sbc_blk_alloc ) 
     147      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), & 
     148         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc ) 
    138149      ! 
    139150      IF( lk_mpp             )   CALL mpp_sum ( sbc_blk_alloc ) 
    … …  
    141152   END FUNCTION sbc_blk_alloc 
    142153 
     154 
    143155   SUBROUTINE sbc_blk_init 
    144156      !!--------------------------------------------------------------------- 
    … …  
    148160      !! 
    149161      !! ** Method  :  
    150       !! 
    151       !!      C A U T I O N : never mask the surface stress fields 
    152       !!                      the stress is assumed to be in the (i,j) mesh referential 
    153       !! 
    154       !! ** Action  :    
    155162      !! 
    156163      !!---------------------------------------------------------------------- 
    … …  
    167174         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm 
    168175         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &  
    169          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12 
     176         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15 
    170177      !!--------------------------------------------------------------------- 
    171178      ! 
    … …  
    262269         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    263270         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
     271         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
     272         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
    264273         ! 
    265274         WRITE(numout,*) 
    … …  
    281290      !! 
    282291      !! ** Purpose :   provide at each time step the surface ocean fluxes 
    283       !!      (momentum, heat, freshwater and runoff) 
     292      !!              (momentum, heat, freshwater and runoff) 
    284293      !! 
    285294      !! ** Method  : (1) READ each fluxes in NetCDF files: 
    … …  
    364373      INTEGER  ::   ji, jj               ! dummy loop indices 
    365374      REAL(wp) ::   zztmp                ! local variable 
    366       REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
    367       REAL(wp), DIMENSION(:,:), POINTER ::   zsq               ! specific humidity at pst 
    368       REAL(wp), DIMENSION(:,:), POINTER ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
    369       REAL(wp), DIMENSION(:,:), POINTER ::   zqla, zevap       ! latent heat fluxes and evaporation 
    370       REAL(wp), DIMENSION(:,:), POINTER ::   Cd                ! transfer coefficient for momentum      (tau) 
    371       REAL(wp), DIMENSION(:,:), POINTER ::   Ch                ! transfer coefficient for sensible heat (Q_sens) 
    372       REAL(wp), DIMENSION(:,:), POINTER ::   Ce                ! tansfert coefficient for evaporation   (Q_lat) 
    373       REAL(wp), DIMENSION(:,:), POINTER ::   zst               ! surface temperature in Kelvin 
    374       REAL(wp), DIMENSION(:,:), POINTER ::   zt_zu             ! air temperature at wind speed height 
    375       REAL(wp), DIMENSION(:,:), POINTER ::   zq_zu             ! air spec. hum.  at wind speed height 
    376       REAL(wp), DIMENSION(:,:), POINTER ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    377       REAL(wp), DIMENSION(:,:), POINTER ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    378       REAL(wp), DIMENSION(:,:), POINTER ::   zrhoa             ! density of air   [kg/m^3] 
    379       !!--------------------------------------------------------------------- 
    380       ! 
    381       IF( nn_timing == 1 )  CALL timing_start('blk_oce') 
    382       ! 
    383       CALL wrk_alloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 
    384       CALL wrk_alloc( jpi,jpj,   Cd, Ch, Ce, zst, zt_zu, zq_zu ) 
    385       CALL wrk_alloc( jpi,jpj,   zU_zu, ztpot, zrhoa ) 
    386       ! 
    387  
     375      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     376      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst 
     377      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes 
     378      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation 
     379      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin 
     380      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
     381      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
     382      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3] 
     383      !!--------------------------------------------------------------------- 
     384      ! 
     385      IF( ln_timing )   CALL timing_start('blk_oce') 
     386      ! 
    388387      ! local scalars ( place there for vector optimisation purposes) 
    389388      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K) 
    … …  
    430429      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave 
    431430 
    432  
    433  
    434431      ! ----------------------------------------------------------------------------- ! 
    435432      !     II    Turbulent FLUXES                                                    ! 
    … …  
    447444      ! 
    448445      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2 
    449          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 
     446         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    450447      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0 
    451          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 
     448         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    452449      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5 
    453          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 
     450         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    454451      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF 
    455          &                                               Cd, Ch, Ce, zt_zu, zq_zu, zU_zu ) 
     452         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce ) 
    456453      CASE DEFAULT 
    457454         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' ) 
    … …  
    460457      !                          ! Compute true air density : 
    461458      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...) 
    462          zrhoa(:,:) = rho_air( zt_zu(:,:)             , zq_zu(:,:)             , sf(jp_slp)%fnow(:,:,1) ) 
     459         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) ) 
    463460      ELSE                                      ! At zt: 
    464461         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    465462      END IF 
    466463 
    467       Cd_oce(:,:) = Cd(:,:)  ! record value of pure ocean-atm. drag (clem) 
     464!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     465!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef. 
    468466 
    469467      DO jj = 1, jpj             ! tau module, i and j component 
    470468         DO ji = 1, jpi 
    471             zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd(ji,jj)   ! using bulk wind speed 
     469            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed 
    472470            taum  (ji,jj) = zztmp * wndm  (ji,jj) 
    473471            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
    … …  
    504502      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN 
    505503         !! q_air and t_air are given at 10m (wind reference height) 
    506          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
    507          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
     504         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed 
     505         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed 
    508506      ELSE 
    509507         !! q_air and t_air are not given at 10m (wind reference height) 
    510508         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!! 
    511          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce(:,:)*(zsq(:,:) - zq_zu(:,:) ) ) ! Evaporation ! using bulk wind speed 
    512          zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch(:,:)*(zst(:,:) - zt_zu(:,:) )   ! Sensible Heat ! using bulk wind speed 
     509         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed 
     510         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed 
    513511      ENDIF 
    514512 
    … …  
    517515 
    518516      IF(ln_ctl) THEN 
    519          CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce , clinfo2=' Ce  : ' ) 
    520          CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch , clinfo2=' Ch  : ' ) 
     517         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' ) 
     518         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' ) 
    521519         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' ) 
    522520         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' ) 
    … …  
    557555         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac   ! output total precipitation [kg/m2/s] 
    558556         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac   ! output solid precipitation [kg/m2/s] 
    559          CALL iom_put( 'snowpre', sprecip * 86400. )        ! Snow 
    560          CALL iom_put( 'precip' , tprecip * 86400. )        ! Total precipitation 
     557         CALL iom_put( 'snowpre', sprecip )                 ! Snow 
     558         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    561559      ENDIF 
    562560      ! 
    … …  
    569567      ENDIF 
    570568      ! 
    571       CALL wrk_dealloc( jpi,jpj,   zwnd_i, zwnd_j, zsq, zqlw, zqsb, zqla, zevap ) 
    572       CALL wrk_dealloc( jpi,jpj,   Cd, Ch, Ce, zst, zt_zu, zq_zu ) 
    573       CALL wrk_dealloc( jpi,jpj,   zU_zu, ztpot, zrhoa ) 
    574       ! 
    575       IF( nn_timing == 1 )  CALL timing_stop('blk_oce') 
     569      IF( ln_timing )   CALL timing_stop('blk_oce') 
    576570      ! 
    577571   END SUBROUTINE blk_oce 
    578572 
    579 #if defined key_lim2 || defined key_lim3 
     573 
     574 
     575   FUNCTION rho_air( ptak, pqa, pslp ) 
     576      !!------------------------------------------------------------------------------- 
     577      !!                           ***  FUNCTION rho_air  *** 
     578      !! 
     579      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
     580      !! 
     581      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
     582      !!------------------------------------------------------------------------------- 
     583      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
     584      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
     585      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
     586      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
     587      !!------------------------------------------------------------------------------- 
     588      ! 
     589      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
     590      ! 
     591   END FUNCTION rho_air 
     592 
     593 
     594   FUNCTION cp_air( pqa ) 
     595      !!------------------------------------------------------------------------------- 
     596      !!                           ***  FUNCTION cp_air  *** 
     597      !! 
     598      !! ** Purpose : Compute specific heat (Cp) of moist air 
     599      !! 
     600      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     601      !!------------------------------------------------------------------------------- 
     602      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
     603      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
     604      !!------------------------------------------------------------------------------- 
     605      ! 
     606      Cp_air = Cp_dry + Cp_vap * pqa 
     607      ! 
     608   END FUNCTION cp_air 
     609 
     610 
     611   FUNCTION q_sat( ptak, pslp ) 
     612      !!---------------------------------------------------------------------------------- 
     613      !!                           ***  FUNCTION q_sat  *** 
     614      !! 
     615      !! ** Purpose : Specific humidity at saturation in [kg/kg] 
     616      !!              Based on accurate estimate of "e_sat" 
     617      !!              aka saturation water vapor (Goff, 1957) 
     618      !! 
     619      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     620      !!---------------------------------------------------------------------------------- 
     621      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
     622      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
     623      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
     624      ! 
     625      INTEGER  ::   ji, jj         ! dummy loop indices 
     626      REAL(wp) ::   ze_sat, ztmp   ! local scalar 
     627      !!---------------------------------------------------------------------------------- 
     628      ! 
     629      DO jj = 1, jpj 
     630         DO ji = 1, jpi 
     631            ! 
     632            ztmp = rt0 / ptak(ji,jj) 
     633            ! 
     634            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
     635            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
     636               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
     637               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     638               ! 
     639            q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
     640            ! 
     641         END DO 
     642      END DO 
     643      ! 
     644   END FUNCTION q_sat 
     645 
     646 
     647   FUNCTION gamma_moist( ptak, pqa ) 
     648      !!---------------------------------------------------------------------------------- 
     649      !!                           ***  FUNCTION gamma_moist  *** 
     650      !! 
     651      !! ** Purpose : Compute the moist adiabatic lapse-rate. 
     652      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
     653      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
     654      !! 
     655      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     656      !!---------------------------------------------------------------------------------- 
     657      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
     658      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
     659      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
     660      ! 
     661      INTEGER  ::   ji, jj         ! dummy loop indices 
     662      REAL(wp) :: zrv, ziRT        ! local scalar 
     663      !!---------------------------------------------------------------------------------- 
     664      ! 
     665      DO jj = 1, jpj 
     666         DO ji = 1, jpi 
     667            zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
     668            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
     669            gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
     670         END DO 
     671      END DO 
     672      ! 
     673   END FUNCTION gamma_moist 
     674 
     675 
     676   FUNCTION L_vap( psst ) 
     677      !!--------------------------------------------------------------------------------- 
     678      !!                           ***  FUNCTION L_vap  *** 
     679      !! 
     680      !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
     681      !! 
     682      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
     683      !!---------------------------------------------------------------------------------- 
     684      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
     685      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
     686      !!---------------------------------------------------------------------------------- 
     687      ! 
     688      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
     689      ! 
     690   END FUNCTION L_vap 
     691 
     692#if defined key_lim3 
     693   !!---------------------------------------------------------------------- 
     694   !!   'key_lim3'                                       ESIM sea-ice model 
     695   !!---------------------------------------------------------------------- 
     696   !!   blk_ice_tau : provide the air-ice stress 
     697   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface 
     698   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler) 
     699   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag 
     700   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag  
     701   !!---------------------------------------------------------------------- 
    580702 
    581703   SUBROUTINE blk_ice_tau 
    … …  
    590712      !!--------------------------------------------------------------------- 
    591713      INTEGER  ::   ji, jj    ! dummy loop indices 
    592       ! 
    593       REAL(wp), DIMENSION(:,:)  , POINTER :: zrhoa 
    594       ! 
    595       REAL(wp) ::   zwnorm_f, zwndi_f , zwndj_f               ! relative wind module and components at F-point 
    596       REAL(wp) ::             zwndi_t , zwndj_t               ! relative wind components at T-point 
    597       REAL(wp), DIMENSION(:,:), POINTER ::   Cd               ! transfer coefficient for momentum      (tau) 
    598       !!--------------------------------------------------------------------- 
    599       ! 
    600       IF( nn_timing == 1 )  CALL timing_start('blk_ice_tau') 
    601       ! 
    602       CALL wrk_alloc( jpi,jpj, zrhoa ) 
    603       CALL wrk_alloc( jpi,jpj, Cd ) 
    604  
    605       Cd(:,:) = Cd_ice 
    606  
    607       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al. 2012) (clem) 
    608 #if defined key_lim3 
    609       IF( ln_Cd_L12 ) THEN 
    610          CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
     714      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point 
     715      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
     716      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau) 
     717      !!--------------------------------------------------------------------- 
     718      ! 
     719      IF( ln_timing )   CALL timing_start('blk_ice_tau') 
     720      ! 
     721      ! set transfer coefficients to default sea-ice values 
     722      Cd_atm(:,:) = Cd_ice 
     723      Ch_atm(:,:) = Cd_ice 
     724      Ce_atm(:,:) = Cd_ice 
     725 
     726      wndm_ice(:,:) = 0._wp      !!gm brutal.... 
     727 
     728      ! ------------------------------------------------------------ ! 
     729      !    Wind module relative to the moving ice ( U10m - U_ice )   ! 
     730      ! ------------------------------------------------------------ ! 
     731      SELECT CASE( cp_ice_msh ) 
     732      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
     733         !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
     734         DO jj = 2, jpjm1 
     735            DO ji = 2, jpim1   ! B grid : NO vector opt 
     736               ! ... scalar wind at T-point (fld being at T-point) 
     737               zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
     738                  &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
     739               zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
     740                  &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
     741               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     742            END DO 
     743         END DO 
     744         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
     745         ! 
     746      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
     747         DO jj = 2, jpjm1 
     748            DO ji = fs_2, fs_jpim1   ! vect. opt. 
     749               zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
     750               zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
     751               wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     752            END DO 
     753         END DO 
     754         CALL lbc_lnk( wndm_ice, 'T',  1. ) 
     755         ! 
     756      END SELECT 
     757 
     758      ! Make ice-atm. drag dependent on ice concentration 
     759      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations 
     760         CALL Cdn10_Lupkes2012( Cd_atm ) 
     761         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical 
     762      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations 
     763         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )  
    611764      ENDIF 
    612 #endif 
     765 
     766!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef. 
     767!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef. 
    613768 
    614769      ! local scalars ( place there for vector optimisation purposes) 
    615770      ! Computing density of air! Way denser that 1.2 over sea-ice !!! 
    616       !! 
    617771      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1)) 
    618772 
    … …  
    620774      utau_ice  (:,:) = 0._wp 
    621775      vtau_ice  (:,:) = 0._wp 
    622       wndm_ice  (:,:) = 0._wp 
    623776      !!gm end 
    624777 
    625       ! ----------------------------------------------------------------------------- ! 
    626       !    Wind components and module relative to the moving ocean ( U10m - U_ice )   ! 
    627       ! ----------------------------------------------------------------------------- ! 
     778      ! ------------------------------------------------------------ ! 
     779      !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
     780      ! ------------------------------------------------------------ ! 
    628781      SELECT CASE( cp_ice_msh ) 
    629782      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation) 
    630          !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked) 
    631783         DO jj = 2, jpjm1 
    632784            DO ji = 2, jpim1   ! B grid : NO vector opt 
    … …  
    636788               zwndj_f = 0.25 * (  sf(jp_wndj)%fnow(ji-1,jj  ,1) + sf(jp_wndj)%fnow(ji  ,jj  ,1)   & 
    637789                  &              + sf(jp_wndj)%fnow(ji-1,jj-1,1) + sf(jp_wndj)%fnow(ji  ,jj-1,1)  ) - rn_vfac * v_ice(ji,jj) 
    638                zwnorm_f = zrhoa(ji,jj) * Cd(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    639790               ! ... ice stress at I-point 
     791               zwnorm_f = zrhoa(ji,jj) * Cd_atm(ji,jj) * SQRT( zwndi_f * zwndi_f + zwndj_f * zwndj_f ) 
    640792               utau_ice(ji,jj) = zwnorm_f * zwndi_f 
    641793               vtau_ice(ji,jj) = zwnorm_f * zwndj_f 
    642                ! ... scalar wind at T-point (fld being at T-point) 
    643                zwndi_t = sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  u_ice(ji,jj+1) + u_ice(ji+1,jj+1)   & 
    644                   &                                                    + u_ice(ji,jj  ) + u_ice(ji+1,jj  )  ) 
    645                zwndj_t = sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.25 * (  v_ice(ji,jj+1) + v_ice(ji+1,jj+1)   & 
    646                   &                                                    + v_ice(ji,jj  ) + v_ice(ji+1,jj  )  ) 
    647                wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    648794            END DO 
    649795         END DO 
    650796         CALL lbc_lnk( utau_ice, 'I', -1. ) 
    651797         CALL lbc_lnk( vtau_ice, 'I', -1. ) 
    652          CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    653798         ! 
    654799      CASE( 'C' )                  ! C-grid ice dynamics :   U & V-points (same as ocean) 
    655          DO jj = 2, jpj 
    656             DO ji = fs_2, jpi   ! vect. opt. 
    657                zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  ) 
    658                zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  ) 
    659                wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
    660             END DO 
    661          END DO 
    662800         DO jj = 2, jpjm1 
    663801            DO ji = fs_2, fs_jpim1   ! vect. opt. 
    664                utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )                          & 
     802               utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            & 
    665803                  &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) ) 
    666                vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )                          & 
     804               vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            & 
    667805                  &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) ) 
    668806            END DO 
    … …  
    670808         CALL lbc_lnk( utau_ice, 'U', -1. ) 
    671809         CALL lbc_lnk( vtau_ice, 'V', -1. ) 
    672          CALL lbc_lnk( wndm_ice, 'T',  1. ) 
    673810         ! 
    674811      END SELECT 
    675  
     812      ! 
    676813      IF(ln_ctl) THEN 
    677814         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ') 
    678815         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ') 
    679816      ENDIF 
    680  
    681       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_tau') 
    682  
     817      ! 
     818      IF( ln_timing )   CALL timing_stop('blk_ice_tau') 
     819      ! 
    683820   END SUBROUTINE blk_ice_tau 
    684821 
    685822 
    686    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     823   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    687824      !!--------------------------------------------------------------------- 
    688825      !!                     ***  ROUTINE blk_ice_flx  *** 
    689826      !! 
    690       !! ** Purpose :   provide the surface boundary condition over sea-ice 
     827      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface 
    691828      !! 
    692829      !! ** Method  :   compute heat and freshwater exchanged 
    … …  
    697834      !!--------------------------------------------------------------------- 
    698835      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     836      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     837      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    699838      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    700839      !! 
    … …  
    703842      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    704843      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
    705       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice 
    706       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice 
    707       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqlw        ! long wave heat sensitivity over ice 
    708       REAL(wp), DIMENSION(:,:,:), POINTER ::   z_dqsb        ! sensible  heat sensitivity over ice 
    709       REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    710       REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
    711       REAL(wp), DIMENSION(:,:)  , POINTER ::   Cd            ! transfer coefficient for momentum      (tau) 
    712       !!--------------------------------------------------------------------- 
    713       ! 
    714       IF( nn_timing == 1 )  CALL timing_start('blk_ice_flx') 
    715       ! 
    716       CALL wrk_alloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    717       CALL wrk_alloc( jpi,jpj,       zrhoa) 
    718       CALL wrk_alloc( jpi,jpj, Cd ) 
    719  
    720       Cd(:,:) = Cd_ice 
    721  
    722       ! Make ice-atm. drag dependent on ice concentration (see Lupkes et al.  2012) (clem) 
    723 #if defined key_lim3 
    724       IF( ln_Cd_L12 ) THEN 
    725          CALL Cdn10_Lupkes2012( Cd ) ! calculate new drag from Lupkes(2012) equations 
    726       ENDIF 
    727 #endif 
    728  
    729       ! 
    730       ! local scalars ( place there for vector optimisation purposes) 
    731       zcoef_dqlw   = 4.0 * 0.95 * Stef 
    732       zcoef_dqla   = -Ls * 11637800. * (-5897.8) 
     844      REAL(wp) ::   zfr1, zfr2               ! local variables 
     845      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
     846      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice 
     847      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice 
     848      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice 
     849      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
     850      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa 
     851      !!--------------------------------------------------------------------- 
     852      ! 
     853      IF( ln_timing )   CALL timing_start('blk_ice_flx') 
     854      ! 
     855      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars 
     856      zcoef_dqla = -Ls * 11637800. * (-5897.8) 
    733857      ! 
    734858      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    … …  
    756880               ! ----------------------------! 
    757881 
    758                ! ... turbulent heat fluxes 
     882               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau 
    759883               ! Sensible Heat 
    760                z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) * ( ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1) ) 
     884               z_qsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1)) 
    761885               ! Latent Heat 
    762                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Cd(ji,jj) * wndm_ice(ji,jj)   & 
    763                   &                         * (  11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
     886               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  & 
     887                  &                ( 11637800. * EXP( -5897.8 / ptsu(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) ) 
    764888               ! Latent heat sensitivity for ice (Dqla/Dt) 
    765889               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN 
    766                   dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Cd(ji,jj) * wndm_ice(ji,jj) / ( zst2 ) * EXP( -5897.8 / ptsu(ji,jj,jl) ) 
     890                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) / zst2 * EXP(-5897.8 / ptsu(ji,jj,jl)) 
    767891               ELSE 
    768892                  dqla_ice(ji,jj,jl) = 0._wp 
    … …  
    770894 
    771895               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice) 
    772                z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Cd(ji,jj) * wndm_ice(ji,jj) 
     896               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj) 
    773897 
    774898               ! ----------------------------! 
    … …  
    787911      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac      ! total precipitation [kg/m2/s] 
    788912      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac      ! solid precipitation [kg/m2/s] 
    789       CALL iom_put( 'snowpre', sprecip * 86400. )                  ! Snow precipitation 
    790       CALL iom_put( 'precip' , tprecip * 86400. )                  ! Total precipitation 
    791  
    792 #if defined  key_lim3 
    793       CALL wrk_alloc( jpi,jpj,   zevap, zsnw ) 
     913      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation 
     914      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation 
    794915 
    795916      ! --- evaporation --- ! 
    … …  
    801922      ! --- evaporation minus precipitation --- ! 
    802923      zsnw(:,:) = 0._wp 
    803       CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing 
    804       emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
     924      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
     925      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    805926      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
    806927      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:) 
    807928 
    808929      ! --- heat flux associated with emp --- ! 
    809       qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst 
     930      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst 
    810931         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair 
    811932         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow) 
    … …  
    815936 
    816937      ! --- total solar and non solar fluxes --- ! 
    817       qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:) 
    818       qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
     938      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  & 
     939         &           + qemp_ice(:,:) + qemp_oce(:,:) 
     940      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 ) 
    819941 
    820942      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- ! 
    … …  
    824946      DO jl = 1, jpl 
    825947         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * cpic * tmask(:,:,1) ) 
    826                                    ! But we do not have Tice => consider it at 0degC => evap=0  
     948         !                         ! But we do not have Tice => consider it at 0degC => evap=0  
    827949      END DO 
    828950 
    829       CALL wrk_dealloc( jpi,jpj,   zevap, zsnw ) 
    830 #endif 
    831  
    832       !-------------------------------------------------------------------- 
    833       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    834       ! thin surface layer and penetrates inside the ice cover 
    835       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    836       ! 
    837       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    838       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    839       ! 
     951      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 
     952      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm 
     953      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
     954      ! 
     955      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     956         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
     957      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
     958         qsr_ice_tr(:,:,:) = qsr_ice(:,:,:) * zfr1 
     959      ELSEWHERE                                                         ! zero when hs>0 
     960         qsr_ice_tr(:,:,:) = 0._wp  
     961      END WHERE 
    840962      ! 
    841963      IF(ln_ctl) THEN 
    … …  
    847969         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ') 
    848970      ENDIF 
    849  
    850       CALL wrk_dealloc( jpi,jpj,jpl,   z_qlw, z_qsb, z_dqlw, z_dqsb ) 
    851       CALL wrk_dealloc( jpi,jpj,       zrhoa ) 
    852       CALL wrk_dealloc( jpi,jpj, Cd ) 
    853       ! 
    854       IF( nn_timing == 1 )  CALL timing_stop('blk_ice_flx') 
    855  
     971      ! 
     972      IF( ln_timing )   CALL timing_stop('blk_ice_flx') 
     973      ! 
    856974   END SUBROUTINE blk_ice_flx 
    857975    
    858 #endif 
    859  
    860    FUNCTION rho_air( ptak, pqa, pslp ) 
    861       !!------------------------------------------------------------------------------- 
    862       !!                           ***  FUNCTION rho_air  *** 
    863       !! 
    864       !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere 
    865       !! 
    866       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)  
    867       !!------------------------------------------------------------------------------- 
    868       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K] 
    869       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg] 
    870       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa] 
    871       REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3] 
    872       !!------------------------------------------------------------------------------- 
    873       ! 
    874       rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  ) 
    875       ! 
    876    END FUNCTION rho_air 
    877  
    878  
    879    FUNCTION cp_air( pqa ) 
    880       !!------------------------------------------------------------------------------- 
    881       !!                           ***  FUNCTION cp_air  *** 
    882       !! 
    883       !! ** Purpose : Compute specific heat (Cp) of moist air 
    884       !! 
    885       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    886       !!------------------------------------------------------------------------------- 
    887       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg] 
    888       REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg] 
    889       !!------------------------------------------------------------------------------- 
    890       ! 
    891       Cp_air = Cp_dry + Cp_vap * pqa 
    892       ! 
    893    END FUNCTION cp_air 
    894  
    895  
    896    FUNCTION q_sat( ptak, pslp ) 
    897       !!---------------------------------------------------------------------------------- 
    898       !!                           ***  FUNCTION q_sat  *** 
    899       !! 
    900       !! ** Purpose : Specific humidity at saturation in [kg/kg] 
    901       !!              Based on accurate estimate of "e_sat" 
    902       !!              aka saturation water vapor (Goff, 1957) 
    903       !! 
    904       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    905       !!---------------------------------------------------------------------------------- 
    906       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K] 
    907       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa] 
    908       REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg] 
    909       ! 
    910       INTEGER  ::   ji, jj         ! dummy loop indices 
    911       REAL(wp) ::   ze_sat, ztmp   ! local scalar 
    912       !!---------------------------------------------------------------------------------- 
    913       ! 
    914       DO jj = 1, jpj 
    915          DO ji = 1, jpi 
    916             ! 
    917             ztmp = rt0 / ptak(ji,jj) 
    918             ! 
    919             ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957) 
    920             ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        & 
    921                &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  & 
    922                &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  ) 
     976 
     977   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     978      !!--------------------------------------------------------------------- 
     979      !!                     ***  ROUTINE blk_ice_qcn  *** 
     980      !! 
     981      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     982      !!                to force sea ice / snow thermodynamics 
     983      !!                in the case JULES coupler is emulated 
     984      !!                 
     985      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     986      !!                following the 0-layer Semtner (1976) approach 
     987      !! 
     988      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     989      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     990      !! 
     991      !!--------------------------------------------------------------------- 
     992      INTEGER                   , INTENT(in   ) ::   k_monocat   ! single-category option 
     993      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu        ! sea ice / snow surface temperature 
     994      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb         ! sea ice base temperature 
     995      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs         ! snow thickness 
     996      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi         ! sea ice thickness 
     997      ! 
     998      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations 
     999      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     1000      ! 
     1001      INTEGER  ::   ji, jj, jl           ! dummy loop indices 
     1002      INTEGER  ::   iter                 ! local integer 
     1003      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars 
     1004      REAL(wp) ::   zkeff_h, ztsu, ztsu0 ! 
     1005      REAL(wp) ::   zqc, zqnet           ! 
     1006      REAL(wp) ::   zhe, zqa0            ! 
     1007      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor 
     1008      !!--------------------------------------------------------------------- 
     1009 
     1010      IF( ln_timing )  CALL timing_start('blk_ice_qcn') 
     1011       
     1012      ! -------------------------------------! 
     1013      !      I   Enhanced conduction factor  ! 
     1014      ! -------------------------------------! 
     1015      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     1016      ! Fichefet and Morales Maqueda, JGR 1997 
     1017      ! 
     1018      zgfac(:,:,:) = 1._wp 
     1019       
     1020      SELECT CASE ( k_monocat )  
     1021      ! 
     1022      CASE ( 1 , 3 ) 
     1023         ! 
     1024         zfac  = 1._wp /  ( rn_cnd_s + rcdic ) 
     1025         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon 
     1026         zfac3 = 2._wp / zepsilon 
     1027         !    
     1028         DO jl = 1, jpl                 
     1029            DO jj = 1 , jpj 
     1030               DO ji = 1, jpi 
     1031                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac                             ! Effective thickness 
     1032                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
     1033               END DO 
     1034            END DO 
     1035         END DO 
     1036         !       
     1037      END SELECT 
     1038       
     1039      ! -------------------------------------------------------------! 
     1040      !      II   Surface temperature and conduction flux            ! 
     1041      ! -------------------------------------------------------------! 
     1042      ! 
     1043      zfac = rcdic * rn_cnd_s 
     1044      ! 
     1045      DO jl = 1, jpl 
     1046         DO jj = 1 , jpj 
     1047            DO ji = 1, jpi 
     1048               !                     
     1049               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
     1050                  &      ( rcdic * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) ) 
     1051               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature 
     1052               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature 
     1053               zqa0    = qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl)  ! Net initial atmospheric heat flux 
    9231054               ! 
    924             q_sat(ji,jj) = reps0 * ze_sat/( 0.01_wp*pslp(ji,jj) - (1._wp - reps0)*ze_sat )   ! 0.01 because SLP is in [Pa] 
    925             ! 
    926          END DO 
    927       END DO 
    928       ! 
    929    END FUNCTION q_sat 
    930  
    931  
    932    FUNCTION gamma_moist( ptak, pqa ) 
    933       !!---------------------------------------------------------------------------------- 
    934       !!                           ***  FUNCTION gamma_moist  *** 
    935       !! 
    936       !! ** Purpose : Compute the moist adiabatic lapse-rate. 
    937       !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate 
    938       !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html 
    939       !! 
    940       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    941       !!---------------------------------------------------------------------------------- 
    942       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K] 
    943       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg] 
    944       REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate 
    945       ! 
    946       INTEGER  ::   ji, jj         ! dummy loop indices 
    947       REAL(wp) :: zrv, ziRT        ! local scalar 
    948       !!---------------------------------------------------------------------------------- 
    949       ! 
    950       DO jj = 1, jpj 
    951          DO ji = 1, jpi 
    952             zrv = pqa(ji,jj) / (1. - pqa(ji,jj)) 
    953             ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT 
    954             gamma_moist(ji,jj) = grav * ( 1. + cevap*zrv*ziRT ) / ( Cp_dry + cevap*cevap*zrv*reps0*ziRT/ptak(ji,jj) ) 
    955          END DO 
    956       END DO 
    957       ! 
    958    END FUNCTION gamma_moist 
    959  
    960  
    961    FUNCTION L_vap( psst ) 
    962       !!--------------------------------------------------------------------------------- 
    963       !!                           ***  FUNCTION L_vap  *** 
    964       !! 
    965       !! ** Purpose : Compute the latent heat of vaporization of water from temperature 
    966       !! 
    967       !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk) 
    968       !!---------------------------------------------------------------------------------- 
    969       REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg] 
    970       REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K] 
    971       !!---------------------------------------------------------------------------------- 
    972       ! 
    973       L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6 
    974       ! 
    975    END FUNCTION L_vap 
    976  
    977  
    978 #if defined key_lim3 
     1055               DO iter = 1, nit     ! --- Iterative loop 
     1056                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards) 
     1057                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget 
     1058                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update 
     1059               END DO 
     1060               ! 
     1061               ptsu   (ji,jj,jl) = MIN( rt0, ztsu ) 
     1062               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     1063               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) 
     1064               qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )   & 
     1065                             &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) ) 
     1066 
     1067            END DO 
     1068         END DO 
     1069         ! 
     1070      END DO  
     1071      !       
     1072      IF( ln_timing )  CALL timing_stop('blk_ice_qcn') 
     1073      ! 
     1074   END SUBROUTINE blk_ice_qcn 
     1075    
     1076 
    9791077   SUBROUTINE Cdn10_Lupkes2012( Cd ) 
    9801078      !!---------------------------------------------------------------------- 
    9811079      !!                      ***  ROUTINE  Cdn10_Lupkes2012  *** 
    9821080      !! 
    983       !! ** Purpose :    Recompute the ice-atm drag at 10m height to make 
    984       !!                 it dependent on edges at leads, melt ponds and flows. 
     1081      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m  
     1082      !!                 to make it dependent on edges at leads, melt ponds and flows. 
    9851083      !!                 After some approximations, this can be resumed to a dependency 
    9861084      !!                 on ice concentration. 
    … …  
    10261124       
    10271125   END SUBROUTINE Cdn10_Lupkes2012 
     1126 
     1127 
     1128   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch ) 
     1129      !!---------------------------------------------------------------------- 
     1130      !!                      ***  ROUTINE  Cdn10_Lupkes2015  *** 
     1131      !! 
     1132      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation 
     1133      !!                 between sea-ice and atmosphere with distinct momentum  
     1134      !!                 and heat coefficients depending on sea-ice concentration  
     1135      !!                 and atmospheric stability (no meltponds effect for now). 
     1136      !!                 
     1137      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015) 
     1138      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme, 
     1139      !!                 it considers specific skin and form drags (Andreas et al. 2010) 
     1140      !!                 to compute neutral transfert coefficients for both heat and  
     1141      !!                 momemtum fluxes. Atmospheric stability effect on transfert 
     1142      !!                 coefficient is also taken into account following Louis (1979). 
     1143      !! 
     1144      !! ** References : Lupkes et al. JGR 2015 (theory) 
     1145      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation) 
     1146      !! 
     1147      !!---------------------------------------------------------------------- 
     1148      ! 
     1149      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd 
     1150      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch 
     1151      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat 
     1152      ! 
     1153      ! ECHAM6 constants 
     1154      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m] 
     1155      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m] 
     1156      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m] 
     1157      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41 
     1158      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41 
     1159      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13 
     1160      REAL(wp), PARAMETER ::   zc2          = zc * zc 
     1161      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14 
     1162      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30 
     1163      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51 
     1164      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56 
     1165      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26 
     1166      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26 
     1167      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma 
     1168      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp 
     1169      ! 
     1170      INTEGER  ::   ji, jj         ! dummy loop indices 
     1171      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu 
     1172      REAL(wp) ::   zrib_o, zrib_i 
     1173      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice 
     1174      REAL(wp) ::   zChn_skin_ice, zChn_form_ice 
     1175      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw 
     1176      REAL(wp) ::   zCdn_form_tmp 
     1177      !!---------------------------------------------------------------------- 
     1178 
     1179      ! Momentum Neutral Transfert Coefficients (should be a constant) 
     1180      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40 
     1181      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7 
     1182      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details) 
     1183      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32) 
     1184 
     1185      ! Heat Neutral Transfert Coefficients 
     1186      zChn_skin_ice = vkarmn**2 / ( LOG( rn_zu / z0_ice + 1._wp ) * LOG( rn_zu * z1_alpha / z0_skin_ice + 1._wp ) )   ! Eq. 50 + Eq. 52 (cf Lupkes email for details) 
     1187      
     1188      ! Atmospheric and Surface Variables 
     1189      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin 
     1190      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg] 
     1191      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
     1192      ! 
     1193      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     1194         DO ji = fs_2, fs_jpim1 
     1195            ! Virtual potential temperature [K] 
     1196            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
     1197            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice 
     1198            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu 
     1199             
     1200            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead) 
     1201            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean 
     1202            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice 
     1203             
     1204            ! Momentum and Heat Neutral Transfert Coefficients 
     1205            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40 
     1206            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53  
     1207                        
     1208            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead) 
     1209            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water 
     1210            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details) 
     1211            IF( zrib_o <= 0._wp ) THEN 
     1212               zfmw = 1._wp - zam * zrib_o / ( 1._wp + 3._wp * zc2 * Cdn_oce(ji,jj) * SQRT( -zrib_o * ( rn_zu / z0w + 1._wp ) ) )  ! Eq. 10 
     1213               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26 
     1214                  &             )**zgamma )**z1_gamma 
     1215            ELSE 
     1216               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12 
     1217               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28 
     1218            ENDIF 
     1219             
     1220            IF( zrib_i <= 0._wp ) THEN 
     1221               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9 
     1222               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25 
     1223            ELSE 
     1224               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11 
     1225               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27 
     1226            ENDIF 
     1227             
     1228            ! Momentum Transfert Coefficients (Eq. 38) 
     1229            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  & 
     1230               &        zCdn_form_ice * ( zfmi * at_i_b(ji,jj) + zfmw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1231             
     1232            ! Heat Transfert Coefficients (Eq. 49) 
     1233            Ch(ji,jj) = zChn_skin_ice *   zfhi +  & 
     1234               &        zChn_form_ice * ( zfhi * at_i_b(ji,jj) + zfhw * ( 1._wp - at_i_b(ji,jj) ) ) / MAX( 1.e-06, at_i_b(ji,jj) ) 
     1235            ! 
     1236         END DO 
     1237      END DO 
     1238      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1239      ! 
     1240   END SUBROUTINE Cdn10_Lupkes2015 
     1241 
    10281242#endif 
    1029     
    10301243 
    10311244   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.