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

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

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

  • Property svn:keywords set to Id
File size: 71.3 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      !PRINT *, '*** LB: allocating tsk!'
156      ALLOCATE( tsk(jpi,jpj), STAT=sbc_blk_cswl_alloc )
157      !PRINT *, '*** 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)) !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) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of d_air and slp !'
513         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
514      CASE( np_humi_rlh )
515         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => computing q_air out of RH, t_air and slp !'
516         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day=', nsec_day
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
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) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
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=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
668         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
669         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
670         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
671         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
672            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
673         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
674         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
675      ENDIF
676
677      ! ----------------------------------------------------------------------------- !
678      !     IV    Total FLUXES                                                       !
679      ! ----------------------------------------------------------------------------- !
680      !
681      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
682         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
683      !
684      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
685         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
686         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
687         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
688         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
689         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
690         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
691      qns(:,:) = qns(:,:) * tmask(:,:,1)
692      !
693#if defined key_si3
694      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
695      qsr_oce(:,:) = qsr(:,:)
696#endif
697      !
698      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
699      CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB
700      CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
701      CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
702      CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
703      CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
704      CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
705      CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
706      CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
707      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
708      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
709      CALL iom_put( 'snowpre', sprecip )                 ! Snow
710      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
711      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
712         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
713         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
714      END IF
715      !!#LB. ENDIF
716      !
717      IF(ln_ctl) THEN
718         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
719         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
720         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
721         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
722            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
723      ENDIF
724      !
725   END SUBROUTINE blk_oce
726
727
728
729#if defined key_si3
730   !!----------------------------------------------------------------------
731   !!   'key_si3'                                       SI3 sea-ice model
732   !!----------------------------------------------------------------------
733   !!   blk_ice_tau : provide the air-ice stress
734   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
735   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
736   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
737   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
738   !!----------------------------------------------------------------------
739
740   SUBROUTINE blk_ice_tau
741      !!---------------------------------------------------------------------
742      !!                     ***  ROUTINE blk_ice_tau  ***
743      !!
744      !! ** Purpose :   provide the surface boundary condition over sea-ice
745      !!
746      !! ** Method  :   compute momentum using bulk formulation
747      !!                formulea, ice variables and read atmospheric fields.
748      !!                NB: ice drag coefficient is assumed to be a constant
749      !!---------------------------------------------------------------------
750      INTEGER  ::   ji, jj    ! dummy loop indices
751      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
752      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
753      !!---------------------------------------------------------------------
754      !
755      ! set transfer coefficients to default sea-ice values
756      Cd_atm(:,:) = rCd_ice
757      Ch_atm(:,:) = rCd_ice
758      Ce_atm(:,:) = rCd_ice
759
760      wndm_ice(:,:) = 0._wp      !!gm brutal....
761
762      ! ------------------------------------------------------------ !
763      !    Wind module relative to the moving ice ( U10m - U_ice )   !
764      ! ------------------------------------------------------------ !
765      ! C-grid ice dynamics :   U & V-points (same as ocean)
766      DO jj = 2, jpjm1
767         DO ji = fs_2, fs_jpim1   ! vect. opt.
768            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
769            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
770            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
771         END DO
772      END DO
773      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
774      !
775      ! Make ice-atm. drag dependent on ice concentration
776      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
777         CALL Cdn10_Lupkes2012( Cd_atm )
778         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
779      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
780         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
781      ENDIF
782
783      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
784      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
785
786      ! local scalars ( place there for vector optimisation purposes)
787
788      !!gm brutal....
789      utau_ice  (:,:) = 0._wp
790      vtau_ice  (:,:) = 0._wp
791      !!gm end
792
793      ! ------------------------------------------------------------ !
794      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
795      ! ------------------------------------------------------------ !
796      ! C-grid ice dynamics :   U & V-points (same as ocean)
797      DO jj = 2, jpjm1
798         DO ji = fs_2, fs_jpim1   ! vect. opt.
799            utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
800               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
801            vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
802               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
803         END DO
804      END DO
805      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
806      !
807      !
808      IF(ln_ctl) THEN
809         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
810         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
811      ENDIF
812      !
813   END SUBROUTINE blk_ice_tau
814
815
816   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
817      !!---------------------------------------------------------------------
818      !!                     ***  ROUTINE blk_ice_flx  ***
819      !!
820      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
821      !!
822      !! ** Method  :   compute heat and freshwater exchanged
823      !!                between atmosphere and sea-ice using bulk formulation
824      !!                formulea, ice variables and read atmmospheric fields.
825      !!
826      !! caution : the net upward water flux has with mm/day unit
827      !!---------------------------------------------------------------------
828      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
829      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
830      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
831      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
832      !!
833      INTEGER  ::   ji, jj, jl               ! dummy loop indices
834      REAL(wp) ::   zst3                     ! local variable
835      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
836      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
837      REAL(wp) ::   zfr1, zfr2               ! local variables
838      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
839      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
840      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
841      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
842      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
843      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
844      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
845      !!---------------------------------------------------------------------
846      !
847      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
848      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
849      !
850
851      !LB:
852      SELECT CASE( nhumi )
853      CASE( np_humi_sph )
854         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
855      CASE( np_humi_dpt )
856         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
857         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
858      CASE( np_humi_rlh )
859         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
860         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
861      END SELECT
862      !LB.
863
864      zztmp = 1. / ( 1. - albo )
865      WHERE( ptsu(:,:,:) /= 0._wp )
866         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
867      ELSEWHERE
868         z1_st(:,:,:) = 0._wp
869      END WHERE
870      !                                     ! ========================== !
871      DO jl = 1, jpl                        !  Loop over ice categories  !
872         !                                  ! ========================== !
873         DO jj = 1 , jpj
874            DO ji = 1, jpi
875               ! ----------------------------!
876               !      I   Radiative FLUXES   !
877               ! ----------------------------!
878               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
879               ! Short Wave (sw)
880               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
881               ! Long  Wave (lw)
882               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
883               ! lw sensitivity
884               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
885
886               ! ----------------------------!
887               !     II    Turbulent FLUXES  !
888               ! ----------------------------!
889
890               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
891               ! Sensible Heat
892               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))
893               ! Latent Heat
894               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
895                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
896               ! Latent heat sensitivity for ice (Dqla/Dt)
897               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
898                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
899                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
900               ELSE
901                  dqla_ice(ji,jj,jl) = 0._wp
902               ENDIF
903
904               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
905               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
906
907               ! ----------------------------!
908               !     III    Total FLUXES     !
909               ! ----------------------------!
910               ! Downward Non Solar flux
911               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
912               ! Total non solar heat flux sensitivity for ice
913               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
914            END DO
915            !
916         END DO
917         !
918      END DO
919      !
920      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
921      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
922      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
923      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
924
925      ! --- evaporation --- !
926      z1_rLsub = 1._wp / rLsub
927      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
928      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
929      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
930
931      ! --- evaporation minus precipitation --- !
932      zsnw(:,:) = 0._wp
933      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
934      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
935      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
936      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
937
938      ! --- heat flux associated with emp --- !
939      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
940         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
941         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
942         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
943      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
944         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
945
946      ! --- total solar and non solar fluxes --- !
947      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
948         &           + qemp_ice(:,:) + qemp_oce(:,:)
949      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
950
951      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
952      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
953
954      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
955      DO jl = 1, jpl
956         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
957         !                         ! But we do not have Tice => consider it at 0degC => evap=0
958      END DO
959
960      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
961      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
962      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
963      !
964      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
965         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
966      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
967         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
968      ELSEWHERE                                                         ! zero when hs>0
969         qtr_ice_top(:,:,:) = 0._wp
970      END WHERE
971      !
972      IF(ln_ctl) THEN
973         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
974         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
975         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
976         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
977         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
978         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
979      ENDIF
980      !
981   END SUBROUTINE blk_ice_flx
982
983
984   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
985      !!---------------------------------------------------------------------
986      !!                     ***  ROUTINE blk_ice_qcn  ***
987      !!
988      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
989      !!                to force sea ice / snow thermodynamics
990      !!                in the case conduction flux is emulated
991      !!
992      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
993      !!                following the 0-layer Semtner (1976) approach
994      !!
995      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
996      !!              - qcn_ice : surface inner conduction flux (W/m2)
997      !!
998      !!---------------------------------------------------------------------
999      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
1000      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
1001      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
1002      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
1003      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
1004      !
1005      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
1006      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
1007      !
1008      INTEGER  ::   ji, jj, jl           ! dummy loop indices
1009      INTEGER  ::   iter                 ! local integer
1010      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
1011      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
1012      REAL(wp) ::   zqc, zqnet           !
1013      REAL(wp) ::   zhe, zqa0            !
1014      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
1015      !!---------------------------------------------------------------------
1016
1017      ! -------------------------------------!
1018      !      I   Enhanced conduction factor  !
1019      ! -------------------------------------!
1020      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1021      ! Fichefet and Morales Maqueda, JGR 1997
1022      !
1023      zgfac(:,:,:) = 1._wp
1024
1025      IF( ld_virtual_itd ) THEN
1026         !
1027         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
1028         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1029         zfac3 = 2._wp / zepsilon
1030         !
1031         DO jl = 1, jpl
1032            DO jj = 1 , jpj
1033               DO ji = 1, jpi
1034                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
1035                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1036               END DO
1037            END DO
1038         END DO
1039         !
1040      ENDIF
1041
1042      ! -------------------------------------------------------------!
1043      !      II   Surface temperature and conduction flux            !
1044      ! -------------------------------------------------------------!
1045      !
1046      zfac = rcnd_i * rn_cnd_s
1047      !
1048      DO jl = 1, jpl
1049         DO jj = 1 , jpj
1050            DO ji = 1, jpi
1051               !
1052               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1053                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1054               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1055               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1056               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1057               !
1058               DO iter = 1, nit     ! --- Iterative loop
1059                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1060                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1061                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1062               END DO
1063               !
1064               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1065               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1066               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1067               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) )  &
1068                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1069
1070               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1071               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)
1072
1073            END DO
1074         END DO
1075         !
1076      END DO
1077      !
1078   END SUBROUTINE blk_ice_qcn
1079
1080
1081   SUBROUTINE Cdn10_Lupkes2012( Cd )
1082      !!----------------------------------------------------------------------
1083      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1084      !!
1085      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1086      !!                 to make it dependent on edges at leads, melt ponds and flows.
1087      !!                 After some approximations, this can be resumed to a dependency
1088      !!                 on ice concentration.
1089      !!
1090      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1091      !!                 with the highest level of approximation: level4, eq.(59)
1092      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1093      !!
1094      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1095      !!
1096      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1097      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1098      !!                    A is the concentration of ice minus melt ponds (if any)
1099      !!
1100      !!                 This new drag has a parabolic shape (as a function of A) starting at
1101      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1102      !!                 and going down to Cdi(say 1.4e-3) for A=1
1103      !!
1104      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1105      !!                 => see Lupkes et al (2013)
1106      !!
1107      !! ** References : Lupkes et al. JGR 2012 (theory)
1108      !!                 Lupkes et al. GRL 2013 (application to GCM)
1109      !!
1110      !!----------------------------------------------------------------------
1111      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1112      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1113      REAL(wp), PARAMETER ::   znu   = 1._wp
1114      REAL(wp), PARAMETER ::   zmu   = 1._wp
1115      REAL(wp), PARAMETER ::   zbeta = 1._wp
1116      REAL(wp)            ::   zcoef
1117      !!----------------------------------------------------------------------
1118      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1119
1120      ! generic drag over a cell partly covered by ice
1121      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1122      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1123      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1124
1125      ! ice-atm drag
1126      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
1127         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1128
1129   END SUBROUTINE Cdn10_Lupkes2012
1130
1131
1132   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1133      !!----------------------------------------------------------------------
1134      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1135      !!
1136      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1137      !!                 between sea-ice and atmosphere with distinct momentum
1138      !!                 and heat coefficients depending on sea-ice concentration
1139      !!                 and atmospheric stability (no meltponds effect for now).
1140      !!
1141      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1142      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1143      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1144      !!                 to compute neutral transfert coefficients for both heat and
1145      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1146      !!                 coefficient is also taken into account following Louis (1979).
1147      !!
1148      !! ** References : Lupkes et al. JGR 2015 (theory)
1149      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1150      !!
1151      !!----------------------------------------------------------------------
1152      !
1153      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1154      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1155      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
1156      !
1157      ! ECHAM6 constants
1158      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1159      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1160      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1161      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1162      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1163      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1164      REAL(wp), PARAMETER ::   zc2          = zc * zc
1165      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1166      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1167      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1168      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1169      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1170      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1171      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1172      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1173      !
1174      INTEGER  ::   ji, jj         ! dummy loop indices
1175      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1176      REAL(wp) ::   zrib_o, zrib_i
1177      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1178      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1179      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1180      REAL(wp) ::   zCdn_form_tmp
1181      !!----------------------------------------------------------------------
1182
1183      ! mean temperature
1184      WHERE( at_i_b(:,:) > 1.e-20 )
1185         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1186      ELSEWHERE
1187         ztm_su(:,:) = rt0
1188      ENDWHERE
1189
1190      ! Momentum Neutral Transfert Coefficients (should be a constant)
1191      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1192      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1193      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1194      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1195
1196      ! Heat Neutral Transfert Coefficients
1197      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)
1198
1199      ! Atmospheric and Surface Variables
1200      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1201      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1202      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
1203      !
1204      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1205         DO ji = fs_2, fs_jpim1
1206            ! Virtual potential temperature [K]
1207            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1208            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1209            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1210
1211            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1212            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1213            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1214
1215            ! Momentum and Heat Neutral Transfert Coefficients
1216            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1217            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1218
1219            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1220            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1221            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1222            IF( zrib_o <= 0._wp ) THEN
1223               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
1224               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1225                  &             )**zgamma )**z1_gamma
1226            ELSE
1227               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1228               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1229            ENDIF
1230
1231            IF( zrib_i <= 0._wp ) THEN
1232               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1233               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1234            ELSE
1235               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1236               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1237            ENDIF
1238
1239            ! Momentum Transfert Coefficients (Eq. 38)
1240            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1241               &        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) )
1242
1243            ! Heat Transfert Coefficients (Eq. 49)
1244            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1245               &        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) )
1246            !
1247         END DO
1248      END DO
1249      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
1250      !
1251   END SUBROUTINE Cdn10_Lupkes2015
1252
1253#endif
1254
1255   !!======================================================================
1256END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.