pkg/obcs/obcs_cost_vol.F

69361556c2 Mart*0001 #include "OBCS_OPTIONS.h"
6805a315c1 Gael*0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
69361556c2 Mart*0005 #ifdef ALLOW_ECCO
                0006 # include "ECCO_OPTIONS.h"
                0007 #endif
5001c65f45 Patr*0008
69361556c2 Mart*0009       subroutine obcs_cost_vol(
                0010      I                          startrec, endrec,
                0011      I                          mytime, myiter, mythid )
5001c65f45 Patr*0012
                0013 c     ==================================================================
69361556c2 Mart*0014 c     SUBROUTINE obcs_cost_vol
5001c65f45 Patr*0015 c     ==================================================================
                0016 c
                0017 c     o cost function contribution obc -- Volume flux imbalance.
                0018 c
                0019 c     ==================================================================
69361556c2 Mart*0020 c     SUBROUTINE obcs_cost_vol
5001c65f45 Patr*0021 c     ==================================================================
                0022
                0023       implicit none
                0024
                0025 c     == global variables ==
                0026
69361556c2 Mart*0027 #ifdef OBCS_VOLFLUX_COST_CONTRIBUTION
5001c65f45 Patr*0028 #include "EEPARAMS.h"
                0029 #include "SIZE.h"
                0030 #include "PARAMS.h"
                0031 #include "GRID.h"
                0032 #ifdef ALLOW_OBCS
46c1d12550 Jean*0033 # include "OBCS_GRID.h"
5001c65f45 Patr*0034 #endif
                0035
c509d7e04a Gael*0036 #ifdef ALLOW_CAL
                0037 # include "cal.h"
                0038 #endif
69361556c2 Mart*0039 c#ifdef ALLOW_ECCO
                0040 c# include "ecco_cost.h"
                0041 c#endif
c509d7e04a Gael*0042 #ifdef ALLOW_CTRL
                0043 # include "CTRL_SIZE.h"
4d72283393 Mart*0044 # include "CTRL.h"
edcd27be69 Mart*0045 # include "CTRL_DUMMY.h"
65754df434 Mart*0046 # include "OPTIMCYCLE.h"
c509d7e04a Gael*0047 #endif
69361556c2 Mart*0048 #endif
5001c65f45 Patr*0049
                0050 c     == routine arguments ==
720a211d38 Ou W*0051       integer startrec, endrec
5001c65f45 Patr*0052       _RL     mytime
720a211d38 Ou W*0053       integer myiter
9f5240b52a Jean*0054       integer mythid
5001c65f45 Patr*0055
c509d7e04a Gael*0056 #if (defined (ALLOW_CTRL) && defined (ALLOW_OBCS))
                0057
9f5240b52a Jean*0058 #ifdef OBCS_VOLFLUX_COST_CONTRIBUTION
                0059 #ifdef BAROTROPIC_OBVEL_CONTROL
                0060 c     == external functions ==
                0061       integer  ilnblnk
                0062       external ilnblnk
5001c65f45 Patr*0063
9f5240b52a Jean*0064 c     == local variables ==
5001c65f45 Patr*0065       integer bi,bj
                0066       integer i,j,k
                0067       integer itlo,ithi
                0068       integer jtlo,jthi
                0069       integer jmin,jmax
                0070       integer imin,imax
                0071       integer irec
                0072       integer iobcs
                0073       integer nrec
                0074       integer ilfld
                0075       integer igg
                0076       _RL fctile
                0077       _RL sumvol
                0078       _RL gg
                0079       _RL tmpx
                0080       _RL tmpy
                0081       _RL wobcsvol
de57a2ec4b Mart*0082       character*(MAX_LEN_FNAM) fnamefldn
                0083       character*(MAX_LEN_FNAM) fnameflds
                0084       character*(MAX_LEN_FNAM) fnamefldw
                0085       character*(MAX_LEN_FNAM) fnameflde
