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

Last change on this file since 11217 was 11217, checked in by laurent, 22 months ago

LB: more readable "cool-skin/warm-layer" organisation + no "cool-skin/warm-layer" scheme in the presence of sea-ice.

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