Index: adw_tricub_lag3d_vec.ftn =================================================================== --- adw_tricub_lag3d_vec.ftn 2006-06-28 14:36:54.000000000 -0400 +++ adw_tricub_lag3d_vec.ftn 2006-06-28 14:31:50.000000000 -0400 @@ -0,0 +1,325 @@ +***s/r adw_tricub_lag3d_vec - Tri-cubic interpolation: Lagrange 3d +* +#include +* + subroutine adw_tricub_lag3d_vec ( F_out, F_in, F_x, F_y, F_z, + % F_num, F_mono_L, i0, in, j0, jn, kn ) +* + implicit none +* + logical F_mono_L +* + integer F_num, i0, in, j0, jn, kn +* + real F_in(*) +* + real F_out (F_num), F_x(F_num), F_y(F_num), F_z(F_num) +* +*authors +* Gravel & Valin & Tanguay +* +* (Based on adw_tricub v_3.1.1) +* +*revision +* v3_20 - Gravel & Valin & Tanguay - initial version +* v3_21 - Desgagne M. - Revision Openmp +* +*object +* see id section +* +*arguments +*______________________________________________________________________ +* | | | +* NAME | DESCRIPTION | I/O | +*--------------|-------------------------------------------------|-----| +* F_out | result of interpolation | o | +* F_in | field to interpolate | i | +* | | | +* F_x | interpolation target X coordinate | i | +* F_y | interpolation target Y coordinate | i | +* F_z | interpolation target Z coordinate | i | +* | | | +* F_num | number of points to interpolate | i | +* | | | +* F_mono_L | switch: .true. : monotonic interpolation | i | +*______________|_________________________________________________|_____| +* +*implicits +#include "glb_ld.cdk" +#include "adw.cdk" +#include "adw_comp.cdk" +* ******************************************************************* +* Internal declarations +* ******************************************************************* +* Statement functions + real *8 triprd,za,zb,zc,zd + triprd(za,zb,zc,zd)=(za-zb)*(za-zc)*(za-zd) +* Standard declarations + integer :: nijk,err,m,cnt,nijag,i,j,k,nij,iimax,jjmax,kkmax + integer, dimension(3) :: timer0,timer + integer, dimension(:), allocatable, save :: ii,jj,kk,n + integer, dimension(:,:), allocatable, save :: o1,o2,o3,o4 + real, dimension(:), allocatable :: prmax, prmin + real(kind=8), dimension(:,:), allocatable :: a,b,c,d,p + real(kind=8), dimension(:,:), allocatable, save :: ra,rb,rc,rd + logical :: init=.true. + logical, dimension(:), allocatable, save :: zcubic_L +* +* ---------------------------------------------------- +* +#if defined (TIMER) +* Initialize system clock (-DTIMER) + call system_clock(count=timer0(1),count_rate=timer0(2),count_max=timer0(3)) +#endif +* +* ******************************************************************* +* Local grid sizing +* ******************************************************************* + nij = l_ni*l_nj + nijag = Adw_nit * Adw_njt + nijk = (kn-1)*nij + ((jn-1)*l_ni) + in + cnt = kn * (jn-j0+1) * (in-i0+1) +* + iimax = G_ni+2*Adw_halox-2 + jjmax = G_nj+Adw_haloy + kkmax = l_nk-1 +* +* ******************************************************************* +* Initialize local positional arrays on startup +* ******************************************************************* + if (init) then +* + allocate(ra(cnt,3),rb(cnt,3),rc(cnt,3),rd(cnt,3),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (D,S)" + allocate(ii(cnt),jj(cnt),kk(cnt),o1(cnt,4),o2(cnt,4),o3(cnt,4),o4(cnt,4),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (I,S)" + allocate(zcubic_L(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (L,S)" + allocate(n(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (I,S)" + m = 1 + do k=1,kn + do j=j0,jn + do i=i0,in + n(m) = (k-1)*nij + ((j-1)*l_ni) + i + m = m+1 + enddo + enddo + enddo +* + init = .false. !execute on startup only + adw_comp_cub_L = .true. !need to initialize values (below) +* + endif +* +* ******************************************************************* +* Set local array sizes and allocate +* ******************************************************************* + allocate(prmax(cnt),prmin(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (R,L)" + allocate(a(cnt,4),b(cnt,4),c(cnt,4),d(cnt,4),p(cnt,4),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (D,L)" +* +* ******************************************************************* +* Set positional parameters if required +* ******************************************************************* + if (adw_comp_cub_L) then + do m=1,cnt +* Prepare x points + ii(m) = ( F_x(n(m)) - Adw_x00_8 ) * Adw_ovdx_8 + ii(m) = Adw_lcx( ii(m)+1 ) + 1 + if (F_x(n(m)) < Adw_bsx_8(ii(m))) ii(m) = ii(m) - 1 + ii(m) = max(2,min(ii(m),iimax)) +* Prepare y points + jj(m) = ( F_y(n(m)) - Adw_y00_8 ) * Adw_ovdy_8 + jj(m) = Adw_lcy( jj(m)+1 ) + 1 + if (F_y(n(m)) < Adw_bsy_8(jj(m))) jj(m) = jj(m) - 1 + jj(m) = max(Adw_haloy,min(jj(m),jjmax)) +* Prepare z points + kk(m) = ( F_z(n(m)) - Adw_z00_8 ) * Adw_ovdz_8 + kk(m) = Adw_lcz( kk(m)+1 ) + if (F_z(n(m)) < Adw_bsz_8(kk(m))) kk(m) = kk(m) - 1 + kk(m) = min(kkmax-1,max(0,kk(m))) +* + zcubic_L(m) = (kk(m) > 0) .and. (kk(m) < kkmax-1) +* + o2(m,1) = (kk(m)-1)*nijag + (jj(m)-Adw_int_j_off-1)*Adw_nit + (ii(m)-Adw_int_i_off) + o1(m,1) = o2(m,1)-Adw_nit + o3(m,1) = o2(m,1)+Adw_nit + o4(m,1) = o3(m,1)+Adw_nit + o1(m,2) = o1(m,1)+nijag !unroll loop for vectorization + o2(m,2) = o2(m,1)+nijag + o3(m,2) = o3(m,1)+nijag + o4(m,2) = o4(m,1)+nijag + o1(m,3) = o1(m,2)+nijag + o2(m,3) = o2(m,2)+nijag + o3(m,3) = o3(m,2)+nijag + o4(m,3) = o4(m,2)+nijag + o1(m,4) = o1(m,3)+nijag + o2(m,4) = o2(m,3)+nijag + o3(m,4) = o3(m,3)+nijag + o4(m,4) = o4(m,3)+nijag +* + ra(m,1) = Adw_bsx_8(ii(m)-1) + rb(m,1) = Adw_bsx_8(ii(m) ) + rc(m,1) = Adw_bsx_8(ii(m)+1) + rd(m,1) = Adw_bsx_8(ii(m)+2) + ra(m,2) = Adw_bsy_8(jj(m)-1) + rb(m,2) = Adw_bsy_8(jj(m) ) + rc(m,2) = Adw_bsy_8(jj(m)+1) + rd(m,2) = Adw_bsy_8(jj(m)+2) + if (zcubic_L(m)) then + ra(m,3) = Adw_bsz_8(kk(m)-1) + rb(m,3) = Adw_bsz_8(kk(m) ) + rc(m,3) = Adw_bsz_8(kk(m)+1) + rd(m,3) = Adw_bsz_8(kk(m)+2) + endif + enddo +* + adw_comp_cub_L = .false. !save position until next request +* + endif +* +* ********************************************************************* +* Begin main interpolation loops +* ********************************************************************* +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + write(6,*) 'PRE-X time (ms): ',timer(1)-timer0(1) +#endif +* + do m=1,cnt +* +* ********************************************************************* +* x interpolation +* ********************************************************************* + p(m,1) = triprd(F_x(n(m)),rb(m,1),rc(m,1),rd(m,1))*Adw_xabcd_8(ii(m)) + p(m,2) = triprd(F_x(n(m)),ra(m,1),rc(m,1),rd(m,1))*Adw_xbacd_8(ii(m)) + p(m,3) = triprd(F_x(n(m)),ra(m,1),rb(m,1),rd(m,1))*Adw_xcabd_8(ii(m)) + p(m,4) = triprd(F_x(n(m)),ra(m,1),rb(m,1),rc(m,1))*Adw_xdabc_8(ii(m)) +* + if (zcubic_L(m)) then + a(m,1) = p(m,1) * F_in (o1(m,1)-1) + p(m,2) * F_in (o1(m,1)) + p(m,3) * F_in (o1(m,1)+1) + p(m,4) * F_in (o1(m,1)+2) + a(m,2) = p(m,1) * F_in (o2(m,1)-1) + p(m,2) * F_in (o2(m,1)) + p(m,3) * F_in (o2(m,1)+1) + p(m,4) * F_in (o2(m,1)+2) + a(m,3) = p(m,1) * F_in (o3(m,1)-1) + p(m,2) * F_in (o3(m,1)) + p(m,3) * F_in (o3(m,1)+1) + p(m,4) * F_in (o3(m,1)+2) + a(m,4) = p(m,1) * F_in (o4(m,1)-1) + p(m,2) * F_in (o4(m,1)) + p(m,3) * F_in (o4(m,1)+1) + p(m,4) * F_in (o4(m,1)+2) +* + d(m,1) = p(m,1) * F_in (o1(m,4)-1) + p(m,2) * F_in (o1(m,4)) + p(m,3) * F_in (o1(m,4)+1) + p(m,4) * F_in (o1(m,4)+2) + d(m,2) = p(m,1) * F_in (o2(m,4)-1) + p(m,2) * F_in (o2(m,4)) + p(m,3) * F_in (o2(m,4)+1) + p(m,4) * F_in (o2(m,4)+2) + d(m,3) = p(m,1) * F_in (o3(m,4)-1) + p(m,2) * F_in (o3(m,4)) + p(m,3) * F_in (o3(m,4)+1) + p(m,4) * F_in (o3(m,4)+2) + d(m,4) = p(m,1) * F_in (o4(m,4)-1) + p(m,2) * F_in (o4(m,4)) + p(m,3) * F_in (o4(m,4)+1) + p(m,4) * F_in (o4(m,4)+2) + endif +* + b(m,1) = p(m,1) * F_in (o1(m,2)-1) + p(m,2) * F_in (o1(m,2)) + p(m,3) * F_in (o1(m,2)+1) + p(m,4) * F_in (o1(m,2)+2) + b(m,2) = p(m,1) * F_in (o2(m,2)-1) + p(m,2) * F_in (o2(m,2)) + p(m,3) * F_in (o2(m,2)+1) + p(m,4) * F_in (o2(m,2)+2) + b(m,3) = p(m,1) * F_in (o3(m,2)-1) + p(m,2) * F_in (o3(m,2)) + p(m,3) * F_in (o3(m,2)+1) + p(m,4) * F_in (o3(m,2)+2) + b(m,4) = p(m,1) * F_in (o4(m,2)-1) + p(m,2) * F_in (o4(m,2)) + p(m,3) * F_in (o4(m,2)+1) + p(m,4) * F_in (o4(m,2)+2) +* + c(m,1) = p(m,1) * F_in (o1(m,3)-1) + p(m,2) * F_in (o1(m,3)) + p(m,3) * F_in (o1(m,3)+1) + p(m,4) * F_in (o1(m,3)+2) + c(m,2) = p(m,1) * F_in (o2(m,3)-1) + p(m,2) * F_in (o2(m,3)) + p(m,3) * F_in (o2(m,3)+1) + p(m,4) * F_in (o2(m,3)+2) + c(m,3) = p(m,1) * F_in (o3(m,3)-1) + p(m,2) * F_in (o3(m,3)) + p(m,3) * F_in (o3(m,3)+1) + p(m,4) * F_in (o3(m,3)+2) + c(m,4) = p(m,1) * F_in (o4(m,3)-1) + p(m,2) * F_in (o4(m,3)) + p(m,3) * F_in (o4(m,3)+1) + p(m,4) * F_in (o4(m,3)+2) +* + enddo +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + print*, 'PRE-Y time (ms): ',timer(1)-timer0(1) +#endif +* +* ********************************************************************* +* y interpolation +* ********************************************************************* + do m=1,cnt + p(m,1) = triprd(F_y(n(m)),rb(m,2),rc(m,2),rd(m,2))*Adw_yabcd_8(jj(m)) + p(m,2) = triprd(F_y(n(m)),ra(m,2),rc(m,2),rd(m,2))*Adw_ybacd_8(jj(m)) + p(m,3) = triprd(F_y(n(m)),ra(m,2),rb(m,2),rd(m,2))*Adw_ycabd_8(jj(m)) + p(m,4) = triprd(F_y(n(m)),ra(m,2),rb(m,2),rc(m,2))*Adw_ydabc_8(jj(m)) +* + if (zcubic_L(m)) then + a(m,1) = p(m,1) * a(m,1) + p(m,2) * a(m,2) + p(m,3) * a(m,3) + p(m,4) * a(m,4) + d(m,1) = p(m,1) * d(m,1) + p(m,2) * d(m,2) + p(m,3) * d(m,3) + p(m,4) * d(m,4) + endif +* + b(m,1) = p(m,1) * b(m,1) + p(m,2) * b(m,2) + p(m,3) * b(m,3) + p(m,4) * b(m,4) + c(m,1) = p(m,1) * c(m,1) + p(m,2) * c(m,2) + p(m,3) * c(m,3) + p(m,4) * c(m,4) + enddo +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + print*, 'PRE-Z time (ms): ',timer(1)-timer0(1) +#endif +* +* ********************************************************************* +* z interpolation +* ********************************************************************* + do m=1,cnt + if (zcubic_L(m)) then + p(m,1) = triprd(F_z(n(m)),rb(m,3),rc(m,3),rd(m,3))*Adw_zabcd_8(kk(m)+1) + p(m,2) = triprd(F_z(n(m)),ra(m,3),rc(m,3),rd(m,3))*Adw_zbacd_8(kk(m)+1) + p(m,3) = triprd(F_z(n(m)),ra(m,3),rb(m,3),rd(m,3))*Adw_zcabd_8(kk(m)+1) + p(m,4) = triprd(F_z(n(m)),ra(m,3),rb(m,3),rc(m,3))*Adw_zdabc_8(kk(m)+1) +* + p(m,1) = p(m,1) * a(m,1) + p(m,2) * b(m,1) + p(m,3) * c(m,1) + p(m,4) * d(m,1) !recycle p1 +* + else + p(m,3) = (F_z(n(m))-Adw_bsz_8(kk(m)))*Adw_zbc_8(kk(m)+1) + p(m,2) = 1. - p(m,3) +* + p(m,1) = p(m,2) * b(m,1) + p(m,3) * c(m,1) !recycle p1 + endif +* +* ********************************************************************* +* End of main loops +* ********************************************************************* + enddo +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + print*, 'PRE-ASSIGNMENT time (ms): ',timer(1)-timer0(1) +#endif +* +* ********************************************************************* +* Final data assignment +* ********************************************************************* + if (.not.F_mono_L) then + do m=1,cnt + F_out(n(m)) = p(m,1) + enddo + else + do m=1,cnt + prmax(m) = max(F_in(o2(m,2)),F_in(o2(m,2)+1),F_in(o3(m,2)),F_in(o3(m,2)+1)) + prmin(m) = min(F_in(o2(m,2)),F_in(o2(m,2)+1),F_in(o3(m,2)),F_in(o3(m,2)+1)) + prmax(m) = max(prmax(m),F_in(o2(m,3)),F_in(o2(m,3)+1),F_in(o3(m,3)),F_in(o3(m,3)+1)) + prmin(m) = min(prmin(m),F_in(o2(m,3)),F_in(o2(m,3)+1),F_in(o3(m,3)),F_in(o3(m,3)+1)) + F_out(n(m)) = max ( prmin(m) , min(prmax(m),p(m,1)) ) + enddo + endif +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + print*, 'POST-ASSIGNMENT time (ms): ',timer(1)-timer0(1) +#endif +* +* ********************************************************************* +* Clear stack +* ********************************************************************* + deallocate(prmin,prmax) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Deallocation error (R,L)" + deallocate(a,b,c,d,p,stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Deallocation error (D,L)" +* +#if defined (TIMER) +* System timing function (-DTIMER) + call system_clock(count=timer(1),count_rate=timer(2),count_max=timer(3)) + print*, 'FINAL WALL-CLOCK time (ms): ',timer(1)-timer0(1) +#endif +* + return + end Index: adw_tritrunc_lag3d.ftn =================================================================== --- adw_tritrunc_lag3d.ftn 2006-06-28 14:36:58.000000000 -0400 +++ adw_tritrunc_lag3d.ftn 2006-06-28 14:31:07.000000000 -0400 @@ -0,0 +1,299 @@ +***s/r adw_tritrunc_lag3d - Tri-Lagrangian (truncated) interpolation. +* +#include +* + subroutine adw_tritrunc_lag3d ( F_out, F_in, F_x, F_y, F_z, + % F_num, F_mono_L, i0, in, j0, jn, kn ) +* + implicit none +* + logical F_mono_L +* + integer F_num, i0, in, j0, jn, kn +* + real F_in(*) +* + real F_out (F_num), F_x(F_num), F_y(F_num), F_z(F_num) +* +*authors +* McTaggart-Cowan +* +* (Based on adw_tricub_lag3d and adw_trilin v_3.2.1 with modifications +* as per ECMWF http://www.ecmwf.int/research/ifsdocs/CY28r1/Dynamics/Dynamics-4-02.html) +* +*revision +* v3_30 - McTaggart-Cowan - initial version +* +*object +* see id section +* +*arguments +*______________________________________________________________________ +* | | | +* NAME | DESCRIPTION | I/O | +*--------------|-------------------------------------------------|-----| +* F_out | result of interpolation | o | +* F_in | field to interpolate | i | +* | | | +* F_x | interpolation target X coordinate | i | +* F_y | interpolation target Y coordinate | i | +* F_z | interpolation target Z coordinate | i | +* | | | +* F_num | number of points to interpolate | i | +* | | | +* F_mono_L | switch: .true. : monotonic interpolation | i | +*______________|_________________________________________________|_____| +* +*implicits +#include "glb_ld.cdk" +#include "adw.cdk" +* +*notes +* This algorithm is a truncated version of the full 3D Lagrangian +* interpolation procedure (adw_trilag_3d). Full Lagrangian interpolation +* requires that the local calculation be done in a 4-point 3D cube, +* thereby needing 64 values for each point. This imposes an enormous +* load on memory access during the gather/scatter operation. In this +* truncation, a 3D diamond rather than a cube is required for the +* interpolation operation (32 values). Interpolation to the points closest to the +* back-trajectory origin is done using cubic function; however, interpolation +* to points further from the origin are done linearly. As a result, +* each 3D truncated interpolation uses only 7 Lagragian interpolations +* and 10 linear interpolations (compared to 21 Lagrangian interpolations +* in the full 3D Lagrangian algorithm). +* +* The organization of the diamond is shown in plan form here, along with +* the order of each interpolation for the 'inner' layers (immediately above +* and below the point of interest, indeces k and k+1), and 'outer' layers +* (indeces k-1 and k+2). The origin of the back-trajectory is denoted with +* {} braces around the interpolation method used in the y-direction to +* obtain it. Both layers are plotted in the horizontal plane. Cube points +* not accessed are denoted with a '0', and those addressed and used by the +* truncated algorithm are denoted with an 'X'. +* +* Inner layers 2 x (2 linear; 3 cubic): +* 0 X -- linear -- X 0 +* | +* | +* | +* X ------------- X --- cubic --- X ------------- X +* | +* {cubic} +* | +* X ------------- X --- cubic --- X ------------- X +* | +* | +* | +* 0 X -- linear -- X 0 +* +* Outer layers 2 x (3 linear): +* 0 0 0 0 +* +* +* +* 0 X -- linear -- X 0 +* | +* {linear} +* | +* 0 X -- linear -- X 0 +* +* +* +* 0 0 0 0 +* +* The vertical interpolation is cubic (Lagrangian) through the +* four points obtained using the layers shown above. This interpolation +* constitues the 7th (and final) higher order interpolation performed +* by the truncated algorithm. +* +* ********************************************************************** + integer n,nijag,i,j,k,nij,iimax,jjmax,kkmax,ii,jj,kk + logical zcubic_L +* + + + integer count + + + real prmin, prmax +* + integer o1, o2, o3, o4 + + real*8 a2, a3 + real*8 b1, b2, b3, b4 + real*8 c1, c2, c3, c4 + real*8 d2, d3 + + real*8 p1, p2, p3, p4 +* + real*8 :: capx,capy,capz + real*8 :: triprd,za,zb,zc,zd,rri,rrj,rrk,ra,rb,rc,rd + triprd(za,zb,zc,zd)=(za-zb)*(za-zc)*(za-zd) +* +* ---------------------------------------------------- +* + nij = l_ni*l_nj + nijag = Adw_nit * Adw_njt +* + iimax = G_ni+2*Adw_halox-2 + jjmax = G_nj+Adw_haloy + kkmax = l_nk-1 +* + + + + count = 0 + + + + +* +c!$omp parallel private(n,i,j,k, +c!$omp& ii,jj,kk,capx,capy,capz,zcubic_L, +c!$omp& o1, o2, o3, o4, a2, a3, b1, b2, b3, b4, +c!$omp& c1, c2, c3, c4, d2, d3, p1, p2, p3, p4, +c!$omp& rri,rrj,rrk,ra,rb,rc,rd,prmin,prmax) +c!$omp do + do 100 k=1,kn + do 90 j=j0,jn + do 80 i=i0,in + + count = count+1 + + + n = (k-1)*nij + ((j-1)*l_ni) + i +* + rri= F_x(n) + ii = ( rri - Adw_x00_8 ) * Adw_ovdx_8 + ii = Adw_lcx( ii+1 ) + 1 + if ( rri .lt. Adw_bsx_8(ii) ) ii = ii - 1 + ii = max(2,min(ii,iimax)) + capx = (rri - Adw_bsx_8(ii)) / (Adw_bsx_8(ii+1) - Adw_bsx_8(ii)) +* + rrj= F_y(n) + jj = ( rrj - Adw_y00_8 ) * Adw_ovdy_8 + jj = Adw_lcy( jj+1 ) + 1 + if ( rrj .lt. Adw_bsy_8(jj) ) jj = jj - 1 + jj = max(Adw_haloy,min(jj,jjmax)) + capy = (rrj - Adw_bsy_8(jj)) / (Adw_bsy_8(jj+1) - Adw_bsy_8(jj)) +* + rrk= F_z(n) + kk = ( rrk - Adw_z00_8 ) * Adw_ovdz_8 + kk = Adw_lcz( kk+1 ) + if ( rrk .lt. Adw_bsz_8(kk) ) kk = kk - 1 + kk = min(kkmax-1,max(0,kk)) + capz = (rrk - Adw_bsz_8(kk)) / (Adw_bsz_8(kk+1) - Adw_bsz_8(kk)) +* + zcubic_L = (kk.gt.0) .and. (kk.lt.kkmax-1) +* +* ********************************************************************* +* x interpolation +* ********************************************************************* + ra = Adw_bsx_8(ii-1) + rb = Adw_bsx_8(ii ) + rc = Adw_bsx_8(ii+1) + rd = Adw_bsx_8(ii+2) + p1 = triprd(rri,rb,rc,rd)*Adw_xabcd_8(ii) + p2 = triprd(rri,ra,rc,rd)*Adw_xbacd_8(ii) + p3 = triprd(rri,ra,rb,rd)*Adw_xcabd_8(ii) + p4 = triprd(rri,ra,rb,rc)*Adw_xdabc_8(ii) +* + o2 = (kk-1)*nijag + (jj-Adw_int_j_off-1)*Adw_nit + (ii-Adw_int_i_off) + o3 = o2+Adw_nit +* + if(zcubic_L) then + a2 = (1.d0 - capx) * F_in(o2) + capx * F_in(o2+1) + a3 = (1.d0 - capx) * F_in(o3) + capx * F_in(o3+1) + else + a2 = 0.d0 + a3 = 0.d0 + endif +* + o2 = o2 + nijag + o3 = o3 + nijag + o1 = o2-Adw_nit + o4 = o3+Adw_nit +* + if (F_mono_L) then + prmax = max(F_in(o2),F_in(o2+1),F_in(o3),F_in(o3+1)) + prmin = min(F_in(o2),F_in(o2+1),F_in(o3),F_in(o3+1)) + else + prmax = 0.d0 + prmin = 0.d0 + endif + b1 = (1.d0 - capx) * F_in(o1) + capx * F_in(o1+1) + b2 = p1 * F_in (o2-1) + p2 * F_in (o2) + p3 * F_in (o2+1) + p4 * F_in (o2+2) + b3 = p1 * F_in (o3-1) + p2 * F_in (o3) + p3 * F_in (o3+1) + p4 * F_in (o3+2) + b4 = (1.d0 - capx) * F_in(o4) + capx * F_in(o4+1) +* + o1 = o1 + nijag + o2 = o2 + nijag + o3 = o3 + nijag + o4 = o4 + nijag +* + if (F_mono_L) then + prmax = max(prmax,F_in(o2),F_in(o2+1),F_in(o3),F_in(o3+1)) + prmin = min(prmin,F_in(o2),F_in(o2+1),F_in(o3),F_in(o3+1)) + endif + c1 = (1.d0 - capx) * F_in(o1) + capx * F_in(o1+1) + c2 = p1 * F_in (o2-1) + p2 * F_in (o2) + p3 * F_in (o2+1) + p4 * F_in (o2+2) + c3 = p1 * F_in (o3-1) + p2 * F_in (o3) + p3 * F_in (o3+1) + p4 * F_in (o3+2) + c4 = (1.d0 - capx) * F_in(o4) + capx * F_in(o4+1) +* + o2 = o2 + nijag + o3 = o3 + nijag +* + if(zcubic_L) then + d2 = (1.d0 - capx) * F_in(o2) + capx * F_in(o2+1) + d3 = (1.d0 - capx) * F_in(o3) + capx * F_in(o3+1) + else + d2 = 0.d0 + d3 = 0.d0 + endif +* ********************************************************************* +* y interpolation +* ********************************************************************* + ra = Adw_bsy_8(jj-1) + rb = Adw_bsy_8(jj ) + rc = Adw_bsy_8(jj+1) + rd = Adw_bsy_8(jj+2) + p1 = triprd(rrj,rb,rc,rd)*Adw_yabcd_8(jj) + p2 = triprd(rrj,ra,rc,rd)*Adw_ybacd_8(jj) + p3 = triprd(rrj,ra,rb,rd)*Adw_ycabd_8(jj) + p4 = triprd(rrj,ra,rb,rc)*Adw_ydabc_8(jj) +* + if (zcubic_L) a2 = (1.d0 - capy) * a2 + capy * a3 + b1 = p1 * b1 + p2 * b2 + p3 * b3 + p4 * b4 + c1 = p1 * c1 + p2 * c2 + p3 * c3 + p4 * c4 + if (zcubic_L) d2 = (1.d0 - capy) * d2 + capy * d3 +* ********************************************************************* +* z interpolation +* ********************************************************************* + if(zcubic_L) then + ra = Adw_bsz_8(kk-1) + rb = Adw_bsz_8(kk ) + rc = Adw_bsz_8(kk+1) + rd = Adw_bsz_8(kk+2) + p1 = triprd(rrk,rb,rc,rd)*Adw_zabcd_8(kk+1) + p2 = triprd(rrk,ra,rc,rd)*Adw_zbacd_8(kk+1) + p3 = triprd(rrk,ra,rb,rd)*Adw_zcabd_8(kk+1) + p4 = triprd(rrk,ra,rb,rc)*Adw_zdabc_8(kk+1) +* + F_out(n) = p1 * a2 + p2 * b1 + p3 * c1 + p4 * d2 +* + else +* + F_out(n) = (1.d0 - capz) * b1 + capz * c1 +* + endif +* + if (F_mono_L) F_out(n) = max ( prmin , min(prmax,F_out(n)) ) +* + 80 continue + 90 continue + 100 continue +c!$omp enddo +c!$omp end parallel +* + return + end Index: adw_comp.cdk =================================================================== --- adw_comp.cdk 2006-06-28 14:37:09.000000000 -0400 +++ adw_comp.cdk 2006-06-28 14:36:17.000000000 -0400 @@ -0,0 +1,29 @@ +#if defined (DOC) +* +* +*revision +* v3_21 - McTaggart-Cowan R. - initial version +* +***comdeck adw_comp.cdk +* +*_______________________________________________________________ +* | +* FLAGS TO HANDLE RE-COMPUTATION OF POSITIONAL COEFFICIENTS | +*_______________________________________________________________ +* | | +* NAME | | +*---------------|-----------------------------------------------| +* adw_comp_cub_L| .true. -> recompute positional coefficients | +* | for cubic interpolation | +* | .false.-> use saved positional coefficients | +* adw_comp_lin_L| .true. -> recompute positional coefficients | +* | for linear interpolation | +*--------------------------------------------------------------- +* +#endif + + logical adw_comp_cub_L, adw_comp_lin_L + + MARK_COMMON_BEG (adw_comp) + common / adw_comp / adw_comp_cub_L, adw_comp_lin_L + MARK_COMMON_END (adw_comp) Index: adw_main_3_int.ftn =================================================================== --- adw_main_3_int.ftn 2006-06-28 14:35:14.000000000 -0400 +++ adw_main_3_int.ftn 2006-06-28 15:28:37.000000000 -0400 @@ -81,7 +81,7 @@ pointer (patr, tr(LDIST_SHAPE,*)),(patr0, tr0(LDIST_SHAPE,*)) real, allocatable, dimension(:,:,:,:) :: tr2 * - if ( Adw_lag3d_L ) then + if (Adw_interp_type(1:3).eq.'lag' .or. Adw_interp_type(1:3).eq.'tru') then call adw_main_3_intlag ( F_u, F_v, F_w ) return endif @@ -155,11 +155,6 @@ * endif * - if( .not. Adw_lag3d_L ) then - call adw_setint ( Adw_n1,Adw_capx1,Adw_xgg1,Adw_xdd1,Adw_capy1,Adw_ygg1, - % Adw_ydd1, Adw_capz1, Adw_cz1, F_u, F_v, F_w, - % .true., .true., .false., nijk,i0,in,j0,jn,l_nk) - else * ---------------------------- * Keep positions in CAP fields * ---------------------------- @@ -175,13 +170,12 @@ end do end do !$omp end parallel do - endif * if ( Adw_fro_a .gt. 0 .and. .not. G_lam) then * if ( Adw_ckbd_L ) call adw_ckbd ( Adw_capy2 ) * - if( .not. Adw_lag3d_L ) then + if( Adw_interp_type(1:3) .eq. 'cub' ) then call adw_setint(Adw_n2,Adw_capx2,Adw_xgg2,Adw_xdd2,Adw_capy2,Adw_ygg2, % Adw_ydd2,Adw_capz2,Adw_cz2,Adw_capx2,Adw_capy2,Adw_capz2, % .true., .true., .false., Adw_fro_a,1,Adw_fro_a,1,1,1) Index: mtn_cfg.ftn =================================================================== --- mtn_cfg.ftn 2006-06-28 14:35:15.000000000 -0400 +++ mtn_cfg.ftn 2006-06-28 14:31:08.000000000 -0400 @@ -106,7 +106,7 @@ c Offc_b1_8 = 0.5 Zblen_L = .true. Zblen_spngtt_L=.false. - Adw_lag3d_L = .true. + Adw_interp_type = 'lagrangian' if ( Theo_case_S.eq.'MTN_SCHAR') then Index: adw_tritrunc_lag3d_vec.ftn =================================================================== --- adw_tritrunc_lag3d_vec.ftn 2006-06-28 14:37:02.000000000 -0400 +++ adw_tritrunc_lag3d_vec.ftn 2006-06-28 14:31:07.000000000 -0400 @@ -0,0 +1,340 @@ +***s/r adw_tritrunc_lag3d_vec - Tri-Lagrangian (truncated) interpolation. +* +#include +* + subroutine adw_tritrunc_lag3d_vec ( F_out, F_in, F_x, F_y, F_z, + % F_num, F_mono_L, i0, in, j0, jn, kn ) +* + implicit none +* + logical F_mono_L +* + integer F_num, i0, in, j0, jn, kn +* + real F_in(*) +* + real F_out (F_num), F_x(F_num), F_y(F_num), F_z(F_num) +* +*authors +* McTaggart-Cowan +* +* (Based on adw_tricub_lag3d and adw_trilin v_3.2.1 with modifications +* as per ECMWF http://www.ecmwf.int/research/ifsdocs/CY28r1/Dynamics/Dynamics-4-02.html) +* +*revision +* v3_30 - McTaggart-Cowan - initial version +* +*object +* see id section +* +*arguments +*______________________________________________________________________ +* | | | +* NAME | DESCRIPTION | I/O | +*--------------|-------------------------------------------------|-----| +* F_out | result of interpolation | o | +* F_in | field to interpolate | i | +* | | | +* F_x | interpolation target X coordinate | i | +* F_y | interpolation target Y coordinate | i | +* F_z | interpolation target Z coordinate | i | +* | | | +* F_num | number of points to interpolate | i | +* | | | +* F_mono_L | switch: .true. : monotonic interpolation | i | +*______________|_________________________________________________|_____| +* +*implicits +#include "glb_ld.cdk" +#include "adw.cdk" +#include "adw_comp.cdk" +* +*notes +* This algorithm is a truncated version of the full 3D Lagrangian +* interpolation procedure (adw_trilag_3d). Full Lagrangian interpolation +* requires that the local calculation be done in a 4-point 3D cube, +* thereby needing 64 values for each point. This imposes an enormous +* load on memory access during the gather/scatter operation. In this +* truncation, a 3D diamond rather than a cube is required for the +* interpolation operation (32 values). Interpolation to the points closest to the +* back-trajectory origin is done using cubic function; however, interpolation +* to points further from the origin are done linearly. As a result, +* each 3D truncated interpolation uses only 7 Lagragian interpolations +* and 10 linear interpolations (compared to 21 Lagrangian interpolations +* in the full 3D Lagrangian algorithm). +* +* The organization of the diamond is shown in plan form here, along with +* the order of each interpolation for the 'inner' layers (immediately above +* and below the point of interest, indeces k and k+1), and 'outer' layers +* (indeces k-1 and k+2). The origin of the back-trajectory is denoted with +* {} braces around the interpolation method used in the y-direction to +* obtain it. Both layers are plotted in the horizontal plane. Cube points +* not accessed are denoted with a '0', and those addressed and used by the +* truncated algorithm are denoted with an 'X'. +* +* Inner layers 2 x (2 linear; 3 cubic): +* 0 X -- linear -- X 0 +* | +* | +* | +* X ------------- X --- cubic --- X ------------- X +* | +* {cubic} +* | +* X ------------- X --- cubic --- X ------------- X +* | +* | +* | +* 0 X -- linear -- X 0 +* +* Outer layers 2 x (3 linear): +* 0 0 0 0 +* +* +* +* 0 X -- linear -- X 0 +* | +* {linear} +* | +* 0 X -- linear -- X 0 +* +* +* +* 0 0 0 0 +* +* The vertical interpolation is cubic (Lagrangian) through the +* four points obtained using the layers shown above. This interpolation +* constitues the 7th (and final) higher order interpolation performed +* by the truncated algorithm. +* +* ********************************************************************** +* Internal declarations +* ********************************************************************** +* Statement functions + real(kind=8) :: triprd,za,zb,zc,zd + triprd(za,zb,zc,zd)=(za-zb)*(za-zc)*(za-zd) +* Standard declarations + integer :: nijk,nijag,i,j,k,nij,iimax,jjmax,kkmax,cnt,err,m + integer, dimension(:), allocatable, save :: ii,jj,kk,n + integer, dimension(:,:), allocatable, save :: o1,o2,o3,o4 + real, dimension(:), allocatable :: prmax,prmin + real(kind=8),dimension(:,:), allocatable, save :: cap + real(kind=8), dimension(:,:), allocatable :: a,b,c,d,p + real(kind=8), dimension(:,:), allocatable, save :: ra,rb,rc,rd + logical :: init=.true. + logical, dimension(:), allocatable, save :: zcubic_L +* +* ---------------------------------------------------- +* +* ********************************************************************** +* Local grid sizing +* ********************************************************************** + nij = l_ni*l_nj + nijag = Adw_nit * Adw_njt + nijk = (kn-1)*nij + ((jn-1)*l_ni) + in + cnt = kn * (jn-j0+1) * (in-i0+1) +* + iimax = G_ni+2*Adw_halox-2 + jjmax = G_nj+Adw_haloy + kkmax = l_nk-1 +* +* ********************************************************************** +* Initialize local positional arrays on startup +* ********************************************************************** + if (init) then +* + allocate(ra(cnt,3),rb(cnt,3),rc(cnt,3),rd(cnt,3),cap(cnt,3),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (D,S)" + allocate(ii(cnt),jj(cnt),kk(cnt),o1(cnt,4),o2(cnt,4),o3(cnt,4),o4(cnt,4),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (I,S)" + allocate(zcubic_L(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (L,S)" + allocate(n(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (I,S)" + m = 1 + do k=1,kn + do j=j0,jn + do i=i0,in + n(m) = (k-1)*nij + ((j-1)*l_ni) + i + m = m+1 + enddo + enddo + enddo +* + init = .false. !execute on startup only + adw_comp_cub_L = .true. !need to initialize values (below) +* + endif +* +* ********************************************************************** +* Set local array sizes and allocate +* ******************************************************************* + allocate(prmax(cnt),prmin(cnt),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (R,L)" + allocate(a(cnt,2:3),b(cnt,4),c(cnt,4),d(cnt,2:3),p(cnt,4),stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Allocation error (D,L)" +* +* ******************************************************************* +* Set positional parameters if required +* ******************************************************************* + if (adw_comp_cub_L) then + do m=1,cnt +* Prepare x points + ii(m) = ( F_x(n(m)) - Adw_x00_8 ) * Adw_ovdx_8 + ii(m) = Adw_lcx( ii(m)+1 ) + 1 + if (F_x(n(m)) < Adw_bsx_8(ii(m))) ii(m) = ii(m) - 1 + ii(m) = max(2,min(ii(m),iimax)) + cap(m,1) = (F_x(n(m)) - Adw_bsx_8(ii(m))) / (Adw_bsx_8(ii(m)+1) - Adw_bsx_8(ii(m))) +* Prepare y points + jj(m) = ( F_y(n(m)) - Adw_y00_8 ) * Adw_ovdy_8 + jj(m) = Adw_lcy( jj(m)+1 ) + 1 + if (F_y(n(m)) < Adw_bsy_8(jj(m))) jj(m) = jj(m) - 1 + jj(m) = max(Adw_haloy,min(jj(m),jjmax)) + cap(m,2) = (F_y(n(m)) - Adw_bsy_8(jj(m))) / (Adw_bsy_8(jj(m)+1) - Adw_bsy_8(jj(m))) +* Prepare z points + kk(m) = ( F_z(n(m)) - Adw_z00_8 ) * Adw_ovdz_8 + kk(m) = Adw_lcz( kk(m)+1 ) + if (F_z(n(m)) < Adw_bsz_8(kk(m))) kk(m) = kk(m) - 1 + kk(m) = min(kkmax-1,max(0,kk(m))) + cap(m,3) = (F_z(n(m)) - Adw_bsz_8(kk(m))) / (Adw_bsz_8(kk(m)+1) - Adw_bsz_8(kk(m))) +* + zcubic_L(m) = (kk(m) > 0) .and. (kk(m) < kkmax-1) +* + o2(m,1) = (kk(m)-1)*nijag + (jj(m)-Adw_int_j_off-1)*Adw_nit + (ii(m)-Adw_int_i_off) + o1(m,1) = o2(m,1)-Adw_nit + o3(m,1) = o2(m,1)+Adw_nit + o4(m,1) = o3(m,1)+Adw_nit + o1(m,2) = o1(m,1)+nijag !unroll loop for vectorization + o2(m,2) = o2(m,1)+nijag + o3(m,2) = o3(m,1)+nijag + o4(m,2) = o4(m,1)+nijag + o1(m,3) = o1(m,2)+nijag + o2(m,3) = o2(m,2)+nijag + o3(m,3) = o3(m,2)+nijag + o4(m,3) = o4(m,2)+nijag + o1(m,4) = o1(m,3)+nijag + o2(m,4) = o2(m,3)+nijag + o3(m,4) = o3(m,3)+nijag + o4(m,4) = o4(m,3)+nijag +* + ra(m,1) = Adw_bsx_8(ii(m)-1) + rb(m,1) = Adw_bsx_8(ii(m) ) + rc(m,1) = Adw_bsx_8(ii(m)+1) + rd(m,1) = Adw_bsx_8(ii(m)+2) + ra(m,2) = Adw_bsy_8(jj(m)-1) + rb(m,2) = Adw_bsy_8(jj(m) ) + rc(m,2) = Adw_bsy_8(jj(m)+1) + rd(m,2) = Adw_bsy_8(jj(m)+2) + if (zcubic_L(m)) then + ra(m,3) = Adw_bsz_8(kk(m)-1) + rb(m,3) = Adw_bsz_8(kk(m) ) + rc(m,3) = Adw_bsz_8(kk(m)+1) + rd(m,3) = Adw_bsz_8(kk(m)+2) + endif + enddo +* + adw_comp_cub_L = .false. !save position until next request +* + endif +* +* ********************************************************************* +* Begin main interpolation loops +* ********************************************************************* + do m=1,cnt +* +* ********************************************************************* +* x interpolation +* ********************************************************************* + p(m,1) = triprd(F_x(n(m)),rb(m,1),rc(m,1),rd(m,1))*Adw_xabcd_8(ii(m)) + p(m,2) = triprd(F_x(n(m)),ra(m,1),rc(m,1),rd(m,1))*Adw_xbacd_8(ii(m)) + p(m,3) = triprd(F_x(n(m)),ra(m,1),rb(m,1),rd(m,1))*Adw_xcabd_8(ii(m)) + p(m,4) = triprd(F_x(n(m)),ra(m,1),rb(m,1),rc(m,1))*Adw_xdabc_8(ii(m)) +* + if (zcubic_L(m)) then + a(m,2) = (1.d0 - cap(m,1)) * F_in(o2(m,1)) + cap(m,1) * F_in(o2(m,1)+1) + a(m,3) = (1.d0 - cap(m,1)) * F_in(o3(m,1)) + cap(m,1) * F_in(o3(m,1)+1) +* + d(m,2) = (1.d0 - cap(m,1)) * F_in(o2(m,4)) + cap(m,1) * F_in(o2(m,4)+1) + d(m,3) = (1.d0 - cap(m,1)) * F_in(o3(m,4)) + cap(m,1) * F_in(o3(m,4)+1) + endif +* + b(m,1) = (1.d0 - cap(m,1)) * F_in(o1(m,2)) + cap(m,1) * F_in(o1(m,2)+1) + b(m,2) = p(m,1) * F_in (o2(m,2)-1) + p(m,2) * F_in (o2(m,2)) + p(m,3) * F_in (o2(m,2)+1) + p(m,4) * F_in (o2(m,2)+2) + b(m,3) = p(m,1) * F_in (o3(m,2)-1) + p(m,2) * F_in (o3(m,2)) + p(m,3) * F_in (o3(m,2)+1) + p(m,4) * F_in (o3(m,2)+2) + b(m,4) = (1.d0 - cap(m,1)) * F_in(o4(m,2)) + cap(m,1) * F_in(o4(m,2)+1) +* + c(m,1) = (1.d0 - cap(m,1)) * F_in(o1(m,3)) + cap(m,1) * F_in(o1(m,3)+1) + c(m,2) = p(m,1) * F_in (o2(m,3)-1) + p(m,2) * F_in (o2(m,3)) + p(m,3) * F_in (o2(m,3)+1) + p(m,4) * F_in (o2(m,3)+2) + c(m,3) = p(m,1) * F_in (o3(m,3)-1) + p(m,2) * F_in (o3(m,3)) + p(m,3) * F_in (o3(m,3)+1) + p(m,4) * F_in (o3(m,3)+2) + c(m,4) = (1.d0 - cap(m,1)) * F_in(o4(m,3)) + cap(m,1) * F_in(o4(m,3)+1) +* + enddo +* +* ********************************************************************* +* y interpolation +* ********************************************************************* + do m=1,cnt + p(m,1) = triprd(F_y(n(m)),rb(m,2),rc(m,2),rd(m,2))*Adw_yabcd_8(jj(m)) + p(m,2) = triprd(F_y(n(m)),ra(m,2),rc(m,2),rd(m,2))*Adw_ybacd_8(jj(m)) + p(m,3) = triprd(F_y(n(m)),ra(m,2),rb(m,2),rd(m,2))*Adw_ycabd_8(jj(m)) + p(m,4) = triprd(F_y(n(m)),ra(m,2),rb(m,2),rc(m,2))*Adw_ydabc_8(jj(m)) +* + if (zcubic_L(m)) then + a(m,2) = (1.d0 - cap(m,2)) * a(m,2) + cap(m,2) * a(m,3) + d(m,2) = (1.d0 - cap(m,2)) * d(m,2) + cap(m,2) * d(m,3) + endif +* + b(m,1) = p(m,1) * b(m,1) + p(m,2) * b(m,2) + p(m,3) * b(m,3) + p(m,4) * b(m,4) + c(m,1) = p(m,1) * c(m,1) + p(m,2) * c(m,2) + p(m,3) * c(m,3) + p(m,4) * c(m,4) + enddo +* +* ********************************************************************* +* z interpolation +* ********************************************************************* + do m=1,cnt + if (zcubic_L(m)) then + p(m,1) = triprd(F_z(n(m)),rb(m,3),rc(m,3),rd(m,3))*Adw_zabcd_8(kk(m)+1) + p(m,2) = triprd(F_z(n(m)),ra(m,3),rc(m,3),rd(m,3))*Adw_zbacd_8(kk(m)+1) + p(m,3) = triprd(F_z(n(m)),ra(m,3),rb(m,3),rd(m,3))*Adw_zcabd_8(kk(m)+1) + p(m,4) = triprd(F_z(n(m)),ra(m,3),rb(m,3),rc(m,3))*Adw_zdabc_8(kk(m)+1) +* + p(m,1) = p(m,1) * a(m,2) + p(m,2) * b(m,1) + p(m,3) * c(m,1) + p(m,4) * d(m,2) !recycle p1 +* + else +* + p(m,1) = (1.d0 - cap(m,3)) * b(m,1) + cap(m,3) * c(m,1) !recycle p1 +* + endif +* +* ********************************************************************* +* End of main loops +* ********************************************************************* + enddo +* +* ********************************************************************* +* Final data assignment +* ********************************************************************* + if (.not.F_mono_L) then + do m=1,cnt + F_out(n(m)) = p(m,1) + enddo + else + do m=1,cnt + prmax(m) = max(F_in(o2(m,2)),F_in(o2(m,2)+1),F_in(o3(m,2)),F_in(o3(m,2)+1)) + prmin(m) = min(F_in(o2(m,2)),F_in(o2(m,2)+1),F_in(o3(m,2)),F_in(o3(m,2)+1)) + prmax(m) = max(prmax(m),F_in(o2(m,3)),F_in(o2(m,3)+1),F_in(o3(m,3)),F_in(o3(m,3)+1)) + prmin(m) = min(prmin(m),F_in(o2(m,3)),F_in(o2(m,3)+1),F_in(o3(m,3)),F_in(o3(m,3)+1)) + F_out(n(m)) = max ( prmin(m) , min(prmax(m),p(m,1)) ) + enddo + endif +* +* ********************************************************************* +* Clear stack +* ********************************************************************* + deallocate(prmin,prmax) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Deallocation error (R,L)" + deallocate(a,b,c,d,p,stat=err) + if (err /= 0) write(6,*) "ADW_TRICUB_LAG3D_VEC: Deallocation error (D,L)" +* + return + end Index: set_world_view1.ftn =================================================================== --- set_world_view1.ftn 2006-06-28 14:35:15.000000000 -0400 +++ set_world_view1.ftn 2006-06-28 15:28:37.000000000 -0400 @@ -58,6 +58,7 @@ * - Move read of namelist coupling to c_init_cplfdl * v3_21 - Lee V. - Added Lctl_reset key, remove TR2D * v3_22 - Lee V. - Add Pres_vtap_L key for vtap +* v3_30 - McTaggart-Cowan R.- Add interpolator option Adw_interp * *object * To initialize model by reading namelist and return the @@ -198,7 +199,7 @@ Adw_exdg_L = .false. Adw_ckbd_L = .false. Adw_mono_L = .true. - Adw_lag3d_L= .false. + Adw_interp_type = 'lagrangian' Adw_nosetint_L = .false. Adw_trunc_traj_L = .false. Adw_halox = 3 @@ -550,7 +551,19 @@ * Lun_sortie_s = Path_dfwmil_S(1:longueur(Path_dfwmil_S)) $ //'/process/output_settings' - +* +* Check for valid back-trajectory interpolation scheme +* + select case (trim(Adw_interp_type)) + case('lagrangian') + continue + case('truncated') + continue + case DEFAULT + write (Lun_out, 9500) trim(Adw_interp_type) + write (Lun_out, 8000) + return + end select * * * Options used for the Off-line mode (MEC) @@ -632,6 +645,9 @@ 9300 format (/,' sensib_L twin_L 4dvar_L sgvc_L when V4dg_conf/100=1 only'/) 9400 format (/,' Step_total lt Init_dfnp-1 and Init_balgm_L when V4dg_conf.ne.0 NOT VALID '/) 9401 format (/,' Schm_theoc_L when V4dg_conf.ne.0 NOT VALID '/) + 9500 format (/,' INVALID back-trajectory interpolator'/ + $ ' (Adw_interp_type = ',a,'). Valid interpolator'/ + $ ' types are "lagrangian", "truncated".') * *------------------------------------------------------------------- * Index: adw_interp2.ftn =================================================================== --- adw_interp2.ftn 2006-06-28 14:35:14.000000000 -0400 +++ adw_interp2.ftn 2006-06-28 14:31:07.000000000 -0400 @@ -1,6 +1,6 @@ ***s/r adw_interp2 - Interpolation of one rhs * -#include "model_macros_f.h" +#include * subroutine adw_interp2 ( F_out, F_in, F_u, F_v, F_wind_L, % F_mono_L, DIST_DIM, Nk,i0,in,j0,jn ) @@ -106,15 +106,53 @@ * Interpolate * ************************************************************************ - call adw_tricub_lag3d ( Adw_wrkc, F_u, - % Adw_capx1, Adw_capy1, Adw_capz1, - % nijk, F_mono_L,i0,in,j0,jn,l_nk) + select case (Adw_interp_type(1:3)) + case ('tru') +#if defined (NEC) + call adw_tritrunc_lag3d_vec ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) +#else + call adw_tritrunc_lag3d ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) +#endif + case DEFAULT +#if defined (NEC) + call adw_tricub_lag3d_vec ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) +#else + call adw_tricub_lag3d ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) +#endif + end select * if (.not.G_lam) then if ( Adw_fro_a .gt. 0 ) then - call adw_tricub_lag3d ( Adw_wrka, F_u, - % Adw_capx2, Adw_capy2, Adw_capz2, - % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) + select case (Adw_interp_type(1:3)) + case ('tru') +#if defined (NEC) + call adw_tritrunc_lag3d_vec ( Adw_wrka, F_u, + % Adw_capx2, Adw_capy2, Adw_capz2, + % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) +#else + call adw_tritrunc_lag3d ( Adw_wrka, F_u, + % Adw_capx2, Adw_capy2, Adw_capz2, + % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) +#endif + case DEFAULT +#if defined (NEC) + call adw_tricub_lag3d_vec ( Adw_wrka, F_u, + % Adw_capx2, Adw_capy2, Adw_capz2, + % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) +#else + call adw_tricub_lag3d ( Adw_wrka, F_u, + % Adw_capx2, Adw_capy2, Adw_capz2, + % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) +#endif + end select endif * !$omp single Index: adw.cdk =================================================================== --- adw.cdk 2006-06-28 14:35:12.000000000 -0400 +++ adw.cdk 2006-06-28 14:31:08.000000000 -0400 @@ -8,6 +8,7 @@ * revision | * v3_20 - Gravel & Valin & Tanguay - Lagrange 3D | * and Optimized SETINT/TRILIN | +* v3_30 - McTaggart-Cowan R. - Add Adw_interp_type options | *______________________________________________________________________| * | * 3 different grids are refered to throughout the advection process: | @@ -39,7 +40,10 @@ * | lation is requested are inside own | * | advection source grid | * Adw_mono_L | switch: .true. : advection of tracers is monotonic | -* Adw_lag3d_L | switch: .true. : cubic Lagrange in 3D | +* Adw_interp_type interpolation scheme for back-trajectories ("cubic"):| +* | "cubic" Lagrangian(h)/cubic(v) interpolation | +* | "truncated" truncated Lagrangian 3D interpolation| +* | "lagrangian" tri-Lagrangian 3D interpolation | * Adw_nosetint_L switch: .true. : no setint for TRILIN (done inside) | * Adw_hor_L | switch: .true. : update horizontal positions | * | for Optimized SETINT/TRILIN | @@ -163,7 +167,7 @@ #endif * logical Adw_nkbz_L, Adw_exdg_L, Adw_ckbd_L, Adw_trunc_traj_L, - % Adw_mono_L,Adw_lag3d_L,Adw_nosetint_L,Adw_hor_L,Adw_ver_L + % Adw_mono_L,Adw_nosetint_L,Adw_hor_L,Adw_ver_L * integer Adw_halox, Adw_haloy, % Adw_nic, Adw_nit, Adw_njc, Adw_njt, @@ -175,9 +179,12 @@ real*8 Adw_ovdx_8, Adw_ovdy_8, Adw_ovdz_8, % Adw_x00_8, Adw_y00_8, Adw_z00_8, Adw_cfl_8(3) * + character(len=64) :: Adw_interp_type +* MARK_COMMON_BEG (adw_l) common / adw_l / Adw_nkbz_L, Adw_exdg_L, Adw_ckbd_L, Adw_trunc_traj_L, - % Adw_mono_L,Adw_lag3d_L,Adw_nosetint_L,Adw_hor_L,Adw_ver_L + % Adw_mono_L,Adw_interp_type,Adw_nosetint_L, + % Adw_hor_L,Adw_ver_L MARK_COMMON_END (adw_l) * MARK_COMMON_BEG (adw_ia) Index: nml.cdk =================================================================== --- nml.cdk 2006-06-28 14:35:15.000000000 -0400 +++ nml.cdk 2006-06-28 14:31:08.000000000 -0400 @@ -16,7 +16,7 @@ % Schm_pcsty_L, Schm_pheat_L, Schm_sfix_L, $ Adw_nkbz_L, Adw_exdg_L, $ Adw_ckbd_L, Adw_halox, Adw_haloy, Adw_trunc_traj_L, - $ Adw_mono_L, Adw_lag3d_L,Adw_nosetint_L, + $ Adw_mono_L, Adw_interp_type, Adw_nosetint_L, $ Glb_pil_n, Glb_pil_s, Glb_pil_w, Glb_pil_e, $ Cstv_dt_8, Cstv_uvdf_8, Cstv_phidf_8, Cori_cornl_L, $ Cstv_pitop_8, Cstv_pisrf_8, Cstv_tstr_8, @@ -49,7 +49,7 @@ % Schm_pcsty_L, Schm_pheat_L, Schm_sfix_L, $ Adw_nkbz_L, Adw_exdg_L, $ Adw_ckbd_L, Adw_halox, Adw_haloy, Adw_trunc_traj_L, - $ Adw_mono_L, Adw_lag3d_L,Adw_nosetint_L, + $ Adw_mono_L, Adw_interp_type, Adw_nosetint_L, $ Glb_pil_n, Glb_pil_s, Glb_pil_w, Glb_pil_e, $ Cstv_dt_8, Cstv_uvdf_8, Cstv_phidf_8, Cori_cornl_L, $ Cstv_pitop_8, Cstv_pisrf_8, Cstv_tstr_8, Index: adw_interp.ftn =================================================================== --- adw_interp.ftn 2006-06-28 14:35:14.000000000 -0400 +++ adw_interp.ftn 2006-06-28 14:31:08.000000000 -0400 @@ -22,7 +22,8 @@ * v2_31 - Tanguay M. - correction parameters adw_vder * v3_00 - Desgagne & Lee - Lam configuration * v3_10 - Corbeil & Desgagne & Lee - AIXport+Opti+OpenMP -* v3_20 - Gravel & Valin & Tanguay - Lagrange 3D +* v3_20 - Gravel & Valin & Tanguay - Lagrange 3D +* v3_20 - McTaggart-Cowan - Add semi-cubic Lagrangian interpolation * *language * fortran 77 @@ -103,30 +104,30 @@ endif ************************************************************************ * -* Compute second derivative for cubic spline in the vertical -* -************************************************************************ - if ( .not. Adw_lag3d_L ) call adw_vder ( F_v, F_u, Adw_nit, Adw_njt, l_nk ) -************************************************************************ -* * Interpolate * ************************************************************************ - if ( .not. Adw_lag3d_L ) then - call adw_tricub ( Adw_wrkc, F_u, F_v, - % Adw_n1, Adw_capx1, Adw_xgg1, - % Adw_xdd1, Adw_capy1, Adw_ygg1, - % Adw_ydd1, Adw_capz1,Adw_cz1, - % nijk, F_mono_L,i0,in,j0,jn,l_nk) - else - call adw_tricub_lag3d ( Adw_wrkc, F_u, - % Adw_capx1, Adw_capy1, Adw_capz1, - % nijk, F_mono_L,i0,in,j0,jn,l_nk) - endif + select case (Adw_interp_type(1:3)) + case ('lag') + call adw_tricub_lag3d ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) + case ('tru') + call adw_tritrunc_lag3d ( Adw_wrkc, F_u, + % Adw_capx1, Adw_capy1, Adw_capz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) + case DEFAULT + call adw_vder ( F_v, F_u, Adw_nit, Adw_njt, l_nk ) !2nd deriv for spline + call adw_tricub ( Adw_wrkc, F_u, F_v, + % Adw_n1, Adw_capx1, Adw_xgg1, + % Adw_xdd1, Adw_capy1, Adw_ygg1, + % Adw_ydd1, Adw_capz1,Adw_cz1, + % nijk, F_mono_L,i0,in,j0,jn,l_nk) + end select * if (.not.G_lam) then if ( Adw_fro_a .gt. 0 ) then - if ( .not. Adw_lag3d_L ) then + if ( Adw_interp_type(1:3) .eq. 'cub' ) then call adw_tricub ( Adw_wrka, F_u, F_v, % Adw_n2, Adw_capx2, Adw_xgg2, % Adw_xdd2, Adw_capy2, Adw_ygg2, Index: theonml.cdk =================================================================== --- theonml.cdk 2006-06-28 14:35:16.000000000 -0400 +++ theonml.cdk 2006-06-28 14:31:08.000000000 -0400 @@ -21,7 +21,7 @@ % Schm_itnlh, Schm_itraj, Schm_alpco_8, Schm_hdlast_L, % Schm_adcub_L, Schm_psadj_L, Schm_difqp_L, Schm_maxcfl_lam, % Schm_sfix_L, - % Adw_nkbz_L, Adw_exdg_L, Adw_mono_L, Adw_lag3d_L, + % Adw_nkbz_L, Adw_exdg_L, Adw_mono_L, Adw_interp_type, $ Adw_ckbd_L, Adw_halox, Adw_haloy, Adw_trunc_traj_L, $ Cstv_dt_8, Cstv_uvdf_8, Cstv_phidf_8, $ Cstv_pitop_8, Cstv_pisrf_8, Cstv_tstr_8, @@ -48,7 +48,7 @@ % Schm_itnlh, Schm_itraj, Schm_alpco_8, Schm_hdlast_L, % Schm_adcub_L, Schm_psadj_L, Schm_difqp_L, Schm_maxcfl_lam, % Schm_sfix_L, - % Adw_nkbz_L, Adw_exdg_L, Adw_mono_L, Adw_lag3d_L, + % Adw_nkbz_L, Adw_exdg_L, Adw_mono_L, Adw_interp_type, $ Adw_ckbd_L, Adw_halox, Adw_haloy, Adw_trunc_traj_L, $ Cstv_dt_8, Cstv_uvdf_8, Cstv_phidf_8, $ Cstv_pitop_8, Cstv_pisrf_8, Cstv_tstr_8, Index: adw_main_3_intlag.ftn =================================================================== --- adw_main_3_intlag.ftn 2006-06-28 14:43:38.000000000 -0400 +++ adw_main_3_intlag.ftn 2006-06-28 14:43:11.000000000 -0400 @@ -1,6 +1,6 @@ ***s/r adw_main_3_intlag - Interpolation of rhs * -#include "model_macros_f.h" +#include * subroutine adw_main_3_intlag ( F_u, F_v, F_w ) * @@ -60,6 +60,7 @@ #include "tr3d.cdk" #include "vt1.cdk" #include "vt0.cdk" +#include "adw_comp.cdk" * *modules integer vmmlod, vmmget, vmmuld, @@ -206,7 +207,7 @@ * ---------------------------- * Keep positions in CAP fields * ---------------------------- -!$omp parallel private(n) +!$omp parallel !$omp do do k=1,l_nk do j=j0,jn @@ -224,6 +225,8 @@ * Perform interpolation *********************************************************************** * + adw_comp_cub_L = .true. !force recomputation of position +* call adw_interp2 (ruw2, ruw1, F_u, F_v, % .true., .false., LDIST_DIM, l_nk,i0,in,j0,jn) *