Index: adw_tritrunc_lag3d.ftn =================================================================== --- adw_tritrunc_lag3d.ftn 2006-06-27 17:49:56.000000000 -0400 +++ adw_tritrunc_lag3d.ftn 2006-06-27 18:14:12.000000000 -0400 @@ -0,0 +1,277 @@ +***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 +* + 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 +* + prmax = 0.; prmin = 0. + a2 = 0.d0; a3 = 0.d0 + d2 = 0.d0; d3 = 0.d0 +* +* +!$omp parallel private(n,i,j,k, +!$omp& ii,jj,kk,capx,capy,capz,zcubic_L, +!$omp& o1, o2, o3, o4, a2, a3, b1, b2, b3, b4, +!$omp& c1, c2, c3, c4, d2, d3, p1, p2, p3, p4, +!$omp& rri,rrj,rrk,ra,rb,rc,rd,prmin,prmax) +!$omp do + do 100 k=1,kn + do 90 j=j0,jn + do 80 i=i0,in + 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) + 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)) + 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) + 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 +!$omp enddo +!$omp end parallel +* + return + end Index: adw_main_3_int.ftn =================================================================== --- adw_main_3_int.ftn 2006-06-27 17:37:16.000000000 -0400 +++ adw_main_3_int.ftn 2006-06-27 17:34:32.000000000 -0400 @@ -150,7 +150,7 @@ * endif * - if( .not. Adw_lag3d_L ) then + if( Adw_interp_type(1:3) .eq. 'cub' ) 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) @@ -176,7 +176,7 @@ * 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-27 17:35:17.000000000 -0400 +++ mtn_cfg.ftn 2006-06-27 17:35:42.000000000 -0400 @@ -108,7 +108,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: set_world_view1.ftn =================================================================== --- set_world_view1.ftn 2006-06-27 15:02:40.000000000 -0400 +++ set_world_view1.ftn 2006-06-27 17:36:32.000000000 -0400 @@ -53,6 +53,7 @@ * - 1d higher order diffusion operator * v3_20 - Zadra A. - Introduce V4dg_sgvc_dt0 * v3_20 - Buehner & Zadra - Introduce V4dg_chum_s +* v3_30 - McTaggart-Cowan R.- Add interpolator option Adw_interp * *object * To initialize model by reading namelist and return the @@ -195,7 +196,7 @@ Adw_exdg_L = .false. Adw_ckbd_L = .false. Adw_mono_L = .true. - Adw_lag3d_L= .false. + Adw_interp_type = 'cubic' Adw_nosetint_L = .false. Adw_trunc_traj_L = .false. Adw_halox = 3 @@ -541,7 +542,21 @@ * 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('cubic') + continue + 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) @@ -626,6 +641,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 "cubic", "lagrangian", "truncated".') * *------------------------------------------------------------------- * Index: Makefile =================================================================== --- Makefile 2006-06-26 16:05:18.000000000 -0400 +++ Makefile 2006-06-27 17:52:09.000000000 -0400 @@ -94,6 +94,7 @@ FTNDECKS= \ adw_pol0.ftn adw_pols.ftn adw_polw.ftn adw_polx.ftn \ adw_set.ftn adw_setint.ftn adw_trajex.ftn adw_trajsp.ftn \ adw_tricub.ftn adw_tricub_lag3d.ftn adw_trilin.ftn adw_trilin_turbo.ftn \ + adw_tritrunc_lag3d.ftn \ adw_vder.ftn bac.ftn bacp_2.ftn bcastc1.ftn \ bin2com.ftn blkcol.ftn blocstat.ftn bmf_perturb.ftn \ bubble.ftn bubble_cfg.ftn c_bloccplg.ftn c_cplg_iml.ftn \ @@ -180,6 +181,7 @@ FDECKS= \ adw_pol0.f adw_pols.f adw_polw.f adw_polx.f \ adw_set.f adw_setint.f adw_trajex.f adw_trajsp.f \ adw_tricub.f adw_tricub_lag3d.f adw_trilin.f adw_trilin_turbo.f \ + adw_tritrunc_lag3d.f \ adw_vder.f bac.f bacp_2.f bcastc1.f \ bin2com.f blkcol.f blocstat.f bmf_perturb.f \ bubble.f bubble_cfg.f c_bloccplg.f c_cplg_iml.f \ @@ -266,6 +268,7 @@ OBJECTS= \ adw_pol0.o adw_pols.o adw_polw.o adw_polx.o \ adw_set.o adw_setint.o adw_trajex.o adw_trajsp.o \ adw_tricub.o adw_tricub_lag3d.o adw_trilin.o adw_trilin_turbo.o \ + adw_tritrunc_lag3d.o \ adw_vder.o bac.o bacp_2.o bcastc1.o \ bin2com.o blkcol.o blocstat.o bmf_perturb.o \ bubble.o bubble_cfg.o c_bloccplg.o c_cplg_iml.o \ @@ -426,6 +429,7 @@ adw_trajex.f: adw_trajex.ftn adw.cdk d adw_trajsp.f: adw_trajsp.ftn adw.cdk dcst.cdk glb_ld.cdk adw_tricub.f: adw_tricub.ftn adw.cdk glb_ld.cdk impnone.cdk adw_tricub_lag3d.f: adw_tricub_lag3d.ftn adw.cdk glb_ld.cdk +adw_tritrunc_lag3d.f: adw_tritrunc_lag3d.ftn adw.cdk glb_ld.cdk adw_trilin.f: adw_trilin.ftn adw.cdk glb_ld.cdk adw_trilin_turbo.f: adw_trilin_turbo.ftn adw.cdk glb_ld.cdk adw_vder.f: adw_vder.ftn adw.cdk geomg.cdk glb_ld.cdk \ @@ -1158,6 +1162,7 @@ adw_trajex.o: adw_trajex.ftn adw.cdk d adw_trajsp.o: adw_trajsp.ftn adw.cdk dcst.cdk glb_ld.cdk adw_tricub.o: adw_tricub.ftn adw.cdk glb_ld.cdk impnone.cdk adw_tricub_lag3d.o: adw_tricub_lag3d.ftn adw.cdk glb_ld.cdk +adw_tritrunc_lag3d.o: adw_tritrunc_lag3d.ftn adw.cdk glb_ld.cdk adw_trilin.o: adw_trilin.ftn adw.cdk glb_ld.cdk adw_trilin_turbo.o: adw_trilin_turbo.ftn adw.cdk glb_ld.cdk adw_vder.o: adw_vder.ftn adw.cdk geomg.cdk glb_ld.cdk \ @@ -1923,7 +1928,7 @@ compil: fn=$${i%.cdk90}; \ $(MAKE) "DEFINE=$(DEFINE)" $$fn.o; \ done \ - fi ; \ + fi ; \ if [ "*.ftn90" != "`echo *.ftn90`" ] ; \ then \ for i in *.ftn90; \ @@ -1976,7 +1981,7 @@ compexp: sortirtout $(MAKE) "DEFINE=$(DEFINE)" "FFLAGS=$(FFLAGS)" $$fn.o; \ $(MAKE) -f make_cdk $$fn.acdk90; \ done \ - fi ; \ + fi ; \ if [ "*.ftn90" != "`echo *.ftn90`" ] ; \ then \ for i in *.ftn90; \ Index: adw.cdk =================================================================== --- adw.cdk 2006-06-27 15:02:40.000000000 -0400 +++ adw.cdk 2006-06-27 17:31:46.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-27 15:02:40.000000000 -0400 +++ nml.cdk 2006-06-27 17:16:57.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, @@ -48,7 +48,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-27 15:02:40.000000000 -0400 +++ adw_interp.ftn 2006-06-28 11:50:44.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,40 +104,45 @@ 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 - call adw_tricub ( Adw_wrka, F_u, F_v, - % Adw_n2, Adw_capx2, Adw_xgg2, - % Adw_xdd2, Adw_capy2, Adw_ygg2, - % Adw_ydd2, Adw_capz2, Adw_cz2, - % 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 + select case (Adw_interp_type(1:3)) + case ('lag') + 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) + case ('tru') + 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) + case DEFAULT + call adw_tricub ( Adw_wrka, F_u, F_v, + % Adw_n2, Adw_capx2, Adw_xgg2, + % Adw_xdd2, Adw_capy2, Adw_ygg2, + % Adw_ydd2, Adw_capz2, Adw_cz2, + % Adw_fro_a, F_mono_L,1,Adw_fro_a,1,1,1) + end select endif * call adw_exch_2 ( Adw_wrkb, dummy, dummy, Index: theonml.cdk =================================================================== --- theonml.cdk 2006-06-27 15:02:40.000000000 -0400 +++ theonml.cdk 2006-06-27 17:15:57.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,