e7f7a6f2e9 Mart*0086 #if (defined ALLOW_OBCSN_CONTROL || defined ALLOW_OBCSS_CONTROL)
9f5240b52a Jean*0087       _RL tmpfldxz (1-OLx:sNx+OLx,Nr,nSx,nSy)
46c1d12550 Jean*0088 #endif
e7f7a6f2e9 Mart*0089 #if (defined ALLOW_OBCSE_CONTROL || defined ALLOW_OBCSW_CONTROL)
9f5240b52a Jean*0090       _RL tmpfldyz (1-OLy:sNy+OLy,Nr,nSx,nSy)
e7f7a6f2e9 Mart*0091 #endif
5001c65f45 Patr*0092       logical doglobalread
                0093       logical ladinit
                0094 #ifdef ECCO_VERBOSE
                0095       character*(MAX_LEN_MBUF) msgbuf
                0096 #endif
                0097 c     == end of interface ==
                0098
69361556c2 Mart*0099       stop 's/r obcs_cost_vol needs to be fixed'
747605eccb Mart*0100
9f5240b52a Jean*0101       jtlo = myByLo(mythid)
                0102       jthi = myByHi(mythid)
                0103       itlo = myBxLo(mythid)
                0104       ithi = myBxHi(mythid)
5001c65f45 Patr*0105       jmin = 1
9f5240b52a Jean*0106       jmax = sNy
5001c65f45 Patr*0107       imin = 1
9f5240b52a Jean*0108       imax = sNx
5001c65f45 Patr*0109
                0110 c--   Read tiled data.
                0111       doglobalread = .false.
                0112       ladinit      = .false.
                0113
                0114 cgg   Assume the number of records is the same for
                0115 cgg   all boundaries.
                0116 c     Number of records to be used.
                0117       nrec = endrec-startrec+1
                0118
                0119       sumvol = 0. _d 0
                0120       wobcsvol = .01
                0121 cgg   Acceptable volume flux is 10^-3. Corresponds to 5 mm change over a year.
951926fb9b Jean*0122 cgg   Added a factor of 1000 because its very important to me.
5001c65f45 Patr*0123       wobcsvol = 1./(wobcsvol * wobcsvol)
                0124
                0125 #ifdef ECCO_VERBOSE
                0126       _BEGIN_MASTER( mythid )
                0127       write(msgbuf,'(a)') ' '
                0128       call print_message( msgbuf, standardmessageunit,
                0129      &                    SQUEEZE_RIGHT , mythid)
                0130       write(msgbuf,'(a)') ' '
                0131       call print_message( msgbuf, standardmessageunit,
                0132      &                    SQUEEZE_RIGHT , mythid)
                0133       write(msgbuf,'(a,i9.8)')
69361556c2 Mart*0134      &  ' obcs_cost_vol: number of records to process: ',nrec
5001c65f45 Patr*0135       call print_message( msgbuf, standardmessageunit,
                0136      &                    SQUEEZE_RIGHT , mythid)
                0137       write(msgbuf,'(a)') ' '
                0138       call print_message( msgbuf, standardmessageunit,
                0139      &                    SQUEEZE_RIGHT , mythid)
                0140       _END_MASTER( mythid )
                0141 #endif
                0142
                0143       if (optimcycle .ge. 0) then
                0144 #ifdef ALLOW_OBCSN_CONTROL
                0145         ilfld=ilnblnk( xx_obcsn_file )
de57a2ec4b Mart*0146         write(fnamefldn,'(2a,i10.10)')
5001c65f45 Patr*0147      &        xx_obcsn_file(1:ilfld),'.', optimcycle
                0148 #endif
                0149 #ifdef ALLOW_OBCSS_CONTROL
                0150         ilfld=ilnblnk( xx_obcss_file )
de57a2ec4b Mart*0151         write(fnameflds,'(2a,i10.10)')
5001c65f45 Patr*0152      &        xx_obcss_file(1:ilfld),'.',optimcycle
                0153 #endif
                0154 #ifdef ALLOW_OBCSW_CONTROL
                0155         ilfld=ilnblnk( xx_obcsw_file )
de57a2ec4b Mart*0156         write(fnamefldw,'(2a,i10.10)')
5001c65f45 Patr*0157      &        xx_obcsw_file(1:ilfld),'.',optimcycle
                0158 #endif
                0159 #ifdef ALLOW_OBCSE_CONTROL
                0160         ilfld=ilnblnk( xx_obcse_file )
