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

Last change on this file since 11615 was 11615, checked in by laurent, 16 months 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.1 KB
Line 
1MODULE sbcblk
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                         Aerodynamic Bulk Formulas
6   !!                        SUCCESSOR OF "sbcblk_core"
7   !!=====================================================================
8   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
9   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
10   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
11   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
12   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
13   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
14   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
15   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
16   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
17   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
18   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)
20   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
24   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
25   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
26   !!             sea-ice case only :
27   !!   blk_ice_tau   : provide the air-ice stress
28   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
29   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
30   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
31   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
32   !!----------------------------------------------------------------------
33   USE oce            ! ocean dynamics and tracers
34   USE dom_oce        ! ocean space and time domain
35   USE phycst         ! physical constants
36   USE fldread        ! read input fields
37   USE sbc_oce        ! Surface boundary condition: ocean fields
38   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
39   USE sbcdcy         ! surface boundary condition: diurnal cycle
40   USE sbcwave , ONLY :   cdn_wave ! wave module
41   USE sbc_ice        ! Surface boundary condition: ice fields
42   USE lib_fortran    ! to use key_nosignedzero
43#if defined key_si3
44   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
45   USE icethd_dh      ! for CALL ice_thd_snwblow
46#endif
47   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
48   USE sbcblk_algo_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            CALL turb_coare3p0 ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                            &  ! COARE v3.0
532               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
533               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
534               &                Tsk_b=tsk(:,:) )
535
536         CASE( np_COARE_3p6 )
537            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITH CSWL options!!!, nsec_day=', nsec_day
538            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
539               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,    &
540               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1), &
541               &                isecday_utc=nsec_day, plong=glamt(:,:) )
542
543         CASE( np_ECMWF     )
544            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITH CSWL options!!!, nsec_day=', nsec_day
545            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
546               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,   &
547               &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
548
549         CASE DEFAULT
550            CALL ctl_stop( 'STOP', 'sbc_oce: unsuported bulk formula selection for "ln_skin_*==.true."' )
551         END SELECT
552
553         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq:
554         WHERE ( fr_i < 0.001_wp )
555            ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep it!!!
556            zst(:,:) = zst(:,:)*tmask(:,:,1)
557            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
558         ELSEWHERE
559            ! we forget about the update...
560            zst(:,:) = zqsb(:,:) !LB: using what we backed up before skin-algo
561            zsq(:,:) = zqla(:,:) !LB:  "   "   "
562         END WHERE
563
564         !LB: Update of tsk, the "official" array for skin temperature
565         tsk(:,:) = zst(:,:)
566
567      ELSE !IF ( ln_skin_cs .OR. ln_skin_wl )
568
569
570         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
571            !
572         CASE( np_NCAR      )
573            CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! NCAR-COREv2
574               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
575
576         CASE( np_COARE_3p0 )
577            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p0" WITHOUT CSWL optional arrays!!!'
578            CALL turb_coare3p0   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,   &  ! COARE v3.0
579               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
580
581         CASE( np_COARE_3p6 )
582            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_coare3p6" WITHOUT CSWL optional arrays!!!'
583            CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,&  ! COARE v3.6
584               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
585
586         CASE( np_ECMWF     )
587            IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => calling "turb_ecmwf" WITHOUT CSWL optional arrays!!!'
588            CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,    &  ! ECMWF
589               &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
590           
591         CASE DEFAULT
592            CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
593         END SELECT
594
595      END IF ! IF ( ln_skin_cs .OR. ln_skin_wl )
596
597
598
599
600      !                          ! Compute true air density :
601      IF( ABS(rn_zu - rn_zqt) > 0.01 ) THEN     ! At zu: (probably useless to remove rhoa*grav*rn_zu from SLP...)
602         rhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
603      ELSE                                      ! At zt:
604         rhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), zqair(:,:), sf(jp_slp)%fnow(:,:,1) )
605      END IF
606
607      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
608      !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
609
610      DO jj = 1, jpj             ! tau module, i and j component
611         DO ji = 1, jpi
612            zztmp = rhoa(ji,jj)  * zU_zu(ji,jj) * Cd_atm(ji,jj)   ! using bulk wind speed
613            taum  (ji,jj) = zztmp * wndm  (ji,jj)
614            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
615            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
616         END DO
617      END DO
618
619      !                          ! add the HF tau contribution to the wind stress module
620      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
621
622      CALL iom_put( "taum_oce", taum )   ! output wind stress module
623
624      ! ... utau, vtau at U- and V_points, resp.
625      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
626      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
627      DO jj = 1, jpjm1
628         DO ji = 1, fs_jpim1
629            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
630               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
631            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
632               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
633         END DO
634      END DO
635      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
636
637      !  Turbulent fluxes over ocean
638      ! -----------------------------
639
640      ! zqla used as temporary array, for rho*U (common term of bulk formulae):
641      zqla(:,:) = rhoa(:,:) * zU_zu(:,:) * tmask(:,:,1)
642
643      IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
644         !! q_air and t_air are given at 10m (wind reference height)
645         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - zqair(:,:)) ) ! Evaporation, using bulk wind speed
646         zqsb (:,:) = cp_air(zqair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
647      ELSE
648         !! q_air and t_air are not given at 10m (wind reference height)
649         ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
650         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
651         zqsb (:,:) = cp_air(zqair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
652      ENDIF
653
654      zqla(:,:) = L_vap(zst(:,:)) * zevap(:,:)     ! Latent Heat flux
655
656
657      ! ----------------------------------------------------------------------------- !
658      !     III    Net longwave radiative FLUX                                        !
659      ! ----------------------------------------------------------------------------- !
660
661      !! 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.)
662      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
663
664
665      IF(ln_ctl) THEN
666         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
667         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
668         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
669         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
670         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
671            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
672         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
673         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
674      ENDIF
675
676      ! ----------------------------------------------------------------------------- !
677      !     IV    Total FLUXES                                                       !
678      ! ----------------------------------------------------------------------------- !
679      !
680      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
681         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
682      !
683      qns(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                &   ! Downward Non Solar
684         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
685         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
686         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
687         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
688         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
689         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
690      qns(:,:) = qns(:,:) * tmask(:,:,1)
691      !
692#if defined key_si3
693      qns_oce(:,:) = zqlw(:,:) - zqsb(:,:) - zqla(:,:)                                ! non solar without emp (only needed by SI3)
694      qsr_oce(:,:) = qsr(:,:)
695#endif
696      !
697      !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN
698      CALL iom_put( "rho_air" ,   rhoa )                 ! output air density (kg/m^3) !#LB
699      CALL iom_put( "qlw_oce" ,   zqlw )                 ! output downward longwave heat over the ocean
700      CALL iom_put( "qsb_oce" , - zqsb )                 ! output downward sensible heat over the ocean
701      CALL iom_put( "qla_oce" , - zqla )                 ! output downward latent   heat over the ocean
702      CALL iom_put( "qemp_oce",   qns-zqlw+zqsb+zqla )   ! output downward heat content of E-P over the ocean
703      CALL iom_put( "qns_oce" ,   qns  )                 ! output downward non solar heat over the ocean
704      CALL iom_put( "qsr_oce" ,   qsr  )                 ! output downward solar heat over the ocean
705      CALL iom_put( "qt_oce"  ,   qns+qsr )              ! output total downward heat over the ocean
706      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
707      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
708      CALL iom_put( 'snowpre', sprecip )                 ! Snow
709      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
710      IF ( ln_skin_cs .OR. ln_skin_wl ) THEN
711         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
712         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
713      END IF
714      !!#LB. ENDIF
715      !
716      IF(ln_ctl) THEN
717         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
718         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
719         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
720         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
721            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
722      ENDIF
723      !
724   END SUBROUTINE blk_oce
725
726
727
728#if defined key_si3
729   !!----------------------------------------------------------------------
730   !!   'key_si3'                                       SI3 sea-ice model
731   !!----------------------------------------------------------------------
732   !!   blk_ice_tau : provide the air-ice stress
733   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
734   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
735   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
736   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
737   !!----------------------------------------------------------------------
738
739   SUBROUTINE blk_ice_tau
740      !!---------------------------------------------------------------------
741      !!                     ***  ROUTINE blk_ice_tau  ***
742      !!
743      !! ** Purpose :   provide the surface boundary condition over sea-ice
744      !!
745      !! ** Method  :   compute momentum using bulk formulation
746      !!                formulea, ice variables and read atmospheric fields.
747      !!                NB: ice drag coefficient is assumed to be a constant
748      !!---------------------------------------------------------------------
749      INTEGER  ::   ji, jj    ! dummy loop indices
750      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
751      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
752      !!---------------------------------------------------------------------
753      !
754      ! set transfer coefficients to default sea-ice values
755      Cd_atm(:,:) = rCd_ice
756      Ch_atm(:,:) = rCd_ice
757      Ce_atm(:,:) = rCd_ice
758
759      wndm_ice(:,:) = 0._wp      !!gm brutal....
760
761      ! ------------------------------------------------------------ !
762      !    Wind module relative to the moving ice ( U10m - U_ice )   !
763      ! ------------------------------------------------------------ !
764      ! C-grid ice dynamics :   U & V-points (same as ocean)
765      DO jj = 2, jpjm1
766         DO ji = fs_2, fs_jpim1   ! vect. opt.
767            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
768            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
769            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
770         END DO
771      END DO
772      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
773      !
774      ! Make ice-atm. drag dependent on ice concentration
775      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
776         CALL Cdn10_Lupkes2012( Cd_atm )
777         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
778      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
779         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
780      ENDIF
781
782      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
783      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
784
785      ! local scalars ( place there for vector optimisation purposes)
786
787      !!gm brutal....
788      utau_ice  (:,:) = 0._wp
789      vtau_ice  (:,:) = 0._wp
790      !!gm end
791
792      ! ------------------------------------------------------------ !
793      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
794      ! ------------------------------------------------------------ !
795      ! C-grid ice dynamics :   U & V-points (same as ocean)
796      DO jj = 2, jpjm1
797         DO ji = fs_2, fs_jpim1   ! vect. opt.
798            utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
799               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
800            vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
801               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
802         END DO
803      END DO
804      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
805      !
806      !
807      IF(ln_ctl) THEN
808         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
809         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
810      ENDIF
811      !
812   END SUBROUTINE blk_ice_tau
813
814
815   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
816      !!---------------------------------------------------------------------
817      !!                     ***  ROUTINE blk_ice_flx  ***
818      !!
819      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
820      !!
821      !! ** Method  :   compute heat and freshwater exchanged
822      !!                between atmosphere and sea-ice using bulk formulation
823      !!                formulea, ice variables and read atmmospheric fields.
824      !!
825      !! caution : the net upward water flux has with mm/day unit
826      !!---------------------------------------------------------------------
827      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
828      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
829      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
830      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
831      !!
832      INTEGER  ::   ji, jj, jl               ! dummy loop indices
833      REAL(wp) ::   zst3                     ! local variable
834      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
835      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
836      REAL(wp) ::   zfr1, zfr2               ! local variables
837      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
838      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
839      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
840      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
841      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
842      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
843      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
844      !!---------------------------------------------------------------------
845      !
846      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
847      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
848      !
849
850      !LB:
851      SELECT CASE( nhumi )
852      CASE( np_humi_sph )
853         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
854      CASE( np_humi_dpt )
855         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of d_air and slp !'
856         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
857      CASE( np_humi_rlh )
858         IF (lwp) PRINT *, ' *** LB(sbcblk.F90) => ICE !!! computing q_air out of RH, t_air and slp !'
859         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
860      END SELECT
861      !LB.
862
863      zztmp = 1. / ( 1. - albo )
864      WHERE( ptsu(:,:,:) /= 0._wp )
865         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
866      ELSEWHERE
867         z1_st(:,:,:) = 0._wp
868      END WHERE
869      !                                     ! ========================== !
870      DO jl = 1, jpl                        !  Loop over ice categories  !
871         !                                  ! ========================== !
872         DO jj = 1 , jpj
873            DO ji = 1, jpi
874               ! ----------------------------!
875               !      I   Radiative FLUXES   !
876               ! ----------------------------!
877               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
878               ! Short Wave (sw)
879               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
880               ! Long  Wave (lw)
881               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
882               ! lw sensitivity
883               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
884
885               ! ----------------------------!
886               !     II    Turbulent FLUXES  !
887               ! ----------------------------!
888
889               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
890               ! Sensible Heat
891               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))
892               ! Latent Heat
893               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
894                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
895               ! Latent heat sensitivity for ice (Dqla/Dt)
896               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
897                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
898                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
899               ELSE
900                  dqla_ice(ji,jj,jl) = 0._wp
901               ENDIF
902
903               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
904               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
905
906               ! ----------------------------!
907               !     III    Total FLUXES     !
908               ! ----------------------------!
909               ! Downward Non Solar flux
910               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
911               ! Total non solar heat flux sensitivity for ice
912               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
913            END DO
914            !
915         END DO
916         !
917      END DO
918      !
919      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
920      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
921      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
922      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
923
924      ! --- evaporation --- !
925      z1_rLsub = 1._wp / rLsub
926      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
927      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
928      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
929
930      ! --- evaporation minus precipitation --- !
931      zsnw(:,:) = 0._wp
932      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
933      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
934      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
935      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
936
937      ! --- heat flux associated with emp --- !
938      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
939         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
940         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
941         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
942      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
943         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
944
945      ! --- total solar and non solar fluxes --- !
946      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
947         &           + qemp_ice(:,:) + qemp_oce(:,:)
948      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
949
950      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
951      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
952
953      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
954      DO jl = 1, jpl
955         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
956         !                         ! But we do not have Tice => consider it at 0degC => evap=0
957      END DO
958
959      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
960      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
961      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
962      !
963      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
964         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
965      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
966         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
967      ELSEWHERE                                                         ! zero when hs>0
968         qtr_ice_top(:,:,:) = 0._wp
969      END WHERE
970      !
971      IF(ln_ctl) THEN
972         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
973         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
974         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
975         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
976         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
977         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
978      ENDIF
979      !
980   END SUBROUTINE blk_ice_flx
981
982
983   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
984      !!---------------------------------------------------------------------
985      !!                     ***  ROUTINE blk_ice_qcn  ***
986      !!
987      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
988      !!                to force sea ice / snow thermodynamics
989      !!                in the case conduction flux is emulated
990      !!
991      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
992      !!                following the 0-layer Semtner (1976) approach
993      !!
994      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
995      !!              - qcn_ice : surface inner conduction flux (W/m2)
996      !!
997      !!---------------------------------------------------------------------
998      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
999      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
1000      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
1001      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
1002      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
1003      !
1004      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
1005      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
1006      !
1007      INTEGER  ::   ji, jj, jl           ! dummy loop indices
1008      INTEGER  ::   iter                 ! local integer
1009      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
1010      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
1011      REAL(wp) ::   zqc, zqnet           !
1012      REAL(wp) ::   zhe, zqa0            !
1013      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
1014      !!---------------------------------------------------------------------
1015
1016      ! -------------------------------------!
1017      !      I   Enhanced conduction factor  !
1018      ! -------------------------------------!
1019      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1020      ! Fichefet and Morales Maqueda, JGR 1997
1021      !
1022      zgfac(:,:,:) = 1._wp
1023
1024      IF( ld_virtual_itd ) THEN
1025         !
1026         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
1027         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1028         zfac3 = 2._wp / zepsilon
1029         !
1030         DO jl = 1, jpl
1031            DO jj = 1 , jpj
1032               DO ji = 1, jpi
1033                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
1034                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
1035               END DO
1036            END DO
1037         END DO
1038         !
1039      ENDIF
1040
1041      ! -------------------------------------------------------------!
1042      !      II   Surface temperature and conduction flux            !
1043      ! -------------------------------------------------------------!
1044      !
1045      zfac = rcnd_i * rn_cnd_s
1046      !
1047      DO jl = 1, jpl
1048         DO jj = 1 , jpj
1049            DO ji = 1, jpi
1050               !
1051               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1052                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1053               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1054               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1055               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1056               !
1057               DO iter = 1, nit     ! --- Iterative loop
1058                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1059                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1060                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1061               END DO
1062               !
1063               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1064               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1065               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1066               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) )  &
1067                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1068
1069               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1070               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)
1071
1072            END DO
1073         END DO
1074         !
1075      END DO
1076      !
1077   END SUBROUTINE blk_ice_qcn
1078
1079
1080   SUBROUTINE Cdn10_Lupkes2012( Cd )
1081      !!----------------------------------------------------------------------
1082      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1083      !!
1084      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1085      !!                 to make it dependent on edges at leads, melt ponds and flows.
1086      !!                 After some approximations, this can be resumed to a dependency
1087      !!                 on ice concentration.
1088      !!
1089      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1090      !!                 with the highest level of approximation: level4, eq.(59)
1091      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1092      !!
1093      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1094      !!
1095      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1096      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1097      !!                    A is the concentration of ice minus melt ponds (if any)
1098      !!
1099      !!                 This new drag has a parabolic shape (as a function of A) starting at
1100      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1101      !!                 and going down to Cdi(say 1.4e-3) for A=1
1102      !!
1103      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1104      !!                 => see Lupkes et al (2013)
1105      !!
1106      !! ** References : Lupkes et al. JGR 2012 (theory)
1107      !!                 Lupkes et al. GRL 2013 (application to GCM)
1108      !!
1109      !!----------------------------------------------------------------------
1110      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1111      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1112      REAL(wp), PARAMETER ::   znu   = 1._wp
1113      REAL(wp), PARAMETER ::   zmu   = 1._wp
1114      REAL(wp), PARAMETER ::   zbeta = 1._wp
1115      REAL(wp)            ::   zcoef
1116      !!----------------------------------------------------------------------
1117      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1118
1119      ! generic drag over a cell partly covered by ice
1120      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1121      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1122      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1123
1124      ! ice-atm drag
1125      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
1126         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1127
1128   END SUBROUTINE Cdn10_Lupkes2012
1129
1130
1131   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1132      !!----------------------------------------------------------------------
1133      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1134      !!
1135      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1136      !!                 between sea-ice and atmosphere with distinct momentum
1137      !!                 and heat coefficients depending on sea-ice concentration
1138      !!                 and atmospheric stability (no meltponds effect for now).
1139      !!
1140      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1141      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1142      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1143      !!                 to compute neutral transfert coefficients for both heat and
1144      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1145      !!                 coefficient is also taken into account following Louis (1979).
1146      !!
1147      !! ** References : Lupkes et al. JGR 2015 (theory)
1148      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1149      !!
1150      !!----------------------------------------------------------------------
1151      !
1152      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1153      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1154      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
1155      !
1156      ! ECHAM6 constants
1157      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1158      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1159      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1160      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1161      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1162      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1163      REAL(wp), PARAMETER ::   zc2          = zc * zc
1164      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1165      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1166      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1167      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1168      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1169      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1170      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1171      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1172      !
1173      INTEGER  ::   ji, jj         ! dummy loop indices
1174      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1175      REAL(wp) ::   zrib_o, zrib_i
1176      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1177      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1178      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1179      REAL(wp) ::   zCdn_form_tmp
1180      !!----------------------------------------------------------------------
1181
1182      ! mean temperature
1183      WHERE( at_i_b(:,:) > 1.e-20 )
1184         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1185      ELSEWHERE
1186         ztm_su(:,:) = rt0
1187      ENDWHERE
1188
1189      ! Momentum Neutral Transfert Coefficients (should be a constant)
1190      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1191      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1192      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1193      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1194
1195      ! Heat Neutral Transfert Coefficients
1196      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)
1197
1198      ! Atmospheric and Surface Variables
1199      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1200      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1201      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
1202      !
1203      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1204         DO ji = fs_2, fs_jpim1
1205            ! Virtual potential temperature [K]
1206            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1207            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1208            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1209
1210            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1211            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1212            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1213
1214            ! Momentum and Heat Neutral Transfert Coefficients
1215            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1216            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1217
1218            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1219            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1220            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1221            IF( zrib_o <= 0._wp ) THEN
1222               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
1223               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1224                  &             )**zgamma )**z1_gamma
1225            ELSE
1226               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1227               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1228            ENDIF
1229
1230            IF( zrib_i <= 0._wp ) THEN
1231               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1232               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1233            ELSE
1234               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1235               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1236            ENDIF
1237
1238            ! Momentum Transfert Coefficients (Eq. 38)
1239            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1240               &        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) )
1241
1242            ! Heat Transfert Coefficients (Eq. 49)
1243            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1244               &        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) )
1245            !
1246         END DO
1247      END DO
1248      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
1249      !
1250   END SUBROUTINE Cdn10_Lupkes2015
1251
1252#endif
1253
1254   !!======================================================================
1255END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.