/[MITgcm]/MITgcm/pkg/ecco/cost_obcs_ageos.F
ViewVC logotype

Diff of /MITgcm/pkg/ecco/cost_obcs_ageos.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.2 by jmc, Tue Oct 9 00:02:50 2007 UTC revision 1.8 by jmc, Tue Sep 18 20:16:34 2012 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "COST_CPPOPTIONS.h"  #include "ECCO_OPTIONS.h"
 #ifdef ALLOW_OBCS  
 # include "OBCS_OPTIONS.h"  
 #endif  
5    
6        subroutine cost_obcs_ageos( myiter, mytime, mythid )        subroutine cost_obcs_ageos( myiter, mytime, mythid )
7    
# Line 33  c     == global variables == Line 30  c     == global variables ==
30  #include "DYNVARS.h"  #include "DYNVARS.h"
31  #include "PARAMS.h"  #include "PARAMS.h"
32  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
33  # include "OBCS.h"  # include "OBCS_GRID.h"
34  #endif  #endif
35    
36  #include "cal.h"  #include "cal.h"
37  #include "ecco_cost.h"  #include "ecco_cost.h"
38    #include "CTRL_SIZE.h"
39  #include "ctrl.h"  #include "ctrl.h"
40  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
41  #include "optim.h"  #include "optim.h"
# Line 50  c     == routine arguments == Line 48  c     == routine arguments ==
48    
49  c     == local variables ==  c     == local variables ==
50    
       _RS        one_rs  
       parameter( one_rs = 1. )  
   
