C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/SOSE/BoxAdj/code_ad/Attic/ctrl_getobcsw.F,v 1.5 2011/01/19 21:08:32 mmazloff Exp $ C $Name: $ #include "CTRL_CPPOPTIONS.h" #ifdef ALLOW_OBCS # include "OBCS_OPTIONS.h" #endif subroutine ctrl_getobcsw( I mytime, I myiter, I mythid & ) c ================================================================== c SUBROUTINE ctrl_getobcsw c ================================================================== c 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 c c modified: gebbie@mit.edu, 18-Mar-2003 c ================================================================== c SUBROUTINE ctrl_getobcsw c ================================================================== implicit none #ifdef ALLOW_OBCSW_CONTROL c == global variables == #include "EEPARAMS.h" #include "SIZE.h" #include "PARAMS.h" #include "GRID.h" #include "OBCS.h" #include "ctrl.h" #include "ctrl_dummy.h" #include "optim.h" c == routine arguments == _RL mytime integer myiter integer mythid c == local variables == integer bi,bj integer i,j,k integer itlo,ithi integer jtlo,jthi integer jmin,jmax integer imin,imax integer ilobcsw integer iobcs integer ip1 _RL dummy _RL obcswfac logical obcswfirst logical obcswchanged integer obcswcount0 integer obcswcount1 integer nk,nz _RL tmpz (nr,nsx,nsy) _RL stmp cgg _RL maskyz (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) c == external functions == integer ilnblnk external ilnblnk c == end of interface == jtlo = mybylo(mythid) jthi = mybyhi(mythid) itlo = mybxlo(mythid) ithi = mybxhi(mythid) jmin = 1-oly jmax = sny+oly imin = 1-olx 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 endif c-- Get the counters, flags, and the interpolation factor. call ctrl_get_gen_rec( I xx_obcswstartdate, xx_obcswperiod, O obcswfac, obcswfirst, obcswchanged, O obcswcount0,obcswcount1, I mytime, myiter, mythid ) do iobcs = 1,nobcs if ( obcswfirst ) then call active_read_yz( 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 cih If open boundary. if ( OB_Iw(j,bi,bj) .ne. 0. ) then i = OB_Iw(J,bi,bj) #ifdef ALLOW_OBCS_CONTROL_MODES cih Determine number of open vertical layers. nz = 0 do k = 1,Nr nz = nz + maskW(i+ip1,j,k,bi,bj) end do cih Compute absolute velocities from the barotropic-baroclinic modes. #ifdef ALLOW_CTRL_OBCS_BALANCE CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are c compatible - to ensure volume conservation can just set barotropic c mode amplitude to 0 c however this means no inflow at every horizontal location.... do k = 1,nr tmpfldyz(k,1,bi,bj)= 0. end do #endif 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 end do do k = 1,Nr tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) & *recip_hFacW(i+ip1,j,k,bi,bj) end do #elif defined (ALLOW_CTRL_OBCS_BALANCE) 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) #endif cih End open boundary. end if 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 cih If open boundary. if ( OB_Iw(j,bi,bj) .ne. 0. ) then i = OB_Iw(J,bi,bj) #ifdef ALLOW_OBCS_CONTROL_MODES cih Determine number of open vertical layers. nz = 0 do k = 1,Nr nz = nz + maskS(i,j,k,bi,bj) end do cih Compute absolute velocities from the barotropic-baroclinic modes. #ifdef ALLOW_CTRL_OBCS_BALANCE CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are c compatible - to ensure volume conservation can just set barotropic c mode amplitude to 0 c however this means no inflow at every horizontal location.... do k = 1,nr tmpfldyz(k,1,bi,bj)= 0. end do #endif 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 end do do k = 1,Nr tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) & *recip_hFacS(i,j,k,bi,bj) end do #elif defined (ALLOW_CTRL_OBCS_BALANCE) 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) #endif cih End if iobcs = 4. end if 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 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( 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 cih If open boundary. if ( OB_Iw(j,bi,bj) .ne. 0. ) then i = OB_Iw(J,bi,bj) #ifdef ALLOW_OBCS_CONTROL_MODES cih Determine number of open vertical layers. nz = 0 do k = 1,Nr nz = nz + maskW(i+ip1,j,k,bi,bj) end do cih Compute absolute velocities from the barotropic-baroclinic modes. #ifdef ALLOW_CTRL_OBCS_BALANCE CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are c compatible - to ensure volume conservation can just set barotropic c mode amplitude to 0 c however this means no inflow at every horizontal location.... do k = 1,nr tmpfldyz(k,1,bi,bj)= 0. end do #endif 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 end do do k = 1,Nr tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) & *recip_hFacW(i+ip1,j,k,bi,bj) end do #elif defined (ALLOW_CTRL_OBCS_BALANCE) 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) #endif cih End open boundary. end if 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 cih If open boundary. if ( OB_Iw(j,bi,bj) .ne. 0. ) then i = OB_Iw(J,bi,bj) #ifdef ALLOW_OBCS_CONTROL_MODES cih Determine number of open vertical layers. nz = 0 do k = 1,Nr nz = nz + maskS(i,j,k,bi,bj) end do cih Compute absolute velocities from the barotropic-baroclinic modes. #ifdef ALLOW_CTRL_OBCS_BALANCE CMM not sure if ALLOW_OBCS_CONTROL_MODES and ALLOW_CTRL_OBCS_BALANCE are c compatible - to ensure volume conservation can just set barotropic c mode amplitude to 0 c however this means no inflow at every horizontal location.... do k = 1,nr tmpfldyz(k,1,bi,bj)= 0. end do #endif 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 end do do k = 1,Nr tmpfldyz(j,k,bi,bj) = tmpz(k,bi,bj) & *recip_hFacS(i,j,k,bi,bj) end do #elif defined (ALLOW_CTRL_OBCS_BALANCE) 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) #endif cih End if iobcs = 4. end if 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 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 */ end