Back to home page

MITgcm
 
 

    

File indexing completed on 2026-03-19 05:08:48 UTC

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