phys/module_sf_exchcoef.F

   1 !  This MODULE holds the routines that calculate air-sea exchange coefficients
   2
   3 MODULE module_sf_exchcoef
   4 CONTAINS
   5
   6   SUBROUTINE znot_m_v1(uref,znotm)
   7   IMPLICIT NONE
   8
   9 ! uref(m/s)   :   Reference level wind
  10 ! znotm(meter):   Roughness scale for momentum
  11 ! Author      :  Biju Thomas on 02/07/2014
  12 !
  13
  14     REAL, INTENT(IN) :: uref
  15     REAL, INTENT(OUT):: znotm
  16     REAL             :: bs0, bs1, bs2, bs3, bs4, bs5, bs6
  17     REAL             :: cf0, cf1, cf2, cf3, cf4, cf5, cf6
  18
  19
  20     bs0 = -8.367276172397277e-12
  21     bs1 = 1.7398510865876079e-09
  22     bs2 = -1.331896578363359e-07
  23     bs3 = 4.507055294438727e-06
  24     bs4 = -6.508676881906914e-05
  25     bs5 = 0.00044745137674732834
  26     bs6 = -0.0010745704660847233
  27
  28     cf0 = 2.1151080765239772e-13
  29     cf1 = -3.2260663894433345e-11
  30     cf2 = -3.329705958751961e-10
  31     cf3 = 1.7648562021709124e-07
  32     cf4 = 7.107636825694182e-06
  33     cf5 = -0.0013914681964973246
  34     cf6 = 0.0406766967657759
  35
  36
  37     IF ( uref .LE. 5.0 ) THEN
  38       znotm = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
  39     ELSEIF (uref .GT. 5.0 .AND. uref .LT. 10.0) THEN
  40       znotm =.00000235*(uref**2 - 25 ) + 3.805129199617346e-05
  41     ELSEIF ( uref .GE. 10.0  .AND. uref .LT. 60.0) THEN
  42       znotm = bs6 + bs5*uref + bs4*uref**2 + bs3*uref**3 + bs2*uref**4 +  &
  43               bs1*uref**5 + bs0*uref**6
  44     ELSE
  45       znotm = cf6 + cf5*uref + cf4*uref**2 + cf3*uref**3 + cf2*uref**4 +  &
  46               cf1*uref**5 + cf0*uref**6
  47
  48     END IF
  49
  50   END SUBROUTINE znot_m_v1
  51
  52   SUBROUTINE znot_m_v0(uref,znotm)
  53   IMPLICIT NONE
  54
  55 ! uref(m/s)   :   Reference level wind
  56 ! znotm(meter):   Roughness scale for momentum
  57 ! Author      :  Biju Thomas on 02/07/2014
  58
  59     REAL, INTENT(IN) :: uref
  60     REAL, INTENT(OUT):: znotm
  61     REAL             :: yz, y1, y2, y3, y4
  62
  63     yz =  0.0001344
  64     y1 =  3.015e-05
  65     y2 =  1.517e-06
  66     y3 = -3.567e-08
  67     y4 =  2.046e-10
  68
  69     IF ( uref .LT. 12.5 ) THEN
  70        znotm  = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
  71     ELSE IF ( uref .GE. 12.5 .AND. uref .LT. 30.0 ) THEN
  72        znotm = (0.0739793 * uref -0.58)/1000.0
  73     ELSE
  74        znotm = yz + uref*y1 + uref**2*y2 + uref**3*y3 + uref**4*y4
  75     END IF
  76
  77   END SUBROUTINE znot_m_v0
  78
  79
  80   SUBROUTINE znot_t_v1(uref,znott)
  81   IMPLICIT NONE
  82
  83 ! uref(m/s)   :   Reference level wind
  84 ! znott(meter):   Roughness scale for temperature/moisture
  85 ! Author      :  Biju Thomas on 02/07/2014
  86
  87     REAL, INTENT(IN) :: uref
  88     REAL, INTENT(OUT):: znott
  89     REAL             :: to0, to1, to2, to3
  90     REAL             :: tr0, tr1, tr2, tr3
  91     REAL             :: tn0, tn1, tn2, tn3, tn4, tn5
  92     REAL             :: ta0, ta1, ta2, ta3, ta4, ta5, ta6
  93     REAL             :: tt0, tt1, tt2, tt3, tt4, tt5, tt6, tt7
  94
  95
  96     tr0 = 6.451939325286488e-08
  97     tr1 = -7.306388137342143e-07
  98     tr2 = -1.3709065148333262e-05
  99     tr3 = 0.00019109962089098182
 100
 101     to0 = 1.4379320027061375e-08
 102     to1 = -2.0674525898850674e-07
 103     to2 = -6.8950970846611e-06
 104     to3 = 0.00012199648268521026
 105
 106     tn0 = 1.4023940955902878e-10
 107     tn1 = -1.4752557214976321e-08
 108     tn2 = 5.90998487691812e-07
 109     tn3 = -1.0920804077770066e-05
 110     tn4 = 8.898205876940546e-05
 111     tn5 = -0.00021123340439418298
 112
 113     tt0 = 1.92409564131838e-12
 114     tt1 = -5.765467086754962e-10
 115     tt2 = 7.276979099726975e-08
 116     tt3 = -5.002261599293387e-06
 117     tt4 = 0.00020220445539973736
 118     tt5 = -0.0048088230565883
 119     tt6 = 0.0623468551971189
 120     tt7 = -0.34019193746967424
 121
 122     ta0 = -1.7787470700719361e-10
 123     ta1 = 4.4691736529848764e-08
 124     ta2 = -3.0261975348463414e-06
 125     ta3 = -0.00011680322286017206
 126     ta4 = 0.024449377821884846
 127     ta5 = -1.1228628619105638
 128     ta6 = 17.358026773905973
 129
 130     IF ( uref .LE. 7.0 ) THEN
 131       znott = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
 132     ELSEIF ( uref  .GE. 7.0 .AND. uref .LT. 12.5 ) THEN
 133       znott =  tr3 + tr2*uref + tr1*uref**2 + tr0*uref**3
 134     ELSEIF ( uref  .GE. 12.5 .AND. uref .LT. 15.0 ) THEN
 135       znott =  to3 + to2*uref + to1*uref**2 + to0*uref**3
 136     ELSEIF ( uref .GE. 15.0  .AND. uref .LT. 30.0) THEN
 137       znott =  tn5 + tn4*uref + tn3*uref**2 + tn2*uref**3 + tn1*uref**4 +   &
 138                                                        tn0*uref**5
 139     ELSEIF ( uref .GE. 30.0  .AND. uref .LT. 60.0) THEN
 140       znott = tt7 + tt6*uref + tt5*uref**2  + tt4*uref**3 + tt3*uref**4 +   &
 141              tt2*uref**5 + tt1*uref**6 + tt0*uref**7
 142     ELSE
 143       znott =  ta6 + ta5*uref + ta4*uref**2  + ta3*uref**3 + ta2*uref**4 +  &
 144               ta1*uref**5 + ta0*uref**6
 145     END IF
 146
 147   END SUBROUTINE znot_t_v1
 148
 149   SUBROUTINE znot_t_v0(uref,znott)
 150   IMPLICIT NONE
 151
 152 ! uref(m/s)   :   Reference level wind
 153 ! znott(meter):   Roughness scale for temperature/moisture
 154 ! Author      :  Biju Thomas on 02/07/2014
 155
 156     REAL, INTENT(IN) :: uref
 157     REAL, INTENT(OUT):: znott
 158
 159     IF ( uref .LT. 7.0 ) THEN
 160        znott = (0.0185 / 9.8*(7.59e-4*uref**2+2.46e-2*uref)**2)
 161     ELSE
 162        znott = (0.2375*exp(-0.5250*uref) + 0.0025*exp(-0.0211*uref))*0.01
 163     END IF
 164
 165   END SUBROUTINE znot_t_v0
 166
 167
 168   SUBROUTINE znot_t_v2(uu,znott)
 169   IMPLICIT NONE
 170
 171 ! uu in MKS
 172 ! znott in m
 173 ! Biju Thomas on 02/12/2015
 174 !
 175
 176     REAL, INTENT(IN) :: uu
 177     REAL, INTENT(OUT):: znott
 178     REAL             :: ta0, ta1, ta2, ta3, ta4, ta5, ta6
 179     REAL             :: tb0, tb1, tb2, tb3, tb4, tb5, tb6
 180     REAL             :: tt0, tt1, tt2, tt3, tt4, tt5, tt6
 181
 182     ta0 = 2.51715926619e-09
 183     ta1 = -1.66917514012e-07
 184     ta2 = 4.57345863551e-06
 185     ta3 = -6.64883696932e-05
 186     ta4 = 0.00054390175125
 187     ta5 = -0.00239645231325
 188     ta6 = 0.00453024927761
 189
 190
 191     tb0 = -1.72935914649e-14
 192     tb1 = 2.50587455802e-12
 193     tb2 = -7.90109676541e-11
 194     tb3 = -4.40976353607e-09
 195     tb4 = 3.68968179733e-07
 196     tb5 = -9.43728336756e-06
 197     tb6 = 8.90731312383e-05
 198
 199     tt0 = 4.68042680888e-14
 200     tt1 = -1.98125754931e-11
 201     tt2 = 3.41357133496e-09
 202     tt3 = -3.05130605309e-07
 203     tt4 = 1.48243563819e-05
 204     tt5 = -0.000367207751936
 205     tt6 = 0.00357204479347
 206
 207     IF ( uu .LE. 7.0 ) THEN
 208        znott = (0.0185 / 9.8*(7.59e-4*uu**2+2.46e-2*uu)**2)
 209     ELSEIF ( uu  .GE. 7.0 .AND. uu .LT. 15. ) THEN
 210        znott = ta6 + ta5*uu + ta4*uu**2  + ta3*uu**3 + ta2*uu**4 +     &
 211                ta1*uu**5 + ta0*uu**6
 212     ELSEIF ( uu .GE. 15.0  .AND. uu .LT. 60.0) THEN
 213        znott = tb6 + tb5*uu + tb4*uu**2 + tb3*uu**3 + tb2*uu**4 +      &
 214                tb1*uu**5 + tb0*uu**6
 215     ELSE
 216        znott = tt6 + tt5*uu + tt4*uu**2  + tt3*uu**3 + tt2*uu**4 +    &
 217                tt1*uu**5 + tt0*uu**6
 218     END IF
 219
 220   END SUBROUTINE znot_t_v2
 221
 222   SUBROUTINE znot_m_v6(uref,znotm)
 223   IMPLICIT NONE
 224 ! Calculate areodynamical roughness over water with input 10-m wind
 225 ! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
 226 ! For high winds, try to fit available observational data
 227 !
 228 ! Bin Liu, NOAA/NCEP/EMC 2017
 229 !
 230 ! uref(m/s)   :   wind speed at 10-m height
 231 ! znotm(meter):   areodynamical roughness scale over water
 232 !
 233
 234     REAL, INTENT(IN) :: uref
 235     REAL, INTENT(OUT):: znotm
 236     REAL             :: p13, p12, p11, p10
 237     REAL             :: p25, p24, p23, p22, p21, p20
 238     REAL             :: p35, p34, p33, p32, p31, p30
 239     REAL             :: p40
 240
 241     p13 = -1.296521881682694e-02
 242     p12 =  2.855780863283819e-01
 243     p11 = -1.597898515251717e+00
 244     p10 = -8.396975715683501e+00
 245
 246     p25 =  3.790846746036765e-10
 247     p24 =  3.281964357650687e-09
 248     p23 =  1.962282433562894e-07
 249     p22 = -1.240239171056262e-06
 250     p21 =  1.739759082358234e-07
 251     p20 =  2.147264020369413e-05
 252
 253     p35 =  1.840430200185075e-07
 254     p34 = -2.793849676757154e-05
 255     p33 =  1.735308193700643e-03
 256     p32 = -6.139315534216305e-02
 257     p31 =  1.255457892775006e+00
 258     p30 = -1.663993561652530e+01
 259
 260     p40 =  4.579369142033410e-04
 261
 262     if (uref >= 0.0 .and.  uref <= 6.5 ) then
 263       znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
 264     elseif (uref > 6.5 .and. uref <= 15.7) then
 265       znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 266     elseif (uref > 15.7 .and. uref <= 53.0) then
 267       znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
 268     elseif ( uref > 53.0) then
 269       znotm = p40
 270     else
 271       print*, 'Wrong input uref value:',uref
 272     endif
 273
 274   END SUBROUTINE znot_m_v6
 275
 276   SUBROUTINE znot_t_v6(uref,znott)
 277   IMPLICIT NONE
 278 ! Calculate scalar roughness over water with input 10-m wind
 279 ! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
 280 ! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
 281 !
 282 ! Bin Liu, NOAA/NCEP/EMC 2017
 283 !
 284 ! uref(m/s)   :   wind speed at 10-m height
 285 ! znott(meter):   scalar roughness scale over water
 286 !
 287
 288     REAL, INTENT(IN) :: uref
 289     REAL, INTENT(OUT):: znott
 290
 291     REAL             :: p00
 292     REAL             :: p15, p14, p13, p12, p11, p10
 293     REAL             :: p25, p24, p23, p22, p21, p20
 294     REAL             :: p35, p34, p33, p32, p31, p30
 295     REAL             :: p45, p44, p43, p42, p41, p40
 296     REAL             :: p56, p55, p54, p53, p52, p51, p50
 297     REAL             :: p60
 298
 299     p00 =  1.100000000000000e-04
 300
 301     p15 = -9.144581627678278e-10
 302     p14 =  7.020346616456421e-08
 303     p13 = -2.155602086883837e-06
 304     p12 =  3.333848806567684e-05
 305     p11 = -2.628501274963990e-04
 306     p10 =  8.634221567969181e-04
 307
 308     p25 = -8.654513012535990e-12
 309     p24 =  1.232380050058077e-09
 310     p23 = -6.837922749505057e-08
 311     p22 =  1.871407733439947e-06
 312     p21 = -2.552246987137160e-05
 313     p20 =  1.428968311457630e-04
 314
 315     p35 =  3.207515102100162e-12
 316     p34 = -2.945761895342535e-10
 317     p33 =  8.788972147364181e-09
 318     p32 = -3.814457439412957e-08
 319     p31 = -2.448983648874671e-06
 320     p30 =  3.436721779020359e-05
 321
 322     p45 = -3.530687797132211e-11
 323     p44 =  3.939867958963747e-09
 324     p43 = -1.227668406985956e-08
 325     p42 = -1.367469811838390e-05
 326     p41 =  5.988240863928883e-04
 327     p40 = -7.746288511324971e-03
 328
 329     p56 = -1.187982453329086e-13
 330     p55 =  4.801984186231693e-11
 331     p54 = -8.049200462388188e-09
 332     p53 =  7.169872601310186e-07
 333     p52 = -3.581694433758150e-05
 334     p51 =  9.503919224192534e-04
 335     p50 = -1.036679430885215e-02
 336
 337     p60 =  4.751256171799112e-05
 338
 339     if (uref >= 0.0 .and. uref < 5.9 ) then
 340       znott = p00
 341     elseif (uref >= 5.9 .and. uref <= 15.4) then
 342       znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
 343     elseif (uref > 15.4 .and. uref <= 21.6) then
 344       znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 345     elseif (uref > 21.6 .and. uref <= 42.2) then
 346       znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
 347     elseif ( uref > 42.2 .and. uref <= 53.3) then
 348       znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
 349     elseif ( uref > 53.3 .and. uref <= 80.0) then
 350       znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
 351     elseif ( uref > 80.0) then
 352       znott = p60
 353     else
 354       print*, 'Wrong input uref value:',uref
 355     endif
 356
 357   END SUBROUTINE znot_t_v6
 358
 359   SUBROUTINE znot_m_v7(uref,znotm)
 360   IMPLICIT NONE
 361 ! Calculate areodynamical roughness over water with input 10-m wind
 362 ! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
 363 ! For high winds, try to fit available observational data
 364 ! Comparing to znot_t_v6, slightly decrease Cd for higher wind speed
 365 !
 366 ! Bin Liu, NOAA/NCEP/EMC 2018
 367 !
 368 ! uref(m/s)   :   wind speed at 10-m height
 369 ! znotm(meter):   areodynamical roughness scale over water
 370 !
 371
 372     REAL, INTENT(IN) :: uref
 373     REAL, INTENT(OUT):: znotm
 374     REAL             :: p13, p12, p11, p10
 375     REAL             :: p25, p24, p23, p22, p21, p20
 376     REAL             :: p35, p34, p33, p32, p31, p30
 377     REAL             :: p40
 378
 379     p13 = -1.296521881682694e-02
 380     p12 =  2.855780863283819e-01
 381     p11 = -1.597898515251717e+00
 382     p10 = -8.396975715683501e+00
 383
 384     p25 =  3.790846746036765e-10
 385     p24 =  3.281964357650687e-09
 386     p23 =  1.962282433562894e-07
 387     p22 = -1.240239171056262e-06
 388     p21 =  1.739759082358234e-07
 389     p20 =  2.147264020369413e-05
 390
 391     p35 =  1.897534489606422e-07
 392     p34 = -3.019495980684978e-05
 393     p33 =  1.931392924987349e-03
 394     p32 = -6.797293095862357e-02
 395     p31 =  1.346757797103756e+00
 396     p30 = -1.707846930193362e+01
 397
 398     p40 =  3.371427455376717e-04
 399
 400     if (uref >= 0.0 .and.  uref <= 6.5 ) then
 401       znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
 402     elseif (uref > 6.5 .and. uref <= 15.7) then
 403       znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 404     elseif (uref > 15.7 .and. uref <= 53.0) then
 405       znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
 406     elseif ( uref > 53.0) then
 407       znotm = p40
 408     else
 409       print*, 'Wrong input uref value:',uref
 410     endif
 411
 412   END SUBROUTINE znot_m_v7
 413
 414   SUBROUTINE znot_t_v7(uref,znott)
 415   IMPLICIT NONE
 416 ! Calculate scalar roughness over water with input 10-m wind
 417 ! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
 418 ! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
 419 ! To be compatible with the slightly decreased Cd for higher wind speed
 420 !
 421 ! Bin Liu, NOAA/NCEP/EMC 2018
 422 !
 423 ! uref(m/s)   :   wind speed at 10-m height
 424 ! znott(meter):   scalar roughness scale over water
 425 !
 426
 427     REAL, INTENT(IN) :: uref
 428     REAL, INTENT(OUT):: znott
 429
 430     REAL             :: p00
 431     REAL             :: p15, p14, p13, p12, p11, p10
 432     REAL             :: p25, p24, p23, p22, p21, p20
 433     REAL             :: p35, p34, p33, p32, p31, p30
 434     REAL             :: p45, p44, p43, p42, p41, p40
 435     REAL             :: p56, p55, p54, p53, p52, p51, p50
 436     REAL             :: p60
 437
 438     p00 =  1.100000000000000e-04
 439
 440     p15 = -9.193764479895316e-10
 441     p14 =  7.052217518653943e-08
 442     p13 = -2.163419217747114e-06
 443     p12 =  3.342963077911962e-05
 444     p11 = -2.633566691328004e-04
 445     p10 =  8.644979973037803e-04
 446
 447     p25 = -9.402722450219142e-12
 448     p24 =  1.325396583616614e-09
 449     p23 = -7.299148051141852e-08
 450     p22 =  1.982901461144764e-06
 451     p21 = -2.680293455916390e-05
 452     p20 =  1.484341646128200e-04
 453
 454     p35 =  7.921446674311864e-12
 455     p34 = -1.019028029546602e-09
 456     p33 =  5.251986927351103e-08
 457     p32 = -1.337841892062716e-06
 458     p31 =  1.659454106237737e-05
 459     p30 = -7.558911792344770e-05
 460
 461     p45 = -2.694370426850801e-10
 462     p44 =  5.817362913967911e-08
 463     p43 = -5.000813324746342e-06
 464     p42 =  2.143803523428029e-04
 465     p41 = -4.588070983722060e-03
 466     p40 =  3.924356617245624e-02
 467
 468     p56 = -1.663918773476178e-13
 469     p55 =  6.724854483077447e-11
 470     p54 = -1.127030176632823e-08
 471     p53 =  1.003683177025925e-06
 472     p52 = -5.012618091180904e-05
 473     p51 =  1.329762020689302e-03
 474     p50 = -1.450062148367566e-02
 475
 476     p60 =  6.840803042788488e-05
 477
 478     if (uref >= 0.0 .and. uref < 5.9 ) then
 479       znott = p00
 480     elseif (uref >= 5.9 .and. uref <= 15.4) then
 481       znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
 482     elseif (uref > 15.4 .and. uref <= 21.6) then
 483       znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 484     elseif (uref > 21.6 .and. uref <= 42.6) then
 485       znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
 486     elseif ( uref > 42.6 .and. uref <= 53.0) then
 487       znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
 488     elseif ( uref > 53.0 .and. uref <= 80.0) then
 489       znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
 490     elseif ( uref > 80.0) then
 491       znott = p60
 492     else
 493       print*, 'Wrong input uref value:',uref
 494     endif
 495
 496   END SUBROUTINE znot_t_v7
 497
 498   SUBROUTINE znot_m_v8(uref,znotm)
 499   IMPLICIT NONE
 500 ! Calculate areodynamical roughness over water with input 10-m wind
 501 ! For low-to-moderate winds, try to match the Cd-U10 relationship from COARE V3.5 (Edson et al. 2013)
 502 ! For high winds, try to fit available observational data
 503 ! Comparing to znot_t_v6, slightly decrease Cd for higher wind speed
 504 ! And this is another variation similar to v7
 505 !
 506 ! Bin Liu, NOAA/NCEP/EMC 2018
 507 !
 508 ! uref(m/s)   :   wind speed at 10-m height
 509 ! znotm(meter):   areodynamical roughness scale over water
 510 !
 511
 512     REAL, INTENT(IN) :: uref
 513     REAL, INTENT(OUT):: znotm
 514     REAL             :: p13, p12, p11, p10
 515     REAL             :: p25, p24, p23, p22, p21, p20
 516     REAL             :: p35, p34, p33, p32, p31, p30
 517     REAL             :: p40
 518
 519     p13 = -1.296521881682694e-02
 520     p12 =  2.855780863283819e-01
 521     p11 = -1.597898515251717e+00
 522     p10 = -8.396975715683501e+00
 523
 524     p25 =  3.790846746036765e-10
 525     p24 =  3.281964357650687e-09
 526     p23 =  1.962282433562894e-07
 527     p22 = -1.240239171056262e-06
 528     p21 =  1.739759082358234e-07
 529     p20 =  2.147264020369413e-05
 530
 531     p35 =  1.897534489606422e-07
 532     p34 = -3.019495980684978e-05
 533     p33 =  1.931392924987349e-03
 534     p32 = -6.797293095862357e-02
 535     p31 =  1.346757797103756e+00
 536     p30 = -1.707846930193362e+01
 537
 538     p40 =  3.886804744928044e-04
 539
 540     if (uref >= 0.0 .and.  uref <= 6.5 ) then
 541       znotm = exp( p10 + p11*uref + p12*uref**2 + p13*uref**3)
 542     elseif (uref > 6.5 .and. uref <= 15.7) then
 543       znotm = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 544     elseif (uref > 15.7 .and. uref <= 51.5) then
 545       znotm = exp( p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30 )
 546     elseif ( uref > 51.5) then
 547       znotm = p40
 548     else
 549       print*, 'Wrong input uref value:',uref
 550     endif
 551
 552   END SUBROUTINE znot_m_v8
 553
 554   SUBROUTINE znot_t_v8(uref,znott)
 555   IMPLICIT NONE
 556 ! Calculate scalar roughness over water with input 10-m wind
 557 ! For low-to-moderate winds, try to match the Ck-U10 relationship from COARE algorithm
 558 ! For high winds, try to retain the Ck-U10 relationship of FY2015 HWRF
 559 ! To be compatible with the slightly decreased Cd for higher wind speed
 560 ! And this is another variation similar to v7
 561 !
 562 ! Bin Liu, NOAA/NCEP/EMC 2018
 563 !
 564 ! uref(m/s)   :   wind speed at 10-m height
 565 ! znott(meter):   scalar roughness scale over water
 566 !
 567
 568     REAL, INTENT(IN) :: uref
 569     REAL, INTENT(OUT):: znott
 570
 571     REAL             :: p00
 572     REAL             :: p15, p14, p13, p12, p11, p10
 573     REAL             :: p25, p24, p23, p22, p21, p20
 574     REAL             :: p35, p34, p33, p32, p31, p30
 575     REAL             :: p45, p44, p43, p42, p41, p40
 576     REAL             :: p56, p55, p54, p53, p52, p51, p50
 577     REAL             :: p60
 578
 579     p00 =  1.100000000000000e-04
 580
 581     p15 = -9.193764479895316e-10
 582     p14 =  7.052217518653943e-08
 583     p13 = -2.163419217747114e-06
 584     p12 =  3.342963077911962e-05
 585     p11 = -2.633566691328004e-04
 586     p10 =  8.644979973037803e-04
 587
 588     p25 = -9.402722450219142e-12
 589     p24 =  1.325396583616614e-09
 590     p23 = -7.299148051141852e-08
 591     p22 =  1.982901461144764e-06
 592     p21 = -2.680293455916390e-05
 593     p20 =  1.484341646128200e-04
 594
 595     p35 =  7.921446674311864e-12
 596     p34 = -1.019028029546602e-09
 597     p33 =  5.251986927351103e-08
 598     p32 = -1.337841892062716e-06
 599     p31 =  1.659454106237737e-05
 600     p30 = -7.558911792344770e-05
 601
 602     p45 = -2.706461188613193e-10
 603     p44 =  5.845859022891930e-08
 604     p43 = -5.027577045502003e-06
 605     p42 =  2.156326523752734e-04
 606     p41 = -4.617267288861201e-03
 607     p40 =  3.951492707214883e-02
 608
 609     p56 = -1.112896580069263e-13
 610     p55 =  4.450334755105140e-11
 611     p54 = -7.375373918500171e-09
 612     p53 =  6.493685149526543e-07
 613     p52 = -3.206421106713471e-05
 614     p51 =  8.407596231678149e-04
 615     p50 = -9.027924333673693e-03
 616
 617     p60 =  5.791179079892191e-05
 618
 619     if (uref >= 0.0 .and. uref < 5.9 ) then
 620       znott = p00
 621     elseif (uref >= 5.9 .and. uref <= 15.4) then
 622       znott = p15*uref**5 + p14*uref**4 + p13*uref**3 + p12*uref**2 + p11*uref + p10
 623     elseif (uref > 15.4 .and. uref <= 21.6) then
 624       znott = p25*uref**5 + p24*uref**4 + p23*uref**3 + p22*uref**2 + p21*uref + p20
 625     elseif (uref > 21.6 .and. uref <= 42.6) then
 626       znott = p35*uref**5 + p34*uref**4 + p33*uref**3 + p32*uref**2 + p31*uref + p30
 627     elseif ( uref > 42.6 .and. uref <= 51.5) then
 628       znott = p45*uref**5 + p44*uref**4 + p43*uref**3 + p42*uref**2 + p41*uref + p40
 629     elseif ( uref > 51.5 .and. uref <= 80.0) then
 630       znott = p56*uref**6 + p55*uref**5 + p54*uref**4 + p53*uref**3 + p52*uref**2 + p51*uref + p50
 631     elseif ( uref > 80.0) then
 632       znott = p60
 633     else
 634       print*, 'Wrong input uref value:',uref
 635     endif
 636
 637   END SUBROUTINE znot_t_v8
 638
 639    SUBROUTINE znot_wind10m(w10m,znott,znotm,icoef_sf)
 640    IMPLICIT NONE
 641
 642 ! w10m(m/s)   :   10-m wind speed
 643 ! znott(meter):   Roughness scale for temperature/moisture, zt
 644 ! znotm(meter):   Roughness scale for momentum, z0
 645 ! Author      :  Weiguo Wang on 02/24/2016
 646 !            convert from icoef=0,1,2 to have 10m level cd, ch match obs
 647      REAL, INTENT(IN) :: w10m
 648      INTEGER, INTENT(IN) :: icoef_sf
 649      REAL, INTENT(OUT):: znott, znotm
 650
 651      real :: zm,zt,windmks, zlev,z10, tmp, zlevt, aaa, zm1,zt1
 652         zlev=20.0
 653         zlevt=10.0
 654         z10=10.0
 655             windmks=w10m
 656             if (windmks > 85.0) windmks=85.0
 657             if (windmks < 1.0) windmks=1.0
 658             if ( icoef_sf .EQ. 1) then
 659               call  znot_m_v1(windmks,zm1)
 660               call  znot_t_v1(windmks,zt1)
 661
 662             else if ( icoef_sf .EQ. 0 ) then
 663               call  znot_m_v0(windmks,zm1)
 664               call  znot_t_v0(windmks,zt1)
 665
 666             else  if( icoef_sf .EQ. 2 ) then
 667               call  znot_m_v1(windmks,zm1)
 668               call  znot_t_v2(windmks,zt1)
 669
 670             else  if( icoef_sf .EQ. 3 ) then
 671               call  znot_m_v1(windmks,zm)
 672               call  znot_t_v2(windmks,zt)
 673 !! adjust a little to match obs at 10m, cd is reduced
 674             tmp=0.4*0.4/(alog(zlev/zm))**2   ! cd at zlev
 675             zm1=z10/exp( sqrt(0.4*0.4/(tmp*0.95-0.0002)) )
 676 !ch
 677             tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
 678             zt1=z10/exp( 0.4*0.4/( 0.95*tmp*alog(z10/zm1) )  )
 679
 680             else if( icoef_sf .EQ. 4 ) then
 681
 682               call  znot_m_v1(windmks,zm)
 683               call  znot_t_v2(windmks,zt)
 684 !!  for wind<20, cd similar to icoef=2 at 10m, then reduced
 685              tmp=0.4*0.4/(alog(10.0/zm))**2   ! cd at zlev
 686              aaa=0.75
 687             if (windmks < 20) then
 688               aaa=0.99
 689             elseif(windmks < 45.0) then
 690               aaa=0.99+(windmks-20)*(0.75-0.99)/(45.0-20.0)
 691             endif
 692             zm1=z10/exp( sqrt(0.4*0.4/(tmp*aaa)) )
 693 !ch
 694           tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
 695             zt1=z10/exp( 0.4*0.4/( 0.95*tmp*alog(z10/zm1) )  )
 696
 697             else if( icoef_sf .EQ. 5 ) then
 698
 699               call  znot_m_v1(windmks,zm)
 700               call  znot_t_v2(windmks,zt)
 701 !!  for wind<20, cd similar to icoef=2 at 10m, then reduced
 702              tmp=0.4*0.4/(alog(10.0/zm))**2   ! cd at zlev
 703              aaa=0.80
 704             if (windmks < 20) then
 705               aaa=1.0
 706             elseif(windmks < 45.0) then
 707               aaa=1.0+(windmks-20)*(0.80-1.0)/(45.0-20.0)
 708             endif
 709             zm1=z10/exp( sqrt(0.4*0.4/(tmp*aaa)) )
 710 !ch
 711           tmp=0.4*0.4/(alog(zlevt/zm)*alog(zlevt/zt))  ! ch at zlev using old formula
 712             zt1=z10/exp( 0.4*0.4/( 1.0*tmp*alog(z10/zm1) )  )
 713
 714             else  if( icoef_sf .EQ. 6 ) then
 715               call  znot_m_v6(windmks,zm1)
 716               call  znot_t_v6(windmks,zt1)
 717             else  if( icoef_sf .EQ. 7 ) then
 718               call  znot_m_v7(windmks,zm1)
 719               call  znot_t_v7(windmks,zt1)
 720             else  if( icoef_sf .EQ. 8 ) then
 721               call  znot_m_v8(windmks,zm1)
 722               call  znot_t_v8(windmks,zt1)
 723            else
 724              CALL wrf_error_fatal('icoef_sf must be one of 0,1,2,3,4,5,6,7,8')
 725
 726           endif
 727           znott=zt1
 728           znotm=zm1
 729
 730   end subroutine znot_wind10m
 731
 732 END MODULE module_sf_exchcoef
 733