diff --git a/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch b/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch index a2dd2540..ad1aa671 100644 --- a/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch +++ b/science/adjoint/patches/kernel/adj_dg_matrix_vector_kernel_mod.patch @@ -27,3 +27,21 @@ contains subroutine adj_dg_matrix_vector_code_r_single(cell, nlayers, lhs, x, ncell_3d, matrix, ndf1, undf1, map1, ndf2, undf2, map2) +@@ -56,7 +56,7 @@ + do df1 = ndf1, 1, -1 + i1 = map1(df1) + do idx = i1 + nl, i1, -1 +- lhs(idx) = 0.0 ++ lhs(idx) = 0.0_r_single + enddo + enddo + +@@ -97,7 +97,7 @@ + do df1 = ndf1, 1, -1 + i1 = map1(df1) + do idx = i1 + nl, i1, -1 +- lhs(idx) = 0.0 ++ lhs(idx) = 0.0_r_double + enddo + enddo + diff --git a/science/adjoint/patches/kernel/adj_horizontal_mass_flux_kernel_mod.patch b/science/adjoint/patches/kernel/adj_horizontal_mass_flux_kernel_mod.patch new file mode 100644 index 00000000..41396902 --- /dev/null +++ b/science/adjoint/patches/kernel/adj_horizontal_mass_flux_kernel_mod.patch @@ -0,0 +1,11 @@ +@@ -38,9 +38,7 @@ + real(kind=r_tran) :: direction + real(kind=r_tran), dimension(nfaces) :: v_dot_n + +- v_dot_n(:) = 1.0_r_tran +- v_dot_n(1) = -1.00000000000000 +- v_dot_n(nfaces) = -1.00000000000000 ++ v_dot_n = (/ -1.0_r_tran, 1.0_r_tran, 1.0_r_tran, -1.0_r_tran /) + do df = nfaces, 1, -1 + do k = nlayers - 1, 0, -1 + direction = v_dot_n(df) * wind(k + map_w2(df)) diff --git a/science/adjoint/patches/kernel/adj_sci_combine_multidata_field_kernel_mod.patch b/science/adjoint/patches/kernel/adj_sci_combine_multidata_field_kernel_mod.patch index 8a7e3191..766fa689 100644 --- a/science/adjoint/patches/kernel/adj_sci_combine_multidata_field_kernel_mod.patch +++ b/science/adjoint/patches/kernel/adj_sci_combine_multidata_field_kernel_mod.patch @@ -34,7 +34,78 @@ contains subroutine adj_combine_multidata_field_code_r_single(nlayers, field_out, n, field1_in, n1, field2_in, n2, ndata_first, ndf, & -@@ -133,4 +133,4 @@ +@@ -49,15 +49,15 @@ + + if (ndata_first) then + do k = nlayers - 1, 0, -1 +- ij = k * n + map(1) + ij = k * n + n1 + map(1) + do df = n2 - 1, 0, -1 + field2_in(map_2(1) + k * n2 + df) = field2_in(map_2(1) + k * n2 + df) + field_out(ij + df) +- field_out(ij + df) = 0.0 ++ field_out(ij + df) = 0.0_r_single + enddo ++ ij = k * n + map(1) + do df = n1 - 1, 0, -1 + field1_in(map_1(1) + k * n1 + df) = field1_in(map_1(1) + k * n1 + df) + field_out(ij + df) +- field_out(ij + df) = 0.0 ++ field_out(ij + df) = 0.0_r_single + enddo + enddo + else +@@ -65,14 +65,14 @@ + ij = df * nlayers + n1 * nlayers - nlayers + map(1) + do k = nlayers - 1, 0, -1 + field2_in(map_2(1) + (df - 1) * nlayers + k) = field2_in(map_2(1) + (df - 1) * nlayers + k) + field_out(ij + k) +- field_out(ij + k) = 0.0 ++ field_out(ij + k) = 0.0_r_single + enddo + enddo + do df = n1, 1, -1 + ij = df * nlayers - nlayers + map(1) + do k = nlayers - 1, 0, -1 + field1_in(map_1(1) + (df - 1) * nlayers + k) = field1_in(map_1(1) + (df - 1) * nlayers + k) + field_out(ij + k) +- field_out(ij + k) = 0.0 ++ field_out(ij + k) = 0.0_r_single + enddo + enddo + end if +@@ -103,15 +103,15 @@ + + if (ndata_first) then + do k = nlayers - 1, 0, -1 +- ij = k * n + map(1) + ij = k * n + n1 + map(1) + do df = n2 - 1, 0, -1 + field2_in(map_2(1) + k * n2 + df) = field2_in(map_2(1) + k * n2 + df) + field_out(ij + df) +- field_out(ij + df) = 0.0 ++ field_out(ij + df) = 0.0_r_double + enddo ++ ij = k * n + map(1) + do df = n1 - 1, 0, -1 + field1_in(map_1(1) + k * n1 + df) = field1_in(map_1(1) + k * n1 + df) + field_out(ij + df) +- field_out(ij + df) = 0.0 ++ field_out(ij + df) = 0.0_r_double + enddo + enddo + else +@@ -119,18 +119,18 @@ + ij = df * nlayers + n1 * nlayers - nlayers + map(1) + do k = nlayers - 1, 0, -1 + field2_in(map_2(1) + (df - 1) * nlayers + k) = field2_in(map_2(1) + (df - 1) * nlayers + k) + field_out(ij + k) +- field_out(ij + k) = 0.0 ++ field_out(ij + k) = 0.0_r_double + enddo + enddo + do df = n1, 1, -1 + ij = df * nlayers - nlayers + map(1) + do k = nlayers - 1, 0, -1 + field1_in(map_1(1) + (df - 1) * nlayers + k) = field1_in(map_1(1) + (df - 1) * nlayers + k) + field_out(ij + k) +- field_out(ij + k) = 0.0 ++ field_out(ij + k) = 0.0_r_double + enddo + enddo + end if end subroutine adj_combine_multidata_field_code_r_double diff --git a/science/adjoint/patches/kernel/adj_w3v_advective_update_kernel_mod.patch b/science/adjoint/patches/kernel/adj_w3v_advective_update_kernel_mod.patch new file mode 100644 index 00000000..350fd603 --- /dev/null +++ b/science/adjoint/patches/kernel/adj_w3v_advective_update_kernel_mod.patch @@ -0,0 +1,42 @@ +@@ -45,12 +45,7 @@ + integer(kind=i_def) :: offset + real(kind=r_tran) :: w + real(kind=r_tran) :: dtdz +- real(kind=r_tran) :: t_u +- real(kind=r_tran) :: t_d + +- dtdz = 0.0_r_tran +- t_u = 0.0_r_tran +- t_d = 0.0_r_tran + if (ndf_w2 == 2) then + df = 1 + offset = 0 +@@ -61,23 +56,16 @@ + do k = nlayers - 1, 0, -1 + w = 0.5 * wind(k + map_w2(df)) + 0.5 * wind(k + map_w2(df) + 1) + ik = cell * nlayers + k - nlayers + 1 +- dtdz = dtdz + advective_increment(map_w3(1) + k) * w * REAL(m3_inv(ik,1,1), r_tran) +- t_d = t_d + (-dtdz) +- t_u = t_u + dtdz +- dtdz = 0.0 ++ dtdz = advective_increment(map_w3(1) + k) * w * REAL(m3_inv(ik,1,1), r_tran) + if (w <= 0.0_r_tran .AND. k < nlayers - 1) then +- tracer(map_md(1) + offset + k + 1) = tracer(map_md(1) + offset + k + 1) + t_u +- t_u = 0.0 ++ tracer(map_md(1) + offset + k + 1) = tracer(map_md(1) + offset + k + 1) + dtdz + else +- tracer(map_md(1) + offset + nlayers + k) = tracer(map_md(1) + offset + nlayers + k) + t_u +- t_u = 0.0 ++ tracer(map_md(1) + offset + nlayers + k) = tracer(map_md(1) + offset + nlayers + k) + dtdz + end if + if (w > 0.0_r_tran .AND. k > 0) then +- tracer(map_md(1) + offset + nlayers + k - 1) = tracer(map_md(1) + offset + nlayers + k - 1) + t_d +- t_d = 0.0 ++ tracer(map_md(1) + offset + nlayers + k - 1) = tracer(map_md(1) + offset + nlayers + k - 1) - dtdz + else +- tracer(map_md(1) + offset + k) = tracer(map_md(1) + offset + k) + t_d +- t_d = 0.0 ++ tracer(map_md(1) + offset + k) = tracer(map_md(1) + offset + k) - dtdz + end if + enddo + diff --git a/science/adjoint/patches/kernel/atl_horizontal_mass_flux_kernel_mod.patch b/science/adjoint/patches/kernel/atl_horizontal_mass_flux_kernel_mod.patch index 9e6a07c0..aa04ceb4 100644 --- a/science/adjoint/patches/kernel/atl_horizontal_mass_flux_kernel_mod.patch +++ b/science/adjoint/patches/kernel/atl_horizontal_mass_flux_kernel_mod.patch @@ -8,7 +8,7 @@ - type, public, extends(kernel_type) :: adj_horizontal_mass_flux_kernel_type + type, public, extends(kernel_type) :: atl_horizontal_mass_flux_kernel_type type(ARG_TYPE) :: META_ARGS(4) = (/ & - arg_type(gh_field, gh_real, gh_inc, any_w2), & + arg_type(gh_field, gh_real, gh_read, any_w2), & arg_type(gh_field, gh_real, gh_read, any_w2), & @@ -11,15 +11,15 @@ arg_type(gh_field, gh_real, gh_read, any_discontinuous_space_1)/) @@ -30,7 +30,18 @@ &undf_md, map_md) integer(kind=i_def), parameter :: nfaces = 4 integer(kind=i_def), intent(in) :: nlayers -@@ -53,6 +53,6 @@ +@@ -39,9 +39,7 @@ + real(kind=r_def) :: direction + real(kind=r_def), dimension(nfaces) :: v_dot_n + +- v_dot_n(:) = 1.0_r_def +- v_dot_n(1) = -1.00000000000000 +- v_dot_n(nfaces) = -1.00000000000000 ++ v_dot_n = (/ -1.0_r_def, 1.0_r_def, 1.0_r_def, -1.0_r_def /) + do df = nfaces, 1, -1 + do k = nlayers - 1, 0, -1 + direction = ls_wind(k + map_w2(df)) * v_dot_n(df) +@@ -52,6 +50,6 @@ enddo enddo diff --git a/science/adjoint/patches/kernel/atl_w3v_advective_update_kernel_mod.patch b/science/adjoint/patches/kernel/atl_w3v_advective_update_kernel_mod.patch index 5f4d873e..7cd2087c 100644 --- a/science/adjoint/patches/kernel/atl_w3v_advective_update_kernel_mod.patch +++ b/science/adjoint/patches/kernel/atl_w3v_advective_update_kernel_mod.patch @@ -32,8 +32,23 @@ &undf_w3, map_w3, ndf_md, undf_md, map_md, ndf_w2, undf_w2, map_w2) integer(kind=i_def), intent(in) :: nlayers integer(kind=i_def), intent(in) :: cell -@@ -79,6 +79,6 @@ - w = 0.0 +@@ -51,7 +51,6 @@ + real(kind=r_def) :: t_u + real(kind=r_def) :: t_d + +- w = 0.0_r_def + if (ndf_w2 == 2) then + df = 1 + offset = 0 +@@ -73,12 +72,11 @@ + end if + dtdz = -t_d + t_u + ik = cell * nlayers + k - nlayers + 1 +- w = w + dtdz * advective_increment(map_w3(1) + k) * m3_inv(ik,1,1) ++ w = dtdz * advective_increment(map_w3(1) + k) * m3_inv(ik,1,1) + wind(k + map_w2(df)) = wind(k + map_w2(df)) + 0.5 * w + wind(k + map_w2(df) + 1) = wind(k + map_w2(df) + 1) + 0.5 * w +- w = 0.0 enddo - end subroutine adj_w3v_advective_update_code diff --git a/science/adjoint/source/algorithm/transport/mol/atl_mol_advective_alg_mod.x90 b/science/adjoint/source/algorithm/transport/mol/atl_mol_advective_alg_mod.x90 index cec84c38..ae4c47cc 100644 --- a/science/adjoint/source/algorithm/transport/mol/atl_mol_advective_alg_mod.x90 +++ b/science/adjoint/source/algorithm/transport/mol/atl_mol_advective_alg_mod.x90 @@ -205,8 +205,7 @@ contains ! Update field: f = f^n - dt_substep*rhs call invoke( inc_X_minus_bY( rhs, dt_mol_substep, field_np1 ), & - inc_X_plus_Y( field_n, field_np1 ), & - setval_c( field_np1, 0.0_r_def ) ) + inc_X_plus_Y( field_n, field_np1 ) ) ! Compute update: rhs = u.grad(rhs_field) call atl_advective_and_flux_alg( & @@ -222,8 +221,8 @@ contains ! Compute the field for this stage: ! rhs_field = sum(s=1,stage): a(stage,s)*field^(s) - call invoke( setval_c( rhs_field, 0.0_r_def ), & - inc_X_plus_Y( field_np1, rk_field(stage) ), & + call invoke( setval_c( rhs_field, 0.0_r_def ), & + setval_x( field_np1, rk_field(stage) ), & setval_c( rk_field(stage), 0.0_r_def ) ) end do ! stage diff --git a/science/adjoint/source/algorithm/transport/mol/atl_mol_conservative_alg_mod.x90 b/science/adjoint/source/algorithm/transport/mol/atl_mol_conservative_alg_mod.x90 index b68fd309..8f2b150b 100644 --- a/science/adjoint/source/algorithm/transport/mol/atl_mol_conservative_alg_mod.x90 +++ b/science/adjoint/source/algorithm/transport/mol/atl_mol_conservative_alg_mod.x90 @@ -381,8 +381,7 @@ module atl_mol_conservative_alg_mod end if if ( do_advective ) then call invoke( inc_X_minus_bY( adv_inc, dt_mol_substep, field_np1 ), & - inc_X_plus_Y( field_n, field_np1 ), & - setval_c( field_np1, 0.0_r_def ) ) + inc_X_plus_Y( field_n, field_np1 ) ) end if call atl_advective_and_flux_alg( & @@ -399,8 +398,8 @@ module atl_mol_conservative_alg_mod ! Compute the field for this stage: ! rhs_field = sum(s=1,stage): a(stage,s)*field^(s) - call invoke( setval_c( rhs_field, 0.0_r_def ), & - inc_X_plus_Y( field_np1, rk_field(stage) ), & + call invoke( setval_c( rhs_field, 0.0_r_def ), & + setval_x( field_np1, rk_field(stage) ), & setval_c( rk_field(stage), 0.0_r_def ) ) end do stage_loop diff --git a/science/adjoint/source/kernel/transport/mol/atl_poly1d_vert_w3_reconstruction_kernel_mod.F90 b/science/adjoint/source/kernel/transport/mol/atl_poly1d_vert_w3_reconstruction_kernel_mod.F90 index d2889501..e330c4ba 100644 --- a/science/adjoint/source/kernel/transport/mol/atl_poly1d_vert_w3_reconstruction_kernel_mod.F90 +++ b/science/adjoint/source/kernel/transport/mol/atl_poly1d_vert_w3_reconstruction_kernel_mod.F90 @@ -129,7 +129,6 @@ subroutine atl_poly1d_vert_w3_reconstruction_code( nlayers, & real(kind=r_def) :: new_tracer real(kind=r_def) :: ls_new_tracer - new_tracer = 0.0_r_def vertical_order = MIN(global_order, nlayers - 1) ij = map_w3(1) if (logspace) then @@ -140,26 +139,24 @@ subroutine atl_poly1d_vert_w3_reconstruction_code( nlayers, & ik = f * global_order + f + k * ndata + p + map_c(1) - 1 ls_new_tracer = ls_new_tracer * MAX(eps, ABS(ls_tracer(ij + stencil(p,k,f)))) ** coeff(ik) end do - new_tracer = new_tracer + ls_new_tracer * reconstruction(map_md(1) + f * nlayers + k) + new_tracer = ls_new_tracer * reconstruction(map_md(1) + f * nlayers + k) reconstruction(map_md(1) + f * nlayers + k) = 0.0_r_def do p = vertical_order + 1, 1, -1 ik = f * global_order + f + k * ndata + p + map_c(1) - 1 tracer(ij + stencil(p,k,f)) = tracer(ij + stencil(p,k,f)) + coeff(ik) * new_tracer / SIGN(MAX(eps, ls_tracer(ij + & &stencil(p,k,f))), ls_tracer(ij + stencil(p,k,f))) enddo - new_tracer = 0.0_r_def enddo enddo else do f = 1, 0, -1 do k = nlayers - 1, 0, -1 - new_tracer = new_tracer + reconstruction(map_md(1) + f * nlayers + k) + new_tracer = reconstruction(map_md(1) + f * nlayers + k) reconstruction(map_md(1) + f * nlayers + k) = 0.0_r_def do p = vertical_order + 1, 1, -1 ik = f * global_order + f + k * ndata + p + map_c(1) - 1 tracer(ij + stencil(p,k,f)) = tracer(ij + stencil(p,k,f)) + coeff(ik) * new_tracer enddo - new_tracer = 0.0_r_def enddo enddo end if diff --git a/science/adjoint/source/kernel/transport/mol/atl_poly_adv_update_kernel_mod.F90 b/science/adjoint/source/kernel/transport/mol/atl_poly_adv_update_kernel_mod.F90 index 972aa52d..e1a89201 100644 --- a/science/adjoint/source/kernel/transport/mol/atl_poly_adv_update_kernel_mod.F90 +++ b/science/adjoint/source/kernel/transport/mol/atl_poly_adv_update_kernel_mod.F90 @@ -152,10 +152,7 @@ subroutine atl_poly_adv_update_code( nlayers, & real(kind=r_tran) :: ls_dtdy - uv = 0.0_r_tran - v_dot_n(:) = 1.0_r_tran - v_dot_n(1) = -1.0_r_tran - v_dot_n(nfaces) = -1.0_r_tran + v_dot_n = (/ -1.0_r_tran, 1.0_r_tran, 1.0_r_tran, -1.0_r_tran /) opposite(:) = -1 missing_neighbour(:) = .false. @@ -186,9 +183,7 @@ subroutine atl_poly_adv_update_code( nlayers, & uv_dir(1,k) = 0.25_r_tran * wind_dir(k + map_w2(1) - 1) + 0.25_r_tran * wind_dir(k + map_w2(3) - 1) uv_dir(2,k) = 0.25_r_tran * wind_dir(k + map_w2(2) - 1) + 0.25_r_tran * wind_dir(k + map_w2(4) - 1) - direction_dofs(:) = 1 - direction_dofs(2) = 2 - direction_dofs(4) = 2 + direction_dofs(:) = (/ 1, 2, 1, 2 /) do df = 1, nfaces, 1 do k = 0, nlayers, 1 @@ -206,39 +201,33 @@ subroutine atl_poly_adv_update_code( nlayers, & do k = nlayers, 0, -1 ls_dtdx = ls_tracer(e,k) - ls_tracer(w,k) ls_dtdy = ls_tracer(n,k) - ls_tracer(s,k) - uv(1,k) = uv(1,k) + ls_dtdx * advective(map_wt(1) + k) - uv(2,k) = uv(2,k) - ls_dtdy * advective(map_wt(1) + k) + uv(1,k) = ls_dtdx * advective(map_wt(1) + k) + uv(2,k) = -ls_dtdy * advective(map_wt(1) + k) advective(map_wt(1) + k) = 0.0_r_tran enddo k = nlayers wind(k + map_w2(2) - 1) = wind(k + map_w2(2) - 1) + 0.25_r_tran * uv(2,k) wind(k + map_w2(4) - 1) = wind(k + map_w2(4) - 1) + 0.25_r_tran * uv(2,k) - uv(2,k) = 0.0_r_tran wind(k + map_w2(1) - 1) = wind(k + map_w2(1) - 1) + 0.25_r_tran * uv(1,k) wind(k + map_w2(3) - 1) = wind(k + map_w2(3) - 1) + 0.25_r_tran * uv(1,k) - uv(1,k) = 0.0_r_tran do k = nlayers - 1, 1, -1 wind(k + map_w2(2)) = wind(k + map_w2(2)) + 0.25_r_tran * uv(2,k) wind(k + map_w2(4)) = wind(k + map_w2(4)) + 0.25_r_tran * uv(2,k) wind(k + map_w2(2) - 1) = wind(k + map_w2(2) - 1) + 0.25_r_tran * uv(2,k) wind(k + map_w2(4) - 1) = wind(k + map_w2(4) - 1) + 0.25_r_tran * uv(2,k) - uv(2,k) = 0.0_r_tran wind(k + map_w2(1)) = wind(k + map_w2(1)) + 0.25_r_tran * uv(1,k) wind(k + map_w2(3)) = wind(k + map_w2(3)) + 0.25_r_tran * uv(1,k) wind(k + map_w2(1) - 1) = wind(k + map_w2(1) - 1) + 0.25_r_tran * uv(1,k) wind(k + map_w2(3) - 1) = wind(k + map_w2(3) - 1) + 0.25_r_tran * uv(1,k) - uv(1,k) = 0.0_r_tran enddo k = 0 wind(map_w2(2)) = wind(map_w2(2)) + 0.25_r_tran * uv(2,k) wind(map_w2(4)) = wind(map_w2(4)) + 0.25_r_tran * uv(2,k) - uv(2,k) = 0.0_r_tran wind(map_w2(1)) = wind(map_w2(1)) + 0.25_r_tran * uv(1,k) wind(map_w2(3)) = wind(map_w2(3)) + 0.25_r_tran * uv(1,k) - uv(1,k) = 0.0_r_tran end subroutine atl_poly_adv_update_code diff --git a/science/adjoint/source/kernel/transport/mol/atl_w3h_advective_update_kernel_mod.F90 b/science/adjoint/source/kernel/transport/mol/atl_w3h_advective_update_kernel_mod.F90 index fca1ec80..021af5e9 100644 --- a/science/adjoint/source/kernel/transport/mol/atl_w3h_advective_update_kernel_mod.F90 +++ b/science/adjoint/source/kernel/transport/mol/atl_w3h_advective_update_kernel_mod.F90 @@ -165,8 +165,6 @@ subroutine atl_w3h_advective_update_code( cell, & ! so if u.n > 0 then we set the field to be the value on this edge from this cell ! and if u.n < 0 then we set the field to be the value on this edge from a ! neighbouring cell - u = 0.0_r_def - v = 0.0_r_def do k = nlayers - 1, 0, -1 ! u * dt/dx @@ -218,16 +216,14 @@ subroutine atl_w3h_advective_update_code( cell, & dtdy = t_N - t_S ik = 1 + k + (cell-1)*nlayers - u = u + m3_inv(ik,1,1) * dtdx * advective_increment(map_w3(1)+k) - v = v + m3_inv(ik,1,1) * dtdy * advective_increment(map_w3(1)+k) + u = m3_inv(ik,1,1) * dtdx * advective_increment(map_w3(1)+k) + v = m3_inv(ik,1,1) * dtdy * advective_increment(map_w3(1)+k) wind(map_w2(2) + k) = wind(map_w2(2) + k) - 0.5_r_def * v wind(map_w2(4) + k) = wind(map_w2(4) + k) - 0.5_r_def * v - v = 0.0_r_def wind(map_w2(1) + k) = wind(map_w2(1) + k) + 0.5_r_def * u wind(map_w2(3) + k) = wind(map_w2(3) + k) + 0.5_r_def * u - u = 0.0_r_def end do end subroutine atl_w3h_advective_update_code