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

Diff of /MITgcm/pkg/ctrl/ctrl_getobcse.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.14 by jmc, Fri Aug 10 19:38:57 2012 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_getobcse(        subroutine ctrl_getobcse(
10       I                             mytime,       I                             mytime,
11       I                             myiter,       I                             myiter,
# Line 26  c     ================================== Line 27  c     ==================================
27    
28        implicit none        implicit none
29    
 #ifdef ALLOW_OBCSE_CONTROL  
   
30  c     == global variables ==  c     == global variables ==
31    #ifdef ALLOW_OBCSE_CONTROL
32  #include "EEPARAMS.h"  #include "EEPARAMS.h"
33  #include "SIZE.h"  #include "SIZE.h"
34  #include "PARAMS.h"  #include "PARAMS.h"
35  #include "GRID.h"  #include "GRID.h"
36  #include "OBCS.h"  c#include "OBCS_PARAMS.h"
37    #include "OBCS_GRID.h"
38    #include "OBCS_FIELDS.h"
39    #include "CTRL_SIZE.h"
40  #include "ctrl.h"  #include "ctrl.h"
41  #include "ctrl_dummy.h"  #include "ctrl_dummy.h"
42  #include "optim.h"  #include "optim.h"
43    #endif /* ALLOW_OBCSE_CONTROL */
44    
45  c     == routine arguments ==  c     == routine arguments ==
   
46        _RL     mytime        _RL     mytime
47        integer myiter        integer myiter
48        integer mythid        integer mythid
49    
50    #ifdef ALLOW_OBCSE_CONTROL
51  c     == local variables ==  c     == local variables ==
52    
53        integer bi,bj        integer bi,bj
# Line 66  c     == local variables == Line 68  c     == local variables ==
68        integer ip1        integer ip1
69    
70  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)  cgg      _RL maskyz   (1-oly:sny+oly,nr,nsx,nsy)
71          _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy)
72    
73        logical doglobalread        logical doglobalread
74        logical ladinit        logical ladinit
75    
76        character*(80) fnameobcse        character*(80) fnameobcse
77    
78  cgg(  Variables for splitting barotropic/baroclinic vels.  #ifdef ALLOW_OBCS_CONTROL_MODES
79        _RL ubaro        integer nk,nz
80        _RL utop        _RL     tmpz (nr,nsx,nsy)
81  cgg)        _RL     stmp
82    #endif
83    
84  c     == external functions ==  c     == external functions ==
85    
# Line 95  c     == end of interface == Line 99  c     == end of interface ==
99        imax = snx+olx        imax = snx+olx
100        ip1  = 0        ip1  = 0
101    
 cgg(  Initialize variables for balancing volume flux.  
       ubaro = 0.d0  
       utop = 0.d0  
 cgg)  
   
