source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk.F90 @ 12015

Last change on this file since 12015 was 12015, checked in by gsamson, 3 months ago

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

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