forked from virginieguemas/CDlib
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathCDlib.py
1277 lines (1128 loc) · 62.6 KB
/
CDlib.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
# This module contains functions to estimate bulk exchange coefficients from
# observations and to compute them from various parameterisations.
#
# List of functions: - CDN - Neutral momentum bulk transfer coefficient
# - CD - Momentum bulk transfer coefficient
# - CSN - Neutral scalar bulk transfer coefficient
# - CS - Scalar bulk transfer coefficient
# for heat or humidity
# - Z0 - Aerodynamic roughness length
# - ZS - Scalar roughness length for
# heat or humidity
# - U - Wind speed
# - UG - Wind speed accounting for gustiness
# - S - Potential temperature or
# specific humidity
# - LMO - Monin-Obukhov length
# - LMOapprox - Monin-Obukhov length approximation
# - RB - Bulk Richardson number
# - Rstar - Roughness Reynolds number
# - Thetavstar - Scaling virtual temperature
# - WEBB - Webb correction to latent heat flux
# - HSR - Sensible heat flux due to rain
# - TAUR - Wind stress due to rain
# - ZETA - Stability parameter
# - PSI - Stability correction (additive)
# - F - Stability correction as a ratio CD/CDN
# - BULK - Bulk algorithms for turbulent fluxes
#
# Author : Virginie Guemas - 2020
# Modified : Sebastien Blein - January 2021 - Version that accepts unumpy arrays
# storing uncertainties as input and propagates those
# uncertainties until the output, still compatible
# with an usage without uncertainties
################################################################################
import numpy as np
import xarray as xr
import meteolib
import sys
from uncertainties import unumpy as unp
g = 9.81 # Gravity - Fairall et al 1996 use 9.72 instead
k = 0.4 # Von Karman constant
cpw = 4210 # J.K-1.kg-1 # specfic heat of liquid water
cp = 1004 # J.K-1.kg-1 # specific heat of dry air
Rv = 461.4 # J.K-1.kg-1
Ra = 287 # J.K-1.kg-1
dv = 0.219 # cm2.s-1 # Water vapour diffusivity
dh = 0.1906 # cm2.s-1 # Heat diffusivtity in air
################################################################################
def CDN(u=None, ustar=None, f=None, z0=None, z=None) :
"""
This function computes the neutral bulk momentum exchange coefficient CDN either as a function of :
- the aerodynamic roughness length z0 (in m),
- the height z (in m),
or as a function of (the relation between observed wind profile and fluxes):
- the horizontal wind speed u (in m/s),
- the friction velocity ustar (in m/s),
- the multiplicative stability correction f.
Author : Virginie Guemas - October 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
if z0 is not None and z is not None:
z0 = np.where(unp.nominal_values(z0)<=0,np.nan,z0)
z0 = np.where(unp.nominal_values(z0)==np.inf,np.nan,z0)
z = np.where(unp.nominal_values(z)<=unp.nominal_values(z0),np.nan,z)
#
CDn = (k/unp.log(z/z0))**2
elif u is not None and ustar is not None and f is not None:
CD = (ustar**2)/(u**2)
CDn = CD/f
else:
sys.exit('CDn can be computed either from (z0,z) or from (u,ustar,f)')
return CDn
################################################################################
def CD (CDN=None, psi=None) :
"""
This function computes the bulk momentum transfer coefficient as a function of the
neutral bulk momentum transfer coefficient CDN and the additive stability correction
psi.
Author : Virginie Guemas - January 2021
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
if CDN is not None and psi is not None :
Cd = CDN / (1-unp.sqrt(CDN)/k*psi)**2
else:
sys.exit('CDN and psi are required to compute CD')
return Cd
################################################################################
def CSN(deltas=None, u=None, sstar=None, ustar=None, f=None, zs=None, z0=None, z=None) :
"""
This function computes the neutral bulk exchange coefficient for either heat or humidity as a function of :
- the scalar roughness length zs for either temperature or humidity (in m),
- the aerodyamic roughness length z0 (in m),
- the height z (in m),
or as function of (the relation between observed meteorological parameters and fluxes):
- the difference in scalar deltas between height z and the surface (either potential temperature in Kelvin or specific humidity in kg/kg),
- the horizontal wind speed u (in m/s),
- the scaling parameter sstar (thetastar for temperature in Kelvin and qstar for humidity in kg/kg),
- the friction velocity ustar (in m/s),
- the multiplicative stability correction f.
Author : Virginie Guemas - October 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
if zs is not None and z0 is not None and z is not None:
zs = np.where(unp.nominal_values(zs)<=0,np.nan,zs)
zs = np.where(unp.nominal_values(zs)==np.inf,np.nan,zs)
z0 = np.where(unp.nominal_values(z0)<=0,np.nan,z0)
z0 = np.where(unp.nominal_values(z0)==np.inf,np.nan,z0)
z = np.where(unp.nominal_values(z)<=unp.nominal_values(z0),np.nan,z)
z = np.where(unp.nominal_values(z)<=unp.nominal_values(zs),np.nan,z)
CSn = k**2/(unp.log(z/zs)*unp.log(z/z0))
elif ustar is not None and sstar is not None and u is not None and deltas is not None and f is not None:
CS = -(ustar*sstar)/(u*deltas)
CSn = CS/f
else:
sys.exit('CSn can be computed either from (zs,z0,z) or from (u,ustar,deltas,sstar,f')
return CSn
################################################################################
def CS (CDN, CSN, psiM, psiH) :
"""
This function computes the bulk heat or humidity transfer coefficient as a function of the
neutral bulk heat or humidity transfer coefficient CSN, the neutral bulk momentum transfer coefficient CDN
and the additive stability correction psiM for momentum and psiH for heat or humidity.
Author : Virginie Guemas - January 2021
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
if CDN is not None and CSN is not None and psiM is not None and psiH is not None :
mask = np.where(CDN == 0.,True,False)
CDN = np.where(CDN == 0., np.nan, CDN)
Cs = CSN / ((1-CSN/(k*unp.sqrt(CDN))*psiH)*(1-unp.sqrt(CDN)/k*psiM))
Cs = np.where(mask,0.,Cs)
else:
sys.exit('CDN, CSN, psiM and psiH are required to compute CH or CE')
return Cs
################################################################################
def Z0(method=None, u=None, ustar=None, psi=None, CDN=None, z=None, T=None, alpha=None, F=None) :
"""
This function returns the aerodynamic roughness length (in m).
With option method = 'CN', (correspondance between neutral exchange coefficient and roughness length),
the function needs :
- the neutral bulk momentum exchange coefficient CDN,
- the height z (in m).
With option method = 'obs' (from the relation between the observed wind profile and flux),
the function needs :
- the horizontal wind speed u (in m/s),
- the friction velocity ustar (in m/s),
- the additive stability correction psi.
- the height z (in m).
With option method = 'coare2.5' (Smith, 1988 formula used in COARE2.5):
the function needs :
- the alpha = 0.011 as in Charnock (1955), Smith (1988) and COARE2.5 or alpha = 0.013 as in Zeng et al (1998) or alpha = 0.018 as in Beljaars (1995)
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
With option method = 'coare3.0' (Fairall et al 2003),
the function needs :
- the 10m horizontal wind speed u (in m/s),
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
With option method = 'tayloryelland' (Taylor and Yelland, 2001) implemented as an option in COARE3.0,
the function needs :
- the 10m horizontal wind speed u (in m/s),
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
With option method = 'oost' (Oost et al, 2002) implemented as an option in COARE3.0,
the function needs :
- the 10m horizontal wind speed u (in m/s),
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
With option method='andreas05' (Andreas et al, 2005),
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
- the factor F before the exponential (5 in Andreas et al (2005), 1 in Andreas et al (2004))
With option method='andreas10' (Andreas et al, 2010),
- the friction velocity ustar (in m/s),
- the temperature T (in Kelvin).
Author : Virginie Guemas - October 2020
Modified : December 2020 - Virginie Guemas - Option Smith (1988) formula used in COARE 2.5 (Fairall, 1996)
Option COARE 3.0 (Fairall et al 2003)
Options Taylor and Yelland (2001) and Oost et al (2002)
January 2021 - Sebastien Blein - Uncertainty propagation
April 2022 - Virginie Guemas - Option Andreas et al (2005), Andreas et al (2010)
"""
if method == 'CN':
if CDN is not None and z is not None:
CDN = np.where(unp.nominal_values(CDN)==0,np.nan,CDN)
CDN = np.where(unp.nominal_values(CDN)==np.inf,np.nan,CDN)
tmp = unp.sqrt(k**2/CDN)
tmp = np.where(tmp>700.,np.nan,tmp)
z0 = z/unp.exp(tmp)
else:
sys.exit('With option method = \'CN\', input CDN and z are required.')
##################################
elif method == 'obs':
if u is not None and ustar is not None and psi is not None and z is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
tmp = u/ustar*k + psi
tmp = np.where(unp.nominal_values(tmp)>3.e+2,np.nan,tmp)
tmp = np.where(unp.nominal_values(tmp)<-3.e+2,np.nan,tmp)
z0 = z/unp.exp(tmp)
else:
sys.exit('With option method = \'obs\', input z, u, ustar and psi are required.')
##################################
elif method == 'coare2.5':
if alpha is not None and ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
#
z0 = alpha*ustar**2/g + 0.11*meteolib.NU(T)/ustar
else:
sys.exit('With option method = \'coare2.5\', input alpha, ustar and T as required')
##################################
elif method == 'coare3.0':
if u is not None and ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
u10n = u
err = 10
while err>0.01:
alpha = np.where(u10n>18,0.018,np.where(u10n<10,0.011,0.011+(0.018-0.011)/(18-10)*(u10n-10)))
z0 = alpha*ustar**2/g + 0.11*meteolib.NU(T)/ustar
tmp = U(z=10, ustar = ustar, z0 = z0)
err = np.nanmax(abs(u10n-tmp))
u10n = tmp
else:
sys.exit('With option method = \'coare3.0\', input u, ustar and T are required')
##################################
elif method == 'tayloryelland':
# What is coded here corresponds to what is reported in Fairall et al (2003) implemented as an option in COARE3.0
# What appears in SURFEX documentation is not exactly the same: u is used instead of u10n and a different formula
# is used for hs.
if u is not None and ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
u = np.where(unp.nominal_values(u)==0,np.nan,u)
u10n = u
err = 10
while err>0.01:
Tp = 0.729*u10n
hs = 0.0248*u10n**2
Lp = (g*Tp**2)/(2*np.pi)
z0 = 1200*hs*(hs/Lp)**4.5 + 0.11*meteolib.NU(T)/ustar
tmp = U(z=10, ustar = ustar, z0 = z0)
err = np.nanmax(abs(u10n-tmp))
u10n = tmp
else:
sys.exit('With option method = \'tayloryelland\', input u, ustar and T are required')
##################################
elif method == 'oost':
# What is coded here corresponds to what is reported in Fairall et al (2003) implemented as an option in COARE3.0.
# What appears in SURFEX documentation is not exactly the same : u is used instead of U10n
if u is not None and ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
u = np.where(unp.nominal_values(u)==0,np.nan,u)
u10n = u
err = 10
while err>0.01:
Tp = 0.729*u10n
Lp = (g*Tp**2)/(2*np.pi)
Cp = g*Tp/(2*np.pi)
z0 = 50/(2*np.pi)*Lp*(ustar/Cp)**4.5 + 0.11*meteolib.NU(T)/ustar
tmp = U(z=10, ustar = ustar, z0 = z0)
err = np.nanmax(abs(u10n-tmp))
u10n = tmp
else:
sys.exit('With option method = \'oost\', input u, ustar and T are required')
##################################
elif method == 'andreas05':
if ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
z0 = 0.135*meteolib.NU(T)/ustar + 0.035*ustar**2/g*(F*unp.exp(-((ustar-0.18)/0.1)**2)+1)
else:
sys.exit('With option method = \'andreas05\', input ustar and T are required')
##################################
elif method == 'andreas10':
if ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
z0 = 0.135*meteolib.NU(T)/ustar + 0.00023*unp.tanh(13*ustar)**3
else:
sys.exit('With option method = \'andreas10\', input ustar and T are required')
else:
sys.exit('Valid methods are \'CN\', \'obs\', \'coare2.5\', \'coare3.0\',\'tayloryelland\',\'oost\',\'andreasà5\',\'andreas10\'.')
return z0
################################################################################
def ZS(method=None, deltas=None, sstar=None, psi=None, CSN=None, z0=None, z=None, rstar=None, ustar=None, T=None, s=None) :
"""
This function returns the scalar roughness length (in m) for either heat or humidity.
With option method = 'CN' (correspondance between neutral exchange coefficients and roughness length),
the function needs :
- the neutral bulk scalar exchange coefficient CHN (for heat) or CQN (for humidity),
- the aerodyamic roughness length z0 (in m),
- the height z (in m).
With option method = 'obs' (from the relation between the observed scalar profiles and fluxes),
the function needs :
- the height z (in m),
- the difference in scalar deltas between height z and the surface (either potential temperature in Kelvin or specific humidity in kg/kg)
- the scaling parameter sstar (for either temperature thetastar or humidity qstar),
- the additive stability correction psi.
With option method = 'LKB' (Liu et al 1979 used in COARE2.5),
the function needs:
- the temperature T (in Kelvin),
- the roughness Reynolds number rstar,
- the friction velocity ustar (in m/s).
- s = 'T'/'Q' for heat/humidity (parameters a and b differ).
With option method = 'andreas' (model from Andreas, 1987),
the function needs:
- the roughness Reynolds number rstar,
- the aerodynamic roughness length z0 (in m),
- s = 'T'/'Q' for heat/humidity (parameters b0, b1 and b2 differ).
With option method = 'brutsaertgarratt' (Garratt, 1992 inspired from Brutsaert 1975 model),
the function needs:
- the roughness Reynolds number rstar,
- the aerodynamic roughness length z0 (in m),
With option method = 'revisedbrutsaertgarratt' (updated coefficients as proposed in Fairall, 2003),
the function needs:
- the roughness Reynolds number rstar,
- the aerodynamic roughness length z0 (in m),
With option method = 'simplebrutsaert' (simplified Brutsaert 1982 model used in several bulk algorithms presented in Fairall 2003),
the function needs:
- the temperature T (in Kelvin),
- the friction velocity ustar (in m/s).
- s = 'T'/'Q' for heat/humidity (parameter r differ).
With option method = 'clayson' (Clayson et al 1996),
the function needs: Not coded yet
With option method = 'mondonredelsperger' (Mondon and Redelsperger, 1998) proposed in SURFEX,
the function needs:
- the temperature T (in Kelvin),
- the friction velocity ustar (in m/s).
- s = 'T'/'Q' for heat/humidity (parameter r differ).
With option method = 'coare3.0' (Fairall et al, 2003),
the function needs:
- the roughness Reynolds number rstar.
Author : Virginie Guemas - October 2020
Modified : December 2020 - Virginie Guemas - Option LKB (Liu et al, 1979) used in COARE2.5 (Fairall et al 1996)
Option Andreas (1987)
Option Brutsaert (1982) simplified model
Option COARE 3.0 (Fairall et al, 2003) and (Garratt, 1992) and revised (Garratt, 1992) as suggested in Fairall et al, 2003
January 2021 - Sebastien Blein - Uncertainty propagation
"""
if method == 'CN':
if CSN is not None and z0 is not None and z is not None:
z0 = np.where(unp.nominal_values(z0)<=0,np.nan,z0)
z0 = np.where(unp.nominal_values(z0)==np.inf,np.nan,z0)
z = np.where(unp.nominal_values(z)<=unp.nominal_values(z0),np.nan,z)
CSN = np.where(unp.nominal_values(CSN)<=0,np.nan,CSN)
CSN = np.where(unp.nominal_values(CSN)==np.inf,np.nan,CSN)
zs = z/unp.exp(k**2/(unp.log(z/z0)*CSN))
else:
sys.exit('With option method = \'CN\', input CSN, z0 and z are required.')
###################################
elif method == 'obs':
if z is not None and deltas is not None and sstar is not None and psi is not None:
sstar = np.where(unp.nominal_values(sstar)==0,np.nan,sstar)
tmp = unp.exp(deltas/sstar*k + psi)
tmp = np.where(unp.nominal_values(tmp)==0,np.nan,tmp)
zs = z/tmp
else:
sys.exit('With option method = \'obs\', input z, deltas, sstar and psi are required.')
###################################
elif method == 'LKB':
if rstar is not None and ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
if s == 'T':
a = np.where(rstar<0, np.nan, np.where(rstar<0.11, 0.177, np.where(rstar<0.825, 1.376, np.where(rstar<3., 1.026, np.where(rstar<10., 1.625, np.where(rstar<30., 4.661, 34.904))))))
b = np.where(rstar<0, np.nan, np.where(rstar<0.11, 0, np.where(rstar<0.825, 0.929, np.where(rstar<3., -0.599, np.where(rstar<10., -1.018, np.where(rstar<30., -1.475, -2.067))))))
# There is a gap in Liu et al (1979) - no values for a/b are given for 0.825 < rstar < 0.925.
# Surely a typo, but between 0.825 and 0.925, I arbitrarily chose 0.825 as boundary between two
# consecutive segments.
a = np.where(np.isnan(rstar.astype(float)),np.nan,a)
b = np.where(np.isnan(rstar.astype(float)),np.nan,b)
elif s == 'Q':
a = np.where(rstar<0, np.nan, np.where(rstar<0.11, 0.292, np.where(rstar<0.825, 1.808, np.where(rstar<3., 1.393, np.where(rstar<10., 1.956, np.where(rstar<30., 4.994, 30.79))))))
b = np.where(rstar<0, np.nan, np.where(rstar<0.11, 0, np.where(rstar<0.825, 0.826, np.where(rstar<3., -0.528, np.where(rstar<10., -0.870, np.where(rstar<30., -1.297, -1.845))))))
a = np.where(np.isnan(rstar.astype(float)),np.nan,a)
b = np.where(np.isnan(rstar.astype(float)),np.nan,b)
else:
sys.exit('s should be T for temperature or Q for humidity')
zs = meteolib.NU(T)/ustar* a*rstar**b
else:
sys.exit('With option method = \'LKB\', input rstar, ustar, T and s are required.')
###################################
elif method == 'andreas':
if rstar is not None and z0 is not None:
rstar = np.where(unp.nominal_values(rstar)<=0,np.nan,rstar)
if s == 'T':
b0 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 1.25, np.where(rstar<2.5, 0.149, 0.317)))
b1 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 0., np.where(rstar<2.5, -0.55, -0.565)))
b2 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 0., np.where(rstar<2.5, 0., -0.183)))
b0 = np.where(np.isnan(rstar.astype(float)),np.nan,b0)
b1 = np.where(np.isnan(rstar.astype(float)),np.nan,b1)
b2 = np.where(np.isnan(rstar.astype(float)),np.nan,b2)
elif s == 'Q':
b0 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 1.61, np.where(rstar<2.5, 0.351, 0.396)))
b1 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 0., np.where(rstar<2.5, -0.628, -0.512)))
b2 = np.where(rstar<0, np.nan, np.where(rstar<=0.135, 0., np.where(rstar<2.5, 0., -0.18)))
b0 = np.where(np.isnan(rstar.astype(float)),np.nan,b0)
b1 = np.where(np.isnan(rstar.astype(float)),np.nan,b1)
b2 = np.where(np.isnan(rstar.astype(float)),np.nan,b2)
else:
sys.exit('s should be T for temperature or Q for humidity')
zs = z0*unp.exp(b0 + b1*unp.log(rstar) + b2*unp.log(rstar)**2)
else:
sys.exit('With option method = \'andreas\', input rstar, z0 and s are required.')
##################################
elif method == 'brutsaertgarratt':
if rstar is not None and z0 is not None:
zs = z0*unp.exp(2-2.28*rstar**0.25)
else:
sys.exit('With option method = \'brutsaertgarratt\', input rstar and z0 are required.')
##################################
elif method == 'revisedbrutsaertgarratt':
if rstar is not None and z0 is not None:
zs = z0*unp.exp(3.4-3.5*rstar**0.25)
else:
sys.exit('With option method = \'revisedbrutsaertgarratt\', input rstar and z0 are required.')
##################################
elif method == 'simplebrutsaert':
if ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
if s == 'T':
r = 0.4
elif s == 'Q':
r = 0.62
else:
sys.exit('s should be T for temperature or Q for humidity')
zs = r*meteolib.NU(T)/ustar
else:
sys.exit('With option method = \'simplebrutsaert\', input ustar, T and s are required.')
##################################
elif method == 'clayson':
sys.exit('Sorry. Option method = \'clayson\' from Clayson et al (1996) is not coded yet. Coming soon')
##################################
elif method == 'mondonredelsperger':
if ustar is not None and T is not None:
ustar = np.where(unp.nominal_values(ustar)==0,np.nan,ustar)
if s == 'T':
zs = np.where(ustar>0.23, 0.14*meteolib.NU(T)/(ustar-0.2) + 7*10**(-6), 0.015*ustar**2/g+0.18*meteolib.NU(T)/ustar)
elif s == 'Q':
zs = np.where(ustar>0.23, 0.2*meteolib.NU(T)/(ustar-0.2) + 9*10**(-6), 0.0205*ustar**2/g+0.294*meteolib.NU(T)/ustar)
else:
sys.exit('s should be T for temperature or Q for humidity')
zs = np.where(np.isnan(ustar.astype(float)),np.nan,zs)
else:
sys.exit('With option method = \'mondonredelsperger\', input ustar, T and s are required.')
##################################
elif method == 'coare3.0':
if rstar is not None:
zs = 0.000055*rstar**(-0.6)
zs = np.where(zs<0.0001,0.0001,zs)
else:
sys.exit('With option method = \'coare3.0\', input rstar is required.')
else:
sys.exit('Valid methods are \'CN\', \'obs\', \'LKB\', \'andreas\',\'brutsaertgarratt\',\'revisedbrutsaertgarratt\',\'simplebrutsaert\',\'mondonredelsperger\',\'coare3.0\'.')
return zs
################################################################################
def U(z, ustar, z0, psi=0) :
"""
This funcion computes the wind speed (in m/s) at height z (in m) as a fonction of the friction velocity in ustar in (m/s), the aerodynamic roughness length z0 (in m), and the additive stability correction psi.
Author : Virginie Guemas - October 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
z0 = np.where(unp.nominal_values(z0)<=0,np.nan,z0)
z0 = np.where(unp.nominal_values(z0)==np.inf,np.nan,z0)
z = np.where(unp.nominal_values(z)<unp.nominal_values(z0),np.nan,z)
u = ustar/k * (unp.log(z/z0) - psi)
return u
################################################################################
def UG(method=None, u=None, h=None, Q0v=None, thetav=None, Q0=None, E0=None, T=None, beta=1.25, zeta=None):
"""
This function computes the corrected wind speed to account for gustiness. It needs :
- the stability parameter zeta to verify the stability conditions to apply a particular gustiness correction
are verified.
With option method = 'godfreybeljaars', (Godfrey and Beljaars, 1991)
the function needs :
- the horizontal wind speed u (in m/s),
- the correction factor beta = 1.25 as in COARE2.5 or beta = 1 as in Zeng et al (1998) or beta = 1.2 as in Beljaars (1995)
- the surface virtual temperature flux Q0v (K.m.s-1)
- the convective boundary layer height h (in m), typically h=600m,
- the virtual potential temperature thetav (in Kelvin).
With option method = 'fairall' (approximation of godfreyBeljaars used in Fairall et al (1996, 2003),
the function needs :
- the horizontal wind speed u (in m/s),
- the correction factor beta = 1.25 as in COARE2.5 or beta = 1 as in Zeng et al (1998) or beta = 1.2 as in Beljaars (1995)
- the temperature T (in Kelvin),
- the potential temperature flux Q0 (K.m.s-1),
- the humidity flux E0 (m.s-1),
- the convective boundary layer height h (in m), typically h=600m.
With option method = 'jordan' (Jordan et al, 1999) used in Andreas et al (2010),
the function needs:
- the horizontal wind speed u (in m/s).
Author : Virginie Guemas - December 2020
Modified : January 2021 - Virginie Guemas - Option fairall which approximates godfreybeljaars
January 2021 - Sebastien Blein - Uncertainty propagation
Verify whether stability allows gustiness correction
"""
if zeta is None:
sys.exit('Please provide zeta so that the function can verify whether the method selected is suitable for the stability')
if method == 'godfreybeljaars':
if u is not None and thetav is not None and h is not None and Q0v is not None:
# Both 'godfreybeljaars' and 'fairall' methods to estimate the wind
# gustiness correction are valid only in unstable cases.
# Both conditions on zeta and Q0v are applied in case zeta is computed
# using some other variables than Q0v (thetastar for example).
# A Q0v_pos array has to be created as the np.where function calculates values everywhere. Otherwise the uncertainties module would rise math errors for the 1/3 power with negative Q0v values.
Q0v_pos = np.where(Q0v>0, Q0v, np.nan)
ug = np.where((zeta<0)&(Q0v>0), beta * (g/thetav*h*Q0v_pos)**(1/3), 0)
ucor = unp.sqrt(u**2+ug**2)
else:
sys.exit('With option method = \'godfreybeljaars\', input u, thetav, h and Q0v are required.')
elif method == 'fairall':
if u is not None and T is not None and h is not None and Q0 is not None and E0 is not None:
Q0v = Q0+0.61*T*E0
Q0v_pos = np.where(Q0v>0, Q0v, np.nan)
ug = np.where((zeta<0)&(Q0v>0), beta * (g/T*h*Q0v_pos)**(1/3), 0.)
ucor = unp.sqrt(u**2+ug**2)
else:
sys.exit('With option method = \'fairall\', input u, T, h, Q0 and E0 are required.')
elif method == 'jordan':
if u is not None:
ucor = np.where(zeta>0, u + 0.5/unp.cosh(u), u)
else:
sys.exit('With option method = \'jordan\', input u is required.')
else:
sys.exit('Valid methods are \'godfreybeljaars\', \'fairall\' and \'jordan\'.')
return ucor
################################################################################
def S(z, s0, sstar, zs, psi=0) :
"""
This function computes either the potential temperature theta (in Kelvin) or the specific humidity (in kg/kg) as a function of :
- the scaling parameter sstar (thetastar for temperature in Kelvin and qstar for humidity in kg/kg),
- the surface parameter s0 (surface potential temperature in Kelvin or surface specific humidity in kg/kg),
- the scalar roughness length for heat or humidity (in m),
- the additive stability correction psi,
- the height z (in m).
Author : Virginie Guemas - October 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
zs = np.where(unp.nominal_values(zs)<=0,np.nan,zs)
zs = np.where(unp.nominal_values(zs)==np.inf,np.nan,zs)
z = np.where(unp.nominal_values(z)<unp.nominal_values(zs),np.nan,z)
s = s0 + sstar/k * (unp.log(z/zs) - psi)
return s
################################################################################
def LMO(ustar, thetav, thetavstar=None, Q0v=None) :
"""
This function computes the Monin-Obukhov length (in meters) either as a function of the virtual temperature scaling parameter thetavstar (in Kelvin) or as a function of the surface virtual temperature flux (K.m.s-1). It also requires the friction velocity (in m.s-1) and the layer-average virtual potential temperature (in Kelvin).
Author : Virginie Guemas - October 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
meteolib.check_T(thetav)
if Q0v is None and thetavstar is None:
sys.exit('At least one of thetavstar or Q0v should be provided')
beta = g/thetav
if Q0v is not None:
Q0v = np.where(unp.nominal_values(Q0v)==0,np.nan,Q0v)
Lmo = -(ustar**3)/(k*beta*Q0v)
if thetavstar is not None:
Qvbis = - ustar*thetavstar
if np.nanmax(np.abs((Qvbis-Q0v)/Q0v)) > 0.001 :
sys.exit('Q0v and -ustar*thetavstar differ')
else:
thetavstar = np.where(unp.nominal_values(thetavstar)==0,np.nan,thetavstar)
Lmo = (ustar**2)/(k*beta*thetavstar)
return Lmo
################################################################################
def LMOapprox(ustar, T, thetastar=None, qstar=None, Q0=None, E0=None) :
"""
This function approximates the Monin-Obukhov length (in meters) as in Fairall et al (1996, 2003) either as of :
- the friction velocity ustar (in m.s-1),
- the temperature scaling parameter thetastar (in Kelvin),
- the humidity scaling parameter qstar (in kg.kg-1),
- the temperature (in Kelvin).
or as a function of :
- the friction velocity ustar (in m.s-1),
- the turbulent heat flux Q0 (in K.m.s-1),
- the turbulent humidity flux E0 (in m.s-1),
- the temperature (in Kelvin).
Author : Virginie Guemas - January 2021
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
meteolib.check_T(T)
if thetastar is None and Q0 is None:
sys.exit('At least one of thetastar or Q0 should be provided')
if qstar is None and E0 is None:
sys.exit('At least one of qstar or E0 should be provided')
beta = g/T
if Q0 is not None and E0 is not None:
Q0v = np.where(unp.nominal_values(Q0+0.61*T*E0)==0,np.nan,Q0+0.61*T*E0)
Lmo = -(ustar**3)/(k*beta*Q0v)
elif thetastar is not None and qstar is not None:
thetavstar = np.where(unp.nominal_values(thetastar+0.61*T*qstar)==0,np.nan,thetastar+0.61*T*qstar)
Lmo = (ustar**2)/(k*beta*thetavstar)
elif Q0 is not None and qstar is not None:
Q0v = np.where(unp.nominal_values(Q0+0.61*T*(-ustar*qstar))==0,np.nan,Q0+0.61*T*(-ustar*qstar))
Lmo = -(ustar**3)/(k*beta*Q0v)
else:
Q0v = np.where(unp.nominal_values((-ustar*thetastar)+0.61*T*E0)==0,np.nan,(-ustar*thetastar)+0.61*T*E0)
Lmo = -(ustar**3)/(k*beta*Q0v)
return Lmo
################################################################################
def RB(thetav, Dthetav, u, v, zt, zu) :
"""
This function computes the Bulk Richardson number as a function of the virtual potential temperature thetav (in Kelvin) at height zt (in m) and wind speed horizontal components u and v (in m/s) at height zu (in m) and the difference in virtual potential temperature between the height zt and the surface Dthetav (in Kelvin).
Author : Virginie Guemas - October 2020
"""
meteolib.check_T(thetav)
wspd2 = np.where(unp.nominal_values(u**2+v**2)==0,np.nan,u**2+v**2)
Rb = g/thetav * (Dthetav*zu**2)/(wspd2*zt)
return Rb
################################################################################
def Rstar(ustar, z0, T) :
"""
This function computes the Roughness Reynolds number as a function of the friction velocity ustar (in m/s), the aerodynamic roughness length (in m) and the temperature (in Kelvin).
Author : Virginie Guemas - October 2020
"""
z0 = np.where(unp.nominal_values(z0)<=0,np.nan,z0)
z0 = np.where(unp.nominal_values(z0)==np.inf,np.nan,z0)
Rs = ustar*z0/meteolib.NU(T)
return Rs
################################################################################
def Thetavstar(thetastar, qstar, theta, q) :
"""
This function estimates the scaling virtual temperature parameter (in Kelvin) from the scaling temperature parameter (in Kelvin), the scaling humidity parameter (in kg/kg), the layer-average temperature (in Kelvin) and the layer-average specific humidity (in kg/kg).
Author : Virginie Guemas - October 2020
"""
meteolib.check_q(q)
meteolib.check_q(qstar)
meteolib.check_T(theta)
thetavstar = thetastar*(1+0.61*q) + 0.61*theta*qstar
return thetavstar
################################################################################
def WEBB(E0,Q0,q,T) :
"""
The total surface humidity flux contains a turbulent component and a component from the mean vertical wind speed. This mean vertical wind speed can be approximated as in Webb et al (1980). This function estimates the total surface humidity flux including the Webb effect as a function of :
- the turbulent humidity flux E0 (in m.s-1),
- the turbulent temperature flux Q0 (in K.m.s-1),
- the specific humidity q (in kg.kg-1),
- the temperature T (in Kelvin).
Author : Virginie Guemas - January 2021
"""
meteolib.check_q(q)
meteolib.check_T(T)
w = 1.61*E0 + (1.61*q)*Q0/T
E0cor = EO + w*q
return E0cor
################################################################################
def HSR(R, T, Ts, deltaq, P, lda=1, mu=dh/dv) :
"""
Precipitation transfers heat to the ocean which should be accounted for in the sensible
heat flux term. This function estimates the correction for liquid precipitation which
should be added to sensible heat flux toward the ocean according to Gosnell et al (1995)
as used in Fairall et al (1996).
This function needs :
- the rain rate R (kg.s-1.m-2),
- the rain (supposed equal to atmospheric) temperature T (in Kelvin),
- the surface temperature Ts (in Kelvin),
- the difference in specific humidity between atmosphere and surface deltaq (in kg.kg-1),
- the pressure (in hPa),
- the lambda parameter lda set to 1 for COARE, Rv/Ra for ECUME and 1 for ECUME6,
- the mu parameter set to 1 for COARE, dh/dv for ECUME and ECUME6. dh/dv seems to be the
correct according to Gosnell et al (1995) derivation but the last line includes an error
which was repeated in COARE. The default here is dh/dv which can be overcritten by setting
mu = 1.
Author : Virginie Guemas - January 2021
"""
deltaT = Ts-T # Difference between atmosphere/rain temperature and surface one
Lv = meteolib.LV(T) # Latent heat of vaporization of rain
qsat = meteolib.ES(T)*Ra/(P*Rv) # Saturation specific humidity
dq = lda * Lv*qsat/(Rv*T**2) # Clausius-Clapeyron relation
alpha = 1/(1 + Lv/cp * dv/dh* dq) # wet-bulb factor
deltaq = np.where(unp.nominal_values(deltaq)==0,np.nan,-deltaq)
B = mu*(cp/Lv)*deltaT/deltaq # Bowen ratio
B = np.where(unp.nominal_values(B)==0,np.nan,B)
Hsrain = - R*cpw*alpha*(1+1/B)*deltaT
return Hsrain
################################################################################
def TAUR(u, R, gamma=0.85) :
"""
Rainfall tends to increase surface drag. This functions computes the correction
to be added to surface wind stress according to Fairall et al (1996) as a
function of :
- the wind speed u at 10m (in m.s-1)
- the precipitation rate (in kg.m-2.s-1),
- the gamma parameter set to 0.85 in COARE and ECUME6, 1 in ECUME.
Author : Virginie Guemas - January 2021
"""
Taurain = gamma*R*u
return Taurain
################################################################################
def ZETA(z,Lmo) :
"""
This function computes the Monin-Obukhov stability parameter zeta as a
function of the height z and the Monin-Obukhov length Lmo.
Author : Sebastien Blein - January 2021
"""
Lmo = np.where(unp.nominal_values(Lmo)==0,np.nan,Lmo)
zeta = z/Lmo
return zeta
################################################################################
def PSI(zetaT, zetaU, gamma=5, stab=None, unstab=None) :
"""
This function computes a stability correction as a function of Monin-Obukhov length. It takes four arguments:
- zetaT = z (height for temperature/humidity) / Lmo (Monin-Obukhov length)
- zetaU = z (height for wind) / Lmo (Monin-Obukhov length)
- stab = formulation for stable regimes : 'dyer-hicks' -- Dyer and Hicks (1970)
'lettau' -- Lettau (1979)
'holtslag-bruin' -- Holtslag and de Bruin (1988)
'beljaars-holtslag' -- Beljaars and Holtslag (1991)
'grachev' -- Grachev (2007)
- unstab = formulation for unstable regimes : 'businger-dyer' -- Paulson (1970)
'fairall1996' -- Fairall et al (1996)
'grachev2000' -- Grachev et al (2000)
- the gamma factor for the 'dyer-hicks' option which can range between about 5
(Webb, 1970; Businger et al., 1971; Dyer, 1974; Large and Pond,1981) and about 7 (Wieringa, 1980;
Large and Pond, 1982; Högström, 1988) - 4.7 seems to be used in COARE2.5
Author : Virginie Guemas - September 2020
Modified : December 2020 - Sebastien Blein - correct Beljaars and Holtslag 1991
December 2020 - Virginie Guemas - add kansas, fairall and holtslag-bruin
January 2021 - Virginie Guemas - include factor gamma to tune the dyer-hicks option. Ex: 4.7 for COARE2.5
January 2021 - Sebastien Blein - Uncertainty propagation
"""
np.seterr(invalid='ignore')
if stab is None or unstab is None:
sys.exit('Stability correction type has to be specified (e.g.: stab=\'beljaars-holtslag\' and unstab=\'businger-dyer\')')
###############################
# zeta discrimination between stab and unstab in order to avoid error with unp function when math domain is not Real.
zeta_stabT = np.where(zetaT>0,zetaT,np.nan)
zeta_unstabT = np.where(zetaT<0,zetaT,np.nan)
zeta_stabU = np.where(zetaU>0,zetaU,np.nan)
zeta_unstabU = np.where(zetaU<0,zetaU,np.nan)
###############################
if unstab == 'businger-dyer':
phiM = (1 - 16*zeta_unstabU)**(-0.25)
phiH = (1 - 16*zeta_unstabT)**(-0.5)
psiM = 2*unp.log((1+phiM**(-1))/2) + unp.log((1+phiM**(-2))/2) - 2*unp.arctan(phiM**(-1)) + np.pi/2
psiH = 2*unp.log((1+phiH**(-1))/2)
# psi is np.nan when zeta >= 0
################################
elif unstab == 'holtslag1990':
# I need to integrate the phi
psiM = None
psiH = None
################################
elif unstab == 'kaderyaglom1990':
# I need to integrate the phi
psiM = None
psiH = None
################################
elif unstab == 'fairall1996':
yU = (1 - 12.87*zeta_unstabU)**(1/3)
yT = (1 - 12.87*zeta_unstabT)**(1/3)
tmpU = 1.5*unp.log((yU**2+yU+1)/3) - unp.sqrt(3)*unp.arctan((2*yU+1)/unp.sqrt(3)) + np.pi/unp.sqrt(3)
tmpT = 1.5*unp.log((yT**2+yT+1)/3) - unp.sqrt(3)*unp.arctan((2*yT+1)/unp.sqrt(3)) + np.pi/unp.sqrt(3)
psi = PSI(zeta_unstabT, zeta_unstabU, stab = stab, unstab = 'businger-dyer')
psiM = 1/(1+zeta_unstabU**2)*psi[0] + zeta_unstabU**2/(1+zeta_unstabU**2)*tmpU
psiH = 1/(1+zeta_unstabT**2)*psi[1] + zeta_unstabT**2/(1+zeta_unstabT**2)*tmpT
# psi is np.nan when zeta >= 0
################################
elif unstab == 'grachev2000':
a={'m':10.15,'h':34.15}
yT = (1 - a['h']*zeta_unstabT)**(1/3)
yU = (1 - a['m']*zeta_unstabU)**(1/3)
tmpT = 1.5*unp.log((yT**2+yT+1)/3) - unp.sqrt(3)*unp.arctan((2*yT+1)/unp.sqrt(3)) + np.pi/unp.sqrt(3)
tmpU = 1.5*unp.log((yU**2+yU+1)/3) - unp.sqrt(3)*unp.arctan((2*yU+1)/unp.sqrt(3)) + np.pi/unp.sqrt(3)
psi = PSI(zeta_unstabT, zeta_unstabU, stab = stab, unstab = 'businger-dyer')
psiM = 1/(1+zeta_unstabU**2)*psi[0] + zeta_unstabU**2/(1+zeta_unstabU**2)*tmpU
psiH = 1/(1+zeta_unstabT**2)*psi[1] + zeta_unstabT**2/(1+zeta_unstabT**2)*tmpT
# psi is np.nan when zeta >= 0
################################
elif unstab == 'beljaars-holtslag':
a,b,c,d = 1.,0.667,5.,0.35
psiM = -(a * zeta_unstabU + b*(zeta_unstabU - c/d) * unp.exp(-d * zeta_unstabU) + b*c/d)
psiH = -(a * zeta_unstabT + b*(zeta_unstabT - c/d) * unp.exp(-d * zeta_unstabT) + b*c/d)
# Beljaars and Holtslag do not propose an option for unstable cases but only for stable cases. This formulation is here
# only to reproduce Elvidge et al (2016).
else:
sys.exit('This option for unstable cases is not coded yet.')
################################
################################
if stab == 'dyer-hicks':
psiM = np.where (zetaU>0, -gamma*zeta_stabU, psiM)
psiH = np.where (zetaT>0, -gamma*zeta_stabT, psiH)
################################
elif stab == 'lettau':
xT = 1 + 4.5*zeta_stabT
xU = 1 + 4.5*zeta_stabU
psiM = np.where (zetaU>0, unp.log((xU**0.5+1)/2) + 2*unp.log((xU**0.25+1)/2) -2*unp.arctan(xU**0.25) + np.pi/2 + 4/3*(1-xU**0.75), psiM)
psiH = np.where (zetaT>0, 2*unp.log((xT**0.25+1)/2) -2*(xT**1.5/3 + xT**0.5 - 4/3), psiH)
################################
elif stab == 'holtslag-bruin':
psiT = -3.75/0.35 - 0.7*zeta_stabT + 3.75/0.35*unp.exp(-0.35*zeta_stabT) -0.75*zeta_stabT*unp.exp(-0.35*zeta_stabT)
psiU = -3.75/0.35 - 0.7*zeta_stabU + 3.75/0.35*unp.exp(-0.35*zeta_stabU) -0.75*zeta_stabU*unp.exp(-0.35*zeta_stabU)
psiM = np.where (zetaU>0, psiU, psiM)
psiH = np.where (zetaT>0, psiT, psiH)
################################
elif stab == 'holtslag1990':
# I need to integrate the phi
psiM = None
psiH = None
################################
elif stab == 'beljaars-holtslag':
a,b,c,d = 1.,0.667,5.,0.35
psiM = np.where (zetaU>0, -(a * zeta_stabU + b*(zeta_stabU - c/d) * unp.exp(-d * zeta_stabU) + b*c/d), psiM)
psiH = np.where (zetaT>0, -((1+2.*a*zeta_stabT/3)**(3./2) + b*(zeta_stabT-c/d)*unp.exp(-d*zeta_stabT) + b*c/d - 1), psiH)
################################
elif stab == 'grachev':
xU = (zeta_stabU+1)**(1/3)
psiM = np.where (zetaU>0, -19.5*(xU-1) + 3.25*0.3**(1/3)*(2*unp.log((xU+0.3**(1/3))/(1+0.3**(1/3))) - unp.log((xU**2-0.3**(1/3)*xU+0.3**(2/3))/(1-0.3**(1/3)+0.3**(2/3))) +2*unp.sqrt(3)*(unp.arctan((2*xU-0.3**(1/3))/(unp.sqrt(3)*0.3**(1/3))) -unp.arctan((2-0.3**(1/3))/(unp.sqrt(3)*0.3**(1/3))))) , psiM)
psiH = np.where (zetaT>0, -2.5*unp.log(1+3*zeta_stabT+zeta_stabT**2) +5/(2*unp.sqrt(5))*(unp.log((2*zeta_stabT+3-unp.sqrt(5))/(2*zeta_stabT+3+unp.sqrt(5))) -unp.log((3-unp.sqrt(5))/(3+unp.sqrt(5)))), psiH)
################################
else:
sys.exit('This option for stable cases is not coded yet.')
################################
################################
# To avoid np.nan for zeta = 0
psiM = np.where (zetaU == 0., 0., psiM)
psiH = np.where (zetaT == 0., 0., psiH)
return (psiM,psiH)
################################################################################
def F(Rb, CDN, z, var='momentum', author='Louis') :
"""
This function computes the CD/CDN ratio following Louis (1979) with either Louis original functions or later adjustments by other authors. It takes five arguments:
- Rb = Richardson bulk number
- CDN = neutral drag coefficient
- z = height
- var = 'momentum' / 'heat'
- author = 'Louis' / 'LupkesGryanik'
Author : Virginie Guemas - September 2020
Modified : January 2021 - Sebastien Blein - Uncertainty propagation
"""
z0=CDlib.z0(CDN)
if author == 'Louis' :
c1 = 9.4
if var == 'momentum' :
alpha = 7.4
elif var == 'heat' or var == 'moisture' :
alpha = 5.3
else:
sys.exit('var can be momentum, heat or moisture only')
z0 = np.where(unp.nominal_values(z0)==0,np.nan,z0)
c2 = c1*alpha*CDN*(z/z0)**0.5
tmp = np.where(unp.nominal_values(1+c1/2*Rb)==0,np.nan,1+c1/2*Rb)
fstab = 1/tmp**2
elif author == 'LupkesGryanik' :
if var == 'momentum' :
c1 = 10
alpha = 7.5
elif var == 'heat' or var == 'moisture' :
c1 = 15
alpha = 5
else:
sys.exit('var can be momentum, heat or moisture only')
z0 = np.where(unp.nominal_values(z0)==0,np.nan,z0)
c2 = c1*alpha*CDN*(z/z0+1)**0.5
Rb = np.where(unp.nominal_values(Rb)==-1,np.nan,Rb)
tmp = np.where(unp.nominal_values(1+10*Rb/unp.sqrt(Rb+1))==0,np.nan,1+10*Rb/unp.sqrt(Rb+1))
fstab = 1/tmp
else:
sys.exit('Only Louis and LupkesGryanik version are coded for now')
dunstab = np.where(unp.nominal_values(1+c2*np.abs(Rb)**0.5)==0,np.nan,1+c2*np.abs(Rb)**0.5)
f = np.where(Rb < 0, 1 - (c1*Rb)/dunstab, fstab)
return f
################################################################################
def BULK(zu, zt, u, theta, thetas, q, qs, T, method='coare2.5') :
"""
This function applies one of the COARE or ECUME algorithms to estimate bulk
turbulent fluxes of velocity, heat and humidity as well as the associated transfer
coefficients above ocean. It can also use a combination of the most recent parametrizations
of scalar and aerodynamic roughness and stability function over sea ice to estimate
those bulk turbulent fluxes and transfer coefficients over sea ice. It needs :
- the atmospheric measurement height zu for wind (in m),
- the atmospheric measurement height zt for temperature and humidity (in m),
- the wind speed at height zu (in m.s-1),
- the potential temperature at height zt (in Kelvin),
- the potential temperature at the surface (in Kelvin),
- the specific humidity at height zt (in kg.kg-1),
- the specific humidity at the surface (in kg.kg-1) - with the 2% reduction over sea
- the layer-averaged temperature T (in Kelvin).
- the method : 'coare2.5', 'coare3.0' or 'seaice'. Default : 'coare2.5'.
Warning: The Webb correction for latent heat flux, the precipitation correction for
sensible heat flux, the warm-layer and cool-skin corrections over ocean for surface
temperature are not included.
Author : Virginie Guemas - January 2021