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 NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcblk.F90 @ 11960

Last change on this file since 11960 was 11960, checked in by acc, 5 years ago

Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. (svn merge -r 11614:11954). Resolved tree conflicts and one actual conflict. Sette tested(these changes alter the ext/AGRIF reference; remember to update). See ticket #2341

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