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

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

Commit all my dev of 2020!

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