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

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

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk.F90 @ 13830

Last change on this file since 13830 was 13830, checked in by laurent, 3 years ago

Sea-ice bulk algorithm "AN05" (Andreas et al. 2005) now operational!

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