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/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90 @ 12040

Last change on this file since 12040 was 12040, checked in by laurent, 4 years ago

Got rid of array "tsk" as it was useless.

  • Property svn:keywords set to Id
File size: 68.5 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   !!----------------------------------------------------------------------
21
22   !!----------------------------------------------------------------------
23   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
24   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
25   !!   blk_oce       : computes momentum, heat and freshwater fluxes over ocean
26   !!             sea-ice case only :
27   !!   blk_ice_tau   : provide the air-ice stress
28   !!   blk_ice_flx   : provide the heat and mass fluxes at air-ice interface
29   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
30   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
31   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
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 sbc_ice        ! Surface boundary condition: ice fields
42   USE lib_fortran    ! to use key_nosignedzero
43#if defined key_si3
44   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif
45   USE icethd_dh      ! for CALL ice_thd_snwblow
46#endif
47   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
48   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
49   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 201 + Edson et al. 2013)
50   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1)
51   !
52   USE iom            ! I/O manager library
53   USE in_out_manager ! I/O manager
54   USE lib_mpp        ! distribued memory computing library
55   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
56   USE prtctl         ! Print control
57
58   USE sbcblk_phy     ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc...
59
60
61   IMPLICIT NONE
62   PRIVATE
63
64   PUBLIC   sbc_blk_init  ! called in sbcmod
65   PUBLIC   sbc_blk       ! called in sbcmod
66#if defined key_si3
67   PUBLIC   blk_ice_tau   ! routine called in icesbc
68   PUBLIC   blk_ice_flx   ! routine called in icesbc
69   PUBLIC   blk_ice_qcn   ! routine called in icesbc
70#endif
71
72   INTEGER , PARAMETER ::   jpfld   =10           ! maximum number of files to read
73   INTEGER , PARAMETER ::   jp_wndi = 1           ! index of 10m wind velocity (i-component) (m/s)    at T-point
74   INTEGER , PARAMETER ::   jp_wndj = 2           ! index of 10m wind velocity (j-component) (m/s)    at T-point
75   INTEGER , PARAMETER ::   jp_tair = 3           ! index of 10m air temperature             (Kelvin)
76   INTEGER , PARAMETER ::   jp_humi = 4           ! index of specific humidity               ( % )
77   INTEGER , PARAMETER ::   jp_qsr  = 5           ! index of solar heat                      (W/m2)
78   INTEGER , PARAMETER ::   jp_qlw  = 6           ! index of Long wave                       (W/m2)
79   INTEGER , PARAMETER ::   jp_prec = 7           ! index of total precipitation (rain+snow) (Kg/m2/s)
80   INTEGER , PARAMETER ::   jp_snow = 8           ! index of snow (solid prcipitation)       (kg/m2/s)
81   INTEGER , PARAMETER ::   jp_slp  = 9           ! index of sea level pressure              (Pa)
82   INTEGER , PARAMETER ::   jp_tdif =10           ! index of tau diff associated to HF tau   (N/m2)   at T-point
83
84   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
85
86   !                           !!* Namelist namsbc_blk : bulk parameters
87   LOGICAL  ::   ln_NCAR        ! "NCAR"      algorithm   (Large and Yeager 2008)
88   LOGICAL  ::   ln_COARE_3p0   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
89   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013)
90   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1)
91   !
92   LOGICAL  ::   ln_taudif      ! logical flag to use the "mean of stress module - module of mean stress" data
93   REAL(wp) ::   rn_pfac        ! multiplication factor for precipitation
94   REAL(wp) ::   rn_efac        ! multiplication factor for evaporation
95   REAL(wp) ::   rn_vfac        ! multiplication factor for ice/ocean velocity in the calculation of wind stress
96   REAL(wp) ::   rn_zqt         ! z(q,t) : height of humidity and temperature measurements
97   REAL(wp) ::   rn_zu          ! z(u)   : height of wind measurements
98   !!gm ref namelist initialize it so remove the setting to false below
99   LOGICAL  ::   ln_Cd_L12 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2012)
100   LOGICAL  ::   ln_Cd_L15 = .FALSE. !  Modify the drag ice-atm depending on ice concentration (from Lupkes et al. JGR2015)
101   !
102   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Cd_atm                    ! transfer coefficient for momentum      (tau)
103   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ch_atm                    ! transfer coefficient for sensible heat (Q_sens)
104   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   Ce_atm                    ! tansfert coefficient for evaporation   (Q_lat)
105   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_zu                      ! air temperature at wind speed height (needed by Lupkes 2015 bulk scheme)
106   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   q_zu                      ! air spec. hum.  at wind speed height (needed by Lupkes 2015 bulk scheme)
107   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
108
109   !LB:
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   !
116   INTEGER  ::   nhumi          ! choice of the bulk algorithm
117   !                            ! associated indices:
118   INTEGER, PARAMETER :: np_humi_sph = 1
119   INTEGER, PARAMETER :: np_humi_dpt = 2
120   INTEGER, PARAMETER :: np_humi_rlh = 3
121   !LB.
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 "vectopt_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( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
144         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
145      !
146      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
147      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
148   END FUNCTION sbc_blk_alloc
149
150
151   SUBROUTINE sbc_blk_init
152      !!---------------------------------------------------------------------
153      !!                    ***  ROUTINE sbc_blk_init  ***
154      !!
155      !! ** Purpose :   choose and initialize a bulk formulae formulation
156      !!
157      !! ** Method  :
158      !!
159      !!----------------------------------------------------------------------
160      INTEGER  ::   ifpr, jfld            ! dummy loop indice and argument
161      INTEGER  ::   ios, ierror, ioptio   ! Local integer
162      !!
163      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
164      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
165      TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read
166      TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        "
167      TYPE(FLD_N) ::   sn_slp , sn_tdif                        !       "                        "
168      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
169         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
170         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
171         &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         &
172         &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           &
173         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh    ! cool-skin / warm-layer !LB
174      !!---------------------------------------------------------------------
175      !
176      !                                      ! allocate sbc_blk_core array
177      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
178      !
179      !                             !** read bulk namelist
180      REWIND( numnam_ref )                !* Namelist namsbc_blk in reference namelist : bulk parameters
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      REWIND( numnam_cfg )                !* Namelist namsbc_blk in configuration namelist : bulk parameters
185      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
186902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
187      !
188      IF(lwm) WRITE( numond, namsbc_blk )
189      !
190      !                             !** initialization of the chosen bulk formulae (+ check)
191      !                                   !* select the bulk chosen in the namelist and check the choice
192      ioptio = 0
193      IF( ln_NCAR      ) THEN
194         nblk =  np_NCAR        ;   ioptio = ioptio + 1
195      ENDIF
196      IF( ln_COARE_3p0 ) THEN
197         nblk =  np_COARE_3p0   ;   ioptio = ioptio + 1
198      ENDIF
199      IF( ln_COARE_3p6 ) THEN
200         nblk =  np_COARE_3p6   ;   ioptio = ioptio + 1
201      ENDIF
202      IF( ln_ECMWF     ) THEN
203         nblk =  np_ECMWF       ;   ioptio = ioptio + 1
204      ENDIF
205      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
206
207
208
209
210      !                             !** initialization of the cool-skin / warm-layer parametrization
211      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
212         !! Some namelist sanity tests:
213         IF( ln_NCAR )      &
214            & CALL ctl_stop( 'sbc_blk_init: Cool-skin/warm-layer param. not compatible with NCAR algorithm' )
215         IF( nn_fsbc /= 1 ) &
216            & CALL ctl_stop( 'sbc_blk_init: Please set "nn_fsbc" to 1 when using cool-skin/warm-layer param.')
217      END IF
218
219      ioptio = 0
220      IF( ln_humi_sph ) THEN
221         nhumi =  np_humi_sph    ;   ioptio = ioptio + 1
222      ENDIF
223      IF( ln_humi_dpt ) THEN
224         nhumi =  np_humi_dpt    ;   ioptio = ioptio + 1
225      ENDIF
226      IF( ln_humi_rlh ) THEN
227         nhumi =  np_humi_rlh    ;   ioptio = ioptio + 1
228      ENDIF
229      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one type of air humidity' )
230      !LB.
231
232      !
233      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
234         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
235         IF( sn_qsr%ln_tint ) THEN
236            CALL ctl_warn( 'sbc_blk_init: ln_dm2dc=T daily qsr time interpolation done by sbcdcy module',   &
237               &           '              ==> We force time interpolation = .false. for qsr' )
238            sn_qsr%ln_tint = .false.
239         ENDIF
240      ENDIF
241      !                                   !* set the bulk structure
242      !                                      !- store namelist information in an array
243      slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
244      slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw
245      slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi
246      slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow
247      slf_i(jp_slp)  = sn_slp    ;   slf_i(jp_tdif) = sn_tdif
248      !
249      lhftau = ln_taudif                     !- add an extra field if HF stress is used
250      jfld = jpfld - COUNT( (/.NOT.lhftau/) )
251      !
252      !                                      !- allocate the bulk structure
253      ALLOCATE( sf(jfld), STAT=ierror )
254      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
255      DO ifpr= 1, jfld
256         ALLOCATE( sf(ifpr)%fnow(jpi,jpj,1) )
257         IF( slf_i(ifpr)%ln_tint )   ALLOCATE( sf(ifpr)%fdta(jpi,jpj,1,2) )
258         IF( slf_i(ifpr)%freqh > 0. .AND. MOD( NINT(3600. * slf_i(ifpr)%freqh), nn_fsbc * NINT(rdt) ) /= 0 )   &
259            &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rdt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   &
260            &                 '               This is not ideal. You should consider changing either rdt or nn_fsbc value...' )
261
262      END DO
263      !                                      !- fill the bulk structure with namelist informations
264      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' )
265      !
266      IF( ln_wave ) THEN
267         !Activated wave module but neither drag nor stokes drift activated
268         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN
269            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
270            !drag coefficient read from wave model definable only with mfs bulk formulae and core
271         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN
272            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')
273         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN
274            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')
275         ENDIF
276      ELSE
277         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                &
278            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    &
279            &                  'with drag coefficient (ln_cdgw =T) '  ,                        &
280            &                  'or Stokes Drift (ln_sdw=T) ' ,                                 &
281            &                  'or ocean stress modification due to waves (ln_tauwoc=T) ',      &
282            &                  'or Stokes-Coriolis term (ln_stcori=T)'  )
283      ENDIF
284      !
285      !
286      IF(lwp) THEN                     !** Control print
287         !
288         WRITE(numout,*)                  !* namelist
289         WRITE(numout,*) '   Namelist namsbc_blk (other than data information):'
290         WRITE(numout,*) '      "NCAR"      algorithm   (Large and Yeager 2008)     ln_NCAR      = ', ln_NCAR
291         WRITE(numout,*) '      "COARE 3.0" algorithm   (Fairall et al. 2003)       ln_COARE_3p0 = ', ln_COARE_3p0
292         WRITE(numout,*) '      "COARE 3.6" algorithm (Fairall 2018 + Edson al 2013)ln_COARE_3p6 = ', ln_COARE_3p6
293         WRITE(numout,*) '      "ECMWF"     algorithm   (IFS cycle 45r1)            ln_ECMWF     = ', ln_ECMWF
294         WRITE(numout,*) '      add High freq.contribution to the stress module     ln_taudif    = ', ln_taudif
295         WRITE(numout,*) '      Air temperature and humidity reference height (m)   rn_zqt       = ', rn_zqt
296         WRITE(numout,*) '      Wind vector reference height (m)                    rn_zu        = ', rn_zu
297         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac
298         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac
299         WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac
300         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))'
301         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12
302         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15
303         !
304         WRITE(numout,*)
305         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
306         CASE( np_NCAR      )   ;   WRITE(numout,*) '   ==>>>   "NCAR" algorithm        (Large and Yeager 2008)'
307         CASE( np_COARE_3p0 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.0" algorithm   (Fairall et al. 2003)'
308         CASE( np_COARE_3p6 )   ;   WRITE(numout,*) '   ==>>>   "COARE 3.6" algorithm (Fairall 2018+Edson et al. 2013)'
309         CASE( np_ECMWF     )   ;   WRITE(numout,*) '   ==>>>   "ECMWF" algorithm       (IFS cycle 45r1)'
310         END SELECT
311         !
312         WRITE(numout,*)
313         WRITE(numout,*) '      use cool-skin  parameterization (SSST)  ln_skin_cs  = ', ln_skin_cs
314         WRITE(numout,*) '      use warm-layer parameterization (SSST)  ln_skin_wl  = ', ln_skin_wl
315         !
316         WRITE(numout,*)
317         SELECT CASE( nhumi )              !* Print the choice of air humidity
318         CASE( np_humi_sph )   ;   WRITE(numout,*) '   ==>>>   air humidity is SPECIFIC HUMIDITY     [kg/kg]'
319         CASE( np_humi_dpt )   ;   WRITE(numout,*) '   ==>>>   air humidity is DEW-POINT TEMPERATURE [K]'
320         CASE( np_humi_rlh )   ;   WRITE(numout,*) '   ==>>>   air humidity is RELATIVE HUMIDITY     [%]'
321         END SELECT
322         !
323      ENDIF
324      !
325   END SUBROUTINE sbc_blk_init
326
327
328   SUBROUTINE sbc_blk( kt )
329      !!---------------------------------------------------------------------
330      !!                    ***  ROUTINE sbc_blk  ***
331      !!
332      !! ** Purpose :   provide at each time step the surface ocean fluxes
333      !!              (momentum, heat, freshwater and runoff)
334      !!
335      !! ** Method  : (1) READ each fluxes in NetCDF files:
336      !!      the 10m wind velocity (i-component) (m/s)    at T-point
337      !!      the 10m wind velocity (j-component) (m/s)    at T-point
338      !!      the 10m or 2m specific humidity     ( % )
339      !!      the solar heat                      (W/m2)
340      !!      the Long wave                       (W/m2)
341      !!      the 10m or 2m air temperature       (Kelvin)
342      !!      the total precipitation (rain+snow) (Kg/m2/s)
343      !!      the snow (solid prcipitation)       (kg/m2/s)
344      !!      the tau diff associated to HF tau   (N/m2)   at T-point   (ln_taudif=T)
345      !!              (2) CALL blk_oce
346      !!
347      !!      C A U T I O N : never mask the surface stress fields
348      !!                      the stress is assumed to be in the (i,j) mesh referential
349      !!
350      !! ** Action  :   defined at each time-step at the air-sea interface
351      !!              - utau, vtau  i- and j-component of the wind stress
352      !!              - taum        wind stress module at T-point
353      !!              - wndm        wind speed  module at T-point over free ocean or leads in presence of sea-ice
354      !!              - qns, qsr    non-solar and solar heat fluxes
355      !!              - emp         upward mass flux (evapo. - precip.)
356      !!              - sfx         salt flux due to freezing/melting (non-zero only if ice is present)
357      !!
358      !! ** References :   Large & Yeager, 2004 / Large & Yeager, 2008
359      !!                   Brodeau et al. Ocean Modelling 2010
360      !!----------------------------------------------------------------------
361      INTEGER, INTENT(in) ::   kt   ! ocean time step
362      !!---------------------------------------------------------------------
363      !
364      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
365      !
366      !                                            ! compute the surface ocean fluxes using bulk formulea
367      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
368
369#if defined key_cice
370      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
371         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
372         IF( ln_dm2dc ) THEN
373            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
374         ELSE
375            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
376         ENDIF
377         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
378         SELECT CASE( nhumi )
379         CASE( np_humi_sph )
380            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1)
381         CASE( np_humi_dpt )
382            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
383         CASE( np_humi_rlh )
384            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
385         END SELECT
386         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
387         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
388         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
389         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
390      ENDIF
391#endif
392      !
393   END SUBROUTINE sbc_blk
394
395
396   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
397      !!---------------------------------------------------------------------
398      !!                     ***  ROUTINE blk_oce  ***
399      !!
400      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
401      !!      the ocean surface at each time step
402      !!
403      !! ** Method  :   bulk formulea for the ocean using atmospheric
404      !!      fields read in sbc_read
405      !!
406      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
407      !!              - vtau    : j-component of the stress at V-point  (N/m2)
408      !!              - taum    : Wind stress module at T-point         (N/m2)
409      !!              - wndm    : Wind speed module at T-point          (m/s)
410      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
411      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
412      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
413      !!
414      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
415      !!---------------------------------------------------------------------
416      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
417      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
418      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
419      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
420      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
421      !
422      INTEGER  ::   ji, jj               ! dummy loop indices
423      REAL(wp) ::   zztmp                ! local variable
424      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
425      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst  [kg/kg]
426      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
427      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
428      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
429      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
430      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
431      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg]
432      !!---------------------------------------------------------------------
433      !
434      ! local scalars ( place there for vector optimisation purposes)
435      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
436
437      ! ----------------------------------------------------------------------------- !
438      !      0   Wind components and module at T-point relative to the moving ocean   !
439      ! ----------------------------------------------------------------------------- !
440
441      ! ... components ( U10m - U_oce ) at T-point (unmasked)
442      !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
443      zwnd_i(:,:) = 0._wp
444      zwnd_j(:,:) = 0._wp
445#if defined key_cyclone
446      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
447      DO jj = 2, jpjm1
448         DO ji = fs_2, fs_jpim1   ! vect. opt.
449            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
450            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
451         END DO
452      END DO
453#endif
454      DO jj = 2, jpjm1
455         DO ji = fs_2, fs_jpim1   ! vect. opt.
456            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
457            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
458         END DO
459      END DO
460      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
461      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
462      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
463         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
464
465      ! ----------------------------------------------------------------------------- !
466      !      I   Solar FLUX                                                           !
467      ! ----------------------------------------------------------------------------- !
468
469      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
470      zztmp = 1. - albo
471      IF( ln_dm2dc ) THEN
472         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
473      ELSE
474         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
475      ENDIF
476
477
478      ! ----------------------------------------------------------------------------- !
479      !     II    Turbulent FLUXES                                                    !
480      ! ----------------------------------------------------------------------------- !
481
482      ! ... specific humidity at SST and IST tmask(
483      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
484
485      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
486         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
487         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
488      ENDIF
489
490      !LB:
491      ! zqair = specific humidity of air at "rn_zqt" m above thes sea:
492      SELECT CASE( nhumi )
493      CASE( np_humi_sph )
494         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
495      CASE( np_humi_dpt )
496         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
497         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
498      CASE( np_humi_rlh )
499         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
500         zqair(:,:) = 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
501      END SELECT
502      !LB.
503
504      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
505      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
506      !!    (since reanalysis products provide absolute temperature "T" at z, not theta !)
507      !!#LB: because AGRIF hates functions that return something else than a scalar, need to
508      !!     use vectorial version of gamma_moist() ...
509      DO jj = 2, jpjm1
510         DO ji = fs_2, fs_jpim1   ! vect. opt.
511            ztpot(ji,jj) = sf(jp_tair)%fnow(ji,jj,1) + gamma_moist( sf(jp_tair)%fnow(ji,jj,1), zqair(ji,jj) ) * rn_zqt
512         END DO
513      END DO
514      CALL lbc_lnk( 'sbcblk', ztpot, 'T',  1. ) ! for CTL-sanity of C1D cases...
515      !ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
516     
517      !! Time to call the user-selected bulk parameterization for
518      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more...
519      SELECT CASE( nblk )
520
521      CASE( np_NCAR      )
522         CALL turb_ncar     (     rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         &
523            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
524
525      CASE( np_COARE_3p0 )
526         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
527            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
528            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
529
530      CASE( np_COARE_3p6 )
531         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
532            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
533            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
534
535      CASE( np_ECMWF     )
536         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  &
537            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
538            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
539
540      CASE DEFAULT
541         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
542      END SELECT
543
544      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
545         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq:
546         WHERE ( fr_i(:,:) < 0.001_wp )
547            ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep them !!!
548            zst(:,:) = zst(:,:)*tmask(:,:,1)
549            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
550         ELSEWHERE
551            ! we forget about the update...
552            zst(:,:) = zqsb(:,:) !#LB: using what we backed up before skin-algo
553            zsq(:,:) = zqla(:,:) !#LB:  "   "   "
554         END WHERE
555      END IF
556
557      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
558      !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
559
560      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
561         !! If zu == zt, then ensuring once for all that:
562         t_zu(:,:) = ztpot(:,:)
563         q_zu(:,:) = zqair(:,:)
564      ENDIF
565
566
567      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
568      ! -------------------------------------------------------------
569
570      CALL BULK_FORMULA( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), &
571         &               Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:),           &
572         &               wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1),   &
573         &               taum(:,:), zqsb(:,:), zqla(:,:),                 &
574         &               pEvap=zevap(:,:), prhoa=rhoa(:,:) )
575
576      zqla(:,:)  =  zqla(:,:) * tmask(:,:,1)
577      zqsb(:,:)  =  zqsb(:,:) * tmask(:,:,1)
578      taum(:,:)  =  taum(:,:) * tmask(:,:,1)
579      zevap(:,:) = zevap(:,:) * tmask(:,:,1)
580
581      ! Tau i and j component on T-grid points, using array "Cd_atm" as a temporary array...
582      Cd_atm(:,:) = 0._wp
583      WHERE ( wndm(:,:) > 0._wp ) Cd_atm(:,:) = taum(:,:) / wndm(:,:)
584      zwnd_i(:,:) = Cd_atm(:,:) * zwnd_i(:,:)
585      zwnd_j(:,:) = Cd_atm(:,:) * zwnd_j(:,:)
586     
587      !                          ! add the HF tau contribution to the wind stress module
588      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
589
590      CALL iom_put( "taum_oce", taum )   ! output wind stress module
591
592      ! ... utau, vtau at U- and V_points, resp.
593      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
594      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
595      DO jj = 1, jpjm1
596         DO ji = 1, fs_jpim1
597            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
598               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
599            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
600               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
601         END DO
602      END DO
603      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
604
605
606      ! ----------------------------------------------------------------------------- !
607      !     III    Net longwave radiative FLUX                                        !
608      ! ----------------------------------------------------------------------------- !
609
610      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
611      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
612      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
613
614
615      IF(ln_ctl) THEN
616         CALL prt_ctl( tab2d_1=zevap , clinfo1=' blk_oce: zevap  : ') !LB
617         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
618         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
619         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
620         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
621         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
622            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
623         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
624         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
625      ENDIF
626
627      ! ----------------------------------------------------------------------------- !
628      !     IV    Total FLUXES                                                       !
629      ! ----------------------------------------------------------------------------- !
630      !
631      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
632         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
633      !
634      qns(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                &   ! Downward Non Solar
635         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
636         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
637         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
638         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
639         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
640         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
641      qns(:,:) = qns(:,:) * tmask(:,:,1)
642      !
643#if defined key_si3
644      qns_oce(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                ! non solar without emp (only needed by SI3)
645      qsr_oce(:,:) = qsr(:,:)
646#endif
647      !
648      CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )    ! output air density [kg/m^3]
649      CALL iom_put( "evap_oce" , zevap )                ! evaporation
650      CALL iom_put( "qlw_oce" ,  zqlw )                 ! output downward longwave heat over the ocean
651      CALL iom_put( "qsb_oce" ,  zqsb )                 ! output downward sensible heat over the ocean
652      CALL iom_put( "qla_oce" ,  zqla )                 ! output downward latent   heat over the ocean
653      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
654      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
655      CALL iom_put( 'snowpre', sprecip )                 ! Snow
656      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
657      !
658      IF ( nn_ice == 0 ) THEN
659         CALL iom_put( "qemp_oce" ,   qns-zqlw-zqsb-zqla )   ! output downward heat content of E-P over the ocean
660         CALL iom_put( "qns_oce"  ,   qns  )                 ! output downward non solar heat over the ocean
661         CALL iom_put( "qsr_oce"  ,   qsr  )                 ! output downward solar heat over the ocean
662         CALL iom_put( "qt_oce"   ,   qns+qsr )              ! output total downward heat over the ocean
663      ENDIF
664      !
665      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
666         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
667         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
668      ENDIF
669      !
670      IF(ln_ctl) THEN
671         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
672         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
673         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
674         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
675            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
676      ENDIF
677      !
678   END SUBROUTINE blk_oce
679
680
681
682#if defined key_si3
683   !!----------------------------------------------------------------------
684   !!   'key_si3'                                       SI3 sea-ice model
685   !!----------------------------------------------------------------------
686   !!   blk_ice_tau : provide the air-ice stress
687   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
688   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
689   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
690   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
691   !!----------------------------------------------------------------------
692
693   SUBROUTINE blk_ice_tau
694      !!---------------------------------------------------------------------
695      !!                     ***  ROUTINE blk_ice_tau  ***
696      !!
697      !! ** Purpose :   provide the surface boundary condition over sea-ice
698      !!
699      !! ** Method  :   compute momentum using bulk formulation
700      !!                formulea, ice variables and read atmospheric fields.
701      !!                NB: ice drag coefficient is assumed to be a constant
702      !!---------------------------------------------------------------------
703      INTEGER  ::   ji, jj    ! dummy loop indices
704      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
705      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
706      !!---------------------------------------------------------------------
707      !
708      ! set transfer coefficients to default sea-ice values
709      Cd_atm(:,:) = rCd_ice
710      Ch_atm(:,:) = rCd_ice
711      Ce_atm(:,:) = rCd_ice
712
713      wndm_ice(:,:) = 0._wp      !!gm brutal....
714
715      ! ------------------------------------------------------------ !
716      !    Wind module relative to the moving ice ( U10m - U_ice )   !
717      ! ------------------------------------------------------------ !
718      ! C-grid ice dynamics :   U & V-points (same as ocean)
719      DO jj = 2, jpjm1
720         DO ji = fs_2, fs_jpim1   ! vect. opt.
721            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
722            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
723            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
724         END DO
725      END DO
726      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
727      !
728      ! Make ice-atm. drag dependent on ice concentration
729      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
730         CALL Cdn10_Lupkes2012( Cd_atm )
731         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
732      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
733         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
734      ENDIF
735
736      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
737      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
738
739      ! local scalars ( place there for vector optimisation purposes)
740
741      !!gm brutal....
742      utau_ice  (:,:) = 0._wp
743      vtau_ice  (:,:) = 0._wp
744      !!gm end
745
746      ! ------------------------------------------------------------ !
747      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
748      ! ------------------------------------------------------------ !
749      ! C-grid ice dynamics :   U & V-points (same as ocean)
750      DO jj = 2, jpjm1
751         DO ji = fs_2, fs_jpim1   ! vect. opt.
752            utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
753               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
754            vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
755               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
756         END DO
757      END DO
758      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
759      !
760      !
761      IF(ln_ctl) THEN
762         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
763         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
764      ENDIF
765      !
766   END SUBROUTINE blk_ice_tau
767
768
769   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
770      !!---------------------------------------------------------------------
771      !!                     ***  ROUTINE blk_ice_flx  ***
772      !!
773      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
774      !!
775      !! ** Method  :   compute heat and freshwater exchanged
776      !!                between atmosphere and sea-ice using bulk formulation
777      !!                formulea, ice variables and read atmmospheric fields.
778      !!
779      !! caution : the net upward water flux has with mm/day unit
780      !!---------------------------------------------------------------------
781      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
782      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
783      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
784      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
785      !!
786      INTEGER  ::   ji, jj, jl               ! dummy loop indices
787      REAL(wp) ::   zst3                     ! local variable
788      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
789      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
790      REAL(wp) ::   zfr1, zfr2               ! local variables
791      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
792      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
793      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
794      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
795      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
796      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
797      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
798      !!---------------------------------------------------------------------
799      !
800      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
801      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
802      !
803
804      !LB:
805      SELECT CASE( nhumi )
806      CASE( np_humi_sph )
807         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
808      CASE( np_humi_dpt )
809         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
810         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
811      CASE( np_humi_rlh )
812         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
813         zqair(:,:) = 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
814      END SELECT
815      !LB.
816
817      zztmp = 1. / ( 1. - albo )
818      WHERE( ptsu(:,:,:) /= 0._wp )
819         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
820      ELSEWHERE
821         z1_st(:,:,:) = 0._wp
822      END WHERE
823      !                                     ! ========================== !
824      DO jl = 1, jpl                        !  Loop over ice categories  !
825         !                                  ! ========================== !
826         DO jj = 1 , jpj
827            DO ji = 1, jpi
828               ! ----------------------------!
829               !      I   Radiative FLUXES   !
830               ! ----------------------------!
831               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
832               ! Short Wave (sw)
833               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
834               ! Long  Wave (lw)
835               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
836               ! lw sensitivity
837               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
838
839               ! ----------------------------!
840               !     II    Turbulent FLUXES  !
841               ! ----------------------------!
842
843               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
844               ! Sensible Heat
845               z_qsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj) * (ptsu(ji,jj,jl) - sf(jp_tair)%fnow(ji,jj,1))
846               ! Latent Heat
847               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
848                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
849               ! Latent heat sensitivity for ice (Dqla/Dt)
850               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
851                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
852                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
853               ELSE
854                  dqla_ice(ji,jj,jl) = 0._wp
855               ENDIF
856
857               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
858               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
859
860               ! ----------------------------!
861               !     III    Total FLUXES     !
862               ! ----------------------------!
863               ! Downward Non Solar flux
864               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
865               ! Total non solar heat flux sensitivity for ice
866               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
867            END DO
868            !
869         END DO
870         !
871      END DO
872      !
873      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
874      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
875      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
876      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
877
878      ! --- evaporation --- !
879      z1_rLsub = 1._wp / rLsub
880      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
881      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
882      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
883
884      ! --- evaporation minus precipitation --- !
885      zsnw(:,:) = 0._wp
886      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
887      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
888      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
889      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
890
891      ! --- heat flux associated with emp --- !
892      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
893         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
894         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
895         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
896      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
897         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
898
899      ! --- total solar and non solar fluxes --- !
900      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
901         &           + qemp_ice(:,:) + qemp_oce(:,:)
902      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
903
904      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
905      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
906
907      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
908      DO jl = 1, jpl
909         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
910         !                         ! But we do not have Tice => consider it at 0degC => evap=0
911      END DO
912
913      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
914      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
915      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
916      !
917      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
918         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
919      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
920         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
921      ELSEWHERE                                                         ! zero when hs>0
922         qtr_ice_top(:,:,:) = 0._wp
923      END WHERE
924      !
925      IF(ln_ctl) THEN
926         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
927         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
928         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
929         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
930         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
931         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
932      ENDIF
933      !
934   END SUBROUTINE blk_ice_flx
935
936
937   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
938      !!---------------------------------------------------------------------
939      !!                     ***  ROUTINE blk_ice_qcn  ***
940      !!
941      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
942      !!                to force sea ice / snow thermodynamics
943      !!                in the case conduction flux is emulated
944      !!
945      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
946      !!                following the 0-layer Semtner (1976) approach
947      !!
948      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
949      !!              - qcn_ice : surface inner conduction flux (W/m2)
950      !!
951      !!---------------------------------------------------------------------
952      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
953      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
954      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
955      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
956      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
957      !
958      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
959      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
960      !
961      INTEGER  ::   ji, jj, jl           ! dummy loop indices
962      INTEGER  ::   iter                 ! local integer
963      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
964      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
965      REAL(wp) ::   zqc, zqnet           !
966      REAL(wp) ::   zhe, zqa0            !
967      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
968      !!---------------------------------------------------------------------
969
970      ! -------------------------------------!
971      !      I   Enhanced conduction factor  !
972      ! -------------------------------------!
973      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
974      ! Fichefet and Morales Maqueda, JGR 1997
975      !
976      zgfac(:,:,:) = 1._wp
977
978      IF( ld_virtual_itd ) THEN
979         !
980         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
981         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
982         zfac3 = 2._wp / zepsilon
983         !
984         DO jl = 1, jpl
985            DO jj = 1 , jpj
986               DO ji = 1, jpi
987                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
988                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
989               END DO
990            END DO
991         END DO
992         !
993      ENDIF
994
995      ! -------------------------------------------------------------!
996      !      II   Surface temperature and conduction flux            !
997      ! -------------------------------------------------------------!
998      !
999      zfac = rcnd_i * rn_cnd_s
1000      !
1001      DO jl = 1, jpl
1002         DO jj = 1 , jpj
1003            DO ji = 1, jpi
1004               !
1005               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1006                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1007               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1008               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1009               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1010               !
1011               DO iter = 1, nit     ! --- Iterative loop
1012                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1013                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1014                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1015               END DO
1016               !
1017               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1018               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1019               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1020               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) )  &
1021                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1022
1023               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1024               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)
1025
1026            END DO
1027         END DO
1028         !
1029      END DO
1030      !
1031   END SUBROUTINE blk_ice_qcn
1032
1033
1034   SUBROUTINE Cdn10_Lupkes2012( Cd )
1035      !!----------------------------------------------------------------------
1036      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1037      !!
1038      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1039      !!                 to make it dependent on edges at leads, melt ponds and flows.
1040      !!                 After some approximations, this can be resumed to a dependency
1041      !!                 on ice concentration.
1042      !!
1043      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1044      !!                 with the highest level of approximation: level4, eq.(59)
1045      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1046      !!
1047      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1048      !!
1049      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1050      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1051      !!                    A is the concentration of ice minus melt ponds (if any)
1052      !!
1053      !!                 This new drag has a parabolic shape (as a function of A) starting at
1054      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1055      !!                 and going down to Cdi(say 1.4e-3) for A=1
1056      !!
1057      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1058      !!                 => see Lupkes et al (2013)
1059      !!
1060      !! ** References : Lupkes et al. JGR 2012 (theory)
1061      !!                 Lupkes et al. GRL 2013 (application to GCM)
1062      !!
1063      !!----------------------------------------------------------------------
1064      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1065      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1066      REAL(wp), PARAMETER ::   znu   = 1._wp
1067      REAL(wp), PARAMETER ::   zmu   = 1._wp
1068      REAL(wp), PARAMETER ::   zbeta = 1._wp
1069      REAL(wp)            ::   zcoef
1070      !!----------------------------------------------------------------------
1071      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1072
1073      ! generic drag over a cell partly covered by ice
1074      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1075      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1076      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1077
1078      ! ice-atm drag
1079      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
1080         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1081
1082   END SUBROUTINE Cdn10_Lupkes2012
1083
1084
1085   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1086      !!----------------------------------------------------------------------
1087      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1088      !!
1089      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1090      !!                 between sea-ice and atmosphere with distinct momentum
1091      !!                 and heat coefficients depending on sea-ice concentration
1092      !!                 and atmospheric stability (no meltponds effect for now).
1093      !!
1094      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1095      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1096      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1097      !!                 to compute neutral transfert coefficients for both heat and
1098      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1099      !!                 coefficient is also taken into account following Louis (1979).
1100      !!
1101      !! ** References : Lupkes et al. JGR 2015 (theory)
1102      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1103      !!
1104      !!----------------------------------------------------------------------
1105      !
1106      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1107      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1108      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
1109      !
1110      ! ECHAM6 constants
1111      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1112      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1113      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1114      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1115      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1116      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1117      REAL(wp), PARAMETER ::   zc2          = zc * zc
1118      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1119      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1120      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1121      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1122      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1123      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1124      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1125      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1126      !
1127      INTEGER  ::   ji, jj         ! dummy loop indices
1128      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1129      REAL(wp) ::   zrib_o, zrib_i
1130      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1131      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1132      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1133      REAL(wp) ::   zCdn_form_tmp
1134      !!----------------------------------------------------------------------
1135
1136      ! mean temperature
1137      WHERE( at_i_b(:,:) > 1.e-20 )
1138         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1139      ELSEWHERE
1140         ztm_su(:,:) = rt0
1141      ENDWHERE
1142
1143      ! Momentum Neutral Transfert Coefficients (should be a constant)
1144      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1145      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1146      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1147      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1148
1149      ! Heat Neutral Transfert Coefficients
1150      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 (cf Lupkes email for details)
1151
1152      ! Atmospheric and Surface Variables
1153      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1154      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1155      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
1156      !
1157      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1158         DO ji = fs_2, fs_jpim1
1159            ! Virtual potential temperature [K]
1160            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1161            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1162            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1163
1164            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1165            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1166            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1167
1168            ! Momentum and Heat Neutral Transfert Coefficients
1169            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1170            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1171
1172            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1173            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1174            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1175            IF( zrib_o <= 0._wp ) THEN
1176               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
1177               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1178                  &             )**zgamma )**z1_gamma
1179            ELSE
1180               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1181               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1182            ENDIF
1183
1184            IF( zrib_i <= 0._wp ) THEN
1185               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1186               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1187            ELSE
1188               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1189               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1190            ENDIF
1191
1192            ! Momentum Transfert Coefficients (Eq. 38)
1193            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1194               &        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) )
1195
1196            ! Heat Transfert Coefficients (Eq. 49)
1197            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1198               &        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) )
1199            !
1200         END DO
1201      END DO
1202      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
1203      !
1204   END SUBROUTINE Cdn10_Lupkes2015
1205
1206#endif
1207
1208   !!======================================================================
1209END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.