18 use data_override_mod,
only : data_override_init, data_override
20 implicit none ;
private
22 #include <MOM_memory.h>
50 logical,
public :: usewaves
51 logical,
public :: lagrangianmixing
64 integer,
public :: stklevelmode=1
71 real,
allocatable,
dimension(:),
public :: &
73 real,
allocatable,
dimension(:),
public :: &
75 real,
allocatable,
dimension(:),
public :: &
77 real,
allocatable,
dimension(:),
public :: &
79 real,
allocatable,
dimension(:,:,:),
public :: &
83 real,
allocatable,
dimension(:,:,:),
public :: &
87 real,
allocatable,
dimension(:,:),
public :: &
88 la_sl,& !< SL Langmuir number (directionality factored later)
92 real,
allocatable,
dimension(:,:),
public :: &
95 real,
allocatable,
dimension(:,:),
public :: &
98 real,
allocatable,
dimension(:,:,:),
public :: &
102 real,
allocatable,
dimension(:,:,:),
public :: &
106 real,
allocatable,
dimension(:,:,:),
public :: &
110 type(time_type),
pointer,
public :: time
119 real :: la_min = 0.05
122 integer,
public :: id_surfacestokes_x = -1 , id_surfacestokes_y = -1
123 integer,
public :: id_3dstokes_x = -1 , id_3dstokes_y = -1
124 integer,
public :: id_la_turb = -1
169 #include "version_variable.h"
171 character(len=40) ::
mdl =
"MOM_wave_interface"
193 type(time_type),
target,
intent(in) :: time
199 type(
diag_ctrl),
target,
intent(inout) :: diag
202 character*(13) :: tmpstring1,tmpstring2
203 character*(5),
parameter :: null_string =
"EMPTY"
204 character*(12),
parameter :: testprof_string =
"TEST_PROFILE"
205 character*(13),
parameter :: surfbands_string =
"SURFACE_BANDS"
206 character*(5),
parameter :: dhh85_string =
"DHH85"
207 character*(4),
parameter :: lf17_string =
"LF17"
208 character*(12),
parameter :: dataovr_string =
"DATAOVERRIDE"
209 character*(7),
parameter :: coupler_string =
"COUPLER"
210 character*(5),
parameter :: input_string =
"INPUT"
213 if (
associated(cs))
then
214 call mom_error(fatal,
"wave_interface_init called with an associated"//&
215 "control structure.")
237 call get_param(param_file,
mdl,
"LAGRANGIAN_MIXING", cs%LagrangianMixing, &
238 "Flag to use Lagrangian Mixing of momentum", units=
"", &
240 if (cs%LagrangianMixing)
then
242 call mom_error(fatal,
"Should you be enabling Lagrangian Mixing? Code not ready.")
244 call get_param(param_file,
mdl,
"STOKES_MIXING", cs%StokesMixing, &
245 "Flag to use Stokes Mixing of momentum", units=
"", &
247 if (cs%StokesMixing)
then
249 call mom_error(fatal,
"Should you be enabling Stokes Mixing? Code not ready.")
251 call get_param(param_file,
mdl,
"CORIOLIS_STOKES", cs%CoriolisStokes, &
252 "Flag to use Coriolis Stokes acceleration", units=
"", &
254 if (cs%CoriolisStokes)
then
256 call mom_error(fatal,
"Should you be enabling Coriolis-Stokes? Code not ready.")
260 call get_param(param_file,
mdl,
"WAVE_METHOD",tmpstring1, &
261 "Choice of wave method, valid options include: \n"// &
262 " TEST_PROFILE - Prescribed from surface Stokes drift \n"// &
263 " and a decay wavelength.\n"// &
264 " SURFACE_BANDS - Computed from multiple surface values \n"// &
265 " and decay wavelengths.\n"// &
266 " DHH85 - Uses Donelan et al. 1985 empirical \n"// &
267 " wave spectrum with prescribed values. \n"// &
268 " LF17 - Infers Stokes drift profile from wind \n"// &
269 " speed following Li and Fox-Kemper 2017.\n", &
270 units=
'', default=null_string)
271 select case (trim(tmpstring1))
273 call mom_error(fatal,
"wave_interface_init called with no specified "//&
275 case (testprof_string)
278 'Surface Stokes (x) for test profile',&
279 units=
'm/s',default=0.1)
281 'Surface Stokes (y) for test profile',&
282 units=
'm/s',default=0.0)
284 units=
'm', default=50.0, scale=us%m_to_Z)
285 case (surfbands_string)
287 call get_param(param_file,
mdl,
"SURFBAND_SOURCE",tmpstring2, &
288 "Choice of SURFACE_BANDS data mode, valid options include: \n"// &
289 " DATAOVERRIDE - Read from NetCDF using FMS DataOverride. \n"// &
290 " COUPLER - Look for variables from coupler pass \n"// &
291 " INPUT - Testing with fixed values.", &
292 units=
'', default=null_string)
293 select case (trim(tmpstring2))
295 call mom_error(fatal,
"wave_interface_init called with SURFACE_BANDS"//&
296 " but no SURFBAND_SOURCE.")
297 case (dataovr_string)
300 "Filename of surface Stokes drift input band data.", default=
"StkSpec.nc")
301 case (coupler_string)
306 "Prescribe number of wavenumber bands for Stokes drift. "// &
307 "Make sure this is consistnet w/ WAVENUMBERS, STOKES_X, and "// &
308 "STOKES_Y, there are no safety checks in the code.", &
310 allocate( cs%WaveNum_Cen(1:
numbands) )
311 cs%WaveNum_Cen(:) = 0.0
312 allocate( cs%PrescribedSurfStkX(1:
numbands))
313 cs%PrescribedSurfStkX(:) = 0.0
314 allocate( cs%PrescribedSurfStkY(1:
numbands))
315 cs%PrescribedSurfStkY(:) = 0.0
316 allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,1:
numbands))
317 cs%STKx0(:,:,:) = 0.0
318 allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,1:
numbands))
319 cs%STKy0(:,:,:) = 0.0
321 call get_param(param_file,
mdl,
"SURFBAND_WAVENUMBERS",cs%WaveNum_Cen, &
322 "Central wavenumbers for surface Stokes drift bands.",units=
'rad/m', &
324 call get_param(param_file,
mdl,
"SURFBAND_STOKES_X",cs%PrescribedSurfStkX, &
325 "X-direction surface Stokes drift for bands.",units=
'm/s', &
327 call get_param(param_file,
mdl,
"SURFBAND_STOKES_Y",cs%PrescribedSurfStkY, &
328 "Y-direction surface Stokes drift for bands.",units=
'm/s', &
331 call mom_error(fatal,
'Check WAVE_METHOD.')
336 call mom_error(warning,
"DHH85 only ever set-up for uniform cases w/"//&
337 " Stokes drift in x-direction.")
339 "Choose true to use waveage in peak frequency.", &
340 units=
'', default=.false.)
342 "Wave Age for DHH85 spectrum.", &
343 units=
'', default=1.2)
345 "Wind speed for DHH85 spectrum.", &
346 units=
'', default=10.0)
348 "Flag to disable updating DHH85 Stokes drift.", &
353 call mom_error(fatal,
'Check WAVE_METHOD.')
358 "The depth (normalized by BLD) to average Stokes drift over in "//&
359 "Langmuir number calculation, where La = sqrt(ust/Stokes).", &
360 units=
"nondim",default=0.04)
362 "Flag (logical) if using misalignment bt shear and waves in LA",&
364 call get_param(param_file,
mdl,
"MIN_LANGMUIR", cs%La_min, &
365 "A minimum value for all Langmuir numbers that is not physical, "//&
366 "but is likely only encountered when the wind is very small and "//&
367 "therefore its effects should be mostly benign.",units=
"nondim",&
372 allocate(cs%Us_x(g%isdB:g%IedB,g%jsd:g%jed,g%ke))
374 allocate(cs%Us_y(g%isd:g%Ied,g%jsdB:g%jedB,g%ke))
377 allocate(cs%US0_x(g%isdB:g%iedB,g%jsd:g%jed))
379 allocate(cs%US0_y(g%isd:g%ied,g%jsdB:g%jedB))
382 allocate(cs%La_SL(g%isc:g%iec,g%jsc:g%jec))
383 allocate(cs%La_turb(g%isc:g%iec,g%jsc:g%jec))
385 cs%La_turb (:,:) = 0.0
387 if (cs%StokesMixing)
then
388 allocate(cs%KvS(g%isd:g%Ied,g%jsd:g%jed,g%ke))
393 cs%id_surfacestokes_y = register_diag_field(
'ocean_model',
'surface_stokes_y', &
394 cs%diag%axesCu1,time,
'Surface Stokes drift (y)',
'm s-1')
395 cs%id_surfacestokes_x = register_diag_field(
'ocean_model',
'surface_stokes_x', &
396 cs%diag%axesCv1,time,
'Surface Stokes drift (x)',
'm s-1')
397 cs%id_3dstokes_y = register_diag_field(
'ocean_model',
'3d_stokes_y', &
398 cs%diag%axesCvL,time,
'3d Stokes drift (y)',
'm s-1')
399 cs%id_3dstokes_x = register_diag_field(
'ocean_model',
'3d_stokes_x', &
400 cs%diag%axesCuL,time,
'3d Stokes drift (y)',
'm s-1')
401 cs%id_La_turb = register_diag_field(
'ocean_model',
'La_turbulent',&
402 cs%diag%axesT1,time,
'Surface (turbulent) Langmuir number',
'nondim')
411 character*(5),
parameter :: null_string =
"EMPTY"
412 character*(4),
parameter :: lf17_string =
"LF17"
413 character*(13) :: tmpstring1
414 logical :: statisticalwaves
418 "The depth (normalized by BLD) to average Stokes drift over in "//&
419 "Langmuir number calculation, where La = sqrt(ust/Stokes).", &
420 units=
"nondim",default=0.04)
423 call get_param(param_file,
mdl,
"USE_LA_LI2016",statisticalwaves, &
424 do_not_log=.true.,default=.false.)
425 if (statisticalwaves)
then
441 type(time_type),
intent(in) :: day
442 type(time_type),
intent(in) :: dt
444 integer :: ii, jj, kk, b
445 type(time_type) :: day_center
448 day_center = day + dt/2
461 cs%STKx0(ii,jj,b) = cs%PrescribedSurfStkX(b)
466 cs%STKY0(ii,jj,b) = cs%PrescribedSurfStkY(b)
483 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
485 real,
dimension(SZI_(G),SZJ_(G)), &
488 real :: top, midpoint, bottom, one_cm
490 real :: cmn_fac, wn, ustokes
492 integer :: ii, jj, kk, b, iim1, jjm1
494 one_cm = 0.01*us%m_to_Z
500 do ii = g%isdB,g%iedB
507 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
508 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
509 cs%Us_x(ii,jj,kk) =
tp_stkx0*exp(midpoint*decayscale)
514 do jj = g%jsdB,g%jedB
520 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
521 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
522 cs%Us_y(ii,jj,kk) =
tp_stky0*exp(midpoint*decayscale)
535 do ii = g%isdB,g%iedB
542 cmn_fac = (1.0-exp(-one_cm*2.*cs%WaveNum_Cen(b))) / &
543 (one_cm*2.*cs%WaveNum_Cen(b))
548 cs%US0_x(ii,jj) = cs%US0_x(ii,jj) + cs%STKx0(ii,jj,b)*cmn_fac
555 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
556 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
560 cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b))-exp(bottom*2.*cs%WaveNum_Cen(b)))&
561 / ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
563 if (cs%StkLevelMode==0)
then
565 cmn_fac = exp(midpoint*2.*(2.*
pi*cs%Freq_Cen(b)*us%T_to_s)**2/(us%L_to_Z**2*gv%g_Earth))
566 elseif (cs%StkLevelMode==1)
then
569 wn = (2.*
pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth)
570 cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
573 cs%US_x(ii,jj,kk) = cs%US_x(ii,jj,kk) + cs%STKx0(ii,jj,b)*cmn_fac
580 do jj = g%jsdB,g%jedB
585 cmn_fac = (1.0-exp(-one_cm*2.*cs%WaveNum_Cen(b))) / &
586 (one_cm*2.*cs%WaveNum_Cen(b))
591 cs%US0_y(ii,jj) = cs%US0_y(ii,jj) + cs%STKy0(ii,jj,b)*cmn_fac
598 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
599 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
603 cmn_fac = (exp(top*2.*cs%WaveNum_Cen(b)) - &
604 exp(bottom*2.*cs%WaveNum_Cen(b))) / &
605 ((top-bottom)*(2.*cs%WaveNum_Cen(b)))
607 if (cs%StkLevelMode==0)
then
609 cmn_fac = exp(midpoint*2.*(2.*
pi*cs%Freq_Cen(b)*us%T_to_s)**2/(us%L_to_Z**2*gv%g_Earth))
610 elseif (cs%StkLevelMode==1)
then
613 wn = (2.*
pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth)
614 cmn_fac = (exp(2.*wn*top)-exp(2.*wn*bottom)) / (2.*wn*(top-bottom))
617 cs%US_y(ii,jj,kk) = cs%US_y(ii,jj,kk) + cs%STKy0(ii,jj,b)*cmn_fac
624 do ii = g%isdB,g%iedB
630 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(iim1,jj,kk))
631 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(iim1,jj,kk))
636 call dhh85_mid(gv, us, midpoint, ustokes)
638 cs%US_x(ii,jj,kk) = ustokes
643 do jj = g%jsdB,g%jedB
648 midpoint = bottom - gv%H_to_Z*0.25*(h(ii,jj,kk)+h(ii,jjm1,kk))
649 bottom = bottom - gv%H_to_Z*0.5*(h(ii,jj,kk)+h(ii,jjm1,kk))
656 cs%US_y(ii,jj,kk) = 0.0
667 do ii = g%isdB,g%iedB
669 cs%Us_x(ii,jj,kk) = 0.
673 do jj = g%jsdB,g%jedB
674 cs%Us_y(ii,jj,kk) = 0.
685 top = h(ii,jj,1)*gv%H_to_Z
687 h(ii,jj,:),override_ma=.false.,waves=cs)
688 cs%La_turb(ii,jj) = la
693 if (cs%id_surfacestokes_y>0) &
694 call post_data(cs%id_surfacestokes_y, cs%us0_y, cs%diag)
695 if (cs%id_surfacestokes_x>0) &
696 call post_data(cs%id_surfacestokes_x, cs%us0_x, cs%diag)
697 if (cs%id_3dstokes_y>0) &
698 call post_data(cs%id_3dstokes_y, cs%us_y, cs%diag)
699 if (cs%id_3dstokes_x>0) &
700 call post_data(cs%id_3dstokes_x, cs%us_x, cs%diag)
701 if (cs%id_La_turb>0) &
702 call post_data(cs%id_La_turb, cs%La_turb, cs%diag)
710 type(time_type),
intent(in) :: day_center
716 real :: temp_x(SZI_(G),SZJ_(G))
717 real :: temp_y(SZI_(G),SZJ_(G))
718 real :: Top, MidPoint
721 integer,
dimension(4) :: start, counter, dims, dim_id
722 character(len=12) :: dim_name(4)
723 character(20) :: varname, varread1, varread2
724 integer :: rcode_fr, rcode_wn, ncid, varid_fr, varid_wn, id, ndims
727 call data_override_init(ocean_domain_in=g%Domain%mpp_domain)
732 varread1 =
'wavenumber'
733 varread2 =
'frequency'
735 if (rcode_wn /= 0)
then
737 " in MOM_wave_interface.")
741 rcode_wn = nf90_inq_varid(ncid, varread1, varid_wn)
742 rcode_fr = nf90_inq_varid(ncid, varread2, varid_fr)
744 if (rcode_wn /= 0 .and. rcode_fr /= 0)
then
745 call mom_error(fatal,
"error finding variable "//trim(varread1)//&
746 " or "//trim(varread2)//
" in file "//trim(
surfbandfilename)//
" in MOM_wave_interface.")
748 elseif (rcode_wn == 0)
then
751 rcode_wn = nf90_inquire_variable(ncid, varid_wn, ndims=ndims, &
753 if (rcode_wn /= 0)
then
755 'error inquiring dimensions MOM_wave_interface.')
757 rcode_wn = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
758 if (rcode_wn /= 0)
then
759 call mom_error(fatal,
"error reading dimension 1 data for "// &
761 " in MOM_wave_interface.")
763 rcode_wn = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
764 if (rcode_wn /= 0)
then
765 call mom_error(fatal,
"error finding variable "//trim(dim_name(1))//&
769 allocate( cs%WaveNum_Cen(1:id) )
770 cs%WaveNum_Cen(:) = 0.0
771 elseif (rcode_fr == 0)
then
774 rcode_fr = nf90_inquire_variable(ncid, varid_fr, ndims=ndims, &
776 if (rcode_fr /= 0)
then
778 'error inquiring dimensions MOM_wave_interface.')
780 rcode_fr = nf90_inquire_dimension(ncid, dims(1), dim_name(1), len=id)
781 if (rcode_fr /= 0)
then
782 call mom_error(fatal,
"error reading dimension 1 data for "// &
784 " in MOM_wave_interface.")
786 rcode_fr = nf90_inq_varid(ncid, dim_name(1), dim_id(1))
787 if (rcode_fr /= 0)
then
788 call mom_error(fatal,
"error finding variable "//trim(dim_name(1))//&
792 allocate( cs%Freq_Cen(1:id) )
795 allocate( cs%WaveNum_Cen(1:id) )
796 cs%WaveNum_Cen(:) = 0.0
797 allocate( cs%STKx0(g%isdB:g%iedB,g%jsd:g%jed,1:id))
798 cs%STKx0(:,:,:) = 0.0
799 allocate( cs%STKy0(g%isd:g%ied,g%jsdB:g%jedB,1:id))
800 cs%STKy0(:,:,:) = 0.0
808 rcode_wn = nf90_get_var(ncid, dim_id(1), cs%WaveNum_Cen, start, counter)
809 if (rcode_wn /= 0)
then
811 "error reading dimension 1 values for var_name "// &
812 trim(varread1)//
",dim_name "//trim(dim_name(1))// &
816 do b = 1,
numbands ; cs%WaveNum_Cen(b) = us%Z_to_m*cs%WaveNum_Cen(b) ;
enddo
818 rcode_fr = nf90_get_var(ncid, dim_id(1), cs%Freq_Cen, start, counter)
819 if (rcode_fr /= 0)
then
821 "error reading dimension 1 values for var_name "// &
822 trim(varread2)//
",dim_name "//trim(dim_name(1))// &
827 cs%WaveNum_Cen(b) = (2.*
pi*cs%Freq_Cen(b)*us%T_to_s)**2 / (us%L_to_Z**2*gv%g_Earth)
837 write(varname,
"(A3,I0)")
'Usx',b
838 call data_override(
'OCN',trim(varname), temp_x, day_center)
840 write(varname,
'(A3,I0)')
'Usy',b
841 call data_override(
'OCN',trim(varname), temp_y, day_center)
847 if (abs(temp_x(i,j)) > 10. .or. abs(temp_y(i,j)) > 10.)
then
858 cs%STKx0(i,j,b) = 0.5 * (temp_x(i,j) + temp_x(i+1,j))
863 cs%STKy0(i,j,b) = 0.5 * (temp_y(i,j) + temp_y(i,j+1))
878 H, U_H, V_H, Override_MA, Waves )
882 integer,
intent(in) :: i
883 integer,
intent(in) :: j
884 real,
intent(in) :: ustar
885 real,
intent(in) :: hbl
886 logical,
optional,
intent(in) :: override_ma
890 real,
dimension(SZK_(GV)),
optional, &
892 real,
dimension(SZK_(GV)),
optional, &
894 real,
dimension(SZK_(GV)),
optional, &
899 real,
intent(out) :: la
902 real :: top, bottom, midpoint
903 real :: dpt_lasl, sheardirection, wavedirection
904 real :: la_stkx, la_stky, la_stk
905 logical :: continueloop, use_ma
906 real,
dimension(SZK_(G)) :: us_h, vs_h
907 real,
dimension(NumBands) :: stkband_x, stkband_y
911 dpt_lasl = min(-0.1*us%m_to_Z, -
la_frachbl*hbl)
914 if (
present(override_ma)) use_ma = override_ma
918 if (.not.(
present(h).and.
present(u_h).and.
present(v_h)))
then
919 call mom_error(fatal,
'Get_LA_waves requested to consider misalignment.')
921 continueloop = .true.
925 midpoint = bottom + gv%H_to_Z*0.5*h(kk)
926 bottom = bottom + gv%H_to_Z*h(kk)
927 if (midpoint > dpt_lasl .and. kk > 1 .and. continueloop)
then
928 sheardirection = atan2(v_h(1)-v_h(kk),u_h(1)-u_h(kk))
929 continueloop = .false.
936 us_h(kk) = 0.5*(waves%US_X(i,j,kk)+waves%US_X(i-1,j,kk))
937 vs_h(kk) = 0.5*(waves%US_Y(i,j,kk)+waves%US_Y(i,j-1,kk))
941 la_stk = sqrt(la_stkx*la_stkx+la_stky*la_stky)
944 stkband_x(bb) = 0.5*(waves%STKx0(i,j,bb)+waves%STKx0(i-1,j,bb))
945 stkband_y(bb) = 0.5*(waves%STKy0(i,j,bb)+waves%STKy0(i,j-1,bb))
949 la_stk = sqrt(la_stkx**2 + la_stky**2)
953 us_h(kk) = 0.5*(waves%US_X(i,j,kk)+waves%US_X(i-1,j,kk))
954 vs_h(kk) = 0.5*(waves%US_Y(i,j,kk)+waves%US_Y(i,j-1,kk))
958 la_stk = sqrt(la_stkx**2 + la_stky**2)
962 call mom_error(fatal,
"Get_Langmuir_number called without defining a WaveMethod. "//&
963 "Suggest to make sure USE_LT is set/overridden to False or "//&
964 "choose a wave method (or set USE_LA_LI2016 to use statistical "//&
975 la = max(waves%La_min, sqrt(us%Z_to_m*us%s_to_T*ustar / (la_stk+1.e-10)))
979 wavedirection = atan2(la_stky, la_stkx)
980 la = la / sqrt(max(1.e-8, cos( wavedirection - sheardirection)))
1003 real,
intent(in) :: ustar
1004 real,
intent(in) :: hbl
1007 real,
intent(out) :: UStokes_SL
1008 real,
intent(out) :: LA
1011 real,
parameter :: &
1013 u19p5_to_u10 = 1.075, &
1016 fm_into_fp = 1.296, &
1018 us_to_u10 = 0.0162, &
1021 real :: UStokes, hm0, fm, fp, vstokes, kphil, kstar
1022 real :: z0, z0i, r1, r2, r3, r4, tmp, lasl_sqr_i
1025 if (ustar > 0.0)
then
1027 call ust_2_u10_coare3p5(us%Z_to_m*us%s_to_T*ustar*sqrt(us%R_to_kg_m3*gv%Rho0/1.225), u10, gv, us)
1029 ustokes = us_to_u10*u10
1033 hm0 = 0.0246 *u10**2
1036 tmp = 2.0 *
pi * u19p5_to_u10 * u10
1037 fp = 0.877 * gv%mks_g_Earth / tmp
1040 fm = fm_into_fp * fp
1046 vstokes = 0.125 *
pi * r_loss * fm * hm0**2
1050 kphil = 0.176 * ustokes / vstokes
1056 kstar = kphil * 2.56
1058 z0 = abs(us%Z_to_m*hbl)
1061 r1 = ( 0.151 / kphil * z0i -0.84 ) * &
1062 ( 1.0 - exp(-2.0 * kphil * z0) )
1063 r2 = -( 0.84 + 0.0591 / kphil * z0i ) * &
1064 sqrt( 2.0 *
pi * kphil * z0 ) * &
1065 erfc( sqrt( 2.0 * kphil * z0 ) )
1066 r3 = ( 0.0632 / kstar * z0i + 0.125 ) * &
1067 (1.0 - exp(-2.0 * kstar * z0) )
1068 r4 = ( 0.125 + 0.0946 / kstar * z0i ) * &
1069 sqrt( 2.0 *
pi *kstar * z0) * &
1070 erfc( sqrt( 2.0 * kstar * z0 ) )
1071 ustokes_sl = ustokes * (0.715 + r1 + r2 + r3 + r4)
1072 la = sqrt(us%Z_to_m*us%s_to_T*ustar / ustokes_sl)
1084 real,
intent(in) :: AvgDepth
1085 real,
dimension(SZK_(GV)), &
1087 real,
dimension(SZK_(GV)), &
1088 intent(in) :: Profile
1090 real,
intent(out) :: Average
1093 real :: top, midpoint, bottom
1104 midpoint = bottom - gv%H_to_Z * 0.5*h(kk)
1105 bottom = bottom - gv%H_to_Z * h(kk)
1106 if (avgdepth < bottom)
then
1107 sum = sum + profile(kk) * (gv%H_to_Z * h(kk))
1108 elseif (avgdepth < top)
then
1109 sum = sum + profile(kk) * (top-avgdepth)
1114 average = sum / abs(avgdepth)
1123 real,
intent(in) :: AvgDepth
1124 integer,
intent(in) :: NB
1125 real,
dimension(NB), &
1126 intent(in) :: WaveNumbers
1127 real,
dimension(NB), &
1128 intent(in) :: SurfStokes
1129 real,
intent(out) :: Average
1139 average = average + surfstokes(bb) * &
1140 (1.-exp(-abs(avgdepth * 2.0 * wavenumbers(bb)))) / &
1141 abs(avgdepth * 2.0 * wavenumbers(bb))
1153 subroutine dhh85_mid(GV, US, zpt, UStokes)
1156 real,
intent(in) :: ZPT
1157 real,
intent(out) :: UStokes
1159 real :: ann, Bnn, Snn, Cnn, Dnn
1160 real :: omega_peak, omega, u10, WA, domega
1161 real :: omega_min, omega_max, wavespec, Stokes
1162 integer :: Nomega, OI
1172 domega = (omega_max-omega_min)/real(nomega)
1176 omega_peak = gv%mks_g_Earth / (wa * u10)
1178 omega_peak = 2. *
pi * 0.13 * gv%mks_g_Earth / u10
1181 ann = 0.006 *
waveage**(-0.55)
1183 snn = 0.08 * (1.0 + 4.0 *
waveage**3)
1186 cnn = cnn - 6.0*log10(wa)
1190 omega = omega_min + 0.5*domega
1192 dnn = exp( -0.5 * (omega-omega_peak)**2 / (snn**2 * omega_peak**2) )
1194 wavespec = (ann * gv%mks_g_Earth**2 / (omega_peak*omega**4 ) ) * &
1195 exp(-bnn*(omega_peak/omega)**4)*cnn**dnn
1197 stokes = 2.0 * wavespec * omega**3 * &
1198 exp( 2.0 * omega**2 * zpt / gv%mks_g_Earth) / gv%mks_g_Earth
1199 ustokes = ustokes + stokes*domega
1200 omega = omega + domega
1213 real,
intent(in) :: dt
1214 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),&
1216 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1218 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1223 real :: dtauup, dtaudn
1232 do i = g%iscB, g%iecB
1233 h_lay = gv%H_to_Z*0.5*(h(i,j,k)+h(i+1,j,k))
1236 dtauup = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i+1,j,k)) * &
1237 (waves%us_x(i,j,k-1)-waves%us_x(i,j,k)) / &
1238 (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k-1)+h(i+1,j,k-1)) ))
1241 dtaudn = 0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i+1,j,k+1)) * &
1242 (waves%us_x(i,j,k)-waves%us_x(i,j,k+1)) / &
1243 (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k+1)+h(i+1,j,k+1)) ))
1244 u(i,j,k) = u(i,j,k) + dt * (dtauup-dtaudn) / h_lay
1250 do j = g%jscB, g%jecB
1252 h_lay = gv%H_to_Z*0.5*(h(i,j,k)+h(i,j+1,k))
1255 dtauup = 0.5*(waves%Kvs(i,j,k)+waves%Kvs(i,j+1,k)) * &
1256 (waves%us_y(i,j,k-1)-waves%us_y(i,j,k)) / &
1257 (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k-1)+h(i,j+1,k-1)) ))
1260 dtaudn =0.5*(waves%Kvs(i,j,k+1)+waves%Kvs(i,j+1,k+1)) * &
1261 (waves%us_y(i,j,k)-waves%us_y(i,j,k+1)) / &
1262 (0.5*(h_lay + gv%H_to_Z*0.5*(h(i,j,k+1)+h(i,j+1,k+1)) ))
1263 v(i,j,k) = v(i,j,k) + dt * (dtauup-dtaudn) / h_lay
1281 real,
intent(in) :: dt
1282 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)), &
1284 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
1286 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)), &
1297 do i = g%iscB, g%iecB
1298 dvel = 0.25*(waves%us_y(i,j+1,k)+waves%us_y(i-1,j+1,k))*g%CoriolisBu(i,j+1) + &
1299 0.25*(waves%us_y(i,j,k)+waves%us_y(i-1,j,k))*g%CoriolisBu(i,j)
1300 u(i,j,k) = u(i,j,k) + dvel*us%s_to_T*dt
1306 do j = g%jscB, g%jecB
1308 dvel = 0.25*(waves%us_x(i+1,j,k)+waves%us_x(i+1,j-1,k))*g%CoriolisBu(i+1,j) + &
1309 0.25*(waves%us_x(i,j,k)+waves%us_x(i,j-1,k))*g%CoriolisBu(i,j)
1310 v(i,j,k) = v(i,j,k) - dvel*us%s_to_T*dt
1321 real,
intent(in) :: USTair
1322 real,
intent(out) :: U10
1327 real,
parameter :: vonkar = 0.4
1328 real,
parameter :: nu=1e-6
1329 real :: z0sm, z0, z0rough, u10a, alpha, CD
1338 z0sm = 0.11 * nu * us%m_to_Z / ustair
1339 u10 = ustair/sqrt(0.001)
1343 do while (abs(u10a/u10-1.) > 0.001)
1346 alpha = min(0.028, 0.0017 * u10 - 0.005)
1347 z0rough = alpha * (us%m_s_to_L_T*ustair)**2 / gv%g_Earth
1349 cd = ( vonkar / log(10.*us%m_to_Z / z0) )**2
1350 u10 = ustair/sqrt(cd)
1357 u10 = ustair/sqrt(0.0015)
1369 if (
allocated(cs%WaveNum_Cen))
deallocate( cs%WaveNum_Cen )
1370 if (
allocated(cs%Freq_Cen))
deallocate( cs%Freq_Cen )
1371 if (
allocated(cs%Us_x))
deallocate( cs%Us_x )
1372 if (
allocated(cs%Us_y))
deallocate( cs%Us_y )
1373 if (
allocated(cs%La_SL))
deallocate( cs%La_SL )
1374 if (
allocated(cs%La_turb))
deallocate( cs%La_turb )
1375 if (
allocated(cs%STKx0))
deallocate( cs%STKx0 )
1376 if (
allocated(cs%STKy0))
deallocate( cs%STKy0 )
1377 if (
allocated(cs%KvS))
deallocate( cs%KvS )
1378 if (
allocated(cs%Us0_y))
deallocate( cs%Us0_y )
1379 if (
allocated(cs%Us0_x))
deallocate( cs%Us0_x )