de57a2ec4b Mart*0161         write(fnameflde,'(2a,i10.10)')
5001c65f45 Patr*0162      &        xx_obcse_file(1:ilfld),'.',optimcycle
                0163 #endif
                0164       else
                0165          print*
ce53c8659f Mart*0166          print*,' obcs_obcsvol: optimcycle has a wrong value.'
5001c65f45 Patr*0167          print*,'                 optimcycle = ',optimcycle
                0168          print*
ce53c8659f Mart*0169          stop   '  ... stopped in obcs_obcsvol.'
5001c65f45 Patr*0170       endif
                0171
                0172       do irec = 1,nrec
                0173 c--   Loop over records. For north boundary, we only need V velocity.
                0174
                0175 cgg    Need to solve for iobcs. Then only keep iobcs=3.or.4.
                0176         gg   = (irec-1)/nobcs
                0177         igg  = int(gg)
                0178         iobcs = irec - igg*nobcs
                0179
                0180 #ifdef ALLOW_OBCSN_CONTROL
2146dab1aa Jean*0181 cgg   Assume that nobcs=4, and V velocity is the 4th record. I cannot
5001c65f45 Patr*0182 cgg   think of a more general way to do this.
                0183         if (iobcs.eq.4) then
720a211d38 Ou W*0184 #ifdef ALLOW_AUTODIFF
5001c65f45 Patr*0185           call active_read_xz( fnamefldn, tmpfldxz, irec, doglobalread,
720a211d38 Ou W*0186      &                         ladinit, optimcycle, mythid,
                0187      &                         xx_obcsn_dummy )
                0188 #else
                0189           CALL READ_REC_XZ_RL( fnamefldn, ctrlprec, Nr,
                0190      &                         tmpfldxz, irec, 1, mythid )
                0191 #endif
5001c65f45 Patr*0192
                0193 cgg  At this point, do not be concerned with the overlap halos.
951926fb9b Jean*0194 cgg  From experience, there is no control contribution in the
5001c65f45 Patr*0195 cgg  velocity points outside the boundaries. This has something
                0196 cgg  to do with the computational stencil, and the fact that we
                0197 cgg  are diagonally offset. Could check later by employing both
                0198 cgg  BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION.
                0199 cgg
                0200 cgg  25-jan-03 --- no idea what i was talking about ^^^^
2146dab1aa Jean*0201 c--     Loop over this thread tiles.
5001c65f45 Patr*0202           do bj = jtlo,jthi
                0203             do bi = itlo,ithi
                0204
                0205 c--         Determine the weights to be used.
                0206               fctile = 0. _d 0
                0207
                0208               do k = 1, Nr
                0209                 do i = imin,imax
                0210                   j = OB_Jn(I,bi,bj)
9fd29e65a3 Jean*0211                   IF ( j.EQ.OB_indexNone ) j = 1
5001c65f45 Patr*0212 cgg    Barotropic velocity is stored in level 1.
                0213                   tmpx = tmpfldxz(i,1,bi,bj)
e7f7a6f2e9 Mart*0214                   if (maskS(i,j,k,bi,bj) .ne. 0.) then
5001c65f45 Patr*0215 cgg -- Positive is flux in.
e7f7a6f2e9 Mart*0216                     fctile = fctile - tmpx* drF(k) *dxg(i,j,bi,bj)
                0217      &                  * _hFacS(i,j,k,bi,bj)
5001c65f45 Patr*0218                   endif
                0219                 enddo
                0220               enddo
                0221
                0222               sumvol         = sumvol + fctile
                0223             enddo
                0224           enddo
                0225         endif
                0226 #endif
                0227
                0228 #ifdef ALLOW_OBCSS_CONTROL
2146dab1aa Jean*0229 cgg   Assume that nobcs=4, and V velocity is the 4th record. I cannot
5001c65f45 Patr*0230 cgg   think of a more general way to do this.
                0231         if (iobcs.eq.4) then
