source: NEMO/branches/2019/dev_r11265_ASINTER-01_Guillaume_ABL1D/src/OCE/SBC/sbcblk.F90 @ 11586

Last change on this file since 11586 was 11586, checked in by gsamson, 7 months ago

dev_r11265_ABL : see #2131

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