Newer
Older
MODULE sbc_phy
!!======================================================================
!! *** MODULE sbc_phy ***
!! A set of functions to compute air themodynamics parameters
!! needed by Aerodynamic Bulk Formulas
!!=====================================================================
!! 4.x ! 2020 L. Brodeau from AeroBulk package (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------
!! virt_temp : virtual (aka sensible) temperature (potential or absolute)
!! rho_air : density of (moist) air (depends on T_air, q_air and SLP
!! visc_air : kinematic viscosity (aka Nu_air) of air from temperature
!! L_vap : latent heat of vaporization of water as a function of temperature
!! cp_air : specific heat of (moist) air (depends spec. hum. q_air)
!! gamma_moist : adiabatic lapse-rate of moist air
!! One_on_L : 1. / ( Obukhov length )
!! Ri_bulk : bulk Richardson number aka BRN
!! q_sat : saturation humidity as a function of SLP and temperature
!! q_air_rh : specific humidity as a function of RH (fraction, not %), t_air and SLP
USE dom_oce ! ocean space and time domain
USE phycst ! physical constants
IMPLICIT NONE
Guillaume Samson
committed
PUBLIC !! Haleluja that was the solution for AGRIF
INTEGER , PARAMETER, PUBLIC :: nb_iter0 = 5 ! Default number of itterations in bulk-param algorithms (can be overriden b.m.o `nb_iter` optional argument)
!! (mainly removed from sbcblk.F90)
Guillaume Samson
committed
REAL(wp), PARAMETER, PUBLIC :: rCp_dry = 1005.0_wp !: Specic heat of dry air, constant pressure [J/K/kg]
REAL(wp), PARAMETER, PUBLIC :: rCp_vap = 1860.0_wp !: Specic heat of water vapor, constant pressure [J/K/kg]
REAL(wp), PARAMETER, PUBLIC :: R_dry = 287.05_wp !: Specific gas constant for dry air [J/K/kg]
REAL(wp), PARAMETER, PUBLIC :: R_vap = 461.495_wp !: Specific gas constant for water vapor [J/K/kg]
REAL(wp), PARAMETER, PUBLIC :: reps0 = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622
REAL(wp), PARAMETER, PUBLIC :: rctv0 = R_vap/R_dry - 1._wp !: for virtual temperature (== (1-eps)/eps) => ~ 0.608
REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp !: specific heat of air (only used for ice fluxes now...)
REAL(wp), PARAMETER, PUBLIC :: albo = 0.066_wp !: ocean albedo assumed to be constant
REAL(wp), PARAMETER, PUBLIC :: R_gas = 8.314510_wp !: Universal molar gas constant [J/mol/K]
REAL(wp), PARAMETER, PUBLIC :: rmm_dryair = 28.9647e-3_wp !: dry air molar mass / molecular weight [kg/mol]
REAL(wp), PARAMETER, PUBLIC :: rmm_water = 18.0153e-3_wp !: water molar mass / molecular weight [kg/mol]
REAL(wp), PARAMETER, PUBLIC :: rmm_ratio = rmm_water / rmm_dryair
REAL(wp), PARAMETER, PUBLIC :: rgamma_dry = R_gas / ( rmm_dryair * rCp_dry ) !: Poisson constant for dry air
REAL(wp), PARAMETER, PUBLIC :: rpref = 1.e5_wp !: reference air pressure for exner function [Pa]
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
!
REAL(wp), PARAMETER, PUBLIC :: rho0_a = 1.2_wp !: Approx. of density of air [kg/m^3]
REAL(wp), PARAMETER, PUBLIC :: rho0_w = 1025._wp !: Density of sea-water (ECMWF->1025) [kg/m^3]
REAL(wp), PARAMETER, PUBLIC :: radrw = rho0_a/rho0_w !: Density ratio
REAL(wp), PARAMETER, PUBLIC :: sq_radrw = SQRT(rho0_a/rho0_w)
REAL(wp), PARAMETER, PUBLIC :: rCp0_w = 4190._wp !: Specific heat capacity of seawater (ECMWF 4190) [J/K/kg]
REAL(wp), PARAMETER, PUBLIC :: rnu0_w = 1.e-6_wp !: kinetic viscosity of water [m^2/s]
REAL(wp), PARAMETER, PUBLIC :: rk0_w = 0.6_wp !: thermal conductivity of water (at 20C) [W/m/K]
!
REAL(wp), PARAMETER, PUBLIC :: emiss_w = 0.98_wp !: Long-wave (thermal) emissivity of sea-water []
!
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 ...
REAL(wp), PARAMETER, PUBLIC :: wspd_thrshld_ice = 0.2_wp !: minimum scalar wind speed accepted over sea-ice... [m/s]
!
REAL(wp), PARAMETER, PUBLIC :: rdct_qsat_salt = 0.98_wp !: reduction factor on specific humidity at saturation (q_sat(T_s)) due to salt
REAL(wp), PARAMETER, PUBLIC :: rtt0 = 273.16_wp !: triple point of temperature [K]
!
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)
! => see eq.(14) in Fairall et al. 1996 (eq.(6) of Zeng aand Beljaars is WRONG! (typo?)
REAL(wp), PARAMETER, PUBLIC :: z0_sea_max = 0.0025_wp !: maximum realistic value for roughness length of sea-surface... [m]
REAL(wp), PUBLIC, SAVE :: pp_cldf = 0.81 !: cloud fraction over sea ice, summer CLIO value [-]
REAL(wp), PARAMETER, PUBLIC :: Cx_min = 0.1E-3_wp ! smallest value allowed for bulk transfer coefficients (usually in stable conditions with now wind)
REAL(wp), PARAMETER :: &
!! Constants for Goff formula in the presence of ice:
& rAg_i = -9.09718_wp, &
& rBg_i = -3.56654_wp, &
& rCg_i = 0.876793_wp, &
& rDg_i = LOG10(6.1071_wp)
REAL(wp), PARAMETER :: rc_louis = 5._wp
REAL(wp), PARAMETER :: rc2_louis = rc_louis * rc_louis
REAL(wp), PARAMETER :: ram_louis = 2. * rc_louis
REAL(wp), PARAMETER :: rah_louis = 3. * rc_louis
INTERFACE virt_temp
MODULE PROCEDURE virt_temp_vctr, virt_temp_sclr
END INTERFACE virt_temp
Guillaume Samson
committed
INTERFACE pres_temp
MODULE PROCEDURE pres_temp_vctr, pres_temp_sclr
END INTERFACE pres_temp
INTERFACE theta_exner
MODULE PROCEDURE theta_exner_vctr, theta_exner_sclr
END INTERFACE theta_exner
INTERFACE visc_air
MODULE PROCEDURE visc_air_vctr, visc_air_sclr
END INTERFACE visc_air
INTERFACE gamma_moist
MODULE PROCEDURE gamma_moist_vctr, gamma_moist_sclr
END INTERFACE gamma_moist
INTERFACE e_sat
MODULE PROCEDURE e_sat_vctr, e_sat_sclr
END INTERFACE e_sat
INTERFACE e_sat_ice
MODULE PROCEDURE e_sat_ice_vctr, e_sat_ice_sclr
END INTERFACE e_sat_ice
Guillaume Samson
committed
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
INTERFACE de_sat_dt_ice
MODULE PROCEDURE de_sat_dt_ice_vctr, de_sat_dt_ice_sclr
END INTERFACE de_sat_dt_ice
INTERFACE Ri_bulk
MODULE PROCEDURE Ri_bulk_vctr, Ri_bulk_sclr
END INTERFACE Ri_bulk
INTERFACE q_sat
MODULE PROCEDURE q_sat_vctr, q_sat_sclr
END INTERFACE q_sat
INTERFACE dq_sat_dt_ice
MODULE PROCEDURE dq_sat_dt_ice_vctr, dq_sat_dt_ice_sclr
END INTERFACE dq_sat_dt_ice
INTERFACE L_vap
MODULE PROCEDURE L_vap_vctr, L_vap_sclr
END INTERFACE L_vap
INTERFACE rho_air
MODULE PROCEDURE rho_air_vctr, rho_air_sclr
END INTERFACE rho_air
INTERFACE cp_air
MODULE PROCEDURE cp_air_vctr, cp_air_sclr
END INTERFACE cp_air
INTERFACE alpha_sw
MODULE PROCEDURE alpha_sw_vctr, alpha_sw_sclr
END INTERFACE alpha_sw
INTERFACE bulk_formula
MODULE PROCEDURE bulk_formula_vctr, bulk_formula_sclr
END INTERFACE bulk_formula
INTERFACE qlw_net
MODULE PROCEDURE qlw_net_vctr, qlw_net_sclr
END INTERFACE qlw_net
INTERFACE f_m_louis
MODULE PROCEDURE f_m_louis_vctr, f_m_louis_sclr
END INTERFACE f_m_louis
INTERFACE f_h_louis
MODULE PROCEDURE f_h_louis_vctr, f_h_louis_sclr
END INTERFACE f_h_louis
PUBLIC virt_temp
Guillaume Samson
committed
PUBLIC pres_temp
PUBLIC theta_exner
PUBLIC rho_air
PUBLIC visc_air
PUBLIC L_vap
PUBLIC cp_air
PUBLIC gamma_moist
PUBLIC One_on_L
PUBLIC Ri_bulk
PUBLIC q_sat
PUBLIC q_air_rh
PUBLIC dq_sat_dt_ice
Guillaume Samson
committed
!
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
PUBLIC update_qnsol_tau
PUBLIC alpha_sw
PUBLIC bulk_formula
PUBLIC qlw_net
!
PUBLIC f_m_louis, f_h_louis
PUBLIC z0_from_Cd
PUBLIC Cd_from_z0
PUBLIC UN10_from_ustar
PUBLIC UN10_from_CD
PUBLIC z0tq_LKB
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sbcblk.F90 10535 2019-01-16 17:36:47Z clem $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
FUNCTION virt_temp_sclr( pta, pqa )
!!------------------------------------------------------------------------
!!
!! Compute the (absolute/potential) VIRTUAL temperature, based on the
!! (absolute/potential) temperature and specific humidity
!!
!! If input temperature is absolute then output virtual temperature is absolute
!! If input temperature is potential then output virtual temperature is potential
!!
!! Author: L. Brodeau, June 2019 / AeroBulk
!! (https://github.com/brodeau/aerobulk/)
!!------------------------------------------------------------------------
REAL(wp) :: virt_temp_sclr !: virtual temperature [K]
REAL(wp), INTENT(in) :: pta !: absolute or potential air temperature [K]
REAL(wp), INTENT(in) :: pqa !: specific humidity of air [kg/kg]
!!-------------------------------------------------------------------
!
virt_temp_sclr = pta * (1._wp + rctv0*pqa)
!!
!! This is exactly the same thing as:
!! virt_temp_sclr = pta * ( pwa + reps0) / (reps0*(1.+pwa))
!! with wpa (mixing ration) defined as : pwa = pqa/(1.-pqa)
!
END FUNCTION virt_temp_sclr
Guillaume Samson
committed
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: virt_temp_vctr !: virtual temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute or potential air temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: specific humidity of air [kg/kg]
Guillaume Samson
committed
virt_temp_vctr(:,:) = pta(:,:) * (1._wp + rctv0*pqa(:,:))
Guillaume Samson
committed
Guillaume Samson
committed
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
FUNCTION pres_temp_sclr( pqspe, pslp, pz, ptpot, pta, l_ice )
!!-------------------------------------------------------------------------------
!! *** FUNCTION pres_temp ***
!!
!! ** Purpose : compute air pressure using barometric equation
!! from either potential or absolute air temperature
!! ** Author: G. Samson, Feb 2021
!!-------------------------------------------------------------------------------
REAL(wp) :: pres_temp_sclr ! air pressure [Pa]
REAL(wp), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
REAL(wp), INTENT(in ) :: pslp ! sea-level pressure [Pa]
REAL(wp), INTENT(in ) :: pz ! height above surface [m]
REAL(wp), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
REAL(wp), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
REAL(wp) :: ztpot, zta, zpa, zxm, zmask, zqsat
INTEGER :: it, niter = 3 ! iteration indice and number
LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
LOGICAL :: lice ! sea-ice presence
IF( PRESENT(ptpot) ) THEN
zmask = 1._wp
ztpot = ptpot
zta = 0._wp
ELSE
zmask = 0._wp
ztpot = 0._wp
zta = pta
ENDIF
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
zpa = pslp ! air pressure first guess [Pa]
DO it = 1, niter
zta = ztpot * ( zpa / rpref )**rgamma_dry * zmask + (1._wp - zmask) * zta
zqsat = q_sat( zta, zpa, l_ice=lice ) ! saturation specific humidity [kg/kg]
zxm = (1._wp - pqspe/zqsat) * rmm_dryair + pqspe/zqsat * rmm_water ! moist air molar mass [kg/mol]
zpa = pslp * EXP( -grav * zxm * pz / ( R_gas * zta ) )
END DO
pres_temp_sclr = zpa
IF(( PRESENT(pta) ).AND.( PRESENT(ptpot) )) pta = zta
END FUNCTION pres_temp_sclr
FUNCTION pres_temp_vctr( pqspe, pslp, pz, ptpot, pta, l_ice )
!!-------------------------------------------------------------------------------
!! *** FUNCTION pres_temp ***
!!
!! ** Purpose : compute air pressure using barometric equation
!! from either potential or absolute air temperature
!! ** Author: G. Samson, Feb 2021
!!-------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: pres_temp_vctr ! air pressure [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pqspe ! air specific humidity [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pslp ! sea-level pressure [Pa]
REAL(wp), INTENT(in ) :: pz ! height above surface [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) , OPTIONAL :: ptpot ! air potential temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout), OPTIONAL :: pta ! air absolute temperature [K]
INTEGER :: ji, jj ! loop indices
LOGICAL , INTENT(in) , OPTIONAL :: l_ice ! sea-ice presence
LOGICAL :: lice ! sea-ice presence
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
IF( PRESENT(ptpot) ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, ptpot=ptpot(ji,jj), pta=pta(ji,jj), l_ice=lice )
END_2D
ELSE
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
pres_temp_vctr(ji,jj) = pres_temp_sclr( pqspe(ji,jj), pslp(ji,jj), pz, pta=pta(ji,jj), l_ice=lice )
END_2D
ENDIF
END FUNCTION pres_temp_vctr
FUNCTION theta_exner_sclr( pta, ppa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION theta_exner ***
!!
!! ** Purpose : compute air/surface potential temperature from absolute temperature
!! and pressure using Exner function
!! ** Author: G. Samson, Feb 2021
!!-------------------------------------------------------------------------------
REAL(wp) :: theta_exner_sclr ! air/surface potential temperature [K]
REAL(wp), INTENT(in) :: pta ! air/surface absolute temperature [K]
REAL(wp), INTENT(in) :: ppa ! air/surface pressure [Pa]
theta_exner_sclr = pta * ( rpref / ppa ) ** rgamma_dry
END FUNCTION theta_exner_sclr
FUNCTION theta_exner_vctr( pta, ppa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION theta_exner ***
!!
!! ** Purpose : compute air/surface potential temperature from absolute temperature
!! and pressure using Exner function
!! ** Author: G. Samson, Feb 2021
!!-------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: theta_exner_vctr ! air/surface potential temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta ! air/surface absolute temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! air/surface pressure [Pa]
INTEGER :: ji, jj ! loop indices
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
theta_exner_vctr(ji,jj) = theta_exner_sclr( pta(ji,jj), ppa(ji,jj) )
END_2D
END FUNCTION theta_exner_vctr
FUNCTION rho_air_vctr( ptak, pqa, ppa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION rho_air_vctr ***
!!
!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
!!
!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! pressure in [Pa]
REAL(wp), DIMENSION(jpi,jpj) :: rho_air_vctr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
Guillaume Samson
committed
rho_air_vctr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
Guillaume Samson
committed
END FUNCTION rho_air_vctr
FUNCTION rho_air_sclr( ptak, pqa, ppa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION rho_air_sclr ***
!!
!! ** Purpose : compute density of (moist) air using the eq. of state of the atmosphere
!!
!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
REAL(wp), INTENT(in) :: ptak ! air temperature [K]
REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp), INTENT(in) :: ppa ! pressure in [Pa]
REAL(wp) :: rho_air_sclr ! density of moist air [kg/m^3]
!!-------------------------------------------------------------------------------
rho_air_sclr = MAX( ppa / (R_dry*ptak * ( 1._wp + rctv0*pqa )) , 0.8_wp )
Guillaume Samson
committed
END FUNCTION rho_air_sclr
FUNCTION visc_air_sclr(ptak)
!!----------------------------------------------------------------------------------
!! Air kinetic viscosity (m^2/s) given from air temperature in Kelvin
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: visc_air_sclr ! kinetic viscosity (m^2/s)
REAL(wp), INTENT(in) :: ptak ! air temperature in (K)
!
REAL(wp) :: ztc, ztc2 ! local scalar
!!----------------------------------------------------------------------------------
!
ztc = ptak - rt0 ! air temp, in deg. C
ztc2 = ztc*ztc
visc_air_sclr = 1.326e-5*(1. + 6.542E-3*ztc + 8.301e-6*ztc2 - 4.84e-9*ztc2*ztc)
!
END FUNCTION visc_air_sclr
FUNCTION visc_air_vctr(ptak)
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: visc_air_vctr ! kinetic viscosity (m^2/s)
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak ! air temperature in (K)
INTEGER :: ji, jj ! dummy loop indices
Guillaume Samson
committed
Guillaume Samson
committed
visc_air_vctr(ji,jj) = visc_air_sclr( ptak(ji,jj) )
Guillaume Samson
committed
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
END FUNCTION visc_air_vctr
FUNCTION L_vap_vctr( psst )
!!---------------------------------------------------------------------------------
!! *** FUNCTION L_vap_vctr ***
!!
!! ** Purpose : Compute the latent heat of vaporization of water from temperature
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: L_vap_vctr ! latent heat of vaporization [J/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
!!----------------------------------------------------------------------------------
!
L_vap_vctr = ( 2.501_wp - 0.00237_wp * ( psst(:,:) - rt0) ) * 1.e6_wp
!
END FUNCTION L_vap_vctr
FUNCTION L_vap_sclr( psst )
!!---------------------------------------------------------------------------------
!! *** FUNCTION L_vap_sclr ***
!!
!! ** Purpose : Compute the latent heat of vaporization of water from temperature
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: L_vap_sclr ! latent heat of vaporization [J/kg]
REAL(wp), INTENT(in) :: psst ! water temperature [K]
!!----------------------------------------------------------------------------------
!
L_vap_sclr = ( 2.501_wp - 0.00237_wp * ( psst - rt0) ) * 1.e6_wp
!
END FUNCTION L_vap_sclr
FUNCTION cp_air_vctr( pqa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION cp_air_vctr ***
!!
!! ** Purpose : Compute specific heat (Cp) of moist air
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp), DIMENSION(jpi,jpj) :: cp_air_vctr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
Guillaume Samson
committed
Guillaume Samson
committed
END FUNCTION cp_air_vctr
FUNCTION cp_air_sclr( pqa )
!!-------------------------------------------------------------------------------
!! *** FUNCTION cp_air_sclr ***
!!
!! ** Purpose : Compute specific heat (Cp) of moist air
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!-------------------------------------------------------------------------------
REAL(wp), INTENT(in) :: pqa ! air specific humidity [kg/kg]
REAL(wp) :: cp_air_sclr ! specific heat of moist air [J/K/kg]
!!-------------------------------------------------------------------------------
Guillaume Samson
committed
Guillaume Samson
committed
END FUNCTION cp_air_sclr
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
FUNCTION gamma_moist_sclr( ptak, pqa )
!!----------------------------------------------------------------------------------
!! ** Purpose : Compute the moist adiabatic lapse-rate.
!! => http://glossary.ametsoc.org/wiki/Moist-adiabatic_lapse_rate
!! => http://www.geog.ucsb.edu/~joel/g266_s10/lecture_notes/chapt03/oh10_3_01/oh10_3_01.html
!!
!! ** Author: L. Brodeau, June 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: gamma_moist_sclr ! [K/m]
REAL(wp), INTENT(in) :: ptak ! absolute air temperature [K] !#LB: double check it's absolute !!!
REAL(wp), INTENT(in) :: pqa ! specific humidity [kg/kg]
!
REAL(wp) :: zta, zqa, zwa, ziRT, zLvap ! local scalars
!!----------------------------------------------------------------------------------
zta = MAX( ptak, 180._wp) ! prevents screw-up over masked regions where field == 0.
zqa = MAX( pqa, 1.E-6_wp) ! " " "
!!
zwa = zqa / (1._wp - zqa) ! w is mixing ratio w = q/(1-q) | q = w/(1+w)
ziRT = 1._wp / (R_dry*zta) ! 1/RT
zLvap = L_vap_sclr( ptak )
!!
gamma_moist_sclr = grav * ( 1._wp + zLvap*zwa*ziRT ) / ( rCp_dry + zLvap*zLvap*zwa*reps0*ziRT/zta )
!!
END FUNCTION gamma_moist_sclr
Guillaume Samson
committed
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: gamma_moist_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa
INTEGER :: ji, jj
Guillaume Samson
committed
Guillaume Samson
committed
gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) )
Guillaume Samson
committed
END FUNCTION gamma_moist_vctr
FUNCTION One_on_L( ptha, pqa, pus, pts, pqs )
!!------------------------------------------------------------------------
!!
!! Evaluates the 1./(Obukhov length) from air temperature,
!! air specific humidity, and frictional scales u*, t* and q*
!!
!! Author: L. Brodeau, June 2019 / AeroBulk
!! (https://github.com/brodeau/aerobulk/)
!!------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: One_on_L !: 1./(Obukhov length) [m^-1]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha !: reference potential temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa !: reference specific humidity of air [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pus !: u*: friction velocity [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pts, pqs !: \theta* and q* friction aka turb. scales for temp. and spec. hum.
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: zqa ! local scalar
!!-------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
zqa = (1._wp + rctv0*pqa(ji,jj))
!
! The main concern is to know whether, the vertical turbulent flux of virtual temperature, < u' theta_v' > is estimated with:
! a/ -u* [ theta* (1 + 0.61 q) + 0.61 theta q* ] => this is the one that seems correct! chose this one!
! or
! b/ -u* [ theta* + 0.61 theta q* ]
!
One_on_L(ji,jj) = grav*vkarmn*( pts(ji,jj)*zqa + rctv0*ptha(ji,jj)*pqs(ji,jj) ) &
& / MAX( pus(ji,jj)*pus(ji,jj)*ptha(ji,jj)*zqa , 1.E-9_wp )
END_2D
!
One_on_L = SIGN( MIN(ABS(One_on_L),200._wp), One_on_L ) ! (prevent FPE from stupid values over masked regions...)
!
END FUNCTION One_on_L
FUNCTION Ri_bulk_sclr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
!!----------------------------------------------------------------------------------
!! Bulk Richardson number according to "wide-spread equation"...
!!
!! Reminder: the Richardson number is the ratio "buoyancy" / "shear"
!!
!! ** Author: L. Brodeau, June 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: Ri_bulk_sclr
REAL(wp), INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
Guillaume Samson
committed
REAL(wp), INTENT(in) :: psst ! potential SST [K]
REAL(wp), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K]
REAL(wp), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg]
REAL(wp), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg]
REAL(wp), INTENT(in) :: pub ! bulk wind speed [m/s]
REAL(wp), INTENT(in), OPTIONAL :: pta_layer ! when possible, a better guess of absolute temperature WITHIN the layer [K]
REAL(wp), INTENT(in), OPTIONAL :: pqa_layer ! when possible, a better guess of specific humidity WITHIN the layer [kg/kg]
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
REAL(wp) :: zqa, zta, zgamma, zdthv, ztv, zsstv ! local scalars
Guillaume Samson
committed
REAL(wp) :: ztptv
!!-------------------------------------------------------------------
Guillaume Samson
committed
IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
Guillaume Samson
committed
zsstv = virt_temp_sclr( psst, pssq ) ! virtual potential SST
ztptv = virt_temp_sclr( ptha, pqa ) ! virtual potential air temperature
zdthv = ztptv - zsstv ! air-sea delta of "virtual potential temperature"
Guillaume Samson
committed
Ri_bulk_sclr = grav * zdthv * pz / ( ztptv * pub * pub ) ! the usual definition of Ri_bulk_sclr
Guillaume Samson
committed
FUNCTION Ri_bulk_vctr( pz, psst, ptha, pssq, pqa, pub, pta_layer, pqa_layer )
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: Ri_bulk_vctr
REAL(wp) , INTENT(in) :: pz ! height above the sea (aka "delta z") [m]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! SST [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptha ! pot. air temp. at height "pz" [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssq ! 0.98*q_sat(SST) [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! air spec. hum. at height "pz" [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pub ! bulk wind speed [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: pta_layer ! when possible, a better guess of absolute temperature WITHIN the layer [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in), OPTIONAL :: pqa_layer ! when possible, a better guess of specific humidity WITHIN the layer [kg/kg]
!!
LOGICAL :: l_ptqa_l_prvd = .FALSE.
INTEGER :: ji, jj
Guillaume Samson
committed
IF( PRESENT(pta_layer) .AND. PRESENT(pqa_layer) ) l_ptqa_l_prvd = .TRUE.
Guillaume Samson
committed
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj), &
& pta_layer=pta_layer(ji,jj ), pqa_layer=pqa_layer(ji,jj ) )
END_2D
Guillaume Samson
committed
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Ri_bulk_vctr(ji,jj) = Ri_bulk_sclr( pz, psst(ji,jj), ptha(ji,jj), pssq(ji,jj), pqa(ji,jj), pub(ji,jj) )
END_2D
Guillaume Samson
committed
Guillaume Samson
committed
FUNCTION e_sat_sclr( ptak )
!!----------------------------------------------------------------------------------
!! *** FUNCTION e_sat_sclr ***
!! < SCALAR argument version >
!! ** Purpose : water vapor at saturation in [Pa]
!! Based on accurate estimate by Goff, 1957
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!
!! Note: what rt0 should be here, is 273.16 (triple point of water) and not 273.15 like here
!!----------------------------------------------------------------------------------
REAL(wp) :: e_sat_sclr ! water vapor at saturation [kg/kg]
REAL(wp), INTENT(in) :: ptak ! air temperature [K]
REAL(wp) :: zta, ztmp ! local scalar
!!----------------------------------------------------------------------------------
zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
ztmp = rt0 / zta !#LB: rt0 or rtt0 ???? (273.15 vs 273.16 )
!
! Vapour pressure at saturation [Pa] : WMO, (Goff, 1957)
e_sat_sclr = 100.*( 10.**( 10.79574*(1. - ztmp) - 5.028*LOG10(zta/rt0) &
& + 1.50475*10.**(-4)*(1. - 10.**(-8.2969*(zta/rt0 - 1.)) ) &
& + 0.42873*10.**(-3)*(10.**(4.76955*(1. - ztmp)) - 1.) + 0.78614) )
!
END FUNCTION e_sat_sclr
Guillaume Samson
committed
FUNCTION e_sat_vctr(ptak)
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_vctr !: vapour pressure at saturation [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: temperature (K)
INTEGER :: ji, jj ! dummy loop indices
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
e_sat_vctr(ji,jj) = e_sat_sclr(ptak(ji,jj))
END_2D
END FUNCTION e_sat_vctr
Guillaume Samson
committed
FUNCTION e_sat_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! Same as "e_sat" but over ice rather than water!
!!---------------------------------------------------------------------------------
REAL(wp) :: e_sat_ice_sclr !: vapour pressure at saturation in presence of ice [Pa]
REAL(wp), INTENT(in) :: ptak
!!
REAL(wp) :: zta, zle, ztmp
!!---------------------------------------------------------------------------------
zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
ztmp = rtt0/zta
!!
zle = rAg_i*(ztmp - 1._wp) + rBg_i*LOG10(ztmp) + rCg_i*(1._wp - zta/rtt0) + rDg_i
!!
e_sat_ice_sclr = 100._wp * 10._wp**zle
Guillaume Samson
committed
Guillaume Samson
committed
FUNCTION e_sat_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: e_sat_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
e_sat_ice_vctr(ji,jj) = e_sat_ice_sclr( ptak(ji,jj) )
Guillaume Samson
committed
Guillaume Samson
committed
FUNCTION de_sat_dt_ice_sclr(ptak)
!!---------------------------------------------------------------------------------
!! d [ e_sat_ice ] / dT (derivative / temperature)
!! Analytical exact formulation: double checked!!!
!! => DOUBLE-check possible / finite-difference version with "./bin/test_phymbl.x"
!!---------------------------------------------------------------------------------
REAL(wp) :: de_sat_dt_ice_sclr !: [Pa/K]
REAL(wp), INTENT(in) :: ptak
!!
REAL(wp) :: zta, zde
!!---------------------------------------------------------------------------------
zta = MAX( ptak , 180._wp ) ! air temp., prevents fpe0 errors dute to unrealistically low values over masked regions...
!!
zde = -(rAg_i*rtt0)/(zta*zta) - rBg_i/(zta*LOG(10._wp)) - rCg_i/rtt0
!!
de_sat_dt_ice_sclr = LOG(10._wp) * zde * e_sat_ice_sclr(zta)
END FUNCTION de_sat_dt_ice_sclr
Guillaume Samson
committed
FUNCTION de_sat_dt_ice_vctr(ptak)
!! Same as "e_sat" but over ice rather than water!
REAL(wp), DIMENSION(jpi,jpj) :: de_sat_dt_ice_vctr !: vapour pressure at saturation in presence of ice [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
de_sat_dt_ice_vctr(ji,jj) = de_sat_dt_ice_sclr( ptak(ji,jj) )
Guillaume Samson
committed
END FUNCTION de_sat_dt_ice_vctr
FUNCTION q_sat_sclr( pta, ppa, l_ice )
!!---------------------------------------------------------------------------------
!! *** FUNCTION q_sat_sclr ***
!!
!! ** Purpose : Conputes specific humidity of air at saturation
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: q_sat_sclr
REAL(wp), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), INTENT(in) :: ppa !: atmospheric pressure [Pa]
LOGICAL, INTENT(in), OPTIONAL :: l_ice !: we are above ice
REAL(wp) :: ze_s
LOGICAL :: lice
!!----------------------------------------------------------------------------------
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
IF( lice ) THEN
ze_s = e_sat_ice( pta )
ELSE
ze_s = e_sat( pta ) ! Vapour pressure at saturation (Goff) :
END IF
q_sat_sclr = reps0*ze_s/(ppa - (1._wp - reps0)*ze_s)
Guillaume Samson
committed
Guillaume Samson
committed
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: q_sat_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
LOGICAL, INTENT(in), OPTIONAL :: l_ice !: we are above ice
LOGICAL :: lice
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
lice = .FALSE.
IF( PRESENT(l_ice) ) lice = l_ice
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
q_sat_vctr(ji,jj) = q_sat_sclr( pta(ji,jj) , ppa(ji,jj), l_ice=lice )
Guillaume Samson
committed
END FUNCTION q_sat_vctr
FUNCTION dq_sat_dt_ice_sclr( pta, ppa )
!!---------------------------------------------------------------------------------
!! *** FUNCTION dq_sat_dt_ice_sclr ***
!! => d [ q_sat_ice(T) ] / dT
!! Analytical exact formulation: double checked!!!
!! => DOUBLE-check possible / finite-difference version with "./bin/test_phymbl.x"
!!----------------------------------------------------------------------------------
REAL(wp) :: dq_sat_dt_ice_sclr
REAL(wp), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), INTENT(in) :: ppa !: atmospheric pressure [Pa]
REAL(wp) :: ze_s, zde_s_dt, ztmp
!!----------------------------------------------------------------------------------
ze_s = e_sat_ice_sclr( pta ) ! Vapour pressure at saturation in presence of ice (Goff)
zde_s_dt = de_sat_dt_ice( pta )
!
ztmp = (reps0 - 1._wp)*ze_s + ppa
!
dq_sat_dt_ice_sclr = reps0*ppa*zde_s_dt / ( ztmp*ztmp )
!
END FUNCTION dq_sat_dt_ice_sclr
Guillaume Samson
committed
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj) :: dq_sat_dt_ice_vctr
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pta !: absolute temperature of air [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
dq_sat_dt_ice_vctr(ji,jj) = dq_sat_dt_ice_sclr( pta(ji,jj) , ppa(ji,jj) )
Guillaume Samson
committed
END FUNCTION dq_sat_dt_ice_vctr
FUNCTION q_air_rh(prha, ptak, ppa)
!!----------------------------------------------------------------------------------
!! Specific humidity of air out of Relative Humidity
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: q_air_rh
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prha !: relative humidity [fraction, not %!!!]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptak !: air temperature [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa !: atmospheric pressure [Pa]
!
INTEGER :: ji, jj ! dummy loop indices
REAL(wp) :: ze ! local scalar
!!----------------------------------------------------------------------------------
!
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj))
q_air_rh(ji,jj) = ze*reps0/(ppa(ji,jj) - (1. - reps0)*ze)
Guillaume Samson
committed
SUBROUTINE UPDATE_QNSOL_TAU( pzu, pTs, pqs, pTa, pqa, pust, ptst, pqst, pwnd, pUb, ppa, prlw, prhoa, &
& pQns, pTau, &
& Qlat)
!!----------------------------------------------------------------------------------
!! Purpose: returns the non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw"
!! and the module of the wind stress => pTau = Tau
!! ** Author: L. Brodeau, Sept. 2019 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pust ! u*
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ptst ! t*
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqst ! q*
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prlw ! downwelling longwave radiative flux [W/m^2]
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! air density [kg/m3]
!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQns ! non-solar heat flux to the ocean aka "Qlat + Qsen + Qlw" [W/m^2]]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
!
REAL(wp), DIMENSION(jpi,jpj), OPTIONAL, INTENT(out) :: Qlat
!
REAL(wp) :: zdt, zdq, zCd, zCh, zCe, zz0, zQlat, zQsen, zQlw
INTEGER :: ji, jj ! dummy loop indices
!!----------------------------------------------------------------------------------
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
zdt = pTa(ji,jj) - pTs(ji,jj) ; zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt )
zdq = pqa(ji,jj) - pqs(ji,jj) ; zdq = SIGN( MAX(ABS(zdq),1.E-9_wp), zdq )
zz0 = pust(ji,jj)/pUb(ji,jj)
zCd = zz0*zz0
zCh = zz0*ptst(ji,jj)/zdt
zCe = zz0*pqst(ji,jj)/zdq
Guillaume Samson
committed
CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), zCd, zCh, zCe, &
& pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
& pTau(ji,jj), zQsen, zQlat )
Guillaume Samson
committed
zQlw = qlw_net_sclr( prlw(ji,jj), pTs(ji,jj) ) ! Net longwave flux
pQns(ji,jj) = zQlat + zQsen + zQlw
IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat
Guillaume Samson
committed
END SUBROUTINE UPDATE_QNSOL_TAU
SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
Guillaume Samson
committed
& pwnd, pUb, ppa, prhoa, &
Guillaume Samson
committed
& pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
Guillaume Samson
committed
REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
REAL(wp), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
REAL(wp), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
REAL(wp), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
REAL(wp), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), INTENT(in) :: pCd
REAL(wp), INTENT(in) :: pCh
REAL(wp), INTENT(in) :: pCe
Guillaume Samson
committed
REAL(wp), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
REAL(wp), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
REAL(wp), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
REAL(wp), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
INTEGER :: jq
!!----------------------------------------------------------------------------------
zfact_evap = 1._wp
IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
Guillaume Samson
committed
zUrho = pUb*MAX(prhoa, 1._wp) ! rho*U10
pTau = zUrho * pCd * pwnd ! Wind stress module
zevap = zUrho * pCe * (pqa - pqs)
pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa)
pQlat = L_vap(pTs) * zevap
IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap
END SUBROUTINE BULK_FORMULA_SCLR
Guillaume Samson
committed
SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, &
& pCd, pCh, pCe, &
Guillaume Samson
committed
& pwnd, pUb, ppa, prhoa, &
Guillaume Samson
committed
& pEvap, pfact_evap )
!!----------------------------------------------------------------------------------
Guillaume Samson
committed
REAL(wp), INTENT(in) :: pzu ! height above the sea-level where all this takes place (normally 10m)
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTs ! water temperature at the air-sea interface [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqs ! satur. spec. hum. at T=pTs [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pTa ! potential air temperature at z=pzu [K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pqa ! specific humidity at z=pzu [kg/kg]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCd
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCh
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pCe
Guillaume Samson
committed
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pwnd ! wind speed module at z=pzu [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pUb ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: ppa ! sea-level atmospheric pressure [Pa]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: prhoa ! Air density at z=pzu [kg/m^3]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pTau ! module of the wind stress [N/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQsen ! [W/m^2]
REAL(wp), DIMENSION(jpi,jpj), INTENT(out) :: pQlat ! [W/m^2]
!!
REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s]
REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent)
!!
REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap
INTEGER :: ji, jj
!!----------------------------------------------------------------------------------
zfact_evap = 1._wp
IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
Guillaume Samson
committed
CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), &
& pCd(ji,jj), pCh(ji,jj), pCe(ji,jj), &
& pwnd(ji,jj), pUb(ji,jj), ppa(ji,jj), prhoa(ji,jj), &
& pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj), &
& pEvap=zevap, pfact_evap=zfact_evap )
Guillaume Samson
committed
IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap
Guillaume Samson
committed
END SUBROUTINE BULK_FORMULA_VCTR
FUNCTION alpha_sw_vctr( psst )
!!---------------------------------------------------------------------------------
!! *** FUNCTION alpha_sw_vctr ***
!!
!! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: alpha_sw_vctr ! thermal expansion coefficient of sea-water [1/K]
REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psst ! water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_vctr = 2.1e-5_wp * MAX(psst(:,:)-rt0 + 3.2_wp, 0._wp)**0.79
Guillaume Samson
committed
END FUNCTION alpha_sw_vctr
FUNCTION alpha_sw_sclr( psst )
!!---------------------------------------------------------------------------------
!! *** FUNCTION alpha_sw_sclr ***
!!
!! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa)
!!
!! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/)
!!----------------------------------------------------------------------------------
REAL(wp) :: alpha_sw_sclr ! thermal expansion coefficient of sea-water [1/K]
REAL(wp), INTENT(in) :: psst ! sea-water temperature [K]
!!----------------------------------------------------------------------------------
alpha_sw_sclr = 2.1e-5_wp * MAX(psst-rt0 + 3.2_wp, 0._wp)**0.79
Guillaume Samson
committed