Mountain glacier adjoint model

From Interactive System for Ice sheet Simulation
Jump to: navigation, search
!                           DISCLAIMER
! Removed by P.H.
!
!   This file was generated by TAF version 1.9.54
!
subroutine adkees1( xx, adxx, fc, adfc )
!******************************************************************
!******************************************************************
!** This routine was generated by Automatic differentiation.     **
!** FastOpt: Transformation of Algorithm in Fortran, TAF 1.9.54  **
!******************************************************************
!******************************************************************
!==============================================
! referencing used modules
!==============================================
use kees1_variables
 
!==============================================
! all entries are defined explicitly
!==============================================
implicit none
 
!==============================================
! declare parameters
!==============================================
real(kind=8) a
parameter ( a = 1.e-16 )
real(kind=8) g
parameter ( g = 9.2 )
real(kind=8) n
parameter ( n = 3 )
real(kind=8) rho
parameter ( rho = 920 )
real(kind=8) c
parameter ( c = 2*a*(rho*g)**n/(n+2)*1000.**n )
real(kind=8) dt
parameter ( dt = 1/12. )
real(kind=8) tend
parameter ( tend = 5000 )
 
!==============================================
! declare arguments
!==============================================
real(kind=8), intent(inout) :: adfc
real(kind=8), intent(inout) :: adxx(int(xend/dx)+1)
real(kind=8), intent(inout) :: fc
real(kind=8), intent(inout) :: xx(int(xend/dx)+1)
 
!==============================================
! declare local variables
!==============================================
real(kind=8) addiff(int(xend/dx))
real(kind=8) addiffi(:)
allocatable addiffi
real(kind=8) adh(int(xend/dx)+1)
real(kind=8) admb(int(xend/dx)+1)
real(kind=8) ads(int(tend/dt)+1,int(xend/dx)+1)
real(kind=8) b(int(xend/dx)+1)
real(kind=8) diff(int(xend/dx))
real(kind=8) diffi(:)
allocatable diffi
real(kind=8) h(int(xend/dx)+1)
integer i
real(kind=8) mb(int(xend/dx)+1)
integer :: nt = int(tend/dt)
integer :: nx = int(xend/dx)
real(kind=8) s(int(tend/dt)+1,int(xend/dx)+1)
integer t
real(kind=8) :: xarr(int(xend/dx)+1) = (/((i-1)*dx,i=1,int(xend/dx)+1)/)
 
!----------------------------------------------
! RESET LOCAL ADJOINT VARIABLES
!----------------------------------------------
addiff(:) = 0._8
adh(:) = 0._8
admb(:) = 0._8
ads(:,:) = 0._8
 
!----------------------------------------------
! ROUTINE BODY
!----------------------------------------------
!----------------------------------------------
! FUNCTION AND TAPE COMPUTATIONS
!----------------------------------------------
b = 1-0.0001*xarr
mb = 0.004-xarr*0.0002
s(1,:) = b
mb = mb+xx
do t = 2, nt
  s(t,1) = 0
  h = s(t-1,:)-b
  diff = c*((h(1:nx)+h(2:nx+1))/2)**(n+2)*abs((s(t-1,2:nx+1)-s(t-1,1:nx))/dx)**(n-1)
  s(t,2:nx) = s(t-1,2:nx)+dt/dx**2*(diff(2:nx)*(s(t-1,3:nx+1)-s(t-1,2:nx))-diff(1:nx-1)*(s(t-1,2:nx)-s(t-1,1:nx-1)))+mb(2:nx)*dt
  where (s(t,:) .lt. b)
    s(t,:) = b
  end where
end do
t = nt+1
s(t,1) = 0
h = s(t-1,:)-b
fc = 0.
do i = 1, size(h)
  fc = fc+h(i)*dx
end do
 
!----------------------------------------------
! ADJOINT COMPUTATIONS
!----------------------------------------------
allocate( addiffi(2:nx+1) )
allocate( diffi(2:nx+1) )
do i = 1, size(h)
  adh(i) = adh(i)+adfc*dx
end do
do t = nt+1, 2, -1
  s(t,1) = 0
  h = s(t-1,:)-b
  diff = c*((h(1:nx)+h(2:nx+1))/2)**(n+2)*abs((s(t-1,2:nx+1)-s(t-1,1:nx))/dx)**(n-1)
  s(t,2:nx) = s(t-1,2:nx)+dt/dx**2*(diff(2:nx)*(s(t-1,3:nx+1)-s(t-1,2:nx))-diff(1:nx-1)*(s(t-1,2:nx)-s(t-1,1:nx-1)))+mb(2:nx)*dt
  where (s(t,:) .lt. b)
    ads(t,:) = 0._8
  end where
  s(t,1) = 0
  addiff(2:nx) = addiff(2:nx)+ads(t,2:nx)*dt/dx**2*(s(t-1,3:nx+1)-s(t-1,2:nx))
  addiff(1:nx-1) = addiff(1:nx-1)-ads(t,2:nx)*dt/dx**2*(s(t-1,2:nx)-s(t-1,1:nx-1))
  admb(2:nx) = admb(2:nx)+ads(t,2:nx)*dt
  ads(t-1,3:nx+1) = ads(t-1,3:nx+1)+ads(t,2:nx)*dt/dx**2*diff(2:nx)
  ads(t-1,2:nx) = ads(t-1,2:nx)+ads(t,2:nx)*(1-dt/dx**2*(diff(2:nx)+diff(1:nx-1)))
  ads(t-1,1:nx-1) = ads(t-1,1:nx-1)+ads(t,2:nx)*dt/dx**2*diff(1:nx-1)
  ads(t,2:nx) = 0._8
  s(t,1) = 0
  diffi = (s(t-1,2:nx+1)-s(t-1,1:nx))/dx
  addiffi = addiff*c*((h(1:nx)+h(2:nx+1))/float(2))**(2+n)*(n-1)*abs(diffi)**((-2)+n)*sign(1._8,diffi)
  adh(2:nx+1) = adh(2:nx+1)+addiff*c/float(2)*(2+n)*((h(1:nx)+h(2:nx+1))/float(2))**(1+n)*abs(diffi)**(n-1)
  adh(1:nx) = adh(1:nx)+addiff*c/float(2)*(2+n)*((h(1:nx)+h(2:nx+1))/float(2))**(1+n)*abs(diffi)**(n-1)
  ads(t-1,2:nx+1) = ads(t-1,2:nx+1)+addiffi/dx
  ads(t-1,1:nx) = ads(t-1,1:nx)-addiffi/dx
  addiff = 0._8
  ads(t-1,:) = ads(t-1,:)+adh
  adh = 0._8
  ads(t,1) = 0._8
end do
adxx = adxx+admb
deallocate( addiffi )
deallocate( diffi )
 
end subroutine adkees1