source: NEMO/trunk/src/OCE/SBC/sbcblk.F90 @ 13208

Last change on this file since 13208 was 13208, checked in by smasson, 19 months ago

trunk: Mid-year merge, merge back dev_r12472_ASINTER-05_Masson_CurrentFeedback

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