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

Last change on this file since 12183 was 12183, checked in by mathiot, 9 months ago

minor changes on sbcblk to compile with Cray

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