C $Header: /u/gcmpack/MITgcm/pkg/seaice/seaice_cost_concentration.F,v 1.11 2015/03/23 21:04:59 gforget Exp $
C $Name: $
#include "SEAICE_OPTIONS.h"
#ifdef ALLOW_ECCO
# include "ECCO_OPTIONS.h"
#endif
#ifdef ALLOW_COST
# include "COST_OPTIONS.h"
#endif
#ifdef ALLOW_CTRL
# include "CTRL_OPTIONS.h"
#endif
#ifdef ALLOW_AUTODIFF
# include "AUTODIFF_OPTIONS.h"
#endif
subroutine SEAICE_COST_CONCENTRATION(
& nnzbar, localbarfile, localbar, xx_localbar_mean_dummy,
& nnzobs, localobsfile, localobs, mult_local,
& nrecloc, localstartdate, localperiod,
& localmask, localweight,
& spminloc, spmaxloc, spzeroloc,
& objf_local, num_local,
& myiter, mytime, mythid )
implicit none
c ian fenty
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "PARAMS.h"
#ifdef ALLOW_CAL
# include "cal.h"
#endif
#ifdef ALLOW_COST
# ifdef ALLOW_CTRL
# include "optim.h"
# endif
# ifdef ALLOW_ECCO
# include "ecco.h"
# endif
# ifdef ALLOW_SEAICE
# include "SEAICE_COST.h"
# include "SEAICE_SIZE.h"
# include "SEAICE_PARAMS.h"
# endif
#endif
c == routine arguments ==
integer nnzbar
integer nnzobs
integer nrecloc
integer myiter
integer mythid
integer localstartdate(4)
_RL localbar (1-olx:snx+olx,1-oly:sny+oly,nnzbar,nsx,nsy)
_RL localobs (1-olx:snx+olx,1-oly:sny+oly,nnzobs,nsx,nsy)
_RL localweight (1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
#ifdef VARIABLE_SMRAREA_WEIGHT
_RL localModWeight(1-olx:snx+olx,1-oly:sny+oly,nsx,nsy)
_RL areaSigma
#endif
_RS localmask (1-olx:snx+olx,1-oly:sny+oly,nr,nsx,nsy)
_RL xx_localbar_mean_dummy
_RL mult_local
_RL mytime
_RL localperiod
_RL spminloc
_RL spmaxloc
_RL spzeroloc
_RL objf_local(nsx,nsy)
_RL num_local(nsx,nsy)
character*(MAX_LEN_FNAM) localbarfile
character*(MAX_LEN_FNAM) localobsfile
#if (defined (ALLOW_ECCO) defined (ALLOW_COST))
#if (defined(ALLOW_SEAICE_COST_SMR_AREA) defined(ALLOW_COST_ICE))
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 il
integer localrec
integer obsrec
logical doglobalread
logical ladinit
_RL spval
parameter (spval = -9999. _d 0 )
_RL localwww
_RL localcost
_RL junk
_RL cmask (1-olx:snx+olx,1-oly:sny+oly,nnzobs)
character*(128) fname1, fname2
cnew(
_RL daytime
_RL diffsecs
integer dayiter
integer daydate(4)
integer difftime(4)
integer middate(4)
integer yday, ymod
integer md, dd, sd, ld, wd
logical exst
cnew)
c == external functions ==
integer ilnblnk
external
c == end of interface ==
jtlo = mybylo(mythid)
jthi = mybyhi(mythid)
itlo = mybxlo(mythid)
ithi = mybxhi(mythid)
jmin = 1
jmax = sny
imin = 1
imax = snx
c-- Initialise local variables.
localwww = 0. _d 0
do bj = jtlo,jthi
do bi = itlo,ithi
objf_local(bi,bj) = 0. _d 0
num_local(bi,bj) = 0. _d 0
enddo
enddo
c-- First, read tiled data.
doglobalread = .false.
ladinit = .false.
write(fname1(1:128),'(80a)') ' '
il=ilnblnk( localbarfile )
write(fname1(1:128),'(2a,i10.10)')
& localbarfile(1:il),'.',optimcycle
cph if ( .NOT. ( mult_local.EQ.0. .OR. localobsfile.EQ.' ' ) ) then
if ( .NOT. ( localobsfile.EQ.' ' ) ) then
c-- Loop over records for the second time.
do irec = 1, nrecloc
if ( nnzbar .EQ. 1 ) then
call ACTIVE_READ_XY( fname1, localbar, irec, doglobalread,
& ladinit, optimcycle, mythid,
& xx_localbar_mean_dummy )
else
call ACTIVE_READ_XYZ( fname1, localbar, irec, doglobalread,
& ladinit, optimcycle, mythid,
& xx_localbar_mean_dummy )
endif
cnew(
obsrec = irec
daytime = FLOAT(secondsperday*(irec-1)) + modelstart
dayiter = hoursperday*(irec-1)+modeliter0
call CAL_GETDATE( dayiter, daytime, daydate, mythid )
call CAL_CONVDATE( daydate,yday,md,dd,sd,ld,wd,mythid )
ymod = smrareastartdate(1)/10000
#ifdef SEAICE_DEBUG
print *,'-- Cost seaice concentration'
print *,'daydate ', daydate
print *,'localobsfile: ', localobsfile
print *,'nrecloc ', nrecloc
print *,'obsrec,daytime ', obsrec,daytime
print *,'dayiter ', dayiter
print *,'yday : ', yday
print *,'md,dd,sd ', md,dd,sd
print *,'ld,wd ', ld,wd
print *,'loclstrtdte(1) ', localstartdate(1)
print *,'ymod, yday ', ymod,yday
print *,'smrarstrtdt(1) ', smrareastartdate(1)
print *,'smrarstartdate ', smrareastartdate
#endif /* SEAICE_DEBUG */
if ( ymod .EQ. yday ) then
middate(1) = smrareastartdate(1)
else
middate(1) = yday*10000+100+1
endif
middate(2) = 0
middate(3) = modelstartdate(3)
middate(4) = modelstartdate(4)
call CAL_TIMEPASSED( middate, daydate, difftime, mythid )
call CAL_TOSECONDS( difftime, diffsecs, mythid )
localrec = int(diffsecs/localperiod) + 1
#ifdef SEAICE_DEBUG
print *,'middate(1,2) ',middate(1),middate(2)
print *,'middate(3,4) ', middate(3),middate(4)
print *,'difftime,diffsecs',difftime,diffsecs
print *,'localrec ',localrec
#endif
il=ilnblnk(localobsfile)
write(fname2(1:128),'(2a,i4)')
& localobsfile(1:il), '_', yday
inquire( file=fname2, exist=exst )
#ifdef SEAICE_DEBUG
print *,'fname2',fname2
#endif
if (.NOT. exst) then
write(fname2(1:128),'(a)') localobsfile(1:il)
localrec = obsrec
#ifdef SEAICE_DEBUG
print *,'localrec ', localrec
print *,'not exist'
#endif
endif
if ( localrec .GT. 0 ) then
#ifdef SEAICE_DEBUG
print *,'calling mdsreadfile',fname2,localrec
#endif
call MDSREADFIELD( fname2, cost_iprec, cost_yftype, nnzobs,
& localobs, localrec, mythid )
else
do bj = jtlo,jthi
do bi = itlo,ithi
do k = 1,nnzobs
do j = jmin,jmax
do i = imin,imax
localobs(i,j,k,bi,bj) = spval
enddo !i
enddo !j
enddo !k
enddo !bi
enddo !bi
endif !obs rec
#ifdef VARIABLE_SMRAREA_WEIGHT
do bj = jtlo,jthi
do bi = itlo,ithi
do k = 1,nnzobs
do j = jmin,jmax
do i = imin,imax
cif set the new weight equal to the old weight
localModWeight(i,j,bi,bj) =
& localweight(i,j,bi,bj)
cif as long we the weight here is not zero we can continue
if (localweight(i,j,bi,bj) .GT. 0. _d 0) THEN
cif back out the original sigma for this location
areaSigma = 1. _d 0/sqrt(localweight(i,j,bi,bj))
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccX
IF (localobs(i,j,k,bi,bj) .eq. 0. _d 0 ) THEN
areaSigma = areaSigma * 0.85 _d 0
ELSEIF ((localobs(i,j,k,bi,bj).gt.0. _d 0 ) .and.
& (localobs(i,j,k,bi,bj).lt.0.15 _d 0)) THEN
areaSigma = areaSigma * 1.2 _d 0
ELSEIF ((localobs(i,j,k,bi,bj).ge.0.15 _d 0) .and.
& (localobs(i,j,k,bi,bj).le.0.25 _d 0)) THEN
areaSigma = areaSigma * 1.1 _d 0
ENDIF
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccX
cif reconstruct the weight = sigma^(-2)
localModWeight(i,j,bi,bj) =
& 1. _d 0 / (areaSigma*areaSigma)
cif if the local weight here is somehow 0,
cif do not divide by its square root but say something
c else
c print *,'costg : localweight <= 0 ',i,j
endif
#ifdef SEAICE_DEBUG
if ((i == SEAICE_debugPointX) .and.
& (j == SEAICE_debugPointY)) then
print '(A,2i4,4(1x,1P2E15.3))',
& 'costg i j obs,locWeight,locModWt,areaigma ',i,j,
& localobs(i,j,k,bi,bj), localweight(i,j,bi,bj),
& localModWeight(i,j,bi,bj),
& areaSigma
endif
#endif
C seaice debug
enddo
enddo
enddo
enddo
enddo
#endif
C variable smrarea weight
#ifdef SEAICE_DEBUG
do bj = jtlo,jthi
do bi = itlo,ithi
do k = 1,nnzobs
do i = imin,imax
do j = jmin,jmax
if (localobs(i,j,k,bi,bj) .LT. -1) THEN
print *,'obs rec date: ', -localobs(i,j,1,bi,bj)
endif
enddo
enddo
enddo
enddo
enddo
#endif
do bj = jtlo,jthi
do bi = itlo,ithi
localcost = 0. _d 0
c-- Determine the mask on weights
do k = 1,nnzobs
do j = jmin,jmax
do i = imin,imax
cmask(i,j,k) = cosphi(i,j,bi,bj)*localmask(i,j,k,bi,bj)
if ( localobs(i,j,k,bi,bj) .lt. spminloc .or.
& localobs(i,j,k,bi,bj) .gt. spmaxloc .or.
& localobs(i,j,k,bi,bj) .eq. spzeroloc ) then
cmask(i,j,k) = 0. _d 0
endif
enddo
enddo
enddo
c--
do k = 1,nnzobs
do j = jmin,jmax
do i = imin,imax
localwww = localweight(i,j,bi,bj)*cmask(i,j,k)
#ifdef VARIABLE_SMRAREA_WEIGHT
localwww = localModWeight(i,j,bi,bj)*cmask(i,j,k)
#endif
junk = ( localbar(i,j,k,bi,bj) -
& localobs(i,j,k,bi,bj) )
localcost = localcost + junk*junk*localwww
#ifdef SEAICE_DEBUG
if ((i == SEAICE_debugPointX) .and.
& (j == SEAICE_debugPointY)) then
print '(A,2i4,2(1x,1P2E15.3))',
& 'costg i j bar, obs ',i,j,
& localbar(i,j,k,bi,bj),
& localobs(i,j,k,bi,bj)
print '(A,2i4,2(1x,1P2E15.3))',
& 'costg i j bar-obs,wgt,loCost ',i,j,
& junk,
& localwww,
& junk*junk*localwww
endif
#endif
if ( localwww .ne. 0. )
& num_local(bi,bj) = num_local(bi,bj) + 1. _d 0
enddo
enddo
enddo
objf_local(bi,bj) = objf_local(bi,bj) + localcost
enddo
enddo
enddo
c-- End of second loop over records.
c-- End of mult_local or localobsfile
endif
#endif /* ifdef ALLOW_SEAICE_COST_SMR_AREA */
#endif /* ifdef ALLOW_ECCO */
return
end