/[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.10 by mlosch, Mon Mar 7 09:24:10 2011 UTC revision 1.17 by gforget, Thu Oct 9 00:49:26 2014 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 "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 76  cgg      _RL maskyz   (1-oly:sny+oly,nr, Line 77  cgg      _RL maskyz   (1-oly:sny+oly,nr,
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 99  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 121  c--   Get the counters, flags, and the i Line 118  c--   Get the counters, flags, and the i
118       O                   obcswcount0,obcswcount1,       O                   obcswcount0,obcswcount1,
119       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
120    
 CML      print *,'ml-getobcs: ',myIter,obcswfirst,obcswchanged,  
 CML     &     obcswcount0,obcswcount1,obcswfac  
121        do iobcs = 1,nobcs        do iobcs = 1,nobcs
122          if ( obcswfirst ) then         if ( obcswfirst ) then
123            call active_read_yz( 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  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
129             do bi = itlo,ithi
130            if ( optimcycle .gt. 0) then  #ifdef ALLOW_OBCS_CONTROL_MODES
131              if (iobcs .eq. 3) then            if (iobcs .gt. 2) then
132  cgg         Special attention is needed for the normal velocity.             do j = jmin,jmax
133  cgg         For the north, this is the v velocity, iobcs = 4.              i = OB_Iw(j,bi,bj)
134  cgg         This is done on a columnwise basis here.              IF ( i.EQ.OB_indexNone ) i = 1
135                do bj = jtlo,jthi  cih    Determine number of open vertical layers.
136                  do bi = itlo, ithi              nz = 0
137                    do j = jmin,jmax              do k = 1,Nr
138                      i = OB_Iw(J,bi,bj)                if (iobcs .eq. 3) then
139                    nz = nz + maskW(i+ip1,j,k,bi,bj)
140  cgg         The barotropic velocity is stored in the level 1.                else
141                      ubaro = tmpfldyz(j,1,bi,bj)                  nz = nz + maskS(i,j,k,bi,bj)
142                      tmpfldyz(j,1,bi,bj) = 0.d0                endif
143                      utop = 0.d0              end do
144    cih    Compute absolute velocities from the barotropic-baroclinic modes.
145                      do k = 1,Nr              do k = 1,Nr
146  cgg    If cells are not full, this should be modified with hFac.               if (k.le.nz) then
147  cgg                stmp = 0.
148  cgg    The xx field (tmpfldyz) does not contain the velocity at the                do nk = 1,nz
149  cgg    surface level. This velocity is not independent; it must                 stmp = stmp +
150  cgg    exactly balance the volume flux, since we are dealing with       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
151  cgg    the baroclinic velocity structure..                end do
152                        utop = tmpfldyz(j,k,bi,bj)*                 tmpz(k,bi,bj) = stmp
153       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop               else
154  cgg    Add the barotropic velocity component.                tmpz(k,bi,bj) = 0.
155                        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 (tmpfldyz) 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  
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            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) = xx_obcsw1(j,k,bi,bj,iobcs)              xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
186                tmpfldyz (j,k,bi,bj)       = 0. _d 0              tmpfldyz (j,k,bi,bj)       = 0. _d 0
              enddo  
             enddo  
187             enddo             enddo
188            enddo            enddo
189             enddo
190            enddo
191    
192            call active_read_yz( 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  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
198             do bi = itlo,ithi
199            if ( optimcycle .gt. 0) then  #ifdef ALLOW_OBCS_CONTROL_MODES
200              if (iobcs .eq. 3) then            if (iobcs .gt. 2) then
201  cgg         Special attention is needed for the normal velocity.             do j = jmin,jmax
202  cgg         For the north, this is the v velocity, iobcs = 4.              i = OB_Iw(j,bi,bj)
203  cgg         This is done on a columnwise basis here.              IF ( i.EQ.OB_indexNone ) i = 1
204                do bj = jtlo,jthi  cih    Determine number of open vertical layers.
205                  do bi = itlo, ithi              nz = 0
206                    do j = jmin,jmax              do k = 1,Nr
207                      i = OB_Iw(J,bi,bj)                if (iobcs .eq. 3) then
208                    nz = nz + maskW(i+ip1,j,k,bi,bj)
209  cgg         The barotropic velocity is stored in the level 1.                else
210                      ubaro = tmpfldyz(j,1,bi,bj)                  nz = nz + maskS(i,j,k,bi,bj)
211                      tmpfldyz(j,1,bi,bj) = 0.d0                endif
212                      utop = 0.d0              end do
213    cih    Compute absolute velocities from the barotropic-baroclinic modes.
214                      do k = 1,Nr              do k = 1,Nr
215  cgg    If cells are not full, this should be modified with hFac.               if (k.le.nz) then
216  cgg                stmp = 0.
217  cgg    The xx field (tmpfldyz) does not contain the velocity at the                do nk = 1,nz
218  cgg    surface level. This velocity is not independent; it must                 stmp = stmp +
219  cgg    exactly balance the volume flux, since we are dealing with       &         modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj)
220  cgg    the baroclinic velocity structure..                end do
221                        utop = tmpfldyz(j,k,bi,bj)*                 tmpz(k,bi,bj) = stmp
222       &                maskW(i+ip1,j,k,bi,bj) * delR(k) + utop               else
223  cgg    Add the barotropic velocity component.                 tmpz(k,bi,bj) = 0.
224                        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 (tmpfldyz) 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  
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.10  
changed lines
  Added in v.1.17
  ViewVC Help
Powered by ViewVC 1.1.22