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

Last change on this file since 11291 was 11291, checked in by laurent, 14 months ago

LB: more appropriated use of constant "emiss_w": emissivity of sea water.

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