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 @ 11111

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

LB: skin temperature now used by ECMWF bulk algorithm, and can be xios-ed ; new module "SBC/sbcblk_skin.F90" ; all physical constants defined in SBC now moved to "DOM/phycst.F90"; all air-properties and other ABL functions gathered in a new module: "SBC/sbcblk_phymbl.F90".

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