! $Id: fsi_lsw.f90 2023-01-20 09:51:28Z ychen $ !****s* bangle/fsi_lsw * ! ! NAME ! fsi_lsw - calculate the bending angle local spectral width ! ! SYNOPSIS ! call fsi_lsw(roc, pp, ee, ppLSW, eeLSW, LSWMG, LSWSS, nout) ! ! DESCRIPTION ! This subroutine calculate the bending angle local spectral width ! ! INPUTS ! REAL(wp), :: roc Radius of curvature (m) ! REAL(wp), DIMENSION(:) :: pp impact parameter (m) ! REAL(wp), DIMENSION(:) :: ee bending angle (rad) ! ! OUTPUT ! REAL(wp), DIMENSION(:) :: ppLSW output impact parameter ! REAL(wp), DIMENSION(:) :: eeLSW output bending angle (rad) ! REAL(wp), DIMENSION(:) :: LSWMG LSW calculated using squared dalp (rad^2) ! REAL(wp), DIMENSION(:) :: LSWSS LSW calculated using absolute dalp (rad^2) ! INTEGER, :: nout array size for impact and bending !! ! AUTHOR ! Yong Chen, NOAA/NESDIS/STAR, yong.chen@noaa.gov ! Loknath Adhikari, UMD/CISESS, loknath.adhikari@noaa.gov ! ! COPYRIGHT ! Copyright (c) 2022-2023 Yong Chen ! For further details please refer to the file COPYRIGHT ! which you should have received as part of this distribution. ! !**** SUBROUTINE fsi_lsw(roc, pp, ee, ppLSW, eeLSW, LSWMG, LSWSS, nout) USE typesizes, ONLY: wp => EightByteReal USE ropp_pp, not_this => fsi_lsw USE ropp_pp_spline USE messages USE ropp_pp_spline IMPLICIT none REAL(wp), INTENT(in) :: roc ! Radius of curvature (m) REAL(wp), DIMENSION(:), INTENT(in) :: pp ! impact parameter (m) REAL(wp), DIMENSION(:), INTENT(in) :: ee ! bending angle (rad) REAL(wp), DIMENSION(:), INTENT(inout) :: ppLSW ! output impact parameter REAL(wp), DIMENSION(:), INTENT(inout) :: eeLSW ! output bending angle (rad) REAL(wp), DIMENSION(:), INTENT(inout) :: LSWMG ! LSW calculated using squared dalp (rad^2) REAL(wp), DIMENSION(:), INTENT(inout) :: LSWSS ! LSW calculated using absolute dalp (rad^2) INTEGER, INTENT(out) :: nout ! array size for impact and bending INTEGER, PARAMETER :: ndalp = 121 ! bin size for bending angle REAL(wp), PARAMETER :: dp = 100.0_wp ! low-pass filetering window size (m) REAL(wp), PARAMETER :: dsw = 500.0_wp ! sliding window size for spectral analysis (m) INTEGER, PARAMETER :: nmax1 = 3000 ! number of data point in the extended signal REAL(wp), PARAMETER :: emin = -0.03_wp ! minimum bending angle departure from mean REAL(wp), PARAMETER :: emax = 0.03_wp ! maximum bending angle departure from mean REAL(wp), PARAMETER :: dplsw = 50.0_wp ! LWS vertical resolution (m) REAL(wp), PARAMETER :: PH = 15550.0_wp ! The highest impact height for LSW (m) INTEGER, PARAMETER :: np_smooth = 3 ! Polynomial degree for smoothing regression CHARACTER(len = 256) :: routine REAL(wp), DIMENSION(:), ALLOCATABLE :: smooth_ba REAL(wp), DIMENSION(:), ALLOCATABLE :: dee REAL(wp), DIMENSION(:), ALLOCATABLE :: cdf REAL(wp), DIMENSION(:), ALLOCATABLE :: dalp REAL(wp), DIMENSION(:), ALLOCATABLE :: de1 REAL(wp), DIMENSION(:), ALLOCATABLE :: de2 REAL(wp), DIMENSION(:, :), ALLOCATABLE :: pdf REAL(wp), DIMENSION(:), ALLOCATABLE :: pp0 REAL(wp), DIMENSION(:), ALLOCATABLE :: ee0 REAL(wp), DIMENSION(:), ALLOCATABLE :: ee0sm REAL(wp), DIMENSION(:), ALLOCATABLE :: dee1 REAL(wp) :: Pmin, Pmax, dp_fine, deltaE REAL(wp) :: dx, p1, p2, p12, dt REAL(wp) :: ang1, ang2, pwr1, ang, uu, vv INTEGER :: nmax, nn, ik, i1, i2, im, nd, nlsw, nbsize INTEGER :: i, j, k, ii, jj, nws ! 1.4 Local arrays !CALL message_get_routine(routine) !CALL message_set_routine('fsi_lsw') nmax = SIZE(pp) Pmax = MAXVAL(pp(1:nmax)) Pmin = MINVAL(pp(1:nmax)) ! only apply for below impact height of 15 km IF (Pmin < roc + PH ) THEN ALLOCATE(smooth_ba(nmax)) dx = (Pmax - Pmin) / (nmax - 1) nws = NINT(dp / dx) ! smoothing the bending angle with 100 meter window size CALL ropp_pp_sliding_polynomial(pp(1:nmax), ee(1:nmax), & nws, np_smooth, smooth_ba(1:nmax)) ! Interpolate to 1 meter uniform vertical grid dp_fine = 1.0_wp nn = NINT((Pmax - Pmin)/dp_fine)+1 ALLOCATE(pp0(nn)) ALLOCATE(ee0(nn)) ALLOCATE(ee0sm(nn)) ALLOCATE(dee1(nn)) ALLOCATE(dee(nmax1)) ALLOCATE(cdf(nmax1)) ALLOCATE(pdf(nmax1, ndalp)) ALLOCATE(dalp(ndalp)) ALLOCATE(de1(ndalp)) ALLOCATE(de2(ndalp)) pp0(1) = Pmin DO i=2,nn pp0(i) = pp0(i-1)+dp_fine ENDDO CALL ropp_pp_interpol(pp(1:nmax), pp0(1:nn), & smooth_ba(1:nmax), ee0sm(1:nn)) CALL ropp_pp_interpol(pp(1:nmax), pp0(1:nn), & ee(1:nmax), ee0(1:nn)) ! bending angle departure from the mean dee1 = ee0sm-ee0 ! bin size of 0.0005 radian deltaE = (emax - emin)/(ndalp - 1) DO ii = 1, ndalp de1(ii) = emin + (ii - 1)*deltaE - 0.5_wp * deltaE de2(ii) = de1(ii) + deltaE dalp(ii) = (de1(ii)+de2(ii))/2.0_wp ENDDO ! redefine Pmin and Pmax for uncertainty within [0 to 15 km] IF (Pmin < roc) Pmin = roc Pmax = roc+PH nlsw = NINT((Pmax-Pmin)/dplsw) ! half sliding window size nd = NINT(dsw/(2.0_wp * dplsw)) DO i = nd, nlsw - nd + 1, 1 k = i - nd + 1 p1 = Pmin + (dplsw * i) - dsw/2.0_wp p2 = p1 + dsw p12 = (p1 + p2)/2.0_wp i1 = MIN(SUM(MINLOC(pp0(:), p2 > pp0(:) .AND. pp0(:) >= p1)), & SUM(MAXLOC(pp0(:), p2 > pp0(:) .AND. pp0(:) >= p1))) i2 = MAX(SUM(MINLOC(pp0(:), p2 > pp0(:) .AND. pp0(:) >= p1)), & SUM(MAXLOC(pp0(:), p2 > pp0(:) .AND. pp0(:) >= p1))) im = SUM(MAXLOC(pp0(:), pp0(:) <= p12 ) ) ik = i2 - i1 + 1 ! bending angle deviation from the mean value dee(1:ik) = dee1(i1:i2) ! sort the data in increasing order DO ii = 1, ik - 1 DO jj = ii + 1, ik IF (dee(ii) >= dee(jj)) THEN dt = dee(ii) dee(ii) = dee(jj) dee(jj) = dt ENDIF ENDDO ENDDO ! Obtain the number of each 0.0005 radian bins cdf(k) = 0 DO ii = 1, ndalp nbsize = COUNT(Mask = (dee(1:ik) >= de1(ii) .AND. dee(1:ik) < de2(ii))) pdf(k, ii) = nbsize cdf(k) = cdf(k) + pdf(k, ii) !write(*, *), 'k, ii, pdf, cdf', k, ii, pdf(k, ii), cdf(k) ENDDO !write(*, *) 'i1, i2, im, ik:', i1, i2, im, ik, p1, p2, p12 ppLSW(k) = pp0(im) eeLSW(k) = ee0(im) ENDDO DO i = 1, k ang1 = 0.0_wp ang2 = 0.0_wp DO j = 1, ndalp pwr1 = pdf(i, j)/cdf(i) ang = ABS(dalp(j)) uu = (pwr1*ang**2) vv = (pwr1*ang) ang1 = ang1+uu ang2 = ang2+vv ENDDO LSWMG(i) = ang1 !SQRT(ang1) LSWSS(i) = ang2**2 !ang2 ENDDO nout = k Pmin = min(ppLSW(1),ppLSW(nout)) Pmax = max(ppLSW(1),ppLSW(nout)) write(*,*)'fsi_lsw: numbers at end: ',nmax, nout, Pmin, Pmax !write(*,*) ppLSW(1:k) !write(*,*) eeLSW(1:k) !write(*,*) sqrt(LSWMG(1:k)) DEALLOCATE(smooth_ba) DEALLOCATE(pp0) DEALLOCATE(ee0) DEALLOCATE(ee0sm) DEALLOCATE(dee1) DEALLOCATE(dee) DEALLOCATE(cdf) DEALLOCATE(pdf) DEALLOCATE(dalp) DEALLOCATE(de1) DEALLOCATE(de2) ELSE nout = 0 ENDIF END SUBROUTINE fsi_lsw