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/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_coare3p6.F90 @ 13655

Last change on this file since 13655 was 13655, checked in by laurent, 4 years ago

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