source: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/SBC/sbcblk_algo_coare3p0.F90

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

dev_ASINTER-01-05_merged: merge dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk branch @rev11988 with dev_r11265_ASINTER-01_Guillaume_ABL1D branch @rev11937 (tickets #2159 and #2131); ORCA2_ICE(_ABL) reproductibility OK

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