source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90 @ 11348

Last change on this file since 11348 was 11348, checked in by gsamson, 2 years ago

dev_r11265_ABL :

  • merge HPC-13_IRRMANN_BDY_optimization branch @ rev11332 with dev_r11265_ABL branch @ rev11334
  • allow ln_dm2dc option with ABL
  • cosmetic change in sbcabl.F90

identical results with rev11334 for bulk and abl orca2

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