/[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.1.2.2 by heimbach, Thu May 30 19:56:44 2002 UTC revision 1.17 by gforget, Thu Oct 9 00:49:26 2014 UTC
# Line 1  Line 1 
1    C $Header$
2    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 20  c       to dyn. fields Line 21  c       to dyn. fields
21  c  c
22  c     started: heimbach@mit.edu, 29-Aug-2001  c     started: heimbach@mit.edu, 29-Aug-2001
23  c  c
24    c     modified: gebbie@mit.edu, 18-Mar-2003
25  c     ==================================================================  c     ==================================================================
26  c     SUBROUTINE ctrl_getobcsw  c     SUBROUTINE ctrl_getobcsw
27  c     ==================================================================  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 "CTRL_OBCS.h"
44  #include "optim.h"  #include "optim.h"
45    #endif /* ALLOW_OBCSW_CONTROL */
46    
47  c     == routine arguments ==  c     == routine arguments ==
   
48        _RL     mytime        _RL     mytime
49        integer myiter        integer myiter
50        integer mythid        integer mythid
51    
52    #ifdef ALLOW_OBCSW_CONTROL
53  c     == local variables ==  c     == local variables ==
54    
55        integer bi,bj        integer bi,bj
# Line 66  c     == local variables == Line 70  c     == local variables ==
70        integer obcswcount1        integer obcswcount1
71    
72  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
73          _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy)
74    
75        logical doglobalread        logical doglobalread
76        logical ladinit        logical ladinit
77    
78        character*(80) fnameobcsw        character*(80) fnameobcsw
79    
80  #ifdef BALANCE_CONTROL_VOLFLUX  #ifdef ALLOW_OBCS_CONTROL_MODES
81  cgg(  Variables for balancing volume flux.        integer nk,nz
82        _RL ubaro        _RL     tmpz (nr,nsx,nsy)
83        _RL ubarocount        _RL     stmp
 cgg)  
84  #endif  #endif
85    
86  c     == external functions ==  c     == external functions ==
# Line 97  c     == end of interface == Line 101  c     == end of interface ==
101        imax = snx+olx        imax = snx+olx
102        ip1  = 1        ip1  = 1
103    
 #ifdef BALANCE_CONTROL_VOLFLUX  
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       ubarocount = 0.d0  
 cgg)  
 #endif  
   
104  c--   Now, read the control vector.  c--   Now, read the control vector.
105        doglobalread = .false.        doglobalread = .false.
106        ladinit      = .false.        ladinit      = .false.
107    
108        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
109          ilobcsw=ilnblnk( xx_obcsw_file )         ilobcsw=ilnblnk( xx_obcsw_file )
110          write(fnameobcsw(1:80),'(2a,i10.10)')         write(fnameobcsw(1:80),'(2a,i10.10)')
111       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle       &      xx_obcsw_file(1:ilobcsw), '.', optimcycle
       else  
         print*  
         print*,' ctrl_getobcsw: optimcycle not set correctly.'  
         print*,' ... stopped in ctrl_getobcsw.'  
