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/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/sbcblk.F90 @ 9727

Last change on this file since 9727 was 9727, checked in by mathiot, 7 years ago

Fix for #2083 (surface flux not masked: bug if isf cavity) in trunk

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