--- MITgcm/pkg/ctrl/ctrl_getobcsw.F 2002/02/05 20:23:58 1.1.2.1 +++ MITgcm/pkg/ctrl/ctrl_getobcsw.F 2002/05/30 19:56:44 1.1.2.2 @@ -15,7 +15,7 @@ c SUBROUTINE ctrl_getobcsw c ================================================================== c -c o Get norhtern obc of the control vector and add it +c o Get western obc of the control vector and add it c to dyn. fields c c started: heimbach@mit.edu, 29-Aug-2001 @@ -56,6 +56,7 @@ integer imin,imax integer ilobcsw integer iobcs + integer ip1 _RL dummy _RL obcswfac @@ -64,14 +65,19 @@ integer obcswcount0 integer obcswcount1 - _RL maskyz (1-oly:sny+oly,nr,nsx,nsy) +cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy) logical doglobalread logical ladinit character*(80) fnameobcsw - +#ifdef BALANCE_CONTROL_VOLFLUX +cgg( Variables for balancing volume flux. + _RL ubaro + _RL ubarocount +cgg) +#endif c == external functions == @@ -89,6 +95,14 @@ jmax = sny+oly imin = 1-olx imax = snx+olx + ip1 = 1 + +#ifdef BALANCE_CONTROL_VOLFLUX +cgg( Initialize variables for balancing volume flux. + ubaro = 0.d0 + ubarocount = 0.d0 +cgg) +#endif c-- Now, read the control vector. doglobalread = .false. @@ -111,41 +125,178 @@ I mytime, myiter, mythid ) do iobcs = 1,nobcs + +cgg if ( (obcswfirst) .or. (obcswchanged) ) then +cgg call active_read_yz( 'maskobcsw', maskyz, +cgg & iobcs, +cgg & doglobalread, ladinit, 0, +cgg & mythid, dummy ) +cgg endif + + if ( obcswfirst ) then + call active_read_yz( fnameobcsw, tmpfldyz, + & (obcswcount0-1)*nobcs+iobcs, + & doglobalread, ladinit, optimcycle, + & mythid, xx_obcsw_dummy ) + +#ifdef BALANCE_CONTROL_VOLFLUX + if ( optimcycle .gt. 0) then + if (iobcs .eq. 3) then +cgg( Make sure that the xx velocity field has a balanced net volume flux. +cgg Find the barotropic flow normal to the boundary. +cgg For the north, this is the v velocity, iobcs = 4. + do bj = jtlo,jthi + do bi = itlo, ithi + do j = jmin,jmax + i = OB_Iw(j,bi,bj) + ubaro = 0.d0 + ubarocount = 0.d0 + do k = 1,nr +cgg( If cells are not full, this should be modified with hFac. + ubaro = tmpfldyz(j,k,bi,bj)*maskW(i+ip1,j,k,bi,bj) + & * delZ(k) + ubaro + ubarocount = maskW(i+ip1,j,k,bi,bj) + & * delZ(k) +ubarocount + enddo + if (ubarocount .eq. 0.) then + ubaro = 0. + else + ubaro = ubaro / ubarocount + endif + do k = 1,nr + tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro + enddo + enddo + enddo + enddo + endif + endif +cgg) +#endif +#ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL + if (optimcycle .gt. 0) then + if ( iobcs .eq. 3) then + do bj = jtlo,jthi + do bi = itlo, ithi + do k = 1,Nr + do j = jmin,jmax + i = OB_Iw(j,bi,bj) + tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) + & - shiftvel(1) *maskW(i+ip1,j,k,bi,bj) + enddo + enddo + enddo + enddo + endif + endif +#endif - call active_read_yz( 'maskobcsw', maskyz, - & iobcs, - & doglobalread, ladinit, 0, - & mythid, dummy ) - - call active_read_yz( fnameobcsw, tmpfldyz, - & (obcswcount0-1)*nobcs+iobcs, - & doglobalread, ladinit, optimcycle, - & mythid, xx_obcsw_dummy ) + 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) + enddo + enddo + enddo + enddo + endif - do bj = jtlo,jthi - do bi = itlo,ithi - do k = 1,nr - do j = jmin,jmax - yz_obcs0(j,k,bi,bj) = tmpfldyz (j,k,bi,bj) + 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 enddo enddo - enddo - call active_read_yz( fnameobcsw, tmpfldyz, - & (obcswcount1-1)*nobcs+iobcs, - & doglobalread, ladinit, optimcycle, - & mythid, xx_obcsw_dummy ) + call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid) - do bj = jtlo,jthi - do bi = itlo,ithi - do k = 1,nr - do j = jmin,jmax - yz_obcs1 (j,k,bi,bj) = tmpfldyz (j,k,bi,bj) + 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 + + call active_read_yz( fnameobcsw, tmpfldyz, + & (obcswcount1-1)*nobcs+iobcs, + & doglobalread, ladinit, optimcycle, + & mythid, xx_obcsw_dummy ) + +#ifdef BALANCE_CONTROL_VOLFLUX + if (optimcycle .gt. 0) then + if (iobcs .eq. 3) then +cgg( Make sure that the xx velocity field has a balanced net volume flux. +cgg Find the barotropic flow normal to the boundary. +cgg For the north, this is the v velocity, iobcs = 4. + do bj = jtlo,jthi + do bi = itlo, ithi + do j = jmin,jmax + i = OB_Iw(j,bi,bj) + ubaro = 0.d0 + ubarocount = 0.d0 + do k = 1,nr +cgg( If cells are not full, this should be modified with hFac. + ubaro = tmpfldyz(j,k,bi,bj)*maskw(i+ip1,j,k,bi,bj) + & * delZ(k) + ubaro + ubarocount = maskW(i+ip1,j,k,bi,bj) + & * delZ(k) + ubarocount + enddo + if (ubarocount .eq. 0.) then + ubaro = 0. + else + ubaro = ubaro / ubarocount + endif + do k = 1,nr + tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) - ubaro + enddo + enddo + enddo + enddo + endif + endif +cgg) +#endif +#ifdef BALANCE_CONTROL_VOLFLUX_GLOBAL + if (optimcycle .gt. 0) then + if ( iobcs .eq. 3) then + do bj = jtlo,jthi + do bi = itlo, ithi + do k = 1,Nr + do j = jmin,jmax + i = OB_Iw(j,bi,bj) + tmpfldyz(j,k,bi,bj) = tmpfldyz(j,k,bi,bj) + & - shiftvel(2) *maskW(i+ip1,j,k,bi,bj) + enddo + enddo + enddo + enddo + endif + 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) + enddo enddo enddo enddo - enddo + endif c-- Add control to model variable. do bj = jtlo,jthi @@ -153,30 +304,31 @@ c-- Calculate mask for tracer cells (0 => land, 1 => water). do k = 1,nr do j = 1,sny + i = OB_Iw(j,bi,bj) if (iobcs .EQ. 1) then OBWt(j,k,bi,bj) = OBWt (j,k,bi,bj) - & + obcswfac *yz_obcs0(j,k,bi,bj) - & + (1. _d 0 - obcswfac)*yz_obcs1(j,k,bi,bj) + & + obcswfac *xx_obcsw0(j,k,bi,bj,iobcs) + & + (1. _d 0 - obcswfac)*xx_obcsw1(j,k,bi,bj,iobcs) OBWt(j,k,bi,bj) = OBWt(j,k,bi,bj) - & *maskyz(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 *yz_obcs0(j,k,bi,bj) - & + (1. _d 0 - obcswfac)*yz_obcs1(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) - & *maskyz(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 *yz_obcs0(j,k,bi,bj) - & + (1. _d 0 - obcswfac)*yz_obcs1(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) - & *maskyz(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 *yz_obcs0(j,k,bi,bj) - & + (1. _d 0 - obcswfac)*yz_obcs1(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) - & *maskyz(j,k,bi,bj) + & *maskS(i,j,k,bi,bj) endif enddo enddo