tortor_newmap Subroutine

    subroutine tortor_newmap(bds, d1, d2, ang1, ang2, v)
        !! Store in module memory the data describing a new torsion-torsion 
        !! map
        use mod_memory, only: mallocate, mfree
        use mod_utils, only: cyclic_spline

        implicit none

        type(ommp_bonded_type), intent(inout) :: bds
        ! Bonded potential data structure
        integer(ip), intent(in) :: d1, d2
        !! Dimensions of the new map to be saved
        real(rp), intent(in) :: ang1(:)
        !! Value of torsion1 for the new map 
        real(rp), intent(in) :: ang2(:)
        !! Value of torsion2 for the new map 
        real(rp), intent(in) :: v(:)
        !! Value of potential for the new map 

        integer :: i, j, ii
        real(rp), allocatable, dimension(:) :: a, b, c, d, dx, dy, dxy, tmpx, tmpy
        
        real(rp), allocatable :: rtmp(:)
        integer(ip), allocatable :: itmp(:,:)
        integer(ip) :: n_data, n_map

        if(allocated(bds%ttmap_ang1)) then
            ! Reallocate the arrays to make space for the new data
            n_data = size(bds%ttmap_ang1)
            call mallocate('torstors_newmap [rtmp]', n_data, rtmp)
            
            rtmp = bds%ttmap_ang1
            call mfree('torstors_newmap [ttmap_ang1]', bds%ttmap_ang1)
            call mallocate('torstors_newmap [ttmap_ang1]', &
                           n_data+d1*d2,  bds%ttmap_ang1)
            bds%ttmap_ang1(:n_data) = rtmp
            
            rtmp = bds%ttmap_ang2
            call mfree('torstors_newmap [ttmap_ang2]', bds%ttmap_ang2)
            call mallocate('torstors_newmap [ttmap_ang2]', &
                           n_data+d1*d2,  bds%ttmap_ang2)
            bds%ttmap_ang2(:n_data) = rtmp
            
            
            call mfree('torstors_newmap [rtmp]', rtmp)
            n_data = size(bds%ttmap_v)
            call mallocate('torstors_newmap [rtmp]', n_data, rtmp)
            
            rtmp = bds%ttmap_v
            call mfree('torstors_newmap [ttmap_v]', bds%ttmap_v)
            call mallocate('torstors_newmap [ttmap_v]', &
                           n_data+d1*d2,  bds%ttmap_v)
            bds%ttmap_v(:n_data) = rtmp
            
            rtmp = bds%ttmap_vx
            call mfree('torstors_newmap [ttmap_vx]', bds%ttmap_vx)
            call mallocate('torstors_newmap [ttmap_vx]', &
                           n_data+d1*d2,  bds%ttmap_vx)
            bds%ttmap_vx(:n_data) = rtmp

            rtmp = bds%ttmap_vy
            call mfree('torstors_newmap [ttmap_vy]', bds%ttmap_vy)
            call mallocate('torstors_newmap [ttmap_vy]', &
                           n_data+d1*d2,  bds%ttmap_vy)
            bds%ttmap_vy(:n_data) = rtmp

            rtmp = bds%ttmap_vxy
            call mfree('torstors_newmap [ttmap_vxy]', bds%ttmap_vxy)
            call mallocate('torstors_newmap [ttmap_vxy]', &
                           n_data+d1*d2,  bds%ttmap_vxy)
            bds%ttmap_vxy(:n_data) = rtmp
            call mfree('torstors_newmap [rtmp]', rtmp)

            n_map = size(bds%ttmap_shape, 2)
            call mallocate('torstors_newmap [itmp]', 2, n_map, itmp)
            itmp = bds%ttmap_shape
            call mfree('torstors_newmap [ttmap_shape]', bds%ttmap_shape)
            call mallocate('torstors_newmap [ttmap_shape]', &
                           2, n_map+1, bds%ttmap_shape)
            bds%ttmap_shape(:,:n_map) = itmp

            call mfree('torstors_newmap [itmp]', itmp)
        else 
            ! First allocation, n_data and n_map are just set for consistency
            n_data = 0
            n_map = 0
            call mallocate('torstors_newmap [ttmap_ang1]', d1*d2,  bds%ttmap_ang1)
            call mallocate('torstors_newmap [ttmap_ang2]', d1*d2,  bds%ttmap_ang2)
            call mallocate('torstors_newmap [ttmap_v]', d1*d2,  bds%ttmap_v)
            call mallocate('torstors_newmap [ttmap_vx]', d1*d2,  bds%ttmap_vx)
            call mallocate('torstors_newmap [ttmap_vy]', d1*d2,  bds%ttmap_vy)
            call mallocate('torstors_newmap [ttmap_vxy]', d1*d2,  bds%ttmap_vxy)
            call mallocate('torstors_newmap [ttmap_shape]', 2, 1, bds%ttmap_shape)
        end if

        call mallocate('tortor_newmap [a]', max(d1,d2), a)
        call mallocate('tortor_newmap [b]', max(d1,d2), b)
        call mallocate('tortor_newmap [c]', max(d1,d2), c)
        call mallocate('tortor_newmap [d]', max(d1,d2), d)
        call mallocate('tortor_newmap [dx]', d1*d2, dx)
        call mallocate('tortor_newmap [dy]', d1*d2, dy)
        call mallocate('tortor_newmap [dxy]', d1*d2, dxy)

        ! This part of the code computes df/dx, df/dy and d^2f/dxdy on the grid.
        ! Since we are basically interpolating on a sphere, we extract the 
        ! coordinate on a meridian, we interpolate it with a cubic spline, and
        ! finally we compute the derivative of this curve at the grid intersection
        ! The same is done in the second direction.
        ! To compute the mixed derivative we apply the same procedure but using
        ! the derivative data (basically we apply the procedure used to compute
        ! df/dx but using  df/dy data instead of actual f values.
        do i=1, d2
            call cyclic_spline(d1, ang1((i-1)*d1+1:i*d1), v((i-1)*d1+1:i*d1), &
                               a(1:d1), b(1:d1), c(1:d1), d(1:d1))
            dx((i-1)*d1+1:i*d1) = b(1:d1)
        end do
        
        ! df/dy since in this direction data are not contiguous, wa allocate 
        ! temporary arrays
        call mallocate('tortor_newmap [tmpx]', d2, tmpx)
        call mallocate('tortor_newmap [tmpy]', d2, tmpy)
        do i=1, d1
            ii = 1
            do j=i, (d2-1)*d1+i, d2
                tmpx(ii) = ang2(j)
                tmpy(ii) = v(j)
                ii = ii + 1
            end do
            call cyclic_spline(d2, tmpx, tmpy, &
                               a(1:d2), b(1:d2), c(1:d2), d(1:d2))
            
            ii = 1
            do j=i, (d2-1)*d1+i, d2
                dy(j) = b(ii)
                ii = ii + 1
            end do
        end do
        
        ! d^2f/dxdy in this case we use df/dx procedure to exploit data contiguity.
        do i=1, d2
            call cyclic_spline(d1, ang1((i-1)*d1+1:i*d1), dy((i-1)*d1+1:i*d1), &
                               a(1:d1), b(1:d1), c(1:d1), d(1:d1))
            dxy((i-1)*d1+1:i*d1) = b(1:d1)
        end do
        call mfree('tortor_newmap [tmpx]', tmpx)
        call mfree('tortor_newmap [tmpy]', tmpy)

        bds%ttmap_ang1(n_data+1:) = ang1
        bds%ttmap_ang2(n_data+1:) = ang2
        bds%ttmap_shape(1,n_map+1) = d1
        bds%ttmap_shape(2,n_map+1) = d2
        bds%ttmap_v(n_data+1:) = v
        bds%ttmap_vx(n_data+1:) = dx
        bds%ttmap_vy(n_data+1:) = dy
        bds%ttmap_vxy(n_data+1:) = dxy
        
        call mfree('tortor_newmap [a]', a)
        call mfree('tortor_newmap [b]', b)
        call mfree('tortor_newmap [c]', c)
        call mfree('tortor_newmap [d]', d)
        call mfree('tortor_newmap [dx]', dx)
        call mfree('tortor_newmap [dy]', dy)
        call mfree('tortor_newmap [dxy]', dxy)

    end subroutine tortor_newmap