New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sbcblk.F90 in NEMO/trunk/src/OCE/SBC – NEMO

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

Last change on this file since 13472 was 13472, checked in by smasson, 4 years ago

trunk: commit changes from r4.0-HEAD from 13284 to 13449, see #2523

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