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

Last change on this file since 13472 was 13472, checked in by smasson, 7 months 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
Line 
1MODULE sbcblk
2   !!======================================================================
3   !!                       ***  MODULE  sbcblk  ***
4   !! Ocean forcing:  momentum, heat and freshwater flux formulation
5   !!                         Aerodynamic Bulk Formulas
6   !!                        SUCCESSOR OF "sbcblk_core"
7   !!=====================================================================
8   !! History :  1.0  !  2004-08  (U. Schweckendiek)  Original CORE code
9   !!            2.0  !  2005-04  (L. Brodeau, A.M. Treguier)  improved CORE bulk and its user interface
10   !!            3.0  !  2006-06  (G. Madec)  sbc rewritting
11   !!             -   !  2006-12  (L. Brodeau)  Original code for turb_core
12   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
13   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
14   !!            3.4  !  2011-11  (C. Harris)  Fill arrays required by CICE
15   !!            3.7  !  2014-06  (L. Brodeau)  simplification and optimization of CORE bulk
16   !!            4.0  !  2016-06  (L. Brodeau)  sbcblk_core becomes sbcblk and is not restricted to the CORE algorithm anymore
17   !!                 !                        ==> based on AeroBulk (https://github.com/brodeau/aerobulk/)
18   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine
19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce conduction flux emulator (M. Vancoppenolle)
20   !!            4.0  !  2019-03  (F. Lemarié & G. Samson)  add ABL compatibility (ln_abl=TRUE)
21   !!----------------------------------------------------------------------
22
23   !!----------------------------------------------------------------------
24   !!   sbc_blk_init  : initialisation of the chosen bulk formulation as ocean surface boundary condition
25   !!   sbc_blk       : bulk formulation as ocean surface boundary condition
26   !!   blk_oce_1     : computes pieces of momentum, heat and freshwater fluxes over ocean for ABL model  (ln_abl=TRUE)
27   !!   blk_oce_2     : finalizes momentum, heat and freshwater fluxes computation over ocean after the ABL step  (ln_abl=TRUE)
28   !!             sea-ice case only :
29   !!   blk_ice_1   : provide the air-ice stress
30   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
31   !!   blk_ice_qcn   : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
32   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
33   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
34   !!----------------------------------------------------------------------
35   USE oce            ! ocean dynamics and tracers
36   USE dom_oce        ! ocean space and time domain
37   USE phycst         ! physical constants
38   USE fldread        ! read input fields
39   USE sbc_oce        ! Surface boundary condition: ocean fields
40   USE cyclone        ! Cyclone 10m wind form trac of cyclone centres
41   USE sbcdcy         ! surface boundary condition: diurnal cycle
42   USE sbcwave , ONLY :   cdn_wave ! wave module
43   USE sbc_ice        ! Surface boundary condition: ice fields
44   USE lib_fortran    ! to use key_nosignedzero
45#if defined key_si3
46   USE ice     , ONLY :   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
48#endif
49   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)
50   USE sbcblk_algo_coare3p0 ! => turb_coare3p0 : COAREv3.0 (Fairall et al. 2003)
51   USE sbcblk_algo_coare3p6 ! => turb_coare3p6 : COAREv3.6 (Fairall et al. 2018 + Edson et al. 2013)
52   USE sbcblk_algo_ecmwf    ! => turb_ecmwf    : ECMWF (IFS cycle 45r1)
53   !
54   USE iom            ! I/O manager library
55   USE in_out_manager ! I/O manager
56   USE lib_mpp        ! distribued memory computing library
57   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
58   USE prtctl         ! Print control
59
60   USE sbcblk_phy     ! a catalog of functions for physical/meteorological parameters in the marine boundary layer, rho_air, q_sat, etc...
61
62
63   IMPLICIT NONE
64   PRIVATE
65
66   PUBLIC   sbc_blk_init  ! called in sbcmod
67   PUBLIC   sbc_blk       ! called in sbcmod
68   PUBLIC   blk_oce_1     ! called in sbcabl
69   PUBLIC   blk_oce_2     ! called in sbcabl
70#if defined key_si3
71   PUBLIC   blk_ice_1     ! routine called in icesbc
72   PUBLIC   blk_ice_2     ! routine called in icesbc
73   PUBLIC   blk_ice_qcn   ! routine called in icesbc
74#endif
75
76   INTEGER , PUBLIC, 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
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
93
94   ! Warning: keep this structure allocatable for Agrif...
95   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read)
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)
100   LOGICAL  ::   ln_COARE_3p6   ! "COARE 3.6" algorithm   (Edson et al. 2013)
101   LOGICAL  ::   ln_ECMWF       ! "ECMWF"     algorithm   (IFS cycle 45r1)
102   !
103   LOGICAL  ::   ln_Cd_L12      ! ice-atm drag = F( ice concentration )                        (Lupkes et al. JGR2012)
104   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015)
105   !
106   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   !
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   !
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)
118
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
124   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature
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
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)
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)
138
139   !! * Substitutions
140#  include "do_loop_substitute.h90"
141   !!----------------------------------------------------------------------
142   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
143   !! $Id$
144   !! Software governed by the CeCILL license (see ./LICENSE)
145   !!----------------------------------------------------------------------
146CONTAINS
147
148   INTEGER FUNCTION sbc_blk_alloc()
149      !!-------------------------------------------------------------------
150      !!             ***  ROUTINE sbc_blk_alloc ***
151      !!-------------------------------------------------------------------
152      ALLOCATE( t_zu(jpi,jpj)   , q_zu(jpi,jpj)   ,                                      &
153         &      Cdn_oce(jpi,jpj), Chn_oce(jpi,jpj), Cen_oce(jpi,jpj),                    &
154         &      Cd_ice (jpi,jpj), Ch_ice (jpi,jpj), Ce_ice (jpi,jpj), STAT=sbc_blk_alloc )
155      !
156      CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
157      IF( sbc_blk_alloc /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_alloc: failed to allocate arrays' )
158   END FUNCTION sbc_blk_alloc
159
160
161   SUBROUTINE sbc_blk_init
162      !!---------------------------------------------------------------------
163      !!                    ***  ROUTINE sbc_blk_init  ***
164      !!
165      !! ** Purpose :   choose and initialize a bulk formulae formulation
166      !!
167      !! ** Method  :
168      !!
169      !!----------------------------------------------------------------------
170      INTEGER  ::   jfpr                  ! dummy loop indice and argument
171      INTEGER  ::   ios, ierror, ioptio   ! Local integer
172      !!
173      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files
174      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read
175      TYPE(FLD_N) ::   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             !       "                        "
178      TYPE(FLD_N) ::   sn_cc, sn_hpgi, sn_hpgj                 !       "                        "
179      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable
180      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
181         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     &
182         &                 sn_cc, sn_hpgi, sn_hpgj,                                   &
183         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm
184         &                 cn_dir , rn_zqt, rn_zu,                                    &
185         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot,           &
186         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback
187         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB
188      !!---------------------------------------------------------------------
189      !
190      !                                      ! allocate sbc_blk_core array
191      IF( sbc_blk_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_blk : unable to allocate standard arrays' )
192      !
193      !                             !** read bulk namelist
194      READ  ( numnam_ref, namsbc_blk, IOSTAT = ios, ERR = 901)
195901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namsbc_blk in reference namelist' )
196      !
197      READ  ( numnam_cfg, namsbc_blk, IOSTAT = ios, ERR = 902 )
198902   IF( ios >  0 )   CALL ctl_nam ( ios , 'namsbc_blk in configuration namelist' )
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
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
217      IF( ioptio /= 1 )   CALL ctl_stop( 'sbc_blk_init: Choose one and only one bulk algorithm' )
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' )
247      !
248      IF( ln_dm2dc ) THEN                 !* check: diurnal cycle on Qsr
249         IF( sn_qsr%freqh /= 24. )   CALL ctl_stop( 'sbc_blk_init: ln_dm2dc=T only with daily short-wave input' )
250         IF( sn_qsr%ln_tint ) THEN
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
258      !
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
263      slf_i(jp_slp  ) = sn_slp     ;   slf_i(jp_cc   ) = sn_cc
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
266      !
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      !
272      !                                      !- allocate the bulk structure
273      ALLOCATE( sf(jpfld), STAT=ierror )
274      IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_init: unable to allocate sf structure' )
275      !
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      !
279      DO jfpr= 1, jpfld
280         !
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)
292            IF(     jfpr == jp_slp ) THEN
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
298            ELSEIF( jfpr == jp_cc  ) THEN
299               sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf
300            ELSE
301               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr
302               CALL ctl_stop( ctmp1 )
303            ENDIF
304         ELSE                                                  !-- used field  --!
305            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation
306            !
307            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   &
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...' )
310         ENDIF
311      END DO
312      !
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
316            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )
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')
322         ENDIF
323      ELSE
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
331      !
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      !
344      IF(lwp) THEN                     !** Control print
345         !
346         WRITE(numout,*)                  !* namelist
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
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
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))'
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
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
365         !
366         WRITE(numout,*)
367         SELECT CASE( nblk )              !* Print the choice of bulk algorithm
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)'
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)'
372         END SELECT
373         !
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         !
385      ENDIF
386      !
387   END SUBROUTINE sbc_blk_init
388
389
390   SUBROUTINE sbc_blk( kt )
391      !!---------------------------------------------------------------------
392      !!                    ***  ROUTINE sbc_blk  ***
393      !!
394      !! ** Purpose :   provide at each time step the surface ocean fluxes
395      !!              (momentum, heat, freshwater and runoff)
396      !!
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
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
425      !!----------------------------------------------------------------------
426      REAL(wp), DIMENSION(jpi,jpj) ::   zssq, zcd_du, zsen, zevp
427      REAL(wp) :: ztmp
428      !!----------------------------------------------------------------------
429      !
430      CALL fld_read( kt, nn_fsbc, sf )             ! input fields provided at the current time-step
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 )
461      !                                            ! compute the surface ocean fluxes using bulk formulea
462      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
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
469
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
474      ENDIF
475      !
476#if defined key_cice
477      IF( MOD( kt - 1, nn_fsbc ) == 0 )   THEN
478         qlw_ice(:,:,1)   = sf(jp_qlw )%fnow(:,:,1)
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
484         tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
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
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)
499      ENDIF
500#endif
501      !
502   END SUBROUTINE sbc_blk
503
504
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
509      !!---------------------------------------------------------------------
510      !!                     ***  ROUTINE blk_oce_1  ***
511      !!
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
514      !!
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
518      !!
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)
523      !!---------------------------------------------------------------------
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]
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]
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]
542      !
543      INTEGER  ::   ji, jj               ! dummy loop indices
544      REAL(wp) ::   zztmp                ! local variable
545      REAL(wp) ::   zstmax, zstau
546#if defined key_cyclone
547      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point
548#endif
549      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point
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]
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
558      !!---------------------------------------------------------------------
559      !
560      ! local scalars ( place there for vector optimisation purposes)
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!)
563
564      ! --- cloud cover --- !
565      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1)
566
567      ! ----------------------------------------------------------------------------- !
568      !      0   Wind components and module at T-point relative to the moving ocean   !
569      ! ----------------------------------------------------------------------------- !
570
571      ! ... components ( U10m - U_oce ) at T-point (unmasked)
572#if defined key_cyclone
573      zwnd_i(:,:) = 0._wp
574      zwnd_j(:,:) = 0._wp
575      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012)
576      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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) )
581      END_2D
582#else
583      ! ... scalar wind module at T-point (not masked)
584      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
585         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  )
586      END_2D
587#endif
588      ! ----------------------------------------------------------------------------- !
589      !      I   Solar FLUX                                                           !
590      ! ----------------------------------------------------------------------------- !
591
592      ! ocean albedo assumed to be constant + modify now Qsr to include the diurnal cycle                    ! Short Wave
593      zztmp = 1. - albo
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)
598      ENDIF
599
600
601      ! ----------------------------------------------------------------------------- !
602      !     II   Turbulent FLUXES                                                     !
603      ! ----------------------------------------------------------------------------- !
604
605      ! specific humidity at SST
606      pssq(:,:) = rdct_qsat_salt * q_sat( ptsk(:,:), pslp(:,:) )
607
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
625      !
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() ...
635         IF( ln_tpot ) THEN
636            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
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
667      CASE DEFAULT
668         CALL ctl_stop( 'STOP', 'sbc_oce: non-existing bulk formula selected' )
669
670      END SELECT
671     
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     
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
692      END IF
693
694      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90
695      ! -------------------------------------------------------------
696
697      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp
698         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
699            zztmp = zU_zu(ji,jj)
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)
704            rhoa(ji,jj)   = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) )
705         END_2D
706      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation
707         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), &
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 )
712
713         zqla(:,:) = zqla(:,:) * tmask(:,:,1)
714         psen(:,:) = psen(:,:) * tmask(:,:,1)
715         taum(:,:) = taum(:,:) * tmask(:,:,1)
716         pevp(:,:) = pevp(:,:) * tmask(:,:,1)
717
718         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
733
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)
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
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
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
746         !     Note that coastal wind stress is not used in the code... so this extra care has no effect
747         DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T
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))
752         END_2D
753
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
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
776
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   : ' )
781      ENDIF
782      !
783   END SUBROUTINE blk_oce_1
784
785
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)
820
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
836      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
844      ENDIF
845
846      ! ----------------------------------------------------------------------------- !
847      !     IV    Total FLUXES                                                       !
848      ! ----------------------------------------------------------------------------- !
849      !
850      emp (:,:) = (  pevp(:,:)                                       &   ! mass flux (evap. - precip.)
851         &         - pprec(:,:) * rn_pfac  ) * tmask(:,:,1)
852      !
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
860      qns(:,:) = qns(:,:) * tmask(:,:,1)
861      !
862#if defined key_si3
863      qns_oce(:,:) = zqlw(:,:) + psen(:,:) + zqla(:,:)                             ! non solar without emp (only needed by SI3)
864      qsr_oce(:,:) = qsr(:,:)
865#endif
866      !
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      !
877      IF ( nn_ice == 0 ) THEN
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
882      ENDIF
883      !
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   : ')
888      ENDIF
889      !
890   END SUBROUTINE blk_oce_2
891
892
893#if defined key_si3
894   !!----------------------------------------------------------------------
895   !!   'key_si3'                                       SI3 sea-ice model
896   !!----------------------------------------------------------------------
897   !!   blk_ice_1   : provide the air-ice stress
898   !!   blk_ice_2   : provide the heat and mass fluxes at air-ice interface
899   !!   blk_ice_qcn : provide ice surface temperature and snow/ice conduction flux (emulating conduction flux)
900   !!   Cdn10_Lupkes2012 : Lupkes et al. (2012) air-ice drag
901   !!   Cdn10_Lupkes2015 : Lupkes et al. (2015) air-ice drag
902   !!----------------------------------------------------------------------
903
904   SUBROUTINE blk_ice_1( pwndi, pwndj, ptair, phumi, pslp , puice, pvice, ptsui,  &   ! inputs
905      &                  putaui, pvtaui, pseni, pevpi, pssqi, pcd_dui             )   ! optional outputs
906      !!---------------------------------------------------------------------
907      !!                     ***  ROUTINE blk_ice_1  ***
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      !!---------------------------------------------------------------------
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      !
930      INTEGER  ::   ji, jj    ! dummy loop indices
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)
934      !!---------------------------------------------------------------------
935      !
936
937      ! ------------------------------------------------------------ !
938      !    Wind module relative to the moving ice ( U10m - U_ice )   !
939      ! ------------------------------------------------------------ !
940      ! C-grid ice dynamics :   U & V-points (same as ocean)
941      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
942         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) )
943      END_2D
944      !
945      ! Make ice-atm. drag dependent on ice concentration
946      IF    ( ln_Cd_L12 ) THEN   ! calculate new drag from Lupkes(2012) equations
947         CALL Cdn10_Lupkes2012( Cd_ice )
948         Ch_ice(:,:) = Cd_ice(:,:)       ! momentum and heat transfer coef. are considered identical
949         Ce_ice(:,:) = Cd_ice(:,:)
950      ELSEIF( ln_Cd_L15 ) THEN   ! calculate new drag from Lupkes(2015) equations
951         CALL Cdn10_Lupkes2015( ptsui, pslp, Cd_ice, Ch_ice )
952         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical
953      ENDIF
954     
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)
958     
959      ! local scalars ( place there for vector optimisation purposes)
960      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:)
961
962      IF( ln_blk ) THEN
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...
967         DO_2D( 0, 1, 0, 1 )    ! at T point
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)
970         END_2D
971         !
972         DO_2D( 0, 0, 0, 0 )    ! U & V-points (same as ocean).
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
979         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp )
980         !
981         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   &
982            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' )
983      ELSE ! ln_abl
984         zztmp1 = 11637800.0_wp
985         zztmp2 =    -5897.8_wp
986         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
994      !
995      IF(sn_cfctl%l_prtctl)  CALL prt_ctl(tab2d_1=wndm_ice  , clinfo1=' blk_ice: wndm_ice : ')
996      !
997   END SUBROUTINE blk_ice_1
998
999
1000   SUBROUTINE blk_ice_2( ptsu, phs, phi, palb, ptair, phumi, pslp, pqlw, pprec, psnow  )
1001      !!---------------------------------------------------------------------
1002      !!                     ***  ROUTINE blk_ice_2  ***
1003      !!
1004      !! ** Purpose :   provide the heat and mass fluxes at air-ice interface
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      !!---------------------------------------------------------------------
1012      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature [K]
1013      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness
1014      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness
1015      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies)
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
1022      !!
1023      INTEGER  ::   ji, jj, jl               ! dummy loop indices
1024      REAL(wp) ::   zst3                     ! local variable
1025      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      -
1026      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      -
1027      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature
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
1032      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (SI3)
1033      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB
1034      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2
1035      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri
1036      !!---------------------------------------------------------------------
1037      !
1038      zcoef_dqlw = 4._wp * 0.95_wp * stefan             ! local scalars
1039      zcoef_dqla = -rLsub * 11637800._wp * (-5897.8_wp) !LB: BAD!
1040      !
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
1049      !
1050      zztmp = 1. / ( 1. - albo )
1051      WHERE( ptsu(:,:,:) /= 0._wp )
1052         z1_st(:,:,:) = 1._wp / ptsu(:,:,:)
1053      ELSEWHERE
1054         z1_st(:,:,:) = 0._wp
1055      END WHERE
1056      !                                     ! ========================== !
1057      DO jl = 1, jpl                        !  Loop over ice categories  !
1058         !                                  ! ========================== !
1059         DO jj = 1 , jpj
1060            DO ji = 1, jpi
1061               ! ----------------------------!
1062               !      I   Radiative FLUXES   !
1063               ! ----------------------------!
1064               zst3 = ptsu(ji,jj,jl) * ptsu(ji,jj,jl) * ptsu(ji,jj,jl)
1065               ! Short Wave (sw)
1066               qsr_ice(ji,jj,jl) = zztmp * ( 1. - palb(ji,jj,jl) ) * qsr(ji,jj)
1067               ! Long  Wave (lw)
1068               z_qlw(ji,jj,jl) = 0.95 * ( pqlw(ji,jj) - stefan * ptsu(ji,jj,jl) * zst3 ) * tmask(ji,jj,1)
1069               ! lw sensitivity
1070               z_dqlw(ji,jj,jl) = zcoef_dqlw * zst3
1071
1072               ! ----------------------------!
1073               !     II    Turbulent FLUXES  !
1074               ! ----------------------------!
1075
1076               ! ... turbulent heat fluxes with Ch_ice recalculated in blk_ice_1
1077               ! Sensible Heat
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))
1079               ! Latent Heat
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) ) )
1083               ! Latent heat sensitivity for ice (Dqla/Dt)
1084               IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
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
1087               ELSE
1088                  dqla_ice(ji,jj,jl) = 0._wp
1089               ENDIF
1090
1091               ! Sensible heat sensitivity (Dqsb_ice/Dtn_ice)
1092               z_dqsb(ji,jj,jl) = rhoa(ji,jj) * rCp_air * Ch_ice(ji,jj) * wndm_ice(ji,jj)
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      !
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
1111
1112      ! --- evaporation --- !
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
1116      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean  !LB: removed rn_efac here, correct???
1117
1118      ! --- evaporation minus precipitation --- !
1119      zsnw(:,:) = 0._wp
1120      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
1121      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
1122      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
1123      emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
1124
1125      ! --- heat flux associated with emp --- !
1126      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
1127         &          + ( tprecip(:,:) - sprecip(:,:) ) * ( ptair(:,:) - rt0 ) * rcp               & ! liquid precip at Tair
1128         &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
1129         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1130      qemp_ice(:,:) =   sprecip(:,:) * zsnw *                                                    & ! solid precip (only)
1131         &              ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1132
1133      ! --- total solar and non solar fluxes --- !
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 )
1137
1138      ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
1139      qprec_ice(:,:) = rhos * ( ( MIN( ptair(:,:), rt0 ) - rt0 ) * rcpi * tmask(:,:,1) - rLfus )
1140
1141      ! --- heat content of evap over ice in W/m2 (to be used in 1D-thermo) ---
1142      DO jl = 1, jpl
1143         qevap_ice(:,:,jl) = 0._wp ! should be -evap_ice(:,:,jl)*( ( Tice - rt0 ) * rcpi * tmask(:,:,1) )
1144         !                         ! But we do not have Tice => consider it at 0degC => evap=0
1145      END DO
1146
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
1169      !
1170      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN
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)
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(:,:) ) )
1177         IF( iom_use('hflx_rain_cea') )  CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * ztmp(:,:) )   ! heat flux from rain (cell average)
1178      ENDIF
1179      IF( iom_use('hflx_snow_cea') .OR. iom_use('hflx_snow_ao_cea') .OR. iom_use('hflx_snow_ai_cea')  )  THEN
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)
1189      ENDIF
1190      !
1191      IF(sn_cfctl%l_prtctl) THEN
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      !
1200   END SUBROUTINE blk_ice_2
1201
1202
1203   SUBROUTINE blk_ice_qcn( ld_virtual_itd, ptsu, ptb, phs, phi )
1204      !!---------------------------------------------------------------------
1205      !!                     ***  ROUTINE blk_ice_qcn  ***
1206      !!
1207      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux
1208      !!                to force sea ice / snow thermodynamics
1209      !!                in the case conduction flux is emulated
1210      !!
1211      !! ** Method  :   compute surface energy balance assuming neglecting heat storage
1212      !!                following the 0-layer Semtner (1976) approach
1213      !!
1214      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K)
1215      !!              - qcn_ice : surface inner conduction flux (W/m2)
1216      !!
1217      !!---------------------------------------------------------------------
1218      LOGICAL                   , INTENT(in   ) ::   ld_virtual_itd  ! single-category option
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
1223      !
1224      INTEGER , PARAMETER ::   nit = 10                  ! number of iterations
1225      REAL(wp), PARAMETER ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction
1226      !
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      !!---------------------------------------------------------------------
1235
1236      ! -------------------------------------!
1237      !      I   Enhanced conduction factor  !
1238      ! -------------------------------------!
1239      ! Emulates the enhancement of conduction by unresolved thin ice (ld_virtual_itd = T)
1240      ! Fichefet and Morales Maqueda, JGR 1997
1241      !
1242      zgfac(:,:,:) = 1._wp
1243
1244      IF( ld_virtual_itd ) THEN
1245         !
1246         zfac  = 1._wp /  ( rn_cnd_s + rcnd_i )
1247         zfac2 = EXP(1._wp) * 0.5_wp * zepsilon
1248         zfac3 = 2._wp / zepsilon
1249         !
1250         DO jl = 1, jpl
1251            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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
1255         END DO
1256         !
1257      ENDIF
1258
1259      ! -------------------------------------------------------------!
1260      !      II   Surface temperature and conduction flux            !
1261      ! -------------------------------------------------------------!
1262      !
1263      zfac = rcnd_i * rn_cnd_s
1264      !
1265      DO jl = 1, jpl
1266         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
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 ) )
1285
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)
1288
1289         END_2D
1290         !
1291      END DO
1292      !
1293   END SUBROUTINE blk_ice_qcn
1294
1295
1296   SUBROUTINE Cdn10_Lupkes2012( pcd )
1297      !!----------------------------------------------------------------------
1298      !!                      ***  ROUTINE  Cdn10_Lupkes2012  ***
1299      !!
1300      !! ** Purpose :    Recompute the neutral air-ice drag referenced at 10m
1301      !!                 to make it dependent on edges at leads, melt ponds and flows.
1302      !!                 After some approximations, this can be resumed to a dependency
1303      !!                 on ice concentration.
1304      !!
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
1316      !!                 Cdw(say 1.5e-3) for A=0, reaching 1.97e-3 for A~0.5
1317      !!                 and going down to Cdi(say 1.4e-3) for A=1
1318      !!
1319      !!                 It is theoretically applicable to all ice conditions (not only MIZ)
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      !!----------------------------------------------------------------------
1326      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pcd
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
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
1339
1340      ! ice-atm drag
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
1344   END SUBROUTINE Cdn10_Lupkes2012
1345
1346
1347   SUBROUTINE Cdn10_Lupkes2015( ptm_su, pslp, pcd, pch )
1348      !!----------------------------------------------------------------------
1349      !!                      ***  ROUTINE  Cdn10_Lupkes2015  ***
1350      !!
1351      !! ** pUrpose :    Alternative turbulent transfert coefficients formulation
1352      !!                 between sea-ice and atmosphere with distinct momentum
1353      !!                 and heat coefficients depending on sea-ice concentration
1354      !!                 and atmospheric stability (no meltponds effect for now).
1355      !!
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)
1359      !!                 to compute neutral transfert coefficients for both heat and
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      !
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
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
1403      zCdn_ice      = zCdn_skin_ice   ! Eq. 7
1404      !zCdn_ice     = 1.89e-3         ! old ECHAM5 value (cf Eq. 32)
1405
1406      ! Heat Neutral Transfert Coefficients
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
1409      ! Atmospheric and Surface Variables
1410      zst(:,:)     = sst_m(:,:) + rt0                                        ! convert SST from Celcius to Kelvin
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]
1413      !
1414      DO_2D( 0, 0, 0, 0 )
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
1457      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1.0_wp, pch, 'T', 1.0_wp )
1458      !
1459   END SUBROUTINE Cdn10_Lupkes2015
1460
1461#endif
1462
1463   !!======================================================================
1464END MODULE sbcblk
Note: See TracBrowser for help on using the repository browser.