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,
!un & wf_zonal_transport, wf_lapldepth, wf_graddepth
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,
!un & weight_lapldepth, weight_graddepth, weight_zonal_transport
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_graddepth(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
!un common /mask/ umask, vmask, etamask
real umask(0:nx+1,0:ny+1)
real vmask(0:nx+1,0:ny+1)
real etamask(0:nx+1,0:ny+1)
!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,
!un & quadfric, suppressio, fullio,
!un & initial_grad, grad_check, optimize, calc_hess
logical xperiodic, yperiodic, spherical, cartesian, quadfric
logical suppressio, fullio
logical initial_grad, grad_check, optimize, calc_hess
!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 gradetau, gradetav
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
cadj init tape_level_1 = 'state'
c
c determine whether to calculate cost function ( can save some time )
c
if ( initial_grad .or. grad_check .or. optimize
& .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
hu(ix,iy) = depth(ix,iy)*umask(ix,iy)
hv(ix,iy) = depth(ix,iy)*vmask(ix,iy)
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
u(ix,iy) = uini(ix,iy)*umask(ix,iy)
v(ix,iy) = vini(ix,iy)*vmask(ix,iy)
#ifndef NT_DO_NOT_BREAK_MFEF90
eta(ix,iy) = etaini(ix,iy)*etamask(ix,iy)
#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 )
c call smoothness_graddepth( cost_gd )
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$$$cadj init tape_level_1 = 'state'
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
cadj store u,v,eta = tape_level_2
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
cadj store u = tape_level_2
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
gradetau = (eta(ix,iy)-eta(ix-1,iy))
& /(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))
u(ix,iy) = umask(ix,iy)/(1+0.5*frictu*dt)*(
& (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
gradetav = (eta(ix,iy)-eta(ix,iy-1))
& /(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))
v(ix,iy) = vmask(ix,iy)/(1+0.5*frictv*dt)*(
& (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
cadj store v = tape_level_2
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
gradetav = (eta(ix,iy)-eta(ix,iy-1))
& /(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))
v(ix,iy) = vmask(ix,iy)/(1+0.5*frictv*dt)*(
& (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
gradetau = (eta(ix,iy)-eta(ix-1,iy))
& /(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))
u(ix,iy) = umask(ix,iy)/(1+0.5*frictu*dt)*(
& (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
eta(ix,iy) = eta(ix,iy) - etamask(ix,iy)*dt*(
& (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 load the data
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 call read_eta_data( time )
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