pkg/ecco/cost_obcse.F

8f7d13d0c9 Jean*0001 #include "ECCO_OPTIONS.h"
6805a315c1 Gael*0002 #ifdef ALLOW_CTRL
                0003 # include "CTRL_OPTIONS.h"
                0004 #endif
5001c65f45 Patr*0005
15f6a15ab5 Mart*0006 CBOP
                0007 C     !ROUTINE: COST_OBCSE
                0008 C     !INTERFACE:
5001c65f45 Patr*0009       subroutine cost_obcse(
720a211d38 Ou W*0010      I                       startrec, endrec,
                0011      I                       myTime, myIter, myThid )
5001c65f45 Patr*0012
15f6a15ab5 Mart*0013 C     !DESCRIPTION: \bv
5001c65f45 Patr*0014 c     ==================================================================
                0015 c     SUBROUTINE cost_obcse
                0016 c     ==================================================================
                0017 c
                0018 c     o cost function contribution obc
                0019 c
                0020 c     ==================================================================
                0021 c     SUBROUTINE cost_obcse
                0022 c     ==================================================================
15f6a15ab5 Mart*0023 C     \ev
                0024
                0025 C     !USES:
5001c65f45 Patr*0026
                0027       implicit none
                0028
                0029 c     == global variables ==
                0030 #include "EEPARAMS.h"
                0031 #include "SIZE.h"
                0032 #include "PARAMS.h"
                0033 #include "GRID.h"
f75c953e99 Jean*0034 c#ifdef ALLOW_OBCS
                0035 c# include "OBCS_GRID.h"
                0036 c#endif
5001c65f45 Patr*0037
c509d7e04a Gael*0038 #ifdef ALLOW_CAL
                0039 # include "cal.h"
                0040 #endif
                0041 #ifdef ALLOW_CTRL
                0042 # include "CTRL_SIZE.h"
4d72283393 Mart*0043 # include "CTRL.h"
edcd27be69 Mart*0044 # include "CTRL_DUMMY.h"
65754df434 Mart*0045 # include "OPTIMCYCLE.h"
c509d7e04a Gael*0046 # include "CTRL_OBCS.h"
                0047 #endif
5001c65f45 Patr*0048
15f6a15ab5 Mart*0049 C     !INPUT/OUTPUT PARAMETERS:
720a211d38 Ou W*0050       integer startrec, endrec
                0051       _RL     myTime
                0052       integer myIter
                0053       integer myThid
5001c65f45 Patr*0054
c509d7e04a Gael*0055 #if (defined (ALLOW_CTRL) && defined (ALLOW_OBCS))
                0056
f75c953e99 Jean*0057 #ifdef ALLOW_OBCSE_COST_CONTRIBUTION
720a211d38 Ou W*0058 C     !FUNCTIONS:
f75c953e99 Jean*0059       integer  ilnblnk
                0060       external ilnblnk
                0061
15f6a15ab5 Mart*0062 C     !LOCAL VARIABLES:
5001c65f45 Patr*0063       integer bi,bj
f75c953e99 Jean*0064       integer j,k
5001c65f45 Patr*0065       integer jmin,jmax
                0066       integer irec
                0067       integer iobcs
15f6a15ab5 Mart*0068       integer nrec
                0069       integer ilfld
                0070       integer igg
5001c65f45 Patr*0071       _RL fctile
                0072       _RL fcthread
                0073       _RL dummy
