/[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.5 by heimbach, Sun Oct 26 00:58:03 2003 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 27  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 "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 67  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  cgg(  Variables for splitting barotropic/baroclinic vels.  #ifdef ALLOW_OBCS_CONTROL_MODES
81        _RL ubaro        integer nk,nz
82        _RL utop        _RL     tmpz (nr,nsx,nsy)
83  cgg)        _RL     stmp
84    #endif
85    
86  c     == external functions ==  c     == external functions ==
87    
# Line 96  c     == end of interface == Line 101  c     == end of interface ==
101        imax = snx+olx        imax = snx+olx
102        ip1  = 1        ip1  = 1
103    
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       utop = 0.d0  
 cgg)  
   
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
112        endif        endif
113    
114  c--   Get the counters, flags, and the interpolation factor.  c--   Get the counters, flags, and the interpolation factor.
# Line 119  c--   Get the counters, flags, and the i Line 119  c--   Get the counters, flags, and the i
119       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
120    
121        do iobcs = 1,nobcs        do iobcs = 1,nobcs
122          if ( obcswfirst ) then         if ( obcswfirst ) then
123            call active_read_yz_loc( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
124       &                         (obcswcount0-1)*nobcs+iobcs,       &                       (obcswcount0-1)*nobcs+iobcs,
125       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
126       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
127    
128            if ( optimcycle .gt. 0) then                    do bj = jtlo,jthi
129              if (iobcs .eq. 3) then           do bi = itlo,ithi
130  cgg         Special attention is needed for the normal velocity.  #ifdef ALLOW_OBCS_CONTROL_MODES
131  cgg         For the north, this is the v velocity, iobcs = 4.            if (iobcs .gt. 2) then
132  cgg         This is done on a columnwise basis here.             do j = jmin,jmax
133                do bj = jtlo,jthi              i = OB_Iw(j,bi,bj)
134                  do bi = itlo, ithi              IF ( i.EQ.OB_indexNone ) i = 1
135                    do j = jmin,jmax  cih    Determine number of open vertical layers.
136                      i = OB_Iw(J,bi,bj)              nz = 0
137                do k = 1,Nr
138  cgg         The barotropic velocity is stored in the level 1.                if (iobcs .eq. 3) then
139                      ubaro = tmpfldyz(j,1,bi,bj)                  nz = nz + maskW(i+ip1,j,k,bi,bj)
140                      tmpfldyz(j,1,bi,bj) = 0.d0                else
141                      utop = 0.d0                  nz = nz + maskS(i,j,k,bi,bj)
142                  endif
143                      do k = 1,Nr              end do
144  cgg    If cells are not full, this should be modified with hFac.  cih    Compute absolute velocities from the barotropic-baroclinic modes.
145  cgg                  do k = 1,Nr
146  cgg    The xx field (tmpfldxz) does not contain the velocity at the               if (k.le.nz) then
147  cgg    surface level. This velocity is not independent; it must                stmp = 0.
148  cgg    exactly balance the volume flux, since we are dealing with                do nk = 1,nz
149  cgg    the baroclinic velocity structure..                 stmp = stmp +
150                        utop = tmpfldyz(j,k,bi,bj)*       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
151       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop                end do
152  cgg    Add the barotropic velocity component.                 tmpz(k,bi,bj) = stmp
153                        if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then               else
154                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro                tmpz(k,bi,bj) = 0.
155                        endif               end if
                     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  
   
           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  
             enddo  
           enddo  
         endif  
   
         if ( (obcswfirst) .or. (obcswchanged)) then  
             
 cgg(    This is a terribly long way to do it. However, the dimensions do not exactly  
 cgg     match up. I will blame Fortran for the ugliness.  
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)  
                 enddo  
               enddo  
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             enddo
176            enddo
177           endif
178    
179            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)         if ( (obcswfirst) .or. (obcswchanged)) then
180    
181            do bj = jtlo,jthi          do bj = jtlo,jthi
182              do bi = itlo,ithi           do bi = itlo,ithi
183                do k = 1,nr            do k = 1,nr
184                  do j = jmin,jmax             do j = jmin,jmax
185                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
186                  enddo              tmpfldyz (j,k,bi,bj)       = 0. _d 0
187                enddo             enddo
             enddo  
188            enddo            enddo
189             enddo
190            enddo
191    
192            call active_read_yz_loc( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
193       &                         (obcswcount1-1)*nobcs+iobcs,       &                       (obcswcount1-1)*nobcs+iobcs,
194       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
195       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
196    
197            if ( optimcycle .gt. 0) then                    do bj = jtlo,jthi
198              if (iobcs .eq. 3) then           do bi = itlo,ithi
199  cgg         Special attention is needed for the normal velocity.  #ifdef ALLOW_OBCS_CONTROL_MODES
200  cgg         For the north, this is the v velocity, iobcs = 4.            if (iobcs .gt. 2) then
201  cgg         This is done on a columnwise basis here.             do j = jmin,jmax
202                do bj = jtlo,jthi              i = OB_Iw(j,bi,bj)
203                  do bi = itlo, ithi              IF ( i.EQ.OB_indexNone ) i = 1
204                    do j = jmin,jmax  cih    Determine number of open vertical layers.
205                      i = OB_Iw(J,bi,bj)              nz = 0
206                do k = 1,Nr
207  cgg         The barotropic velocity is stored in the level 1.                if (iobcs .eq. 3) then
208                      ubaro = tmpfldyz(j,1,bi,bj)                  nz = nz + maskW(i+ip1,j,k,bi,bj)
209                      tmpfldyz(j,1,bi,bj) = 0.d0                else
210                      utop = 0.d0                  nz = nz + maskS(i,j,k,bi,bj)
211                  endif
212                      do k = 1,Nr              end do
213  cgg    If cells are not full, this should be modified with hFac.  cih    Compute absolute velocities from the barotropic-baroclinic modes.
214  cgg                  do k = 1,Nr
215  cgg    The xx field (tmpfldxz) does not contain the velocity at the               if (k.le.nz) then
216  cgg    surface level. This velocity is not independent; it must                stmp = 0.
217  cgg    exactly balance the volume flux, since we are dealing with                do nk = 1,nz
218  cgg    the baroclinic velocity structure..                 stmp = stmp +
219                        utop = tmpfldyz(j,k,bi,bj)*       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
220       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop                end do
221  cgg    Add the barotropic velocity component.                 tmpz(k,bi,bj) = stmp
222                        if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then               else
223                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro                 tmpz(k,bi,bj) = 0.
224                        endif               end if
                     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  
   
           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  
225              enddo              enddo
226                do k = 1,Nr
227                  if (iobcs .eq. 3) then
228                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
229         &            *recip_hFacW(i+ip1,j,k,bi,bj)
230                  else
231                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
232         &            *recip_hFacS(i,j,k,bi,bj)
233                  endif
234                end do
235               enddo
236              endif
237    #endif
238              do k = 1,nr
239               do j = jmin,jmax
240                xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
241    cgg   &                                        *   maskyz (j,k,bi,bj)
242               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.5  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22