--- MITgcm/pkg/ctrl/ctrl_getobcsw.F 2003/10/26 00:58:03 1.5 +++ MITgcm/pkg/ctrl/ctrl_getobcsw.F 2011/03/07 09:24:10 1.10 @@ -1,3 +1,5 @@ +C $Header: /home/ubuntu/mnt/e9_copy/MITgcm/pkg/ctrl/ctrl_getobcsw.F,v 1.10 2011/03/07 09:24:10 mlosch Exp $ +C $Name: $ #include "CTRL_CPPOPTIONS.h" #ifdef ALLOW_OBCS @@ -67,6 +69,7 @@ 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 @@ -107,7 +110,7 @@ if (optimcycle .ge. 0) then ilobcsw=ilnblnk( xx_obcsw_file ) - write(fnameobcsw(1:80),'(2a,i10.10)') + write(fnameobcsw(1:80),'(2a,i10.10)') & xx_obcsw_file(1:ilobcsw), '.', optimcycle endif @@ -118,14 +121,18 @@ O obcswcount0,obcswcount1, I mytime, myiter, mythid ) +CML print *,'ml-getobcs: ',myIter,obcswfirst,obcswchanged, +CML & obcswcount0,obcswcount1,obcswfac do iobcs = 1,nobcs if ( obcswfirst ) then - call active_read_yz_loc( fnameobcsw, tmpfldyz, + call active_read_yz( fnameobcsw, tmpfldyz, & (obcswcount0-1)*nobcs+iobcs, & doglobalread, ladinit, optimcycle, & mythid, xx_obcsw_dummy ) - if ( optimcycle .gt. 0) then +#ifdef ALLOW_CTRL_OBCS_BALANCE + + 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. @@ -142,20 +149,20 @@ 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 +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.. +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. +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) + tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj) & - utop / delR(1) enddo enddo @@ -177,20 +184,20 @@ 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 +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.. +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. +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) + tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj) & - utop / delR(1) enddo enddo @@ -198,6 +205,8 @@ endif endif +#endif /* ALLOW_CTRL_OBCS_BALANCE */ + do bj = jtlo,jthi do bi = itlo,ithi do k = 1,nr @@ -211,38 +220,26 @@ 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 - enddo - - call exf_swapffields_yz( tmpfldyz2, tmpfldyz, mythid) - - 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 + 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, + call active_read_yz( fnameobcsw, tmpfldyz, & (obcswcount1-1)*nobcs+iobcs, & doglobalread, ladinit, optimcycle, & mythid, xx_obcsw_dummy ) - if ( optimcycle .gt. 0) then +#ifdef ALLOW_CTRL_OBCS_BALANCE + + 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. @@ -259,20 +256,20 @@ 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 +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.. +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. +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) + tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj) & - utop / delR(1) enddo enddo @@ -294,20 +291,20 @@ 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 +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.. +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. +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) + tmpfldyz(j,1,bi,bj) = tmpfldyz(j,1,bi,bj) & - utop / delR(1) enddo enddo @@ -315,6 +312,8 @@ endif endif +#endif /* ALLOW_CTRL_OBCS_BALANCE */ + do bj = jtlo,jthi do bi = itlo,ithi do k = 1,nr @@ -327,41 +326,41 @@ enddo endif -c-- Add control to model variable. +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 + 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