```
module force
use size
!c     forcing fields for momentum equations
!un      common /force/ uforce, vforce
real uforce(nx,ny)
real vforce(nx,ny)
!c
end module

module data
use size

!c
!c     data arrays
!c
!un      common /data/ eta_data, u_data, v_data, depth_data,
!un     &     eta_data_time, zonal_transport_data, nedt
integer nedt, nedtmax
parameter ( nedtmax = 1000 )
real eta_data(nx,ny), u_data(nx,ny), v_data(nx,ny)
real depth_data(nx,ny)
real eta_data_time(nedtmax)
real zonal_transport_data

end module
module weights
use size
!c
!c     weight matices are all diagonal for now
!c
!un      common /weight_factors/ wf_depth, wf_eta, wf_u, wf_v,
real wf_depth, wf_eta, wf_u, wf_v, wf_lapldepth, wf_graddepth
real wf_zonal_transport
!c
!un      common /weights/ weight_depth, weight_eta, weight_u, weight_v,
real weight_depth(nx,nx,ny,ny)
real weight_eta(nx,ny)
real weight_u(nx,ny)
real weight_v(nx,ny)
real weight_lapldepth(nx,ny)
real weight_zonal_transport
!c
end module
module pfields
use size
!c     parameter fields, like depth, grid etc.
!un      common /pfields/ fcoriu, fcoriv, depth, frict,
!un     &     hu, hv, invhu, invhv
!c     coriolis parameter, depth, and friction coefficient
!c     defined at eta grid point
real fcoriu(nx,ny)
real fcoriv(nx,ny)
real depth(0:nx+1,0:ny+1)
real hu(0:nx+1,0:ny+1), hv(0:nx+1,0:ny+1)
real invhu(0:nx+1,0:ny+1), invhv(0:nx+1,0:ny+1)
real frict(0:nx+1,0:ny+1)
!c
!un      common /geom/ x, y, dx, dy, rx, ry, hy
real x(0:nx+1), y(0:ny+1)
real dx(0:nx), dy(0:ny)
real rx(0:ny+1), ry, hy(0:ny+1)
!c
!c
!un      common /ini/ inidepth, uini, vini, etaini
real inidepth(nx,ny)
real uini(nx,ny)
real vini(nx,ny)
real etaini(nx,ny)
!c
!un      common /scales/ scaledepth, scaleu, scalev, scaleeta
real scaledepth(nx,ny)
real scaleu(nx,ny), scalev(nx,ny), scaleeta(nx,ny)

end module

module vars
use size
!c     dynamic variables velocity and sea-surface height
!un      common /vars/ u, v, eta
real u(0:nx+1,0:ny+1)
real v(0:nx+1,0:ny+1)
#ifndef NT_DO_NOT_BREAK_MFEF90
real eta(0:nx+1,0:ny+1)
#endif
!c
end module

module size
!c
!c     size of the domain
!c
integer nx, ny
parameter ( nx = 20, ny = 20 )
end module

module parms
!c
!c     basic parameters
!c
real g, earth, pi, deg2rad, rho0, invrho0
parameter ( g = 9.81, earth = 6371000 )
parameter ( pi = 3.14159265358979323844, deg2rad = pi/180. )
parameter ( rho0 = 1028., invrho0 = 1./rho0 )
!c
!c     ``variable'' parameters
!c
!un      common /timeparameter/ dt, start_time, dt_dump, nt, ntspinup
real dt, start_time, dt_dump
integer nt, ntspinup
!c
!un      common /iterpar/ iteration
integer iteration
!c
!un      common /parms/ om, f0, beta, rini, ah, xstart, ystart
real om
real f0, beta
real rini
real ah
real xstart, ystart
!c
!c     flags
!c
!un      common /flags/ xperiodic, yperiodic, spherical, cartesian,
logical xperiodic, yperiodic, spherical, cartesian, quadfric
logical suppressio, fullio
!c
!c     names
!c
!un      common /names/ foutname, runname, depthfile, forcingfile,
!un     &     uinifile, vinifile, etainifile, ncdatafile, ncrestartfile
character*(80) foutname
character*(80) runname
character*(80) depthfile, forcingfile
character*(80) uinifile, vinifile, etainifile
character*(80) ncdatafile, ncrestartfile
!c
end module

subroutine forward_model( nctrl, xc, cost_final )

use size
use parms
use vars
use pfields
use weights
use data
use force

implicit none

c ===============================================================================

c     interface
integer nctrl
real xc(nctrl)
real cost_final
c     end interface

c     locals
integer time_index, ntotal
integer it, jt, nio
logical calc_cost
real cost_d, cost_sd, cost_gd, cost, time
c     zonal volume transport
real zonal_transport
c     checkpoint frequency for tamc code
integer nouter, ninner
parameter ( ninner = 1000 )
real routin
cph new from inlining
integer ix, iy, jx, jy
integer k
real frictu, frictv
real fv, fu
cnt      logical is_eta_data_time
cnt      external is_eta_data_time
cnt
logical testResult
real one
parameter ( one = 1. )

c     we will need to store the dynamic variables u, v, and eta, so they
c     do not have to be recomputed during the adjoint calculations, so far
c     use two level checkpointing, intialize the outer tape (first tape
c     level) here to avoid unecessary recomputation of the forward model.
c
c
c     determine whether to calculate cost function ( can save some time )
c
&     .or. calc_hess ) then
calc_cost = .true.
else
calc_cost = .false.
end if
c
c     initial local variables
c
cost_d = 0.
cost_sd = 0.
cost_gd = 0.
cost = 0.

c
c     initialise other stuff
do iy = 1, ny
do ix = 1, nx
depth(ix,iy) = 0.
invhu(ix,iy) = 0.
invhv(ix,iy) = 0.
end do
end do

c
cph(
cph      call map_from_control_vector( nctrl, xc )
k = 0
do iy = 1, ny
do ix = 1, nx
if ( etamask(ix,iy) .ne. 0 ) then
k = k+1
depth(ix,iy) = scaledepth(ix,iy)*( one + xc(k) )
c               print *, ix, iy, xc(k), depth(ix,iy)
c               depth(ix,iy) = scaledepth(ix,iy) + xc(k)
c               depth(ix,iy) = xc(k)
end if
end do
end do
if ( k .ne. nctrl ) then
print *, 'map_from_control_vector: ',
&        'dimensions of control vector are wrong'
print *, k, ' should be ', nctrl
stop
end if
cph)
c
c     calculate depth at the faces of the grid points ( u and v points )
c
cph(
cph      call calc_depth_uv
c     create depth at u and v points
do iy = 1, ny+1
do ix = 1, nx+1

if ( hu(ix,iy) .ne. 0. ) then
invhu(ix,iy) = 1./hu(ix,iy)
end if
if ( hv(ix,iy) .ne. 0. ) then
invhv(ix,iy) = 1./hv(ix,iy)
end if
end do
end do

c\$\$\$      open(20,file='huhv.dat',form='formatted',recl=102*20)
c\$\$\$      write(20,'(102E20.12)') ((hu(ix,iy),ix=0,nx+1),iy=ny+1,0,-1)
c\$\$\$      write(20,'(102E20.12)') ((hu(ix,iy),ix=0,nx+1),iy=ny+1,0,-1)
c\$\$\$      close(20)
cph)
c
c     initialize dynamic variables
c
cph(
cph      call initial_values
do iy = 1, ny
do ix = 1, nx
#ifndef NT_DO_NOT_BREAK_MFEF90
#endif
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
u(nx+1,iy) = u(1,iy)
v(0,iy)    = v(nx,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
eta(0,iy)  = eta(nx,iy)
#endif
end do
end if
if ( yperiodic ) then
do ix = 0, nx+1
u(ix,0)    = u(ix,ny)
v(ix,ny+1) = v(ix,1)
#ifndef NT_DO_NOT_BREAK_MFEF90
eta(ix,0)  = eta(ix,0)
#endif
end do
end if
cph)
c
c     initialize I/O
if ( .not. suppressio )  then
cph(
cph these contains netcdf stuff. Can be avoided for OpenAD???
call ini_io
call pfields_io
cph)
end if
c     first cost function contribution is the time independend variable
c     depth
if ( calc_cost ) then
cph(
cph         call cost_depth( cost_d )
do iy = 1, ny
do jy = 1, ny
do ix = 1, nx
do jx = 1, nx
cost_d = cost_d
&                 + .5*(depth(ix,iy)-depth_data(ix,iy))
&                 *weight_depth(ix,jx,iy,jy)
&                 *(depth(jx,jy)-depth_data(jx,jy))
end do
end do
end do
end do
cph)
c         call smoothness_lapldepth( cost_sd )
end if
c
c     start time stepping
c
nio = 0
time_index = 0
time = 0.
c
c     total number of time steps include spin up
c
ntotal = nt+ntspinup
c\$\$\$c
c\$\$\$c     store the dynamic variables u, v, and eta, so they do not have
c\$\$\$c     to be recomputed during the adjoint calculations, so far use
c\$\$\$c     two level checkpointing
c\$\$\$c
c
c     start integration over assimilation interval, determine the number
c     of cycles for the outer loop of check pointing, round to the nearest
c     integer larger than routin (ratio of total to inner loop cycles)
c
routin = real(ntotal)/real(ninner)
if ( routin .ne. ntotal/ninner ) then
nouter = int( routin ) + 1
else
nouter = int( routin )
end if
if ( fullio .and. .not. suppressio ) then
print *, 'number of outer loops',
&        ' = number of tape records = ', nouter
end if
time_index = 0
time = start_time+time_index*dt
if ( .not. suppressio ) then
nio = nio + 1
print *, 'Writing Time Step ', time_index
cph(
cph this contains netcdf stuff
cph         call state_io( time, nio )
cph)
end if
do it = 1, nouter
c
cadj store u,v,eta = tape_level_1, key=it
cadj init tape_level_2 = common, ninner
c
do jt = 1, ninner
time_index = (it-1)*ninner+jt
time = start_time+time_index*dt
if ( time_index .le. ntotal ) then
c

cph(
cph               call time_step( time_index )
!cadj init tape_level_aux = common, 1
c     alternate evaluation of momentum equations, because coriolis is
c     treated explicitly
if ( mod(time_index,2) .ne. 0 ) then

cph(
cph 2nd level inlining
cph                  call umomentum
do iy = 1, ny
do ix = 1, nx
frictu = .5*(frict(ix,iy)+frict(ix-1,iy))*invhu(ix,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
&           /(rx(iy)*.5*(dx(ix)+dx(ix-1)))
#endif
c     factor .25 is contained in fcoriu
fv = fcoriu(ix,iy)*( v(ix,iy)  + v(ix-1,iy)
&                         + v(ix,iy+1)+ v(ix-1,iy+1))
&           (1-.5*frictu*dt)*u(ix,iy) - g*dt*gradetau + dt*fv
&           + dt*uforce(ix,iy)*invhu(ix,iy)
&           )
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
u(nx+1,iy) = u(1,iy)
end do
end if
if ( yperiodic ) then
c     for coriolis term in v equation
do ix = 0, nx+1
u(ix,0) = u(ix,ny)
end do
end if
cph)
cph(
cph 2nd level inlining
cph                  call vmomentum
do iy = 1, ny
do ix = 1, nx
frictv = .5*(frict(ix,iy)+frict(ix,iy-1))*invhv(ix,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
&           /(ry*.5*(dy(iy)+dy(iy-1)))
#endif
c     factor .25 is contained in fcoriv
fu = fcoriv(ix,iy)*( u(ix,iy)  + u(ix+1,iy)
&                         + u(ix,iy-1)+ u(ix+1,iy-1))
&           (1-.5*frictv*dt)*v(ix,iy) - g*dt*gradetav - dt*fu
&           + dt*vforce(ix,iy)*invhv(ix,iy)
&           )
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
c     for coriolis term in u equation
v(0,iy) = v(nx,iy)
end do
end if
if ( yperiodic ) then
do ix = 0, nx+1
v(ix,ny+1) = v(ix,1)
end do
end if
cph)
else

cph(
cph 2nd level inlining
cph                  call vmomentum
do iy = 1, ny
do ix = 1, nx
frictv = .5*(frict(ix,iy)+frict(ix,iy-1))*invhv(ix,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
&           /(ry*.5*(dy(iy)+dy(iy-1)))
#endif
c     factor .25 is contained in fcoriv
fu = fcoriv(ix,iy)*( u(ix,iy)  + u(ix+1,iy)
&                         + u(ix,iy-1)+ u(ix+1,iy-1))
&           (1-.5*frictv*dt)*v(ix,iy) - g*dt*gradetav - dt*fu
&           + dt*vforce(ix,iy)*invhv(ix,iy)
&           )
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
c     for coriolis term in u equation
v(0,iy) = v(nx,iy)
end do
end if
if ( yperiodic ) then
do ix = 0, nx+1
v(ix,ny+1) = v(ix,1)
end do
end if
cph)
cph(
cph 2nd level inlining
cph                  call umomentum
do iy = 1, ny
do ix = 1, nx
frictu = .5*(frict(ix,iy)+frict(ix-1,iy))*invhu(ix,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
&           /(rx(iy)*.5*(dx(ix)+dx(ix-1)))
#endif
c     factor .25 is contained in fcoriu
fv = fcoriu(ix,iy)*( v(ix,iy)  + v(ix-1,iy)
&                         + v(ix,iy+1)+ v(ix-1,iy+1))
&           (1-.5*frictu*dt)*u(ix,iy) - g*dt*gradetau + dt*fv
&           + dt*uforce(ix,iy)*invhu(ix,iy)
&           )
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
u(nx+1,iy) = u(1,iy)
end do
end if
if ( yperiodic ) then
c     for coriolis term in v equation
do ix = 0, nx+1
u(ix,0) = u(ix,ny)
end do
end if
cph)
end if
c
c     continuity equation calculates eta
c
cph(
cph 2nd level inlining
cph               call continuity
c     eta equation
do iy = 1, ny
do ix = 1, nx
#ifndef NT_DO_NOT_BREAK_MFEF90
&           (hu(ix+1,iy)*u(ix+1,iy) - hu(ix,iy)*u(ix,iy))
&           /(rx(iy)*dx(ix))
&           + ( hy(iy+1)* hv(ix,iy+1)*v(ix,iy+1)
&             - hy(iy)  * hv(ix,iy)  *v(ix,iy)  )
&           /(ry*dy(iy))
&           )
#endif
end do
end do
c     handle domain boundaries
if ( xperiodic ) then
do iy = 0, ny+1
#ifndef NT_DO_NOT_BREAK_MFEF90
eta(0,iy) = eta(nx,iy)
#endif
end do
end if
if ( yperiodic ) then
do ix = 0, nx+1
#ifndef NT_DO_NOT_BREAK_MFEF90
eta(ix,0) = eta(ix,ny)
#endif
end do
end if
cph)

if ( ( time_index .gt. ntspinup ) .and. calc_cost ) then
cph(
cph                  call cost_function( time , cost )
c
c
cnt: convert to sub call
call is_eta_data_time( time, testResult )
if ( testResult ) then
cph(
cph no observed eta are used for this experiment
cph i.e. simplified cost function!
cph)
end if
c
c     calculate cost function, if there is only one data slab, or if
c     the time step matches a data time
c
cnt: convert to sub call
call is_eta_data_time( time, testResult )
if ( nedt .eq. -1 .or. testResult ) then
c
c     calculate cost function terms
c
c     zonal transport
cph(
cph 2nd level inlining
cph         call calc_zonal_transport( zonal_transport )
ix = 6
zonal_transport = 0.
do iy = 1, ny
zonal_transport = zonal_transport
&        + u(ix,iy)*dy(iy)*hu(ix,iy)
c     &       + u(ix,iy)*dy(iy)*( eta(ix-1,iy)+eta(ix,iy) )
end do
cph)
c\$\$\$         cf = cf + zonal_transport*wf_zonal_transport
cost = cost + .5*(zonal_transport-zonal_transport_data)**2
&        *weight_zonal_transport
c
do iy = 1, ny
do ix = 1, nx
c     sea surface height
cph               cost = cost + .5*(eta(ix,iy)-eta_data(ix,iy))**2
cph     &              *weight_eta(ix,iy)
c     velocity
cph               cost = cost + .5*(u(ix,iy)-u_data(ix,iy))**2
cph     &              *weight_u(ix,iy)
cph               cost = cost + .5*(v(ix,iy)-v_data(ix,iy))**2
cph     &              *weight_v(ix,iy)
end do
end do
end if
cph)
end if
c
c     after stepping the model forward, store the dynamic variables
c
if ( fullio .and. .not. suppressio .and.
&              ( mod(real(time_index)*dt,dt_dump) .eq. 0. ) ) then
nio = nio + 1
print *, 'Writing Time Step ', time_index
cph(
cph this contains netcdf stuff
cph                  call state_io( time, nio )
cph)
end if
end if
end do
c
end do
c\$\$\$c     overturning
c\$\$\$      zonal_transport = 0.
c\$\$\$      call calc_overturning( zonal_transport )
c\$\$\$      print *, 'overturning transport = ', zonal_transport*1e-6, ' Sv'
c\$\$\$      cost = zonal_transport
c     zonal transport
zonal_transport = 0.
cph(
cph      call calc_zonal_transport( zonal_transport )
ix = 6
zonal_transport = 0.
do iy = 1, ny
zonal_transport = zonal_transport
&        + u(ix,iy)*dy(iy)*hu(ix,iy)
c     &       + u(ix,iy)*dy(iy)*( eta(ix-1,iy)+eta(ix,iy) )
end do
cph)
print *, 'zonal volume transport = ', zonal_transport*1e-6, ' Sv'
c
c     save cost function value (to trick TAMC)
c
cph(
cph cost function is just zonal transport
cph      cost_final = cost + cost_d + cost_sd + cost_gd
cost_final = zonal_transport
cph)

cph      if ( iteration .ge. 0 ) then
open(10,file='cost.txt',form='formatted',position='append')
write(10,'(I5,5E15.8)') iteration, cost_d, cost_sd, cost_gd,
&        cost, cost_final
close(10)
iteration = iteration + 1
cph      end if

return
end !subroutine forward_model
```