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

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

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

Last change on this file since 13305 was 13305, checked in by acc, 4 years ago

Trunk changes required to avoid issues with the outer halo in ORCA2_ICE_PISCES,
REPRO_8_4 tests with nn_hls=2. These changes ensure that tmask and output
from sbc_blk are set correctly in the outer halo. Failure to set valid
values in the outer halo can generate NaNs? which lead to OOB errors in the
XRGB lookup table used for the TRC optics.See #2366 for details. With these
changes all variants of the ORCA2_ICE_PISCES SETTE test will complete. There
are still differences between the 1 and 2 halo width runs but running with:
no land suppression; partial land suppression or full land suppression does
not alter either set of results. Likewise setting ln_nnogather either true
or false does not alter results. Differences in run.stat start after 140
timesteps and differences in tracer.stat start after 60 timesteps between
the different halo width sets. Equivalent tests with ln_icebergs = F show no
differences in run.stat but halo-width dependent differences in tracer.stat
persist (now after 64 timesteps).

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