Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

sbc_phy.F90 @ 15540

Last change on this file since 15540 was 15540, checked in by sparonuz, 3 years ago
Mixed precision version, tested up to 30 years on ORCA2.
File size: 58.7 KB

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

Note: See TracBrowser for help on using the repository browser.

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: NEMO/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/SBC/sbc_phy.F90 @ 15540

Download in other formats: