--- MITgcm/pkg/ctrl/ctrl_getobcsw.F 2003/07/18 21:10:16 1.3 +++ MITgcm/pkg/ctrl/ctrl_getobcsw.F 2012/07/31 16:05:56 1.14 @@ -1,10 +1,11 @@ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ctrl/ctrl_getobcsw.F,v 1.14 2012/07/31 16:05:56 heimbach Exp $ +C $Name: $ #include "CTRL_CPPOPTIONS.h" #ifdef ALLOW_OBCS # include "OBCS_OPTIONS.h" #endif - subroutine ctrl_getobcsw( I mytime, I myiter, @@ -27,26 +28,27 @@ implicit none -#ifdef ALLOW_OBCSW_CONTROL - c == global variables == - +#ifdef ALLOW_OBCSW_CONTROL #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" -#include "OBCS.h" - +c#include "OBCS_PARAMS.h" +#include "OBCS_GRID.h" +#include "OBCS_FIELDS.h" +#include "CTRL_SIZE.h" #include "ctrl.h" #include "ctrl_dummy.h" #include "optim.h" +#endif /* ALLOW_OBCSW_CONTROL */ c == routine arguments == - _RL mytime integer myiter integer mythid +#ifdef ALLOW_OBCSW_CONTROL c == local variables == integer bi,bj @@ -67,16 +69,18 @@ integer obcswcount1 cgg _RL maskyz (1-oly:sny+oly,nr,nsx,nsy) + _RL tmpfldyz (1-oly:sny+oly,nr,nsx,nsy) logical doglobalread logical ladinit character*(80) fnameobcsw -cgg( Variables for splitting barotropic/baroclinic vels. - _RL ubaro - _RL utop -cgg) +#ifdef ALLOW_OBCS_CONTROL_MODES + integer nk,nz + _RL tmpz (nr,nsx,nsy) + _RL stmp +#endif c == external functions == @@ -96,19 +100,14 @@ imax = snx+olx ip1 = 1 -cgg( Initialize variables for balancing volume flux. - ubaro = 0.d0 - utop = 0.d0 -cgg) - c-- Now, read the control vector. doglobalread = .false. ladinit = .false. if (optimcycle .ge. 0) then - ilobcsw=ilnblnk( xx_obcsw_file ) - write(fnameobcsw(1:80),'(2a,i10.10)') - & xx_obcsw_file(1:ilobcsw), '.', optimcycle + ilobcsw=ilnblnk( xx_obcsw_file ) + write(fnameobcsw(1:80),'(2a,i10.10)') + & xx_obcsw_file(1:ilobcsw), '.', optimcycle endif c-- Get the counters, flags, and the interpolation factor. @@ -119,265 +118,171 @@ I mytime, myiter, mythid ) do iobcs = 1,nobcs - if ( obcswfirst ) then - call active_read_yz_loc( fnameobcsw, tmpfldyz, - & (obcswcount0-1)*nobcs+iobcs, - & doglobalread, ladinit, optimcycle, - & mythid, xx_obcsw_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_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) * 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_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) * 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_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 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 + if ( obcswfirst ) then + 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 +#ifdef ALLOW_OBCS_CONTROL_MODES + if (iobcs .gt. 2) then + do j = jmin,jmax + i = OB_Iw(j,bi,bj) +cih Determine number of open vertical layers. + nz = 0 + do k = 1,Nr + if (iobcs .eq. 3) then + nz = nz + maskW(i+ip1,j,k,bi,bj) + else + nz = nz + maskS(i,j,k,bi,bj) + endif + end do +cih Compute absolute velocities from the barotropic-baroclinic modes. + do k = 1,Nr + if (k.le.nz) then + stmp = 0. + do nk = 1,nz + stmp = stmp + + & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj) + end do + tmpz(k,bi,bj) = stmp + else + tmpz(k,bi,bj) = 0. + end if enddo + do k = 1,Nr + if (iobcs .eq. 3) then + tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) + & *recip_hFacW(i+ip1,j,k,bi,bj) + else + tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) + & *recip_hFacS(i,j,k,bi,bj) + endif + end do + enddo + endif +#endif + 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 - 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 + do bj = jtlo,jthi + do bi = itlo,ithi + do k = 1,nr + do j = jmin,jmax + xx_obcsw0(j,k,bi,bj,iobcs) = xx_obcsw1(j,k,bi,bj,iobcs) + tmpfldyz (j,k,bi,bj) = 0. _d 0 + enddo enddo + enddo + enddo - call active_read_yz_loc( fnameobcsw, tmpfldyz, - & (obcswcount1-1)*nobcs+iobcs, - & doglobalread, ladinit, optimcycle, - & mythid, xx_obcsw_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_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) * 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_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) * 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_obcsw1 (j,k,bi,bj,iobcs) = tmpfldyz (j,k,bi,bj) -cgg & * maskyz (j,k,bi,bj) - enddo - enddo + call active_read_yz( fnameobcsw, tmpfldyz, + & (obcswcount1-1)*nobcs+iobcs, + & doglobalread, ladinit, optimcycle, + & mythid, xx_obcsw_dummy ) + + do bj = jtlo,jthi + do bi = itlo,ithi +#ifdef ALLOW_OBCS_CONTROL_MODES + if (iobcs .gt. 2) then + do j = jmin,jmax + i = OB_Iw(j,bi,bj) +cih Determine number of open vertical layers. + nz = 0 + do k = 1,Nr + if (iobcs .eq. 3) then + nz = nz + maskW(i+ip1,j,k,bi,bj) + else + nz = nz + maskS(i,j,k,bi,bj) + endif + end do +cih Compute absolute velocities from the barotropic-baroclinic modes. + do k = 1,Nr + if (k.le.nz) then + stmp = 0. + do nk = 1,nz + stmp = stmp + + & modesv(k,nk,nz)*tmpfldyz(j,nk,bi,bj) + end do + tmpz(k,bi,bj) = stmp + else + tmpz(k,bi,bj) = 0. + end if enddo + do k = 1,Nr + if (iobcs .eq. 3) then + tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) + & *recip_hFacW(i+ip1,j,k,bi,bj) + else + tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) + & *recip_hFacS(i,j,k,bi,bj) + endif + end do + enddo + endif +#endif + 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 - endif + enddo + enddo + endif -c-- Add control to model variable. - do bj = jtlo, jthi - do bi = itlo, ithi -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 *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) - & *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 - enddo +c-- Add control to model variable. + do bj = jtlo, jthi + do bi = itlo, ithi +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 *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) + & *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 enddo + enddo C-- End over iobcs loop enddo -#else /* ALLOW_OBCSW_CONTROL undefined */ - -c == routine arguments == - - _RL mytime - integer myiter - integer mythid - -c-- CPP flag ALLOW_OBCSW_CONTROL undefined. - #endif /* ALLOW_OBCSW_CONTROL */ + return end -