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_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC – NEMO

source: branches/2017/dev_r7881_ENHANCE09_RK3/NEMOGCM/NEMO/RK3_SRC/SBC/sbcblk.F90 @ 8637

Last change on this file since 8637 was 8637, checked in by gm, 5 years ago

#1911 (ENHANCE-09): PART I.3 - phasing with updated branch dev_r8183_ICEMODEL revision 8626

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