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/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC – NEMO

source: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90 @ 11845

Last change on this file since 11845 was 11845, checked in by laurent, 4 years ago

Improving syntax consistency

File size: 26.0 KB
RevLine 
[11284]1MODULE sbcblk_algo_coare3p0
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p0  ***
[11615]4   !!
5   !!                After Fairall et al, 2003
[11284]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
[11631]40   USE sbcblk_skin_coare ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
[11284]41
42   IMPLICIT NONE
43   PRIVATE
44
[11615]45   PUBLIC :: COARE3P0_INIT, TURB_COARE3P0
[11284]46
[11785]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
[11284]50
[11772]51   INTEGER , PARAMETER ::   nb_itt = 10             ! number of itterations
[11284]52
53   !!----------------------------------------------------------------------
54CONTAINS
55
[11615]56
57   SUBROUTINE coare3p0_init(l_use_cs, l_use_wl)
58      !!---------------------------------------------------------------------
59      !!                  ***  FUNCTION 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
[11816]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( ' COARE3P0_INIT => allocation of Tau_ac, Qnt_ac, dT_wl & Hz_wl failed!' )
[11615]74         Tau_ac(:,:) = 0._wp
[11772]75         Qnt_ac(:,:) = 0._wp
[11785]76         dT_wl(:,:)  = 0._wp
[11772]77         Hz_wl(:,:)  = Hwl_max
[11615]78      END IF
79      IF ( l_use_cs ) THEN
80         ierr = 0
[11772]81         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr )
[11816]82         IF( ierr > 0 ) CALL ctl_stop( ' COARE3P0_INIT => allocation of dT_cs failed!' )
[11772]83         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction
[11615]84      END IF
85   END SUBROUTINE coare3p0_init
86
87
88
[11672]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)
[11772]93      &                      pdT_wl, pHz_wl )                                                 ! optionals for warm-layer only
[11284]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
[11329]100      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
[11284]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      !! -------
[11615]109      !!    *  kt   : current time step (starts at 1)
[11284]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]
[11615]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
[11284]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]
[11615]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)
[11284]127      !!
[11615]128      !! OPTIONAL INPUT:
[11284]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]
[11785]133      !!
134      !! OPTIONAL OUTPUT:
135      !! ----------------
[11615]136      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
137      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
[11772]138      !!    * pHz_wl  : thickness of warm-layer                               [m]
[11284]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      !!----------------------------------------------------------------------------------
[11615]152      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
[11284]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]
[11615]160      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
161      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
[11284]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]
[11615]173      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
174      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
[11772]175      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pHz_wl   !             [m]
[11615]176      !
[11284]177      INTEGER :: j_itt
178      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
179      !
[11845]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
[11284]186      !
[11615]187      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
[11845]188         &              zeta_t,   &  ! stability parameter at height zt
[11615]189         &                zsst,   &  ! to back up the initial bulk SST
190         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
191         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K]
192         &                zHwl       ! depth of warm-layer [m]
[11772]193      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0'
[11284]194      !!----------------------------------------------------------------------------------
[11772]195      IF ( kt == nit000 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl)
[11284]196
197      l_zt_equal_zu = .FALSE.
198      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
199      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
200
[11615]201      !! Initializations for cool skin and warm layer:
[11816]202      IF ( l_use_cs .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 cool-skin param!' )
[11817]204
205      IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) &
206         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' )
207
[11672]208      IF ( l_use_cs .OR. l_use_wl ) THEN
209         ALLOCATE ( zsst(jpi,jpj) )
210         zsst = T_s ! backing up the bulk SST
[11772]211         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction
[11785]212         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s
[11672]213      END IF
[11615]214
215
[11284]216      !! First guess of temperature and humidity at height zu:
217      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
218      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
219
220      !! Pot. temp. difference (and we don't want it to be 0!)
221      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
222      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
223
224      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
225
226      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
227
228      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
229      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
230      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
231
232      z0     = alfa_charn_3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
[11785]233      z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
234
[11284]235      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
[11785]236      z0t    = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
[11284]237
[11615]238      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
[11284]239
240      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
241
242      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
243
244      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
245      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
246      ztmp0 = ztmp0*ztmp2
247      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
248         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
249      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
250
251      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
252      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
253
[11785]254      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)
[11284]255      t_star = dt_zu*ztmp0
256      q_star = dq_zu*ztmp0
257
258      ! What needs to be done if zt /= zu:
259      IF( .NOT. l_zt_equal_zu ) THEN
260         !! First update of values at zu (or zt for wind)
261         zeta_t = zt*zeta_u/zu
262         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
263         ztmp1 = LOG(zt/zu) + ztmp0
264         t_zu = t_zt - t_star/vkarmn*ztmp1
265         q_zu = q_zt - q_star/vkarmn*ztmp1
266         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
267         !
268         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
269         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
270      END IF
271
272      !! ITERATION BLOCK
273      DO j_itt = 1, nb_itt
274
275         !!Inverse of Monin-Obukov length (1/L) :
276         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
[11785]277         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...)
[11284]278
279         ztmp1 = u_star*u_star   ! u*^2
280
[11615]281         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
282         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
283         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
284         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
[11284]285         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
286
287         !! Stability parameters:
288         zeta_u = zu*ztmp0
289         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
290         IF( .NOT. l_zt_equal_zu ) THEN
291            zeta_t = zt*ztmp0
292            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
293         END IF
294
295         !! Adjustment the wind at 10m (not needed in the current algo form):
296         !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u))
[11615]297
[11284]298         !! Roughness lengthes z0, z0t (z0q = z0t) :
299         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
[11845]300         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
[11785]301         z0     = MIN( MAX(ABS(z0), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
302
[11631]303         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
[11615]304         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!!
[11785]305         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on)
[11284]306
307         !! Turbulent scales at zu :
308         ztmp0   = psi_h_coare(zeta_u)
[11785]309         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LB: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
[11284]310
311         t_star = dt_zu*ztmp1
312         q_star = dq_zu*ztmp1
[11785]313         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)
[11284]314
315         IF( .NOT. l_zt_equal_zu ) THEN
[11615]316            !! Re-updating temperature and humidity at zu if zt /= zu :
317            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
[11284]318            t_zu = t_zt - t_star/vkarmn*ztmp1
319            q_zu = q_zt - q_star/vkarmn*ztmp1
320         END IF
321
[11615]322
323         IF( l_use_cs ) THEN
324            !! Cool-skin contribution
325
[11772]326            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, &
[11615]327               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
328
[11772]329            CALL CS_COARE( Qsw, ztmp1, u_star, zsst, ztmp2 )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
[11615]330
[11772]331            T_s(:,:) = zsst(:,:) + dT_cs(:,:)*tmask(:,:,1)
332            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + dT_wl(:,:)*tmask(:,:,1)
[11615]333            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
334
[11284]335         END IF
336
[11615]337         IF( l_use_wl ) THEN
338            !! Warm-layer contribution
[11772]339            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, &
[11615]340               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
[11626]341            !! In WL_COARE or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
[11772]342            CALL WL_COARE( Qsw, ztmp1, zeta_u, zsst, MOD(nb_itt,j_itt) )
343
[11615]344            !! Updating T_s and q_s !!!
[11772]345            T_s(:,:) = zsst(:,:) + dT_wl(:,:)*tmask(:,:,1)
346            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1)
[11615]347            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
348         END IF
349
350         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
[11284]351            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
352            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
353         END IF
354
355      END DO !DO j_itt = 1, nb_itt
356
357      ! compute transfer coefficients at zu :
358      ztmp0 = u_star/U_blk
359      Cd   = ztmp0*ztmp0
360      Ch   = ztmp0*t_star/dt_zu
361      Ce   = ztmp0*q_star/dq_zu
[11329]362
[11284]363      ztmp1 = zu + z0
364      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
365      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
366      Cen = Chn
[11329]367
[11284]368      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
[11329]369
[11772]370      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs
371      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl
372      IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl
[11615]373
374      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
375
[11284]376   END SUBROUTINE turb_coare3p0
377
378
379   FUNCTION alfa_charn_3p0( pwnd )
380      !!-------------------------------------------------------------------
381      !! Compute the Charnock parameter as a function of the wind speed
382      !!
383      !! (Fairall et al., 2003 p.577-578)
384      !!
385      !! Wind below 10 m/s :  alfa = 0.011
386      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
387      !! Wind greater than 18 m/s :  alfa = 0.018
388      !!
389      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
390      !!-------------------------------------------------------------------
391      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0
392      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
393      !
394      INTEGER  ::   ji, jj         ! dummy loop indices
395      REAL(wp) :: zw, zgt10, zgt18
396      !!-------------------------------------------------------------------
397      !
398      DO jj = 1, jpj
399         DO ji = 1, jpi
400            !
401            zw = pwnd(ji,jj)   ! wind speed
402            !
403            ! Charnock's constant, increases with the wind :
[11631]404            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1
[11284]405            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
406            !
407            alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
408               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
409               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
410            !
411         END DO
412      END DO
413      !
414   END FUNCTION alfa_charn_3p0
415
416   FUNCTION psi_m_coare( pzeta )
417      !!----------------------------------------------------------------------------------
418      !! ** Purpose: compute the universal profile stability function for momentum
419      !!             COARE 3.0, Fairall et al. 2003
420      !!             pzeta : stability paramenter, z/L where z is altitude
421      !!                     measurement and L is M-O length
422      !!       Stability function for wind speed and scalars matching Kansas and free
423      !!       convection forms with weighting f convective form, follows Fairall et
424      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
425      !!       form from Beljaars and Holtslag (1991)
426      !!
427      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
428      !!----------------------------------------------------------------------------------
429      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
430      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
431      !
432      INTEGER  ::   ji, jj    ! dummy loop indices
433      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
434      !!----------------------------------------------------------------------------------
435      !
436      DO jj = 1, jpj
437         DO ji = 1, jpi
438            !
439            zta = pzeta(ji,jj)
440            !
441            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
442            !
443            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
444               & - 2.*ATAN(zphi_m) + 0.5*rpi
445            !
446            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
447            !
448            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
449               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
450            !
451            zf = zta*zta
452            zf = zf/(1. + zf)
453            zc = MIN(50._wp, 0.35_wp*zta)
454            zstab = 0.5 + SIGN(0.5_wp, zta)
455            !
456            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
457               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
458               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
459            !
460         END DO
461      END DO
462      !
463   END FUNCTION psi_m_coare
464
465
466   FUNCTION psi_h_coare( pzeta )
467      !!---------------------------------------------------------------------
468      !! Universal profile stability function for temperature and humidity
469      !! COARE 3.0, Fairall et al. 2003
470      !!
471      !! pzeta : stability paramenter, z/L where z is altitude measurement
472      !!         and L is M-O length
473      !!
474      !! Stability function for wind speed and scalars matching Kansas and free
475      !! convection forms with weighting f convective form, follows Fairall et
476      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
477      !! form from Beljaars and Holtslag (1991)
478      !!
479      !! Author: L. Brodeau, June 2016 / AeroBulk
480      !!         (https://github.com/brodeau/aerobulk/)
481      !!----------------------------------------------------------------
482      !!
483      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
484      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
485      !
486      INTEGER  ::   ji, jj     ! dummy loop indices
487      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
488      !
489      DO jj = 1, jpj
490         DO ji = 1, jpi
491            !
492            zta = pzeta(ji,jj)
493            !
494            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
495            !
496            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
497            !
498            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
499            !
500            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
501               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
502            !
503            zf = zta*zta
504            zf = zf/(1. + zf)
505            zc = MIN(50._wp,0.35_wp*zta)
506            zstab = 0.5 + SIGN(0.5_wp, zta)
507            !
508            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
509               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
510               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
511            !
512         END DO
513      END DO
514      !
515   END FUNCTION psi_h_coare
516
517   !!======================================================================
518END MODULE sbcblk_algo_coare3p0
Note: See TracBrowser for help on using the repository browser.