720a211d38 Ou W*0232 #ifdef ALLOW_AUTODIFF
5001c65f45 Patr*0233           call active_read_xz( fnameflds, tmpfldxz, irec, doglobalread,
720a211d38 Ou W*0234      &                         ladinit, optimcycle, mythid,
                0235      &                         xx_obcss_dummy )
                0236 #else
                0237           CALL READ_REC_XZ_RL( fnameflds, ctrlprec, Nr,
                0238      &                         tmpfldxz, irec, 1, mythid )
                0239 #endif
5001c65f45 Patr*0240
                0241 cgg  At this point, do not be concerned with the overlap halos.
951926fb9b Jean*0242 cgg  From experience, there is no control contribution in the
5001c65f45 Patr*0243 cgg  velocity points outside the boundaries. This has something
                0244 cgg  to do with the computational stencil, and the fact that we
                0245 cgg  are diagonally offset. Could check later by employing both
                0246 cgg  BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION.
                0247
2146dab1aa Jean*0248 c--     Loop over this thread tiles.
5001c65f45 Patr*0249           do bj = jtlo,jthi
                0250             do bi = itlo,ithi
                0251
                0252 c--         Determine the weights to be used.
                0253               fctile = 0. _d 0
                0254
                0255               do k = 1, Nr
                0256                 do i = imin,imax
                0257                   j = OB_Js(I,bi,bj)
9fd29e65a3 Jean*0258                   IF ( j.EQ.OB_indexNone ) j = 1
5001c65f45 Patr*0259 cgg    Barotropic velocity is stored in level 1.
                0260                   tmpx = tmpfldxz(i,1,bi,bj)
e7f7a6f2e9 Mart*0261                   if (maskS(i,j+1,k,bi,bj) .ne. 0.) then
5001c65f45 Patr*0262 cgg -- Positive is flux in.
e7f7a6f2e9 Mart*0263                     fctile = fctile + tmpx* drF(k) *dxg(i,j+1,bi,bj)
                0264      &                  * _hFacS(i,j+1,k,bi,bj)
5001c65f45 Patr*0265                   endif
                0266                 enddo
                0267               enddo
                0268
                0269               sumvol         = sumvol + fctile
                0270             enddo
                0271           enddo
                0272         endif
                0273
                0274 #endif
                0275
                0276 #ifdef ALLOW_OBCSW_CONTROL
2146dab1aa Jean*0277 cgg   Assume that nobcs=4, and V velocity is the 4th record. I cannot
5001c65f45 Patr*0278 cgg   think of a more general way to do this.
                0279         if (iobcs.eq.3) then
720a211d38 Ou W*0280 #ifdef ALLOW_AUTODIFF
5001c65f45 Patr*0281           call active_read_yz( fnamefldw, tmpfldyz, irec, doglobalread,
720a211d38 Ou W*0282      &                         ladinit, optimcycle, mythid,
                0283      &                         xx_obcsw_dummy )
                0284 #else
                0285           CALL READ_REC_YZ_RL( fnamefldw, ctrlprec, Nr,
                0286      &                         tmpfldyz, irec, 1, mythid )
                0287 #endif
5001c65f45 Patr*0288
                0289 cgg  At this point, do not be concerned with the overlap halos.
951926fb9b Jean*0290 cgg  From experience, there is no control contribution in the
5001c65f45 Patr*0291 cgg  velocity points outside the boundaries. This has something
                0292 cgg  to do with the computational stencil, and the fact that we
                0293 cgg  are diagonally offset. Could check later by employing both
                0294 cgg  BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION.
                0295
2146dab1aa Jean*0296 c--     Loop over this thread tiles.
5001c65f45 Patr*0297           do bj = jtlo,jthi
                0298             do bi = itlo,ithi
                0299
                0300 c--         Determine the weights to be used.
                0301               fctile = 0. _d 0
                0302
                0303               do k = 1, Nr
                0304                 do j = jmin,jmax
                0305                   i = OB_Iw(j,bi,bj)
9fd29e65a3 Jean*0306                   IF ( i.EQ.OB_indexNone ) i = 1
5001c65f45 Patr*0307 cgg    Barotropic velocity is stored in the level 1.
                0308                   tmpy = tmpfldyz(j,1,bi,bj)
