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

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

Last change on this file since 15520 was 15520, checked in by sparonuz, 7 months ago

fix name of variables mispelled in subprogram declaration. Removed file added for error.

File size: 26.1 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   !!            4.2  ! 2020-12  (L. Brodeau) Introduction of various air-ice bulk parameterizations + improvements
19   !!----------------------------------------------------------------------
20
21   !!----------------------------------------------------------------------
22   !!   turb_coare3p0  : 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 oce             ! ocean dynamics and tracers
27   USE dom_oce         ! ocean space and time domain
28   USE phycst          ! physical constants
29   USE iom             ! I/O manager library
30   USE lib_mpp         ! distribued memory computing library
31   USE in_out_manager  ! I/O manager
32   USE prtctl          ! Print control
33   USE sbcwave, ONLY   :  cdn_wave ! wave module
34#if defined key_si3 || defined key_cice
35   USE sbc_ice         ! Surface boundary condition: ice fields
36#endif
37   USE lib_fortran     ! to use key_nosignedzero
38
39   USE sbc_oce         ! Surface boundary condition: ocean fields
40   USE sbc_phy         ! Catalog of functions for physical/meteorological parameters in the marine boundary layer
41   USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
42
43   IMPLICIT NONE
44   PRIVATE
45
46   PUBLIC :: sbcblk_algo_coare3p0_init, TURB_COARE3P0
47   !! * Substitutions
48#  include "do_loop_substitute.h90"
49
50   !! COARE own values for given constants:
51   REAL(wp), PARAMETER :: zi0   = 600._wp     ! scale height of the atmospheric boundary layer...
52   REAL(wp), PARAMETER :: Beta0 =  1.25_wp    ! gustiness parameter
53   REAL(wp), PARAMETER :: zeta_abs_max = 50._wp
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, Ubzu,                               &
93      &                      nb_iter, CdN, ChN, CeN,                                     & ! optional output
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      !!    *  Ubzu  : 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) ::   Ubzu    ! bulk wind speed at zu                     [m/s]
170      !
171      INTEGER , INTENT(in   ), OPTIONAL                     :: nb_iter  ! number of iterations
172      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CdN
173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   ChN
174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   CeN
175      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
176      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
177      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
178      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
179      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
181      !
182      INTEGER :: nbit, jit
183      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
184      !
185      REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star
186      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu
187      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air
188      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t
189      REAL(wp), DIMENSION(jpi,jpj) :: zeta_u        ! stability parameter at height zu
190      REAL(wp), DIMENSION(jpi,jpj) :: ztmp0, ztmp1, ztmp2
191      !
192      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zeta_t  ! stability parameter at height zt
193      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zsst    ! to back up the initial bulk SST
194      !
195      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0'
196      !!----------------------------------------------------------------------------------
197      IF( kt == nit000 ) CALL sbcblk_algo_coare3p0_init(l_use_cs, l_use_wl)
198
199      nbit = nb_iter0
200      IF( PRESENT(nb_iter) ) nbit = nb_iter
201
202      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision
203      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
204
205      !! Initializations for cool skin and warm layer:
206      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
207         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
208
209      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
210         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
211
212      IF( l_use_cs .OR. l_use_wl ) THEN
213         ALLOCATE ( zsst(jpi,jpj) )
214         zsst = T_s ! backing up the bulk SST
215         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
216         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
217      ENDIF
218
219      !! First guess of temperature and humidity at height zu:
220      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
221      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
222
223      !! Pot. temp. difference (and we don't want it to be 0!)
224      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
225      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
226
227      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
228
229      Ubzu = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
230
231      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
232      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
233      u_star = 0.035_wp*Ubzu*ztmp1/ztmp0       ! (u* = 0.035*Un10)
234
235      z0     = charn_coare3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
236      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
237
238      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
239      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
240
241      Cd     = MAX( (vkarmn/ztmp0)**2 , Cx_min )    ! first guess of Cd
242
243      ztmp0 = vkarmn2/LOG(zt/z0t)/Cd
244
245      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, Ubzu ) ! Bulk Richardson Number (BRN)
246
247      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
248      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
249      zeta_u = (1._wp - ztmp1) *   ztmp0*ztmp2 / (1._wp - ztmp2*zi0*0.004_wp*Beta0**3/zu) & !  BRN < 0
250         &  +       ztmp1      * ( ztmp0*ztmp2 + 27._wp/9._wp*ztmp2*ztmp2 )                 !  BRN > 0
251
252      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
253      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
254
255      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)
256      t_star = dt_zu*ztmp0
257      q_star = dq_zu*ztmp0
258
259      ! What needs to be done if zt /= zu:
260      IF( .NOT. l_zt_equal_zu ) THEN
261         !! First update of values at zu (or zt for wind)
262         zeta_t = zt*zeta_u/zu
263         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
264         ztmp1 = LOG(zt/zu) + ztmp0
265         t_zu = t_zt - t_star/vkarmn*ztmp1
266         q_zu = q_zt - q_star/vkarmn*ztmp1
267         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
268         !
269         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
270         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
271      ENDIF
272
273      !! ITERATION BLOCK
274      DO jit = 1, nbit
275
276         !!Inverse of Obukov length (1/L) :
277         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Obukhov length]
278         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! 1/L (prevents FPE from stupid values from masked region later on...)
279
280         ztmp1 = u_star*u_star   ! u*^2
281
282         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
283         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
284         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
285         Ubzu = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
286         ! => 0.2 prevents Ubzu to be 0 in stable case when U_zu=0.
287
288         !! Stability parameters:
289         zeta_u = zu*ztmp0
290         zeta_u = SIGN( MIN(ABS(zeta_u),zeta_abs_max), zeta_u )
291         IF( .NOT. l_zt_equal_zu ) THEN
292            zeta_t = zt*ztmp0
293            zeta_t = SIGN( MIN(ABS(zeta_t),zeta_abs_max), zeta_t )
294         ENDIF
295
296         !! Adjustment the wind at 10m (not needed in the current algo form):
297         !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))
298
299         !! Roughness lengthes z0, z0t (z0q = z0t) :
300         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
301         z0    = charn_coare3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
302         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
303
304         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
305         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 !!!
306         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
307
308         !! Turbulent scales at zu :
309         ztmp0   = psi_h_coare(zeta_u)
310         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LB: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
311
312         t_star = dt_zu*ztmp1
313         q_star = dq_zu*ztmp1
314         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)
315
316         IF( .NOT. l_zt_equal_zu ) THEN
317            !! Re-updating temperature and humidity at zu if zt /= zu :
318            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
319            t_zu = t_zt - t_star/vkarmn*ztmp1
320            q_zu = q_zt - q_star/vkarmn*ztmp1
321         ENDIF
322
323         IF( l_use_cs ) THEN
324            !! Cool-skin contribution
325
326            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, &
327               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
328
329            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
330
331            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
332            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
333            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
334         ENDIF
335
336         IF( l_use_wl ) THEN
337            !! Warm-layer contribution
338            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, &
339               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
340            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
341            CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nbit,jit) )
342
343            !! Updating T_s and q_s !!!
344            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1)
345            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
346            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
347         ENDIF
348
349         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
350            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
351            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
352         ENDIF
353
354      END DO !DO jit = 1, nbit
355
356      ! compute transfer coefficients at zu :
357      ztmp0 = u_star/Ubzu
358      Cd   = MAX( ztmp0*ztmp0        , Cx_min )
359      Ch   = MAX( ztmp0*t_star/dt_zu , Cx_min )
360      Ce   = MAX( ztmp0*q_star/dq_zu , Cx_min )
361
362      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
363
364      IF(PRESENT(Cdn)) Cdn = MAX( vkarmn2 / (LOG(zu/z0 )*LOG(zu/z0 )) , Cx_min )
365      IF(PRESENT(Chn)) Chn = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
366      IF(PRESENT(Cen)) Cen = MAX( vkarmn2 / (LOG(zu/z0t)*LOG(zu/z0t)) , Cx_min )
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 charn_coare3p0( 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) :: charn_coare3p0
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      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
396            !
397            zw = pwnd(ji,jj)   ! wind speed
398            !
399            ! Charnock's constant, increases with the wind :
400            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1
401            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18._wp)) ! If zw<18. --> 0, else --> 1
402            !
403            charn_coare3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
404               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
405               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
406            !
407      END_2D
408   END FUNCTION charn_coare3p0
409
410   FUNCTION psi_m_coare( pzeta )
411      !!----------------------------------------------------------------------------------
412      !! ** Purpose: compute the universal profile stability function for momentum
413      !!             COARE 3.0, Fairall et al. 2003
414      !!             pzeta : stability paramenter, z/L where z is altitude
415      !!                     measurement and L is M-O length
416      !!       Stability function for wind speed and scalars matching Kansas and free
417      !!       convection forms with weighting f convective form, follows Fairall et
418      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
419      !!       form from Beljaars and Holtslag (1991)
420      !!
421      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
422      !!----------------------------------------------------------------------------------
423      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
424      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
425      !
426      INTEGER  ::   ji, jj    ! dummy loop indices
427      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
428      !!----------------------------------------------------------------------------------
429      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
430            !
431            zta = pzeta(ji,jj)
432            !
433            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
434            !
435            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
436               & - 2.*ATAN(zphi_m) + 0.5*rpi
437            !
438            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
439            !
440            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
441               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
442            !
443            zf = zta*zta
444            zf = zf/(1. + zf)
445            zc = MIN(50._wp, 0.35_wp*zta)
446            zstab = 0.5 + SIGN(0.5_wp, zta)
447            !
448            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
449               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
450               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )  !     "
451      END_2D
452   END FUNCTION psi_m_coare
453
454
455   FUNCTION psi_h_coare( pzeta )
456      !!---------------------------------------------------------------------
457      !! Universal profile stability function for temperature and humidity
458      !! COARE 3.0, Fairall et al. 2003
459      !!
460      !! pzeta : stability paramenter, z/L where z is altitude measurement
461      !!         and L is M-O length
462      !!
463      !! Stability function for wind speed and scalars matching Kansas and free
464      !! convection forms with weighting f convective form, follows Fairall et
465      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
466      !! form from Beljaars and Holtslag (1991)
467      !!
468      !! Author: L. Brodeau, June 2016 / AeroBulk
469      !!         (https://github.com/brodeau/aerobulk/)
470      !!----------------------------------------------------------------
471      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
472      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
473      !
474      INTEGER  ::   ji, jj     ! dummy loop indices
475      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
476      !!----------------------------------------------------------------
477      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
478            !
479            zta = pzeta(ji,jj)
480            !
481            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
482            !
483            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
484            !
485            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
486            !
487            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
488               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
489            !
490            zf = zta*zta
491            zf = zf/(1. + zf)
492            zc = MIN(50._wp,0.35_wp*zta)
493            zstab = 0.5 + SIGN(0.5_wp, zta)
494            !
495            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
496               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
497               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
498      END_2D
499   END FUNCTION psi_h_coare
500
501   !!======================================================================
502END MODULE sbcblk_algo_coare3p0
Note: See TracBrowser for help on using the repository browser.