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 @ 14036

Last change on this file since 14036 was 14036, checked in by laurent, 3 years ago

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