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/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 11215

Last change on this file since 11215 was 11215, checked in by laurent, 5 years ago

LB: inclusion of "cool-skin/warm-layer" support (a la ECMWF) into COARE 3.0 algorithm.

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