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

Last change on this file since 13501 was 13501, checked in by gsamson, 7 months ago

define default value (=0.) for sbc ABL hpgi and hpgj fields when 'not used'

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