pkg/ecco/cost_obcsvol.F

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