Context Navigation

sbc_phy.F90 @ 15551

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

Download in other formats: