Back to home page

MITgcm
 
 

    

File indexing completed on 2024-03-19 05:10:16 UTC

view on githubraw file Latest commit 720a211d on 2024-03-18 17:10:23 UTC
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_obcs_ageos( myiter, mytime, mythid )
                0007 
                0008 c     ==================================================================
                0009 c     SUBROUTINE cost_obcs_ageos
                0010 c     ==================================================================
                0011 c
                0012 c     o cost function contribution obc -- Ageostrophic boundary flow.
                0013 c
                0014 c     started: G. Gebbie gebbie@mit.edu 4-Feb-2003
                0015 c
                0016 c     warning: masks redundantly assume no gradient of topography at
                0017 c              boundary.
                0018 c
                0019 c     ==================================================================
                0020 c     SUBROUTINE cost_obcs_ageos
                0021 c     ==================================================================
                0022 
                0023       implicit none
                0024 
                0025 c     == global variables ==
                0026 
                0027 #include "EEPARAMS.h"
                0028 #include "SIZE.h"
                0029 #include "GRID.h"
                0030 #include "DYNVARS.h"
                0031 #include "PARAMS.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
                0039 #ifdef ALLOW_ECCO
                0040 # include "ecco_cost.h"
                0041 #endif
                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 # include "CTRL_OBCS.h"
                0048 #endif
5001c65f45 Patr*0049 
                0050 c     == routine arguments ==
                0051 
                0052       integer myiter
                0053       _RL     mytime
                0054       integer mythid
                0055 
c509d7e04a Gael*0056 #if (defined (ALLOW_CTRL) && \
                0057      defined (ALLOW_OBCS) && \
                0058      defined (ALLOW_ECCO))
                0059 
9f5240b52a Jean*0060 #ifdef OBCS_AGEOS_COST_CONTRIBUTION
5001c65f45 Patr*0061 c     == local variables ==
                0062 
                0063       integer bi,bj
                0064       integer i,j,k
                0065       integer itlo,ithi
                0066       integer jtlo,jthi
                0067       integer jmin,jmax
                0068       integer imin,imax
                0069       integer irec
                0070       integer levmon
                0071       integer levoff
                0072       integer iltheta
                0073       integer ilsalt
                0074       integer iluvel
                0075       integer ilvvel
                0076       integer ip1, jp1
                0077 
                0078       _RL fctile
                0079       _RL fcthread
                0080 
