/[MITgcm]/MITgcm/pkg/ctrl/ctrl_getobcsw.F
ViewVC logotype

Diff of /MITgcm/pkg/ctrl/ctrl_getobcsw.F

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

revision 1.9 by mlosch, Wed Jan 19 08:42:06 2011 UTC revision 1.15 by jmc, Fri Aug 10 19:38:57 2012 UTC
# Line 1  Line 1 
1  C $Header$  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CTRL_CPPOPTIONS.h"  #include "CTRL_OPTIONS.h"
5  #ifdef ALLOW_OBCS  #ifdef ALLOW_OBCS
6  # include "OBCS_OPTIONS.h"  # include "OBCS_OPTIONS.h"
7  #endif  #endif
8    
   
9        subroutine ctrl_getobcsw(        subroutine ctrl_getobcsw(
10       I                             mytime,       I                             mytime,
11       I                             myiter,       I                             myiter,
# Line 29  c     ================================== Line 28  c     ==================================
28    
29        implicit none        implicit none
30    
 #ifdef ALLOW_OBCSW_CONTROL  
   
31  c     == global variables ==  c     == global variables ==
32    #ifdef ALLOW_OBCSW_CONTROL
33  #include "EEPARAMS.h"  #include "EEPARAMS.h"
34  #include "SIZE.h"  #include "SIZE.h"
35  #include "PARAMS.h"  #include "PARAMS.h"
36  #include "GRID.h"  #include "GRID.h"
37  #include "OBCS.h"  c#include "OBCS_PARAMS.h"
38    #include "OBCS_GRID.h"
39    #include "OBCS_FIELDS.h"
40    #include "CTRL_SIZE.h"
41  #include "ctrl.h"  #include "ctrl.h"
42  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
43  #include "optim.h"  #include "optim.h"
44    #endif /* ALLOW_OBCSW_CONTROL */
45    
46  c     == routine arguments ==  c     == routine arguments ==
   
47        _RL     mytime        _RL     mytime
48        integer myiter        integer myiter
49        integer mythid        integer mythid
50    
51    #ifdef ALLOW_OBCSW_CONTROL
52  c     == local variables ==  c     == local variables ==
53    
54        integer bi,bj        integer bi,bj
# Line 69  c     == local variables == Line 69  c     == local variables ==
69        integer obcswcount1        integer obcswcount1
70    
71  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
72          _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy)
73    
74        logical doglobalread        logical doglobalread
75        logical ladinit        logical ladinit
76    
77        character*(80) fnameobcsw        character*(80) fnameobcsw
78    
79  cgg(  Variables for splitting barotropic/baroclinic vels.  #ifdef ALLOW_OBCS_CONTROL_MODES
80        _RL ubaro        integer nk,nz
81        _RL utop        _RL     tmpz (nr,nsx,nsy)
82  cgg)        _RL     stmp
83    #endif
84    
85  c     == external functions ==  c     == external functions ==
86    
# Line 98  c     == end of interface == Line 100  c     == end of interface ==
100        imax = snx+olx        imax = snx+olx
101        ip1  = 1        ip1  = 1
102    
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       utop = 0.d0  
 cgg)  
   
