source: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 @ 8585

Last change on this file since 8585 was 8585, checked in by clem, 3 years ago

add option Lupkes2015 for ice-atm drag

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