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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk.F90 in branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

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

Last change on this file since 8321 was 8321, checked in by clem, 7 years ago

STEP3 (1): clean separation between sea-ice and ocean

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