C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_ctdsclim.F,v 1.4 2004/10/11 16:38:53 heimbach Exp $
#include "COST_CPPOPTIONS.h"
subroutine COST_CTDSCLIM(
I myiter,
I mytime,
I mythid
& )
c ==================================================================
c SUBROUTINE cost_Ctdsclim
c ==================================================================
c
c o Evaluate cost function contribution of salinity.
c
c started: Christian Eckert eckert@mit.edu 30-Jun-1999
c
c changed: Christian Eckert eckert@mit.edu 25-Feb-2000
c
c - Restructured the code in order to create a package
c for the MITgcmUV.
c
c changed: Patrick Heimbach heimbach@mit.edu 27-May-2000
c
c - set ladinit to .true. to initialise adsbar file
c
c ==================================================================
c SUBROUTINE cost_Ctdsclim
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 ==
_RS one_rs
parameter( one_rs = 1. )
integer bi,bj
integer i,j,k
integer itlo,ithi
integer jtlo,jthi
integer jmin,jmax
integer imin,imax
integer irec
integer levmon
integer levoff
integer ilctdsclim
_RL fctile
_RL fcthread
_RL cmask (1-olx:snx+olx,1-oly:sny+oly)
_RL spval
character*(80) fnamesalt
logical doglobalread
logical ladinit
character*(MAX_LEN_MBUF) msgbuf
#ifdef GENERIC_BAR_MONTH
integer mrec, nyears, iyear
#endif
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 = -9990.
c-- Read tiled data.
doglobalread = .false.
ladinit = .false.
#ifdef ALLOW_CTDSCLIM_COST_CONTRIBUTION
if (optimcycle .ge. 0) then
ilctdsclim = ilnblnk( sbarfile )
write(fnamesalt(1:80),'(2a,i10.10)')
& sbarfile(1:ilctdsclim),'.',optimcycle
endif
fcthread = 0. _d 0
#ifdef GENERIC_BAR_MONTH
c-- Loop over month
do irec = 1,12
nyears=int((nmonsrec-irec)/12)+1
if(nyears.gt.0) then
do iyear=1,nyears
mrec=irec+(iyear-1)*12
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnamesalt, sbar, mrec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_sbar_mean_dummy )
do bj = jtlo,jthi
do bi = itlo,ithi
do k = 1,nr
do j = jmin,jmax
do i = imin,imax
if(iyear.eq.1) then
sbar_gen(i,j,k,bi,bj) =sbar(i,j,k,bi,bj)
elseif(iyear.eq.nyears) then
sbar(i,j,k,bi,bj) =(sbar_gen(i,j,k,bi,bj)
$ +sbar(i,j,k,bi,bj))/float(nyears)
else
sbar_gen(i,j,k,bi,bj) =sbar_gen(i,j,k,bi,bj)
$ +sbar(i,j,k,bi,bj)
endif
enddo
enddo
enddo
enddo
enddo
enddo
#else
c-- Loop over records.
do irec = 1,nmonsrec
c-- Read time averages and the monthly mean data.
call ACTIVE_READ_XYZ( fnamesalt, sbar, irec,
& doglobalread, ladinit,
& optimcycle, mythid,
& xx_sbar_mean_dummy )
#endif
c-- Determine the month to be read.
levoff = mod(modelstartdate(1)/100,100)
levmon = (irec-1) + levoff
levmon = mod(levmon-1,12)+1
call MDSREADFIELD( ctdsclimfile, cost_iprec, cost_yftype,
& nr, sdat, levmon, mythid)
do bj = jtlo,jthi
do bi = itlo,ithi
c-- Loop over the model layers
fctile = 0. _d 0
do k = 1,nr
c-- Determine the mask or weights
do j = jmin,jmax
do i = imin,imax
cmask(i,j) = 1. _d 0
if (sdat(i,j,k,bi,bj) .eq. 0.) then
cmask(i,j) = 0. _d 0
endif
if (sdat(i,j,k,bi,bj) .lt. spval) 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) .eq. 0.) then
cmask(i,j) = 0. _d 0
endif
enddo
enddo
c-- Compute model data misfit and cost function term for
c the salinity field.
do j = jmin,jmax
do i = imin,imax
if (_hFacC(i,j,k,bi,bj) .ne. 0.) then
fctile = fctile +
& (wsalt2(i,j,k,bi,bj)*cosphi(i,j,bi,bj)*
& cmask(i,j)*
& (sbar(i,j,k,bi,bj) - sdat(i,j,k,bi,bj))*
& (sbar(i,j,k,bi,bj) - sdat(i,j,k,bi,bj)) )
endif
enddo
enddo
enddo
c-- End of loop over layers.
fcthread = fcthread + fctile
objf_ctdsclim(bi,bj) = objf_ctdsclim(bi,bj) + fctile
#ifdef ECCO_VERBOSE
c-- Print cost function for each tile in each thread.
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
& ' cost_Ctdsclim: irec,bi,bj = ',irec,bi,bj
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,d22.15)')
& ' cost function (salinity) = ',
& fctile
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 , myThid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,i8.8)')
& ' cost_Ctdsclim: irec = ',irec
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a,d22.15)')
& ' global cost function value',
& ' (salinity) = ',fcthread
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
#endif
#ifdef GENERIC_BAR_MONTH
endif
#endif
enddo
c-- End of loop over records.
#else
c-- Do not enter the calculation of the salinity contribution to
c-- the final cost function.
_BEGIN_MASTER( mythid )
write(msgbuf,'(a)') ' '
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a)')
& ' cost_Ctdsclim: no contribution of salinity field ',
& 'to cost function.'
call PRINT_MESSAGE( msgbuf, standardmessageunit,
& SQUEEZE_RIGHT , mythid)
write(msgbuf,'(a,a,i9.8)')
& ' cost_Ctdsclim: number of records that would have',
& ' been processed: ',nmonsrec
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