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

Last change on this file since 11665 was 11665, checked in by laurent, 4 years ago

LB: syntax improvement.

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