112        endif        endif
113    
114  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
115        call ctrl_GetRec( 'xx_obcsw',        call ctrl_get_gen_rec(
116         I                   xx_obcswstartdate, xx_obcswperiod,
117       O                   obcswfac, obcswfirst, obcswchanged,       O                   obcswfac, obcswfirst, obcswchanged,
118       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
119       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
120    
121        do iobcs = 1,nobcs        do iobcs = 1,nobcs
122                 if ( obcswfirst ) then
123  cgg        if ( (obcswfirst) .or. (obcswchanged) ) then          call active_read_yz( fnameobcsw, tmpfldyz,
124  cgg          call active_read_yz( 'maskobcsw', maskyz,       &                       (obcswcount0-1)*nobcs+iobcs,
125  cgg     &                         iobcs,       &                       doglobalread, ladinit, optimcycle,
126  cgg     &                         doglobalread, ladinit, 0,       &                       mythid, xx_obcsw_dummy )
 cgg     &                         mythid, dummy )  
 cgg        endif  
   
         if ( obcswfirst ) then  
           call active_read_yz( fnameobcsw, tmpfldyz,  
      &                         (obcswcount0-1)*nobcs+iobcs,  
      &                         doglobalread, ladinit, optimcycle,  
      &                         mythid, xx_obcsw_dummy )  
   
 #ifdef BALANCE_CONTROL_VOLFLUX  
           if ( optimcycle .gt. 0) then            
             if (iobcs .eq. 3) then  
 cgg(        Make sure that the xx velocity field has a balanced net volume flux.  
 cgg         Find the barotropic flow normal to the boundary.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(j,bi,bj)  
                     ubaro = 0.d0  
                     ubarocount = 0.d0  
                     do k = 1,nr  
 cgg(   If cells are not full, this should be modified with hFac.  
                       ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj)  
      &                      * delZ(k) + ubaro  
                       ubarocount = maskW(i+ip1,j,k,bi,bj)  
      &                           * delZ(k) +ubarocount  
                     enddo  
                     if (ubarocount .eq. 0.) then  
                       ubaro = 0.  
                     else  
                       ubaro = ubaro / ubarocount  
                     endif  
                     do k = 1,nr  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 cgg)  
 #endif  
 #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do k = 1,Nr  
                     do j = jmin,jmax  
                       i = OB_Iw(j,bi,bj)  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  
      &                   - shiftvel(1) *maskW(i+ip1,j,k,bi,bj)  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 #endif  
127    
128            do bj = jtlo,jthi          do bj = jtlo,jthi
129              do bi = itlo,ithi           do bi = itlo,ithi
130                do k = 1,nr  #ifdef ALLOW_OBCS_CONTROL_MODES
131                  do j = jmin,jmax            if (iobcs .gt. 2) then
132                    xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)             do j = jmin,jmax
133                   enddo              i = OB_Iw(j,bi,bj)
134                enddo              IF ( i.EQ.OB_indexNone ) i = 1
135    cih    Determine number of open vertical layers.
136                nz = 0
137                do k = 1,Nr
138                  if (iobcs .eq. 3) then
139                    nz = nz + maskW(i+ip1,j,k,bi,bj)
140                  else
141                    nz = nz + maskS(i,j,k,bi,bj)
142                  endif
143                end do
144    cih    Compute absolute velocities from the barotropic-baroclinic modes.
145                do k = 1,Nr
146                 if (k.le.nz) then
147                  stmp = 0.
148                  do nk = 1,nz
149                   stmp = stmp +
150         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
151                  end do
152                   tmpz(k,bi,bj) = stmp
153                 else
154                  tmpz(k,bi,bj) = 0.
155                 end if
156              enddo              enddo
157                do k = 1,Nr
158                  if (iobcs .eq. 3) then
159                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
160         &            *recip_hFacW(i+ip1,j,k,bi,bj)
161                  else
162                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
163         &            *recip_hFacS(i,j,k,bi,bj)
164                  endif
165                end do
166               enddo
167              endif
168    #endif
169              do k = 1,nr
170               do j = jmin,jmax
171                xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
172    cgg   &                                        *   maskyz (j,k,bi,bj)
173               enddo
174            enddo            enddo
175          endif           enddo
176            enddo
177           endif
178    
179          if ( (obcswfirst) .or. (obcswchanged)) then         if ( (obcswfirst) .or. (obcswchanged)) then
180              
181  cgg(    This is a terribly long way to do it. However, the dimensions don't exactly          do bj = jtlo,jthi
182  cgg     match up. I will blame Fortran for the ugliness.           do bi = itlo,ithi
183              do k = 1,nr
184            do bj = jtlo,jthi             do j = jmin,jmax
185              do bi = itlo,ithi              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
186                do k = 1,nr              tmpfldyz (j,k,bi,bj)       = 0. _d 0
187                  do j = jmin,jmax             enddo
   
                   tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)  
                 enddo  
               enddo  
             enddo  
188            enddo            enddo
189             enddo
190            enddo
191    
192            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)          call active_read_yz( fnameobcsw, tmpfldyz,
193         &                       (obcswcount1-1)*nobcs+iobcs,
194         &                       doglobalread, ladinit, optimcycle,
195         &                       mythid, xx_obcsw_dummy )
196    
197            do bj = jtlo,jthi          do bj = jtlo,jthi
198              do bi = itlo,ithi           do bi = itlo,ithi
199                do k = 1,nr  #ifdef ALLOW_OBCS_CONTROL_MODES
200                  do j = jmin,jmax            if (iobcs .gt. 2) then
201                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)                     do j = jmin,jmax
202                  enddo              i = OB_Iw(j,bi,bj)
203                enddo              IF ( i.EQ.OB_indexNone ) i = 1
204    cih    Determine number of open vertical layers.
205                nz = 0
206                do k = 1,Nr
207                  if (iobcs .eq. 3) then
208                    nz = nz + maskW(i+ip1,j,k,bi,bj)
209                  else
210                    nz = nz + maskS(i,j,k,bi,bj)
211                  endif
212                end do
213    cih    Compute absolute velocities from the barotropic-baroclinic modes.
214                do k = 1,Nr
215                 if (k.le.nz) then
216                  stmp = 0.
217                  do nk = 1,nz
218                   stmp = stmp +
219         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
220                  end do
221                   tmpz(k,bi,bj) = stmp
222                 else
223                   tmpz(k,bi,bj) = 0.
224                 end if
225              enddo              enddo
226            enddo              do k = 1,Nr
227                  if (iobcs .eq. 3) then
228            call active_read_yz( fnameobcsw, tmpfldyz,                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
229       &                         (obcswcount1-1)*nobcs+iobcs,       &            *recip_hFacW(i+ip1,j,k,bi,bj)
230       &                         doglobalread, ladinit, optimcycle,                else
231       &                         mythid, xx_obcsw_dummy )                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
232         &            *recip_hFacS(i,j,k,bi,bj)
233  #ifdef BALANCE_CONTROL_VOLFLUX                endif
234            if (optimcycle .gt. 0) then              end do
235              if (iobcs .eq. 3) then             enddo
 cgg(        Make sure that the xx velocity field has a balanced net volume flux.  
 cgg         Find the barotropic flow normal to the boundary.  
 cgg         For the north, this is the v velocity, iobcs = 4.  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do j = jmin,jmax  
                     i = OB_Iw(j,bi,bj)  
                     ubaro = 0.d0  
                     ubarocount = 0.d0  
                     do k = 1,nr  
 cgg(   If cells are not full, this should be modified with hFac.  
                       ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj)  
      &                      * delZ(k) + ubaro  
                       ubarocount = maskW(i+ip1,j,k,bi,bj)  
      &                           * delZ(k) + ubarocount  
                     enddo  
                     if (ubarocount .eq. 0.) then  
                       ubaro = 0.  
                     else  
                       ubaro = ubaro / ubarocount  
                     endif  
                     do k = 1,nr  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
           endif  
 cgg)  
 #endif  
 #ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL  
           if (optimcycle .gt. 0) then  
             if ( iobcs .eq. 3) then  
               do bj = jtlo,jthi  
                 do bi = itlo, ithi  
                   do k = 1,Nr  
                     do j = jmin,jmax  
                       i = OB_Iw(j,bi,bj)  
                       tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)  
      &                   - shiftvel(2) *maskW(i+ip1,j,k,bi,bj)  
                     enddo  
                   enddo  
                 enddo  
               enddo  
             endif  