15f6a15ab5 Mart*0074       _RL gg
                0075       _RL tmpx
                0076 cgg(
9f5240b52a Jean*0077       _RL tmpfield (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0078       _RL maskyz   (1-OLy:sNy+OLy,Nr,nSx,nSy)
de57a2ec4b Mart*0079       character*(MAX_LEN_FNAM) fnamefld
5001c65f45 Patr*0080       logical doglobalread
                0081       logical ladinit
                0082 #ifdef ECCO_VERBOSE
                0083       character*(MAX_LEN_MBUF) msgbuf
                0084 #endif
15f6a15ab5 Mart*0085 CEOP
5001c65f45 Patr*0086
                0087       jmin = 1
9f5240b52a Jean*0088       jmax = sNy
5001c65f45 Patr*0089
                0090 c--   Read tiled data.
                0091       doglobalread = .false.
                0092       ladinit      = .false.
                0093
15f6a15ab5 Mart*0094 c     Number of records to be used.
                0095       nrec = endrec-startrec+1
                0096
f75c953e99 Jean*0097 c     ip1 = 1
5001c65f45 Patr*0098       fcthread = 0. _d 0
                0099
15f6a15ab5 Mart*0100 #ifdef ECCO_VERBOSE
720a211d38 Ou W*0101       _BEGIN_MASTER( myThid )
15f6a15ab5 Mart*0102       write(msgbuf,'(a)') ' '
720a211d38 Ou W*0103       call print_message( msgbuf, standardMessageUnit,
                0104      &                    SQUEEZE_RIGHT, myThid )
15f6a15ab5 Mart*0105       write(msgbuf,'(a)') ' '
720a211d38 Ou W*0106       call print_message( msgbuf, standardMessageUnit,
                0107      &                    SQUEEZE_RIGHT, myThid )
15f6a15ab5 Mart*0108       write(msgbuf,'(a,i9.8)')
                0109      &  ' cost_obcse: number of records to process: ',nrec
720a211d38 Ou W*0110       call print_message( msgbuf, standardMessageUnit,
                0111      &                    SQUEEZE_RIGHT, myThid )
15f6a15ab5 Mart*0112       write(msgbuf,'(a)') ' '
720a211d38 Ou W*0113       call print_message( msgbuf, standardMessageUnit,
                0114      &                    SQUEEZE_RIGHT, myThid )
                0115       _END_MASTER( myThid )
15f6a15ab5 Mart*0116 #endif
                0117
                0118       if (optimcycle .ge. 0) then
                0119         ilfld=ilnblnk( xx_obcse_file )
de57a2ec4b Mart*0120         write(fnamefld,'(2a,i10.10)')
15f6a15ab5 Mart*0121      &       xx_obcse_file(1:ilfld), '.', optimcycle
                0122       endif
                0123
5001c65f45 Patr*0124 c--   Loop over records.
15f6a15ab5 Mart*0125       do irec = 1,nrec
5001c65f45 Patr*0126
720a211d38 Ou W*0127 #ifdef ALLOW_AUTODIFF
                0128         call active_read_yz( fnamefld, tmpfield, irec,
                0129      &                       doglobalread, ladinit, optimcycle,
                0130      &                       myThid, xx_obcse_dummy )
                0131 #else
                0132         CALL READ_REC_YZ_RL( fnamefld, ctrlprec, Nr,
                0133      &                       tmpfield, irec, 1, myThid )
                0134 #endif
c9dc83bee0 Patr*0135
720a211d38 Ou W*0136 cgg     Need to solve for iobcs would have been.
                0137         gg    = (irec-1)/nobcs
                0138         igg   = int(gg)
                0139         iobcs = irec - igg*nobcs
                0140
                0141 #ifdef ALLOW_AUTODIFF
                0142         call active_read_yz( 'maskobcse', maskyz, iobcs,
                0143      &                       doglobalread, ladinit, 0,
                0144      &                       myThid, dummy )
                0145 #else
                0146         CALL READ_REC_YZ_RL( 'maskobcse', ctrlprec, Nr,
                0147      &                       maskyz, iobcs, 1, myThid )
                0148 #endif
5001c65f45 Patr*0149
80451941d6 Jean*0150 c--     Loop over this thread s tiles.
9f5240b52a Jean*0151         do bj = myByLo(myThid), myByHi(myThid)
                0152           do bi = myBxLo(myThid), myBxHi(myThid)
5001c65f45 Patr*0153
15f6a15ab5 Mart*0154 c--         Determine the weights to be used.
5001c65f45 Patr*0155             fctile = 0. _d 0
c9dc83bee0 Patr*0156
15f6a15ab5 Mart*0157             do k = 1, Nr
                0158               do j = jmin,jmax
f75c953e99 Jean*0159 c               i = OB_Iw(j,bi,bj)
15f6a15ab5 Mart*0160 cgg                if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then
                0161                   tmpx = tmpfield(j,k,bi,bj)
                0162 CMM                  fctile = fctile + wobcse2(j,k,bi,bj,iobcs)
                0163                   fctile = fctile + wobcse(k,iobcs)
                0164      &                        *tmpx*tmpx*maskyz(j,k,bi,bj)
                0165 cgg                endif
                0166 CMM                  if (wobcsw2(j,k,bi,bj,iobcs)*maskyz(j,k,bi,bj).ne.0.)
                0167                   if (wobcse(k,iobcs)*maskyz(j,k,bi,bj).ne.0.)
                0168      &                    num_obcse(bi,bj) = num_obcse(bi,bj) + 1. _d 0
                0169               enddo
c9dc83bee0 Patr*0170             enddo
15f6a15ab5 Mart*0171
5001c65f45 Patr*0172             objf_obcse(bi,bj) = objf_obcse(bi,bj) + fctile
15f6a15ab5 Mart*0173             fcthread         = fcthread + fctile
5001c65f45 Patr*0174           enddo
                0175         enddo
                0176
                0177 #ifdef ECCO_VERBOSE
                0178 c--     Print cost function for all tiles.
6637358eea Jean*0179         _GLOBAL_SUM_RL( fcthread , myThid )
5001c65f45 Patr*0180         write(msgbuf,'(a)') ' '
720a211d38 Ou W*0181         call print_message( msgbuf, standardMessageUnit,
                0182      &                      SQUEEZE_RIGHT, myThid )
5001c65f45 Patr*0183         write(msgbuf,'(a,i8.8)')
                0184      &    ' cost_obcse: irec = ',irec
720a211d38 Ou W*0185         call print_message( msgbuf, standardMessageUnit,
                0186      &                      SQUEEZE_RIGHT, myThid )
5001c65f45 Patr*0187         write(msgbuf,'(a,a,d22.15)')
                0188      &    ' global cost function value',
                0189      &    ' (obcse) = ',fcthread
720a211d38 Ou W*0190         call print_message( msgbuf, standardMessageUnit,
                0191      &                      SQUEEZE_RIGHT, myThid )
5001c65f45 Patr*0192         write(msgbuf,'(a)') ' '
720a211d38 Ou W*0193         call print_message( msgbuf, standardMessageUnit,
                0194      &                      SQUEEZE_RIGHT, myThid )
5001c65f45 Patr*0195 #endif
                0196
                0197       enddo
                0198 c--   End of loop over records.
                0199
f75c953e99 Jean*0200 #endif /* ALLOW_OBCSE_COST_CONTRIBUTION */
5001c65f45 Patr*0201
c509d7e04a Gael*0202 #endif /* ALLOW_CTRL and ALLOW_OBCS */
                0203
5001c65f45 Patr*0204       return
                0205       end