New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk.F90 in NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/SBC/sbcblk.F90 @ 12342

Last change on this file since 12342 was 12342, checked in by acc, 2 years ago

Branch 2019/dev_r12072_MERGE_OPTION2_2019. Implement recommended changes to fix ticket #2370

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