C $Header: /u/gcmpack/MITgcm/pkg/ctrl/ctrl_cost_gen.F,v 1.9 2015/10/29 03:43:59 gforget Exp $
C $Name: $
c ----------------------------------------------------------------
c --- ctrl_cost_gen2d
c --- ctrl_cost_gen3d
c ----------------------------------------------------------------
c ----------------------------------------------------------------
#include "CTRL_OPTIONS.h"
#ifdef ALLOW_ECCO
# include "ECCO_OPTIONS.h"
#endif
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: ctrl_cost_gen2d
C !INTERFACE:
subroutine CTRL_COST_GEN2D(
I startrec,
I endrec,
I xx_gen_file,
I xx_gen_dummy,
I xx_gen_period,
I xx_gen_weight,
I dodimensionalcost,
O num_gen_anom,
O objf_gen_anom,
#ifdef ECCO_CTRL_DEPRECATED
I xx_gen_wmean,
O num_gen_mean,
O objf_gen_mean,
O objf_gen_smoo,
I xx_gen_remo_intercept,
I xx_gen_remo_slope,
#endif /* ECCO_CTRL_DEPRECATED */
I xx_gen_mask,
I myThid
& )
C !DESCRIPTION: \bv
C Generic routine for all 2D control penalty terms
C \ev
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
#ifdef ALLOW_CTRL
# include "ctrl.h"
# include "optim.h"
#endif
c == routine arguments ==
integer startrec
integer endrec
character*(MAX_LEN_FNAM) xx_gen_file
_RL xx_gen_dummy
_RL xx_gen_period
_RL xx_gen_weight(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
logical dodimensionalcost
_RL num_gen_anom(nsx,nsy)
_RL objf_gen_anom(nsx,nsy)
#ifdef ECCO_CTRL_DEPRECATED
_RL xx_gen_wmean
_RL num_gen_mean(nsx,nsy)
_RL num_gen_smoo(nsx,nsy)
_RL objf_gen_mean(nsx,nsy)
_RL objf_gen_smoo(nsx,nsy)
_RL xx_gen_remo_intercept
_RL xx_gen_remo_slope
#endif /* ECCO_CTRL_DEPRECATED */
_RS xx_gen_mask(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
integer myThid
#ifdef ALLOW_CTRL
c == local variables ==
integer bi,bj
integer i,j,kk
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer nrec
integer irec
integer ilfld
_RL fctile
_RL fctilem
_RL tmpx
_RL lengthscale
#ifdef ECCO_CTRL_DEPRECATED
_RL xx_mean(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
#endif /* ECCO_CTRL_DEPRECATED */
logical doglobalread
logical ladinit
character*(80) fnamefld
c == external functions ==
integer ilnblnk
external
CEOP
jtlo = mybylo(myThid)
jthi = mybyhi(myThid)
itlo = mybxlo(myThid)
ithi = mybxhi(myThid)
jmin = 1
jmax = sny
imin = 1
imax = snx
lengthscale = 1. _d 0
c-- Read state record from global file.
doglobalread = .false.
ladinit = .false.
c Number of records to be used.
nrec = endrec-startrec+1
if (optimcycle .ge. 0) then
ilfld=ilnblnk( xx_gen_file )
write(fnamefld(1:80),'(2a,i10.10)')
& xx_gen_file(1:ilfld),'.',optimcycle
endif
c-- >>> Loop 1 to compute mean forcing:
do bj = jtlo,jthi
do bi = itlo,ithi
num_gen_anom(bi,bj) = 0. _d 0
objf_gen_anom(bi,bj) = 0. _d 0
#ifdef ECCO_CTRL_DEPRECATED
do j = jmin,jmax
do i = imin,imax
xx_mean(i,j,bi,bj) = 0. _d 0
enddo
enddo
num_gen_mean(bi,bj) = 0. _d 0
num_gen_smoo(bi,bj) = 0. _d 0
objf_gen_mean(bi,bj) = 0. _d 0
objf_gen_smoo(bi,bj) = 0. _d 0
#endif /* ECCO_CTRL_DEPRECATED */
enddo
enddo
#ifndef ECCO_CTRL_DEPRECATED
c-- >>> Loop over records.
do irec = 1,nrec
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY(
& fnamefld, tmpfld2d, irec, doglobalread,
& ladinit, optimcycle, myThid, xx_gen_dummy )
#else
CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif
c-- Loop over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Determine the weights to be used.
kk = 1
fctile = 0. _d 0
do j = jmin,jmax
do i = imin,imax
if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
tmpx = tmpfld2d(i,j,bi,bj)
IF ( dodimensionalcost ) THEN
fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
ELSE
fctile = fctile + tmpx*tmpx
ENDIF
if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
& num_gen_anom(bi,bj) = num_gen_anom(bi,bj)
& + 1. _d 0
endif
enddo
enddo
objf_gen_anom(bi,bj) = objf_gen_anom(bi,bj) + fctile
enddo
enddo
c-- End of loop over records.
enddo
#else /* ECCO_CTRL_DEPRECATED */
IF ( dodimensionalcost ) THEN
do irec = 1,nrec
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY(
& fnamefld, tmpfld2d, irec, doglobalread,
& ladinit, optimcycle, myThid, xx_gen_dummy )
#else
CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif
c-- Loop over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
do j = jmin,jmax
do i = imin,imax
xx_mean(i,j,bi,bj) = xx_mean(i,j,bi,bj)
& + tmpfld2d(i,j,bi,bj)
& - ( xx_gen_remo_intercept +
& xx_gen_remo_slope*(irec-1)*xx_gen_period )
enddo
enddo
enddo
enddo
enddo
if ( xx_gen_wmean .NE. 0. ) then
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Determine the weights to be used.
kk = 1
fctilem = 0. _d 0
do j = jmin,jmax
do i = imin,imax
xx_mean(i,j,bi,bj)
& = xx_mean(i,j,bi,bj)/float(nrec)
tmpx = xx_mean(i,j,bi,bj)/xx_gen_wmean
if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
#ifdef ALLOW_ECCO
if ( ABS(R_low(i,j,bi,bj)) .LT. 100. _d 0)
& tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100. _d 0
fctilem = fctilem + cosphi(i,j,bi,bj)*tmpx*tmpx
if ( cosphi(i,j,bi,bj) .ne. 0. _d 0)
& num_gen_mean(bi,bj) = num_gen_mean(bi,bj) + 1. _d 0
#else
fctilem = fctilem + tmpx*tmpx
num_gen_mean(bi,bj) = num_gen_mean(bi,bj) + 1. _d 0
#endif
endif
enddo
enddo
objf_gen_mean(bi,bj) = objf_gen_mean(bi,bj) + fctilem
enddo
enddo
endif
ENDIF !IF ( dodimensionalcost ) THEN
c-- >>> Loop 2 over records.
do irec = 1,nrec
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY(
& fnamefld, tmpfld2d, irec, doglobalread,
& ladinit, optimcycle, myThid, xx_gen_dummy )
#else
CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif
c-- Loop over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Determine the weights to be used.
kk = 1
fctile = 0. _d 0
do j = jmin,jmax
do i = imin,imax
if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
IF ( dodimensionalcost ) THEN
tmpx = tmpfld2d(i,j,bi,bj)
& - xx_mean(i,j,bi,bj)
& - ( xx_gen_remo_intercept +
& xx_gen_remo_slope*(irec-1)*xx_gen_period )
#ifdef ALLOW_ECCO
if ( ABS(R_low(i,j,bi,bj)) .LT. 100. _d 0 )
& tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100. _d 0
fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
& * cosphi(i,j,bi,bj)
#else
fctile = fctile + xx_gen_weight(i,j,bi,bj)*tmpx*tmpx
#endif
ELSE !IF ( dodimensionalcost ) THEN
tmpx = tmpfld2d(i,j,bi,bj)
fctile = fctile + tmpx*tmpx
ENDIF !IF ( dodimensionalcost ) THEN
#ifdef ALLOW_ECCO
if ( xx_gen_weight(i,j,bi,bj)
& *cosphi(i,j,bi,bj) .ne. 0. _d 0 )
#else
if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
#endif
& num_gen_anom(bi,bj) = num_gen_anom(bi,bj)
& + 1. _d 0
endif
enddo
enddo
objf_gen_anom(bi,bj) = objf_gen_anom(bi,bj) + fctile
enddo
enddo
c-- End of loop over records.
enddo
IF ( dodimensionalcost ) THEN
#ifdef ALLOW_SMOOTH_BC_COST_CONTRIBUTION
c-- >>> Loop 2 over records.
do irec = 1,nrec
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XY(
& fnamefld, tmpfld2d, irec, doglobalread,
& ladinit, optimcycle, myThid, xx_gen_dummy )
#else
CALL READ_REC_XY_RL( fnamefld, tmpfld2d, iRec, 1, myThid )
#endif
_EXCH_XY_RL(tmpfld2d, myThid)
c-- Loop over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Determine the weights to be used.
kk = 1
fctile = 0. _d 0
do j = jmin,jmax
do i = imin,imax
if (xx_gen_mask(i,j,kk,bi,bj) .ne. 0. _d 0) then
tmpx =
& ( tmpfld2d(i+2,j,bi,bj)-tmpfld2d(i+1,j,bi,bj) )
& *maskW(i+1,j,kk,bi,bj)*maskW(i+2,j,kk,bi,bj)
& + ( tmpfld2d(i+1,j,bi,bj)-tmpfld2d(i,j,bi,bj) )
& *maskW(i+1,j,kk,bi,bj)
& + ( tmpfld2d(i,j+2,bi,bj)-tmpfld2d(i,j+1,bi,bj) )
& *maskS(i,j+1,kk,bi,bj)*maskS(i,j+2,kk,bi,bj)
& + ( tmpfld2d(i,j+1,bi,bj)-tmpfld2d(i,j,bi,bj) )
& *maskS(i,j+1,kk,bi,bj)
#ifdef ALLOW_ECCO
if ( ABS(R_low(i,j,bi,bj)) .LT. 100. _d 0 )
& tmpx = tmpx*ABS(R_low(i,j,bi,bj))/100. _d 0
fctile = fctile
& + xx_gen_weight(i,j,bi,bj)*cosphi(i,j,bi,bj)
#else
fctile = fctile
& + xx_gen_weight(i,j,bi,bj)
#endif
* *0.0161 _d 0*lengthscale/4.0 _d 0
& *tmpx*tmpx
#ifdef ALLOW_ECCO
if ( xx_gen_weight(i,j,bi,bj)*cosphi(i,j,bi,bj)
& .ne. 0. _d 0 )
#else
if ( xx_gen_weight(i,j,bi,bj) .ne. 0. _d 0 )
#endif
& num_gen_smoo(bi,bj) = num_gen_smoo(bi,bj)
& + 1. _d 0
endif
enddo
enddo
objf_gen_smoo(bi,bj) = objf_gen_smoo(bi,bj) + fctile
enddo
enddo
c-- End of loop over records.
enddo
#endif /* ALLOW_SMOOTH_BC_COST_CONTRIBUTION */
ENDIF !IF ( dodimensionalcost ) THEN
#endif /* ECCO_CTRL_DEPRECATED */
#endif /* ALLOW_CTRL */
return
end
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
CBOP
C !ROUTINE: ctrl_cost_gen3d
C !INTERFACE:
subroutine CTRL_COST_GEN3D(
I xx_gen_file,
I xx_gen_dummy,
I xx_gen_weight,
I dodimensionalcost,
O num_gen,
O objf_gen,
I xx_gen_mask,
I myThid
& )
C !DESCRIPTION: \bv
C Generic routine for all 3D control penalty terms
C \ev
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#include "GRID.h"
#ifdef ALLOW_ECCO
# include "ecco.h"
#endif
#ifdef ALLOW_CTRL
# include "ctrl.h"
# include "optim.h"
#endif
c == routine arguments ==
character*(MAX_LEN_FNAM) xx_gen_file
_RL xx_gen_dummy
_RL xx_gen_weight(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
logical dodimensionalcost
_RL num_gen(nsx,nsy)
_RL objf_gen(nsx,nsy)
_RS xx_gen_mask(1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
INTEGER myThid
#ifdef ALLOW_CTRL
c == local variables ==
integer bi,bj
integer i,j,k
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer irec
integer ilfld
_RL tmpx
logical doglobalread
logical ladinit
character*(80) fnamefld
c == external functions ==
integer ilnblnk
external
CEOP
jtlo = mybylo(myThid)
jthi = mybyhi(myThid)
itlo = mybxlo(myThid)
ithi = mybxhi(myThid)
jmin = 1
jmax = sny
imin = 1
imax = snx
c-- Read state record from global file.
doglobalread = .false.
ladinit = .false.
if (optimcycle .ge. 0) then
ilfld = ilnblnk( xx_gen_file )
write(fnamefld(1:80),'(2a,i10.10)')
& xx_gen_file(1:ilfld),'.',optimcycle
endif
c-- >>> Loop 1 to compute mean forcing:
do bj = jtlo,jthi
do bi = itlo,ithi
num_gen(bi,bj) = 0. _d 0
objf_gen(bi,bj) = 0. _d 0
enddo
enddo
irec = 1
#ifdef ALLOW_AUTODIFF
call ACTIVE_READ_XYZ( fnamefld, tmpfld3d, irec, doglobalread,
& ladinit, optimcycle, myThid
& , xx_gen_dummy )
#else
CALL READ_REC_XYZ_RL( fnamefld, tmpfld3d, iRec, 1, myThid )
#endif
c-- Loop over this thread tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
num_gen(bi,bj) = 0. _d 0
objf_gen(bi,bj) = 0. _d 0
do k = 1,nr
do j = jmin,jmax
do i = imin,imax
if (xx_gen_mask(i,j,k,bi,bj) .ne. 0. _d 0) then
tmpx = tmpfld3d(i,j,k,bi,bj)
IF ( dodimensionalcost ) THEN
objf_gen(bi,bj) = objf_gen(bi,bj)
& + xx_gen_weight(i,j,k,bi,bj)
& *tmpx*tmpx
ELSE
objf_gen(bi,bj) = objf_gen(bi,bj) + tmpx*tmpx
ENDIF
if ( xx_gen_weight(i,j,k,bi,bj) .ne. 0. _d 0 )
& num_gen(bi,bj) = num_gen(bi,bj) + 1. _d 0
endif
enddo
enddo
enddo
enddo
enddo
#endif
return
end