102  c--   Now, read the control vector.  c--   Now, read the control vector.
103        doglobalread = .false.        doglobalread = .false.
104        ladinit      = .false.        ladinit      = .false.
105    
106        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
107          ilobcse=ilnblnk( xx_obcse_file )         ilobcse=ilnblnk( xx_obcse_file )
108          write(fnameobcse(1:80),'(2a,i10.10)')         write(fnameobcse(1:80),'(2a,i10.10)')
109       &       xx_obcse_file(1:ilobcse), '.', optimcycle       &      xx_obcse_file(1:ilobcse), '.', optimcycle
110        endif        endif
111    
112  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    
119        do iobcs = 1,nobcs        do iobcs = 1,nobcs
120    
121          if ( obcsefirst ) then         if ( obcsefirst ) then
122            call active_read_yz( fnameobcse, tmpfldyz,          call active_read_yz( fnameobcse, tmpfldyz,
123       &                         (obcsecount0-1)*nobcs+iobcs,       &                       (obcsecount0-1)*nobcs+iobcs,
124       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
125       &                         mythid, xx_obcse_dummy )       &                       mythid, xx_obcse_dummy )
   
           if ( optimcycle .gt. 0) then            
             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_Ie(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) * delZ(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 / 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_Ie(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_obcse1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)  
 cgg     &                                        *   maskyz (j,k,bi,bj)  
                  enddo  
               enddo  
             enddo  
           enddo  
         endif  
126    
127          if ( (obcsefirst) .or. (obcsechanged)) then          do bj = jtlo,jthi
128                       do bi = itlo,ithi
129            do bj = jtlo,jthi  #ifdef ALLOW_OBCS_CONTROL_MODES
130              do bi = itlo,ithi            if (iobcs .gt. 2) then
131                do j = jmin,jmax             do j = jmin,jmax
132                  do k = 1,nr              i = OB_Ie(j,bi,bj)
133                    tmpfldyz(j,k,bi,bj) = xx_obcse1(j,k,bi,bj,iobcs)  cih    Determine number of open vertical layers.
134                  enddo              nz = 0
135                enddo              do k = 1,Nr
136                  if (iobcs .eq. 3) then
137                    nz = nz + maskW(i+ip1,j,k,bi,bj)
138                  else
139                    nz = nz + maskS(i,j,k,bi,bj)
140                  endif
141                end do
142    cih    Compute absolute velocities from the barotropic-baroclinic modes.
143                do k = 1,Nr
144                 if (k.le.nz) then
145                  stmp = 0.
146                  do nk = 1,nz
147                   stmp = stmp +
148         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
149                  end do
150                   tmpz(k,bi,bj) = stmp
151                 else
152                  tmpz(k,bi,bj) = 0.
153                 end if
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_obcse1(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 ( (obcsefirst) .or. (obcsechanged)) then
178    
179            do bj = jtlo,jthi          do bj = jtlo,jthi
180              do bi = itlo,ithi           do bi = itlo,ithi
181                do j = jmin,jmax            do j = jmin,jmax
182                  do k = 1,nr             do k = 1,nr
183                    xx_obcse0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)              xx_obcse0(j,k,bi,bj,iobcs) = xx_obcse1(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( fnameobcse, tmpfldyz,          call active_read_yz( fnameobcse, tmpfldyz,
191       &                         (obcsecount1-1)*nobcs+iobcs,       &                       (obcsecount1-1)*nobcs+iobcs,
192       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
193       &                         mythid, xx_obcse_dummy )       &                       mythid, xx_obcse_dummy )
   
           if ( optimcycle .gt. 0) then            
             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_Ie(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) * delZ(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 / 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_Ie(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  
194    
195            do bj = jtlo,jthi          do bj = jtlo,jthi
196              do bi = itlo,ithi           do bi = itlo,ithi
197                do k = 1,nr  #ifdef ALLOW_OBCS_CONTROL_MODES
198                  do j = jmin,jmax            if (iobcs .gt. 2) then
199                    xx_obcse1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)             do j = jmin,jmax
200  cgg     &                                        *   maskyz (j,k,bi,bj)              i = OB_Ie(j,bi,bj)
201                   enddo  cih    Determine number of open vertical layers.
202                enddo              nz = 0
203                do k = 1,Nr
204                  if (iobcs .eq. 3) then
205                    nz = nz + maskW(i+ip1,j,k,bi,bj)
206                  else
207                    nz = nz + maskS(i,j,k,bi,bj)
208                  endif
209                end do
210    cih    Compute absolute velocities from the barotropic-baroclinic modes.
211                do k = 1,Nr
212                 if (k.le.nz) then
213                  stmp = 0.
214                  do nk = 1,nz
215                   stmp = stmp +
216         &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
217                  end do
218                   tmpz(k,bi,bj) = stmp
219                 else
220                  tmpz(k,bi,bj) = 0.
221                 endif
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_obcse1 (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_Ie(j,bi,bj)             i = OB_Ie(j,bi,bj)
252                      if (iobcs .EQ. 1) then             if (iobcs .EQ. 1) then
253                         OBEt(j,k,bi,bj) = OBEt (j,k,bi,bj)              OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)
254       &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)       &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
255       &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
256                         OBEt(j,k,bi,bj) = OBEt(j,k,bi,bj)              OBEt(j,k,bi,bj) = OBEt(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                         OBEs(j,k,bi,bj) = OBEs (j,k,bi,bj)              OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)
260       &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)       &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
261       &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
262                         OBEs(j,k,bi,bj) = OBEs(j,k,bi,bj)              OBEs(j,k,bi,bj) = OBEs(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                         OBEu(j,k,bi,bj) = OBEu (j,k,bi,bj)              OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)
266       &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)       &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
267       &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
268                         OBEu(j,k,bi,bj) = OBEu(j,k,bi,bj)              OBEu(j,k,bi,bj) = OBEu(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                         OBEv(j,k,bi,bj) = OBEv (j,k,bi,bj)              OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)
272       &                 + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)       &           + obcsefac            *xx_obcse0(j,k,bi,bj,iobcs)
273       &                 + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)       &           + (1. _d 0 - obcsefac)*xx_obcse1(j,k,bi,bj,iobcs)
274                         OBEv(j,k,bi,bj) = OBEv(j,k,bi,bj)              OBEv(j,k,bi,bj) = OBEv(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_OBCSE_CONTROL undefined */  
   
 c     == routine arguments ==  
   
       _RL     mytime  
       integer myiter  
       integer mythid  
   
 c--   CPP flag ALLOW_OBCSE_CONTROL undefined.  
   
285  #endif /* ALLOW_OBCSE_CONTROL */  #endif /* ALLOW_OBCSE_CONTROL */
286    
287          return
288        end        end
   
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.22