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/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ice_an05.F90 @ 15149

Last change on this file since 15149 was 15081, checked in by clem, 3 years ago

debug algo andreas2005 in the trunk
File size: 18.9 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            zsmoot = 0._wp ; ztrans = 0._wp ; zrough = 0._wp
270            IF    ( zre <= 0.135_wp ) THEN   ! Smooth flow condition (R* <= 0.135):
271               zsmoot = 1._wp
272            ELSEIF( zre < 2.5_wp )    THEN   ! Transition (0.135 < R* < 2.5)
273               ztrans = 1._wp
274            ELSE                             ! Rough ( R* > 2.5)
275               zrough = 1._wp
276            ENDIF
277               
278            zlog  = LOG(zre)
279            zlog2 = zlog*zlog
280
281            !! z0t:
282            zb0 = zsmoot*1.25_wp + ztrans*0.149_wp + zrough*0.317_wp
283            zb1 =                - ztrans*0.550_wp - zrough*0.565_wp
284            zb2 =                                  - zrough*0.183_wp
285            zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2
286            rough_leng_tq(ji,jj,1) = zz0 * EXP( zlog_z0s_on_z0 )
287
288            !! z0q:
289            zb0 = zsmoot*1.61_wp + ztrans*0.351_wp + zrough*0.396_wp
290            zb1 =                - ztrans*0.628_wp - zrough*0.512_wp
291            zb2 =                                  - zrough*0.180_wp
292            zlog = LOG(zre)
293            zlog_z0s_on_z0 = zb0 + zb1*zlog + zb2*zlog2
294            rough_leng_tq(ji,jj,2) = zz0 * EXP( zlog_z0s_on_z0 )
295
296      END_2D
297      !!
298   END FUNCTION rough_leng_tq
299
300
301
302   FUNCTION psi_m_ice( pzeta )
303      !!----------------------------------------------------------------------------------
304      !! ** Purpose: compute the universal profile stability function for momentum
305      !!
306      !!
307      !!     Andreas et al 2005 == Jordan et al. 1999
308      !!
309      !!     Psi:
310      !!     Unstable => Paulson 1970
311      !!     Stable   => Holtslag & De Bruin 1988
312      !!
313      !!             pzeta : stability paramenter, z/L where z is altitude
314      !!                     measurement and L is M-O length
315      !!
316      !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
317      !!----------------------------------------------------------------------------------
318      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_ice
319      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
320      !
321      INTEGER  ::   ji, jj    ! dummy loop indices
322      REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab
323      !!----------------------------------------------------------------------------------
324      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )            !
325            zta = pzeta(ji,jj)
326            !
327            ! Unstable stratification:
328            zx = ABS(1._wp - 16._wp*zta)**.25              !  (16 here, not 15!)
329
330            zpsi_u =      LOG( (1._wp + zx*zx)/2. ) &  ! Eq.(30) Jordan et al. 1999
331               &     + 2.*LOG( (1._wp + zx       )/2. ) &
332               &    - 2.*ATAN( zx ) + 0.5*rpi
333
334            ! Stable stratification:
335            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
336
337            !! Combine:
338            zstab = 0.5 + SIGN(0.5_wp, zta)
339            psi_m_ice(ji,jj) = (1._wp - zstab) * zpsi_u   & ! Unstable (zta < 0)
340               &                   + zstab     * zpsi_s     ! Stable (zta > 0)
341            !
342      END_2D
343   END FUNCTION psi_m_ice
344
345
346   FUNCTION psi_h_ice( pzeta )
347      !!----------------------------------------------------------------------------------
348      !! ** Purpose: compute the universal profile stability function for
349      !!             temperature and humidity
350      !!
351      !!
352      !!     Andreas et al 2005 == Jordan et al. 1999
353      !!
354      !!     Psi:
355      !!     Unstable => Paulson 1970
356      !!     Stable   => Holtslag & De Bruin 1988
357      !!
358      !!             pzeta : stability paramenter, z/L where z is altitude
359      !!                     measurement and L is M-O length
360      !!
361      !! ** Author: L. Brodeau, 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
362      !!----------------------------------------------------------------------------------
363      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_ice
364      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
365      !
366      INTEGER  ::   ji, jj    ! dummy loop indices
367      REAL(wp) :: zta, zx, zpsi_u, zpsi_s, zstab
368      !!----------------------------------------------------------------------------------
369      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )            !
370            zta = pzeta(ji,jj)
371            !
372            ! Unstable stratification:
373            zx = ABS(1._wp - 16._wp*zta)**.25              !  (16 here, not 15!)
374
375            zpsi_u =   2.*LOG( (1._wp + zx*zx)/2. )  ! Eq.(31) Jordan et al. 1999
376
377            ! Stable stratification (identical to Psi_m!):
378            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
379
380            !! Combine:
381            zstab = 0.5 + SIGN(0.5_wp, zta)
382            psi_h_ice(ji,jj) = (1._wp - zstab) * zpsi_u   & ! Unstable (zta < 0)
383               &                   + zstab     * zpsi_s     ! Stable (zta > 0)
384            !
385      END_2D
386   END FUNCTION psi_h_ice
387
388   !!======================================================================
389END MODULE sbcblk_algo_ice_an05
Note: See TracBrowser for help on using the repository browser.