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

Last change on this file since 11266 was 11266, checked in by laurent, 15 months ago

LB: CSWL param now uses NEMO time step (rdt) and previous-t value of t_skin as first guess.

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