103  c--   Now, read the control vector.  c--   Now, read the control vector.
104        doglobalread = .false.        doglobalread = .false.
105        ladinit      = .false.        ladinit      = .false.
106    
107        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
108          ilobcsw=ilnblnk( xx_obcsw_file )         ilobcsw=ilnblnk( xx_obcsw_file )
109          write(fnameobcsw(1:80),'(2a,i10.10)')         write(fnameobcsw(1:80),'(2a,i10.10)')
110       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle       &      xx_obcsw_file(1:ilobcsw), '.', optimcycle
111        endif        endif
112    
113  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
# Line 121  c--   Get the counters, flags, and the i Line 118  c--   Get the counters, flags, and the i
118       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
119    
120        do iobcs = 1,nobcs        do iobcs = 1,nobcs
121          if ( obcswfirst ) then         if ( obcswfirst ) then
122            call active_read_yz( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
123       &                         (obcswcount0-1)*nobcs+iobcs,       &                       (obcswcount0-1)*nobcs+iobcs,
124       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
125       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
126    
127  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
128             do bi = itlo,ithi
129            if ( optimcycle .gt. 0) then  #ifdef ALLOW_OBCS_CONTROL_MODES
130              if (iobcs .eq. 3) then            if (iobcs .gt. 2) then
131  cgg         Special attention is needed for the normal velocity.             do j = jmin,jmax
132  cgg         For the north, this is the v velocity, iobcs = 4.              i = OB_Iw(j,bi,bj)
133  cgg         This is done on a columnwise basis here.  cih    Determine number of open vertical layers.
134                do bj = jtlo,jthi              nz = 0
135                  do bi = itlo, ithi              do k = 1,Nr
136                    do j = jmin,jmax                if (iobcs .eq. 3) then
137                      i = OB_Iw(J,bi,bj)                  nz = nz + maskW(i+ip1,j,k,bi,bj)
138                  else
139  cgg         The barotropic velocity is stored in the level 1.                  nz = nz + maskS(i,j,k,bi,bj)
140                      ubaro = tmpfldyz(j,1,bi,bj)                endif
141                      tmpfldyz(j,1,bi,bj) = 0.d0              end do
142                      utop = 0.d0  cih    Compute absolute velocities from the barotropic-baroclinic modes.
143                do k = 1,Nr
144                      do k = 1,Nr               if (k.le.nz) then
145  cgg    If cells are not full, this should be modified with hFac.                stmp = 0.
146  cgg                do nk = 1,nz
147  cgg    The xx field (tmpfldxz) does not contain the velocity at the                 stmp = stmp +
148  cgg    surface level. This velocity is not independent; it must       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
149  cgg    exactly balance the volume flux, since we are dealing with                end do
150  cgg    the baroclinic velocity structure..                 tmpz(k,bi,bj) = stmp
151                        utop = tmpfldyz(j,k,bi,bj)*               else
152       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop                tmpz(k,bi,bj) = 0.
153  cgg    Add the barotropic velocity component.               end if
                       if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                   tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
             if (iobcs .eq. 4) then  
 cgg         Special attention is needed for the normal velocity.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
 cgg         This is done on a columnwise basis here.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(J,bi,bj)  
   
 cgg         The barotropic velocity is stored in the level 1.  
                     ubaro = tmpfldyz(j,1,bi,bj)  
                     tmpfldyz(j,1,bi,bj) = 0.d0  
                     utop = 0.d0  
   
                     do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg  
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    surface level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                       utop = tmpfldyz(j,k,bi,bj)*  
      &                maskS(i,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                   tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
   
 #endif /* ALLOW_CTRL_OBCS_BALANCE */  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
154              enddo              enddo
155                do k = 1,Nr
156                  if (iobcs .eq. 3) then
157                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
158         &            *recip_hFacW(i+ip1,j,k,bi,bj)
159                  else
160                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
161         &            *recip_hFacS(i,j,k,bi,bj)
162                  endif
163                end do
164               enddo
165              endif
166    #endif
167              do k = 1,nr
168               do j = jmin,jmax
169                xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
170    cgg   &                                        *   maskyz (j,k,bi,bj)
171               enddo
172            enddo            enddo
173          endif           enddo
174            enddo
175           endif
176    
177          if ( (obcswfirst) .or. (obcswchanged)) then         if ( (obcswfirst) .or. (obcswchanged)) then
178    
179            do bj = jtlo,jthi          do bj = jtlo,jthi
180             do bi = itlo,ithi           do bi = itlo,ithi
181              do k = 1,nr            do k = 1,nr
182               do j = jmin,jmax             do j = jmin,jmax
183                xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
184                tmpfldyz (j,k,bi,bj)       = 0. _d 0              tmpfldyz (j,k,bi,bj)       = 0. _d 0
              enddo  
             enddo  
185             enddo             enddo
186            enddo            enddo
187             enddo
188            enddo
189    
190            call active_read_yz( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
191       &                         (obcswcount1-1)*nobcs+iobcs,       &                       (obcswcount1-1)*nobcs+iobcs,
192       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
193       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
194    
195  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
196             do bi = itlo,ithi
197            if ( optimcycle .gt. 0) then  #ifdef ALLOW_OBCS_CONTROL_MODES
198              if (iobcs .eq. 3) then            if (iobcs .gt. 2) then
199  cgg         Special attention is needed for the normal velocity.             do j = jmin,jmax
200  cgg         For the north, this is the v velocity, iobcs = 4.              i = OB_Iw(j,bi,bj)
201  cgg         This is done on a columnwise basis here.  cih    Determine number of open vertical layers.
202                do bj = jtlo,jthi              nz = 0
203                  do bi = itlo, ithi              do k = 1,Nr
204                    do j = jmin,jmax                if (iobcs .eq. 3) then
205                      i = OB_Iw(J,bi,bj)                  nz = nz + maskW(i+ip1,j,k,bi,bj)
206                  else
207  cgg         The barotropic velocity is stored in the level 1.                  nz = nz + maskS(i,j,k,bi,bj)
208                      ubaro = tmpfldyz(j,1,bi,bj)                endif
209                      tmpfldyz(j,1,bi,bj) = 0.d0              end do
210                      utop = 0.d0  cih    Compute absolute velocities from the barotropic-baroclinic modes.
211                do k = 1,Nr
212                      do k = 1,Nr               if (k.le.nz) then
213  cgg    If cells are not full, this should be modified with hFac.                stmp = 0.
214  cgg                do nk = 1,nz
215  cgg    The xx field (tmpfldxz) does not contain the velocity at the                 stmp = stmp +
216  cgg    surface level. This velocity is not independent; it must       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
217  cgg    exactly balance the volume flux, since we are dealing with                end do
218  cgg    the baroclinic velocity structure..                 tmpz(k,bi,bj) = stmp
219                        utop = tmpfldyz(j,k,bi,bj)*               else
220       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop                 tmpz(k,bi,bj) = 0.
221  cgg    Add the barotropic velocity component.               end if
                       if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
             if (iobcs .eq. 4) then  
 cgg         Special attention is needed for the normal velocity.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
 cgg         This is done on a columnwise basis here.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(J,bi,bj)  
   
 cgg         The barotropic velocity is stored in the level 1.  
                     ubaro = tmpfldyz(j,1,bi,bj)  
                     tmpfldyz(j,1,bi,bj) = 0.d0  
                     utop = 0.d0  
   
                     do k = 1,Nr  
 cgg    If cells are not full, this should be modified with hFac.  
 cgg  
 cgg    The xx field (tmpfldxz) does not contain the velocity at the  
 cgg    surface level. This velocity is not independent; it must  
 cgg    exactly balance the volume flux, since we are dealing with  
 cgg    the baroclinic velocity structure..  
                       utop = tmpfldyz(j,k,bi,bj)*  
      &                maskS(i,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       if (maskS(i,j,k,bi,bj) .ne. 0.) then  
                         tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro  
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delR(1)  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
   
 #endif /* ALLOW_CTRL_OBCS_BALANCE */  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
222              enddo              enddo
223                do k = 1,Nr
224                  if (iobcs .eq. 3) then
225                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
226         &            *recip_hFacW(i+ip1,j,k,bi,bj)
227                  else
228                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
229         &            *recip_hFacS(i,j,k,bi,bj)
230                  endif
231                end do
232               enddo
233              endif
234    #endif
235              do k = 1,nr
236               do j = jmin,jmax
237                xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
238    cgg   &                                        *   maskyz (j,k,bi,bj)
239               enddo
240            enddo            enddo
241          endif           enddo
242            enddo
243           endif
244    
245  c--     Add control to model variable.  c--   Add control to model variable.
246          do bj = jtlo, jthi         do bj = jtlo, jthi
247             do bi = itlo, ithi          do bi = itlo, ithi
248  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--   Calculate mask for tracer cells (0 => land, 1 => water).
249                do k = 1,nr           do k = 1,nr
250                   do j = 1,sny            do j = 1,sny
251                      i = OB_Iw(j,bi,bj)             i = OB_Iw(j,bi,bj)
252                      if (iobcs .EQ. 1) then             if (iobcs .EQ. 1) then
253                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
254       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
255       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
256                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
257       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
258                      else if (iobcs .EQ. 2) then             else if (iobcs .EQ. 2) then
259                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
260       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
261       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
262                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
263       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
264                      else if (iobcs .EQ. 3) then             else if (iobcs .EQ. 3) then
265                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
266       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
267       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
268                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
269       &                      *maskW(i+ip1,j,k,bi,bj)       &           *maskW(i+ip1,j,k,bi,bj)
270                      else if (iobcs .EQ. 4) then             else if (iobcs .EQ. 4) then
271                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
272       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
273       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
274                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
275       &                      *maskS(i,j,k,bi,bj)       &           *maskS(i,j,k,bi,bj)
276                      endif             endif
277                   enddo            enddo
278                enddo           enddo
            enddo  
279          enddo          enddo
280           enddo
281    
282  C--   End over iobcs loop  C--   End over iobcs loop
283        enddo        enddo
284    
 #else /* ALLOW_OBCSW_CONTROL undefined */  
   
 c     == routine arguments ==  
   
       _RL     mytime  
       integer myiter  
       integer mythid  
   
 c--   CPP flag ALLOW_OBCSW_CONTROL undefined.  
   
285  #endif /* ALLOW_OBCSW_CONTROL */  #endif /* ALLOW_OBCSW_CONTROL */
286    
287          return
288        end        end
   
Legend:
Removed from v.1.9  
changed lines
  Added in v.1.15
  ViewVC Help
Powered by ViewVC 1.1.22