/[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.2 by heimbach, Tue Jun 24 16:07:06 2003 UTC revision 1.13 by jmc, Tue May 24 20:48:28 2011 UTC
# Line 1  Line 1 
1    C $Header$
2    C $Name$
3    
4  #include "CTRL_CPPOPTIONS.h"  #include "CTRL_CPPOPTIONS.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    
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 67  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 96  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 119  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            if ( optimcycle .gt. 0) then                    do bj = jtlo,jthi
128              if (iobcs .eq. 3) then           do bi = itlo,ithi
129  cgg         Special attention is needed for the normal velocity.  #ifdef ALLOW_OBCS_CONTROL_MODES
130  cgg         For the north, this is the v velocity, iobcs = 4.            if (iobcs .gt. 2) then
131  cgg         This is done on a columnwise basis here.             do j = jmin,jmax
132                do bj = jtlo,jthi              i = OB_Iw(j,bi,bj)
133                  do bi = itlo, ithi  cih    Determine number of open vertical layers.
134                    do j = jmin,jmax              nz = 0
135                      i = OB_Iw(J,bi,bj)              do k = 1,Nr
136                  if (iobcs .eq. 3) then
137  cgg         The barotropic velocity is stored in the level 1.                  nz = nz + maskW(i+ip1,j,k,bi,bj)
138                      ubaro = tmpfldyz(j,1,bi,bj)                else
139                      tmpfldyz(j,1,bi,bj) = 0.d0                  nz = nz + maskS(i,j,k,bi,bj)
140                      utop = 0.d0                endif
141                end do
142                      do k = 1,Nr  cih    Compute absolute velocities from the barotropic-baroclinic modes.
143  cgg    If cells are not full, this should be modified with hFac.              do k = 1,Nr
144  cgg                   if (k.le.nz) then
145  cgg    The xx field (tmpfldxz) does not contain the velocity at the                stmp = 0.
146  cgg    surface level. This velocity is not independent; it must                do nk = 1,nz
147  cgg    exactly balance the volume flux, since we are dealing with                 stmp = stmp +
148  cgg    the baroclinic velocity structure..       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
149                        utop = tmpfldyz(j,k,bi,bj)*                end do
150       &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop                 tmpz(k,bi,bj) = stmp
151  cgg    Add the barotropic velocity component.               else
152                        if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then                tmpz(k,bi,bj) = 0.
153                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro               end if
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                   tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delZ(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) * delZ(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 / delZ(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 don't 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  
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             enddo
174            enddo
175           endif
176    
177            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)         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) = tmpfldyz2(j,k,bi,bj)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
184                  enddo              tmpfldyz (j,k,bi,bj)       = 0. _d 0
185                enddo             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            if ( optimcycle .gt. 0) then                    do bj = jtlo,jthi
196              if (iobcs .eq. 3) then           do bi = itlo,ithi
197  cgg         Special attention is needed for the normal velocity.  #ifdef ALLOW_OBCS_CONTROL_MODES
198  cgg         For the north, this is the v velocity, iobcs = 4.            if (iobcs .gt. 2) then
199  cgg         This is done on a columnwise basis here.             do j = jmin,jmax
200                do bj = jtlo,jthi              i = OB_Iw(j,bi,bj)
201                  do bi = itlo, ithi  cih    Determine number of open vertical layers.
202                    do j = jmin,jmax              nz = 0
203                      i = OB_Iw(J,bi,bj)              do k = 1,Nr
204                  if (iobcs .eq. 3) then
205  cgg         The barotropic velocity is stored in the level 1.                  nz = nz + maskW(i+ip1,j,k,bi,bj)
206                      ubaro = tmpfldyz(j,1,bi,bj)                else
207                      tmpfldyz(j,1,bi,bj) = 0.d0                  nz = nz + maskS(i,j,k,bi,bj)
208                      utop = 0.d0                endif
209                end do
210                      do k = 1,Nr  cih    Compute absolute velocities from the barotropic-baroclinic modes.
211  cgg    If cells are not full, this should be modified with hFac.              do k = 1,Nr
212  cgg                   if (k.le.nz) then
213  cgg    The xx field (tmpfldxz) does not contain the velocity at the                stmp = 0.
214  cgg    surface level. This velocity is not independent; it must                do nk = 1,nz
215  cgg    exactly balance the volume flux, since we are dealing with                 stmp = stmp +
216  cgg    the baroclinic velocity structure..       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
217                        utop = tmpfldyz(j,k,bi,bj)*                end do
218       &                maskW(i+ip1,j,k,bi,bj) * delZ(k) + utop                 tmpz(k,bi,bj) = stmp
219  cgg    Add the barotropic velocity component.               else
220                        if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then                 tmpz(k,bi,bj) = 0.
221                          tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj)+ ubaro               end if
                       endif  
                     enddo  
 cgg    Compute the baroclinic velocity at level 1. Should balance flux.  
                     tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj)  
      &                                      - utop / delZ(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) * delZ(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 / delZ(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  
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.2  
changed lines
  Added in v.1.13

  ViewVC Help
Powered by ViewVC 1.1.22