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_lg15.F90 in NEMO/trunk/src/OCE/SBC – NEMO

source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_ice_lg15.F90 @ 14072

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

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369
File size: 15.1 KB
Line 
1MODULE sbcblk_algo_ice_lg15
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_ice_lg15  ***
4   !!       Computes turbulent components of surface fluxes over sea-ice
5   !!
6   !!
7   !!  Lüpkes, C., and Gryanik, V. M. ( 2015), A stability‐dependent parametrization
8   !!  of transfer coefficients for momentum and heat over polar sea ice to be used in climate models,
9   !!  J. Geophys. Res. Atmos., 120, 552– 581, doi:10.1002/2014JD022418.
10   !!
11   !!       => Despite the fact that the sea-ice concentration (frice) must be provided,
12   !!          only transfer coefficients, and air temp. + hum. height adjustement
13   !!          over ice are returned/performed.
14   !!        ==> 'frice' is only here to estimate the form drag caused by sea-ice...
15   !!
16   !!       Routine turb_ice_lg15 maintained and developed in AeroBulk
17   !!                     (https://github.com/brodeau/aerobulk/)
18   !!
19   !!            Author: Laurent Brodeau, Summer 2020
20   !!
21   !!----------------------------------------------------------------------
22   USE par_kind, ONLY: wp
23   USE par_oce,  ONLY: jpi, jpj
24   USE phycst          ! physical constants
25   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
26   USE sbcblk_algo_ice_cdn
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC :: turb_ice_lg15
32
33   REAL(wp), PARAMETER ::   ralpha_0  = 0.2_wp     ! (Eq.12) (ECHAM6 value)
34
35   !! To be namelist parameters in NEMO:
36   REAL(wp), PARAMETER :: rz0_i_s_0  = 0.69e-3_wp  !           Eq. 43 [m]
37   REAL(wp), PARAMETER :: rz0_i_f_0  = 4.54e-4_wp  ! bottom p.562 MIZ [m]
38
39   LOGICAL,  PARAMETER :: l_add_form_drag = .TRUE.
40   LOGICAL,  PARAMETER :: l_use_pond_info = .FALSE.
41   LOGICAL,  PARAMETER :: l_dbg_print     = .FALSE.
42
43   INTEGER , PARAMETER ::   nbit = 8        ! number of itterations
44
45   !!----------------------------------------------------------------------
46CONTAINS
47
48   SUBROUTINE turb_ice_lg15( zt, zu, Ts_i, t_zt, qs_i, q_zt, U_zu, frice, &
49      &                      Cd_i, Ch_i, Ce_i, t_zu_i, q_zu_i,            &
50      &                      CdN, ChN, CeN, xz0, xu_star, xL, xUN10 )
51      !!----------------------------------------------------------------------
52      !!                      ***  ROUTINE  turb_ice_lg15  ***
53      !!
54      !! ** Purpose :   Computes turbulent transfert coefficients of surface
55      !!                fluxes according to:
56      !!            Lüpkes, C., and Gryanik, V. M. ( 2015), A stability‐dependent
57      !!            parametrization of transfer coefficients for momentum and heat
58      !!            over polar sea ice to be used in climate models,
59      !!            J. Geophys. Res. Atmos., 120, 552– 581, doi:10.1002/2014JD022418.
60      !!
61      !!           If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
62      !!           Returns the effective bulk wind speed at zu to be used in the bulk formulas
63      !!
64      !! INPUT :
65      !! -------
66      !!    *  zt   : height for temperature and spec. hum. of air            [m]
67      !!    *  zu   : height for wind speed (usually 10m)                     [m]
68      !!    *  Ts_i  : surface temperature of sea-ice                         [K]
69      !!    *  t_zt : potential air temperature at zt                         [K]
70      !!    *  qs_i  : saturation specific humidity at temp. Ts_i over ice    [kg/kg]
71      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
72      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
73      !!    * frice : sea-ice concentration        (fraction)
74      !!
75      !! OUTPUT :
76      !! --------
77      !!    *  Cd_i   : drag coefficient over sea-ice
78      !!    *  Ch_i   : sensible heat coefficient over sea-ice
79      !!    *  Ce_i   : sublimation coefficient over sea-ice
80      !!    *  t_zu_i : pot. air temp. adjusted at zu over sea-ice             [K]
81      !!    *  q_zu_i : spec. hum. of air adjusted at zu over sea-ice          [kg/kg]
82      !!
83      !! OPTIONAL OUTPUT:
84      !! ----------------
85      !!    * CdN     : neutral-stability drag coefficient
86      !!    * ChN     : neutral-stability sensible heat coefficient
87      !!    * CeN     : neutral-stability evaporation coefficient
88      !!    * xz0     : return the aerodynamic roughness length (integration constant for wind stress) [m]
89      !!    * xu_star : return u* the friction velocity                    [m/s]
90      !!    * xL      : return the Obukhov length                          [m]
91      !!    * xUN10   : neutral wind speed at 10m                          [m/s]
92      !!
93      !! ** Author: L. Brodeau, January 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
94      !!----------------------------------------------------------------------------------
95      REAL(wp), INTENT(in )                     :: zt    ! height for t_zt and q_zt                    [m]
96      REAL(wp), INTENT(in )                     :: zu    ! height for U_zu                             [m]
97      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: Ts_i  ! ice surface temperature                [Kelvin]
98      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: t_zt  ! potential air temperature              [Kelvin]
99      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: qs_i  ! sat. spec. hum. at ice/air interface    [kg/kg]
100      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: q_zt  ! spec. air humidity at zt               [kg/kg]
101      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: U_zu  ! relative wind module at zu                [m/s]
102      REAL(wp), INTENT(in ), DIMENSION(jpi,jpj) :: frice ! sea-ice concentration        (fraction)
103      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Cd_i  ! drag coefficient over sea-ice
104      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ch_i  ! transfert coefficient for heat over ice
105      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: Ce_i  ! transfert coefficient for sublimation over ice
106      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: t_zu_i ! pot. air temp. adjusted at zu               [K]
107      REAL(wp), INTENT(out), DIMENSION(jpi,jpj) :: q_zu_i ! spec. humidity adjusted at zu           [kg/kg]
108      !!----------------------------------------------------------------------------------
109      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CdN
110      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: ChN
111      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: CeN
112      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xz0  ! Aerodynamic roughness length   [m]
113      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xu_star  ! u*, friction velocity
114      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xL  ! zeta (zu/L)
115      REAL(wp), INTENT(out), DIMENSION(jpi,jpj), OPTIONAL :: xUN10  ! Neutral wind at zu
116      !!----------------------------------------------------------------------------------
117      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: Ubzu
118      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: ztmp1, ztmp2      ! temporary stuff
119      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: dt_zu, dq_zu
120      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zz0_s, zz0_f, RiB ! third dimensions (size=2):
121      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zCdN_s, zChN_s, zCdN_f, zChN_f
122      !!
123      INTEGER :: jit
124      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
125      !!
126      LOGICAL :: lreturn_cdn=.FALSE., lreturn_chn=.FALSE., lreturn_cen=.FALSE.
127      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE.
128      !!
129      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ice_lg15@sbcblk_algo_ice_lg15.f90'
130      !!----------------------------------------------------------------------------------
131      ALLOCATE ( Ubzu(jpi,jpj) )
132      ALLOCATE ( ztmp1(jpi,jpj),  ztmp2(jpi,jpj) )
133      ALLOCATE ( dt_zu(jpi,jpj),  dq_zu(jpi,jpj) )
134      ALLOCATE ( zz0_s(jpi,jpj),  zz0_f(jpi,jpj),    RiB(jpi,jpj), &
135         &      zCdN_s(jpi,jpj), zChN_s(jpi,jpj), zCdN_f(jpi,jpj), zChN_f(jpi,jpj) )
136
137      lreturn_cdn   = PRESENT(CdN)
138      lreturn_chn   = PRESENT(ChN)
139      lreturn_cen   = PRESENT(CeN)
140      lreturn_z0    = PRESENT(xz0)
141      lreturn_ustar = PRESENT(xu_star)
142      lreturn_L     = PRESENT(xL)
143      lreturn_UN10  = PRESENT(xUN10)
144
145      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp )
146
147      !! Scalar wind speed cannot be below 0.2 m/s
148      Ubzu = MAX( U_zu, wspd_thrshld_ice )
149
150      !! First guess of temperature and humidity at height zu:
151      t_zu_i = MAX( t_zt ,   100._wp )   ! who knows what's given on masked-continental regions...
152      q_zu_i = MAX( q_zt , 0.1e-6_wp )   !               "
153
154      !! Air-Ice & Air-Sea differences (and we don't want them to be 0!)
155      dt_zu = t_zu_i - Ts_i ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
156      dq_zu = q_zu_i - qs_i ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
157
158      !! Very crude first guess:
159      Cd_i(:,:) = 1.4e-3_wp
160      Ch_i(:,:) = 1.4e-3_wp
161      Ce_i(:,:) = 1.4e-3_wp
162     
163      !! For skin drag :
164      zz0_s(:,:)  = rz0_i_s_0        !#LB/RFI! ! Room for improvement. We use the same z0_skin everywhere (= rz0_i_s_0)...
165      zCdN_s(:,:) = Cd_from_z0( zu, zz0_s(:,:) )
166      zChN_s(:,:) = vkarmn2 / ( LOG( zu / zz0_s(:,:) ) * LOG( zu / (ralpha_0*zz0_s(:,:)) ) )     ! (Eq.11,12)  [ "" ]
167     
168      !! For form drag in MIZ:
169      zz0_f(:,:)  = 0._wp
170      zCdN_f(:,:) = 0._wp
171      zChN_f(:,:) = 0._wp
172      IF ( l_add_form_drag ) THEN
173         zz0_f(:,:)  = rz0_i_f_0        !#LB/RFI! ! Room for improvement. We use the same z0_form everywhere !!!
174         zCdN_f(:,:) = CdN_f_LG15_light( zu, frice(:,:), zz0_f(:,:) )
175         zChN_f(:,:) = zCdN_f(:,:)/( 1._wp + LOG(1._wp/ralpha_0)/vkarmn*SQRT(zCdN_f(:,:)) ) ! (Eq.60,61)   [ "" ]
176      END IF
177
178      !! Some other first guess values, needed to compute wind at zt:
179      Cd_i(:,:) = zCdN_s(:,:) + zCdN_f(:,:)
180      Ch_i(:,:) = zChN_s(:,:) + zChN_f(:,:)
181      RiB(:,:) = Ri_bulk( zt, Ts_i(:,:), t_zt(:,:), qs_i(:,:), q_zt(:,:), Ubzu(:,:) )  ! over ice (index=1)
182
183
184      !! ITERATION BLOCK
185      DO jit = 1, nbit
186
187         IF(l_dbg_print) PRINT *, 'LOLO: LOOP #', INT(jit,1)
188         IF(l_dbg_print) PRINT *, 'LOLO: theta_zu, Ts_i, Ubzu =', REAL(t_zu_i(:,:),4), REAL(Ts_i(:,:),4), REAL(Ubzu(:,:),4)
189         IF(l_dbg_print) PRINT *, 'LOLO:     q_zu =', REAL(q_zu_i(:,:),4)
190         IF(l_dbg_print) PRINT *, 'LOLO:  CdN_s, zCdN_f   =', REAL(zCdN_s(:,:),4), REAL(zCdN_f(:,:),4)
191
192
193         !! Bulk Richardson Number
194         !! ======================
195         !! PROBLEM: when computed at z=zu, with adjusted theta and q, it is numerically unstable in some rare events (unstable)
196         !!          => fix: compute RiB at zt, with ajusted wind at zt... => seems to be more stable
197         IF( .NOT. l_zt_equal_zu ) THEN
198            ! U_zt = U_zu + u_star/vkarmn*(LOG(zt/zu) + psi_m_coare(zu/L) - psi_m_coare(zt/L))
199            ztmp1(:,:) = zCdN_s(:,:) + zCdN_f(:,:)    ! total neutral drag coeff!
200            ztmp2(:,:) = zz0_s(:,:) + zz0_f(:,:)      ! total roughness length z0
201            ztmp1 = LOG(zt/zu) + f_h_louis( zu, RiB(:,:), ztmp1(:,:), ztmp2(:,:) ) &
202               &               - f_h_louis( zt, RiB(:,:), ztmp1(:,:), ztmp2(:,:) )
203            ztmp2(:,:) = MAX( Ubzu(:,:) + (SQRT(Cd_i(:,:))*Ubzu)*ztmp1 , wspd_thrshld_ice ) ! wind at zt ( SQRT(Cd_i(:,:))*Ubzu == u* !)
204            ztmp2(:,:) = MIN( ztmp2(:,:) , Ubzu(:,:) )
205            IF(l_dbg_print) PRINT *, 'LOLO: ADJUSTED WIND AT ZT =', ztmp2
206         ELSE
207            ztmp2(:,:) = Ubzu(:,:)
208         END IF
209         RiB(:,:) = Ri_bulk( zt, Ts_i(:,:), t_zt(:,:), qs_i(:,:), q_zt(:,:), ztmp2(:,:) )  ! over ice (index=1)
210         IF(l_dbg_print) PRINT *, 'LOLO: RiB_zt =', RiB(:,:)
211
212
213         ! Momentum and Heat transfer coefficients WITHOUT FORM DRAG / (Eq.6) and (Eq.10):
214         Cd_i(:,:) = zCdN_s(:,:) * f_m_louis( zu, RiB(:,:), zCdN_s(:,:), zz0_s(:,:) ) ! (Eq.6)
215         Ch_i(:,:) = zChN_s(:,:) * f_h_louis( zu, RiB(:,:), zCdN_s(:,:), zz0_s(:,:) ) ! (Eq.10) / LOLO: why "zCdN_s" (ztmp1) and not "zChn" ???
216         IF(l_dbg_print) PRINT *, 'LOLO: f_m_louis_s =', f_m_louis( zu, RiB(:,:), zCdN_s(:,:), zz0_s(:,:) )
217         IF(l_dbg_print) PRINT *, 'LOLO: f_h_louis_s =', f_h_louis( zu, RiB(:,:), zCdN_s(:,:), zz0_s(:,:) )
218         IF(l_dbg_print) PRINT *, 'LOLO: Cd / skin only / ice   =', REAL(Cd_i(:,:),4)
219
220         
221         IF ( l_add_form_drag ) THEN
222            !! Form-drag-related NEUTRAL momentum and Heat transfer coefficients:
223            !!   MIZ:
224            Cd_i(:,:) = Cd_i(:,:) + zCdN_f(:,:) * f_m_louis( zu, RiB(:,:), zCdN_f(:,:), zz0_f(:,:) ) ! (Eq.6)
225            Ch_i(:,:) = Ch_i(:,:) + zChN_f(:,:) * f_h_louis( zu, RiB(:,:), zCdN_f(:,:), zz0_f(:,:) ) ! (Eq.10) / LOLO: why "zCdN_f" and not "zChn" ???
226            IF(l_dbg_print) PRINT *, 'LOLO: f_m_louis_f =', f_m_louis( zu, RiB(:,:), zCdN_f(:,:), zz0_f(:,:) )
227            IF(l_dbg_print) PRINT *, 'LOLO: f_h_louis_f =', f_h_louis( zu, RiB(:,:), zCdN_f(:,:), zz0_f(:,:) )
228
229            IF(l_dbg_print) PRINT *, 'LOLO: Cd / form only / ice   =', REAL(zCdN_f(:,:) * f_m_louis( zu, RiB(:,:), zCdN_f(:,:), zz0_f(:,:) ),4)
230
231         END IF
232
233         IF(l_dbg_print) PRINT *, 'LOLO: Cd, Ch / TOTAL / ice   =', REAL(Cd_i(:,:),4), REAL(Ch_i(:,:),4)
234
235
236         !! Adjusting temperature and humidity from zt to zu:
237         IF( .NOT. l_zt_equal_zu ) THEN
238
239            !! Over ice:
240            ztmp1(:,:) = zCdN_s(:,:) + zCdN_f(:,:)    ! total neutral drag coeff!
241            ztmp2(:,:) = zz0_s(:,:) + zz0_f(:,:)      ! total roughness length z0
242            ztmp1 = LOG(zt/zu) + f_h_louis( zu, RiB(:,:), ztmp1(:,:), ztmp2(:,:) ) &
243               &               - f_h_louis( zt, RiB(:,:), ztmp1(:,:), ztmp2(:,:) )
244            ztmp2 = 1._wp/SQRT(Cd_i(:,:))
245
246            t_zu_i(:,:) = t_zt - (Ch_i(:,:) * dt_zu(:,:) * ztmp2) / vkarmn * ztmp1   ! t_star = Ch * dt_zu / SQRT(Cd)
247            q_zu_i(:,:) = q_zt - (Ch_i(:,:) * dq_zu(:,:) * ztmp2) / vkarmn * ztmp1   ! q_star = Ce * dq_zu / SQRT(Cd)
248            q_zu_i(:,:) = MAX(0._wp, q_zu_i(:,:))
249
250            dt_zu(:,:) = t_zu_i(:,:) - Ts_i
251            dq_zu(:,:) = q_zu_i(:,:) - qs_i
252
253            dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
254            dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
255         END IF
256
257         IF(l_dbg_print) PRINT *, ''!LOLO
258
259      END DO !DO jit = 1, nbit
260
261      Ce_i(:,:)   =  Ch_i(:,:)
262
263      IF( lreturn_cdn )   CdN = zCdN_s(:,:)+zCdN_f(:,:)
264      IF( lreturn_chn )   ChN = zChN_s(:,:)+zChN_f(:,:)
265      IF( lreturn_cen )   CeN = zChN_s(:,:)+zChN_f(:,:)
266
267      IF( lreturn_z0 ) xz0   = z0_from_Cd( zu, zCdN_s(:,:)+zCdN_f(:,:) )
268
269      IF( lreturn_ustar ) xu_star = SQRT(Cd_i) * Ubzu
270
271      IF( lreturn_L ) THEN
272         ztmp1 = SQRT(Cd_i)
273         xL    = 1./One_on_L( t_zu_i, q_zu_i, ztmp1*Ubzu, Ch_i*dt_zu(:,:)/ztmp1, Ce_i*dq_zu(:,:)/ztmp1 )
274      END IF
275
276      IF( lreturn_UN10 ) THEN
277         ztmp1 = zCdN_s(:,:) + zCdN_f(:,:)  ! => CdN
278         xUN10 = SQRT(Cd_i) * Ubzu/vkarmn * LOG( 10._wp / z0_from_Cd(zu, ztmp1) )
279      END IF
280
281      DEALLOCATE ( Ubzu )
282      DEALLOCATE ( ztmp1, ztmp2 )
283      DEALLOCATE ( dt_zu, dq_zu )
284      DEALLOCATE ( zz0_s, zz0_f, RiB, zCdN_s, zChN_s, zCdN_f, zChN_f )
285
286   END SUBROUTINE turb_ice_lg15
287
288   !!======================================================================
289END MODULE sbcblk_algo_ice_lg15
Note: See TracBrowser for help on using the repository browser.