source: NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/OCE/SBC/sbcblk.F90 @ 12811

Last change on this file since 12811 was 12811, checked in by clem, 13 months ago

debug a restartability issue. All sette tests passed now

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