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

Last change on this file since 11304 was 11304, checked in by smasson, 22 months ago

dev_r11265_ABL : si3 compatibility, see #2131

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