e7f7a6f2e9 Mart*0309                   if (maskW(i+1,j,k,bi,bj) .ne. 0.) then
5001c65f45 Patr*0310 cgg -- Positive is flux in.
e7f7a6f2e9 Mart*0311                     fctile = fctile + tmpy* drF(k) *dyg(i+1,j,bi,bj)
                0312      &                  * _hFacW(i+1,j,k,bi,bj)
5001c65f45 Patr*0313                   endif
                0314                 enddo
                0315               enddo
                0316
                0317               sumvol         = sumvol + fctile
                0318             enddo
                0319           enddo
                0320         endif
                0321
                0322 #endif
                0323
                0324 #ifdef ALLOW_OBCSE_CONTROL
2146dab1aa Jean*0325 cgg   Assume that nobcs=4, and V velocity is the 4th record. I cannot
5001c65f45 Patr*0326 cgg   think of a more general way to do this.
                0327         if (iobcs.eq.3) then
720a211d38 Ou W*0328 #ifdef ALLOW_AUTODIFF
5001c65f45 Patr*0329           call active_read_yz( fnameflde, tmpfldyz, irec, doglobalread,
720a211d38 Ou W*0330      &                         ladinit, optimcycle, mythid,
                0331      &                         xx_obcse_dummy )
                0332 #else
                0333           CALL READ_REC_YZ_RL( fnameflde, ctrlprec, Nr,
                0334      &                         tmpfldyz, irec, 1, mythid )
                0335 #endif
5001c65f45 Patr*0336
                0337 cgg  At this point, do not be concerned with the overlap halos.
951926fb9b Jean*0338 cgg  From experience, there is no control contribution in the
5001c65f45 Patr*0339 cgg  velocity points outside the boundaries. This has something
                0340 cgg  to do with the computational stencil, and the fact that we
                0341 cgg  are diagonally offset. Could check later by employing both
                0342 cgg  BALANCE_CONTROL_VOLFLUX and VOLFLUX_COST_CONTRIBUTION.
                0343
2146dab1aa Jean*0344 c--     Loop over this thread tiles.
5001c65f45 Patr*0345           do bj = jtlo,jthi
                0346             do bi = itlo,ithi
                0347
                0348 c--         Determine the weights to be used.
                0349               fctile = 0. _d 0
                0350
                0351               do k = 1, Nr
                0352                 do j = jmin,jmax
                0353                   i = OB_Ie(j,bi,bj)
9fd29e65a3 Jean*0354                   IF ( i.EQ.OB_indexNone ) i = 1
5001c65f45 Patr*0355 cgg    Barotropic velocity stored in level 1.
                0356                   tmpy = tmpfldyz(j,1,bi,bj)
e7f7a6f2e9 Mart*0357                   if (maskW(i,j,k,bi,bj) .ne. 0.) then
5001c65f45 Patr*0358 cgg -- Positive is flux in.
e7f7a6f2e9 Mart*0359                     fctile = fctile - tmpy* drF(k) *dyg(i,j,bi,bj)
                0360      &                  * _hFacW(i,j,k,bi,bj)
5001c65f45 Patr*0361                   endif
                0362                 enddo
                0363               enddo
                0364
                0365               sumvol         = sumvol + fctile
                0366             enddo
                0367           enddo
                0368         endif
                0369
                0370 #endif
                0371
                0372       enddo
                0373 c--   End of loop over records.
                0374
                0375 c--   Do the global summation.
6637358eea Jean*0376       _GLOBAL_SUM_RL( sumvol, mythid )
951926fb9b Jean*0377       objf_obcsvol =  wobcsvol * sumvol* sumvol
5001c65f45 Patr*0378
9f5240b52a Jean*0379 #endif /* BAROTROPIC_OBVEL_CONTROL */
                0380 #endif /* OBCS_VOLFLUX_COST_CONTRIBUTION */
5001c65f45 Patr*0381
c509d7e04a Gael*0382 #endif /* ALLOW_CTRL and ALLOW_OBCS */
                0383
5001c65f45 Patr*0384       return
                0385       end