51        integer bi,bj        integer bi,bj
52        integer i,j,k        integer i,j,k
53        integer itlo,ithi        integer itlo,ithi
# Line 206  c--     Read time averages and the month Line 201  c--     Read time averages and the month
201    
202  cgg     Minor problem : grad T,S needs overlap values.  cgg     Minor problem : grad T,S needs overlap values.
203          _BARRIER          _BARRIER
204          _EXCH_XYZ_R8(tbar,myThid)          _EXCH_XYZ_RL(tbar,myThid)
205          _EXCH_XYZ_R8(sbar,myThid)          _EXCH_XYZ_RL(sbar,myThid)
206    
207    
208          do bj = jtlo,jthi          do bj = jtlo,jthi
# Line 219  cgg     Minor problem : grad T,S needs o Line 214  cgg     Minor problem : grad T,S needs o
214  cgg    Make a mask for the velocity shear comparison.  cgg    Make a mask for the velocity shear comparison.
215              do k = 1,nr-1              do k = 1,nr-1
216                do i = imin, imax                do i = imin, imax
217                  j = ob_jn(i,bi,bj)                  j = OB_Jn(i,bi,bj)
218  cgg    All these points need to be wet.  cgg    All these points need to be wet.
219                  if (j .eq. 0) then                  if ( j.eq.OB_indexNone ) then
220                    maskxzageos(i,k,bi,bj) = 0.                    maskxzageos(i,k,bi,bj) = 0.
221                  else                  else
222                    maskxzageos(i,k,bi,bj) =                    maskxzageos(i,k,bi,bj) =
# Line 234  cgg    All these points need to be wet. Line 229  cgg    All these points need to be wet.
229              enddo              enddo
230    
231              do k = 1,nr              do k = 1,nr
232                call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,  
233       &                      tbar,sbar,rholoc,mythid)  C-- jmc: both calls below are wrong if more than 1 tile => stop here
234                _EXCH_XY_R8(rholoc , myThid)         IF ( bi.NE.1 .OR. bj.NE.1 )
235         &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
236                  call find_rho_2d(
237         I                         iMin, iMax, jMin, jMax, k,
238         I                         tbar(1-OLx,1-OLy,k,bi,bj),
239         I                         sbar(1-OLx,1-OLy,k,bi,bj),
240         O                         rholoc,
241         I                         k, bi, bj, myThid )
242                  _EXCH_XY_RL(rholoc , myThid)
243    
244  cgg   Compute centered difference horizontal gradient on bdy.  cgg   Compute centered difference horizontal gradient on bdy.
245                do i = imin, imax                do i = imin, imax
246                  j = ob_jn(i,bi,bj)                  j = OB_Jn(i,bi,bj)
247                    if ( j.eq.OB_indexNone ) j = 1
248                    xzgrdrho(i,k,bi,bj) =                    xzgrdrho(i,k,bi,bj) =
249       &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))       &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
250       &            /(2.*dxc(i,j+jp1,bi,bj))       &            /(2.*dxc(i,j+jp1,bi,bj))
# Line 252  cgg         Above level 4 needs not be g Line 256  cgg         Above level 4 needs not be g
256  cgg         Please get rid of the "4" ASAP. Ridiculous.  cgg         Please get rid of the "4" ASAP. Ridiculous.
257              do k = 4,Nr-1              do k = 4,Nr-1
258                do i = imin,imax                do i = imin,imax
259                  j = ob_jn(i,bi,bj)                  j = OB_Jn(i,bi,bj)
260                    if ( j.eq.OB_indexNone ) j = 1
261                    xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)                    xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)
262       &                               - vbar(i,j+jp1,k+1,bi,bj)       &                               - vbar(i,j+jp1,k+1,bi,bj)
263                   xzdvel2(i,k,bi,bj)=((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+                   xzdvel2(i,k,bi,bj)=((xzgrdrho(i,k,bi,bj)*delz(k)/2.)+
# Line 278  c--         End of loop over layers. Line 283  c--         End of loop over layers.
283  cgg    Make a mask for the velocity shear comparison.  cgg    Make a mask for the velocity shear comparison.
284              do k = 1,nr-1              do k = 1,nr-1
285                do i = imin, imax                do i = imin, imax
286                  j = ob_js(i,bi,bj)                  j = OB_Js(i,bi,bj)
287                  if (j .eq. 0) then                  if ( j.eq.OB_indexNone ) then
288                    maskxzageos(i,k,bi,bj) = 0.                    maskxzageos(i,k,bi,bj) = 0.
289                  else                  else
290  cgg    All these points need to be wet.  cgg    All these points need to be wet.
# Line 294  cgg    All these points need to be wet. Line 299  cgg    All these points need to be wet.
299    
300              do k = 1,nr              do k = 1,nr
301    
302                call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,  C-- jmc: both calls below are wrong if more than 1 tile => stop here
303       &                      tbar,sbar,rholoc,mythid)         IF ( bi.NE.1 .OR. bj.NE.1 )
304                _EXCH_XY_R8(rholoc , myThid)       &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
305                  call find_rho_2d(
306         I                         iMin, iMax, jMin, jMax, k,
307         I                         tbar(1-OLx,1-OLy,k,bi,bj),
308         I                         sbar(1-OLx,1-OLy,k,bi,bj),
309         O                         rholoc,
310         I                         k, bi, bj, myThid )
311    
312                  _EXCH_XY_RL(rholoc , myThid)
313    
314  cgg   Compute centered difference horizontal gradient on bdy.  cgg   Compute centered difference horizontal gradient on bdy.
315                 do i = imin, imax                 do i = imin, imax
316                   j = ob_js(i,bi,bj)                   j = OB_Js(i,bi,bj)
317                     if ( j.eq.OB_indexNone ) j = 1
318                   xzgrdrho(i,k,bi,bj) =                   xzgrdrho(i,k,bi,bj) =
319       &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))       &            (rholoc(i-1,j+jp1,bi,bj)-rholoc(i+1,j+jp1,bi,bj))
320       &            /(2.*dxc(i,j+jp1,bi,bj))       &            /(2.*dxc(i,j+jp1,bi,bj))
# Line 310  cgg   Compute centered difference horizo Line 324  cgg   Compute centered difference horizo
324  cgg         Compute vertical shear from geostrophy/thermal wind.  cgg         Compute vertical shear from geostrophy/thermal wind.
325               do k = 4,Nr-1               do k = 4,Nr-1
326                 do i = imin,imax                 do i = imin,imax
327                   j = ob_js(i,bi,bj)                   j = OB_Js(i,bi,bj)
328                     if ( j.eq.OB_indexNone ) j = 1
329  cgg         Retrieve the model vertical shear.  cgg         Retrieve the model vertical shear.
330                   xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)                   xzdvel1(i,k,bi,bj) = vbar(i,j+jp1,k  ,bi,bj)
331       &                               - vbar(i,j+jp1,k+1,bi,bj)       &                               - vbar(i,j+jp1,k+1,bi,bj)
# Line 337  c--         End of loop over layers. Line 352  c--         End of loop over layers.
352  cgg    Make a mask for the velocity shear comparison.  cgg    Make a mask for the velocity shear comparison.
353              do k = 1,nr-1              do k = 1,nr-1
354                do j = jmin, jmax                do j = jmin, jmax
355                  i = ob_iw(j,bi,bj)                  i = OB_Iw(j,bi,bj)
356  cgg    All these points need to be wet.  cgg    All these points need to be wet.
357                  if (i.eq. 0) then                  if ( i.eq.OB_indexNone ) then
358                    maskyzageos(j,k,bi,bj) = 0.                    maskyzageos(j,k,bi,bj) = 0.
359                  else                  else
360                    maskyzageos(j,k,bi,bj) =                    maskyzageos(j,k,bi,bj) =
# Line 353  cgg    All these points need to be wet. Line 368  cgg    All these points need to be wet.
368    
369              do k = 1,nr              do k = 1,nr
370    
371                call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,         IF ( bi.NE.1 .OR. bj.NE.1 )
372       &                      tbar,sbar,rholoc,mythid)       &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
373                _EXCH_XY_R8(rholoc , myThid)                call find_rho_2d(
374         I                         iMin, iMax, jMin, jMax, k,
375         I                         tbar(1-OLx,1-OLy,k,bi,bj),
376         I                         sbar(1-OLx,1-OLy,k,bi,bj),
377         O                         rholoc,
378         I                         k, bi, bj, myThid )
379                  _EXCH_XY_RL(rholoc , myThid)
380    
381  cgg   Compute centered difference horizontal gradient on bdy.  cgg   Compute centered difference horizontal gradient on bdy.
382                do j = jmin, jmax                do j = jmin, jmax
383                  i = ob_iw(j,bi,bj)                  i = OB_Iw(j,bi,bj)
384                    if ( i.eq.OB_indexNone ) i = 1
385  cgg             Negative sign due to geostrophy.  cgg             Negative sign due to geostrophy.
386                  yzgrdrho(j,k,bi,bj) =                  yzgrdrho(j,k,bi,bj) =
387       &            (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))       &            (rholoc(i+ip1,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
# Line 370  cgg             Negative sign due to geo Line 392  cgg             Negative sign due to geo
392  cgg         Compute vertical shear from geostrophy/thermal wind.  cgg         Compute vertical shear from geostrophy/thermal wind.
393              do k = 4,Nr-1              do k = 4,Nr-1
394                do j = jmin,jmax                do j = jmin,jmax
395                  i = ob_iw(j,bi,bj)                  i = OB_Iw(j,bi,bj)
396                    if ( i.eq.OB_indexNone ) i = 1
397  cgg         Retrieve the model vertical shear.  cgg         Retrieve the model vertical shear.
398                  yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)                  yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)
399       &                               - ubar(i+ip1,j,k+1,bi,bj)       &                               - ubar(i+ip1,j,k+1,bi,bj)
# Line 396  c--         End of loop over layers. Line 419  c--         End of loop over layers.
419  cgg    Make a mask for the velocity shear comparison.  cgg    Make a mask for the velocity shear comparison.
420              do k = 1,nr-1              do k = 1,nr-1
421                do j = jmin, jmax                do j = jmin, jmax
422                  i = ob_ie(j,bi,bj)                  i = OB_Ie(j,bi,bj)
423                  if (i.eq.0) then                  if ( i.eq.OB_indexNone ) then
424                    maskyzageos(j,k,bi,bj) =0.                    maskyzageos(j,k,bi,bj) =0.
425                  else                  else
426  cgg    All these points need to be wet.  cgg    All these points need to be wet.
# Line 412  cgg    All these points need to be wet. Line 435  cgg    All these points need to be wet.
435    
436              do k = 1,nr              do k = 1,nr
437    
438                call find_rho(bi,bj,imin,imax,jmin,jmax,k,k,         IF ( bi.NE.1 .OR. bj.NE.1 )
439       &                      tbar,sbar,rholoc,mythid)       &  STOP 'COST_OBCS_AGEOS wrong with more than 1 tile/proc'
440                _EXCH_XY_R8(rholoc , myThid)                call find_rho_2d(
441         I                         iMin, iMax, jMin, jMax, k,
442         I                         tbar(1-OLx,1-OLy,k,bi,bj),
443         I                         sbar(1-OLx,1-OLy,k,bi,bj),
444         O                         rholoc,
445         I                         k, bi, bj, myThid )
446                  _EXCH_XY_RL(rholoc , myThid)
447    
448  cgg   Compute centered difference horizontal gradient on bdy.  cgg   Compute centered difference horizontal gradient on bdy.
449                do j = jmin, jmax                do j = jmin, jmax
450                  i = ob_ie(j,bi,bj)                  i = OB_Ie(j,bi,bj)
451                    if ( i.eq.OB_indexNone ) i = 1
452  cgg             Negative sign due to geostrophy.  cgg             Negative sign due to geostrophy.
453                  yzgrdrho(j,k,bi,bj) =                  yzgrdrho(j,k,bi,bj) =
454       &            (rholoc(i+ip1,,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))       &            (rholoc(i+ip1,,j+1,bi,bj)-rholoc(i+ip1,j-1,bi,bj))
# Line 429  cgg             Negative sign due to geo Line 459  cgg             Negative sign due to geo
459  cgg         Compute vertical shear from geostrophy/thermal wind.  cgg         Compute vertical shear from geostrophy/thermal wind.
460              do k = 4,Nr-1              do k = 4,Nr-1
461                do j = jmin,jmax                do j = jmin,jmax
462                  i = ob_ie(j,bi,bj)                  i = OB_Ie(j,bi,bj)
463                    if ( i.eq.OB_indexNone ) i = 1
464  cgg         Retrieve the model vertical shear.  cgg         Retrieve the model vertical shear.
465                  yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)                  yzdvel1(j,k,bi,bj) = ubar(i+ip1,j,k  ,bi,bj)
466       &                             - ubar(i+ip1,j,k+1,bi,bj)       &                             - ubar(i+ip1,j,k+1,bi,bj)
# Line 477  c--         Print cost function for each Line 508  c--         Print cost function for each
508    
509  #ifdef ECCO_VERBOSE  #ifdef ECCO_VERBOSE
510  c--     Print cost function for all tiles.  c--     Print cost function for all tiles.
511          _GLOBAL_SUM_R8( fcthread , myThid )          _GLOBAL_SUM_RL( fcthread , myThid )
512          write(msgbuf,'(a)') ' '          write(msgbuf,'(a)') ' '
513          call print_message( msgbuf, standardmessageunit,          call print_message( msgbuf, standardmessageunit,
514       &                      SQUEEZE_RIGHT , mythid)       &                      SQUEEZE_RIGHT , mythid)
# Line 502  c--   End of loop over records. Line 533  c--   End of loop over records.
533    
534        return        return
535        end        end
   
   
   
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.22