Return the asfc of a surface. The different grids are obtained by Hermite 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 |
|
logical(kind=I4), | intent(in) | :: | omp |
with openmp ? |
subroutine calcul_asfc_hermite(tab_in, scal, asfc_res, omp) !================================================================================================ !! Return the *asfc* of a surface. The different grids are obtained by Hermite 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* logical(kind=I4), intent(in ) :: omp !! *with openmp ?* 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(:,:) :: tab_in_dx ! new grid function evaluations real(kind=R8), allocatable, dimension(:,:) :: tab_in_dy ! new grid function evaluations real(kind=R8), allocatable, dimension(:,:) :: tab_in_xy ! new grid function evaluations real(kind=R8), allocatable, dimension(:) :: gx real(kind=R8), allocatable, dimension(:) :: gy real(kind=R8), allocatable, dimension(:) :: gw real(kind=R8), allocatable, dimension(:) :: tab_dnq real(kind=R8) :: rr integer(kind=I4) :: i, ii, j, jj, k, long_tmp, larg_tmp real(kind=R8), dimension(1:nb_beta) :: beta real(kind=R8), dimension(1:npp) :: 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) :: asfc1, asfc2, aire, hx, hy, hhx, hhy, width, height real (kind=R8) :: xi, yi, eps_x integer(kind=I4) :: it, ng, nbpt integer(kind=I4) :: long, larg integer(kind=I4) :: nb_th long = scal%xres larg = scal%yres if (out_her) call get_unit(unit_out_her) if (out_her) open(unit=unit_out_her, file="out/asfc_her_her_all.txt") width = scal%lx height = scal%ly hx = width/(long -1) hy = height/(larg -1) eps_x = min(hx/10**3, hy/10**3) ! 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 allocate( x_new(1:long), & ! y_new(1:larg), & ! tab_ou(1:long, 1:larg), & ! tab_in_dx(1:long, 1:larg), & ! tab_in_dy(1:long, 1:larg), & ! tab_in_xy(1:long, 1:larg) ) ! 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 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) call init_aire_hermite(gx = gx, gy = gy, gw = gw, tab_dnq = tab_dnq, ng = ng) ! 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 call calcul_aire_hermite( tab_in = tab_ou(1:long_new, 1:larg_new), & ! long = long_new, & ! larg = larg_new, & ! gw = gw(1:ng), & ! tab_dnq = tab_dnq(1:ng), & ! ng = ng, & ! hx = hhx, & ! hy = hhy, & ! width = width, & ! height = height, & ! aire = aire) ! vec_s(1) = log( aire ) vec_x(1) = log( (hhx*1e6)*(hhy*1e6)/2 ) call calcul_tabd_hermite( tab_in = tab_in (1:long, 1:larg), & ! tab_dx = tab_in_dx(1:long, 1:larg), & ! tab_dy = tab_in_dy(1:long, 1:larg), & ! tab_xy = tab_in_xy(1:long, 1:larg), & ! long = long, & ! larg = larg, & ! hx = hx, & ! hy = hy ) ! nb_th = 1 if (omp) then nb_th = omp_get_num_procs() endif !$OMP PARALLEL DEFAULT(SHARED) NUM_THREADS(nb_th) IF (omp) !$OMP DO ORDERED SCHEDULE (STATIC,1) PRIVATE(it, long_tmp, larg_tmp, long_new, larg_new, hhx, hhy, i, x_new, j, y_new, yi, jj, xi, ii, tab_ou, aire) do it = 1, npp - 1 long_tmp = nint(long*(rr**it)) ! nb points en suite géométrique larg_tmp = nint(larg*(rr**it)) if ( long_new == long_tmp .or. larg_new == larg_tmp) then vec_s(it + 1) = 0 cycle ! à découper trop fin, on peut tomber sur les mêmes entiers endif long_new = long_tmp larg_new = larg_tmp hhx = width/(long_new -1) hhy = height/(larg_new -1) vec_x(it + 1) = log( (hhx*1e6)*(hhy*1e6)/2 ) 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 do j = 1, larg_new yi = y_new(j) jj = locate2(n = larg, xx = y(1:larg), x = yi, eps = eps_x) yi = (yi -(y(jj)+y(jj+1))/2)/(hy/2) do i = 1, long_new xi = x_new(i) ii = locate2(n = long, xx = x(1:long), x = xi, eps = eps_x) xi = (xi -(x(ii)+x(ii+1))/2)/(hx/2) tab_ou(i, j) = & ! nq_i(xi, yi, 1, 1)*tab_in(ii , jj ) + & ! u1 nq_i(xi, yi, 2, 1)*tab_in(ii +1, jj ) + & ! u2 nq_i(xi, yi, 3, 1)*tab_in(ii +1, jj +1) + & ! u3 nq_i(xi, yi, 4, 1)*tab_in(ii , jj +1) + & ! u4 nq_i(xi, yi, 1, 2)*tab_in_dx(ii , jj )*hx/2 + & ! du1/dx nq_i(xi, yi, 2, 2)*tab_in_dx(ii +1, jj )*hx/2 + & ! du2/dx nq_i(xi, yi, 3, 2)*tab_in_dx(ii +1, jj +1)*hx/2 + & ! du3/dx nq_i(xi, yi, 4, 2)*tab_in_dx(ii , jj +1)*hx/2 + & ! du4/dx nq_i(xi, yi, 1, 3)*tab_in_dy(ii , jj )*hy/2 + & ! du1/dy nq_i(xi, yi, 2, 3)*tab_in_dy(ii +1, jj )*hy/2 + & ! du2/dy nq_i(xi, yi, 3, 3)*tab_in_dy(ii +1, jj +1)*hy/2 + & ! du3/dy nq_i(xi, yi, 4, 3)*tab_in_dy(ii , jj +1)*hy/2 + & ! du4/dy nq_i(xi, yi, 1, 4)*tab_in_xy(ii , jj )*hx*hy/4 + & ! du1/dxdy nq_i(xi, yi, 2, 4)*tab_in_xy(ii +1, jj )*hx*hy/4 + & ! du2/dxdy nq_i(xi, yi, 3, 4)*tab_in_xy(ii +1, jj +1)*hx*hy/4 + & ! du3/dxdy nq_i(xi, yi, 4, 4)*tab_in_xy(ii , jj +1)*hx*hy/4 ! du4/dxdy enddo enddo call calcul_aire_hermite( tab_in = tab_ou(1:long_new, 1:larg_new), & ! long = long_new, & ! larg = larg_new, & ! gw = gw(1:ng), & ! tab_dnq = tab_dnq(1:ng), & ! ng = ng, & ! hx = hhx, & ! hy = hhy, & ! width = width, & ! height = height, & ! aire = aire) ! vec_s(it + 1) = log( aire ) !~ write(unit_out_her, *) vec_x(it+1), vec_s(it+1) enddo !$OMP END DO !$OMP END PARALLEL deallocate( x_new, y_new, x, y, tab_ou, tab_in_dx, tab_in_dy, tab_in_xy, tab_dnq, gx, gy, gw ) if (out_lin) close(unit_out_lin) if (out_her) close(unit_out_her) !............................................................. k = 1 do i = 1, npp if ( abs(vec_s(i)) > EPS_R8 ) then vec_y(k) = vec_s(i) vec_x(k) = vec_x(i) k = k + 1 endif enddo nbpt = k - 1 call interpolate_asfc( bt = beta(1:nb_beta), & ! n_bt = nb_beta, & ! n_pt = nbpt, & ! 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 subroutine init_aire_hermite(gx, gy, gw, tab_dnq, ng) !================================================================================================ !< @note ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none real (kind=R8), intent(out), allocatable, dimension(:) :: gx real (kind=R8), intent(out), allocatable, dimension(:) :: gy real (kind=R8), intent(out), allocatable, dimension(:) :: gw real (kind=R8), intent(out), allocatable, dimension(:) :: tab_dnq integer(kind=I4), intent(out) :: ng real (kind=R8) :: x1, x2, y1, y2 integer(kind=I4) :: i, k, nb_gauss_1d nb_gauss_1d = 2 ng = nb_gauss_1d**2 allocate( gx(1:ng), gy(1:ng), gw(1:ng) ) select case(nb_gauss_1d) case(1) x1 = 0._R8 y1 = 2._R8 gx(1:ng) = [ x1] gy(1:ng) = [ x1] gw(1:ng) = [y1*y1] case(2) x1 = sqrt(1._R8/3._R8) y1 = UN gx(1:ng) = [ -x1, -x1, +x1, +x1] gy(1:ng) = [ -x1, +x1, -x1, +x1] gw(1:ng) = [ y1*y1, y1*y1, y1*y1, y1*y1] case(3) x1 = sqrt(3._R8/5.0_R8) x2 = 0._R8 y1 = 5._R8/9._R8 y2 = 8._R8/9._R8 gx(1:ng) = [ -x1, -x1, -x1, x2, x2, x2, +x1, +x1, +x1] gy(1:ng) = [ -x1, x2, +x1, -x1, x2, +x1, -x1, x2, +x1] gw(1:ng) = [ y1*y1, y1*y2, y1*y1, y1*y2, y2*y2, y1*y2, y1*y1, y1*y2, y1*y1] endselect allocate( tab_dnq(1:32*ng) ) i = 1 do k = 1, ng tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 1) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 1) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 1) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 1) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 2) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 2) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 2) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 2) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 3) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 3) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 3) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 3) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 1, 4) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 2, 4) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 3, 4) ; i = i +1 tab_dnq(i) = dnq_xi_i(gx(k), gy(k), 4, 4) ; i = i +1 !! tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 1) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 1) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 1) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 1) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 2) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 2) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 2) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 2) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 3) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 3) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 3) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 3) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 1, 4) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 2, 4) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 3, 4) ; i = i +1 tab_dnq(i) = dnq_et_i(gx(k), gy(k), 4, 4) ; i = i +1 enddo return endsubroutine init_aire_hermite subroutine calcul_tabd_hermite( tab_in, tab_dx, tab_dy, tab_xy, long, larg, hx, hy) !================================================================================================ !< @note ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in ) :: long integer(kind=I4), intent(in ) :: larg real (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in real (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_dx real (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_dy real (kind=R8), intent(out), dimension(1:long, 1:larg) :: tab_xy real (kind=R8), intent(in ) :: hx real (kind=R8), intent(in ) :: hy integer(kind=I4) :: i, im, ip, j, jm, jp real (kind=R8) :: ui, uim, uip, ujm, ujp, upp, ump, upm, umm do j = 1, larg jm = max(j -1, 1) jp = min(j +1, larg) do i = 1, long im = max(i -1, 1) ip = min(i +1, long) ui = tab_in(i , j ) uim = tab_in(im, j ) uip = tab_in(ip, j ) ujm = tab_in(i , jm) ujp = tab_in(i , jp) upp = tab_in(ip, jp) ump = tab_in(im, jp) upm = tab_in(ip, jm) umm = tab_in(im, jm) tab_dx(i, j) = (uip -uim)/(2*hx) tab_dy(i, j) = (ujp -ujm)/(2*hy) tab_xy(i, j) = (upp -ump -upm +umm)/(4*hx*hy) enddo enddo tab_dx( 1, 1:larg) = ( tab_in( 2, 1:larg) - & ! tab_in( 1, 1:larg) )/hx ! tab_dx(long, 1:larg) = ( tab_in(long , 1:larg) - & ! tab_in(long -1, 1:larg) )/hx ! tab_dy(1:long, 1) = ( tab_in(1:long, 2) - & ! tab_in(1:long, 1) )/hy ! tab_dy(1:long, larg) = ( tab_in(1:long, larg ) - & ! tab_in(1:long, larg -1) )/hy ! tab_xy( 1, 1:larg) = ( tab_dy( 2, 1:larg) - & ! tab_dy( 1, 1:larg) )/hx ! tab_xy(long, 1:larg) = ( tab_dy(long , 1:larg) - & ! tab_dy(long -1, 1:larg) )/hx ! tab_xy(1:long, 1) = ( tab_dx(1:long, 2) - & ! tab_dx(1:long, 1) )/hy ! tab_xy(1:long, larg) = ( tab_dx(1:long, larg ) - & ! tab_dx(1:long, larg -1) )/hy ! return endsubroutine calcul_tabd_hermite subroutine calcul_aire_hermite(tab_in, long, larg, gw, tab_dnq, ng, hx, hy, width, height, aire) !================================================================================================ !< @note ! ! @endnote !------------------------------------------------------------------------------------------------ implicit none integer(kind=I4), intent(in ) :: long integer(kind=I4), intent(in ) :: larg integer(kind=I4), intent(in ) :: ng real (kind=R8), intent(in ), dimension(1:long, 1:larg) :: tab_in real (kind=R8), intent(in ), dimension(1:ng) :: gw real (kind=R8), intent(in ), dimension(1:32*ng) :: tab_dnq real (kind=R8), intent(in ) :: hx real (kind=R8), intent(in ) :: hy real (kind=R8), intent(in ) :: width real (kind=R8), intent(in ) :: height real (kind=R8), intent(out) :: aire integer(kind=I4) :: i, j, k, i1, i2, j1, j2 real (kind=R8) :: aire_tmp real(kind=R8), allocatable, dimension(:) :: dfx real(kind=R8), allocatable, dimension(:) :: dfy real(kind=R8), allocatable, dimension(:,:) :: tab_dx ! new grid function evaluations real(kind=R8), allocatable, dimension(:,:) :: tab_dy ! new grid function evaluations real(kind=R8), allocatable, dimension(:,:) :: tab_xy ! new grid function evaluations allocate( dfx(1:ng), & ! dfy(1:ng) ) ! allocate( tab_dx (1:long, 1:larg), & ! tab_dy (1:long, 1:larg), & ! tab_xy (1:long, 1:larg) ) ! call calcul_tabd_hermite( tab_in = tab_in(1:long, 1:larg), & ! tab_dx = tab_dx(1:long, 1:larg), & ! tab_dy = tab_dy(1:long, 1:larg), & ! tab_xy = tab_xy(1:long, 1:larg), & ! long = long, & ! larg = larg, & ! hx = hx, & ! hy = hy ) ! aire_tmp = 0._R8 do j = 1, larg -1 j1 = j ; j2 = j +1 do i = 1, long -1 i1 = i ; i2 = i +1 do k = 1, ng dfx(k) = (2._R8/hx)*( & ! tab_dnq(32*(k -1) +01)*tab_in(i1, j1) + & ! u1 tab_dnq(32*(k -1) +02)*tab_in(i2, j1) + & ! u2 tab_dnq(32*(k -1) +03)*tab_in(i2, j2) + & ! u3 tab_dnq(32*(k -1) +04)*tab_in(i1, j2) + & ! u4 tab_dnq(32*(k -1) +05)*tab_dx(i1, j1)*hx/2 + & ! du1/dx tab_dnq(32*(k -1) +06)*tab_dx(i2, j1)*hx/2 + & ! du2/dx tab_dnq(32*(k -1) +07)*tab_dx(i2, j2)*hx/2 + & ! du3/dx tab_dnq(32*(k -1) +08)*tab_dx(i1, j2)*hx/2 + & ! du4/dx tab_dnq(32*(k -1) +09)*tab_dy(i1, j1)*hy/2 + & ! du1/dy tab_dnq(32*(k -1) +10)*tab_dy(i2, j1)*hy/2 + & ! du2/dy tab_dnq(32*(k -1) +11)*tab_dy(i2, j2)*hy/2 + & ! du3/dy tab_dnq(32*(k -1) +12)*tab_dy(i1, j2)*hy/2 + & ! du4/dy tab_dnq(32*(k -1) +13)*tab_xy(i1, j1)*hx*hy/4 + & ! du1/dxdy tab_dnq(32*(k -1) +14)*tab_xy(i2, j1)*hx*hy/4 + & ! du2/dxdy tab_dnq(32*(k -1) +15)*tab_xy(i2, j2)*hx*hy/4 + & ! du3/dxdy tab_dnq(32*(k -1) +16)*tab_xy(i1, j2)*hx*hy/4 ) ! du4/dxdy dfy(k) = (2._R8/hy)*( & ! tab_dnq(32*(k -1) +17)*tab_in(i1, j1) + & ! u1 tab_dnq(32*(k -1) +18)*tab_in(i2, j1) + & ! u2 tab_dnq(32*(k -1) +19)*tab_in(i2, j2) + & ! u3 tab_dnq(32*(k -1) +20)*tab_in(i1, j2) + & ! u4 tab_dnq(32*(k -1) +21)*tab_dx(i1, j1)*hx/2 + & ! du1/dx tab_dnq(32*(k -1) +22)*tab_dx(i2, j1)*hx/2 + & ! du2/dx tab_dnq(32*(k -1) +23)*tab_dx(i2, j2)*hx/2 + & ! du3/dx tab_dnq(32*(k -1) +24)*tab_dx(i1, j2)*hx/2 + & ! du4/dx tab_dnq(32*(k -1) +25)*tab_dy(i1, j1)*hy/2 + & ! du1/dy tab_dnq(32*(k -1) +26)*tab_dy(i2, j1)*hy/2 + & ! du2/dy tab_dnq(32*(k -1) +27)*tab_dy(i2, j2)*hy/2 + & ! du3/dy tab_dnq(32*(k -1) +28)*tab_dy(i1, j2)*hy/2 + & ! du4/dy tab_dnq(32*(k -1) +29)*tab_xy(i1, j1)*hx*hy/4 + & ! du1/dxdy tab_dnq(32*(k -1) +30)*tab_xy(i2, j1)*hx*hy/4 + & ! du2/dxdy tab_dnq(32*(k -1) +31)*tab_xy(i2, j2)*hx*hy/4 + & ! du3/dxdy tab_dnq(32*(k -1) +32)*tab_xy(i1, j2)*hx*hy/4 ) ! du4/dxdy enddo do k = 1, ng aire_tmp = aire_tmp + gw(k) * sqrt( UN + dfx(k)**2 + dfy(k)**2 ) enddo enddo enddo aire = aire_tmp*( hx/2 )*( hy/2 ) aire = aire/(width*height) deallocate( tab_dx , & ! tab_dy , & ! tab_xy ) ! deallocate( dfx, dfy ) return endsubroutine calcul_aire_hermite endsubroutine calcul_asfc_hermite