source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90 @ 9098

Last change on this file since 9098 was 9098, checked in by flavoni, 4 years ago

change lbc_lnk in lbc_lnk_multi

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