/[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.8 by jmc, Tue Oct 9 00:00:00 2007 UTC revision 1.16 by jmc, Tue Sep 18 20:21:23 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.              IF ( i.EQ.OB_indexNone ) i = 1
134                do bj = jtlo,jthi  cih    Determine number of open vertical layers.
135                  do bi = itlo, ithi              nz = 0
136                    do j = jmin,jmax              do k = 1,Nr
137                      i = OB_Iw(J,bi,bj)                if (iobcs .eq. 3) then
138                    nz = nz + maskW(i+ip1,j,k,bi,bj)
139  cgg         The barotropic velocity is stored in the level 1.                else
140                      ubaro = tmpfldyz(j,1,bi,bj)                  nz = nz + maskS(i,j,k,bi,bj)
141                      tmpfldyz(j,1,bi,bj) = 0.d0                endif
142                      utop = 0.d0              end do
143    cih    Compute absolute velocities from the barotropic-baroclinic modes.
144                      do k = 1,Nr              do k = 1,Nr
145  cgg    If cells are not full, this should be modified with hFac.               if (k.le.nz) then
146  cgg                stmp = 0.
147  cgg    The xx field (tmpfldxz) does not contain the velocity at the                do nk = 1,nz
148  cgg    surface level. This velocity is not independent; it must                 stmp = stmp +
149  cgg    exactly balance the volume flux, since we are dealing with       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
150  cgg    the baroclinic velocity structure..                end do
151                        utop = tmpfldyz(j,k,bi,bj)*                 tmpz(k,bi,bj) = stmp
152       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop               else
153  cgg    Add the barotropic velocity component.                tmpz(k,bi,bj) = 0.
154                        if (maskW(i+ip1,j,k,bi,bj) .ne. 0.) then               end if
                         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  
155              enddo              enddo
156                do k = 1,Nr
157                  if (iobcs .eq. 3) then
158                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
159         &            *recip_hFacW(i+ip1,j,k,bi,bj)
160                  else
161                    tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
162         &            *recip_hFacS(i,j,k,bi,bj)
163                  endif
164                end do
165               enddo
166              endif
167    #endif
168              do k = 1,nr
169               do j = jmin,jmax
170                xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
171    cgg   &                                        *   maskyz (j,k,bi,bj)
172               enddo
173            enddo            enddo
174          endif           enddo
175            enddo
176          if ( (obcswfirst) .or. (obcswchanged)) then         endif
177    
178  cgg(    This is a terribly long way to do it. However, the dimensions do not exactly         if ( (obcswfirst) .or. (obcswchanged)) then
 cgg     match up. I will blame Fortran for the ugliness.  
179    
180            do bj = jtlo,jthi          do bj = jtlo,jthi
181              do bi = itlo,ithi           do bi = itlo,ithi
182                do k = 1,nr            do k = 1,nr
183                  do j = jmin,jmax             do j = jmin,jmax
184                    tmpfldyz(j,k,bi,bj) = xx_obcsw1(j,k,bi,bj,iobcs)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
185                  enddo              tmpfldyz (j,k,bi,bj)       = 0. _d 0
186                enddo             enddo
             enddo  
187            enddo            enddo
188             enddo
189            enddo
190    
191            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)          call active_read_yz( fnameobcsw, tmpfldyz,
192         &                       (obcswcount1-1)*nobcs+iobcs,
193            do bj = jtlo,jthi       &                       doglobalread, ladinit, optimcycle,
194              do bi = itlo,ithi       &                       mythid, xx_obcsw_dummy )
195                do k = 1,nr  
196                  do j = jmin,jmax          do bj = jtlo,jthi
197                    xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)           do bi = itlo,ithi
198                  enddo  #ifdef ALLOW_OBCS_CONTROL_MODES
199                enddo            if (iobcs .gt. 2) then
200               do j = jmin,jmax
201                i = OB_Iw(j,bi,bj)
202                IF ( i.EQ.OB_indexNone ) i = 1
203    cih    Determine number of open vertical layers.
204                nz = 0
205                do k = 1,Nr
206                  if (iobcs .eq. 3) then
207                    nz = nz + maskW(i+ip1,j,k,bi,bj)
208                  else
209                    nz = nz + maskS(i,j,k,bi,bj)
210                  endif
211                end do
212    cih    Compute absolute velocities from the barotropic-baroclinic modes.
213                do k = 1,Nr
214                 if (k.le.nz) then
215                  stmp = 0.
216                  do nk = 1,nz
217                   stmp = stmp +
218         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
219                  end do
220                   tmpz(k,bi,bj) = stmp
221                 else
222                   tmpz(k,bi,bj) = 0.
223                 end if
224              enddo              enddo
225            enddo              do k = 1,Nr
226                  if (iobcs .eq. 3) then
227            call active_read_yz( fnameobcsw, tmpfldyz,                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
228       &                         (obcswcount1-1)*nobcs+iobcs,       &            *recip_hFacW(i+ip1,j,k,bi,bj)
229       &                         doglobalread, ladinit, optimcycle,                else
230       &                         mythid, xx_obcsw_dummy )                  tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj)
231         &            *recip_hFacS(i,j,k,bi,bj)
232  #ifdef ALLOW_CTRL_OBCS_BALANCE                endif
233                end do
234            if ( optimcycle .gt. 0) then             enddo
             if (iobcs .eq. 3) 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)*  
      &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop  
 cgg    Add the barotropic velocity component.  
                       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  
