source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk.F90 @ 13135

Last change on this file since 13135 was 13135, checked in by orioltp, 3 months ago

dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation: merge with trunk@13134, see #2364

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