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

source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

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

Merging branch "2020/dev_r13648_ASINTER-04_laurent_bulk_ice", ticket #2369

File size: 26.1 KB
Line 
1MODULE sbcblk_algo_coare3p6
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p6  ***
4   !!
5   !!                After Fairall et al 2018 & Edson et al 2013
6   !! Computes:
7   !!   * bulk transfer coefficients C_D, C_E and C_H
8   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
9   !!   * the effective bulk wind speed at 10m Ubzu
10   !!   => all these are used in bulk formulas in sbcblk.F90
11   !!
12   !!       Routine turb_coare3p6 maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk)
14   !!
15   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
16   !!----------------------------------------------------------------------
17   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
18   !!            4.2  !  2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   turb_coare3p6  : computes the bulk turbulent transfer coefficients
23   !!                   adjusts t_air and q_air from zt to zu m
24   !!                   returns the effective bulk wind speed at 10m
25   !!----------------------------------------------------------------------
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE lib_mpp,        ONLY: ctl_stop         ! distribued memory computing library
29   USE in_out_manager, ONLY: nit000  ! I/O manager
30   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
31   USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC :: SBCBLK_ALGO_COARE3P6_INIT, TURB_COARE3P6
37   !! * Substitutions
38#  include "do_loop_substitute.h90"
39
40   !! COARE own values for given constants:
41   REAL(wp), PARAMETER :: zi0   = 600._wp     ! scale height of the atmospheric boundary layer...
42   REAL(wp), PARAMETER :: Beta0 =  1.2_wp     ! gustiness parameter
43   REAL(wp), PARAMETER :: zeta_abs_max = 50._wp
44
45   !!----------------------------------------------------------------------
46CONTAINS
47
48
49   SUBROUTINE sbcblk_algo_coare3p6_init(l_use_cs, l_use_wl)
50      !!---------------------------------------------------------------------
51      !!                  ***  FUNCTION sbcblk_algo_coare3p6_init  ***
52      !!
53      !! INPUT :
54      !! -------
55      !!    * l_use_cs : use the cool-skin parameterization
56      !!    * l_use_wl : use the warm-layer parameterization
57      !!---------------------------------------------------------------------
58      LOGICAL , INTENT(in) ::   l_use_cs ! use the cool-skin parameterization
59      LOGICAL , INTENT(in) ::   l_use_wl ! use the warm-layer parameterization
60      INTEGER :: ierr
61      !!---------------------------------------------------------------------
62      IF( l_use_wl ) THEN
63         ierr = 0
64         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr )
65         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P6_INIT => allocation of Tau_ac, Qnt_ac, dT_wl & Hz_wl failed!' )
66         Tau_ac(:,:) = 0._wp
67         Qnt_ac(:,:) = 0._wp
68         dT_wl(:,:)  = 0._wp
69         Hz_wl(:,:)  = Hwl_max
70      ENDIF
71      IF( l_use_cs ) THEN
72         ierr = 0
73         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
74         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P6_INIT => allocation of dT_cs failed!' )
75         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
76      ENDIF
77   END SUBROUTINE sbcblk_algo_coare3p6_init
78
79
80
81   SUBROUTINE turb_coare3p6( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl, &
82      &                      Cd, Ch, Ce, t_zu, q_zu, Ubzu,                               &
83      &                      nb_iter, Cdn, Chn, Cen,                                     & ! optional output
84      &                      Qsw, rad_lw, slp, pdT_cs,                                   & ! optionals for cool-skin (and warm-layer)
85      &                      pdT_wl, pHz_wl )                                              ! optionals for warm-layer only
86      !!----------------------------------------------------------------------
87      !!                      ***  ROUTINE  turb_coare3p6  ***
88      !!
89      !! ** Purpose :   Computes turbulent transfert coefficients of surface
90      !!                fluxes according to Fairall et al. (2003)
91      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
92      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
93      !!
94      !!                Applies the cool-skin warm-layer correction of the SST to T_s
95      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
96      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
97      !!                are provided as (optional) arguments!
98      !!
99      !! INPUT :
100      !! -------
101      !!    *  kt   : current time step (starts at 1)
102      !!    *  zt   : height for temperature and spec. hum. of air            [m]
103      !!    *  zu   : height for wind speed (usually 10m)                     [m]
104      !!    *  t_zt : potential air temperature at zt                         [K]
105      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
106      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
107      !!    * l_use_cs : use the cool-skin parameterization
108      !!    * l_use_wl : use the warm-layer parameterization
109      !!
110      !! INPUT/OUTPUT:
111      !! -------------
112      !!    *  T_s  : always "bulk SST" as input                              [K]
113      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
114      !!              -> skin temperature as output if CSWL used              [K]
115      !!
116      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
117      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
118      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
119      !!
120      !! OPTIONAL INPUT:
121      !! ---------------
122      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
123      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
124      !!    *  slp    : sea-level pressure                                    [Pa]
125      !!
126      !! OPTIONAL OUTPUT:
127      !! ----------------
128      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
129      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
130      !!    * pHz_wl  : thickness of warm-layer                               [m]
131      !!
132      !! OUTPUT :
133      !! --------
134      !!    *  Cd     : drag coefficient
135      !!    *  Ch     : sensible heat coefficient
136      !!    *  Ce     : evaporation coefficient
137      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
138      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
139      !!    *  Ubzu  : bulk wind speed at zu                                 [m/s]
140      !!
141      !!
142      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
143      !!----------------------------------------------------------------------------------
144      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
145      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
146      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
147      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
148      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
149      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
150      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
151      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
152      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
153      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
154      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
155      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
156      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
157      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
158      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
159      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
160      !
161      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
162      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
163      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
164      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
165      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
166      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
167      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
168      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
169      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
170      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
171      !
172      INTEGER :: nbit, jit
173      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
174      !
175      REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
176      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
177      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air
178      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
179      REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu
180      REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
181      !
182      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt
183      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst    ! to back up the initial bulk SST
184      !
185      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6'
186      !!----------------------------------------------------------------------------------
187      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl)
188
189      nbit = nb_iter0
190      IF( PRESENT(nb_iter) ) nbit = nb_iter
191
192      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
193      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
194
195      !! Initializations for cool skin and warm layer:
196      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
197         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
198
199      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
200         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
201
202      IF( l_use_cs .OR. l_use_wl ) THEN
203         ALLOCATE ( zsst(jpi,jpj) )
204         zsst = T_s ! backing up the bulk SST
205         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
206         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
207      ENDIF
208
209      !! First guess of temperature and humidity at height zu:
210      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
211      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
212
213      !! Pot. temp. difference (and we don't want it to be 0!)
214      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
215      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
216
217      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
218
219      Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
220
221      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
222      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
223      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10)
224
225      z0     = charn_coare3p6(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
226      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
227
228      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
229      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
230
231      Cd     = MAX( (vkarmn/ztmp0)**2 , Cx_min )    ! first guess of Cd
232
233      ztmp0 = vkarmn2/LOG(zt/z0t)/Cd
234
235      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN)
236
237      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
238      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
239      zeta_u = (1._wp - ztmp1) *   ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & !  BRN < 0
240         &  +       ztmp1      * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 )                 !  BRN > 0
241
242      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
243      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
244
245      u_star = MAX ( Ubzu*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on)
246      t_star = dt_zu*ztmp0
247      q_star = dq_zu*ztmp0
248
249      ! What needs to be done if zt /= zu:
250      IF( .NOT. l_zt_equal_zu ) THEN
251         !! First update of values at zu (or zt for wind)
252         zeta_t = zt*zeta_u/zu
253         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
254         ztmp1 = LOG(zt/zu) + ztmp0
255         t_zu = t_zt - t_star/vkarmn*ztmp1
256         q_zu = q_zt - q_star/vkarmn*ztmp1
257         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
258         !
259         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
260         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
261      ENDIF
262
263      !! ITERATION BLOCK
264      DO jit = 1, nbit
265
266         !!Inverse of Obukov length (1/L) :
267         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Obukhov length]
268         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! 1/L (prevents FPE from stupid values from masked region later on...)
269
270         ztmp1 = u_star*u_star   ! u*^2
271
272         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
273         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
274         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
275         Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
276         ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0.
277
278         !! Stability parameters:
279         zeta_u = zu*ztmp0
280         zeta_u = SIGN( MIN(ABS(zeta_u),zeta_abs_max), zeta_u )
281         IF( .NOT. l_zt_equal_zu ) THEN
282            zeta_t = zt*ztmp0
283            zeta_t = SIGN( MIN(ABS(zeta_t),zeta_abs_max), zeta_t )
284         ENDIF
285
286         !! Adjustment the wind at 10m (not needed in the current algo form):
287         !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
288
289         !! Roughness lengthes z0, z0t (z0q = z0t) :
290         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
291         z0    = charn_coare3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
292         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
293
294         ztmp1 = ( znu_a / (z0*u_star) )**0.72_wp     ! COARE3.6-specific! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
295         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1 )   ! COARE3.6-specific!
296         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
297
298         !! Turbulent scales at zu :
299         ztmp0   = psi_h_coare(zeta_u)
300         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LB: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
301
302         t_star = dt_zu*ztmp1
303         q_star = dq_zu*ztmp1
304         u_star = MAX( Ubzu*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u)) , 1.E-9 )  !  (MAX => prevents FPE from stupid values from masked region later on)
305
306         IF( .NOT. l_zt_equal_zu ) THEN
307            !! Re-updating temperature and humidity at zu if zt /= zu :
308            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
309            t_zu = t_zt - t_star/vkarmn*ztmp1
310            q_zu = q_zt - q_star/vkarmn*ztmp1
311         ENDIF
312
313         IF( l_use_cs ) THEN
314            !! Cool-skin contribution
315
316            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
317               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
318
319            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
320
321            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
322            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
323            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
324         ENDIF
325
326         IF( l_use_wl ) THEN
327            !! Warm-layer contribution
328            CALL UPDATE_QNSOL_TAU( zu, T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_zu, Ubzu, slp, rad_lw, &
329               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
330            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
331            CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
332
333            !! Updating T_s and q_s !!!
334            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1)
335            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
336            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
337         ENDIF
338
339         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
340            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
341            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
342         ENDIF
343
344      END DO !DO jit = 1, nbit
345
346      ! compute transfer coefficients at zu :
347      ztmp0 = u_star/Ubzu
348      Cd   = MAX( ztmp0*ztmp0        , Cx_min )
349      Ch   = MAX( ztmp0*t_star/dt_zu , Cx_min )
350      Ce   = MAX( ztmp0*q_star/dq_zu , Cx_min )
351
352      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
353
354      IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min )
355      IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
356      IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
357
358      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
359      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
360      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
361
362      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
363
364   END SUBROUTINE turb_coare3p6
365
366
367   FUNCTION charn_coare3p6( pwnd )
368      !!-------------------------------------------------------------------
369      !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m
370      !!  "wind speed dependent formulation"
371      !!  (Eq. 13 in Edson et al., 2013)
372      !!
373      !! Author: L. Brodeau, July 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/)
374      !!-------------------------------------------------------------------
375      REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6
376      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! neutral wind speed at 10m
377      !
378      REAL(wp), PARAMETER :: charn0_max = 0.028  !: value above which the Charnock parameter levels off for winds > 18 m/s
379      !!-------------------------------------------------------------------
380      charn_coare3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp )
381      !!
382   END FUNCTION charn_coare3p6
383
384   FUNCTION charn_coare3p6_wave( pus, pwsh, pwps )
385      !!-------------------------------------------------------------------
386      !! Computes the Charnock parameter as a function of wave information and u*
387      !!
388      !!  (COARE 3.6, Fairall et al., 2018)
389      !!
390      !! Author: L. Brodeau, October 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/)
391      !!-------------------------------------------------------------------
392      REAL(wp), DIMENSION(jpi,jpj) :: charn_coare3p6_wave
393      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus   ! friction velocity             [m/s]
394      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwsh  ! significant wave height       [m]
395      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwps  ! phase speed of dominant waves [m/s]
396      !!-------------------------------------------------------------------
397      charn_coare3p6_wave = ( pwsh*0.2_wp*(pus/pwps)**2.2_wp ) * grav/(pus*pus)
398      !!
399   END FUNCTION charn_coare3p6_wave
400
401
402   FUNCTION psi_m_coare( pzeta )
403      !!----------------------------------------------------------------------------------
404      !! ** Purpose: compute the universal profile stability function for momentum
405      !!             COARE 3.0, Fairall et al. 2003
406      !!             pzeta : stability paramenter, z/L where z is altitude
407      !!                     measurement and L is M-O length
408      !!       Stability function for wind speed and scalars matching Kansas and free
409      !!       convection forms with weighting f convective form, follows Fairall et
410      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
411      !!       form from Beljaars and Holtslag (1991)
412      !!
413      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
414      !!----------------------------------------------------------------------------------
415      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
416      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
417      !
418      INTEGER  ::   ji, jj    ! dummy loop indices
419      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
420      !!----------------------------------------------------------------------------------
421      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
422            !
423            zta = pzeta(ji,jj)
424            !
425            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
426            !
427            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
428               & - 2.*ATAN(zphi_m) + 0.5*rpi
429            !
430            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
431            !
432            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
433               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
434            !
435            zf = zta*zta
436            zf = zf/(1. + zf)
437            zc = MIN(50._wp, 0.35_wp*zta)
438            zstab = 0.5 + SIGN(0.5_wp, zta)
439            !
440            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
441               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
442               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )  !     "
443      END_2D
444   END FUNCTION psi_m_coare
445
446
447   FUNCTION psi_h_coare( pzeta )
448      !!---------------------------------------------------------------------
449      !! Universal profile stability function for temperature and humidity
450      !! COARE 3.0, Fairall et al. 2003
451      !!
452      !! pzeta : stability paramenter, z/L where z is altitude measurement
453      !!         and L is M-O length
454      !!
455      !! Stability function for wind speed and scalars matching Kansas and free
456      !! convection forms with weighting f convective form, follows Fairall et
457      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
458      !! form from Beljaars and Holtslag (1991)
459      !!
460      !! Author: L. Brodeau, June 2016 / AeroBulk
461      !!         (https://github.com/brodeau/aerobulk/)
462      !!----------------------------------------------------------------
463      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
464      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
465      !
466      INTEGER  ::   ji, jj     ! dummy loop indices
467      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
468      !!----------------------------------------------------------------
469      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
470            !
471            zta = pzeta(ji,jj)
472            !
473            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
474            !
475            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
476            !
477            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
478            !
479            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
480               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
481            !
482            zf = zta*zta
483            zf = zf/(1. + zf)
484            zc = MIN(50._wp,0.35_wp*zta)
485            zstab = 0.5 + SIGN(0.5_wp, zta)
486            !
487            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
488               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
489               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
490      END_2D
491   END FUNCTION psi_h_coare
492
493   !!======================================================================
494END MODULE sbcblk_algo_coare3p6
Note: See TracBrowser for help on using the repository browser.