235            endif            endif
236    #endif
237  #endif /* ALLOW_CTRL_OBCS_BALANCE */            do k = 1,nr
238               do j = jmin,jmax
239            do bj = jtlo,jthi              xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
240              do bi = itlo,ithi  cgg   &                                        *   maskyz (j,k,bi,bj)
241                do k = 1,nr             enddo
                 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  
242            enddo            enddo
243          endif           enddo
244            enddo
245           endif
246    
247  c--     Add control to model variable.  c--   Add control to model variable.
248          do bj = jtlo, jthi         do bj = jtlo, jthi
249             do bi = itlo, ithi          do bi = itlo, ithi
250  c--        Calculate mask for tracer cells (0 => land, 1 => water).  c--   Calculate mask for tracer cells (0 => land, 1 => water).
251                do k = 1,nr           do k = 1,nr
252                   do j = 1,sny            do j = 1,sny
253                      i = OB_Iw(j,bi,bj)             i = OB_Iw(j,bi,bj)
254                      if (iobcs .EQ. 1) then             IF ( i.EQ.OB_indexNone ) i = 1
255                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)             if (iobcs .EQ. 1) then
256       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
257       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
258                         OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
259       &                      *maskW(i+ip1,j,k,bi,bj)              OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
260                      else if (iobcs .EQ. 2) then       &           *maskW(i+ip1,j,k,bi,bj)
261                         OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)             else if (iobcs .EQ. 2) then
262       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
263       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
264                         OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
265       &                      *maskW(i+ip1,j,k,bi,bj)              OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
266                      else if (iobcs .EQ. 3) then       &           *maskW(i+ip1,j,k,bi,bj)
267                         OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)             else if (iobcs .EQ. 3) then
268       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
269       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
270                         OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
271       &                      *maskW(i+ip1,j,k,bi,bj)              OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
272                      else if (iobcs .EQ. 4) then       &           *maskW(i+ip1,j,k,bi,bj)
273                         OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)             else if (iobcs .EQ. 4) then
274       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
275       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)       &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
276                         OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)       &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
277       &                      *maskS(i,j,k,bi,bj)              OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
278                      endif       &           *maskS(i,j,k,bi,bj)
279                   enddo             endif
280                enddo            enddo
281             enddo           enddo
282          enddo          enddo
283           enddo
284    
285  C--   End over iobcs loop  C--   End over iobcs loop
286        enddo        enddo
287    
 #else /* ALLOW_OBCSW_CONTROL undefined */  
   
 c     == routine arguments ==  
   
       _RL     mytime  
       integer myiter  
       integer mythid  
   
 c--   CPP flag ALLOW_OBCSW_CONTROL undefined.  
   
288  #endif /* ALLOW_OBCSW_CONTROL */  #endif /* ALLOW_OBCSW_CONTROL */
289    
290          return
291        end        end
   
Legend:
Removed from v.1.8  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22