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.
sbcblk.F90 in NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC – NEMO

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/sbcblk.F90 @ 10297

Last change on this file since 10297 was 10297, checked in by smasson, 5 years ago

dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of mppmin/max/sum, see #2133

  • Property svn:keywords set to Id
File size: 70.0 KB
RevLine 
[6723]1MODULE sbcblk
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                         Aerodynamic Bulk Formulas
6   !!                        SUCCESSOR OF "sbcblk_core"
7   !!=====================================================================
[7163]8   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
9   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
10   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
11   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
[6723]12   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
13   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
[7163]14   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
15   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
16   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
[9019]17   !!                 !                        ==> based on AeroBulk (http://aerobulk.sourceforge.net/)
[7163]18   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
[9019]19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)
[6723]20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
[7163]23   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
24   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
[6727]25   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
26   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
27   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
28   !!   q_sat         : saturation humidity as a function of SLP and temperature
29   !!   L_vap         : latent heat of vaporization of water as a function of temperature
[9019]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
[6723]36   !!----------------------------------------------------------------------
37   USE oce            ! ocean dynamics and tracers
38   USE dom_oce        ! ocean space and time domain
39   USE phycst         ! physical constants
40   USE fldread        ! read input fields
41   USE sbc_oce        ! Surface boundary condition: ocean fields
42   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
43   USE sbcdcy         ! surface boundary condition: diurnal cycle
44   USE sbcwave , ONLY :   cdn_wave ! wave module
45   USE sbc_ice        ! Surface boundary condition: ice fields
46   USE lib_fortran    ! to use key_nosignedzero
[9570]47#if defined key_si3
[9938]48   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su, rn_cnd_s, hfx_err_dif
[9019]49   USE icethd_dh      ! for CALL ice_thd_snwblow
[6723]50#endif
51   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
52   USE sbcblk_algo_coare    ! => turb_coare    : COAREv3.0 (Fairall et al. 2003)
53   USE sbcblk_algo_coare3p5 ! => turb_coare3p5 : COAREv3.5 (Edson et al. 2013)
54   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 31)
[6727]55   !
[6723]56   USE iom            ! I/O manager library
57   USE in_out_manager ! I/O manager
58   USE lib_mpp        ! distribued memory computing library
59   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
60   USE prtctl         ! Print control
61
62   IMPLICIT NONE
63   PRIVATE
64
[7163]65   PUBLIC   sbc_blk_init  ! called in sbcmod
66   PUBLIC   sbc_blk       ! called in sbcmod
[9570]67#if defined key_si3
[9019]68   PUBLIC   blk_ice_tau   ! routine called in iceforcing
69   PUBLIC   blk_ice_flx   ! routine called in iceforcing
70   PUBLIC   blk_ice_qcn   ! routine called in iceforcing
71#endif 
[6723]72
[6727]73!!Lolo: should ultimately be moved in the module with all physical constants ?
74!!gm  : In principle, yes.
[6723]75   REAL(wp), PARAMETER ::   Cp_dry = 1005.0       !: Specic heat of dry air, constant pressure      [J/K/kg]
76   REAL(wp), PARAMETER ::   Cp_vap = 1860.0       !: Specic heat of water vapor, constant pressure  [J/K/kg]
77   REAL(wp), PARAMETER ::   R_dry = 287.05_wp     !: Specific gas constant for dry air              [J/K/kg]
78   REAL(wp), PARAMETER ::   R_vap = 461.495_wp    !: Specific gas constant for water vapor          [J/K/kg]
79   REAL(wp), PARAMETER ::   reps0 = R_dry/R_vap   !: ratio of gas constant for dry air and water vapor => ~ 0.622
80   REAL(wp), PARAMETER ::   rctv0 = R_vap/R_dry   !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
81
82   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
83   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
84   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
[7163]85   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
86   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
87   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
88   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
[6723]89   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
90   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
[7421]91   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
[6723]92   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
93
94   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
95
96   !                                             !!! Bulk parameters
[7355]97   REAL(wp), PARAMETER ::   cpa    = 1000.5         ! specific heat of air (only used for ice fluxes now...)
98   REAL(wp), PARAMETER ::   Ls     =    2.839e6     ! latent heat of sublimation
99   REAL(wp), PARAMETER ::   Stef   =    5.67e-8     ! Stefan Boltzmann constant
[9019]100   REAL(wp), PARAMETER ::   Cd_ice =    1.4e-3      ! transfer coefficient over ice
[7355]101   REAL(wp), PARAMETER ::   albo   =    0.066       ! ocean albedo assumed to be constant
[6723]102   !
103   !                           !!* Namelist namsbc_blk : bulk parameters
104   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
105   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
106   LOGICAL  ::   ln_COARE_3p5   ! "COARE 3.5" algorithm   (Edson et al. 2013)
107   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 31)
108   !
109   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
110   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
[9019]111   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation
112   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress
[6723]113   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
114   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
[9019]115!!gm ref namelist initialize it so remove the setting to false below
116   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
117   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
[7355]118   !
[9019]119   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
120   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
121   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
122   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
123   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
124   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
[6723]125
126   INTEGER  ::   nblk           ! choice of the bulk algorithm
127   !                            ! associated indices:
128   INTEGER, PARAMETER ::   np_NCAR      = 1   ! "NCAR" algorithm        (Large and Yeager 2008)
129   INTEGER, PARAMETER ::   np_COARE_3p0 = 2   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
130   INTEGER, PARAMETER ::   np_COARE_3p5 = 3   ! "COARE 3.5" algorithm   (Edson et al. 2013)
131   INTEGER, PARAMETER ::   np_ECMWF     = 4   ! "ECMWF" algorithm       (IFS cycle 31)
132
133   !! * Substitutions
134#  include "vectopt_loop_substitute.h90"
135   !!----------------------------------------------------------------------
[9598]136   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
[10069]137   !! $Id$
[10068]138   !! Software governed by the CeCILL license (see ./LICENSE)
[6723]139   !!----------------------------------------------------------------------
140CONTAINS
141
[7355]142   INTEGER FUNCTION sbc_blk_alloc()
143      !!-------------------------------------------------------------------
144      !!             ***  ROUTINE sbc_blk_alloc ***
145      !!-------------------------------------------------------------------
[9019]146      ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
147         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
[7355]148      !
[10297]149      IF( lk_mpp             )   CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
[7355]150      IF( sbc_blk_alloc /= 0 )   CALL ctl_warn('sbc_blk_alloc: failed to allocate arrays')
151   END FUNCTION sbc_blk_alloc
152
[9019]153
[7163]154   SUBROUTINE sbc_blk_init
155      !!---------------------------------------------------------------------
156      !!                    ***  ROUTINE sbc_blk_init  ***
157      !!
158      !! ** Purpose :   choose and initialize a bulk formulae formulation
159      !!
160      !! ** Method  :
161      !!
162      !!----------------------------------------------------------------------
163      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
164      INTEGER  ::   ios, ierror, ioptio   ! Local integer
165      !!
166      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
167      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
168      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
169      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
170      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
171      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
172         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
173         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
[7355]174         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
[9019]175         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
[7163]176      !!---------------------------------------------------------------------
177      !
[7355]178      !                                      ! allocate sbc_blk_core array
179      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
180      !
[7163]181      !                             !** read bulk namelist 
182      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
183      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
[9168]184901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist', lwp )
[7163]185      !
186      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters
187      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
[9168]188902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist', lwp )
[7163]189      !
190      IF(lwm) WRITE( numond, namsbc_blk )
191      !
192      !                             !** initialization of the chosen bulk formulae (+ check)
193      !                                   !* select the bulk chosen in the namelist and check the choice
[8666]194                                                               ioptio = 0
[7163]195      IF( ln_NCAR      ) THEN   ;   nblk =  np_NCAR        ;   ioptio = ioptio + 1   ;   ENDIF
196      IF( ln_COARE_3p0 ) THEN   ;   nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1   ;   ENDIF
197      IF( ln_COARE_3p5 ) THEN   ;   nblk =  np_COARE_3p5   ;   ioptio = ioptio + 1   ;   ENDIF
198      IF( ln_ECMWF     ) THEN   ;   nblk =  np_ECMWF       ;   ioptio = ioptio + 1   ;   ENDIF
199      !
200      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
201      !
202      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
203         IF( sn_qsr%nfreqh /= 24 )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
204         IF( sn_qsr%ln_tint ) THEN
205            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
206               &           '              ==> We force time interpolation = .false. for qsr' )
207            sn_qsr%ln_tint = .false.
208         ENDIF
209      ENDIF
210      !                                   !* set the bulk structure
211      !                                      !- store namelist information in an array
212      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
213      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
214      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
215      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
216      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
217      !
218      lhftau = ln_taudif                     !- add an extra field if HF stress is used
219      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
220      !
221      !                                      !- allocate the bulk structure
222      ALLOCATE( sf(jfld), STAT=ierror )
223      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
224      DO ifpr= 1, jfld
225         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
226         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
[9168]227         IF( slf_i(ifpr)%nfreqh > 0. .AND. MOD( 3600. * slf_i(ifpr)%nfreqh , REAL(nn_fsbc) * rdt) /= 0. )   &
228            &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
229            &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
[8530]230
[7163]231      END DO
232      !                                      !- fill the bulk structure with namelist informations
233      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
234      !
[7431]235      IF ( ln_wave ) THEN
236      !Activated wave module but neither drag nor stokes drift activated
[9033]237         IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
238            CALL ctl_warn( 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F')
[7431]239      !drag coefficient read from wave model definable only with mfs bulk formulae and core
240         ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN       
241             CALL ctl_stop( 'drag coefficient read from wave model definable only with mfs bulk formulae and core')
242         ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN
243             CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
244         ENDIF
245      ELSE
[9033]246      IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
[7431]247         &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
248         &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
249         &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
[9033]250         &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      & 
[7431]251         &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
252      ENDIF 
253      !
[7163]254      !           
255      IF(lwp) THEN                     !** Control print
256         !
257         WRITE(numout,*)                  !* namelist
258         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
259         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
260         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
261         WRITE(numout,*) '      "COARE 3.5" algorithm   (Edson et al. 2013)         ln_COARE_3p5 = ', ln_COARE_3p0
262         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 31)              ln_ECMWF     = ', ln_ECMWF
263         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
264         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
265         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
266         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
267         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
268         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
269         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
[9019]270         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
271         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
[7163]272         !
273         WRITE(numout,*)
274         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
[9190]275         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
276         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
277         CASE( np_COARE_3p5 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.5" algorithm   (Edson et al. 2013)'
278         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 31)'
[7163]279         END SELECT
280         !
281      ENDIF
282      !
283   END SUBROUTINE sbc_blk_init
284
285
[6723]286   SUBROUTINE sbc_blk( kt )
287      !!---------------------------------------------------------------------
288      !!                    ***  ROUTINE sbc_blk  ***
289      !!
290      !! ** Purpose :   provide at each time step the surface ocean fluxes
[9019]291      !!              (momentum, heat, freshwater and runoff)
[6723]292      !!
293      !! ** Method  : (1) READ each fluxes in NetCDF files:
294      !!      the 10m wind velocity (i-component) (m/s)    at T-point
295      !!      the 10m wind velocity (j-component) (m/s)    at T-point
296      !!      the 10m or 2m specific humidity     ( % )
297      !!      the solar heat                      (W/m2)
298      !!      the Long wave                       (W/m2)
299      !!      the 10m or 2m air temperature       (Kelvin)
300      !!      the total precipitation (rain+snow) (Kg/m2/s)
301      !!      the snow (solid prcipitation)       (kg/m2/s)
302      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
303      !!              (2) CALL blk_oce
304      !!
305      !!      C A U T I O N : never mask the surface stress fields
306      !!                      the stress is assumed to be in the (i,j) mesh referential
307      !!
308      !! ** Action  :   defined at each time-step at the air-sea interface
309      !!              - utau, vtau  i- and j-component of the wind stress
310      !!              - taum        wind stress module at T-point
311      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
312      !!              - qns, qsr    non-solar and solar heat fluxes
313      !!              - emp         upward mass flux (evapo. - precip.)
314      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
315      !!
316      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
317      !!                   Brodeau et al. Ocean Modelling 2010
318      !!----------------------------------------------------------------------
319      INTEGER, INTENT(in) ::   kt   ! ocean time step
320      !!---------------------------------------------------------------------
321      !
322      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
[7163]323      !
[6723]324      !                                            ! compute the surface ocean fluxes using bulk formulea
325      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
326
327#if defined key_cice
328      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
[7753]329         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
330         IF( ln_dm2dc ) THEN ; qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
331         ELSE                ; qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1) 
[7282]332         ENDIF
[7753]333         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
334         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
335         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
336         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
337         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
338         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
[6723]339      ENDIF
340#endif
341      !
342   END SUBROUTINE sbc_blk
343
344
345   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
346      !!---------------------------------------------------------------------
347      !!                     ***  ROUTINE blk_oce  ***
348      !!
349      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
350      !!      the ocean surface at each time step
351      !!
352      !! ** Method  :   bulk formulea for the ocean using atmospheric
353      !!      fields read in sbc_read
354      !!
355      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
356      !!              - vtau    : j-component of the stress at V-point  (N/m2)
357      !!              - taum    : Wind stress module at T-point         (N/m2)
358      !!              - wndm    : Wind speed module at T-point          (m/s)
359      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
360      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
361      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
362      !!
363      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
364      !!---------------------------------------------------------------------
365      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
366      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
367      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
368      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
369      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
370      !
371      INTEGER  ::   ji, jj               ! dummy loop indices
372      REAL(wp) ::   zztmp                ! local variable
[9019]373      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
374      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst
375      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
376      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
377      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
378      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
379      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
380      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa             ! density of air   [kg/m^3]
[6723]381      !!---------------------------------------------------------------------
382      !
[7753]383      ! local scalars ( place there for vector optimisation purposes)
384      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
[6723]385
386      ! ----------------------------------------------------------------------------- !
387      !      0   Wind components and module at T-point relative to the moving ocean   !
388      ! ----------------------------------------------------------------------------- !
389
[7753]390      ! ... components ( U10m - U_oce ) at T-point (unmasked)
391!!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
392      zwnd_i(:,:) = 0._wp
393      zwnd_j(:,:) = 0._wp
[6723]394#if defined key_cyclone
395      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
396      DO jj = 2, jpjm1
397         DO ji = fs_2, fs_jpim1   ! vect. opt.
398            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
399            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
400         END DO
401      END DO
402#endif
403      DO jj = 2, jpjm1
404         DO ji = fs_2, fs_jpim1   ! vect. opt.
405            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
406            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
407         END DO
408      END DO
[10170]409      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
[6723]410      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
[7753]411      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
412         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
[6723]413
414      ! ----------------------------------------------------------------------------- !
415      !      I   Radiative FLUXES                                                     !
416      ! ----------------------------------------------------------------------------- !
417
418      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
419      zztmp = 1. - albo
420      IF( ln_dm2dc ) THEN   ;   qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
[7753]421      ELSE                  ;   qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
[6723]422      ENDIF
423
[7753]424      zqlw(:,:) = (  sf(jp_qlw)%fnow(:,:,1) - Stef * zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:)  ) * tmask(:,:,1)   ! Long  Wave
[6723]425
426      ! ----------------------------------------------------------------------------- !
427      !     II    Turbulent FLUXES                                                    !
428      ! ----------------------------------------------------------------------------- !
429
430      ! ... specific humidity at SST and IST tmask(
[6727]431      zsq(:,:) = 0.98 * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
[6723]432      !!
433      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
434      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
435      !!    (since reanalysis products provide T at z, not theta !)
[6727]436      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
[6723]437
438      SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
[6727]439      !
440      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
[9019]441         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]442      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
[9019]443         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]444      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
[9019]445         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6727]446      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
[9019]447         &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
[6723]448      CASE DEFAULT
[7163]449         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
[6723]450      END SELECT
451
[6727]452      !                          ! Compute true air density :
453      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove zrho*grav*rn_zu from SLP...)
[9019]454         zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
[6727]455      ELSE                                      ! At zt:
456         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[6723]457      END IF
458
[9019]459!!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
460!!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
[7355]461
[6727]462      DO jj = 1, jpj             ! tau module, i and j component
[6723]463         DO ji = 1, jpi
[9019]464            zztmp = zrhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
[6723]465            taum  (ji,jj) = zztmp * wndm  (ji,jj)
466            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
467            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
468         END DO
469      END DO
470
[6727]471      !                          ! add the HF tau contribution to the wind stress module
[7753]472      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
[6727]473
[6723]474      CALL iom_put( "taum_oce", taum )   ! output wind stress module
475
476      ! ... utau, vtau at U- and V_points, resp.
477      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
478      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
479      DO jj = 1, jpjm1
480         DO ji = 1, fs_jpim1
481            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
482               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
483            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
484               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
485         END DO
486      END DO
[10170]487      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
[6723]488
489      !  Turbulent fluxes over ocean
490      ! -----------------------------
491
492      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
[9727]493      zqla(:,:) = zrhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
[6723]494
495      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
496         !! q_air and t_air are given at 10m (wind reference height)
[9019]497         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
498         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
[6723]499      ELSE
500         !! q_air and t_air are not given at 10m (wind reference height)
501         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
[9019]502         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
503         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
[6723]504      ENDIF
505
[6727]506      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
[6723]507
508
509      IF(ln_ctl) THEN
[9019]510         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
511         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
[6723]512         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
513         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
514         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
515            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
516         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
517         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
518      ENDIF
519
520      ! ----------------------------------------------------------------------------- !
521      !     III    Total FLUXES                                                       !
522      ! ----------------------------------------------------------------------------- !
523      !
[7753]524      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
525         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
526      !
527      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
[9935]528         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
[7753]529         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST
530         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
531         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
532         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
[9935]533         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
[9727]534      qns(:,:) = qns(:,:) * tmask(:,:,1)
[7753]535      !
[9570]536#if defined key_si3
[9656]537      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
[7753]538      qsr_oce(:,:) = qsr(:,:)
[6723]539#endif
540      !
541      IF ( nn_ice == 0 ) THEN
542         CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
543         CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
544         CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
545         CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
546         CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
547         CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
548         CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
[9727]549         tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
550         sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
[9019]551         CALL iom_put( 'snowpre', sprecip )                 ! Snow
552         CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
[6723]553      ENDIF
554      !
555      IF(ln_ctl) THEN
556         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
557         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
558         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
559         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
560            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
561      ENDIF
562      !
563   END SUBROUTINE blk_oce
564
565
[9019]566
567   FUNCTION rho_air( ptak, pqa, pslp )
568      !!-------------------------------------------------------------------------------
569      !!                           ***  FUNCTION rho_air  ***
570      !!
571      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
572      !!
573      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
574      !!-------------------------------------------------------------------------------
575      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
576      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
577      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp      ! pressure in                [Pa]
578      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air   ! density of moist air   [kg/m^3]
579      !!-------------------------------------------------------------------------------
580      !
581      rho_air = pslp / (  R_dry*ptak * ( 1._wp + rctv0*pqa )  )
582      !
583   END FUNCTION rho_air
584
585
586   FUNCTION cp_air( pqa )
587      !!-------------------------------------------------------------------------------
588      !!                           ***  FUNCTION cp_air  ***
589      !!
590      !! ** Purpose : Compute specific heat (Cp) of moist air
591      !!
592      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
593      !!-------------------------------------------------------------------------------
594      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
595      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air   ! specific heat of moist air   [J/K/kg]
596      !!-------------------------------------------------------------------------------
597      !
598      Cp_air = Cp_dry + Cp_vap * pqa
599      !
600   END FUNCTION cp_air
601
602
603   FUNCTION q_sat( ptak, pslp )
604      !!----------------------------------------------------------------------------------
605      !!                           ***  FUNCTION q_sat  ***
606      !!
607      !! ** Purpose : Specific humidity at saturation in [kg/kg]
608      !!              Based on accurate estimate of "e_sat"
609      !!              aka saturation water vapor (Goff, 1957)
610      !!
611      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
612      !!----------------------------------------------------------------------------------
613      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak    ! air temperature                       [K]
614      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pslp    ! sea level atmospheric pressure       [Pa]
615      REAL(wp), DIMENSION(jpi,jpj)             ::   q_sat   ! Specific humidity at saturation   [kg/kg]
616      !
617      INTEGER  ::   ji, jj         ! dummy loop indices
618      REAL(wp) ::   ze_sat, ztmp   ! local scalar
619      !!----------------------------------------------------------------------------------
620      !
621      DO jj = 1, jpj
622         DO ji = 1, jpi
623            !
624            ztmp = rt0 / ptak(ji,jj)
625            !
626            ! Vapour pressure at saturation [hPa] : WMO, (Goff, 1957)
627            ze_sat = 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(ptak(ji,jj)/rt0)        &
628               &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(ptak(ji,jj)/rt0 - 1.)) )  &
629               &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614  )
630               !
631            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]
632            !
633         END DO
634      END DO
635      !
636   END FUNCTION q_sat
637
638
639   FUNCTION gamma_moist( ptak, pqa )
640      !!----------------------------------------------------------------------------------
641      !!                           ***  FUNCTION gamma_moist  ***
642      !!
643      !! ** Purpose : Compute the moist adiabatic lapse-rate.
644      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
645      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
646      !!
647      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
648      !!----------------------------------------------------------------------------------
649      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak          ! air temperature       [K]
650      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa           ! specific humidity [kg/kg]
651      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist   ! moist adiabatic lapse-rate
652      !
653      INTEGER  ::   ji, jj         ! dummy loop indices
654      REAL(wp) :: zrv, ziRT        ! local scalar
655      !!----------------------------------------------------------------------------------
656      !
657      DO jj = 1, jpj
658         DO ji = 1, jpi
659            zrv = pqa(ji,jj) / (1. - pqa(ji,jj))
660            ziRT = 1. / (R_dry*ptak(ji,jj))    ! 1/RT
[9935]661            gamma_moist(ji,jj) = grav * ( 1. + rLevap*zrv*ziRT ) / ( Cp_dry + rLevap*rLevap*zrv*reps0*ziRT/ptak(ji,jj) )
[9019]662         END DO
663      END DO
664      !
665   END FUNCTION gamma_moist
666
667
668   FUNCTION L_vap( psst )
669      !!---------------------------------------------------------------------------------
670      !!                           ***  FUNCTION L_vap  ***
671      !!
672      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
673      !!
674      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://sourceforge.net/p/aerobulk)
675      !!----------------------------------------------------------------------------------
676      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap   ! latent heat of vaporization   [J/kg]
677      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
678      !!----------------------------------------------------------------------------------
679      !
680      L_vap = (  2.501 - 0.00237 * ( psst(:,:) - rt0)  ) * 1.e6
681      !
682   END FUNCTION L_vap
683
[9570]684#if defined key_si3
[9019]685   !!----------------------------------------------------------------------
[9570]686   !!   'key_si3'                                       SI3 sea-ice model
[9019]687   !!----------------------------------------------------------------------
688   !!   blk_ice_tau : provide the air-ice stress
689   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
690   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating JULES coupler)
691   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
692   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
693   !!----------------------------------------------------------------------
694
[6723]695   SUBROUTINE blk_ice_tau
696      !!---------------------------------------------------------------------
697      !!                     ***  ROUTINE blk_ice_tau  ***
698      !!
699      !! ** Purpose :   provide the surface boundary condition over sea-ice
700      !!
701      !! ** Method  :   compute momentum using bulk formulation
702      !!                formulea, ice variables and read atmospheric fields.
703      !!                NB: ice drag coefficient is assumed to be a constant
704      !!---------------------------------------------------------------------
705      INTEGER  ::   ji, jj    ! dummy loop indices
[9019]706      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
707      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
708      REAL(wp), DIMENSION(jpi,jpj) ::   zrhoa     ! transfer coefficient for momentum      (tau)
[6723]709      !!---------------------------------------------------------------------
710      !
[9019]711      ! set transfer coefficients to default sea-ice values
712      Cd_atm(:,:) = Cd_ice
713      Ch_atm(:,:) = Cd_ice
714      Ce_atm(:,:) = Cd_ice
[6723]715
[9019]716      wndm_ice(:,:) = 0._wp      !!gm brutal....
[7355]717
[9019]718      ! ------------------------------------------------------------ !
719      !    Wind module relative to the moving ice ( U10m - U_ice )   !
720      ! ------------------------------------------------------------ !
[9767]721      ! C-grid ice dynamics :   U & V-points (same as ocean)
722      DO jj = 2, jpjm1
723         DO ji = fs_2, fs_jpim1   ! vect. opt.
724            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
725            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
726            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
[9019]727         END DO
[9767]728      END DO
[10170]729      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
[9767]730      !
[9019]731      ! Make ice-atm. drag dependent on ice concentration
732      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
733         CALL Cdn10_Lupkes2012( Cd_atm )
734         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
735      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
736         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm ) 
[7355]737      ENDIF
738
[9019]739!!      CALL iom_put( "Cd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
740!!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
741
[6723]742      ! local scalars ( place there for vector optimisation purposes)
743      ! Computing density of air! Way denser that 1.2 over sea-ice !!!
[7355]744      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
[6723]745
[7753]746      !!gm brutal....
747      utau_ice  (:,:) = 0._wp
748      vtau_ice  (:,:) = 0._wp
749      !!gm end
[6723]750
[9019]751      ! ------------------------------------------------------------ !
752      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
753      ! ------------------------------------------------------------ !
[9767]754      ! C-grid ice dynamics :   U & V-points (same as ocean)
755      DO jj = 2, jpjm1
756         DO ji = fs_2, fs_jpim1   ! vect. opt.
757            utau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
758               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
759            vtau_ice(ji,jj) = 0.5 * zrhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
760               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
[6723]761         END DO
[9767]762      END DO
[10170]763      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
[9019]764      !
[9767]765      !
[6723]766      IF(ln_ctl) THEN
767         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
768         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
769      ENDIF
[9019]770      !
[6723]771   END SUBROUTINE blk_ice_tau
772
773
[9019]774   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
[6723]775      !!---------------------------------------------------------------------
776      !!                     ***  ROUTINE blk_ice_flx  ***
777      !!
[9019]778      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
[6723]779      !!
780      !! ** Method  :   compute heat and freshwater exchanged
781      !!                between atmosphere and sea-ice using bulk formulation
782      !!                formulea, ice variables and read atmmospheric fields.
783      !!
784      !! caution : the net upward water flux has with mm/day unit
785      !!---------------------------------------------------------------------
[6727]786      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
[9019]787      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
788      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
[6727]789      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
[6723]790      !!
[6727]791      INTEGER  ::   ji, jj, jl               ! dummy loop indices
[9454]792      REAL(wp) ::   zst3                     ! local variable
[6727]793      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
[9935]794      REAL(wp) ::   zztmp, z1_rLsub           !   -      -
[9019]795      REAL(wp) ::   zfr1, zfr2               ! local variables
[9454]796      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
[9019]797      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
798      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
799      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
800      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
[9656]801      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
[9019]802      REAL(wp), DIMENSION(jpi,jpj)     ::   zrhoa
[6723]803      !!---------------------------------------------------------------------
804      !
[9019]805      zcoef_dqlw = 4.0 * 0.95 * Stef             ! local scalars
806      zcoef_dqla = -Ls * 11637800. * (-5897.8)
[6723]807      !
[6727]808      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
[6723]809      !
810      zztmp = 1. / ( 1. - albo )
[9454]811      WHERE( ptsu(:,:,:) /= 0._wp )   ;   z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
812      ELSEWHERE                       ;   z1_st(:,:,:) = 0._wp
813      END WHERE
[7753]814      !                                     ! ========================== !
815      DO jl = 1, jpl                        !  Loop over ice categories  !
816         !                                  ! ========================== !
[6723]817         DO jj = 1 , jpj
818            DO ji = 1, jpi
819               ! ----------------------------!
820               !      I   Radiative FLUXES   !
821               ! ----------------------------!
[9454]822               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
[6723]823               ! Short Wave (sw)
824               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
825               ! Long  Wave (lw)
826               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - Stef * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
827               ! lw sensitivity
828               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
829
830               ! ----------------------------!
831               !     II    Turbulent FLUXES  !
832               ! ----------------------------!
833
[9019]834               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
[6723]835               ! Sensible Heat
[9019]836               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))
[6723]837               ! Latent Heat
[9019]838               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
[9454]839                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
[6723]840               ! Latent heat sensitivity for ice (Dqla/Dt)
841               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
[9454]842                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
843                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
[6723]844               ELSE
845                  dqla_ice(ji,jj,jl) = 0._wp
846               ENDIF
847
848               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
[9019]849               z_dqsb(ji,jj,jl) = zrhoa(ji,jj) * cpa * Ch_atm(ji,jj) * wndm_ice(ji,jj)
[6723]850
851               ! ----------------------------!
852               !     III    Total FLUXES     !
853               ! ----------------------------!
854               ! Downward Non Solar flux
855               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
856               ! Total non solar heat flux sensitivity for ice
857               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
858            END DO
859            !
860         END DO
861         !
862      END DO
863      !
[9727]864      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
865      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
[9019]866      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
867      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
[6723]868
869      ! --- evaporation --- !
[9935]870      z1_rLsub = 1._wp / rLsub
871      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
872      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
873      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
[6723]874
[7753]875      ! --- evaporation minus precipitation --- !
876      zsnw(:,:) = 0._wp
[9019]877      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
878      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
[7753]879      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
880      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
[6723]881
[7753]882      ! --- heat flux associated with emp --- !
[9019]883      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
[7753]884         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
885         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
[9935]886         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[7753]887      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
[9935]888         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]889
[7753]890      ! --- total solar and non solar fluxes --- !
[9019]891      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
892         &           + qemp_ice(:,:) + qemp_oce(:,:)
893      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
[6723]894
[7753]895      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
[9935]896      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
[6723]897
[7504]898      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
899      DO jl = 1, jpl
[9935]900         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
[9019]901         !                         ! But we do not have Tice => consider it at 0degC => evap=0
[7504]902      END DO
903
[9019]904      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
905      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
906      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
[6723]907      !
[9019]908      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
[9910]909         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
[9019]910      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
[9910]911         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
[9019]912      ELSEWHERE                                                         ! zero when hs>0
[9910]913         qtr_ice_top(:,:,:) = 0._wp 
[9019]914      END WHERE
[6723]915      !
916      IF(ln_ctl) THEN
917         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
918         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
919         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
920         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
921         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
922         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
923      ENDIF
924      !
925   END SUBROUTINE blk_ice_flx
[6727]926   
[6723]927
[9076]928   SUBROUTINE blk_ice_qcn( k_virtual_itd, ptsu, ptb, phs, phi )
[9019]929      !!---------------------------------------------------------------------
930      !!                     ***  ROUTINE blk_ice_qcn  ***
[6723]931      !!
[9019]932      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
933      !!                to force sea ice / snow thermodynamics
934      !!                in the case JULES coupler is emulated
935      !!               
936      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
937      !!                following the 0-layer Semtner (1976) approach
[6727]938      !!
[9019]939      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
940      !!              - qcn_ice : surface inner conduction flux (W/m2)
941      !!
942      !!---------------------------------------------------------------------
[9076]943      INTEGER                   , INTENT(in   ) ::   k_virtual_itd   ! single-category option
944      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
945      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
946      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
947      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
[6723]948      !
[9019]949      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
950      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
[6723]951      !
[9019]952      INTEGER  ::   ji, jj, jl           ! dummy loop indices
953      INTEGER  ::   iter                 ! local integer
954      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
955      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
956      REAL(wp) ::   zqc, zqnet           !
957      REAL(wp) ::   zhe, zqa0            !
958      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
959      !!---------------------------------------------------------------------
960     
961      ! -------------------------------------!
962      !      I   Enhanced conduction factor  !
963      ! -------------------------------------!
[9076]964      ! Emulates the enhancement of conduction by unresolved thin ice (k_virtual_itd = 1/2)
[9019]965      ! Fichefet and Morales Maqueda, JGR 1997
[6723]966      !
[9019]967      zgfac(:,:,:) = 1._wp
968     
[9076]969      SELECT CASE ( k_virtual_itd )
[6723]970      !
[9076]971      CASE ( 1 , 2 )
[9019]972         !
[9935]973         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
[9019]974         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
975         zfac3 = 2._wp / zepsilon
976         !   
977         DO jl = 1, jpl               
978            DO jj = 1 , jpj
979               DO ji = 1, jpi
[9935]980                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
[9019]981                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
982               END DO
983            END DO
984         END DO
985         !     
986      END SELECT
987     
988      ! -------------------------------------------------------------!
989      !      II   Surface temperature and conduction flux            !
990      ! -------------------------------------------------------------!
[6723]991      !
[9935]992      zfac = rcnd_i * rn_cnd_s
[6723]993      !
[9019]994      DO jl = 1, jpl
995         DO jj = 1 , jpj
996            DO ji = 1, jpi
997               !                   
998               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
[9935]999                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
[9019]1000               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1001               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
[9910]1002               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
[6727]1003               !
[9019]1004               DO iter = 1, nit     ! --- Iterative loop
1005                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1006                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1007                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1008               END DO
1009               !
1010               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1011               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1012               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
[9910]1013               qml_ice(ji,jj,jl) = ( qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) - qcn_ice(ji,jj,jl) )  &
[9019]1014                             &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
[6723]1015
[9938]1016               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1017               hfx_err_dif(ji,jj) = hfx_err_dif(ji,jj) - ( dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 ) ) * a_i_b(ji,jj,jl) 
1018
[9019]1019            END DO
[6723]1020         END DO
[9019]1021         !
1022      END DO 
1023      !     
1024   END SUBROUTINE blk_ice_qcn
1025   
[6723]1026
[7355]1027   SUBROUTINE Cdn10_Lupkes2012( Cd )
1028      !!----------------------------------------------------------------------
1029      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1030      !!
[9019]1031      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1032      !!                 to make it dependent on edges at leads, melt ponds and flows.
[7355]1033      !!                 After some approximations, this can be resumed to a dependency
1034      !!                 on ice concentration.
1035      !!               
1036      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1037      !!                 with the highest level of approximation: level4, eq.(59)
1038      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1039      !!
1040      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1041      !!
1042      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1043      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1044      !!                    A is the concentration of ice minus melt ponds (if any)
1045      !!
1046      !!                 This new drag has a parabolic shape (as a function of A) starting at
1047      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1048      !!                 and going down to Cdi(say 1.4e-3) for A=1
1049      !!
[7507]1050      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
[7355]1051      !!                 => see Lupkes et al (2013)
1052      !!
1053      !! ** References : Lupkes et al. JGR 2012 (theory)
1054      !!                 Lupkes et al. GRL 2013 (application to GCM)
1055      !!
1056      !!----------------------------------------------------------------------
1057      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1058      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1059      REAL(wp), PARAMETER ::   znu   = 1._wp
1060      REAL(wp), PARAMETER ::   zmu   = 1._wp
1061      REAL(wp), PARAMETER ::   zbeta = 1._wp
1062      REAL(wp)            ::   zcoef
1063      !!----------------------------------------------------------------------
1064      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1065
1066      ! generic drag over a cell partly covered by ice
[7507]1067      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1068      !!   &      Cd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1069      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
[7355]1070
1071      ! ice-atm drag
[7507]1072      Cd(:,:) = Cd_ice +  &                                                         ! pure ice drag
1073         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
[7355]1074     
1075   END SUBROUTINE Cdn10_Lupkes2012
[9019]1076
1077
1078   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1079      !!----------------------------------------------------------------------
1080      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1081      !!
1082      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1083      !!                 between sea-ice and atmosphere with distinct momentum
1084      !!                 and heat coefficients depending on sea-ice concentration
1085      !!                 and atmospheric stability (no meltponds effect for now).
1086      !!               
1087      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1088      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1089      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1090      !!                 to compute neutral transfert coefficients for both heat and
1091      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1092      !!                 coefficient is also taken into account following Louis (1979).
1093      !!
1094      !! ** References : Lupkes et al. JGR 2015 (theory)
1095      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1096      !!
1097      !!----------------------------------------------------------------------
1098      !
1099      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1100      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1101      REAL(wp), DIMENSION(jpi,jpj)            ::   zst, zqo_sat, zqi_sat
1102      !
1103      ! ECHAM6 constants
1104      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1105      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1106      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1107      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1108      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1109      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1110      REAL(wp), PARAMETER ::   zc2          = zc * zc
1111      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1112      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1113      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1114      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1115      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1116      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1117      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1118      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1119      !
1120      INTEGER  ::   ji, jj         ! dummy loop indices
1121      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1122      REAL(wp) ::   zrib_o, zrib_i
1123      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1124      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1125      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1126      REAL(wp) ::   zCdn_form_tmp
1127      !!----------------------------------------------------------------------
1128
1129      ! Momentum Neutral Transfert Coefficients (should be a constant)
1130      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1131      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1132      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1133      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1134
1135      ! Heat Neutral Transfert Coefficients
1136      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)
1137     
1138      ! Atmospheric and Surface Variables
1139      zst(:,:)     = sst_m(:,:) + rt0                                       ! convert SST from Celcius to Kelvin
1140      zqo_sat(:,:) = 0.98_wp * q_sat( zst(:,:)  , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1141      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg]
1142      !
1143      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1144         DO ji = fs_2, fs_jpim1
1145            ! Virtual potential temperature [K]
1146            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1147            zthetav_is = tm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1148            zthetav_zu = t_zu (ji,jj) * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1149           
1150            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1151            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1152            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1153           
1154            ! Momentum and Heat Neutral Transfert Coefficients
1155            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1156            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1157                       
1158            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1159            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1160            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1161            IF( zrib_o <= 0._wp ) THEN
1162               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
1163               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1164                  &             )**zgamma )**z1_gamma
1165            ELSE
1166               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1167               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1168            ENDIF
1169           
1170            IF( zrib_i <= 0._wp ) THEN
1171               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1172               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1173            ELSE
1174               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1175               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1176            ENDIF
1177           
1178            ! Momentum Transfert Coefficients (Eq. 38)
1179            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1180               &        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) )
1181           
1182            ! Heat Transfert Coefficients (Eq. 49)
1183            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1184               &        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) )
1185            !
1186         END DO
1187      END DO
[10170]1188      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
[9019]1189      !
1190   END SUBROUTINE Cdn10_Lupkes2015
1191
[7355]1192#endif
1193
[6723]1194   !!======================================================================
1195END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.