source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk.F90 @ 13159

Last change on this file since 13159 was 13159, checked in by gsamson, 6 weeks ago

merge trunk@r13136 into ASINTER-06 branch; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

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