C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ctdt.F,v 1.5 2005/03/28 23:49:49 heimbach Exp $
#include "COST_CPPOPTIONS.h"
subroutine COST_CTDT(
I myiter,
I mytime,
I mythid
& )
c ==================================================================
c SUBROUTINE cost_CTDT
c ==================================================================
c
c o Evaluate cost function contribution of CTD temperature data.
c
c started: Elisabeth Remy eremy@ucsd.edu 30-Aug-2000
c
c
c ==================================================================
c SUBROUTINE cost_CTDT
c ==================================================================
implicit none
c == global variables ==
#include "EEPARAMS.h"
#include "SIZE.h"
#include "GRID.h"
#include "DYNVARS.h"
#include "cal.h"
#include "ecco_cost.h"
#include "ctrl.h"
#include "ctrl_dummy.h"
#include "optim.h"
c == routine arguments ==
integer myiter
_RL mytime
integer mythid
c == local variables ==
integer bi,bj
integer i,j,k
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer nrec
integer irec
integer ilu
_RL fctile_ctdt
_RL fcthread_ctdt
_RL www (1-olx:snx+olx,1-oly:sny+oly)
_RL wtmp (1-olx:snx+olx,1-oly:sny+oly)
_RL tmpobs (1-olx:snx+olx,1-oly:sny+oly)
_RL tmpbar (1-olx:snx+olx,1-oly:sny+oly)
_RL cmask (1-olx:snx+olx,1-oly:sny+oly)
_RL spval
_RL spmax
character*(80) fnametheta
logical doglobalread
logical ladinit
character*(MAX_LEN_MBUF) msgbuf
cnew(
integer il
integer mody, modm
integer iyear, imonth
character*(80) fnametmp
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
spval = -1.8
spmax = 40.
c-- Read state record from global file.
doglobalread = .false.
ladinit = .false.
#ifdef ALLOW_CTDT_COST_CONTRIBUTION
if (optimcycle .ge. 0) then
ilu=ilnblnk( tbarfile )
write(fnametheta(1:80),'(2a,i10.10)')
& tbarfile(1:ilu),'.',optimcycle
endif
fcthread_ctdt = 0. _d 0
cnew(
mody = modelstartdate(1)/10000
modm = modelstartdate(1)/100 - mody*100
cnew)
c-- Loop over records.
do irec = 1,nmonsrec
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnametheta, tbar, irec,
& doglobalread, ladinit,
& optimcycle, mythid, xx_tbar_mean_dummy )
cnew(
iyear = mody + INT((modm-1+irec-1)/12)
imonth = 1 + MOD(modm-1+irec-1,12)
il=ilnblnk(ctdtfile)
write(fnametmp(1:80),'(2a,i4)')
& ctdtfile(1:il), '_', iyear
inquire( file=fnametmp, exist=exst )
if (.NOT. exst) then
write(fnametmp(1:80),'(a)') ctdtfile(1:il)
imonth = irec
endif
call MDSREADFIELD( fnametmp, cost_iprec, 'RL', nr, ctdtobs,
& imonth, mythid)
cnew)
c-- Loop over this thread's tiles.
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Loop over the model layers
fctile_ctdt = 0. _d 0
do k = 1,nr
c-- Determine the weights to be used.
do j = jmin,jmax
do i = imin,imax
cmask(i,j) = 1. _d 0
if (ctdtobs(i,j,k,bi,bj) .lt. spval .or.
& ctdtobs(i,j,k,bi,bj) .gt. spmax .or.
& ctdtobs(i,j,k,bi,bj) .eq. 0. ) then
cmask(i,j) = 0. _d 0
endif
cph(
cph print *, 'WARNING: SPECIFIC SETUP FOR ECCO'
cph below statement could be replaced by following
cph to make it independnet of Nr:
cph
cph if ( rC(K) .GT. -1000. ) then
cph)
c set cmask=0 in areas shallower than 1000m
if (
& (_hFacC(i,j,13,bi,bj) .ne. 0.).and.
& (_hFacC(i,j,k,bi,bj) .ne. 0.)) then
www(i,j) = cosphi(i,j,bi,bj)*cmask(i,j)
tmpobs(i,j) = ctdtobs(i,j,k,bi,bj)
tmpbar(i,j) = tbar(i,j,k,bi,bj)
wtmp(i,j) = wtheta2(i,j,k,bi,bj)
c-- The array ctdtobs contains CTD temperature.
fctile_ctdt = fctile_ctdt +
& (wtmp(i,j)*www(i,j))*
& (tmpbar(i,j)-tmpobs(i,j))*
& (tmpbar(i,j)-tmpobs(i,j))
if ( wtmp(i,j)*www(i,j) .ne. 0. )
& num_ctdt(bi,bj) = num_ctdt(bi,bj) + 1. _d 0
endif
enddo
enddo
enddo
c-- End of loop over layers.
fcthread_ctdt = fcthread_ctdt + fctile_ctdt
objf_ctdt(bi,bj) = objf_ctdt(bi,bj) + fctile_ctdt
#ifdef ECCO_VERBOSE
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
& ' COST_CTDT: irec,bi,bj = ',irec,bi,bj
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,d22.15)')
& ' COST_CTDT: cost function = ', fctile_ctdt
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
#endif
enddo
enddo
#ifdef ECCO_VERBOSE
c-- Print cost function for all tiles.
_GLOBAL_SUM_R8( fcthread_ctdt , myThid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8)')
& ' cost_CTDT: irec = ',irec
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a,d22.15)')
& ' global cost function value',
& ' ( CTD temp. ) = ',fcthread_ctdt
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
#endif
enddo
c-- End of second loop over records.
#else
c-- Do not enter the calculation of the CTD temperature contribution
c-- to the final cost function.
fctile_ctdt = 0. _d 0
fcthread_ctdt = 0. _d 0
crg
nrec = 1
crg
_BEGIN_MASTER( mythid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a)')
& ' cost_CTDT: no contribution of CTD temperature ',
& ' to cost function.'
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a,i9.8)')
& ' cost_CTDT: number of records that would have',
& ' been processed: ',nrec
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
_END_MASTER( mythid )
#endif
return
end