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/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC – NEMO

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbcblk_algo_ice_lg15.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago

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