source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/SBC/sbcblk_algo_coare3p0.F90

Last change on this file was 13197, checked in by gsamson, 4 months ago

merge with trunk@r13136 with a more recent svn version; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

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