Return the asfc of a surface. The different grids are obtained by spline of degree 3
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| real(kind=R8), | intent(in), | dimension(1:scal%xres, 1:scal%yres) | :: | tab_in |
input surface |
|
| type(SCALE_SURF), | intent(in) | :: | scal |
surface characteristics |
||
| real(kind=R8), | intent(out), | dimension(1:2) | :: | asfc_res |
result: asfc, adjustment factor |
subroutine calcul_asfc_spl_all(tab_in, scal, asfc_res) !================================================================================================ !! Return the *asfc* of a surface. The different grids are obtained by spline of degree 3 implicit none type(SCALE_SURF), intent(in ) :: scal !! *surface characteristics* real (kind=R8), intent(in ), dimension(1:scal%xres, 1:scal%yres) :: tab_in !! *input surface* real (kind=R8), intent(out), dimension(1:2) :: asfc_res !! *result: asfc, adjustment factor* integer(kind=I4), parameter :: idx = 0 ! [[db2val]] input integer(kind=I4), parameter :: idy = 0 ! [[db2val]] input real(kind=R8), allocatable, dimension(:) :: x ! x points in original grid real(kind=R8), allocatable, dimension(:) :: y ! y points in original grid integer(kind=I4) :: long_new ! number of points in x dimension for new grid integer(kind=I4) :: larg_new ! number of points in y dimension for new grid real(kind=R8), allocatable, dimension(:) :: x_new ! new grid x points real(kind=R8), allocatable, dimension(:) :: y_new ! new grid y points real(kind=R8), allocatable, dimension(:,:) :: tab_ou ! new grid function evaluations real(kind=R8), allocatable, dimension(:,:) :: coeff, coeff_tmp real(kind=R8), allocatable, dimension(:) :: tx, tx_tmp ! x knots real(kind=R8), allocatable, dimension(:) :: ty, ty_tmp ! y knots real(kind=R8) :: val, rr integer(kind=I4) :: i, j, k, ip, long_tmp, larg_tmp integer(kind=I4) :: iflag integer(kind=I4) :: inbvx, inbvy, iloy integer(kind=I4) :: nb_pt integer(kind=I4), dimension(1:8) :: d1, d2 real(kind=R8), dimension(1:nb_beta) :: beta real(kind=R8), dimension(1:npp) :: vec_l, vec_s real(kind=R8), dimension(1:npp) :: vec_x ! points coordinates real(kind=R8), dimension(1:npp) :: vec_y ! points coordinates real(kind=R8), dimension(1:8) :: gx, gy, vf real(kind=R8) :: asfc1, asfc2, aire, aire_lin, aire_tmp, x1, x2, y1, y2, val_x, val_y, width, height, hx, hy, hhx, hhy logical(kind=I4) :: new_it integer(kind=I4) :: long, larg long = scal%xres larg = scal%yres if (out_lin) open(unit=unit_out_lin, file="out/asfc_lin_spl_all.txt") if (out_spl) open(unit=unit_out_spl, file="out/asfc_spl_spl_all.txt") width = scal%lx height = scal%ly hx = width/(long -1) hy = height/(larg -1) ! définition d'abscisses pour l'interpolation par splines allocate( x(1:long), y(1:larg) ) do i = 1, long x(i) = hx*real(i-1, kind=R8) enddo do j = 1, larg y(j) = hy*real(j-1, kind=R8) enddo rr = (real(long0, kind=R8)/long)**(UN/npp) ! facteur de réduction pour aller du maillage initial au maillage minimal avec npp points rr = max( rr, (real(larg0, kind=R8)/larg)**(UN/npp) ) ! facteur de réduction pour aller du maillage initial au maillage minimal avec npp points allocate( x_new(1:long), & y_new(1:larg), & ! nouvelles abscisses tab_ou(1:long, 1:larg) ) x_new(1:long) = x(1:long) y_new(1:larg) = y(1:larg) tab_ou(1:long, 1:larg) = tab_in(1:long, 1:larg) ! pour chaque réduction de maillage, calcul du maillage résultant et de l'aire relative associée !............................................................. long_new = long larg_new = larg hhx = hx hhy = hy nb_pt = 1 new_it = .true. do ip = 1, npp +1 if (new_it) then long_tmp = long_new larg_tmp = larg_new ! calcul de l'aire relative call calcul_aire(tab_ou(1:long_tmp, 1:larg_tmp), long_tmp, larg_tmp, hhx, hhy, aire_lin)!, rx=real(long-1, kind=R8)/(long_tmp-1), & ! ry=real(larg-1, kind=R8)/(larg_tmp-1)) allocate( coeff_tmp(1:long_tmp, 1:larg_tmp) ) allocate( tx_tmp(1:(long_tmp +kx)), & ty_tmp(1:(larg_tmp +ky)) ) iflag = 0 call db2ink( x = x_new(1:long_tmp), & nx = long_tmp, & y = y_new(1:larg_tmp), & ny = larg_tmp, & fcn = tab_ou(1:long_tmp, 1:larg_tmp), & kx = kx, & ky = ky, & tx = tx_tmp(1:(long_tmp +kx)), & ty = ty_tmp(1:(larg_tmp +ky)), & bcoef = coeff_tmp(1:long_tmp, 1:larg_tmp), & iflag = iflag) if (iflag/=1) then write(*,*) iflag, long_tmp, larg_tmp, kx, ky stop 'error calling db2ink' endif if (ip==1) then allocate( coeff(1:long, 1:larg) ) allocate( tx(1:(long +kx)), & ty(1:(larg +ky)) ) coeff(1:long, 1:larg) = coeff_tmp(1:long, 1:larg) tx(1:(long +kx)) = tx_tmp(1:(long +kx)) ty(1:(larg +ky)) = ty_tmp(1:(larg +ky)) endif inbvx = 1 inbvy = 1 iloy = 1 aire_tmp = 0._R8 !~ !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(NB_PROCS) IF(MULTI_PROCS2) !~ !$OMP DO SCHEDULE (STATIC,(larg_tmp-1)/NB_PROCS) PRIVATE(i, k, iflag, val, y1, y2, x1, x2, d1, d2, gx, gy, vf) FIRSTPRIVATE(inbvx, inbvy, iloy) REDUCTION(+:aire_tmp) do j = 1, larg_tmp -1 y1 = y_new(j) +(hhy/2)*(UN -UN/sqrt(3._R8)) y2 = y_new(j) +(hhy/2)*(UN +UN/sqrt(3._R8)) do i = 1, long_tmp -1 x1 = x_new(i) +(hhx/2)*(UN -UN/sqrt(3._R8)) x2 = x_new(i) +(hhx/2)*(UN +UN/sqrt(3._R8)) d1(1:8) = [ 1, 0, 1, 0, 1, 0, 1, 0] d2(1:8) = [ 0, 1, 0, 1, 0, 1, 0, 1] gx(1:8) = [x1, x1, x1, x1, x2, x2, x2, x2] gy(1:8) = [y1, y1, y2, y2, y1, y1, y2, y2] do k = 1, 8 call db2val(xval = gx(k), & yval = gy(k), & idx = d1(k), & idy = d2(k), & tx = tx_tmp(1:(long_tmp +kx)), & ty = ty_tmp(1:(larg_tmp +ky)), & nx = long_tmp, & ny = larg_tmp, & kx = kx, & ky = ky, & bcoef = coeff_tmp(1:long_tmp, 1:larg_tmp), & f = vf(k), & iflag = iflag, & inbvx = inbvx, & inbvy = inbvy, & iloy = iloy) if (iflag/=0) then write(*,*) iflag stop 'error calling db2val' endif enddo do k = 1, 4 val_x = vf(2*k -1) val_y = vf(2*k ) aire_tmp = aire_tmp +sqrt( UN +val_x**2 +val_y**2 ) enddo enddo enddo !~ !$OMP END DO !~ !$OMP END PARALLEL aire = aire_tmp*( hhx/2 )*( hhy/2 ) aire = aire/(width*height) vec_x(nb_pt) = log( (hhx*1e6)*(hhy*1e6)/2 ) vec_l(nb_pt) = log(aire_lin) vec_s(nb_pt) = log(aire ) if (out_lin .and. nb_pt>1) write(unit_out_lin, *) vec_x(nb_pt -1), vec_l(nb_pt -1) if (out_spl .and. nb_pt>1) write(unit_out_spl, *) vec_x(nb_pt -1), vec_s(nb_pt -1) if (out_ter) write(*, *) hhx, hhy, aire, ip endif if (ip == npp +1) then deallocate( x_new, y_new, x, y, tab_ou, coeff, tx, ty, coeff_tmp, tx_tmp, ty_tmp ) exit endif long_new = nint(long*(rr**ip)) ! nb points en suite géométrique larg_new = nint(larg*(rr**ip)) new_it = .true. if (long_new==long_tmp .or. larg_new==larg_tmp) then new_it = .false. cycle ! à découper trop fin, on peut tomber sur les mêmes entiers endif hhx = width/(long_new -1) hhy = height/(larg_new -1) nb_pt = nb_pt +1 deallocate( x_new, y_new, tab_ou ) allocate( x_new(1:long_new), & y_new(1:larg_new), & ! nouvelles abscisses tab_ou(1:long_new, 1:larg_new) ) do i = 1, long_new x_new(i) = hhx*real(i-1, kind=R8) enddo do j = 1, larg_new y_new(j) = hhy*real(j-1, kind=R8) enddo inbvx = 1 inbvy = 1 iloy = 1 ! calcul des hauteurs de la surface interpolée !~ !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(NB_PROCS) IF(MULTI_PROCS2) !~ !$OMP DO SCHEDULE (STATIC,larg_new/NB_PROCS) PRIVATE(i, iflag, val) FIRSTPRIVATE(inbvx, inbvy, iloy) do j = 1, larg_new do i = 1, long_new call db2val(xval = x_new(i), & yval = y_new(j), & idx = idx, & idy = idy, & tx = tx(1:(long +kx)), & ty = ty(1:(larg +ky)), & nx = long, & ny = larg, & kx = kx, & ky = ky, & bcoef = coeff(1:long, 1:larg), & f = val, & iflag = iflag, & inbvx = inbvx, & inbvy = inbvy, & iloy = iloy) if (iflag/=0) error stop 'error calling db2val' iflag = 0 tab_ou(i, j) = val enddo enddo !~ !$OMP END DO !~ !$OMP END PARALLEL deallocate( coeff_tmp, tx_tmp, ty_tmp ) enddo if (out_lin) close(unit_out_lin) if (out_spl) close(unit_out_spl) !............................................................. nb_pt = nb_pt -1 vec_y(1:nb_pt) = vec_s(1:nb_pt) call interpolate_asfc( bt = beta(1:nb_beta), & n_bt = nb_beta, & n_pt = nb_pt, & asf = asfc1, & r2 = asfc2) asfc_res = [-1000*asfc1, asfc2] return contains subroutine interpolate_asfc(bt, n_bt, n_pt, asf, r2) !================================================================================================ !< @note Function that fits the data points for the asfc determination.<br/> ! \(f\_boltz(x_i)=\beta_2 + \beta_3 \tanh \left( \dfrac{x_i -\beta_1}{\beta_4} \right) \) ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in ) :: n_bt !! *number of parameters* integer(kind=I4), intent(in ) :: n_pt !! *data vector length* real (kind=R8), intent(out), dimension(1:n_bt) :: bt !! *vector \(\beta\) of parameters* real (kind=R8), intent(out) :: asf !! *Asfc number* real (kind=R8), intent(out) :: r2 !! *correlation number to assess validity of the Asfc calculus* real (kind=R8), dimension(1:n_pt, 1:n_bt) :: jacob real (kind=R8), dimension(1:n_pt) :: v_x, v_y, res_y, pentes integer(kind=I4) :: i0, i, j, info real (kind=R8) :: delta1, delta2, y_mean v_x(1:n_pt) = vec_x(1:n_pt) ! smoothing v_y(1) = vec_y(1) do i = 1 + 1, n_pt - 1 v_y(i) = 0.25_R8 * ( vec_y(i - 1) + 2 * vec_y(i) + vec_y(i + 1) ) enddo v_y(n_pt) = vec_y(n_pt) call init_beta_boltz( bt = bt(1:n_bt), & ! n_bt = n_bt, & ! v_x = v_x(1:n_pt), & ! v_y = v_y(1:n_pt), & ! n_pt = n_pt) ! res_y(1:n_pt) = 0._R8 jacob(1:n_pt, 1:n_bt) = 0._R8 call lmder1( fcn = lmder1_f, & ! m = n_pt, & ! n = n_bt, & ! x = bt(1:n_bt), & ! fvec = res_y(1:n_pt), & ! fjac = jacob(1:n_pt, 1:n_bt), & ! ldfjac = n_pt, & ! tol = 1.0e-8_R8, & ! info = info) ! delta1 = 0._R8 delta2 = 0._R8 y_mean = sum( vec_y(1:n_pt) )/n_pt do i = 1, n_pt delta1 = delta1 +( vec_y(i) -f_boltz(xi = vec_x(i), & ! beta = bt(1:n_bt), & ! n_beta = n_bt) & ! )**2 delta2 = delta2 +( vec_y(i) -y_mean )**2 enddo r2 = 1._R8 -delta1/delta2 !~ if (r2<0) then !~ do i=1,n_pt !~ write(99,*) v_x(i), v_y(i) !~ enddo !~ stop 'error' !~ endif i0 = locate(n = n_pt, & ! xx = v_x(1:n_pt), & ! x = bt(1)) ! j = 0 do i = i0 -5, i0 +5 if (i<1 .or. i>n_pt) cycle j = j +1 pentes(j) = +df_boltz( xi = v_x(i), & ! beta = bt(1:n_bt), & ! n_beta = n_bt, & ! ivar = 0) ! enddo asf = sum( pentes(1:j)/j ) return endsubroutine interpolate_asfc subroutine lmder1_f(m, n, x, fvec, fjac, ldfjac, iflag) !================================================================================================ !< @note Function called by [[lmder1]] as part of **minpack**, modified by ! [Burkardt](https://people.sc.fsu.edu/~jburkardt/f_src/minpack/minpack.html) <br/> ! According *iflag* value it calculates the function [[f_boltz]] at the data points ! or the jacobian. ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in ) :: m !! *number of points* integer(kind=I4), intent(in ) :: n !! *number of parameters* integer(kind=I4), intent(in ) :: ldfjac !! *leading dimension of fjac, which must be at least n* integer(kind=I4), intent(in ) :: iflag !! *which calculus to perform* real (kind=R8), intent(out ), dimension(1:m) :: fvec !! *vector of f_boltz(xi)* real (kind=R8), intent(out ), dimension(1:ldfjac, 1:n) :: fjac !! *jacobian* real (kind=R8), intent(inout), dimension(1:n) :: x !! *parameter values* integer(kind=I4) :: i, k select case (iflag) case(0) continue case(1) do i = 1, m fvec(i) = f_boltz( xi = vec_x(i), & ! beta = x(1:n), & ! n_beta = n) & ! -vec_y(i) enddo case(2) do i = 1, m do k = 1, n fjac(i, k) = df_boltz(xi = vec_x(i), & ! beta = x(1:n), & ! ivar = k, & ! n_beta = n) ! enddo enddo case default write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'LMDER1_F - Fatal error!' write ( *, '(a,i6)' ) ' Called with unexpected value of IFLAG = ', iflag stop endselect return endsubroutine lmder1_f endsubroutine calcul_asfc_spl_all