source: NEMO/trunk/src/OCE/SBC/sbcblk_algo_coare3p0.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 8 months ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge —ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The —ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

File size: 25.6 KB
RevLine 
[12015]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
[12340]46   !! * Substitutions
47#  include "do_loop_substitute.h90"
[12015]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 = .FALSE.
197      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
198      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
199
200      !! Initializations for cool skin and warm layer:
201      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
202         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' )
203
204      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
205         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
206
207      IF( l_use_cs .OR. l_use_wl ) THEN
208         ALLOCATE ( zsst(jpi,jpj) )
209         zsst = T_s ! backing up the bulk SST
210         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
211         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
212      ENDIF
213
214
215      !! First guess of temperature and humidity at height zu:
216      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
217      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
218
219      !! Pot. temp. difference (and we don't want it to be 0!)
220      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
221      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
222
223      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
224
225      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
226
227      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
228      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
229      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
230
231      z0     = alfa_charn_3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
232      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
233
234      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
235      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
236
237      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
238
239      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
240
241      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
242
243      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
244      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
245      ztmp0 = ztmp0*ztmp2
246      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
247         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
248      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
249
250      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
251      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
252
253      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)
254      t_star = dt_zu*ztmp0
255      q_star = dq_zu*ztmp0
256
257      ! What needs to be done if zt /= zu:
258      IF( .NOT. l_zt_equal_zu ) THEN
259         !! First update of values at zu (or zt for wind)
260         zeta_t = zt*zeta_u/zu
261         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
262         ztmp1 = LOG(zt/zu) + ztmp0
263         t_zu = t_zt - t_star/vkarmn*ztmp1
264         q_zu = q_zt - q_star/vkarmn*ztmp1
265         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
266         !
267         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
268         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
269      ENDIF
270
271      !! ITERATION BLOCK
272      DO j_itt = 1, nb_itt
273
274         !!Inverse of Monin-Obukov length (1/L) :
275         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
276         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...)
277
278         ztmp1 = u_star*u_star   ! u*^2
279
280         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
281         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
282         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
283         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
284         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
285
286         !! Stability parameters:
287         zeta_u = zu*ztmp0
288         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
289         IF( .NOT. l_zt_equal_zu ) THEN
290            zeta_t = zt*ztmp0
291            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
292         ENDIF
293
294         !! Adjustment the wind at 10m (not needed in the current algo form):
295         !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))
296
297         !! Roughness lengthes z0, z0t (z0q = z0t) :
298         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
299         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
300         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
301
302         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
303         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 !!!
304         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
305
306         !! Turbulent scales at zu :
307         ztmp0   = psi_h_coare(zeta_u)
308         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LB: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
309
310         t_star = dt_zu*ztmp1
311         q_star = dq_zu*ztmp1
312         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)
313
314         IF( .NOT. l_zt_equal_zu ) THEN
315            !! Re-updating temperature and humidity at zu if zt /= zu :
316            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
317            t_zu = t_zt - t_star/vkarmn*ztmp1
318            q_zu = q_zt - q_star/vkarmn*ztmp1
319         ENDIF
320
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, U_blk, 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
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, U_blk, 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(nb_itt,j_itt) )
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 j_itt = 1, nb_itt
355
356      ! compute transfer coefficients at zu :
357      ztmp0 = u_star/U_blk
358      Cd   = ztmp0*ztmp0
359      Ch   = ztmp0*t_star/dt_zu
360      Ce   = ztmp0*q_star/dq_zu
361
362      ztmp1 = zu + z0
363      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
364      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
365      Cen = Chn
366
367      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
368
369      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
370      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
371      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
372
373      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
374
375   END SUBROUTINE turb_coare3p0
376
377
378   FUNCTION alfa_charn_3p0( pwnd )
379      !!-------------------------------------------------------------------
380      !! Compute the Charnock parameter as a function of the wind speed
381      !!
382      !! (Fairall et al., 2003 p.577-578)
383      !!
384      !! Wind below 10 m/s :  alfa = 0.011
385      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
386      !! Wind greater than 18 m/s :  alfa = 0.018
387      !!
388      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
389      !!-------------------------------------------------------------------
390      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0
391      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
392      !
393      INTEGER  ::   ji, jj         ! dummy loop indices
394      REAL(wp) :: zw, zgt10, zgt18
395      !!-------------------------------------------------------------------
396      !
[12340]397      DO_2D_11_11
398         !
399         zw = pwnd(ji,jj)   ! wind speed
400         !
401         ! Charnock's constant, increases with the wind :
402         zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1
403         zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
404         !
405         alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
406            &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
407            &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
408         !
409      END_2D
[12015]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      !
[12340]433      DO_2D_11_11
434         !
435         zta = pzeta(ji,jj)
436         !
437         zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
438         !
439         zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
440            & - 2.*ATAN(zphi_m) + 0.5*rpi
441         !
442         zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
443         !
444         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
445            &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
446         !
447         zf = zta*zta
448         zf = zf/(1. + zf)
449         zc = MIN(50._wp, 0.35_wp*zta)
450         zstab = 0.5 + SIGN(0.5_wp, zta)
451         !
452         psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
453            &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
454            &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
455         !
456      END_2D
[12015]457      !
458   END FUNCTION psi_m_coare
459
460
461   FUNCTION psi_h_coare( pzeta )
462      !!---------------------------------------------------------------------
463      !! Universal profile stability function for temperature and humidity
464      !! COARE 3.0, Fairall et al. 2003
465      !!
466      !! pzeta : stability paramenter, z/L where z is altitude measurement
467      !!         and L is M-O length
468      !!
469      !! Stability function for wind speed and scalars matching Kansas and free
470      !! convection forms with weighting f convective form, follows Fairall et
471      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
472      !! form from Beljaars and Holtslag (1991)
473      !!
474      !! Author: L. Brodeau, June 2016 / AeroBulk
475      !!         (https://github.com/brodeau/aerobulk/)
476      !!----------------------------------------------------------------
477      !!
478      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
479      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
480      !
481      INTEGER  ::   ji, jj     ! dummy loop indices
482      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
483      !
[12340]484      DO_2D_11_11
485         !
486         zta = pzeta(ji,jj)
487         !
488         zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
489         !
490         zpsi_k = 2.*LOG((1. + zphi_h)/2.)
491         !
492         zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
493         !
494         zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
495            &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
496         !
497         zf = zta*zta
498         zf = zf/(1. + zf)
499         zc = MIN(50._wp,0.35_wp*zta)
500         zstab = 0.5 + SIGN(0.5_wp, zta)
501         !
502         psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
503            &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
504            &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
505         !
506      END_2D
[12015]507      !
508   END FUNCTION psi_h_coare
509
510   !!======================================================================
511END MODULE sbcblk_algo_coare3p0
Note: See TracBrowser for help on using the repository browser.