source: NEMO/branches/2020/dev_r12702_ASINTER-02_emanuelaclementi_Waves/src/OCE/SBC/sbcblk.F90 @ 12991

Last change on this file since 12991 was 12991, checked in by emanuelaclementi, 5 months ago

Included wave-current processes following Couvelard et al 2019 and Gurvan code -ticket #2155

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