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_coare3p6.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_coare3p6.F90 @ 11631

Last change on this file since 11631 was 11631, checked in by laurent, 5 years ago

LB: syntax improved...

File size: 26.7 KB
Line 
1MODULE sbcblk_algo_coare3p6
2   !!======================================================================
3   !!                   ***  MODULE  sbcblk_algo_coare3p6  ***
4   !!
5   !!                After Fairall et al 2018 & Edson et al 2013
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_coare3p6 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_coare3p6  : 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 :: COARE3P6_INIT, TURB_COARE3P6
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.2_wp   ! gustiness parameter
50
51   INTEGER , PARAMETER ::   nb_itt = 5             ! number of itterations
52
53   !!----------------------------------------------------------------------
54CONTAINS
55
56
57   SUBROUTINE coare3p6_init(l_use_cs, l_use_wl)
58      !!---------------------------------------------------------------------
59      !!                  ***  FUNCTION coare3p6_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 *, ' *** coare3p6_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 ' COARE3P6_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 *, ' *** coare3p6_init: CS => allocating delta_vl :', jpi,jpj
84         ALLOCATE ( delta_vl(jpi,jpj), STAT=ierr )
85         !IF( ierr > 0 ) STOP ' COARE3P6_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 coare3p6_init
90
91
92
93   SUBROUTINE turb_coare3p6( 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_coare3p6  ***
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_skin=True)
130      !!              -> MUST be given the correct value if not computing skint temp. (in case l_use_skin=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_coare3p6@sbcblk_algo_coare3p6.F90'
204      CHARACTER(len=128) :: cf_tmp
205      !!----------------------------------------------------------------------------------
206
207      IF ( kt == 1 ) CALL COARE3P6_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_3p6(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_3p6(ztmp2)*ztmp1/grav + 0.11_wp*znu_a/u_star   ! Roughness length (eq.6) [ ztmp1==u*^2 ]
325         ztmp1 = ( znu_a / (z0*u_star) )**0.72_wp     ! COARE3.6-specific! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific!
326         z0t   = MIN( 1.6E-4_wp , 5.8E-5_wp*ztmp1 )   ! COARE3.6-specific!
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_COARE( 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_COARE 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_COARE( kt, Qsw, ztmp1, zeta_u, zsst, plong, isecday_utc, MOD(nb_itt,j_itt),  pdTw,  Hwl=zHwl )
365            ELSE
366               CALL WL_COARE( 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_coare3p6
414
415
416   FUNCTION alfa_charn_3p6( pwnd )
417      !!-------------------------------------------------------------------
418      !! Computes the Charnock parameter as a function of the Neutral wind speed at 10m
419      !!
420      !!  (Eq. 13 in Edson et al., 2013)
421      !!
422      !! Author: L. Brodeau, July 2019 / AeroBulk  (https://github.com/brodeau/aerobulk/)
423      !!-------------------------------------------------------------------
424      REAL(wp), DIMENSION(jpi,jpj) :: alfa_charn_3p6
425      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd   ! neutral wind speed at 10m
426      !
427      REAL(wp), PARAMETER :: charn0_max = 0.028  !: value above which the Charnock parameter levels off for winds > 18 m/s
428      !!-------------------------------------------------------------------
429      !
430      alfa_charn_3p6 = MAX( MIN( 0.0017_wp*pwnd - 0.005_wp , charn0_max) , 0._wp )
431      !
432   END FUNCTION alfa_charn_3p6
433
434   FUNCTION psi_m_coare( pzeta )
435      !!----------------------------------------------------------------------------------
436      !! ** Purpose: compute the universal profile stability function for momentum
437      !!             COARE 3.0, Fairall et al. 2003
438      !!             pzeta : stability paramenter, z/L where z is altitude
439      !!                     measurement and L is M-O length
440      !!       Stability function for wind speed and scalars matching Kansas and free
441      !!       convection forms with weighting f convective form, follows Fairall et
442      !!       al (1996) with profile constants from Grachev et al (2000) BLM stable
443      !!       form from Beljaars and Holtslag (1991)
444      !!
445      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
446      !!----------------------------------------------------------------------------------
447      REAL(wp), DIMENSION(jpi,jpj) :: psi_m_coare
448      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
449      !
450      INTEGER  ::   ji, jj    ! dummy loop indices
451      REAL(wp) :: zta, zphi_m, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
452      !!----------------------------------------------------------------------------------
453      !
454      DO jj = 1, jpj
455         DO ji = 1, jpi
456            !
457            zta = pzeta(ji,jj)
458            !
459            zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable
460            !
461            zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   &
462               & - 2.*ATAN(zphi_m) + 0.5*rpi
463            !
464            zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective
465            !
466            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
467               &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
468            !
469            zf = zta*zta
470            zf = zf/(1. + zf)
471            zc = MIN(50._wp, 0.35_wp*zta)
472            zstab = 0.5 + SIGN(0.5_wp, zta)
473            !
474            psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0)
475               &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0)
476               &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     "
477            !
478         END DO
479      END DO
480      !
481   END FUNCTION psi_m_coare
482
483
484   FUNCTION psi_h_coare( pzeta )
485      !!---------------------------------------------------------------------
486      !! Universal profile stability function for temperature and humidity
487      !! COARE 3.0, Fairall et al. 2003
488      !!
489      !! pzeta : stability paramenter, z/L where z is altitude measurement
490      !!         and L is M-O length
491      !!
492      !! Stability function for wind speed and scalars matching Kansas and free
493      !! convection forms with weighting f convective form, follows Fairall et
494      !! al (1996) with profile constants from Grachev et al (2000) BLM stable
495      !! form from Beljaars and Holtslag (1991)
496      !!
497      !! Author: L. Brodeau, June 2016 / AeroBulk
498      !!         (https://github.com/brodeau/aerobulk/)
499      !!----------------------------------------------------------------
500      !!
501      REAL(wp), DIMENSION(jpi,jpj) :: psi_h_coare
502      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pzeta
503      !
504      INTEGER  ::   ji, jj     ! dummy loop indices
505      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab
506      !
507      DO jj = 1, jpj
508         DO ji = 1, jpi
509            !
510            zta = pzeta(ji,jj)
511            !
512            zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable)
513            !
514            zpsi_k = 2.*LOG((1. + zphi_h)/2.)
515            !
516            zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective
517            !
518            zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) &
519               &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447
520            !
521            zf = zta*zta
522            zf = zf/(1. + zf)
523            zc = MIN(50._wp,0.35_wp*zta)
524            zstab = 0.5 + SIGN(0.5_wp, zta)
525            !
526            psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) &
527               &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     &
528               &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 )
529            !
530         END DO
531      END DO
532      !
533   END FUNCTION psi_h_coare
534
535   !!======================================================================
536END MODULE sbcblk_algo_coare3p6
Note: See TracBrowser for help on using the repository browser.