Speed-up of adjoint transport

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
+

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))

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
 

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
+

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
 

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

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

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

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

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
 

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