236            endif            endif
237  #endif  #endif
238              do k = 1,nr
239            do bj = jtlo,jthi             do j = jmin,jmax
240              do bi = itlo,ithi              xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
241                do k = 1,nr  cgg   &                                        *   maskyz (j,k,bi,bj)
242                  do j = jmin,jmax             enddo
                   xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)  
                  enddo  
               enddo  
             enddo  
243            enddo            enddo
244          endif           enddo
245            enddo
246           endif
247    
248  c--     Add control to model variable.  c--   Add control to model variable.
249          do bj = jtlo,jthi         do bj = jtlo, jthi
250             do bi = itlo,ithi          do bi = itlo, ithi
251  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--   Calculate mask for tracer cells (0 => land, 1 => water).
252                do k = 1,nr           do k = 1,nr
253                   do j = 1,sny            do j = 1,sny
254                      i = OB_Iw(j,bi,bj)             i = OB_Iw(j,bi,bj)
255                      if (iobcs .EQ. 1) then             IF ( i.EQ.OB_indexNone ) i = 1
256                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)             if (iobcs .EQ. 1) then
257       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
258       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
259                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
260       &                      *maskW(i+ip1,j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
261                      else if (iobcs .EQ. 2) then       &           *maskW(i+ip1,j,k,bi,bj)
262                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)             else if (iobcs .EQ. 2) then
263       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
264       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
265                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
266       &                      *maskW(i+ip1,j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
267                      else if (iobcs .EQ. 3) then       &           *maskW(i+ip1,j,k,bi,bj)
268                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)             else if (iobcs .EQ. 3) then
269       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
270       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
271                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
272       &                      *maskW(i+ip1,j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
273                      else if (iobcs .EQ. 4) then       &           *maskW(i+ip1,j,k,bi,bj)
274                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)             else if (iobcs .EQ. 4) then
275       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
276       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
277                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
278       &                      *maskS(i,j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
279                      endif       &           *maskS(i,j,k,bi,bj)
280                   enddo             endif
281                enddo            enddo
282             enddo           enddo
283          enddo          enddo
284           enddo
285    
286  C--   End over iobcs loop  C--   End over iobcs loop
287        enddo        enddo
288    
 #else /* ALLOW_OBCSW_CONTROL undefined */  
   
 c     == routine arguments ==  
   
       _RL     mytime  
       integer myiter  
       integer mythid  
   
 c--   CPP flag ALLOW_OBCSW_CONTROL undefined.  
   
289  #endif /* ALLOW_OBCSW_CONTROL */  #endif /* ALLOW_OBCSW_CONTROL */
290    
291          return
292        end        end
   
Legend:
Removed from v.1.1.2.2  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22