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_algo_ice_an05.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbcblk_algo_ice_an05.F90 @ 14021

Last change on this file since 14021 was 14021, checked in by laurent, 3 years ago

Caught up with trunk rev 14020...
File size: 19.1 KB
Line 
1MODULE sbcblk_algo_ice_an05
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_ice_an05  ***
4   !!       Computes turbulent components of surface fluxes over sea-ice
5   !!
6   !!   Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results.
7   !!   Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7
8   !!
9   !!   * bulk transfer coefficients C_D, C_E and C_H
10   !!   * air temp. and spec. hum. adjusted from zt (usually 2m) to zu (usually 10m) if needed
11   !!   * the "effective" bulk wind speed at zu: Ub (including gustiness contribution in unstable conditions)
12   !!   => all these are used in bulk formulas in sbcblk.F90
13   !!
14   !!       Routine turb_ice_an05 maintained and developed in AeroBulk
15   !!                     (https://github.com/brodeau/aerobulk/)
16   !!
17   !!            Author: Laurent Brodeau, Summer 2020
18   !!
19   !!----------------------------------------------------------------------
20   USE par_kind, ONLY: wp
21   USE par_oce,  ONLY: jpi, jpj, Nis0, Nie0, Njs0, Nje0, nn_hls, ntsi, ntsj, ntei, ntej
22   USE lib_mpp,  ONLY: ctl_stop         ! distribued memory computing library
23   USE phycst          ! physical constants
24   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
25
26   IMPLICIT NONE
27   PRIVATE
28
29   PUBLIC :: turb_ice_an05
30
31   INTEGER , PARAMETER ::   nbit = 8        ! number of itterations
32
33   !! * Substitutions
34#  include "do_loop_substitute.h90"
35   !!----------------------------------------------------------------------
36CONTAINS
37
38   SUBROUTINE turb_ice_an05( zt, zu, Ts_i, t_zt, qs_i, q_zt, U_zu,        &
39      &                      Cd_i, Ch_i, Ce_i, t_zu_i, q_zu_i,            &
40      &                      CdN, ChN, CeN, xz0, xu_star, xL, xUN10 )
41      !!----------------------------------------------------------------------
42      !!                      ***  ROUTINE  turb_ice_an05  ***
43      !!
44      !! ** Purpose :   Computes turbulent transfert coefficients of surface
45      !!                fluxes according to:
46      !!   Andreas, E.L., Jordan, R.E. & Makshtas, A.P. Parameterizing turbulent exchange over sea ice: the ice station weddell results.
47      !!   Boundary-Layer Meteorology 114, 439–460 (2005). https://doi.org/10.1007/s10546-004-1414-7
48      !!
49      !!           If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
50      !!           Returns the effective bulk wind speed at zu to be used in the bulk formulas
51      !!
52      !! INPUT :
53      !! -------
54      !!    *  zt   : height for temperature and spec. hum. of air            [m]
55      !!    *  zu   : height for wind speed (usually 10m)                     [m]
56      !!    *  Ts_i  : surface temperature of sea-ice                         [K]
57      !!    *  t_zt : potential air temperature at zt                         [K]
58      !!    *  qs_i  : saturation specific humidity at temp. Ts_i over ice    [kg/kg]
59      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
60      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
61      !!
62      !! OUTPUT :
63      !! --------
64      !!    *  Cd_i   : drag coefficient over sea-ice
65      !!    *  Ch_i   : sensible heat coefficient over sea-ice
66      !!    *  Ce_i   : sublimation coefficient over sea-ice
67      !!    *  t_zu_i : pot. air temp. adjusted at zu over sea-ice             [K]
68      !!    *  q_zu_i : spec. hum. of air adjusted at zu over sea-ice          [kg/kg]
69      !!
70      !! OPTIONAL OUTPUT:
71      !! ----------------
72      !!    * CdN     : neutral-stability drag coefficient
73      !!    * ChN     : neutral-stability sensible heat coefficient
74      !!    * CeN     : neutral-stability evaporation coefficient
75      !!    * xz0     : return the aerodynamic roughness length (integration constant for wind stress) [m]
76      !!    * xu_star : return u* the friction velocity                    [m/s]
77      !!    * xL      : return the Obukhov length                          [m]
78      !!    * xUN10   : neutral wind speed at 10m                          [m/s]
79      !!
80      !! ** Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
81      !!----------------------------------------------------------------------------------
82      REAL(wp), INTENT(in )                     :: zt    ! height for t_zt and q_zt                    [m]
83      REAL(wp), INTENT(in )                     :: zu    ! height for U_zu                             [m]
84      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: Ts_i  ! ice surface temperature                [Kelvin]
85      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt  ! potential air temperature              [Kelvin]
86      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: qs_i  ! sat. spec. hum. at ice/air interface    [kg/kg]
87      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt  ! spec. air humidity at zt               [kg/kg]
88      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu  ! relative wind module at zu                [m/s]
89      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Cd_i  ! drag coefficient over sea-ice
90      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ch_i  ! transfert coefficient for heat over ice
91      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ce_i  ! transfert coefficient for sublimation over ice
92      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: t_zu_i ! pot. air temp. adjusted at zu               [K]
93      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: q_zu_i ! spec. humidity adjusted at zu           [kg/kg]
94      !!----------------------------------------------------------------------------------
95      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CdN
96      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: ChN
97      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CeN
98      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xz0  ! Aerodynamic roughness length   [m]
99      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xu_star  ! u*, friction velocity
100      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xL  ! zeta (zu/L)
101      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xUN10  ! Neutral wind at zu
102      !!----------------------------------------------------------------------------------
103      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Ubzu
104      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp0, ztmp1, ztmp2      ! temporary stuff
105      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: z0, dt_zu, dq_zu
106      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: u_star, t_star, q_star
107      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: znu_a                    !: Nu_air = kinematic viscosity of air
108      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_u, zeta_t           ! stability parameter at height zu
109      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: z0tq
110      !!
111      INTEGER :: jit
112      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
113      !!
114      LOGICAL :: lreturn_cdn=.FALSE., lreturn_chn=.FALSE., lreturn_cen=.FALSE.
115      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE.
116      !!
117      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ice_an05@sbcblk_algo_ice_an05.f90'
118      !!----------------------------------------------------------------------------------
119      ALLOCATE (  Ubzu(jpi,jpj), u_star(jpi,jpj),  t_star(jpi,jpj),  q_star(jpi,jpj),  &
120         &      zeta_u(jpi,jpj),  dt_zu(jpi,jpj),   dq_zu(jpi,jpj),  &
121         &       znu_a(jpi,jpj),  ztmp1(jpi,jpj),   ztmp2(jpi,jpj),  &
122         &          z0(jpi,jpj),   z0tq(jpi,jpj,2), ztmp0(jpi,jpj)   )
123
124      lreturn_cdn   = PRESENT(CdN)
125      lreturn_chn   = PRESENT(ChN)
126      lreturn_cen   = PRESENT(CeN)
127      lreturn_z0    = PRESENT(xz0)
128      lreturn_ustar = PRESENT(xu_star)
129      lreturn_L     = PRESENT(xL)
130      lreturn_UN10  = PRESENT(xUN10)
131
132      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp )
133      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
134
135      !! Scalar wind speed cannot be below 0.2 m/s
136      Ubzu = MAX( U_zu, wspd_thrshld_ice )
137
138      !! First guess of temperature and humidity at height zu:
139      t_zu_i = MAX( t_zt ,   100._wp )   ! who knows what's given on masked-continental regions...
140      q_zu_i = MAX( q_zt , 0.1e-6_wp )   !               "
141
142      !! Air-Ice differences (and we don't want it to be 0!)
143      dt_zu = t_zu_i - Ts_i ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
144      dq_zu = q_zu_i - qs_i ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
145
146      znu_a = visc_air(t_zu_i) ! Air viscosity (m^2/s) at zt given from temperature in (K)
147
148      !! Very crude first guesses of z0:
149      z0 = 8.0E-4_wp
150
151      !! Crude first guess of turbulent scales
152      u_star = 0.035_wp*Ubzu*LOG( 10._wp/z0 )/LOG( zu/z0 )
153      z0 = rough_leng_m( u_star , znu_a )
154
155      DO jit = 1, 2
156         u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)) , 1.E-9 )
157         z0 = rough_leng_m( u_star , znu_a )
158      END DO
159
160      z0tq = rough_leng_tq( z0, u_star , znu_a )
161      t_star = dt_zu*vkarmn/(LOG(zu/z0tq(:,:,1)))
162      q_star = dq_zu*vkarmn/(LOG(zu/z0tq(:,:,2)))
163     
164     
165      !! ITERATION BLOCK
166      DO jit = 1, nbit
167
168         !!Inverse of Obukov length (1/L) :
169         ztmp0 = One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star)  ! 1/L == 1/[Obukhov length]
170         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...)
171
172         !! Stability parameters "zeta" :
173         zeta_u = zu*ztmp0
174         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
175         IF( .NOT. l_zt_equal_zu ) THEN
176            zeta_t = zt*ztmp0
177            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
178         END IF
179
180         !! Roughness lengthes z0, z0t, & z0q :
181         z0   = rough_leng_m (     u_star , znu_a )
182         z0tq = rough_leng_tq( z0, u_star , znu_a )
183
184         !! Turbulent scales at zu :
185         ztmp0   = psi_h_ice(zeta_u)
186         t_star =      dt_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,1)) - ztmp0)
187         q_star =      dq_zu*vkarmn/(LOG(zu) - LOG(z0tq(:,:,2)) - ztmp0)
188         u_star = MAX( Ubzu*vkarmn/(LOG(zu) - LOG(z0(:,:)) - psi_m_ice(zeta_u)) , 1.E-9 )
189
190         IF( .NOT. l_zt_equal_zu ) THEN
191            !! Re-updating temperature and humidity at zu if zt /= zu :
192            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_ice(zeta_t)
193            t_zu_i = t_zt - t_star/vkarmn*ztmp1
194            q_zu_i = q_zt - q_star/vkarmn*ztmp1
195            dt_zu = t_zu_i - Ts_i ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
196            dq_zu = q_zu_i - qs_i ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
197         END IF
198
199      END DO !DO jit = 1, nbit
200
201      ! compute transfer coefficients at zu :
202      ztmp0 = u_star/Ubzu
203      Cd_i   = ztmp0*ztmp0
204      Ch_i   = ztmp0*t_star/dt_zu
205      Ce_i   = ztmp0*q_star/dq_zu
206
207      IF( lreturn_cdn .OR. lreturn_chn .OR. lreturn_cen ) ztmp0 = 1._wp/LOG( zu/z0(:,:) )
208      IF( lreturn_cdn )   CdN = vkarmn2*ztmp0*ztmp0
209      IF( lreturn_chn )   ChN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,1))
210      IF( lreturn_cen )   CeN = vkarmn2*ztmp0/LOG(zu/z0tq(:,:,2))
211
212      IF( lreturn_z0 )    xz0     = z0
213      IF( lreturn_ustar ) xu_star = u_star
214      IF( lreturn_L )     xL      = 1./One_on_L(t_zu_i, q_zu_i, u_star, t_star, q_star)
215      IF( lreturn_UN10 )  xUN10   = u_star/vkarmn*LOG(10./z0)
216
217      DEALLOCATE ( Ubzu, u_star, t_star, q_star, zeta_u, dt_zu, dq_zu, z0, z0tq, znu_a, ztmp0, ztmp1, ztmp2 )
218      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
219
220   END SUBROUTINE turb_ice_an05
221
222
223
224   FUNCTION rough_leng_m( pus , pnua )
225      !!----------------------------------------------------------------------------------
226      !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 19)
227      !!
228      !! Author: L. Brodeau, January 2020 / AeroBulk  (https://github.com/brodeau/aerobulk/)
229      !!----------------------------------------------------------------------------------
230      REAL(wp), DIMENSION(jpi,jpj) :: rough_leng_m      ! roughness length over sea-ice [m]
231      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus   ! u* = friction velocity    [m/s]
232      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua  ! kinematic viscosity of air [m^2/s]
233      !!
234      INTEGER  :: ji, jj    ! dummy loop indices
235      REAL(wp) :: zus, zz
236      !!----------------------------------------------------------------------------------
237      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
238            zus = MAX( pus(ji,jj) , 1.E-9_wp )
239
240            zz = (zus - 0.18_wp) / 0.1_wp
241
242            rough_leng_m(ji,jj) = 0.135*pnua(ji,jj)/zus + 0.035*zus*zus/grav*( 5.*EXP(-zz*zz) + 1._wp ) ! Eq.(19) Andreas et al., 2005
243      END_2D
244      !!
245   END FUNCTION rough_leng_m
246
247   FUNCTION rough_leng_tq( pz0, pus , pnua )
248      !!----------------------------------------------------------------------------------
249      !! Computes the roughness length of sea-ice according to Andreas et al. 2005, (eq. 22)
250      !!    => which still relies on Andreas 1987 !
251      !!
252      !! Author: L. Brodeau, January 2020 / AeroBulk  (https://github.com/brodeau/aerobulk/)
253      !!----------------------------------------------------------------------------------
254      REAL(wp), DIMENSION(jpi,jpj,2)           :: rough_leng_tq     ! temp.,hum. roughness lengthes over sea-ice [m]
255      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0   ! roughness length            [m]
256      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus   ! u* = friction velocity    [m/s]
257      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pnua  ! kinematic viscosity of air [m^2/s]
258      !!
259      INTEGER  :: ji, jj    ! dummy loop indices
260      REAL(wp) :: zz0, zus, zre, zsmoot, ztrans, zrough
261      REAL(wp) :: zb0, zb1, zb2, zlog, zlog2, zlog_z0s_on_z0
262      !!----------------------------------------------------------------------------------
263      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
264            zz0 = pz0(ji,jj)
265            zus = MAX( pus(ji,jj) , 1.E-9_wp )
266            zre = MAX( zus*zz0/pnua(ji,jj) , 0._wp ) ! Roughness Reynolds number
267
268            !! *** TABLE 1 of Andreas et al. 2005 ***
269            !! Smooth flow condition (R* <= 0.135):
270            zsmoot = 0.5_wp + SIGN( 0.5_wp, (0.135_wp   - zre) ) ! zre <= 0.135: zsmoot==1 ; otherwize: zsmoot==0
271            !! Transition (0.135 < R* < 2.5):
272            ztrans = 0.5_wp + SIGN( 0.5_wp, (2.49999_wp - zre) ) - zsmoot
273            !! Rough ( R* > 2.5):
274            zrough = 0.5_wp + SIGN( 0.5_wp, (zre - 2.5_wp) )
275
276            IF( (zsmoot+ztrans+zrough > 1.001_wp).OR.(zsmoot+ztrans+zrough < 0.999_wp) ) &
277               CALL ctl_stop( ' rough_leng_tq@mod_blk_ice_an05.f90 => something wrong with zsmoot, ztrans, zrough!' )
278
279            zlog  = LOG(zre)
280            zlog2 = zlog*zlog
281
282            !! z0t:
283            zb0 = zsmoot*1.25_wp + ztrans*0.149_wp + zrough*0.317_wp
284            zb1 =                - ztrans*0.550_wp - zrough*0.565_wp
285            zb2 =                                  - zrough*0.183_wp
286            zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2
287            rough_leng_tq(ji,jj,1) = zz0 * EXP( zlog_z0s_on_z0 )
288
289            !! z0q:
290            zb0 = zsmoot*1.61_wp + ztrans*0.351_wp + zrough*0.396_wp
291            zb1 =                - ztrans*0.628_wp - zrough*0.512_wp
292            zb2 =                                  - zrough*0.180_wp
293            zlog = LOG(zre)
294            zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2
295            rough_leng_tq(ji,jj,2) = zz0 * EXP( zlog_z0s_on_z0 )
296
297      END_2D
298      !!
299   END FUNCTION rough_leng_tq
300
301
302
303   FUNCTION psi_m_ice( pzeta )
304      !!----------------------------------------------------------------------------------
305      !! ** Purpose: compute the universal profile stability function for momentum
306      !!
307      !!
308      !!     Andreas et al 2005 == Jordan et al. 1999
309      !!
310      !!     Psi:
311      !!     Unstable => Paulson 1970
312      !!     Stable   => Holtslag & De Bruin 1988
313      !!
314      !!             pzeta : stability paramenter, z/L where z is altitude
315      !!                     measurement and L is M-O length
316      !!
317      !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
318      !!----------------------------------------------------------------------------------
319      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ice
320      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
321      !
322      INTEGER  ::   ji, jj    ! dummy loop indices
323      REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab
324      !!----------------------------------------------------------------------------------
325      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )            !
326            zta = pzeta(ji,jj)
327            !
328            ! Unstable stratification:
329            zx = ABS(1._wp - 16._wp*zta)**.25              !  (16 here, not 15!)
330
331            zpsi_u =      LOG( (1._wp + zx*zx)/2. ) &  ! Eq.(30) Jordan et al. 1999
332               &     + 2.*LOG( (1._wp + zx       )/2. ) &
333               &    - 2.*ATAN( zx ) + 0.5*rpi
334
335            ! Stable stratification:
336            zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp )  ! Eq.(33) Jordan et al. 1999
337
338            !! Combine:
339            zstab = 0.5 + SIGN(0.5_wp, zta)
340            psi_m_ice(ji,jj) = (1._wp - zstab) * zpsi_u   & ! Unstable (zta < 0)
341               &                   + zstab     * zpsi_s     ! Stable (zta > 0)
342            !
343      END_2D
344   END FUNCTION psi_m_ice
345
346
347   FUNCTION psi_h_ice( pzeta )
348      !!----------------------------------------------------------------------------------
349      !! ** Purpose: compute the universal profile stability function for
350      !!             temperature and humidity
351      !!
352      !!
353      !!     Andreas et al 2005 == Jordan et al. 1999
354      !!
355      !!     Psi:
356      !!     Unstable => Paulson 1970
357      !!     Stable   => Holtslag & De Bruin 1988
358      !!
359      !!             pzeta : stability paramenter, z/L where z is altitude
360      !!                     measurement and L is M-O length
361      !!
362      !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
363      !!----------------------------------------------------------------------------------
364      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ice
365      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
366      !
367      INTEGER  ::   ji, jj    ! dummy loop indices
368      REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab
369      !!----------------------------------------------------------------------------------
370      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )            !
371            zta = pzeta(ji,jj)
372            !
373            ! Unstable stratification:
374            zx = ABS(1._wp - 16._wp*zta)**.25              !  (16 here, not 15!)
375
376            zpsi_u =   2.*LOG( (1._wp + zx*zx)/2. )  ! Eq.(31) Jordan et al. 1999
377
378            ! Stable stratification (identical to Psi_m!):
379            zpsi_s = - ( 0.7_wp*zta + 0.75_wp*(zta - 14.3_wp)*EXP( -0.35*zta) + 10.7_wp )  ! Eq.(33) Jordan et al. 1999
380
381            !! Combine:
382            zstab = 0.5 + SIGN(0.5_wp, zta)
383            psi_h_ice(ji,jj) = (1._wp - zstab) * zpsi_u   & ! Unstable (zta < 0)
384               &                   + zstab     * zpsi_s     ! Stable (zta > 0)
385            !
386      END_2D
387   END FUNCTION psi_h_ice
388
389   !!======================================================================
390END MODULE sbcblk_algo_ice_an05
Note: See TracBrowser for help on using the repository browser.