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

Last change on this file since 11615 was 11615, checked in by laurent, 16 months ago

LB: new cool-skin and warm-layer parameterizations for ECMWF and COARE3.6, COARE3.0 uses same CSWL param as for COARE3.6.

File size: 27.4 KB
Line 
1MODULE sbcblk_algo_coare3p0
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p0  ***
4   !!
5   !!                After Fairall et al, 2003
6   !! Computes:
7   !!   * bulk transfer coefficients C_D, C_E and C_H
8   !!   * air temp. and spec. hum. adjusted from zt (2m) to zu (10m) if needed
9   !!   * the effective bulk wind speed at 10m U_blk
10   !!   => all these are used in bulk formulas in sbcblk.F90
11   !!
12   !!       Routine turb_coare3p0 maintained and developed in AeroBulk
13   !!                     (https://github.com/brodeau/aerobulk)
14   !!
15   !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk)
16   !!----------------------------------------------------------------------
17   !! History :  4.0  !  2016-02  (L.Brodeau)   Original code
18   !!----------------------------------------------------------------------
19
20   !!----------------------------------------------------------------------
21   !!   turb_coare3p0  : computes the bulk turbulent transfer coefficients
22   !!                   adjusts t_air and q_air from zt to zu m
23   !!                   returns the effective bulk wind speed at 10m
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE phycst          ! physical constants
28   USE iom             ! I/O manager library
29   USE lib_mpp         ! distribued memory computing library
30   USE in_out_manager  ! I/O manager
31   USE prtctl          ! Print control
32   USE sbcwave, ONLY   :  cdn_wave ! wave module
33#if defined key_si3 || defined key_cice
34   USE sbc_ice         ! Surface boundary condition: ice fields
35#endif
36   USE lib_fortran     ! to use key_nosignedzero
37
38   USE sbc_oce         ! Surface boundary condition: ocean fields
39   USE sbcblk_phy      ! all thermodynamics functions, rho_air, q_sat, etc... !LB
40   USE sbcblk_skin_coare3p6  ! cool-skin/warm layer scheme (CSWL_ECMWF) !LB
41
42   IMPLICIT NONE
43   PRIVATE
44
45   PUBLIC :: 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 = 5             ! number of itterations
52
53   !!----------------------------------------------------------------------
54CONTAINS
55
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
72         PRINT *, ' *** coare3p0_init: WL => allocating Tau_ac, Qnt_ac, and H_wl :', jpi,jpj
73         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), H_wl(jpi,jpj), STAT=ierr )
74         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of Tau_ac and Qnt_ac failed!'
75         Tau_ac(:,:) = 0._wp
76         Qnt_ac(:,:)   = 0._wp
77         H_wl(:,:)    = H_wl_max
78         PRINT *, ' *** Tau_ac , Qnt_ac, and H_wl allocated!'
79      END IF
80      !!
81      IF ( l_use_cs ) THEN
82         ierr = 0
83         PRINT *, ' *** coare3p0_init: CS => allocating delta_vl :', jpi,jpj
84         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr )
85         !IF( ierr > 0 ) STOP ' COARE3P0_INIT => allocation of delta_vl and Qnt_ac failed!'
86         delta_vl(:,:) = 0.001_wp      ! First guess of zdelta [m]
87         PRINT *, ' *** delta_vl allocated!'
88      END IF
89   END SUBROUTINE coare3p0_init
90
91
92
93   SUBROUTINE turb_coare3p0( kt, zt, zu, T_s, t_zt, q_s, q_zt, U_zu, l_use_cs, l_use_wl,  &
94      &                      Cd, Ch, Ce, t_zu, q_zu, U_blk,                               &
95      &                      Cdn, Chn, Cen,                      &
96      &                      Qsw, rad_lw, slp, pdT_cs,                                    & ! optionals for cool-skin (and warm-layer)
97      &                      isecday_utc, plong, pdT_wl, Hwl )                             ! optionals for warm-layer only
98      !!----------------------------------------------------------------------
99      !!                      ***  ROUTINE  turb_coare3p0  ***
100      !!
101      !! ** Purpose :   Computes turbulent transfert coefficients of surface
102      !!                fluxes according to Fairall et al. (2003)
103      !!                If relevant (zt /= zu), adjust temperature and humidity from height zt to zu
104      !!                Returns the effective bulk wind speed at zu to be used in the bulk formulas
105      !!
106      !!                Applies the cool-skin warm-layer correction of the SST to T_s
107      !!                if the net shortwave flux at the surface (Qsw), the downwelling longwave
108      !!                radiative fluxes at the surface (rad_lw), and the sea-leve pressure (slp)
109      !!                are provided as (optional) arguments!
110      !!
111      !! INPUT :
112      !! -------
113      !!    *  kt   : current time step (starts at 1)
114      !!    *  zt   : height for temperature and spec. hum. of air            [m]
115      !!    *  zu   : height for wind speed (usually 10m)                     [m]
116      !!    *  t_zt : potential air temperature at zt                         [K]
117      !!    *  q_zt : specific humidity of air at zt                          [kg/kg]
118      !!    *  U_zu : scalar wind speed at zu                                 [m/s]
119      !!    * l_use_cs : use the cool-skin parameterization
120      !!    * l_use_wl : use the warm-layer parameterization
121      !!
122      !! INPUT/OUTPUT:
123      !! -------------
124      !!    *  T_s  : always "bulk SST" as input                              [K]
125      !!              -> unchanged "bulk SST" as output if CSWL not used      [K]
126      !!              -> skin temperature as output if CSWL used              [K]
127      !!
128      !!    *  q_s  : SSQ aka saturation specific humidity at temp. T_s       [kg/kg]
129      !!              -> doesn't need to be given a value if skin temp computed (in case l_use_cs=True or l_use_wl=True)
130      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_cs=False or l_use_wl=False)
131      !!
132      !! OPTIONAL INPUT:
133      !! ---------------
134      !!    *  Qsw    : net solar flux (after albedo) at the surface (>0)     [W/m^2]
135      !!    *  rad_lw : downwelling longwave radiation at the surface  (>0)   [W/m^2]
136      !!    *  slp    : sea-level pressure                                    [Pa]
137      !!    * pdT_cs  : SST increment "dT" for cool-skin correction           [K]
138      !!    * isecday_utc:
139      !!    *  plong  : longitude array                                       [deg.E]
140      !!    * pdT_wl  : SST increment "dT" for warm-layer correction          [K]
141      !!    * Hwl     : depth of warm layer                                   [m]
142      !!
143      !! OUTPUT :
144      !! --------
145      !!    *  Cd     : drag coefficient
146      !!    *  Ch     : sensible heat coefficient
147      !!    *  Ce     : evaporation coefficient
148      !!    *  t_zu   : pot. air temperature adjusted at wind height zu       [K]
149      !!    *  q_zu   : specific humidity of air        //                    [kg/kg]
150      !!    *  U_blk  : bulk wind speed at zu                                 [m/s]
151      !!
152      !!
153      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
154      !!----------------------------------------------------------------------------------
155      INTEGER,  INTENT(in   )                     ::   kt       ! current time step
156      REAL(wp), INTENT(in   )                     ::   zt       ! height for t_zt and q_zt                    [m]
157      REAL(wp), INTENT(in   )                     ::   zu       ! height for U_zu                             [m]
158      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   T_s      ! sea surface temperature                [Kelvin]
159      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   t_zt     ! potential air temperature              [Kelvin]
160      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj) ::   q_s      ! sea surface specific humidity           [kg/kg]
161      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   q_zt     ! specific air humidity at zt             [kg/kg]
162      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   U_zu     ! relative wind module at zu                [m/s]
163      LOGICAL , INTENT(in   )                     ::   l_use_cs ! use the cool-skin parameterization
164      LOGICAL , INTENT(in   )                     ::   l_use_wl ! use the warm-layer parameterization
165      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cd       ! transfer coefficient for momentum         (tau)
166      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ch       ! transfer coefficient for sensible heat (Q_sens)
167      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Ce       ! transfert coefficient for evaporation   (Q_lat)
168      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   t_zu     ! pot. air temp. adjusted at zu               [K]
169      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   q_zu     ! spec. humidity adjusted at zu           [kg/kg]
170      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   U_blk    ! bulk wind speed at zu                     [m/s]
171      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   Cdn, Chn, Cen ! neutral transfer coefficients
172      !
173      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   Qsw      !             [W/m^2]
174      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   rad_lw   !             [W/m^2]
175      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   slp      !             [Pa]
176      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_cs
177      !
178      INTEGER,  INTENT(in   ), OPTIONAL                     ::   isecday_utc ! current UTC time, counted in second since 00h of the current day
179      REAL(wp), INTENT(in   ), OPTIONAL, DIMENSION(jpi,jpj) ::   plong    !             [deg.E]
180      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   pdT_wl   !             [K]
181      REAL(wp), INTENT(  out), OPTIONAL, DIMENSION(jpi,jpj) ::   Hwl      !             [m]
182      !
183      INTEGER :: j_itt
184      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U
185      !
186      REAL(wp), DIMENSION(jpi,jpj) ::  &
187         &  u_star, t_star, q_star, &
188         &  dt_zu, dq_zu,    &
189         &  znu_a,           & !: Nu_air, Viscosity of air
190         &  z0, z0t
191      REAL(wp), DIMENSION(jpi,jpj) ::   zeta_u        ! stability parameter at height zu
192      REAL(wp), DIMENSION(jpi,jpj) ::   ztmp0, ztmp1, ztmp2
193      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zeta_t        ! stability parameter at height zt
194      !
195      REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
196         &                zsst,   &  ! to back up the initial bulk SST
197         &                pdTc,   &  ! SST increment "dT" for cool-skin correction           [K]
198         &                pdTw,   &  ! SST increment "dT" for warm layer correction          [K]
199         &                zHwl       ! depth of warm-layer [m]
200
201      !
202      LOGICAL :: lreturn_z0=.FALSE., lreturn_ustar=.FALSE., lreturn_L=.FALSE., lreturn_UN10=.FALSE.
203      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0.F90'
204      CHARACTER(len=128) :: cf_tmp
205      !!----------------------------------------------------------------------------------
206
207      IF ( kt == 1 ) CALL COARE3P0_INIT(l_use_cs, l_use_wl) ! allocation of accumulation arrays
208
209
210      l_zt_equal_zu = .FALSE.
211      IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision
212      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) )
213
214      !! Initializations for cool skin and warm layer:
215      IF ( l_use_cs ) THEN
216         IF( .NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp)) ) THEN
217            PRINT *, ' * PROBLEM ('//trim(crtnm)//'): you need to provide Qsw, rad_lw & slp to use cool-skin param!'
218            STOP
219         END IF
220         ALLOCATE ( pdTc(jpi,jpj) )
221         pdTc(:,:) = -0.25_wp  ! First guess of skin correction
222      END IF
223
224      IF ( l_use_wl ) THEN
225         IF(.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp) .AND. PRESENT(isecday_utc) .AND. PRESENT(plong))) THEN
226            PRINT *, ' * PROBLEM ('//TRIM(crtnm)//'): you need to provide Qsw, rad_lw, slp, isecday_utc & plong to use warm-layer param!'
227            STOP
228         END IF
229         ALLOCATE ( pdTw(jpi,jpj) )
230         IF (PRESENT(Hwl)) ALLOCATE ( zHwl(jpi,jpj) )
231      END IF
232
233      IF ( l_use_cs .OR. l_use_wl ) THEN
234         ALLOCATE ( zsst(jpi,jpj) )
235         zsst = T_s ! backing up the bulk SST
236         IF( l_use_cs ) T_s = T_s - 0.25   ! First guess of correction
237         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s !LOLO WL too!!!
238      END IF
239
240
241      !! First guess of temperature and humidity at height zu:
242      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions...
243      q_zu = MAX( q_zt , 1.e-6_wp )   !               "
244
245      !! Pot. temp. difference (and we don't want it to be 0!)
246      dt_zu = t_zu - T_s ;   dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
247      dq_zu = q_zu - q_s ;   dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
248
249      znu_a = visc_air(t_zu) ! Air viscosity (m^2/s) at zt given from temperature in (K)
250
251      U_blk = SQRT(U_zu*U_zu + 0.5_wp*0.5_wp) ! initial guess for wind gustiness contribution
252
253      ztmp0   = LOG(    zu*10000._wp) ! optimization: 10000. == 1/z0 (with z0 first guess == 0.0001)
254      ztmp1   = LOG(10._wp*10000._wp) !       "                    "               "
255      u_star = 0.035_wp*U_blk*ztmp1/ztmp0       ! (u* = 0.035*Un10)
256
257      z0     = alfa_charn_3p0(U_zu)*u_star*u_star/grav + 0.11_wp*znu_a/u_star
258      z0     = MIN(ABS(z0), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
259      z0t    = 1._wp / ( 0.1_wp*EXP(vkarmn/(0.00115/(vkarmn/ztmp1))) )
260      z0t    = MIN(ABS(z0t), 0.001_wp)  ! (prevent FPE from stupid values from masked region later on...) !#LOLO
261
262      Cd     = (vkarmn/ztmp0)**2    ! first guess of Cd
263
264      ztmp0 = vkarmn*vkarmn/LOG(zt/z0t)/Cd
265
266      ztmp2 = Ri_bulk( zu, T_s, t_zu, q_s, q_zu, U_blk ) ! Bulk Richardson Number (BRN)
267
268      !! First estimate of zeta_u, depending on the stability, ie sign of BRN (ztmp2):
269      ztmp1 = 0.5 + SIGN( 0.5_wp , ztmp2 )
270      ztmp0 = ztmp0*ztmp2
271      zeta_u = (1._wp-ztmp1) * (ztmp0/(1._wp+ztmp2/(-zu/(zi0*0.004_wp*Beta0**3)))) & !  BRN < 0
272         &  +     ztmp1   * (ztmp0*(1._wp + 27._wp/9._wp*ztmp2/ztmp0))               !  BRN > 0
273      !#LB: should make sure that the "ztmp0" of "27./9.*ztmp2/ztmp0" is "ztmp0*ztmp2" and not "ztmp0==vkarmn*vkarmn/LOG(zt/z0t)/Cd" !
274
275      !! First guess M-O stability dependent scaling params.(u*,t*,q*) to estimate z0 and z/L
276      ztmp0  = vkarmn/(LOG(zu/z0t) - psi_h_coare(zeta_u))
277
278      u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0)  - psi_m_coare(zeta_u))
279      t_star = dt_zu*ztmp0
280      q_star = dq_zu*ztmp0
281
282      ! What needs to be done if zt /= zu:
283      IF( .NOT. l_zt_equal_zu ) THEN
284         !! First update of values at zu (or zt for wind)
285         zeta_t = zt*zeta_u/zu
286         ztmp0 = psi_h_coare(zeta_u) - psi_h_coare(zeta_t)
287         ztmp1 = LOG(zt/zu) + ztmp0
288         t_zu = t_zt - t_star/vkarmn*ztmp1
289         q_zu = q_zt - q_star/vkarmn*ztmp1
290         q_zu = (0.5_wp + SIGN(0.5_wp,q_zu))*q_zu !Makes it impossible to have negative humidity :
291         !
292         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
293         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
294      END IF
295
296      !! ITERATION BLOCK
297      DO j_itt = 1, nb_itt
298
299         !!Inverse of Monin-Obukov length (1/L) :
300         ztmp0 = One_on_L(t_zu, q_zu, u_star, t_star, q_star)  ! 1/L == 1/[Monin-Obukhov length]
301         ztmp0 = SIGN( MIN(ABS(ztmp0),200._wp), ztmp0 ) ! (prevents FPE from stupid values from masked region later on...) !#LOLO
302
303         ztmp1 = u_star*u_star   ! u*^2
304
305         !! Update wind at zu with convection-related wind gustiness in unstable conditions (Fairall et al. 2003, Eq.8):
306         ztmp2 = Beta0*Beta0*ztmp1*(MAX(-zi0*ztmp0/vkarmn,0._wp))**(2._wp/3._wp) ! square of wind gustiness contribution, ztmp2 == Ug^2
307         !!   ! Only true when unstable (L<0) => when ztmp0 < 0 => explains "-" before zi0
308         U_blk = MAX(SQRT(U_zu*U_zu + ztmp2), 0.2_wp)        ! include gustiness in bulk wind speed
309         ! => 0.2 prevents U_blk to be 0 in stable case when U_zu=0.
310
311         !! Stability parameters:
312         zeta_u = zu*ztmp0
313         zeta_u = SIGN( MIN(ABS(zeta_u),50.0_wp), zeta_u )
314         IF( .NOT. l_zt_equal_zu ) THEN
315            zeta_t = zt*ztmp0
316            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t )
317         END IF
318
319         !! Adjustment the wind at 10m (not needed in the current algo form):
320         !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))
321
322         !! Roughness lengthes z0, z0t (z0q = z0t) :
323         ztmp2 = u_star/vkarmn*LOG(10./z0)                                 ! Neutral wind speed at 10m
324         z0    = alfa_charn_3p0(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6)
325         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp     ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
326         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 !!!
327
328         !! Turbulent scales at zu :
329         ztmp0   = psi_h_coare(zeta_u)
330         ztmp1   = vkarmn/(LOG(zu) - LOG(z0t) - ztmp0) ! #LOLO: in ztmp0, some use psi_h_coare(zeta_t) rather than psi_h_coare(zeta_t) ???
331
332         t_star = dt_zu*ztmp1
333         q_star = dq_zu*ztmp1
334         u_star = U_blk*vkarmn/(LOG(zu) - LOG(z0) - psi_m_coare(zeta_u))
335
336         IF( .NOT. l_zt_equal_zu ) THEN
337            !! Re-updating temperature and humidity at zu if zt /= zu :
338            ztmp1 = LOG(zt/zu) + ztmp0 - psi_h_coare(zeta_t)
339            t_zu = t_zt - t_star/vkarmn*ztmp1
340            q_zu = q_zt - q_star/vkarmn*ztmp1
341         END IF
342
343
344         IF( l_use_cs ) THEN
345            !! Cool-skin contribution
346
347            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
348               &                   ztmp1, zeta_u,  Qlat=ztmp2)  ! Qnsol -> ztmp1 / Tau -> zeta_u
349
350            CALL CS_COARE3P6( Qsw, ztmp1, u_star, zsst, ztmp2,  pdTc )  ! ! Qnsol -> ztmp1 / Qlat -> ztmp2
351
352            T_s(:,:) = zsst(:,:) + pdTc(:,:)*tmask(:,:,1)
353            IF( l_use_wl ) T_s(:,:) = T_s(:,:) + pdTw(:,:)*tmask(:,:,1)
354            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
355
356         END IF
357
358         IF( l_use_wl ) THEN
359            !! Warm-layer contribution
360            CALL UPDATE_QNSOL_TAU( T_s, q_s, t_zu, q_zu, u_star, t_star, q_star, U_blk, slp, rad_lw, &
361               &                   ztmp1, zeta_u)  ! Qnsol -> ztmp1 / Tau -> zeta_u
362            !! In WL_COARE3P6 or , Tau_ac and Qnt_ac must be updated at the final itteration step => add a flag to do this!
363            IF (PRESENT(Hwl)) THEN
364               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
365            ELSE
366               CALL WL_COARE3P6( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw )
367            END IF
368            !! Updating T_s and q_s !!!
369            T_s(:,:) = zsst(:,:) + pdTw(:,:)*tmask(:,:,1)
370            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + pdTc(:,:)*tmask(:,:,1)
371            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:))
372            !LOLO:
373            !IF ( j_itt == nb_itt) THEN
374            !   WRITE(cf_tmp,'("Qnt_ac_k",i5.5,"_p",i4.4,".nc")') kt, narea
375            !   CALL DUMP_FIELD(REAL( Qnt_ac*tmask(:,:,1) , 4), TRIM(cf_tmp), 'Qnt_ac')
376            !   WRITE(cf_tmp,  '("pdTw_k",i5.5,"_p",i4.4,".nc")') kt, narea
377            !   CALL DUMP_FIELD(REAL( pdTw*tmask(:,:,1) , 4), TRIM(cf_tmp), 'pdTw')
378            !END IF
379            !LOLO.
380         END IF
381
382
383         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN
384            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu )
385            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu )
386         END IF
387
388      END DO !DO j_itt = 1, nb_itt
389
390      ! compute transfer coefficients at zu :
391      ztmp0 = u_star/U_blk
392      Cd   = ztmp0*ztmp0
393      Ch   = ztmp0*t_star/dt_zu
394      Ce   = ztmp0*q_star/dq_zu
395
396      ztmp1 = zu + z0
397      Cdn = vkarmn*vkarmn / (log(ztmp1/z0 )*log(ztmp1/z0 ))
398      Chn = vkarmn*vkarmn / (log(ztmp1/z0t)*log(ztmp1/z0t))
399      Cen = Chn
400
401      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t )
402
403      IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = pdTc*tmask(:,:,1)
404      IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = pdTw*tmask(:,:,1)
405
406      IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst )
407      IF (          l_use_cs      ) DEALLOCATE ( pdTc )
408      IF (          l_use_wl      ) THEN
409         DEALLOCATE ( pdTw )
410         IF (PRESENT(Hwl)) DEALLOCATE ( zHwl )
411      END IF
412
413   END SUBROUTINE turb_coare3p0
414
415
416   FUNCTION alfa_charn_3p0( pwnd )
417      !!-------------------------------------------------------------------
418      !! Compute the Charnock parameter as a function of the wind speed
419      !!
420      !! (Fairall et al., 2003 p.577-578)
421      !!
422      !! Wind below 10 m/s :  alfa = 0.011
423      !! Wind between 10 and 18 m/s : linear increase from 0.011 to 0.018
424      !! Wind greater than 18 m/s :  alfa = 0.018
425      !!
426      !! Author: L. Brodeau, June 2016 / AeroBulk  (https://github.com/brodeau/aerobulk/)
427      !!-------------------------------------------------------------------
428      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p0
429      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! wind speed
430      !
431      INTEGER  ::   ji, jj         ! dummy loop indices
432      REAL(wp) :: zw, zgt10, zgt18
433      !!-------------------------------------------------------------------
434      !
435      DO jj = 1, jpj
436         DO ji = 1, jpi
437            !
438            zw = pwnd(ji,jj)   ! wind speed
439            !
440            ! Charnock's constant, increases with the wind :
441            zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10)) ! If zw<10. --> 0, else --> 1
442            zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1
443            !
444            alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s
445               &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) &
446               &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999)
447            !
448         END DO
449      END DO
450      !
451   END FUNCTION alfa_charn_3p0
452
453   FUNCTION psi_m_coare( pzeta )
454      !!----------------------------------------------------------------------------------
455      !! ** Purpose: compute the universal profile stability function for momentum
456      !!             COARE 3.0, Fairall et al. 2003
457      !!             pzeta : stability paramenter, z/L where z is altitude
458      !!                     measurement and L is M-O length
459      !!       Stability function for wind speed and scalars matching Kansas and free
460      !!       convection forms with weighting f convective form, follows Fairall et
461      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
462      !!       form from Beljaars and Holtslag (1991)
463      !!
464      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
465      !!----------------------------------------------------------------------------------
466      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
467      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
468      !
469      INTEGER  ::   ji, jj    ! dummy loop indices
470      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
471      !!----------------------------------------------------------------------------------
472      !
473      DO jj = 1, jpj
474         DO ji = 1, jpi
475            !
476            zta = pzeta(ji,jj)
477            !
478            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
479            !
480            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
481               & - 2.*ATAN(zphi_m) + 0.5*rpi
482            !
483            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
484            !
485            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
486               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
487            !
488            zf = zta*zta
489            zf = zf/(1. + zf)
490            zc = MIN(50._wp, 0.35_wp*zta)
491            zstab = 0.5 + SIGN(0.5_wp, zta)
492            !
493            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
494               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
495               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
496            !
497         END DO
498      END DO
499      !
500   END FUNCTION psi_m_coare
501
502
503   FUNCTION psi_h_coare( pzeta )
504      !!---------------------------------------------------------------------
505      !! Universal profile stability function for temperature and humidity
506      !! COARE 3.0, Fairall et al. 2003
507      !!
508      !! pzeta : stability paramenter, z/L where z is altitude measurement
509      !!         and L is M-O length
510      !!
511      !! Stability function for wind speed and scalars matching Kansas and free
512      !! convection forms with weighting f convective form, follows Fairall et
513      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
514      !! form from Beljaars and Holtslag (1991)
515      !!
516      !! Author: L. Brodeau, June 2016 / AeroBulk
517      !!         (https://github.com/brodeau/aerobulk/)
518      !!----------------------------------------------------------------
519      !!
520      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
521      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
522      !
523      INTEGER  ::   ji, jj     ! dummy loop indices
524      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
525      !
526      DO jj = 1, jpj
527         DO ji = 1, jpi
528            !
529            zta = pzeta(ji,jj)
530            !
531            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
532            !
533            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
534            !
535            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
536            !
537            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
538               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
539            !
540            zf = zta*zta
541            zf = zf/(1. + zf)
542            zc = MIN(50._wp,0.35_wp*zta)
543            zstab = 0.5 + SIGN(0.5_wp, zta)
544            !
545            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
546               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
547               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
548            !
549         END DO
550      END DO
551      !
552   END FUNCTION psi_h_coare
553
554   !!======================================================================
555END MODULE sbcblk_algo_coare3p0
Note: See TracBrowser for help on using the repository browser.