9f5240b52a Jean*0081       _RL rholoc (1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy)
                0082       _RL xzgrdrho(1-OLx:sNx+OLx,Nr,nSx,nSy)
                0083       _RL yzgrdrho(1-OLy:sNy+OLy,Nr,nSx,nSy)
                0084       _RL xzdvel1   (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0085       _RL xzdvel2   (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0086       _RL yzdvel1   (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0087       _RL yzdvel2   (1-OLy:sNy+OLy,Nr,nSx,nSy)
                0088       _RL maskxzageos   (1-OLx:sNx+OLx,Nr,nSx,nSy)
                0089       _RL maskyzageos   (1-OLy:sNy+OLy,Nr,nSx,nSy)
5001c65f45 Patr*0090       _RL dummy
                0091 
de57a2ec4b Mart*0092       character*(MAX_LEN_FNAM) fnametheta
                0093       character*(MAX_LEN_FNAM) fnamesalt
                0094       character*(MAX_LEN_FNAM) fnameuvel
                0095       character*(MAX_LEN_FNAM) fnamevvel
5001c65f45 Patr*0096 
                0097       logical doglobalread
                0098       logical ladinit
                0099 
                0100       character*(MAX_LEN_MBUF) msgbuf
                0101 
                0102 c     == external functions ==
                0103 
                0104       integer  ilnblnk
                0105       external ilnblnk
                0106 
                0107 c     == end of interface ==
                0108 
9f5240b52a Jean*0109       jtlo = myByLo(mythid)
                0110       jthi = myByHi(mythid)
                0111       itlo = myBxLo(mythid)
                0112       ithi = myBxHi(mythid)
5001c65f45 Patr*0113       jmin = 1
9f5240b52a Jean*0114       jmax = sNy
5001c65f45 Patr*0115       imin = 1
9f5240b52a Jean*0116       imax = sNx
5001c65f45 Patr*0117 
                0118 c--   Read tiled data.
                0119       doglobalread = .false.
                0120       ladinit      = .false.
                0121 
                0122 #ifdef ECCO_VERBOSE
                0123       _BEGIN_MASTER( mythid )
                0124         write(msgbuf,'(a)') ' '
                0125         call print_message( msgbuf, standardmessageunit,
                0126      &                      SQUEEZE_RIGHT , mythid)
                0127         write(msgbuf,'(a,i8.8)')
                0128      &    ' cost_obcs_ageos: number of records to process = ',nmonsrec
                0129         call print_message( msgbuf, standardmessageunit,
                0130      &                      SQUEEZE_RIGHT , mythid)
                0131         write(msgbuf,'(a)') ' '
                0132         call print_message( msgbuf, standardmessageunit,
                0133      &                      SQUEEZE_RIGHT , mythid)
                0134       _END_MASTER( mythid )
9f5240b52a Jean*0135 #endif /* ECCO_VERBOSE */
5001c65f45 Patr*0136 
                0137       if (optimcycle .ge. 0) then
                0138         iltheta = ilnblnk( tbarfile )
de57a2ec4b Mart*0139         write(fnametheta,'(2a,i10.10)')
5001c65f45 Patr*0140      &    tbarfile(1:iltheta),'.',optimcycle
                0141         ilsalt = ilnblnk( sbarfile )
de57a2ec4b Mart*0142         write(fnamesalt,'(2a,i10.10)')
5001c65f45 Patr*0143      &    sbarfile(1:ilsalt),'.',optimcycle
                0144         iluvel = ilnblnk( ubarfile )
de57a2ec4b Mart*0145         write(fnameuvel,'(2a,i10.10)')
5001c65f45 Patr*0146      &    ubarfile(1:iluvel),'.',optimcycle
                0147         ilvvel = ilnblnk( vbarfile )
de57a2ec4b Mart*0148         write(fnamevvel,'(2a,i10.10)')
5001c65f45 Patr*0149      &    vbarfile(1:ilvvel),'.',optimcycle
                0150       endif
                0151 
                0152       fcthread = 0. _d 0
                0153       fctile   = 0. _d 0
                0154 
                0155 cgg   Code safe: always initialize to zero.
                0156       do bj = jtlo,jthi
                0157         do bi = itlo,ithi
9f5240b52a Jean*0158           do k = 1,Nr
                0159             do i = 1-OLx,sNx+OLx
5001c65f45 Patr*0160               maskxzageos(i,k,bi,bj) = 0. _d 0
                0161               xzdvel1(i,k,bi,bj)     = 0. _d 0
                0162               xzdvel2(i,k,bi,bj)     = 0. _d 0
                0163               xzgrdrho(i,k,bi,bj)    = 0. _d 0
                0164             enddo
9f5240b52a Jean*0165             do j = 1-OLy,sNy+OLy
5001c65f45 Patr*0166               maskyzageos(j,k,bi,bj) = 0. _d 0
                0167               yzdvel1(j,k,bi,bj)     = 0. _d 0
                0168               yzdvel2(j,k,bi,bj)     = 0. _d 0
                0169               yzgrdrho(j,k,bi,bj)    = 0. _d 0
                0170             enddo
                0171           enddo
                0172         enddo
                0173       enddo
                0174 
                0175       do bj = jtlo,jthi
                0176         do bi = itlo,ithi
9f5240b52a Jean*0177           do j = 1-OLy,sNy+OLy
                0178             do i = 1-OLx,sNx+OLx
5001c65f45 Patr*0179               rholoc(i,j,bi,bj) = 0. _d 0
                0180             enddo
                0181           enddo
                0182         enddo
                0183       enddo
                0184 
                0185 c--   Loop over records.
                0186       do irec = 1,nmonsrec
                0187 
                0188 c--     Read time averages and the monthly mean data.
720a211d38 Ou W*0189 #ifdef ALLOW_AUTODIFF
5001c65f45 Patr*0190         call active_read_xyz( fnametheta, tbar, irec,
                0191      &                        doglobalread, ladinit,
                0192      &                        optimcycle, mythid,
                0193      &                        xx_tbar_mean_dummy )
                0194 
                0195         call active_read_xyz( fnamesalt, sbar, irec,
                0196      &                        doglobalread, ladinit,
                0197      &                        optimcycle, mythid,
                0198      &                        xx_sbar_mean_dummy )
                0199 
                0200         call active_read_xyz( fnameuvel, ubar, irec,
                0201      &                        doglobalread, ladinit,
                0202      &                        optimcycle, mythid,
                0203      &                        xx_ubar_mean_dummy )
                0204 
                0205         call active_read_xyz( fnamevvel, vbar, irec,
                0206      &                        doglobalread, ladinit,
                0207      &                        optimcycle, mythid,
                0208      &                        xx_vbar_mean_dummy )
720a211d38 Ou W*0209 #else
                0210         call read_rec_xyz_rl( fnametheta, tbar, irec, 1, myThid )
                0211 
                0212         call read_rec_xyz_rl( fnamesalt, sbar, irec, 1, myThid )
                0213 
                0214         call read_rec_xyz_rl( fnameuvel, ubar, irec, 1, myThid )
                0215 
                0216         call read_rec_xyz_rl( fnamevvel, vbar, irec, 1, myThid )
                0217 #endif
5001c65f45 Patr*0218 
                0219 cgg     Minor problem : grad T,S needs overlap values.
                0220         _BARRIER
6637358eea Jean*0221         _EXCH_XYZ_RL(tbar,myThid)
                0222         _EXCH_XYZ_RL(sbar,myThid)
5001c65f45 Patr*0223 
                0224         do bj = jtlo,jthi
                0225           do bi = itlo,ithi
                0226 
                0227 #ifdef ALLOW_OBCSN_CONTROL
                0228             jp1 = 0
                0229 
                0230 cgg    Make a mask for the velocity shear comparison.
9f5240b52a Jean*0231             do k = 1,Nr-1
5001c65f45 Patr*0232               do i = imin, imax
9fd29e65a3 Jean*0233                 j = OB_Jn(i,bi,bj)
951926fb9b Jean*0234 cgg    All these points need to be wet.
9fd29e65a3 Jean*0235                 if ( j.eq.OB_indexNone ) then
5001c65f45 Patr*0236                   maskxzageos(i,k,bi,bj) = 0.
                0237                 else
951926fb9b Jean*0238                   maskxzageos(i,k,bi,bj) =
5001c65f45 Patr*0239      &       hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
                0240      &       hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
                0241      &       hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
                0242      &       hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
                0243                 endif
                0244               enddo
                0245             enddo
951926fb9b Jean*0246 
9f5240b52a Jean*0247             do k = 1,Nr
852c7dfbb9 Jean*0248 
                0249 C-- jmc: both calls below are wrong if more than 1 tile => stop here
46c1d12550 Jean*0250        IF ( bi.NE.1 .OR. bj.NE.1 )
852c7dfbb9 Jean*0251      &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
                0252               call find_rho_2d(
                0253      I                         iMin, iMax, jMin, jMax, k,
46c1d12550 Jean*0254      I                         tbar(1-OLx,1-OLy,k,bi,bj),
852c7dfbb9 Jean*0255      I                         sbar(1-OLx,1-OLy,k,bi,bj),
                0256      O                         rholoc,
                0257      I                         k, bi, bj, myThid )
6637358eea Jean*0258               _EXCH_XY_RL(rholoc , myThid)
5001c65f45 Patr*0259 
                0260 cgg   Compute centered difference horizontal gradient on bdy.
                0261               do i = imin, imax
9fd29e65a3 Jean*0262                 j = OB_Jn(i,bi,bj)
                0263                 if ( j.eq.OB_indexNone ) j = 1
951926fb9b Jean*0264                   xzgrdrho(i,k,bi,bj) =
5001c65f45 Patr*0265      &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
                0266      &            /(2.*dxc(i,j+jp1,bi,bj))
                0267               enddo
                0268             enddo
                0269 
                0270 cgg         Compute vertical shear from geostrophy/thermal wind.
                0271 cgg         Above level 4 needs not be geostrophic.
                0272 cgg         Please get rid of the "4" ASAP. Ridiculous.
                0273             do k = 4,Nr-1
                0274               do i = imin,imax
9fd29e65a3 Jean*0275                 j = OB_Jn(i,bi,bj)
                0276                 if ( j.eq.OB_indexNone ) j = 1
5001c65f45 Patr*0277                   xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)
                0278      &                               - vbar(i,j+jp1,k+1,bi,bj)
                0279                  xzdvel2(i,k,bi,bj)=((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
                0280      &                             (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
                0281      &                        * gravity / f0 / rhonil
                0282 
                0283                   fctile = fctile + 100*wcurrent(k,bi,bj) *
                0284      &            maskxzageos(i,k,bi,bj)*
                0285      &          (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
                0286      &          (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
                0287       if (maskxzageos(i,k,bi,bj) .ne. 0) then
                0288 cgg      print*,'XX shear,geos shear',xzdvel1(i,k,bi,bj),xzdvel2(i,k,bi,bj)
                0289 cgg      print*,'i,j,k,fctile N',i,j,k,fctile
                0290       endif
                0291               enddo
                0292             enddo
                0293 c--         End of loop over layers.
                0294 #endif /* ALLOW_OBCSN_CONTROL */
                0295 
                0296 #ifdef ALLOW_OBCSS_CONTROL
                0297             jp1 = 1
                0298 
                0299 cgg    Make a mask for the velocity shear comparison.
9f5240b52a Jean*0300             do k = 1,Nr-1
5001c65f45 Patr*0301               do i = imin, imax
9fd29e65a3 Jean*0302                 j = OB_Js(i,bi,bj)
                0303                 if ( j.eq.OB_indexNone ) then
5001c65f45 Patr*0304                   maskxzageos(i,k,bi,bj) = 0.
                0305                 else
951926fb9b Jean*0306 cgg    All these points need to be wet.
                0307                   maskxzageos(i,k,bi,bj) =
5001c65f45 Patr*0308      &       hfacC(i,j+jp1,k,bi,bj)*hfacC(i+1,j+jp1,k,bi,bj) *
                0309      &       hfacC(i-1,j+jp1,k,bi,bj)*hfacC(i,j+jp1,k+1,bi,bj)*
                0310      &       hfacC(i-1,j+jp1,k+1,bi,bj)*hfacC(i+1,j+jp1,k+1,bi,bj)*
                0311      &       hfacS(i,j+jp1,k,bi,bj)*hfacS(i,j+jp1,k+1,bi,bj)
951926fb9b Jean*0312                 endif
5001c65f45 Patr*0313               enddo
                0314             enddo
                0315 
9f5240b52a Jean*0316             do k = 1,Nr
5001c65f45 Patr*0317 
46c1d12550 Jean*0318 C-- jmc: both calls below are wrong if more than 1 tile => stop here
                0319        IF ( bi.NE.1 .OR. bj.NE.1 )
                0320      &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
852c7dfbb9 Jean*0321               call find_rho_2d(
                0322      I                         iMin, iMax, jMin, jMax, k,
46c1d12550 Jean*0323      I                         tbar(1-OLx,1-OLy,k,bi,bj),
852c7dfbb9 Jean*0324      I                         sbar(1-OLx,1-OLy,k,bi,bj),
                0325      O                         rholoc,
                0326      I                         k, bi, bj, myThid )
                0327 
6637358eea Jean*0328               _EXCH_XY_RL(rholoc , myThid)
5001c65f45 Patr*0329 
                0330 cgg   Compute centered difference horizontal gradient on bdy.
                0331                do i = imin, imax
9fd29e65a3 Jean*0332                  j = OB_Js(i,bi,bj)
                0333                  if ( j.eq.OB_indexNone ) j = 1
951926fb9b Jean*0334                  xzgrdrho(i,k,bi,bj) =
5001c65f45 Patr*0335      &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
                0336      &            /(2.*dxc(i,j+jp1,bi,bj))
                0337                enddo
                0338              enddo
                0339 
                0340 cgg         Compute vertical shear from geostrophy/thermal wind.
                0341              do k = 4,Nr-1
                0342                do i = imin,imax
9fd29e65a3 Jean*0343                  j = OB_Js(i,bi,bj)
                0344                  if ( j.eq.OB_indexNone ) j = 1
5001c65f45 Patr*0345 cgg         Retrieve the model vertical shear.
                0346                  xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)
                0347      &                               - vbar(i,j+jp1,k+1,bi,bj)
                0348 
                0349 cgg         Compute vertical shear from geostrophy/thermal wind.
                0350                  xzdvel2(i,k,bi,bj) =((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
                0351      &                             (xzgrdrho(i,k+1,bi,bj)*delz(k+1)/2.))
                0352      &                        * gravity / f0 /rhonil
                0353 
                0354 cgg   Make a comparison.
                0355                   fctile = fctile + 100*wcurrent(k,bi,bj) *
                0356      &          maskxzageos(i,k,bi,bj)*
                0357      &          (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))*
                0358      &          (xzdvel2(i,k,bi,bj) - xzdvel1(i,k,bi,bj))
                0359 cgg        print*,'fctile S',fctile
                0360               enddo
                0361             enddo
                0362 c--         End of loop over layers.
9f5240b52a Jean*0363 #endif /* ALLOW_OBCSS_CONTROL */
5001c65f45 Patr*0364 
                0365 #ifdef ALLOW_OBCSW_CONTROL
                0366             ip1 = 1
                0367 
                0368 cgg    Make a mask for the velocity shear comparison.
9f5240b52a Jean*0369             do k = 1,Nr-1
5001c65f45 Patr*0370               do j = jmin, jmax
9fd29e65a3 Jean*0371                 i = OB_Iw(j,bi,bj)
5001c65f45 Patr*0372 cgg    All these points need to be wet.
9fd29e65a3 Jean*0373                 if ( i.eq.OB_indexNone ) then
5001c65f45 Patr*0374                   maskyzageos(j,k,bi,bj) = 0.
951926fb9b Jean*0375                 else
                0376                   maskyzageos(j,k,bi,bj) =
5001c65f45 Patr*0377      &       hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
                0378      &       hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
                0379      &       hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
                0380      &       hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
                0381                 endif
                0382               enddo
                0383             enddo
                0384 
9f5240b52a Jean*0385             do k = 1,Nr
5001c65f45 Patr*0386 
46c1d12550 Jean*0387        IF ( bi.NE.1 .OR. bj.NE.1 )
                0388      &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
852c7dfbb9 Jean*0389               call find_rho_2d(
                0390      I                         iMin, iMax, jMin, jMax, k,
46c1d12550 Jean*0391      I                         tbar(1-OLx,1-OLy,k,bi,bj),
852c7dfbb9 Jean*0392      I                         sbar(1-OLx,1-OLy,k,bi,bj),
                0393      O                         rholoc,
                0394      I                         k, bi, bj, myThid )
6637358eea Jean*0395               _EXCH_XY_RL(rholoc , myThid)
5001c65f45 Patr*0396 
                0397 cgg   Compute centered difference horizontal gradient on bdy.
                0398               do j = jmin, jmax
9fd29e65a3 Jean*0399                 i = OB_Iw(j,bi,bj)
                0400                 if ( i.eq.OB_indexNone ) i = 1
5001c65f45 Patr*0401 cgg             Negative sign due to geostrophy.
951926fb9b Jean*0402                 yzgrdrho(j,k,bi,bj) =
5001c65f45 Patr*0403      &            (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
                0404      &            /(2.*dyc(i+ip1,j,bi,bj))
                0405               enddo
                0406             enddo
                0407 
                0408 cgg         Compute vertical shear from geostrophy/thermal wind.
                0409             do k = 4,Nr-1
                0410               do j = jmin,jmax
9fd29e65a3 Jean*0411                 i = OB_Iw(j,bi,bj)
                0412                 if ( i.eq.OB_indexNone ) i = 1
5001c65f45 Patr*0413 cgg         Retrieve the model vertical shear.
                0414                 yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)
                0415      &                               - ubar(i+ip1,j,k+1,bi,bj)
                0416 
                0417 cgg         Compute vertical shear from geostrophy/thermal wind.
                0418                 yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k  ,bi,bj)*delz(k)/2.)+
                0419      &                             (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
                0420      &                        * gravity / f0 / rhonil
                0421 
                0422 cgg   Make a comparison.
                0423                 fctile = fctile + 100*wcurrent(k,bi,bj) *
                0424      &           maskyzageos(j,k,bi,bj) *
                0425      &          (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
                0426      &          (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
                0427               enddo
                0428             enddo
                0429 c--         End of loop over layers.
9f5240b52a Jean*0430 #endif /* ALLOW_OBCSW_CONTROL */
5001c65f45 Patr*0431 
                0432 #ifdef ALLOW_OBCSE_CONTROL
                0433             ip1 = 0
                0434 
                0435 cgg    Make a mask for the velocity shear comparison.
9f5240b52a Jean*0436             do k = 1,Nr-1
5001c65f45 Patr*0437               do j = jmin, jmax
9fd29e65a3 Jean*0438                 i = OB_Ie(j,bi,bj)
                0439                 if ( i.eq.OB_indexNone ) then
5001c65f45 Patr*0440                   maskyzageos(j,k,bi,bj) =0.
                0441                 else
951926fb9b Jean*0442 cgg    All these points need to be wet.
                0443                   maskyzageos(j,k,bi,bj) =
5001c65f45 Patr*0444      &       hfacC(i+ip1,j,k,bi,bj)*hfacC(i+ip1,j+1,k,bi,bj) *
                0445      &       hfacC(i+ip1,j-1,k,bi,bj)*hfacC(i+ip1,j,k+1,bi,bj)*
                0446      &       hfacC(i+ip1,j-1,k+1,bi,bj)*hfacC(i+ip1,j+1,k+1,bi,bj)*
                0447      &       hfacW(i+ip1,j,k,bi,bj)*hfacW(i+ip1,j,k+1,bi,bj)
                0448                 endif
                0449               enddo
                0450             enddo
                0451 
9f5240b52a Jean*0452             do k = 1,Nr
5001c65f45 Patr*0453 
46c1d12550 Jean*0454        IF ( bi.NE.1 .OR. bj.NE.1 )
                0455      &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
852c7dfbb9 Jean*0456               call find_rho_2d(
                0457      I                         iMin, iMax, jMin, jMax, k,
46c1d12550 Jean*0458      I                         tbar(1-OLx,1-OLy,k,bi,bj),
852c7dfbb9 Jean*0459      I                         sbar(1-OLx,1-OLy,k,bi,bj),
                0460      O                         rholoc,
                0461      I                         k, bi, bj, myThid )
6637358eea Jean*0462               _EXCH_XY_RL(rholoc , myThid)
5001c65f45 Patr*0463 
                0464 cgg   Compute centered difference horizontal gradient on bdy.
                0465               do j = jmin, jmax
9fd29e65a3 Jean*0466                 i = OB_Ie(j,bi,bj)
                0467                 if ( i.eq.OB_indexNone ) i = 1
5001c65f45 Patr*0468 cgg             Negative sign due to geostrophy.
951926fb9b Jean*0469                 yzgrdrho(j,k,bi,bj) =
720a211d38 Ou W*0470      &            (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
5001c65f45 Patr*0471      &            /(2.*dyc(i+ip1,j,bi,bj))
                0472               enddo
                0473             enddo
                0474 
                0475 cgg         Compute vertical shear from geostrophy/thermal wind.
                0476             do k = 4,Nr-1
                0477               do j = jmin,jmax
9fd29e65a3 Jean*0478                 i = OB_Ie(j,bi,bj)
                0479                 if ( i.eq.OB_indexNone ) i = 1
5001c65f45 Patr*0480 cgg         Retrieve the model vertical shear.
                0481                 yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)
                0482      &                             - ubar(i+ip1,j,k+1,bi,bj)
                0483 
                0484 cgg         Compute vertical shear from geostrophy/thermal wind.
                0485                 yzdvel2(j,k,bi,bj) =((yzgrdrho(j,k  ,bi,bj)*delz(k)/2.)+
                0486      &                             (yzgrdrho(j,k+1,bi,bj)*delz(k+1)/2.))
                0487      &                        * gravity / f0 /rhonil
                0488 
                0489 cgg   Make a comparison.
                0490                 fctile = fctile + 100*wcurrent(k,bi,bj) *
                0491      &           maskyzageos(j,k,bi,bj) *
                0492      &          (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))*
                0493      &          (yzdvel2(j,k,bi,bj) - yzdvel1(j,k,bi,bj))
                0494 
                0495               enddo
                0496             enddo
                0497 c--         End of loop over layers.
9f5240b52a Jean*0498 #endif /* ALLOW_OBCSE_CONTROL */
5001c65f45 Patr*0499 
                0500             fcthread          = fcthread          + fctile
                0501             objf_ageos(bi,bj) = objf_ageos(bi,bj) + fctile
                0502 
                0503 #ifdef ECCO_VERBOSE
                0504 c--         Print cost function for each tile in each thread.
                0505             write(msgbuf,'(a)') ' '
                0506             call print_message( msgbuf, standardmessageunit,
                0507      &                          SQUEEZE_RIGHT , mythid)
                0508             write(msgbuf,'(a,i8.8,1x,i3.3,1x,i3.3)')
                0509      &        ' cost_Theta: irec,bi,bj          =  ',irec,bi,bj
                0510             call print_message( msgbuf, standardmessageunit,
                0511      &                          SQUEEZE_RIGHT , mythid)
                0512             write(msgbuf,'(a,d22.15)')
                0513      &        '     cost function (temperature) = ',
                0514      &        fctile
                0515             call print_message( msgbuf, standardmessageunit,
                0516      &                      SQUEEZE_RIGHT , mythid)
                0517             write(msgbuf,'(a)') ' '
                0518             call print_message( msgbuf, standardmessageunit,
                0519      &                          SQUEEZE_RIGHT , mythid)
9f5240b52a Jean*0520 #endif /* ECCO_VERBOSE */
5001c65f45 Patr*0521 
                0522           enddo
                0523         enddo
                0524 
                0525 #ifdef ECCO_VERBOSE
                0526 c--     Print cost function for all tiles.
6637358eea Jean*0527         _GLOBAL_SUM_RL( fcthread , myThid )
5001c65f45 Patr*0528         write(msgbuf,'(a)') ' '
                0529         call print_message( msgbuf, standardmessageunit,
                0530      &                      SQUEEZE_RIGHT , mythid)
                0531         write(msgbuf,'(a,i8.8)')
                0532      &    ' cost_Theta: irec = ',irec
                0533         call print_message( msgbuf, standardmessageunit,
                0534      &                      SQUEEZE_RIGHT , mythid)
                0535         write(msgbuf,'(a,a,d22.15)')
                0536      &    ' global cost function value',
                0537      &    ' (temperature) = ',fcthread
                0538         call print_message( msgbuf, standardmessageunit,
                0539      &                      SQUEEZE_RIGHT , mythid)
                0540         write(msgbuf,'(a)') ' '
                0541         call print_message( msgbuf, standardmessageunit,
                0542      &                      SQUEEZE_RIGHT , mythid)
9f5240b52a Jean*0543 #endif /* ECCO_VERBOSE */
5001c65f45 Patr*0544 
                0545       enddo
                0546 c--   End of loop over records.
                0547 
9f5240b52a Jean*0548 #endif /* OBCS_AGEOS_COST_CONTRIBUTION */
5001c65f45 Patr*0549 
c509d7e04a Gael*0550 #endif /* ALLOW_CTRL, ALLOW_OBCS and ALLOW_ECCO */
                0551 
5001c65f45 Patr*0552       return
                0553       end
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!