/[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.6 by mlosch, Fri Dec 3 00:48:57 2004 UTC revision 1.11 by mlosch, Mon Mar 14 17:08:00 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
# 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
# Line 106  c--   Now, read the control vector. Line 109  c--   Now, read the control vector.
109        ladinit      = .false.        ladinit      = .false.
110    
111        if (optimcycle .ge. 0) then        if (optimcycle .ge. 0) then
112          ilobcsw=ilnblnk( xx_obcsw_file )         ilobcsw=ilnblnk( xx_obcsw_file )
113          write(fnameobcsw(1:80),'(2a,i10.10)')         write(fnameobcsw(1:80),'(2a,i10.10)')
114       &       xx_obcsw_file(1:ilobcsw), '.', optimcycle       &      xx_obcsw_file(1:ilobcsw), '.', optimcycle
115        endif        endif
116    
117  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 122  c--   Get the counters, flags, and the i
122       I                   mytime, myiter, mythid )       I                   mytime, myiter, mythid )
123    
124        do iobcs = 1,nobcs        do iobcs = 1,nobcs
125          if ( obcswfirst ) then         if ( obcswfirst ) then
126            call active_read_yz_loc( fnameobcsw, tmpfldyz,          call active_read_yz( fnameobcsw, tmpfldyz,
127       &                         (obcswcount0-1)*nobcs+iobcs,       &                       (obcswcount0-1)*nobcs+iobcs,
128       &                         doglobalread, ladinit, optimcycle,       &                       doglobalread, ladinit, optimcycle,
129       &                         mythid, xx_obcsw_dummy )       &                       mythid, xx_obcsw_dummy )
130    
131  #ifdef ALLOW_CTRL_OBCS_BALANCE          do bj = jtlo,jthi
132             do bi = itlo,ithi
133            if ( optimcycle .gt. 0) then                      do k = 1,nr
134              if (iobcs .eq. 3) then             do j = jmin,jmax
135  cgg         Special attention is needed for the normal velocity.              xx_obcsw1(j,k,bi,bj,iobcs)  = tmpfldyz (j,k,bi,bj)
136  cgg         For the north, this is the v velocity, iobcs = 4.  cgg   &                                        *   maskyz (j,k,bi,bj)
137  cgg         This is done on a columnwise basis here.             enddo
               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  
           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  
             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  
             enddo  
138            enddo            enddo
139             enddo
140            enddo
141           endif
142    
143            call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid)         if ( (obcswfirst) .or. (obcswchanged)) then
   
           do bj = jtlo,jthi  
             do bi = itlo,ithi  
               do k = 1,nr  
                 do j = jmin,jmax  
                   xx_obcsw0(j,k,bi,bj,iobcs) = tmpfldyz2(j,k,bi,bj)  
                 enddo  
               enddo  
             enddo  
           enddo  
144    
145            call active_read_yz_loc( fnameobcsw, tmpfldyz,          do bj = jtlo,jthi
146       &                         (obcswcount1-1)*nobcs+iobcs,           do bi = itlo,ithi
147       &                         doglobalread, ladinit, optimcycle,            do k = 1,nr
148       &                         mythid, xx_obcsw_dummy )             do j = jmin,jmax
149                xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs)
150  #ifdef ALLOW_CTRL_OBCS_BALANCE              tmpfldyz (j,k,bi,bj)       = 0. _d 0
151               enddo
           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_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  
           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  
             enddo  
152            enddo            enddo
153          endif           enddo
154            enddo
155    
156  c--     Add control to model variable.          call active_read_yz( fnameobcsw, tmpfldyz,
157          do bj = jtlo, jthi       &                       (obcswcount1-1)*nobcs+iobcs,
158             do bi = itlo, ithi       &                       doglobalread, ladinit, optimcycle,
159  c--        Calculate mask for tracer cells (0 => land, 1 => water).       &                       mythid, xx_obcsw_dummy )
160                do k = 1,nr  
161                   do j = 1,sny          do bj = jtlo,jthi
162                      i = OB_Iw(j,bi,bj)           do bi = itlo,ithi
163                      if (iobcs .EQ. 1) then            do k = 1,nr
164                         OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)             do j = jmin,jmax
165       &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)              xx_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj)
166       &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  cgg   &                                        *   maskyz (j,k,bi,bj)
                        OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 2) then  
                        OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 3) then  
                        OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)  
      &                      *maskW(i+ip1,j,k,bi,bj)  
                     else if (iobcs .EQ. 4) then  
                        OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)  
      &                 + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)  
      &                 + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)  
                        OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)  
      &                      *maskS(i,j,k,bi,bj)  
                     endif  
                  enddo  
               enddo  
167             enddo             enddo
168              enddo
169             enddo
170          enddo          enddo
171           endif
172    
173    c--   Add control to model variable.
174           do bj = jtlo, jthi
175            do bi = itlo, ithi
176    c--   Calculate mask for tracer cells (0 => land, 1 => water).
177             do k = 1,nr
178              do j = 1,sny
179               i = OB_Iw(j,bi,bj)
180               if (iobcs .EQ. 1) then
181                OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj)
182         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
183         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
184                OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj)
185         &           *maskW(i+ip1,j,k,bi,bj)
186               else if (iobcs .EQ. 2) then
187                OBWs(j,k,bi,bj) = OBWs (j,k,bi,bj)
188         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
189         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
190                OBWs(j,k,bi,bj) = OBWs(j,k,bi,bj)
191         &           *maskW(i+ip1,j,k,bi,bj)
192               else if (iobcs .EQ. 3) then
193                OBWu(j,k,bi,bj) = OBWu (j,k,bi,bj)
194         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
195         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
196                OBWu(j,k,bi,bj) = OBWu(j,k,bi,bj)
197         &           *maskW(i+ip1,j,k,bi,bj)
198               else if (iobcs .EQ. 4) then
199                OBWv(j,k,bi,bj) = OBWv (j,k,bi,bj)
200         &           + obcswfac            *xx_obcsw0(j,k,bi,bj,iobcs)
201         &           + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs)
202                OBWv(j,k,bi,bj) = OBWv(j,k,bi,bj)
203         &           *maskS(i,j,k,bi,bj)
204               endif
205              enddo
206             enddo
207            enddo
208           enddo
209          
210  C--   End over iobcs loop  C--   End over iobcs loop
211        enddo        enddo
212    
Legend:
Removed from v.1.6  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.22