Return the asfc of a surface. The different grids are obtained by linear interpolation
| 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_lin_all(tab_in, scal, asfc_res) !================================================================================================ !! Return the *asfc* of a surface. The different grids are obtained by linear interpolation 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* 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) :: rr integer(kind=I4) :: i, ii, j, jj, ip, long_tmp, larg_tmp integer(kind=I4) :: nb_pt real(kind=R8), dimension(1:nb_beta) :: beta real(kind=R8), dimension(1:npp) :: vec_l real(kind=R8), dimension(1:npp) :: vec_x ! points coordinates real(kind=R8), dimension(1:npp) :: vec_y ! points coordinates real (kind=R8) :: asfc1, asfc2, aire_lin, xm, xp, ym, yp, hx, hy, hhx, hhy, h1, h2, h3, h4, width, height 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_lin_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 new_it = .true. nb_pt = 1 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) if (nb_pt>1) then vec_x(nb_pt -1) = log( (hhx*1e6)*(hhy*1e6)/2 ) vec_l(nb_pt -1) = log(aire_lin) endif if (out_lin .and. nb_pt>1) write(unit_out_lin, *) vec_x(nb_pt -1), vec_l(nb_pt -1) if (out_ter) write(*, *) aire_lin, long_tmp, larg_tmp endif if (ip == npp +1) then deallocate( x_new, y_new, x, y, tab_ou ) 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) 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) ) nb_pt = nb_pt +1 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 do j = 1, larg_new jj = locate(n=larg, xx=y(1:larg), x=y_new(j)) ym = y_new(j) -y(jj) ; yp = y_new(j) -y(jj +1) ym = ym/hy ; yp = yp/hy do i = 1, long_new ii = locate(n=long, xx=x(1:long), x=x_new(i)) xm = x_new(i) -x(ii) ; xp = x_new(i) -x(ii +1) xm = xm/hx ; xp = xp/hx h1 = tab_in(ii , jj ) h2 = tab_in(ii +1, jj ) h3 = tab_in(ii +1, jj +1) h4 = tab_in(ii , jj +1) tab_ou(i, j) = h1*xp*yp -h2*xm*yp +h3*xm*ym -h4*xp*ym enddo enddo enddo if (out_lin) close(unit_out_lin) !call system('python pyt/filetoplot.py &') !............................................................. nb_pt = nb_pt -1 vec_y(1:nb_pt) = vec_l(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_lin_all