C $Header: /u/gcmpack/MITgcm/pkg/ecco/cost_obcsvol.F,v 1.13 2014/10/20 03:16:12 gforget Exp $ C $Name: $ #include "ECCO_OPTIONS.h" #ifdef ALLOW_CTRL # include "CTRL_OPTIONS.h" #endif subroutine COST_OBCSVOL( I myiter, I mytime, I startrec, I endrec, I mythid & ) c ================================================================== c SUBROUTINE cost_obcsvol c ================================================================== c c o cost function contribution obc -- Volume flux imbalance. c c ================================================================== c SUBROUTINE cost_obcsvol c ================================================================== implicit none c == global variables == #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" #ifdef ALLOW_OBCS # include "OBCS_GRID.h" #endif #ifdef ALLOW_CAL # include "cal.h" #endif #ifdef ALLOW_ECCO # include "ecco_cost.h" #endif #ifdef ALLOW_CTRL # include "CTRL_SIZE.h" # include "ctrl.h" # include "ctrl_dummy.h" # include "optim.h" #endif c == routine arguments == integer myiter _RL mytime integer mythid integer startrec integer endrec #if (defined (ALLOW_CTRL) defined (ALLOW_OBCS)) 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 iobcs integer nrec integer ilfld integer igg _RL fctile _RL sumvol _RL gg _RL tmpx _RL tmpy _RL wobcsvol character*(80) fnamefldn character*(80) fnameflds character*(80) fnamefldw character*(80) fnameflde #if (defined ALLOW_OBCSN_CONTROL defined ALLOW_OBCSS_CONTROL) _RL tmpfldxz (1-olx:snx+olx,nr,nsx,nsy) #endif #if (defined ALLOW_OBCSE_CONTROL defined ALLOW_OBCSW_CONTROL) _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy) #endif logical doglobalread logical ladinit #ifdef ECCO_VERBOSE character*(MAX_LEN_MBUF) msgbuf #endif c == external functions == integer ilnblnk external c == end of interface == #ifdef OBCS_VOLFLUX_COST_CONTRIBUTION #ifdef BAROTROPIC_OBVEL_CONTROL stop 's/r cost_obcsvol needs to be fixed' jtlo = mybylo(mythid) jthi = mybyhi(mythid) itlo = mybxlo(mythid) ithi = mybxhi(mythid) jmin = 1 jmax = sny imin = 1 imax = snx c-- Read tiled data. doglobalread = .false. ladinit = .false. cgg Assume the number of records is the same for cgg all boundaries. c Number of records to be used. nrec = endrec-startrec+1 sumvol = 0. _d 0 wobcsvol = .01 cgg Acceptable volume flux is 10^-3. Corresponds to 5 mm change over a year. cgg Added a factor of 1000 because its very important to me. wobcsvol = 1./(wobcsvol * wobcsvol) #ifdef ECCO_VERBOSE _BEGIN_MASTER( mythid ) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a,i9.8)') & ' cost_obcsvol: number of records to process: ',nrec call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) write(msgbuf,'(a)') ' ' call PRINT_MESSAGE( msgbuf, standardmessageunit, & SQUEEZE_RIGHT , mythid) _END_MASTER( mythid ) #endif if (optimcycle .ge. 0) then #ifdef ALLOW_OBCSN_CONTROL ilfld=ilnblnk( xx_obcsn_file ) write(fnamefldn(1:80),'(2a,i10.10)') & xx_obcsn_file(1:ilfld),'.', optimcycle #endif #ifdef ALLOW_OBCSS_CONTROL ilfld=ilnblnk( xx_obcss_file ) write(fnameflds(1:80),'(2a,i10.10)') & xx_obcss_file(1:ilfld),'.',optimcycle #endif #ifdef ALLOW_OBCSW_CONTROL ilfld=ilnblnk( xx_obcsw_file ) write(fnamefldw(1:80),'(2a,i10.10)') & xx_obcsw_file(1:ilfld),'.',optimcycle #endif #ifdef ALLOW_OBCSE_CONTROL ilfld=ilnblnk( xx_obcse_file ) write(fnameflde(1:80),'(2a,i10.10)') & xx_obcse_file(1:ilfld),'.',optimcycle #endif else print* print*,' obcs_obcsvol: optimcycle has a wrong value.' print*,' optimcycle = ',optimcycle print* stop ' ... stopped in obcs_obcsvol.' endif do irec = 1,nrec c-- Loop over records. For north boundary, we only need V velocity. cgg Need to solve for iobcs. Then only keep iobcs=3.or.4. gg = (irec-1)/nobcs igg = int(gg) iobcs = irec - igg*nobcs #ifdef ALLOW_OBCSN_CONTROL cgg Assume that nobcs=4, and V velocity is the 4th record. I cannot cgg think of a more general way to do this. if (iobcs.eq.4) then call ACTIVE_READ_XZ( fnamefldn, tmpfldxz, irec, doglobalread, & ladinit, optimcycle, mythid & , xx_obcsn_dummy ) cgg At this point, do not be concerned with the overlap halos. cgg From experience, there is no control contribution in the cgg velocity points outside the boundaries. This has something cgg to do with the computational stencil, and the fact that we cgg are diagonally offset. Could check later by employing both cgg BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION. cgg cgg 25-jan-03 --- no idea what i was talking about ^^^^ c-- Loop over this thread tiles. do bj = jtlo,jthi do bi = itlo,ithi c-- Determine the weights to be used. fctile = 0. _d 0 do k = 1, Nr do i = imin,imax j = OB_Jn(I,bi,bj) IF ( j.EQ.OB_indexNone ) j = 1 cgg Barotropic velocity is stored in level 1. tmpx = tmpfldxz(i,1,bi,bj) if (maskS(i,j,k,bi,bj) .ne. 0.) then cgg -- Positive is flux in. fctile = fctile - tmpx* drF(k) *dxg(i,j,bi,bj) & * _hFacS(i,j,k,bi,bj) endif enddo enddo sumvol = sumvol + fctile enddo enddo endif #endif #ifdef ALLOW_OBCSS_CONTROL cgg Assume that nobcs=4, and V velocity is the 4th record. I cannot cgg think of a more general way to do this. if (iobcs.eq.4) then call ACTIVE_READ_XZ( fnameflds, tmpfldxz, irec, doglobalread, & ladinit, optimcycle, mythid & , xx_obcss_dummy ) cgg At this point, do not be concerned with the overlap halos. cgg From experience, there is no control contribution in the cgg velocity points outside the boundaries. This has something cgg to do with the computational stencil, and the fact that we cgg are diagonally offset. Could check later by employing both cgg BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION. c-- Loop over this thread tiles. do bj = jtlo,jthi do bi = itlo,ithi c-- Determine the weights to be used. fctile = 0. _d 0 do k = 1, Nr do i = imin,imax j = OB_Js(I,bi,bj) IF ( j.EQ.OB_indexNone ) j = 1 cgg Barotropic velocity is stored in level 1. tmpx = tmpfldxz(i,1,bi,bj) if (maskS(i,j+1,k,bi,bj) .ne. 0.) then cgg -- Positive is flux in. fctile = fctile + tmpx* drF(k) *dxg(i,j+1,bi,bj) & * _hFacS(i,j+1,k,bi,bj) endif enddo enddo sumvol = sumvol + fctile enddo enddo endif #endif #ifdef ALLOW_OBCSW_CONTROL cgg Assume that nobcs=4, and V velocity is the 4th record. I cannot cgg think of a more general way to do this. if (iobcs.eq.3) then call ACTIVE_READ_YZ( fnamefldw, tmpfldyz, irec, doglobalread, & ladinit, optimcycle, mythid & , xx_obcsw_dummy ) cgg At this point, do not be concerned with the overlap halos. cgg From experience, there is no control contribution in the cgg velocity points outside the boundaries. This has something cgg to do with the computational stencil, and the fact that we cgg are diagonally offset. Could check later by employing both cgg BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION. c-- Loop over this thread tiles. do bj = jtlo,jthi do bi = itlo,ithi c-- Determine the weights to be used. fctile = 0. _d 0 do k = 1, Nr do j = jmin,jmax i = OB_Iw(j,bi,bj) IF ( i.EQ.OB_indexNone ) i = 1 cgg Barotropic velocity is stored in the level 1. tmpy = tmpfldyz(j,1,bi,bj) if (maskW(i+1,j,k,bi,bj) .ne. 0.) then cgg -- Positive is flux in. fctile = fctile + tmpy* drF(k) *dyg(i+1,j,bi,bj) & * _hFacW(i+1,j,k,bi,bj) endif enddo enddo sumvol = sumvol + fctile enddo enddo endif #endif #ifdef ALLOW_OBCSE_CONTROL cgg Assume that nobcs=4, and V velocity is the 4th record. I cannot cgg think of a more general way to do this. if (iobcs.eq.3) then call ACTIVE_READ_YZ( fnameflde, tmpfldyz, irec, doglobalread, & ladinit, optimcycle, mythid & , xx_obcse_dummy ) cgg At this point, do not be concerned with the overlap halos. cgg From experience, there is no control contribution in the cgg velocity points outside the boundaries. This has something cgg to do with the computational stencil, and the fact that we cgg are diagonally offset. Could check later by employing both cgg BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION. c-- Loop over this thread tiles. do bj = jtlo,jthi do bi = itlo,ithi c-- Determine the weights to be used. fctile = 0. _d 0 do k = 1, Nr do j = jmin,jmax i = OB_Ie(j,bi,bj) IF ( i.EQ.OB_indexNone ) i = 1 cgg Barotropic velocity stored in level 1. tmpy = tmpfldyz(j,1,bi,bj) if (maskW(i,j,k,bi,bj) .ne. 0.) then cgg -- Positive is flux in. fctile = fctile - tmpy* drF(k) *dyg(i,j,bi,bj) & * _hFacW(i,j,k,bi,bj) endif enddo enddo sumvol = sumvol + fctile enddo enddo endif #endif enddo c-- End of loop over records. c-- Do the global summation. _GLOBAL_SUM_RL( sumvol, mythid ) objf_obcsvol = wobcsvol * sumvol* sumvol #endif #endif #endif /* ALLOW_CTRL and ALLOW_OBCS */ return end