source: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/SPITZ12/MY_SRC/sbcblk.F90 @ 9552

Last change on this file since 9552 was 9552, checked in by clem, 3 years ago

add a new configuration SPITZ12 which is a regional simulation around the Spitzbergen (Svalbard archipelago) which includes tides and sea-ice with open boundaries.

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