source: NEMO/branches/UKMO/NEMO_4.0.2_FKOSM/src/OCE/SBC/sbcblk.F90 @ 12912

Last change on this file since 12912 was 12912, checked in by cguiavarch, 9 months ago

OSMOSIS and Fox-Kemper changes (merged from NEMO_4.0.1_FKOSM)

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