New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk.F90 in NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC – NEMO

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

Last change on this file since 11360 was 11360, checked in by gsamson, 5 years ago

dev_r11265_ABL : see #2131

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