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

Last change on this file since 12038 was 12038, checked in by laurent, 10 months ago

Tiny bug fix…

  • Property svn:keywords set to Id
File size: 68.7 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), tsk (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      IF( kt == nit000 ) tsk(:,:) = sst_m(:,:)*tmask(:,:,1)  ! no previous estimate of skin temperature => using bulk SST (use restart?)
367      !                                                      ! compute the surface ocean fluxes using bulk formulea
368      IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
369
370#if defined key_cice
371      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
372         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
373         IF( ln_dm2dc ) THEN
374            qsr_ice(:,:,1) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )
375         ELSE
376            qsr_ice(:,:,1) =          sf(jp_qsr)%fnow(:,:,1)
377         ENDIF
378         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
379         SELECT CASE( nhumi )
380         CASE( np_humi_sph )
381            qatm_ice(:,:) =           sf(jp_humi)%fnow(:,:,1)
382         CASE( np_humi_dpt )
383            qatm_ice(:,:) = q_sat(    sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
384         CASE( np_humi_rlh )
385            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
386         END SELECT
387         tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
388         sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
389         wndi_ice(:,:)    = sf(jp_wndi)%fnow(:,:,1)
390         wndj_ice(:,:)    = sf(jp_wndj)%fnow(:,:,1)
391      ENDIF
392#endif
393      !
394   END SUBROUTINE sbc_blk
395
396
397   SUBROUTINE blk_oce( kt, sf, pst, pu, pv )
398      !!---------------------------------------------------------------------
399      !!                     ***  ROUTINE blk_oce  ***
400      !!
401      !! ** Purpose :   provide the momentum, heat and freshwater fluxes at
402      !!      the ocean surface at each time step
403      !!
404      !! ** Method  :   bulk formulea for the ocean using atmospheric
405      !!      fields read in sbc_read
406      !!
407      !! ** Outputs : - utau    : i-component of the stress at U-point  (N/m2)
408      !!              - vtau    : j-component of the stress at V-point  (N/m2)
409      !!              - taum    : Wind stress module at T-point         (N/m2)
410      !!              - wndm    : Wind speed module at T-point          (m/s)
411      !!              - qsr     : Solar heat flux over the ocean        (W/m2)
412      !!              - qns     : Non Solar heat flux over the ocean    (W/m2)
413      !!              - emp     : evaporation minus precipitation       (kg/m2/s)
414      !!
415      !!  ** Nota  :   sf has to be a dummy argument for AGRIF on NEC
416      !!---------------------------------------------------------------------
417      INTEGER  , INTENT(in   )                 ::   kt    ! time step index
418      TYPE(fld), INTENT(inout), DIMENSION(:)   ::   sf    ! input data
419      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pst   ! surface temperature                      [Celcius]
420      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pu    ! surface current at U-point (i-component) [m/s]
421      REAL(wp) , INTENT(in)   , DIMENSION(:,:) ::   pv    ! surface current at V-point (j-component) [m/s]
422      !
423      INTEGER  ::   ji, jj               ! dummy loop indices
424      REAL(wp) ::   zztmp                ! local variable
425      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
426      REAL(wp), DIMENSION(jpi,jpj) ::   zsq               ! specific humidity at pst  [kg/kg]
427      REAL(wp), DIMENSION(jpi,jpj) ::   zqlw, zqsb        ! long wave and sensible heat fluxes
428      REAL(wp), DIMENSION(jpi,jpj) ::   zqla, zevap       ! latent heat fluxes and evaporation
429      REAL(wp), DIMENSION(jpi,jpj) ::   zst               ! surface temperature in Kelvin
430      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s]
431      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K]
432      REAL(wp), DIMENSION(jpi,jpj) ::   zqair             ! specific humidity     of air at z=rn_zqt [kg/kg]
433      !!---------------------------------------------------------------------
434      !
435      ! local scalars ( place there for vector optimisation purposes)
436      zst(:,:) = pst(:,:) + rt0      ! convert SST from Celcius to Kelvin (and set minimum value far above 0 K)
437
438      ! ----------------------------------------------------------------------------- !
439      !      0   Wind components and module at T-point relative to the moving ocean   !
440      ! ----------------------------------------------------------------------------- !
441
442      ! ... components ( U10m - U_oce ) at T-point (unmasked)
443      !!gm    move zwnd_i (_j) set to zero  inside the key_cyclone ???
444      zwnd_i(:,:) = 0._wp
445      zwnd_j(:,:) = 0._wp
446#if defined key_cyclone
447      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
448      DO jj = 2, jpjm1
449         DO ji = fs_2, fs_jpim1   ! vect. opt.
450            sf(jp_wndi)%fnow(ji,jj,1) = sf(jp_wndi)%fnow(ji,jj,1) + zwnd_i(ji,jj)
451            sf(jp_wndj)%fnow(ji,jj,1) = sf(jp_wndj)%fnow(ji,jj,1) + zwnd_j(ji,jj)
452         END DO
453      END DO
454#endif
455      DO jj = 2, jpjm1
456         DO ji = fs_2, fs_jpim1   ! vect. opt.
457            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
458            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
459         END DO
460      END DO
461      CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. )
462      ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked)
463      wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   &
464         &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1)
465
466      ! ----------------------------------------------------------------------------- !
467      !      I   Solar FLUX                                                           !
468      ! ----------------------------------------------------------------------------- !
469
470      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
471      zztmp = 1. - albo
472      IF( ln_dm2dc ) THEN
473         qsr(:,:) = zztmp * sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(:,:,1)
474      ELSE
475         qsr(:,:) = zztmp *          sf(jp_qsr)%fnow(:,:,1)   * tmask(:,:,1)
476      ENDIF
477
478
479      ! ----------------------------------------------------------------------------- !
480      !     II    Turbulent FLUXES                                                    !
481      ! ----------------------------------------------------------------------------- !
482
483      ! ... specific humidity at SST and IST tmask(
484      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) )
485
486      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
487         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action
488         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action
489      ENDIF
490
491      !LB:
492      ! zqair = specific humidity of air at "rn_zqt" m above thes sea:
493      SELECT CASE( nhumi )
494      CASE( np_humi_sph )
495         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
496      CASE( np_humi_dpt )
497         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm
498         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
499      CASE( np_humi_rlh )
500         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm
501         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
502      END SELECT
503      !LB.
504
505      !! Estimate of potential temperature at z=rn_zqt, based on adiabatic lapse-rate
506      !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
507      !!    (since reanalysis products provide absolute temperature "T" at z, not theta !)
508      !!#LB: because AGRIF hates functions that return something else than a scalar, need to
509      !!     use vectorial version of gamma_moist() ...
510      DO jj = 2, jpjm1
511         DO ji = fs_2, fs_jpim1   ! vect. opt.
512            ztpot(ji,jj) = sf(jp_tair)%fnow(ji,jj,1) + gamma_moist( sf(jp_tair)%fnow(ji,jj,1), zqair(ji,jj) ) * rn_zqt
513         END DO
514      END DO
515      CALL lbc_lnk( 'sbcblk', ztpot, 'T',  1. ) ! for CTL-sanity of C1D cases...
516      !ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), zqair(:,:) ) * rn_zqt
517     
518      !! Time to call the user-selected bulk parameterization for
519      !!  ==  transfer coefficients  ==!   Cd, Ch, Ce at T-point, and more...
520      SELECT CASE( nblk )
521
522      CASE( np_NCAR      )
523         CALL turb_ncar     (     rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm,                         &
524            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
525
526      CASE( np_COARE_3p0 )
527         CALL turb_coare3p0 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
528            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
529            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
530
531      CASE( np_COARE_3p6 )
532         CALL turb_coare3p6 ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl, &
533            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
534            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
535
536      CASE( np_ECMWF     )
537         CALL turb_ecmwf   ( kt, rn_zqt, rn_zu, zst, ztpot, zsq, zqair, wndm, ln_skin_cs, ln_skin_wl,  &
538            &                Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce,     &
539            &                Qsw=qsr(:,:), rad_lw=sf(jp_qlw)%fnow(:,:,1), slp=sf(jp_slp)%fnow(:,:,1) )
540
541      CASE DEFAULT
542         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
543      END SELECT
544
545      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
546         !! In the presence of sea-ice we forget about the cool-skin/warm-layer update of zst and zsq:
547         WHERE ( fr_i(:,:) < 0.001_wp )
548            ! zst and zsq have been updated by cool-skin/warm-layer scheme and we keep them !!!
549            zst(:,:) = zst(:,:)*tmask(:,:,1)
550            zsq(:,:) = zsq(:,:)*tmask(:,:,1)
551         ELSEWHERE
552            ! we forget about the update...
553            zst(:,:) = zqsb(:,:) !#LB: using what we backed up before skin-algo
554            zsq(:,:) = zqla(:,:) !#LB:  "   "   "
555         END WHERE
556         tsk(:,:) = zst(:,:) !#LB: Update of tsk, the "official" array for skin temperature
557      END IF
558
559      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef.
560      !!      CALL iom_put( "Ch_oce", Ch_atm)  ! output value of pure ocean-atm. transfer coef.
561
562      IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN
563         !! If zu == zt, then ensuring once for all that:
564         t_zu(:,:) = ztpot(:,:)
565         q_zu(:,:) = zqair(:,:)
566      ENDIF
567
568
569      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
570      ! -------------------------------------------------------------
571
572      CALL BULK_FORMULA( rn_zu, zst(:,:), zsq(:,:), t_zu(:,:), q_zu(:,:), &
573         &               Cd_atm(:,:), Ch_atm(:,:), Ce_atm(:,:),           &
574         &               wndm(:,:), zU_zu(:,:), sf(jp_slp)%fnow(:,:,1),   &
575         &               taum(:,:), zqsb(:,:), zqla(:,:),                 &
576         &               pEvap=zevap(:,:), prhoa=rhoa(:,:) )
577
578      zqla(:,:)  =  zqla(:,:) * tmask(:,:,1)
579      zqsb(:,:)  =  zqsb(:,:) * tmask(:,:,1)
580      taum(:,:)  =  taum(:,:) * tmask(:,:,1)
581      zevap(:,:) = zevap(:,:) * tmask(:,:,1)
582
583      ! Tau i and j component on T-grid points, using array "Cd_atm" as a temporary array...
584      Cd_atm(:,:) = 0._wp
585      WHERE ( wndm(:,:) > 0._wp ) Cd_atm(:,:) = taum(:,:) / wndm(:,:)
586      zwnd_i(:,:) = Cd_atm(:,:) * zwnd_i(:,:)
587      zwnd_j(:,:) = Cd_atm(:,:) * zwnd_j(:,:)
588     
589      !                          ! add the HF tau contribution to the wind stress module
590      IF( lhftau )   taum(:,:) = taum(:,:) + sf(jp_tdif)%fnow(:,:,1)
591
592      CALL iom_put( "taum_oce", taum )   ! output wind stress module
593
594      ! ... utau, vtau at U- and V_points, resp.
595      !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines
596      !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves
597      DO jj = 1, jpjm1
598         DO ji = 1, fs_jpim1
599            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) &
600               &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1))
601            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) &
602               &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1))
603         END DO
604      END DO
605      CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. )
606
607
608      ! ----------------------------------------------------------------------------- !
609      !     III    Net longwave radiative FLUX                                        !
610      ! ----------------------------------------------------------------------------- !
611
612      !! LB: now moved after Turbulent fluxes because must use the skin temperature rather that the SST
613      !! (zst is skin temperature if ln_skin_cs==.TRUE. .OR. ln_skin_wl==.TRUE.)
614      zqlw(:,:) = emiss_w * ( sf(jp_qlw)%fnow(:,:,1) - stefan*zst(:,:)*zst(:,:)*zst(:,:)*zst(:,:) ) * tmask(:,:,1)   ! Net radiative longwave flux
615
616
617      IF(ln_ctl) THEN
618         CALL prt_ctl( tab2d_1=zevap , clinfo1=' blk_oce: zevap  : ') !LB
619         CALL prt_ctl( tab2d_1=zqla  , clinfo1=' blk_oce: zqla   : ', tab2d_2=Ce_atm , clinfo2=' Ce_oce  : ' )
620         CALL prt_ctl( tab2d_1=zqsb  , clinfo1=' blk_oce: zqsb   : ', tab2d_2=Ch_atm , clinfo2=' Ch_oce  : ' )
621         CALL prt_ctl( tab2d_1=zqlw  , clinfo1=' blk_oce: zqlw   : ', tab2d_2=qsr, clinfo2=' qsr : ' )
622         CALL prt_ctl( tab2d_1=zsq   , clinfo1=' blk_oce: zsq    : ', tab2d_2=zst, clinfo2=' zst : ' )
623         CALL prt_ctl( tab2d_1=utau  , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
624            &          tab2d_2=vtau  , clinfo2=           ' vtau : ', mask2=vmask )
625         CALL prt_ctl( tab2d_1=wndm  , clinfo1=' blk_oce: wndm   : ')
626         CALL prt_ctl( tab2d_1=zst   , clinfo1=' blk_oce: zst    : ')
627      ENDIF
628
629      ! ----------------------------------------------------------------------------- !
630      !     IV    Total FLUXES                                                       !
631      ! ----------------------------------------------------------------------------- !
632      !
633      emp (:,:) = (  zevap(:,:)                                          &   ! mass flux (evap. - precip.)
634         &         - sf(jp_prec)%fnow(:,:,1) * rn_pfac  ) * tmask(:,:,1)
635      !
636      qns(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                &   ! Downward Non Solar
637         &     - sf(jp_snow)%fnow(:,:,1) * rn_pfac * rLfus                        &   ! remove latent melting heat for solid precip
638         &     - zevap(:,:) * pst(:,:) * rcp                                      &   ! remove evap heat content at SST !LB??? pst is Celsius !?
639         &     + ( sf(jp_prec)%fnow(:,:,1) - sf(jp_snow)%fnow(:,:,1) ) * rn_pfac  &   ! add liquid precip heat content at Tair
640         &     * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp                          &
641         &     + sf(jp_snow)%fnow(:,:,1) * rn_pfac                                &   ! add solid  precip heat content at min(Tair,Tsnow)
642         &     * ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi
643      qns(:,:) = qns(:,:) * tmask(:,:,1)
644      !
645#if defined key_si3
646      qns_oce(:,:) = zqlw(:,:) + zqsb(:,:) + zqla(:,:)                                ! non solar without emp (only needed by SI3)
647      qsr_oce(:,:) = qsr(:,:)
648#endif
649      !
650      CALL iom_put( "rho_air"  , rhoa*tmask(:,:,1) )    ! output air density [kg/m^3]
651      CALL iom_put( "evap_oce" , zevap )                ! evaporation
652      CALL iom_put( "qlw_oce" ,  zqlw )                 ! output downward longwave heat over the ocean
653      CALL iom_put( "qsb_oce" ,  zqsb )                 ! output downward sensible heat over the ocean
654      CALL iom_put( "qla_oce" ,  zqla )                 ! output downward latent   heat over the ocean
655      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output total precipitation [kg/m2/s]
656      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1) ! output solid precipitation [kg/m2/s]
657      CALL iom_put( 'snowpre', sprecip )                 ! Snow
658      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation
659      !
660      IF ( nn_ice == 0 ) THEN
661         CALL iom_put( "qemp_oce" ,   qns-zqlw-zqsb-zqla )   ! output downward heat content of E-P over the ocean
662         CALL iom_put( "qns_oce"  ,   qns  )                 ! output downward non solar heat over the ocean
663         CALL iom_put( "qsr_oce"  ,   qsr  )                 ! output downward solar heat over the ocean
664         CALL iom_put( "qt_oce"   ,   qns+qsr )              ! output total downward heat over the ocean
665      ENDIF
666      !
667      IF( ln_skin_cs .OR. ln_skin_wl ) THEN
668         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius
669         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference...
670      ENDIF
671      !
672      IF(ln_ctl) THEN
673         CALL prt_ctl(tab2d_1=zqsb , clinfo1=' blk_oce: zqsb   : ', tab2d_2=zqlw , clinfo2=' zqlw  : ')
674         CALL prt_ctl(tab2d_1=zqla , clinfo1=' blk_oce: zqla   : ', tab2d_2=qsr  , clinfo2=' qsr   : ')
675         CALL prt_ctl(tab2d_1=pst  , clinfo1=' blk_oce: pst    : ', tab2d_2=emp  , clinfo2=' emp   : ')
676         CALL prt_ctl(tab2d_1=utau , clinfo1=' blk_oce: utau   : ', mask1=umask,   &
677            &         tab2d_2=vtau , clinfo2=              ' vtau  : ' , mask2=vmask )
678      ENDIF
679      !
680   END SUBROUTINE blk_oce
681
682
683
684#if defined key_si3
685   !!----------------------------------------------------------------------
686   !!   'key_si3'                                       SI3 sea-ice model
687   !!----------------------------------------------------------------------
688   !!   blk_ice_tau : provide the air-ice stress
689   !!   blk_ice_flx : provide the heat and mass fluxes at air-ice interface
690   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
691   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
692   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
693   !!----------------------------------------------------------------------
694
695   SUBROUTINE blk_ice_tau
696      !!---------------------------------------------------------------------
697      !!                     ***  ROUTINE blk_ice_tau  ***
698      !!
699      !! ** Purpose :   provide the surface boundary condition over sea-ice
700      !!
701      !! ** Method  :   compute momentum using bulk formulation
702      !!                formulea, ice variables and read atmospheric fields.
703      !!                NB: ice drag coefficient is assumed to be a constant
704      !!---------------------------------------------------------------------
705      INTEGER  ::   ji, jj    ! dummy loop indices
706      REAL(wp) ::   zwndi_f , zwndj_f, zwnorm_f   ! relative wind module and components at F-point
707      REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point
708      !!---------------------------------------------------------------------
709      !
710      ! set transfer coefficients to default sea-ice values
711      Cd_atm(:,:) = rCd_ice
712      Ch_atm(:,:) = rCd_ice
713      Ce_atm(:,:) = rCd_ice
714
715      wndm_ice(:,:) = 0._wp      !!gm brutal....
716
717      ! ------------------------------------------------------------ !
718      !    Wind module relative to the moving ice ( U10m - U_ice )   !
719      ! ------------------------------------------------------------ !
720      ! C-grid ice dynamics :   U & V-points (same as ocean)
721      DO jj = 2, jpjm1
722         DO ji = fs_2, fs_jpim1   ! vect. opt.
723            zwndi_t = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( u_ice(ji-1,jj  ) + u_ice(ji,jj) )  )
724            zwndj_t = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( v_ice(ji  ,jj-1) + v_ice(ji,jj) )  )
725            wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1)
726         END DO
727      END DO
728      CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. )
729      !
730      ! Make ice-atm. drag dependent on ice concentration
731      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
732         CALL Cdn10_Lupkes2012( Cd_atm )
733         Ch_atm(:,:) = Cd_atm(:,:)       ! momentum and heat transfer coef. are considered identical
734      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
735         CALL Cdn10_Lupkes2015( Cd_atm, Ch_atm )
736      ENDIF
737
738      !!      CALL iom_put( "rCd_ice", Cd_atm)  ! output value of pure ice-atm. transfer coef.
739      !!      CALL iom_put( "Ch_ice", Ch_atm)  ! output value of pure ice-atm. transfer coef.
740
741      ! local scalars ( place there for vector optimisation purposes)
742
743      !!gm brutal....
744      utau_ice  (:,:) = 0._wp
745      vtau_ice  (:,:) = 0._wp
746      !!gm end
747
748      ! ------------------------------------------------------------ !
749      !    Wind stress relative to the moving ice ( U10m - U_ice )   !
750      ! ------------------------------------------------------------ !
751      ! C-grid ice dynamics :   U & V-points (same as ocean)
752      DO jj = 2, jpjm1
753         DO ji = fs_2, fs_jpim1   ! vect. opt.
754            utau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji+1,jj  ) + wndm_ice(ji,jj) )            &
755               &          * ( 0.5 * (sf(jp_wndi)%fnow(ji+1,jj,1) + sf(jp_wndi)%fnow(ji,jj,1) ) - rn_vfac * u_ice(ji,jj) )
756            vtau_ice(ji,jj) = 0.5 * rhoa(ji,jj) * Cd_atm(ji,jj) * ( wndm_ice(ji,jj+1  ) + wndm_ice(ji,jj) )            &
757               &          * ( 0.5 * (sf(jp_wndj)%fnow(ji,jj+1,1) + sf(jp_wndj)%fnow(ji,jj,1) ) - rn_vfac * v_ice(ji,jj) )
758         END DO
759      END DO
760      CALL lbc_lnk_multi( 'sbcblk', utau_ice, 'U', -1., vtau_ice, 'V', -1. )
761      !
762      !
763      IF(ln_ctl) THEN
764         CALL prt_ctl(tab2d_1=utau_ice  , clinfo1=' blk_ice: utau_ice : ', tab2d_2=vtau_ice  , clinfo2=' vtau_ice : ')
765         CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
766      ENDIF
767      !
768   END SUBROUTINE blk_ice_tau
769
770
771   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb )
772      !!---------------------------------------------------------------------
773      !!                     ***  ROUTINE blk_ice_flx  ***
774      !!
775      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
776      !!
777      !! ** Method  :   compute heat and freshwater exchanged
778      !!                between atmosphere and sea-ice using bulk formulation
779      !!                formulea, ice variables and read atmmospheric fields.
780      !!
781      !! caution : the net upward water flux has with mm/day unit
782      !!---------------------------------------------------------------------
783      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature
784      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
785      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
786      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
787      !!
788      INTEGER  ::   ji, jj, jl               ! dummy loop indices
789      REAL(wp) ::   zst3                     ! local variable
790      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
791      REAL(wp) ::   zztmp, z1_rLsub          !   -      -
792      REAL(wp) ::   zfr1, zfr2               ! local variables
793      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
794      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice
795      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qsb         ! sensible  heat flux over ice
796      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqlw        ! long wave heat sensitivity over ice
797      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_dqsb        ! sensible  heat sensitivity over ice
798      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
799      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
800      !!---------------------------------------------------------------------
801      !
802      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
803      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
804      !
805
806      !LB:
807      SELECT CASE( nhumi )
808      CASE( np_humi_sph )
809         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity!
810      CASE( np_humi_dpt )
811         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm
812         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
813      CASE( np_humi_rlh )
814         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm
815         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
816      END SELECT
817      !LB.
818
819      zztmp = 1. / ( 1. - albo )
820      WHERE( ptsu(:,:,:) /= 0._wp )
821         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
822      ELSEWHERE
823         z1_st(:,:,:) = 0._wp
824      END WHERE
825      !                                     ! ========================== !
826      DO jl = 1, jpl                        !  Loop over ice categories  !
827         !                                  ! ========================== !
828         DO jj = 1 , jpj
829            DO ji = 1, jpi
830               ! ----------------------------!
831               !      I   Radiative FLUXES   !
832               ! ----------------------------!
833               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
834               ! Short Wave (sw)
835               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
836               ! Long  Wave (lw)
837               z_qlw(ji,jj,jl) = 0.95 * ( sf(jp_qlw)%fnow(ji,jj,1) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
838               ! lw sensitivity
839               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
840
841               ! ----------------------------!
842               !     II    Turbulent FLUXES  !
843               ! ----------------------------!
844
845               ! ... turbulent heat fluxes with Ch_atm recalculated in blk_ice_tau
846               ! Sensible Heat
847               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))
848               ! Latent Heat
849               qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, rhoa(ji,jj) * rLsub  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
850                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / rhoa(ji,jj) - zqair(ji,jj) ) )
851               ! Latent heat sensitivity for ice (Dqla/Dt)
852               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
853                  dqla_ice(ji,jj,jl) = rn_efac * zcoef_dqla * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
854                     &                 z1_st(ji,jj,jl)*z1_st(ji,jj,jl) * EXP(-5897.8 * z1_st(ji,jj,jl))
855               ELSE
856                  dqla_ice(ji,jj,jl) = 0._wp
857               ENDIF
858
859               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
860               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_atm(ji,jj) * wndm_ice(ji,jj)
861
862               ! ----------------------------!
863               !     III    Total FLUXES     !
864               ! ----------------------------!
865               ! Downward Non Solar flux
866               qns_ice (ji,jj,jl) =     z_qlw (ji,jj,jl) - z_qsb (ji,jj,jl) - qla_ice (ji,jj,jl)
867               ! Total non solar heat flux sensitivity for ice
868               dqns_ice(ji,jj,jl) = - ( z_dqlw(ji,jj,jl) + z_dqsb(ji,jj,jl) + dqla_ice(ji,jj,jl) )
869            END DO
870            !
871         END DO
872         !
873      END DO
874      !
875      tprecip(:,:) = sf(jp_prec)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! total precipitation [kg/m2/s]
876      sprecip(:,:) = sf(jp_snow)%fnow(:,:,1) * rn_pfac * tmask(:,:,1)  ! solid precipitation [kg/m2/s]
877      CALL iom_put( 'snowpre', sprecip )                    ! Snow precipitation
878      CALL iom_put( 'precip' , tprecip )                    ! Total precipitation
879
880      ! --- evaporation --- !
881      z1_rLsub = 1._wp / rLsub
882      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation
883      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT
884      zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean
885
886      ! --- evaporation minus precipitation --- !
887      zsnw(:,:) = 0._wp
888      CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
889      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
890      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
891      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
892
893      ! --- heat flux associated with emp --- !
894      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
895         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
896         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
897         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
898      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
899         &              ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
900
901      ! --- total solar and non solar fluxes --- !
902      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
903         &           + qemp_ice(:,:) + qemp_oce(:,:)
904      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
905
906      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
907      qprec_ice(:,:) = rhos * ( ( MIN( sf(jp_tair)%fnow(:,:,1), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
908
909      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
910      DO jl = 1, jpl
911         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
912         !                         ! But we do not have Tice => consider it at 0degC => evap=0
913      END DO
914
915      ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- !
916      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm
917      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1
918      !
919      WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm
920         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) )
921      ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm
922         qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1
923      ELSEWHERE                                                         ! zero when hs>0
924         qtr_ice_top(:,:,:) = 0._wp
925      END WHERE
926      !
927      IF(ln_ctl) THEN
928         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl)
929         CALL prt_ctl(tab3d_1=z_qlw   , clinfo1=' blk_ice: z_qlw    : ', tab3d_2=dqla_ice, clinfo2=' dqla_ice : ', kdim=jpl)
930         CALL prt_ctl(tab3d_1=z_dqsb  , clinfo1=' blk_ice: z_dqsb   : ', tab3d_2=z_dqlw  , clinfo2=' z_dqlw   : ', kdim=jpl)
931         CALL prt_ctl(tab3d_1=dqns_ice, clinfo1=' blk_ice: dqns_ice : ', tab3d_2=qsr_ice , clinfo2=' qsr_ice  : ', kdim=jpl)
932         CALL prt_ctl(tab3d_1=ptsu    , clinfo1=' blk_ice: ptsu     : ', tab3d_2=qns_ice , clinfo2=' qns_ice  : ', kdim=jpl)
933         CALL prt_ctl(tab2d_1=tprecip , clinfo1=' blk_ice: tprecip  : ', tab2d_2=sprecip , clinfo2=' sprecip  : ')
934      ENDIF
935      !
936   END SUBROUTINE blk_ice_flx
937
938
939   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
940      !!---------------------------------------------------------------------
941      !!                     ***  ROUTINE blk_ice_qcn  ***
942      !!
943      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
944      !!                to force sea ice / snow thermodynamics
945      !!                in the case conduction flux is emulated
946      !!
947      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
948      !!                following the 0-layer Semtner (1976) approach
949      !!
950      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
951      !!              - qcn_ice : surface inner conduction flux (W/m2)
952      !!
953      !!---------------------------------------------------------------------
954      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
955      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptsu            ! sea ice / snow surface temperature
956      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   ptb             ! sea ice base temperature
957      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phs             ! snow thickness
958      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   phi             ! sea ice thickness
959      !
960      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
961      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
962      !
963      INTEGER  ::   ji, jj, jl           ! dummy loop indices
964      INTEGER  ::   iter                 ! local integer
965      REAL(wp) ::   zfac, zfac2, zfac3   ! local scalars
966      REAL(wp) ::   zkeff_h, ztsu, ztsu0 !
967      REAL(wp) ::   zqc, zqnet           !
968      REAL(wp) ::   zhe, zqa0            !
969      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zgfac   ! enhanced conduction factor
970      !!---------------------------------------------------------------------
971
972      ! -------------------------------------!
973      !      I   Enhanced conduction factor  !
974      ! -------------------------------------!
975      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
976      ! Fichefet and Morales Maqueda, JGR 1997
977      !
978      zgfac(:,:,:) = 1._wp
979
980      IF( ld_virtual_itd ) THEN
981         !
982         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
983         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
984         zfac3 = 2._wp / zepsilon
985         !
986         DO jl = 1, jpl
987            DO jj = 1 , jpj
988               DO ji = 1, jpi
989                  zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness
990                  IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor
991               END DO
992            END DO
993         END DO
994         !
995      ENDIF
996
997      ! -------------------------------------------------------------!
998      !      II   Surface temperature and conduction flux            !
999      ! -------------------------------------------------------------!
1000      !
1001      zfac = rcnd_i * rn_cnd_s
1002      !
1003      DO jl = 1, jpl
1004         DO jj = 1 , jpj
1005            DO ji = 1, jpi
1006               !
1007               zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness
1008                  &      ( rcnd_i * phs(ji,jj,jl) + rn_cnd_s * MAX( 0.01, phi(ji,jj,jl) ) )
1009               ztsu    = ptsu(ji,jj,jl)                                                ! Store current iteration temperature
1010               ztsu0   = ptsu(ji,jj,jl)                                                ! Store initial surface temperature
1011               zqa0    = qsr_ice(ji,jj,jl) - qtr_ice_top(ji,jj,jl) + qns_ice(ji,jj,jl) ! Net initial atmospheric heat flux
1012               !
1013               DO iter = 1, nit     ! --- Iterative loop
1014                  zqc   = zkeff_h * ( ztsu - ptb(ji,jj) )                              ! Conduction heat flux through snow-ice system (>0 downwards)
1015                  zqnet = zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc  ! Surface energy budget
1016                  ztsu  = ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h )              ! Temperature update
1017               END DO
1018               !
1019               ptsu   (ji,jj,jl) = MIN( rt0, ztsu )
1020               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) )
1021               qns_ice(ji,jj,jl) = qns_ice(ji,jj,jl) + dqns_ice(ji,jj,jl) * ( ptsu(ji,jj,jl) - ztsu0 )
1022               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) )  &
1023                  &   * MAX( 0._wp , SIGN( 1._wp, ptsu(ji,jj,jl) - rt0 ) )
1024
1025               ! --- Diagnose the heat loss due to changing non-solar flux (as in icethd_zdf_bl99) --- !
1026               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)
1027
1028            END DO
1029         END DO
1030         !
1031      END DO
1032      !
1033   END SUBROUTINE blk_ice_qcn
1034
1035
1036   SUBROUTINE Cdn10_Lupkes2012( Cd )
1037      !!----------------------------------------------------------------------
1038      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1039      !!
1040      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1041      !!                 to make it dependent on edges at leads, melt ponds and flows.
1042      !!                 After some approximations, this can be resumed to a dependency
1043      !!                 on ice concentration.
1044      !!
1045      !! ** Method :     The parameterization is taken from Lupkes et al. (2012) eq.(50)
1046      !!                 with the highest level of approximation: level4, eq.(59)
1047      !!                 The generic drag over a cell partly covered by ice can be re-written as follows:
1048      !!
1049      !!                 Cd = Cdw * (1-A) + Cdi * A + Ce * (1-A)**(nu+1/(10*beta)) * A**mu
1050      !!
1051      !!                    Ce = 2.23e-3       , as suggested by Lupkes (eq. 59)
1052      !!                    nu = mu = beta = 1 , as suggested by Lupkes (eq. 59)
1053      !!                    A is the concentration of ice minus melt ponds (if any)
1054      !!
1055      !!                 This new drag has a parabolic shape (as a function of A) starting at
1056      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1057      !!                 and going down to Cdi(say 1.4e-3) for A=1
1058      !!
1059      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
1060      !!                 => see Lupkes et al (2013)
1061      !!
1062      !! ** References : Lupkes et al. JGR 2012 (theory)
1063      !!                 Lupkes et al. GRL 2013 (application to GCM)
1064      !!
1065      !!----------------------------------------------------------------------
1066      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1067      REAL(wp), PARAMETER ::   zCe   = 2.23e-03_wp
1068      REAL(wp), PARAMETER ::   znu   = 1._wp
1069      REAL(wp), PARAMETER ::   zmu   = 1._wp
1070      REAL(wp), PARAMETER ::   zbeta = 1._wp
1071      REAL(wp)            ::   zcoef
1072      !!----------------------------------------------------------------------
1073      zcoef = znu + 1._wp / ( 10._wp * zbeta )
1074
1075      ! generic drag over a cell partly covered by ice
1076      !!Cd(:,:) = Cd_oce(:,:) * ( 1._wp - at_i_b(:,:) ) +  &                        ! pure ocean drag
1077      !!   &      rCd_ice      *           at_i_b(:,:)   +  &                        ! pure ice drag
1078      !!   &      zCe         * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**zmu   ! change due to sea-ice morphology
1079
1080      ! ice-atm drag
1081      Cd(:,:) = rCd_ice +  &                                                         ! pure ice drag
1082         &      zCe    * ( 1._wp - at_i_b(:,:) )**zcoef * at_i_b(:,:)**(zmu-1._wp)  ! change due to sea-ice morphology
1083
1084   END SUBROUTINE Cdn10_Lupkes2012
1085
1086
1087   SUBROUTINE Cdn10_Lupkes2015( Cd, Ch )
1088      !!----------------------------------------------------------------------
1089      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1090      !!
1091      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1092      !!                 between sea-ice and atmosphere with distinct momentum
1093      !!                 and heat coefficients depending on sea-ice concentration
1094      !!                 and atmospheric stability (no meltponds effect for now).
1095      !!
1096      !! ** Method :     The parameterization is adapted from Lupkes et al. (2015)
1097      !!                 and ECHAM6 atmospheric model. Compared to Lupkes2012 scheme,
1098      !!                 it considers specific skin and form drags (Andreas et al. 2010)
1099      !!                 to compute neutral transfert coefficients for both heat and
1100      !!                 momemtum fluxes. Atmospheric stability effect on transfert
1101      !!                 coefficient is also taken into account following Louis (1979).
1102      !!
1103      !! ** References : Lupkes et al. JGR 2015 (theory)
1104      !!                 Lupkes et al. ECHAM6 documentation 2015 (implementation)
1105      !!
1106      !!----------------------------------------------------------------------
1107      !
1108      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Cd
1109      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   Ch
1110      REAL(wp), DIMENSION(jpi,jpj)            ::   ztm_su, zst, zqo_sat, zqi_sat
1111      !
1112      ! ECHAM6 constants
1113      REAL(wp), PARAMETER ::   z0_skin_ice  = 0.69e-3_wp  ! Eq. 43 [m]
1114      REAL(wp), PARAMETER ::   z0_form_ice  = 0.57e-3_wp  ! Eq. 42 [m]
1115      REAL(wp), PARAMETER ::   z0_ice       = 1.00e-3_wp  ! Eq. 15 [m]
1116      REAL(wp), PARAMETER ::   zce10        = 2.80e-3_wp  ! Eq. 41
1117      REAL(wp), PARAMETER ::   zbeta        = 1.1_wp      ! Eq. 41
1118      REAL(wp), PARAMETER ::   zc           = 5._wp       ! Eq. 13
1119      REAL(wp), PARAMETER ::   zc2          = zc * zc
1120      REAL(wp), PARAMETER ::   zam          = 2. * zc     ! Eq. 14
1121      REAL(wp), PARAMETER ::   zah          = 3. * zc     ! Eq. 30
1122      REAL(wp), PARAMETER ::   z1_alpha     = 1._wp / 0.2_wp  ! Eq. 51
1123      REAL(wp), PARAMETER ::   z1_alphaf    = z1_alpha    ! Eq. 56
1124      REAL(wp), PARAMETER ::   zbetah       = 1.e-3_wp    ! Eq. 26
1125      REAL(wp), PARAMETER ::   zgamma       = 1.25_wp     ! Eq. 26
1126      REAL(wp), PARAMETER ::   z1_gamma     = 1._wp / zgamma
1127      REAL(wp), PARAMETER ::   r1_3         = 1._wp / 3._wp
1128      !
1129      INTEGER  ::   ji, jj         ! dummy loop indices
1130      REAL(wp) ::   zthetav_os, zthetav_is, zthetav_zu
1131      REAL(wp) ::   zrib_o, zrib_i
1132      REAL(wp) ::   zCdn_skin_ice, zCdn_form_ice, zCdn_ice
1133      REAL(wp) ::   zChn_skin_ice, zChn_form_ice
1134      REAL(wp) ::   z0w, z0i, zfmi, zfmw, zfhi, zfhw
1135      REAL(wp) ::   zCdn_form_tmp
1136      !!----------------------------------------------------------------------
1137
1138      ! mean temperature
1139      WHERE( at_i_b(:,:) > 1.e-20 )
1140         ztm_su(:,:) = SUM( t_su(:,:,:) * a_i_b(:,:,:) , dim=3 ) / at_i_b(:,:)
1141      ELSEWHERE
1142         ztm_su(:,:) = rt0
1143      ENDWHERE
1144
1145      ! Momentum Neutral Transfert Coefficients (should be a constant)
1146      zCdn_form_tmp = zce10 * ( LOG( 10._wp / z0_form_ice + 1._wp ) / LOG( rn_zu / z0_form_ice + 1._wp ) )**2   ! Eq. 40
1147      zCdn_skin_ice = ( vkarmn                                      / LOG( rn_zu / z0_skin_ice + 1._wp ) )**2   ! Eq. 7
1148      zCdn_ice      = zCdn_skin_ice   ! Eq. 7 (cf Lupkes email for details)
1149      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1150
1151      ! Heat Neutral Transfert Coefficients
1152      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)
1153
1154      ! Atmospheric and Surface Variables
1155      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
1156      zqo_sat(:,:) = rdct_qsat_salt * q_sat( zst(:,:)   , sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ocean [kg/kg]
1157      zqi_sat(:,:) =                  q_sat( ztm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] !LB: no 0.98 !!(rdct_qsat_salt)
1158      !
1159      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility
1160         DO ji = fs_2, fs_jpim1
1161            ! Virtual potential temperature [K]
1162            zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean
1163            zthetav_is = ztm_su(ji,jj) * ( 1._wp + rctv0 * zqi_sat(ji,jj) )   ! ocean ice
1164            zthetav_zu = t_zu (ji,jj)  * ( 1._wp + rctv0 * q_zu(ji,jj)    )   ! at zu
1165
1166            ! Bulk Richardson Number (could use Ri_bulk function from aerobulk instead)
1167            zrib_o = grav / zthetav_os * ( zthetav_zu - zthetav_os ) * rn_zu / MAX( 0.5, wndm(ji,jj)     )**2   ! over ocean
1168            zrib_i = grav / zthetav_is * ( zthetav_zu - zthetav_is ) * rn_zu / MAX( 0.5, wndm_ice(ji,jj) )**2   ! over ice
1169
1170            ! Momentum and Heat Neutral Transfert Coefficients
1171            zCdn_form_ice = zCdn_form_tmp * at_i_b(ji,jj) * ( 1._wp - at_i_b(ji,jj) )**zbeta  ! Eq. 40
1172            zChn_form_ice = zCdn_form_ice / ( 1._wp + ( LOG( z1_alphaf ) / vkarmn ) * SQRT( zCdn_form_ice ) )               ! Eq. 53
1173
1174            ! Momentum and Heat Stability functions (possibility to use psi_m_ecmwf instead)
1175            z0w = rn_zu * EXP( -1._wp * vkarmn / SQRT( Cdn_oce(ji,jj) ) ) ! over water
1176            z0i = z0_skin_ice                                             ! over ice (cf Lupkes email for details)
1177            IF( zrib_o <= 0._wp ) THEN
1178               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
1179               zfhw = ( 1._wp + ( zbetah * ( zthetav_os - zthetav_zu )**r1_3 / ( Chn_oce(ji,jj) * MAX(0.01, wndm(ji,jj)) )   &     ! Eq. 26
1180                  &             )**zgamma )**z1_gamma
1181            ELSE
1182               zfmw = 1._wp / ( 1._wp + zam * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 12
1183               zfhw = 1._wp / ( 1._wp + zah * zrib_o / SQRT( 1._wp + zrib_o ) )   ! Eq. 28
1184            ENDIF
1185
1186            IF( zrib_i <= 0._wp ) THEN
1187               zfmi = 1._wp - zam * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq.  9
1188               zfhi = 1._wp - zah * zrib_i / (1._wp + 3._wp * zc2 * zCdn_ice * SQRT( -zrib_i * ( rn_zu / z0i + 1._wp)))   ! Eq. 25
1189            ELSE
1190               zfmi = 1._wp / ( 1._wp + zam * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 11
1191               zfhi = 1._wp / ( 1._wp + zah * zrib_i / SQRT( 1._wp + zrib_i ) )   ! Eq. 27
1192            ENDIF
1193
1194            ! Momentum Transfert Coefficients (Eq. 38)
1195            Cd(ji,jj) = zCdn_skin_ice *   zfmi +  &
1196               &        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) )
1197
1198            ! Heat Transfert Coefficients (Eq. 49)
1199            Ch(ji,jj) = zChn_skin_ice *   zfhi +  &
1200               &        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) )
1201            !
1202         END DO
1203      END DO
1204      CALL lbc_lnk_multi( 'sbcblk', Cd, 'T',  1., Ch, 'T', 1. )
1205      !
1206   END SUBROUTINE Cdn10_Lupkes2015
1207
1208#endif
1209
1210   !!======================================================================
1211END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.