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.
sbc_phy.F90 in NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC – NEMO

source: NEMO/branches/2020/dev_r13648_ASINTER-04_laurent_bulk_ice/src/OCE/SBC/sbc_phy.F90 @ 13655

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

Commit all my dev of 2020!
File size: 58.6 KB
Line 
1MODULE sbc_phy
2   !!======================================================================
3   !!                       ***  MODULE  sbc_phy  ***
4   !! A set of functions to compute air themodynamics parameters
5   !!                     needed by Aerodynamic Bulk Formulas
6   !!=====================================================================
7   !!            4.x  !  2020 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
8   !!----------------------------------------------------------------------
9
10   !!   virt_temp     : virtual (aka sensible) temperature (potential or absolute)
11   !!   rho_air       : density of (moist) air (depends on T_air, q_air and SLP
12   !!   visc_air      : kinematic viscosity (aka Nu_air) of air from temperature
13   !!   L_vap         : latent heat of vaporization of water as a function of temperature
14   !!   cp_air        : specific heat of (moist) air (depends spec. hum. q_air)
15   !!   gamma_moist   : adiabatic lapse-rate of moist air
16   !!   One_on_L      : 1. / ( Obukhov length )
17   !!   Ri_bulk       : bulk Richardson number aka BRN
18   !!   q_sat         : saturation humidity as a function of SLP and temperature
19   !!   q_air_rh      : specific humidity as a function of RH (fraction, not %), t_air and SLP
20
21   USE dom_oce        ! ocean space and time domain
22   USE phycst         ! physical constants
23
24   IMPLICIT NONE
25   PRIVATE
26
27   INTEGER , PARAMETER, PUBLIC  :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument)
28
29   !!  (mainly removed from sbcblk.F90)
30   REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp   !: Specic heat of dry air, constant pressure      [J/K/kg]
31   REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp   !: Specic heat of water vapor, constant pressure  [J/K/kg]
32   REAL(wp), PARAMETER, PUBLIC :: R_dry   = 287.05_wp   !: Specific gas constant for dry air              [J/K/kg]
33   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg]
34   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
35   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp  !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
36   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...)
37   REAL(wp), PARAMETER, PUBLIC :: albo    = 0.066_wp    !: ocean albedo assumed to be constant
38   !
39   REAL(wp), PARAMETER, PUBLIC :: rho0_a  = 1.2_wp      !: Approx. of density of air                       [kg/m^3]
40   REAL(wp), PARAMETER, PUBLIC :: rho0_w  = 1025._wp    !: Density of sea-water  (ECMWF->1025)             [kg/m^3]
41   REAL(wp), PARAMETER, PUBLIC :: radrw = rho0_a/rho0_w !: Density ratio
42   REAL(wp), PARAMETER, PUBLIC :: sq_radrw = SQRT(rho0_a/rho0_w)
43   REAL(wp), PARAMETER, PUBLIC :: rCp0_w  = 4190._wp    !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg]
44   REAL(wp), PARAMETER, PUBLIC :: rnu0_w  = 1.e-6_wp    !: kinetic viscosity of water                      [m^2/s]
45   REAL(wp), PARAMETER, PUBLIC :: rk0_w   = 0.6_wp      !: thermal conductivity of water (at 20C)          [W/m/K]
46   !
47   REAL(wp), PARAMETER, PUBLIC :: emiss_w = 0.98_wp     !: Long-wave (thermal) emissivity of sea-water []
48   !
49   REAL(wp), PARAMETER, PUBLIC :: emiss_i = 0.996_wp    !:  "   for ice and snow => but Rees 1993 suggests can be lower in winter on fresh snow... 0.72 ...
50
51   REAL(wp), PARAMETER, PUBLIC :: wspd_thrshld_ice = 0.2_wp !: minimum scalar wind speed accepted over sea-ice... [m/s]
52
53   !
54   REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp  !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt
55   REAL(wp), PARAMETER, PUBLIC :: rtt0 = 273.16_wp        !: triple point of temperature    [K]
56   !
57   REAL(wp), PARAMETER, PUBLIC :: rcst_cs = -16._wp*9.80665_wp*rho0_w*rCp0_w*rnu0_w*rnu0_w*rnu0_w/(rk0_w*rk0_w) !: for cool-skin parameterizations... (grav = 9.80665_wp)
58   !                              => see eq.(14) in Fairall et al. 1996   (eq.(6) of Zeng aand Beljaars is WRONG! (typo?)
59
60   REAL(wp), PARAMETER, PUBLIC :: z0_sea_max = 0.0025_wp   !: maximum realistic value for roughness length of sea-surface... [m]
61
62   REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-]
63
64
65   REAL(wp), PARAMETER, PUBLIC :: Cx_min = 0.1E-3_wp ! smallest value allowed for bulk transfer coefficients (usually in stable conditions with now wind)
66
67   REAL(wp), PARAMETER :: &
68                                !! Constants for Goff formula in the presence of ice:
69      &      rAg_i = -9.09718_wp, &
70      &      rBg_i = -3.56654_wp, &
71      &      rCg_i = 0.876793_wp, &
72      &      rDg_i = LOG10(6.1071_wp)
73
74   REAL(wp), PARAMETER :: rc_louis  = 5._wp
75   REAL(wp), PARAMETER :: rc2_louis = rc_louis * rc_louis
76   REAL(wp), PARAMETER :: ram_louis = 2. * rc_louis
77   REAL(wp), PARAMETER :: rah_louis = 3. * rc_louis
78
79
80   INTERFACE virt_temp
81      MODULE PROCEDURE virt_temp_vctr, virt_temp_sclr
82   END INTERFACE virt_temp
83
84   INTERFACE visc_air
85      MODULE PROCEDURE visc_air_vctr, visc_air_sclr
86   END INTERFACE visc_air
87
88   INTERFACE gamma_moist
89      MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
90   END INTERFACE gamma_moist
91
92   INTERFACE e_sat
93      MODULE PROCEDURE e_sat_vctr, e_sat_sclr
94   END INTERFACE e_sat
95
96   INTERFACE e_sat_ice
97      MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr
98   END INTERFACE e_sat_ice
99   INTERFACE de_sat_dt_ice
100      MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr
101   END INTERFACE de_sat_dt_ice
102
103   INTERFACE Ri_bulk
104      MODULE PROCEDURE Ri_bulk_vctr, Ri_bulk_sclr
105   END INTERFACE Ri_bulk
106
107   INTERFACE q_sat
108      MODULE PROCEDURE q_sat_vctr, q_sat_sclr
109   END INTERFACE q_sat
110
111   INTERFACE dq_sat_dt_ice
112      MODULE PROCEDURE dq_sat_dt_ice_vctr, dq_sat_dt_ice_sclr
113   END INTERFACE dq_sat_dt_ice
114
115   INTERFACE L_vap
116      MODULE PROCEDURE L_vap_vctr, L_vap_sclr
117   END INTERFACE L_vap
118
119   INTERFACE rho_air
120      MODULE PROCEDURE rho_air_vctr, rho_air_sclr
121   END INTERFACE rho_air
122
123   INTERFACE cp_air
124      MODULE PROCEDURE cp_air_vctr, cp_air_sclr
125   END INTERFACE cp_air
126
127   INTERFACE alpha_sw
128      MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr
129   END INTERFACE alpha_sw
130
131   INTERFACE bulk_formula
132      MODULE PROCEDURE bulk_formula_vctr, bulk_formula_sclr
133   END INTERFACE bulk_formula
134
135   INTERFACE qlw_net
136      MODULE PROCEDURE qlw_net_vctr, qlw_net_sclr
137   END INTERFACE qlw_net
138
139   INTERFACE f_m_louis
140      MODULE PROCEDURE f_m_louis_vctr, f_m_louis_sclr
141   END INTERFACE f_m_louis
142
143   INTERFACE f_h_louis
144      MODULE PROCEDURE f_h_louis_vctr, f_h_louis_sclr
145   END INTERFACE f_h_louis
146
147
148   PUBLIC virt_temp
149   PUBLIC rho_air
150   PUBLIC visc_air
151   PUBLIC L_vap
152   PUBLIC cp_air
153   PUBLIC gamma_moist
154   PUBLIC One_on_L
155   PUBLIC Ri_bulk
156   PUBLIC q_sat
157   PUBLIC q_air_rh
158   PUBLIC dq_sat_dt_ice
159   !:
160   PUBLIC update_qnsol_tau
161   PUBLIC alpha_sw
162   PUBLIC bulk_formula
163   PUBLIC qlw_net
164   !
165   PUBLIC f_m_louis, f_h_louis
166   PUBLIC z0_from_Cd
167   PUBLIC Cd_from_z0
168   PUBLIC UN10_from_ustar
169   PUBLIC UN10_from_CD
170   PUBLIC z0tq_LKB
171
172   !! * Substitutions
173#  include "do_loop_substitute.h90"
174   !!----------------------------------------------------------------------
175   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
176   !! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
177   !! Software governed by the CeCILL license (see ./LICENSE)
178   !!----------------------------------------------------------------------
179CONTAINS
180
181
182   FUNCTION virt_temp_sclr( pta, pqa )
183      !!------------------------------------------------------------------------
184      !!
185      !! Compute the (absolute/potential) VIRTUAL temperature, based on the
186      !! (absolute/potential) temperature and specific humidity
187      !!
188      !! If input temperature is absolute then output virtual temperature is absolute
189      !! If input temperature is potential then output virtual temperature is potential
190      !!
191      !! Author: L. Brodeau, June 2019 / AeroBulk
192      !!         (https://github.com/brodeau/aerobulk/)
193      !!------------------------------------------------------------------------
194      REAL(wp)             :: virt_temp_sclr !: virtual temperature [K]
195      REAL(wp), INTENT(in) :: pta       !: absolute or potential air temperature [K]
196      REAL(wp), INTENT(in) :: pqa       !: specific humidity of air   [kg/kg]
197      !!-------------------------------------------------------------------
198      !
199      virt_temp_sclr = pta * (1._wp + rctv0*pqa)
200      !!
201      !! This is exactly the same thing as:
202      !! virt_temp_sclr = pta * ( pwa + reps0) / (reps0*(1.+pwa))
203      !! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
204      !
205   END FUNCTION virt_temp_sclr
206   !!
207   FUNCTION virt_temp_vctr( pta, pqa )
208      REAL(wp), DIMENSION(jpi,jpj)             :: virt_temp_vctr !: virtual temperature [K]
209      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
210      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air   [kg/kg]
211      virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
212   END FUNCTION virt_temp_vctr
213   !===============================================================================================
214
215
216   FUNCTION rho_air_vctr( ptak, pqa, ppa )
217      !!-------------------------------------------------------------------------------
218      !!                           ***  FUNCTION rho_air_vctr  ***
219      !!
220      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
221      !!
222      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
223      !!-------------------------------------------------------------------------------
224      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak      ! air temperature             [K]
225      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa       ! air specific humidity   [kg/kg]
226      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ppa      ! pressure in                [Pa]
227      REAL(wp), DIMENSION(jpi,jpj)             ::   rho_air_vctr   ! density of moist air   [kg/m^3]
228      !!-------------------------------------------------------------------------------
229      rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
230   END FUNCTION rho_air_vctr
231
232   FUNCTION rho_air_sclr( ptak, pqa, ppa )
233      !!-------------------------------------------------------------------------------
234      !!                           ***  FUNCTION rho_air_sclr  ***
235      !!
236      !! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
237      !!
238      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
239      !!-------------------------------------------------------------------------------
240      REAL(wp), INTENT(in) :: ptak           ! air temperature             [K]
241      REAL(wp), INTENT(in) :: pqa            ! air specific humidity   [kg/kg]
242      REAL(wp), INTENT(in) :: ppa           ! pressure in                [Pa]
243      REAL(wp)             :: rho_air_sclr   ! density of moist air   [kg/m^3]
244      !!-------------------------------------------------------------------------------
245      rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
246   END FUNCTION rho_air_sclr
247
248
249
250
251   FUNCTION visc_air_sclr(ptak)
252      !!----------------------------------------------------------------------------------
253      !! Air kinetic viscosity (m^2/s) given from air temperature in Kelvin
254      !!
255      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
256      !!----------------------------------------------------------------------------------
257      REAL(wp)             :: visc_air_sclr   ! kinetic viscosity (m^2/s)
258      REAL(wp), INTENT(in) :: ptak       ! air temperature in (K)
259      !
260      REAL(wp) ::   ztc, ztc2   ! local scalar
261      !!----------------------------------------------------------------------------------
262      !
263      ztc  = ptak - rt0   ! air temp, in deg. C
264      ztc2 = ztc*ztc
265      visc_air_sclr = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
266      !
267   END FUNCTION visc_air_sclr
268
269   FUNCTION visc_air_vctr(ptak)
270      REAL(wp), DIMENSION(jpi,jpj)             ::   visc_air_vctr   ! kinetic viscosity (m^2/s)
271      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak       ! air temperature in (K)
272      INTEGER  ::   ji, jj      ! dummy loop indices
273      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
274      visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
275      END_2D
276   END FUNCTION visc_air_vctr
277
278
279   FUNCTION L_vap_vctr( psst )
280      !!---------------------------------------------------------------------------------
281      !!                           ***  FUNCTION L_vap_vctr  ***
282      !!
283      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
284      !!
285      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
286      !!----------------------------------------------------------------------------------
287      REAL(wp), DIMENSION(jpi,jpj)             ::   L_vap_vctr   ! latent heat of vaporization   [J/kg]
288      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
289      !!----------------------------------------------------------------------------------
290      !
291      L_vap_vctr = (  2.501_wp - 0.00237_wp * ( psst(:,:) - rt0)  ) * 1.e6_wp
292      !
293   END FUNCTION L_vap_vctr
294
295   FUNCTION L_vap_sclr( psst )
296      !!---------------------------------------------------------------------------------
297      !!                           ***  FUNCTION L_vap_sclr  ***
298      !!
299      !! ** Purpose : Compute the latent heat of vaporization of water from temperature
300      !!
301      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
302      !!----------------------------------------------------------------------------------
303      REAL(wp)             ::   L_vap_sclr   ! latent heat of vaporization   [J/kg]
304      REAL(wp), INTENT(in) ::   psst         ! water temperature                [K]
305      !!----------------------------------------------------------------------------------
306      !
307      L_vap_sclr = (  2.501_wp - 0.00237_wp * ( psst - rt0)  ) * 1.e6_wp
308      !
309   END FUNCTION L_vap_sclr
310
311
312   FUNCTION cp_air_vctr( pqa )
313      !!-------------------------------------------------------------------------------
314      !!                           ***  FUNCTION cp_air_vctr  ***
315      !!
316      !! ** Purpose : Compute specific heat (Cp) of moist air
317      !!
318      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
319      !!-------------------------------------------------------------------------------
320      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa      ! air specific humidity         [kg/kg]
321      REAL(wp), DIMENSION(jpi,jpj)             ::   cp_air_vctr   ! specific heat of moist air   [J/K/kg]
322      !!-------------------------------------------------------------------------------
323      cp_air_vctr = rCp_dry + rCp_vap * pqa
324   END FUNCTION cp_air_vctr
325
326   FUNCTION cp_air_sclr( pqa )
327      !!-------------------------------------------------------------------------------
328      !!                           ***  FUNCTION cp_air_sclr  ***
329      !!
330      !! ** Purpose : Compute specific heat (Cp) of moist air
331      !!
332      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
333      !!-------------------------------------------------------------------------------
334      REAL(wp), INTENT(in) :: pqa           ! air specific humidity         [kg/kg]
335      REAL(wp)             :: cp_air_sclr   ! specific heat of moist air   [J/K/kg]
336      !!-------------------------------------------------------------------------------
337      cp_air_sclr = rCp_dry + rCp_vap * pqa
338   END FUNCTION cp_air_sclr
339
340
341
342   !===============================================================================================
343   FUNCTION gamma_moist_sclr( ptak, pqa )
344      !!----------------------------------------------------------------------------------
345      !! ** Purpose : Compute the moist adiabatic lapse-rate.
346      !!     => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
347      !!     => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
348      !!
349      !! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
350      !!----------------------------------------------------------------------------------
351      REAL(wp)             :: gamma_moist_sclr !                           [K/m]
352      REAL(wp), INTENT(in) ::   ptak           ! absolute air temperature  [K] !#LB: double check it's absolute !!!
353      REAL(wp), INTENT(in) ::   pqa            ! specific humidity     [kg/kg]
354      !
355      REAL(wp) :: zta, zqa, zwa, ziRT, zLvap        ! local scalars
356      !!----------------------------------------------------------------------------------
357      zta = MAX( ptak,  180._wp) ! prevents screw-up over masked regions where field == 0.
358      zqa = MAX( pqa,  1.E-6_wp) !    "                   "                     "
359      !!
360      zwa = zqa / (1._wp - zqa)   ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
361      ziRT = 1._wp / (R_dry*zta)    ! 1/RT
362      zLvap = L_vap_sclr( ptak )
363      !!
364      gamma_moist_sclr = grav * ( 1._wp + zLvap*zwa*ziRT ) / ( rCp_dry + zLvap*zLvap*zwa*reps0*ziRT/zta )
365      !!
366   END FUNCTION gamma_moist_sclr
367   !!
368   FUNCTION gamma_moist_vctr( ptak, pqa )
369      REAL(wp), DIMENSION(jpi,jpj)             ::   gamma_moist_vctr
370      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   ptak
371      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pqa
372      INTEGER  :: ji, jj
373      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
374      gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
375      END_2D
376   END FUNCTION gamma_moist_vctr
377   !===============================================================================================
378
379
380   FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
381      !!------------------------------------------------------------------------
382      !!
383      !! Evaluates the 1./(Obukhov length) from air temperature,
384      !! air specific humidity, and frictional scales u*, t* and q*
385      !!
386      !! Author: L. Brodeau, June 2019 / AeroBulk
387      !!         (https://github.com/brodeau/aerobulk/)
388      !!------------------------------------------------------------------------
389      REAL(wp), DIMENSION(jpi,jpj)             :: One_on_L     !: 1./(Obukhov length) [m^-1]
390      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha         !: reference potential temperature of air [K]
391      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa          !: reference specific humidity of air   [kg/kg]
392      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus          !: u*: friction velocity [m/s]
393      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs     !: \theta* and q* friction aka turb. scales for temp. and spec. hum.
394      !
395      INTEGER  ::   ji, jj         ! dummy loop indices
396      REAL(wp) ::     zqa          ! local scalar
397      !!-------------------------------------------------------------------
398      !
399      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
400      !
401      zqa = (1._wp + rctv0*pqa(ji,jj))
402      !
403      ! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
404      !  a/  -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
405      !                      or
406      !  b/  -u* [ theta*              + 0.61 theta q* ]
407      !
408      One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
409         &               / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
410      !
411      END_2D
412      !
413      One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
414      !
415   END FUNCTION One_on_L
416
417
418   !===============================================================================================
419   FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub,  pta_layer, pqa_layer )
420      !!----------------------------------------------------------------------------------
421      !! Bulk Richardson number according to "wide-spread equation"...
422      !!
423      !!    Reminder: the Richardson number is the ratio "buoyancy" / "shear"
424      !!
425      !! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
426      !!----------------------------------------------------------------------------------
427      REAL(wp)             :: Ri_bulk_sclr
428      REAL(wp), INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
429      REAL(wp), INTENT(in) :: psst  ! SST                                   [K]
430      REAL(wp), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
431      REAL(wp), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
432      REAL(wp), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
433      REAL(wp), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
434      REAL(wp), INTENT(in), OPTIONAL :: pta_layer ! when possible, a better guess of absolute temperature WITHIN the layer [K]
435      REAL(wp), INTENT(in), OPTIONAL :: pqa_layer ! when possible, a better guess of specific humidity    WITHIN the layer [kg/kg]
436      !!
437      LOGICAL  :: l_ptqa_l_prvd = .FALSE.
438      REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv  ! local scalars
439      !!-------------------------------------------------------------------
440      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
441      !
442      zsstv = virt_temp_sclr( psst, pssq )          ! virtual SST (absolute==potential because z=0!)
443      !
444      zdthv = virt_temp_sclr( ptha, pqa  ) - zsstv  ! air-sea delta of "virtual potential temperature"
445      !
446      !! ztv: estimate of the ABSOLUTE virtual temp. within the layer
447      IF( l_ptqa_l_prvd ) THEN
448         ztv = virt_temp_sclr( pta_layer, pqa_layer )
449      ELSE
450         zqa = 0.5_wp*( pqa  + pssq )                             ! ~ mean q within the layer...
451         zta = 0.5_wp*( psst + ptha - gamma_moist(ptha, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
452         zta = 0.5_wp*( psst + ptha - gamma_moist( zta, zqa)*pz ) ! ~ mean absolute temperature of air within the layer
453         zgamma =  gamma_moist(zta, zqa)                          ! Adiabatic lapse-rate for moist air within the layer
454         ztv = 0.5_wp*( zsstv + virt_temp_sclr( ptha-zgamma*pz, pqa ) )
455      END IF
456      !
457      Ri_bulk_sclr = grav*zdthv*pz / ( ztv*pub*pub )      ! the usual definition of Ri_bulk_sclr
458      !
459   END FUNCTION Ri_bulk_sclr
460   !!
461   FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub,  pta_layer, pqa_layer )
462      REAL(wp), DIMENSION(jpi,jpj)             :: Ri_bulk_vctr
463      REAL(wp)                    , INTENT(in) :: pz    ! height above the sea (aka "delta z")  [m]
464      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst  ! SST                                   [K]
465      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha  ! pot. air temp. at height "pz"         [K]
466      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq  ! 0.98*q_sat(SST)                   [kg/kg]
467      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa   ! air spec. hum. at height "pz"     [kg/kg]
468      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub   ! bulk wind speed                     [m/s]
469      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: pta_layer ! when possible, a better guess of absolute temperature WITHIN the layer [K]
470      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: pqa_layer ! when possible, a better guess of specific humidity    WITHIN the layer [kg/kg]
471      !!
472      LOGICAL  :: l_ptqa_l_prvd = .FALSE.
473      INTEGER  ::   ji, jj
474      IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd=.TRUE.
475      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
476      IF( l_ptqa_l_prvd ) THEN
477         Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
478            &                                pta_layer=pta_layer(ji,jj ),  pqa_layer=pqa_layer(ji,jj ) )
479      ELSE
480         Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
481      END IF
482      END_2D
483   END FUNCTION Ri_bulk_vctr
484   !===============================================================================================
485
486   !===============================================================================================
487   FUNCTION e_sat_sclr( ptak )
488      !!----------------------------------------------------------------------------------
489      !!                   ***  FUNCTION e_sat_sclr  ***
490      !!                  < SCALAR argument version >
491      !! ** Purpose : water vapor at saturation in [Pa]
492      !!              Based on accurate estimate by Goff, 1957
493      !!
494      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
495      !!
496      !!    Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
497      !!----------------------------------------------------------------------------------
498      REAL(wp)             ::   e_sat_sclr   ! water vapor at saturation   [kg/kg]
499      REAL(wp), INTENT(in) ::   ptak    ! air temperature                  [K]
500      REAL(wp) ::   zta, ztmp   ! local scalar
501      !!----------------------------------------------------------------------------------
502      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
503      ztmp = rt0 / zta   !#LB: rt0 or rtt0 ???? (273.15 vs 273.16 )
504      !
505      ! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
506      e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0)        &
507         &    + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) )  &
508         &    + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
509      !
510   END FUNCTION e_sat_sclr
511   !!
512   FUNCTION e_sat_vctr(ptak)
513      REAL(wp), DIMENSION(jpi,jpj)             :: e_sat_vctr !: vapour pressure at saturation  [Pa]
514      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak    !: temperature (K)
515      INTEGER  ::   ji, jj         ! dummy loop indices
516      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
517      e_sat_vctr(ji,jj) = e_sat_sclr(ptak(ji,jj))
518      END_2D
519   END FUNCTION e_sat_vctr
520   !===============================================================================================
521
522   !===============================================================================================
523   FUNCTION e_sat_ice_sclr(ptak)
524      !!---------------------------------------------------------------------------------
525      !! Same as "e_sat" but over ice rather than water!
526      !!---------------------------------------------------------------------------------
527      REAL(wp)             :: e_sat_ice_sclr !: vapour pressure at saturation in presence of ice [Pa]
528      REAL(wp), INTENT(in) :: ptak
529      !!
530      REAL(wp) :: zta, zle, ztmp
531      !!---------------------------------------------------------------------------------
532      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
533      ztmp = rtt0/zta
534      !!
535      zle  = rAg_i*(ztmp - 1._wp) + rBg_i*LOG10(ztmp) + rCg_i*(1._wp - zta/rtt0) + rDg_i
536      !!
537      e_sat_ice_sclr = 100._wp * 10._wp**zle
538   END FUNCTION e_sat_ice_sclr
539   !!
540   FUNCTION e_sat_ice_vctr(ptak)
541      !! Same as "e_sat" but over ice rather than water!
542      REAL(wp), DIMENSION(jpi,jpj) :: e_sat_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
543      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
544      INTEGER  :: ji, jj
545      !!----------------------------------------------------------------------------------
546      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
547      e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
548      END_2D
549   END FUNCTION e_sat_ice_vctr
550   !!
551   FUNCTION de_sat_dt_ice_sclr(ptak)
552      !!---------------------------------------------------------------------------------
553      !! d [ e_sat_ice ] / dT   (derivative / temperature)
554      !! Analytical exact formulation: double checked!!!
555      !!  => DOUBLE-check possible / finite-difference version with "./bin/test_phymbl.x"
556      !!---------------------------------------------------------------------------------
557      REAL(wp)             :: de_sat_dt_ice_sclr !:  [Pa/K]
558      REAL(wp), INTENT(in) :: ptak
559      !!
560      REAL(wp) :: zta, zde
561      !!---------------------------------------------------------------------------------
562      zta = MAX( ptak , 180._wp )   ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
563      !!
564      zde = -(rAg_i*rtt0)/(zta*zta) - rBg_i/(zta*LOG(10._wp)) - rCg_i/rtt0
565      !!
566      de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta)
567   END FUNCTION de_sat_dt_ice_sclr
568   !!
569   FUNCTION de_sat_dt_ice_vctr(ptak)
570      !! Same as "e_sat" but over ice rather than water!
571      REAL(wp), DIMENSION(jpi,jpj) :: de_sat_dt_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
572      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
573      INTEGER  :: ji, jj
574      !!----------------------------------------------------------------------------------
575      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
576      de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
577      END_2D
578   END FUNCTION de_sat_dt_ice_vctr
579
580
581
582   !===============================================================================================
583
584   !===============================================================================================
585   FUNCTION q_sat_sclr( pta, ppa,  l_ice )
586      !!---------------------------------------------------------------------------------
587      !!                           ***  FUNCTION q_sat_sclr  ***
588      !!
589      !! ** Purpose : Conputes specific humidity of air at saturation
590      !!
591      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
592      !!----------------------------------------------------------------------------------
593      REAL(wp) :: q_sat_sclr
594      REAL(wp), INTENT(in) :: pta  !: absolute temperature of air [K]
595      REAL(wp), INTENT(in) :: ppa  !: atmospheric pressure        [Pa]
596      LOGICAL,  INTENT(in), OPTIONAL :: l_ice  !: we are above ice
597      REAL(wp) :: ze_s
598      LOGICAL  :: lice
599      !!----------------------------------------------------------------------------------
600      lice = .FALSE.
601      IF( PRESENT(l_ice) ) lice = l_ice
602      IF( lice ) THEN
603         ze_s = e_sat_ice( pta )
604      ELSE
605         ze_s = e_sat( pta ) ! Vapour pressure at saturation (Goff) :
606      END IF
607      q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s)
608   END FUNCTION q_sat_sclr
609   !!
610   FUNCTION q_sat_vctr( pta, ppa,  l_ice )
611      REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr
612      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta  !: absolute temperature of air [K]
613      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa  !: atmospheric pressure        [Pa]
614      LOGICAL,  INTENT(in), OPTIONAL :: l_ice  !: we are above ice
615      LOGICAL  :: lice
616      INTEGER  :: ji, jj
617      !!----------------------------------------------------------------------------------
618      lice = .FALSE.
619      IF( PRESENT(l_ice) ) lice = l_ice
620      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
621      q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
622      END_2D
623   END FUNCTION q_sat_vctr
624   !===============================================================================================
625
626
627   !===============================================================================================
628   FUNCTION dq_sat_dt_ice_sclr( pta, ppa )
629      !!---------------------------------------------------------------------------------
630      !!     ***  FUNCTION dq_sat_dt_ice_sclr  ***
631      !!    => d [ q_sat_ice(T) ] / dT
632      !! Analytical exact formulation: double checked!!!
633      !!  => DOUBLE-check possible / finite-difference version with "./bin/test_phymbl.x"
634      !!----------------------------------------------------------------------------------
635      REAL(wp) :: dq_sat_dt_ice_sclr
636      REAL(wp), INTENT(in) :: pta  !: absolute temperature of air [K]
637      REAL(wp), INTENT(in) :: ppa  !: atmospheric pressure        [Pa]
638      REAL(wp) :: ze_s, zde_s_dt, ztmp
639      !!----------------------------------------------------------------------------------
640      ze_s     =  e_sat_ice_sclr( pta ) ! Vapour pressure at saturation  in presence of ice (Goff)
641      zde_s_dt = de_sat_dt_ice(   pta )
642      !
643      ztmp = (reps0 - 1._wp)*ze_s + ppa
644      !
645      dq_sat_dt_ice_sclr = reps0*ppa*zde_s_dt / ( ztmp*ztmp )
646      !
647   END FUNCTION dq_sat_dt_ice_sclr
648   !!
649   FUNCTION dq_sat_dt_ice_vctr( pta, ppa )
650      REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr
651      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta  !: absolute temperature of air [K]
652      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa  !: atmospheric pressure        [Pa]
653      INTEGER  :: ji, jj
654      !!----------------------------------------------------------------------------------
655      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
656      dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
657      END_2D
658   END FUNCTION dq_sat_dt_ice_vctr
659   !===============================================================================================
660
661
662   !===============================================================================================
663   FUNCTION q_air_rh(prha, ptak, ppa)
664      !!----------------------------------------------------------------------------------
665      !! Specific humidity of air out of Relative Humidity
666      !!
667      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
668      !!----------------------------------------------------------------------------------
669      REAL(wp), DIMENSION(jpi,jpj)             :: q_air_rh
670      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha        !: relative humidity      [fraction, not %!!!]
671      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak        !: air temperature        [K]
672      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa        !: atmospheric pressure   [Pa]
673      !
674      INTEGER  ::   ji, jj      ! dummy loop indices
675      REAL(wp) ::   ze      ! local scalar
676      !!----------------------------------------------------------------------------------
677      !
678      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
679      ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
680      q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
681      END_2D
682      !
683   END FUNCTION q_air_rh
684
685
686   SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, &
687      &                         pQns, pTau,  &
688      &                         Qlat)
689      !!----------------------------------------------------------------------------------
690      !! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
691      !!          and the module of the wind stress => pTau = Tau
692      !! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
693      !!----------------------------------------------------------------------------------
694      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
695      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
696      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
697      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
698      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
699      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pust ! u*
700      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ptst ! t*
701      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqst ! q*
702      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
703      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
704      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa]
705      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: prlw ! downwelling longwave radiative flux [W/m^2]
706      !
707      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
708      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
709      !
710      REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
711      !
712      REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zz0, zQlat, zQsen, zQlw
713      INTEGER  ::   ji, jj     ! dummy loop indices
714      !!----------------------------------------------------------------------------------
715      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
716      zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
717      zdq = pqa(ji,jj) - pqs(ji,jj) ;  zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
718      zz0 = pust(ji,jj)/pUb(ji,jj)
719      zCd = zz0*zz0
720      zCh = zz0*ptst(ji,jj)/zdt
721      zCe = zz0*pqst(ji,jj)/zdq
722
723      CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
724         &              pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), &
725         &              pTau(ji,jj), zQsen, zQlat )
726
727      zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
728
729      pQns(ji,jj) = zQlat + zQsen + zQlw
730
731      IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
732      END_2D
733   END SUBROUTINE UPDATE_QNSOL_TAU
734
735
736   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
737      &                          pCd, pCh, pCe,           &
738      &                          pwnd, pUb, ppa,         &
739      &                          pTau, pQsen, pQlat,      &
740      &                          pEvap, prhoa, pfact_evap )
741      !!----------------------------------------------------------------------------------
742      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
743      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
744      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
745      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
746      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
747      REAL(wp), INTENT(in)  :: pCd
748      REAL(wp), INTENT(in)  :: pCh
749      REAL(wp), INTENT(in)  :: pCe
750      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
751      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
752      REAL(wp), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa]
753      !!
754      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
755      REAL(wp), INTENT(out) :: pQsen !  [W/m^2]
756      REAL(wp), INTENT(out) :: pQlat !  [W/m^2]
757      !!
758      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
759      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
760      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
761      !!
762      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
763      INTEGER  :: jq
764      !!----------------------------------------------------------------------------------
765      zfact_evap = 1._wp
766      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
767
768      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa")
769      ztaa = pTa ! first guess...
770      DO jq = 1, 4
771         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )  !#LB: why not "0.5*(pqs+pqa)" rather then "pqa" ???
772         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder...
773      END DO
774      zrho = rho_air(ztaa, pqa, ppa)
775      zrho = rho_air(ztaa, pqa, ppa-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given!
776
777      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10
778
779      pTau = zUrho * pCd * pwnd ! Wind stress module
780
781      zevap = zUrho * pCe * (pqa - pqs)
782      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
783      pQlat = L_vap(pTs) * zevap
784
785      IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
786      IF( PRESENT(prhoa) ) prhoa = zrho
787
788   END SUBROUTINE BULK_FORMULA_SCLR
789   !!
790   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
791      &                          pCd, pCh, pCe,           &
792      &                          pwnd, pUb, ppa,         &
793      &                          pTau, pQsen, pQlat,      &
794      &                          pEvap, prhoa, pfact_evap )
795      !!----------------------------------------------------------------------------------
796      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m)
797      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K]
798      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg]
799      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K]
800      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg]
801      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCd
802      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCh
803      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pCe
804      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s]
805      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
806      REAL(wp), DIMENSION(jpi,jpj), INTENT(in)  :: ppa ! sea-level atmospheric pressure [Pa]
807      !!
808      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau  ! module of the wind stress [N/m^2]
809      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen !  [W/m^2]
810      REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat !  [W/m^2]
811      !!
812      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
813      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3]
814      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
815      !!
816      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
817      INTEGER  :: ji, jj
818      !!----------------------------------------------------------------------------------
819      zfact_evap = 1._wp
820      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
821
822      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
823
824      CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
825         &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  &
826         &                    pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj),                &
827         &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             &
828         &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       )
829
830      IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
831      IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho
832      END_2D
833   END SUBROUTINE BULK_FORMULA_VCTR
834
835
836   FUNCTION alpha_sw_vctr( psst )
837      !!---------------------------------------------------------------------------------
838      !!                           ***  FUNCTION alpha_sw_vctr  ***
839      !!
840      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
841      !!
842      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
843      !!----------------------------------------------------------------------------------
844      REAL(wp), DIMENSION(jpi,jpj)             ::   alpha_sw_vctr   ! thermal expansion coefficient of sea-water [1/K]
845      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   psst   ! water temperature                [K]
846      !!----------------------------------------------------------------------------------
847      alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
848   END FUNCTION alpha_sw_vctr
849
850   FUNCTION alpha_sw_sclr( psst )
851      !!---------------------------------------------------------------------------------
852      !!                           ***  FUNCTION alpha_sw_sclr  ***
853      !!
854      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
855      !!
856      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
857      !!----------------------------------------------------------------------------------
858      REAL(wp)             ::   alpha_sw_sclr   ! thermal expansion coefficient of sea-water [1/K]
859      REAL(wp), INTENT(in) ::   psst   ! sea-water temperature                   [K]
860      !!----------------------------------------------------------------------------------
861      alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
862   END FUNCTION alpha_sw_sclr
863
864
865   !===============================================================================================
866   FUNCTION qlw_net_sclr( pdwlw, pts,  l_ice )
867      !!---------------------------------------------------------------------------------
868      !!                           ***  FUNCTION qlw_net_sclr  ***
869      !!
870      !! ** Purpose : Estimate of the net longwave flux at the surface
871      !!----------------------------------------------------------------------------------
872      REAL(wp) :: qlw_net_sclr
873      REAL(wp), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2]
874      REAL(wp), INTENT(in) :: pts   !: surface temperature [K]
875      LOGICAL,  INTENT(in), OPTIONAL :: l_ice  !: we are above ice
876      REAL(wp) :: zemiss, zt2
877      LOGICAL  :: lice
878      !!----------------------------------------------------------------------------------
879      lice = .FALSE.
880      IF( PRESENT(l_ice) ) lice = l_ice
881      IF( lice ) THEN
882         zemiss = emiss_i
883      ELSE
884         zemiss = emiss_w
885      END IF
886      zt2 = pts*pts
887      qlw_net_sclr = zemiss*( pdwlw - stefan*zt2*zt2)  ! zemiss used both as the IR albedo and IR emissivity...
888   END FUNCTION qlw_net_sclr
889   !!
890   FUNCTION qlw_net_vctr( pdwlw, pts,  l_ice )
891      REAL(wp), DIMENSION(jpi,jpj) :: qlw_net_vctr
892      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pdwlw !: downwelling longwave (aka infrared, aka thermal) radiation [W/m^2]
893      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts   !: surface temperature [K]
894      LOGICAL,  INTENT(in), OPTIONAL :: l_ice  !: we are above ice
895      LOGICAL  :: lice
896      INTEGER  :: ji, jj
897      !!----------------------------------------------------------------------------------
898      lice = .FALSE.
899      IF( PRESENT(l_ice) ) lice = l_ice
900      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
901      qlw_net_vctr(ji,jj) = qlw_net_sclr( pdwlw(ji,jj) , pts(ji,jj), l_ice=lice )
902      END_2D
903   END FUNCTION qlw_net_vctr
904   !===============================================================================================
905
906   FUNCTION z0_from_Cd( pzu, pCd,  ppsi )
907      REAL(wp), DIMENSION(jpi,jpj) :: z0_from_Cd        !: roughness length [m]
908      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m]
909      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd   !: (neutral or non-neutral) drag coefficient []
910      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
911      !!
912      !! If pCd is the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given
913      !! If pCd is the drag coefficient (in stable or unstable conditions) then pssi must be provided
914      !!----------------------------------------------------------------------------------
915      IF( PRESENT(ppsi) ) THEN
916         !! Cd provided is the actual Cd (not the neutral-stability CdN) :
917         z0_from_Cd = pzu * EXP( - ( vkarmn/SQRT(pCd(:,:)) + ppsi(:,:) ) ) !LB: ok, double-checked!
918      ELSE
919         !! Cd provided is the neutral-stability Cd, aka CdN :
920         z0_from_Cd = pzu * EXP( - vkarmn/SQRT(pCd(:,:)) )            !LB: ok, double-checked!
921      END IF
922   END FUNCTION z0_from_Cd
923
924   FUNCTION Cd_from_z0( pzu, pz0,  ppsi )
925      REAL(wp), DIMENSION(jpi,jpj) :: Cd_from_z0        !: (neutral or non-neutral) drag coefficient []
926      REAL(wp)                    , INTENT(in) :: pzu   !: reference height zu [m]
927      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0   !: roughness length [m]
928      REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
929      !!
930      !! If we want to return the NEUTRAL-STABILITY drag coefficient then ppsi must be 0 or not given
931      !! If we want to return the stability-corrected Cd (i.e. in stable or unstable conditions) then pssi must be provided
932      !!----------------------------------------------------------------------------------
933      IF( PRESENT(ppsi) ) THEN
934         !! The Cd we return is the actual Cd (not the neutral-stability CdN) :
935         Cd_from_z0 = 1._wp / ( LOG( pzu / pz0(:,:) ) - ppsi(:,:) )
936      ELSE
937         !! The Cd we return is the neutral-stability Cd, aka CdN :
938         Cd_from_z0 = 1._wp /   LOG( pzu / pz0(:,:) )
939      END IF
940      Cd_from_z0 = vkarmn2 * Cd_from_z0 * Cd_from_z0
941   END FUNCTION Cd_from_z0
942
943
944   FUNCTION f_m_louis_sclr( pzu, pRib, pCdn, pz0 )
945      !!----------------------------------------------------------------------------------
946      !!  Stability correction function for MOMENTUM
947      !!                 Louis (1979)
948      !!----------------------------------------------------------------------------------
949      REAL(wp)             :: f_m_louis_sclr ! term "f_m" in Eq.(6) when option "Louis" rather than "Psi(zeta) is chosen, Lupkes & Gryanik (2015),
950      REAL(wp), INTENT(in) :: pzu     ! reference height (height for pwnd)  [m]
951      REAL(wp), INTENT(in) :: pRib    ! Bulk Richardson number
952      REAL(wp), INTENT(in) :: pCdn    ! neutral drag coefficient
953      REAL(wp), INTENT(in) :: pz0     ! roughness length                      [m]
954      !!----------------------------------------------------------------------------------
955      REAL(wp) :: ztu, zts, zstab
956      !!----------------------------------------------------------------------------------
957      zstab = 0.5 + SIGN(0.5_wp, pRib) ; ! Unstable (Ri<0) => zstab = 0 | Stable (Ri>0) => zstab = 1
958      !
959      ztu = pRib / ( 1._wp + 3._wp * rc2_louis * pCdn * SQRT( ABS( -pRib * ( pzu / pz0 + 1._wp) ) ) ) ! ABS is just here for when it's stable conditions and ztu is not used anyways
960      zts = pRib / SQRT( ABS( 1._wp + pRib ) ) ! ABS is just here for when it's UNstable conditions and zts is not used anyways
961      !
962      f_m_louis_sclr = (1._wp - zstab) *         ( 1._wp - ram_louis * ztu )  &  ! Unstable Eq.(A6)
963         &               +      zstab  * 1._wp / ( 1._wp + ram_louis * zts )     ! Stable   Eq.(A7)
964      !
965   END FUNCTION f_m_louis_sclr
966   !!
967   FUNCTION f_m_louis_vctr( pzu, pRib, pCdn, pz0 )
968      REAL(wp), DIMENSION(jpi,jpj)             :: f_m_louis_vctr
969      REAL(wp),                     INTENT(in) :: pzu
970      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
971      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCdn
972      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
973      INTEGER  :: ji, jj
974      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
975      f_m_louis_vctr(ji,jj) = f_m_louis_sclr( pzu, pRib(ji,jj), pCdn(ji,jj), pz0(ji,jj) )
976      END_2D
977   END FUNCTION f_m_louis_vctr
978
979
980   FUNCTION f_h_louis_sclr( pzu, pRib, pChn, pz0 )
981      !!----------------------------------------------------------------------------------
982      !!  Stability correction function for HEAT
983      !!                 Louis (1979)
984      !!----------------------------------------------------------------------------------
985      REAL(wp)             :: f_h_louis_sclr ! term "f_h" in Eq.(6) when option "Louis" rather than "Psi(zeta) is chosen, Lupkes & Gryanik (2015),
986      REAL(wp), INTENT(in) :: pzu     ! reference height (height for pwnd)  [m]
987      REAL(wp), INTENT(in) :: pRib    ! Bulk Richardson number
988      REAL(wp), INTENT(in) :: pChn    ! neutral heat transfer coefficient
989      REAL(wp), INTENT(in) :: pz0     ! roughness length                      [m]
990      !!----------------------------------------------------------------------------------
991      REAL(wp) :: ztu, zts, zstab
992      !!----------------------------------------------------------------------------------
993      zstab = 0.5 + SIGN(0.5_wp, pRib) ; ! Unstable (Ri<0) => zstab = 0 | Stable (Ri>0) => zstab = 1
994      !
995      ztu = pRib / ( 1._wp + 3._wp * rc2_louis * pChn * SQRT( ABS(-pRib * ( pzu / pz0 + 1._wp) ) ) )
996      zts = pRib / SQRT( ABS( 1._wp + pRib ) )
997      !
998      f_h_louis_sclr = (1._wp - zstab) *         ( 1._wp - rah_louis * ztu )  &  ! Unstable Eq.(A6)
999         &              +       zstab  * 1._wp / ( 1._wp + rah_louis * zts )     ! Stable   Eq.(A7)  !#LB: in paper it's "ram_louis" and not "rah_louis" typo or what????
1000      !
1001   END FUNCTION f_h_louis_sclr
1002   !!
1003   FUNCTION f_h_louis_vctr( pzu, pRib, pChn, pz0 )
1004      REAL(wp), DIMENSION(jpi,jpj)             :: f_h_louis_vctr
1005      REAL(wp),                     INTENT(in) :: pzu
1006      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRib
1007      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pChn
1008      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0
1009      INTEGER  :: ji, jj
1010      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1011      f_h_louis_vctr(ji,jj) = f_h_louis_sclr( pzu, pRib(ji,jj), pChn(ji,jj), pz0(ji,jj) )
1012      END_2D
1013   END FUNCTION f_h_louis_vctr
1014
1015   FUNCTION UN10_from_ustar( pzu, pUzu, pus, ppsi )
1016      !!----------------------------------------------------------------------------------
1017      !!  Provides the neutral-stability wind speed at 10 m
1018      !!----------------------------------------------------------------------------------
1019      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_ustar  !: neutral stability wind speed at 10m [m/s]
1020      REAL(wp),                     INTENT(in) :: pzu   !: measurement heigh of wind speed   [m]
1021      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUzu  !: bulk wind speed at height pzu m   [m/s]
1022      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus   !: friction velocity                 [m/s]
1023      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
1024      !!----------------------------------------------------------------------------------
1025      UN10_from_ustar(:,:) = pUzu(:,:) - pus(:,:)/vkarmn * ( LOG(pzu/10._wp) - ppsi(:,:) )
1026      !!
1027   END FUNCTION UN10_from_ustar
1028
1029
1030   FUNCTION UN10_from_CD( pzu, pUb, pCd, ppsi )
1031      !!----------------------------------------------------------------------------------
1032      !!  Provides the neutral-stability wind speed at 10 m
1033      !!----------------------------------------------------------------------------------
1034      REAL(wp), DIMENSION(jpi,jpj)             :: UN10_from_CD  !: [m/s]
1035      REAL(wp),                     INTENT(in) :: pzu  !: measurement heigh of bulk wind speed
1036      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb  !: bulk wind speed at height pzu m   [m/s]
1037      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd  !: drag coefficient
1038      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppsi !: "Psi_m(pzu/L)" stability correction profile for momentum []
1039      !!----------------------------------------------------------------------------------
1040      !! Reminder: UN10 = u*/vkarmn * log(10/z0)
1041      !!     and: u* = sqrt(Cd) * Ub
1042      !!                                  u*/vkarmn * log(   10   /       z0    )
1043      UN10_from_CD(:,:) = SQRT(pCd(:,:))*pUb/vkarmn * LOG( 10._wp / z0_from_Cd( pzu, pCd(:,:), ppsi=ppsi(:,:) ) )
1044      !!
1045   END FUNCTION UN10_from_CD
1046
1047
1048   FUNCTION z0tq_LKB( iflag, pRer, pz0 )
1049      !!---------------------------------------------------------------------------------
1050      !!       ***  FUNCTION z0tq_LKB  ***
1051      !!
1052      !! ** Purpose : returns the "temperature/humidity roughness lengths"
1053      !!              * iflag==1 => temperature => returns: z_{0t}
1054      !!              * iflag==2 => humidity    => returns: z_{0q}
1055      !!              from roughness reynold number "pRer" (i.e. [z_0 u*]/Nu_{air})
1056      !!              between 0 and 1000. Out of range "pRer" indicated by prt=-999.
1057      !!              and roughness length (for momentum)
1058      !!
1059      !!              Based on Liu et al. (1979) JAS 36 1722-1723s
1060      !!
1061      !!              Note: this is what is used into COARE 2.5 to estimate z_{0t} and z_{0q}
1062      !!
1063      !! ** Author: L. Brodeau, April 2020 / AeroBulk (https://github.com/brodeau/aerobulk/)
1064      !!----------------------------------------------------------------------------------
1065      REAL(wp), DIMENSION(jpi,jpj)             :: z0tq_LKB
1066      INTEGER,                      INTENT(in) :: iflag     !: 1 => dealing with temperature; 2 => dealing with humidity
1067      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pRer      !: roughness Reynolds number  [z_0 u*]/Nu_{air}
1068      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pz0       !: roughness length (for momentum) [m]
1069      !-------------------------------------------------------------------
1070      ! Scalar Re_r relation from Liu et al.
1071      REAL(wp), DIMENSION(8,2), PARAMETER :: &
1072         & XA = (/ 0.177, 1.376, 1.026, 1.625, 4.661, 34.904, 1667.19, 5.88e5,  &
1073         &         0.292, 1.808, 1.393, 1.956, 4.994, 30.709, 1448.68, 2.98e5 /)
1074      !!
1075      REAL(wp), DIMENSION(8,2), PARAMETER :: &
1076         & XB = (/ 0., 0.929, -0.599, -1.018, -1.475, -2.067, -2.907, -3.935,  &
1077         &         0., 0.826, -0.528, -0.870, -1.297, -1.845, -2.682, -3.616 /)
1078      !!
1079      REAL(wp), DIMENSION(0:8),   PARAMETER :: &
1080         & XRAN = (/ 0., 0.11, 0.825, 3.0, 10.0, 30.0, 100., 300., 1000. /)
1081      !-------------------------------------------------------------------
1082      !
1083      !-------------------------------------------------------------------
1084      ! Scalar Re_r relation from Moana Wave data.
1085      !
1086      !      real*8 A(9,2),B(9,2),RAN(9),pRer,prt
1087      !      integer iflag
1088      !      DATA A/0.177,2.7e3,1.03,1.026,1.625,4.661,34.904,1667.19,5.88E5,
1089      !     &       0.292,3.7e3,1.4,1.393,1.956,4.994,30.709,1448.68,2.98E5/
1090      !      DATA B/0.,4.28,0,-0.599,-1.018,-1.475,-2.067,-2.907,-3.935,
1091      !     &       0.,4.28,0,-0.528,-0.870,-1.297,-1.845,-2.682,-3.616/
1092      !      DATA RAN/0.11,.16,1.00,3.0,10.0,30.0,100.,300.,1000./
1093      !-------------------------------------------------------------------
1094
1095      LOGICAL  :: lfound=.FALSE.
1096      REAL(wp) :: zrr
1097      INTEGER  :: ji, jj, jm
1098
1099      z0tq_LKB(:,:) = -999._wp
1100
1101      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
1102
1103      zrr    = pRer(ji,jj)
1104      lfound = .FALSE.
1105
1106      IF( (zrr > 0.).AND.(zrr < 1000.) ) THEN
1107         jm = 0
1108         DO WHILE ( .NOT. lfound )
1109            jm = jm + 1
1110            lfound = ( (zrr > XRAN(jm-1)) .AND. (zrr <= XRAN(jm)) )
1111         END DO
1112
1113         z0tq_LKB(ji,jj) = XA(jm,iflag)*zrr**XB(jm,iflag) * pz0(ji,jj)/zrr
1114
1115      END IF
1116
1117      END_2D
1118
1119      z0tq_LKB(:,:) = MIN( MAX(ABS(z0tq_LKB(:,:)), 1.E-9) , 0.05_wp )
1120
1121   END FUNCTION z0tq_LKB
1122
1123
1124   !!======================================================================
1125END MODULE sbc_phy
Note: See TracBrowser for help on using the repository browser.