#if !defined usrdiagnostic | !defined nesttime subroutine clinic(j) #else subroutine clinic(j,dcpg,dco) #endif c c======================================================================= c === c CLINIC computes, for one row, the internal mode component of === c the U- and V-velocities, as well as the vorticity === c driving function for use by RELAX later in determining === c the external modes, where: === c === c J = the row number === c === c======================================================================= c c----------------------------------------------------------------------- c Define global data. c----------------------------------------------------------------------- c #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #include #if defined bndy_rlx | defined cstfrc | defined botfrc # include #endif #if defined ext_tide & defined tide_zero # include #endif c c----------------------------------------------------------------------- c Define local data. c----------------------------------------------------------------------- c integer i,j,k #ifndef bioadjloc integer iend,im1,istr,kz,ll,sofar #endif #ifndef islands integer l #endif #ifndef bioadjloc FLOAT * boxfac,boxvol,detmr,diag1,diag2,fx FLOAT * ueng(imt,km),umean(imt),veng(imt,km),vmean(imt) #endif #ifndef explicitvmix FLOAT * r2dtuv,twodt(km) #endif #if defined usrdiagnostic & defined nesttime FLOAT * dcpg(2),dco(2),ddco(2) #endif c c----------------------------------------------------------------------- c Begin executable code. c----------------------------------------------------------------------- c c======================================================================= c Begin introductory section, preparing various =================== c arrays for the computation of the internal modes =================== c======================================================================= c c----------------------------------------------------------------------- c Save internal mode (baroclinic) velocities at rows J-1, J, and J+1. c----------------------------------------------------------------------- c do 100 k=1,km do 100 i=1,imt uclm(i,k)=ucl(i,k) vclm(i,k)=vcl(i,k) ucl(i,k)=uclp(i,k) vcl(i,k)=vclp(i,k) uclb(i,k)=uclbp(i,k) vclb(i,k)=vclbp(i,k) uclp(i,k)=up(i,k) vclp(i,k)=vp(i,k) uclbp(i,k)=ubp(i,k) vclbp(i,k)=vbp(i,k) 100 continue c c----------------------------------------------------------------------- c Save ubarm and ubarbm for bnd conditions. Since ubar is reset to ubar c at row j, at the call to clinic.F, it contains ubar at row j-1 c----------------------------------------------------------------------- c if (j.gt.2) then do 105 i=1,imt ubarm(i)=ubar(i) vbarm(i)=vbar(i) ubarbm(i)=ubarb(i) vbarbm(i)=vbarb(i) 105 continue endif #ifndef bioadjloc c c----------------------------------------------------------------------- c Calculate advective coefficients FUW at west face of U,V box and c FVN at north face of U,V box by combining external an internal c modes with a grid factor. c----------------------------------------------------------------------- c c 1st, calculate external mode U at west face of U,V box c and V at north face of U,V box. c do 110 i=1,imtm1 im1=max(1,i-1) ubar(i)=-(p(i ,j+1)-p(i,j ))*min(c1,fkmu(im1),fkmu (i))* * dyur(j)/(p5*(hdv(im1,j)+hdv(i,j ))) vbar(i)= (p(i+1,j+1)-p(i,j+1))*min(c1,fkmu(i ),fkmup(i))* * dxur(i)/(p5*(hdv(i ,j)+hdv(i,j+1)))*cstr(j+1) 110 continue #ifdef cyclic ubar(1 )=ubar(imtm1) vbar(1 )=vbar(imtm1) ubar(imt)=ubar(2 ) vbar(imt)=vbar(2 ) #endif c c 2nd, calculate internal mode U at west face of U,V box c and V at north face of U,V box. c fx=p5 do 120 k=1,km do 120 i=2,imu fuw(i,k)=(ucl(i,k)*dzvqz(i,k,0)+ucl(i-1,k)*dzvqz(i-1,k,0))*fx fvn(i,k)=(vclp(i,k)*xzvqz(i,k)+vcl (i ,k)*dzvqz(i,k,0))*fx 120 continue c c 3rd, combine external and internal modes and add grid weight factor. c fx=dyu2r(j)*csr(j)*cst(j+1) do 130 k=1,km do 130 i=2,imu fuw(i,k)=(fuw(i,k)+ * p5*ubar(i)*(dzvqz(i-1,k,0)+dzvqz(i,k,0)))*csr(j) fvn(i,k)=(fvn(i,k)+p5*vbar(i)*(dzvqz(i,k,0)+xzvqz(i,k)))*fx 130 continue #ifdef cyclic c c Impose periodicity. c do 138 k=1,km fuw(imt,k)=fuw(2,k) fvn(imt,k)=fvn(2,k) 138 continue c #endif c----------------------------------------------------------------------- c Compute external mode velocities for row J+1. c----------------------------------------------------------------------- c c 1st, compute for TAU-1 time level c do 140 i=1,imtm1 diag1=pb(i+1,j+2)-pb(i ,j+1) diag2=pb(i ,j+2)-pb(i+1,j+1) ubarb(i)=-(diag1+diag2)*dyu2r(j+1)*min(c1,fkmup(i))*hv(i,j+1) vbarb(i)= (diag1-diag2)*dxu2r(i )*min(c1,fkmup(i))*hv(i,j+1)* * csr(j+1) 140 continue #ifdef cyclic ubarb(1 )=ubarb(imtm1) vbarb(1 )=vbarb(imtm1) ubarb(imt)=ubarb(2 ) vbarb(imt)=vbarb(2 ) #endif c c 2nd, compute for TAU time level. c do 150 i=1,imtm1 diag1=p (i+1,j+2)-p (i ,j+1) diag2=p (i ,j+2)-p (i+1,j+1) ubar (i)=-(diag1+diag2)*dyu2r(j+1)*min(c1,fkmup(i))*hv(i,j+1) vbar (i)= (diag1-diag2)*dxu2r(i )*min(c1,fkmup(i))*hv(i,j+1)* * csr(j+1) 150 continue #ifdef cyclic ubar(1 )=ubar(imtm1) vbar(1 )=vbar(imtm1) ubar(imt)=ubar(2 ) vbar(imt)=vbar(2 ) #endif c c----------------------------------------------------------------------- c Add external mode to internal mode for row J+1 (ocean points only). c----------------------------------------------------------------------- c do 160 k=1,km do 160 i=1,imu if(int(fkmup(i)).ge.kar(k)) then ubp(i,k)=ubp(i,k)+ubarb(i) vbp(i,k)=vbp(i,k)+vbarb(i) up (i,k)=up (i,k)+ubar (i) vp (i,k)=vp (i,k)+vbar (i) endif 160 continue c c----------------------------------------------------------------------- c Compute external mode velocities for row J (for bdy conditions). c----------------------------------------------------------------------- c c Save ubarp and ubarbp for bnd conditions. c do 163 i=1,imt ubarp(i)=ubar(i) vbarp(i)=vbar(i) ubarbp(i)=ubarb(i) vbarbp(i)=vbarb(i) 163 continue c c 1st, compute for TAU-1 time level c do 165 i=1,imtm1 diag1=pb(i+1,j+1)-pb(i ,j) diag2=pb(i ,j+1)-pb(i+1,j) ubarb(i)=-(diag1+diag2)*dyu2r(j)*min(c1,fkmu(i))*hv(i,j) vbarb(i)= (diag1-diag2)*dxu2r(i )*min(c1,fkmu(i))*hv(i,j)* * csr(j) 165 continue #ifdef cyclic ubarb(1 )=ubarb(imtm1) vbarb(1 )=vbarb(imtm1) ubarb(imt)=ubarb(2 ) vbarb(imt)=vbarb(2 ) #endif c c 2nd, compute for TAU time level. c do 168 i=1,imtm1 diag1=p (i+1,j+1)-p (i ,j) diag2=p (i ,j+1)-p (i+1,j) ubar (i)=-(diag1+diag2)*dyu2r(j)*min(c1,fkmu(i))*hv(i,j) vbar (i)= (diag1-diag2)*dxu2r(i )*min(c1,fkmup(i))*hv(i,j)* * csr(j) 168 continue #ifdef cyclic ubar(1 )=ubar(imtm1) vbar(1 )=vbar(imtm1) ubar(imt)=ubar(2 ) vbar(imt)=vbar(2 ) #endif c c----------------------------------------------------------------------- c Accumulate kinetic energy from row J+1 every NTSI timesteps or c in energy timesteps. c----------------------------------------------------------------------- c if(prntsi.or.(diagts.and.eots)) then do 170 k=1,km do 170 i=1,imu eke(i,k,1)=p5*(up(i,k)*up(i,k)+vp(i,k)*vp(i,k)) 170 continue fx=p5*cs(j+1)*dyu(j+1) do 180 k=1,km do 180 i=2,imum1 ektot(k,j)=ektot(k,j)+fx*eke(i,k,1)*c2*dzvqz(i,k,0)*dxuq(i,k) 180 continue endif c c----------------------------------------------------------------------- c Compute "omega" vertical velocity at the U,V points. c----------------------------------------------------------------------- c c Set "omega" vertical velocity at the surface to zero (rigid-lid). c do 190 i=2,imum1 wu(i,1)=c0 190 continue c c Compute change of W between levels. c do 200 k=1,km do 200 i=2,imum1 wu(i,k+1)= c2*( (fuw(i+1,k)-fuw (i ,k))*dxu2rq(i,k)+ * fvn(i ,k)-fvsu(i ,k) ) 200 continue c c Integrate downward from the surface. c do 210 k=1,km do 210 i=2,imum1 wu(i,k+1)=wu(i,k)+wu(i,k+1) 210 continue c c----------------------------------------------------------------------- c Compute standard vertical velocity from the "omega" vertical c velocity. c----------------------------------------------------------------------- c if((diagts.or.wrtts).and.eots) then do 220 k=1,km wvelu(1,k)=c0 wvelu(imu,k)=c0 do 220 i=2,imum1 wvelu(i,k)=p5*(wu(i,k+1)+wu(i,k))- * u(i,k)*dxu2r(i)* * ( (tdepth(i+1,k,jrn)- * tdepth(i ,k,jrn))*cstr(j+1)+ * (tdepth(i+1,k,jrs)- * tdepth(i ,k,jrs))*cstr(j ) )- * v(i,k)*dyu2r(j)* * ( (tdepth(i+1,k,jrn)- * tdepth(i+1,k,jrs))+ * (tdepth(i ,k,jrn)- * tdepth(i ,k,jrs)) ) 220 continue endif c c----------------------------------------------------------------------- c Compute hydrostatic pressure gradient. c----------------------------------------------------------------------- c #if defined usrdiagnostic & defined nesttime call dtime (dco) #endif call grad24_p(dpdx,dpdy,j) #if defined pressbias # if defined pressinbias if (itt.eq.1) then call press_bias0(dpdx,dpdy,j) endif # else call press_bias(dpdx,dpdy,j) # endif #endif #if defined usrdiagnostic & defined nesttime call dtime (dcpg) #endif #if defined ext_tide & defined tide_zero call tide_stress0(j) #endif c c----------------------------------------------------------------------- c Compute quantities for the the computation of vertical diffusion. c----------------------------------------------------------------------- c #ifndef barotropic do 230 k=1,kmm1 do 230 i=2,imum1 vmf(i,k,1)=vvc(i,k)*(ub(i,k)-ub(i,k+1))/dzzvqz(i,k+1,0) vmf(i,k,2)=vvc(i,k)*(vb(i,k)-vb(i,k+1))/dzzvqz(i,k+1,0) 230 continue #endif c c Set the K=0 elements of VMF to reflect wind stress and set the c K=KZ elements of VMF for bottom drag condition. c do 240 i=2,imum1 kz=kmu(i) vmf(i,0,1)=smf(i,1) vmf(i,0,2)=smf(i,2) vmf(i,kz,1)=bmf(i,1) vmf(i,kz,2)=bmf(i,2) 240 continue c c======================================================================= c End introductory section =========================================== c======================================================================= c c======================================================================= c Begin computation of the internal modes. ============ c The new values "UA" and "VA", will first be loaded with ============ c the time rate of change, and then updated. ============ c======================================================================= c #if !defined notadvm & !defined linear_physics c c----------------------------------------------------------------------- c Compute total advection of momentum. c----------------------------------------------------------------------- c c Compute zonal advection of momentum. c do 250 k=1,km do 250 i=2,imum1 UUx(i,k)=(fuw(i ,k)*(u(i-1,k)+u(i ,k))- * fuw(i+1,k)*(u(i ,k)+u(i+1,k)) * )*dxu2rq(i,k)/dzvqz(i,k,0) UVx(i,k)=(fuw(i ,k)*(v(i-1,k)+v(i ,k))- * fuw(i+1,k)*(v(i ,k)+v(i+1,k)) * )*dxu2rq(i,k)/dzvqz(i,k,0) 250 continue c c Compute meridional advection of momentum. c do 260 k=1,km do 260 i=2,imum1 VUy(i,k)= (fvsu(i,k)*(u (i,k)+um(i,k))- * fvn (i,k)*(up(i,k)+u (i,k)))/dzvqz(i,k,0) VVy(i,k)= (fvsu(i,k)*(v (i,k)+vm(i,k))- * fvn (i,k)*(vp(i,k)+v (i,k)))/dzvqz(i,k,0) 260 continue c c Compute vertical advection of momentum. c #ifndef barotropic do 270 k=2,km do 270 i=2,imum1 tempa(i,k)=wu(i,k)*(u(i,k-1)+u(i,k)) tempb(i,k)=wu(i,k)*(v(i,k-1)+v(i,k)) 270 continue #endif do 280 i=2,imum1 tempa(i, 1)=wu(i, 1)*u(i, 1) tempa(i,kmp1)=wu(i,kmp1)*u(i,km) tempb(i, 1)=wu(i, 1)*v(i, 1) tempb(i,kmp1)=wu(i,kmp1)*v(i,km) 280 continue c do 290 k=1,km do 290 i=2,imum1 WUz(i,k)=(tempa(i,k+1)-tempa(i,k))*dzv2rqz(i,k) WVz(i,k)=(tempb(i,k+1)-tempb(i,k))*dzv2rqz(i,k) 290 continue #endif c c----------------------------------------------------------------------- c Compute horizontal diffusion of momentum (evaluate at tau-1 tstep). c----------------------------------------------------------------------- c if((mixvel.eq.2).or.(mixvel.eq.3)) then if(mixvel.eq.2) then call lapv_depth(j,ubm,ub,ubp,Uxx,Uyy,Vmet) call lapv_depth(j,vbm,vb,vbp,Vxx,Vyy,Umet) else call lapv_lev(j,km,ubm,ub,ubp,Uxx,Uyy,Vmet) call lapv_lev(j,km,vbm,vb,vbp,Vxx,Vyy,Umet) endif do 300 k=1,km do 300 i=2,imum1 Uxx(i,k)=am*Uxx(i,k) Uyy(i,k)=am*Uyy(i,k) Vxx(i,k)=am*Vxx(i,k) Vyy(i,k)=am*Vyy(i,k) Umet(i,k)=am*(lpmtvl(j)*ub(i,k)-Umet(i,k)) Vmet(i,k)=am*(lpmtvl(j)*vb(i,k)+Vmet(i,k)) 300 continue else do 305 k=1,km do 305 i=2,imum1 Uxx(i,k)=c0 Uyy(i,k)=c0 Vxx(i,k)=c0 Vyy(i,k)=c0 Umet(i,k)=c0 Vmet(i,k)=c0 305 continue endif c c----------------------------------------------------------------------- c Compute vertical diffusion of momentum. c----------------------------------------------------------------------- c do 310 k=1,km do 310 i=2,imum1 Uzz(i,k)=(vmf(i,k-1,1)-vmf(i,k,1))/dzvqz(i,k,0) Vzz(i,k)=(vmf(i,k-1,2)-vmf(i,k,2))/dzvqz(i,k,0) 310 continue c c----------------------------------------------------------------------- c Compute Coriolis force (evaluate on TAU timestep for explicit c treatment, or evaluate on TAU-1 timestep for implicit treatment c with remainder of term to be added later). c----------------------------------------------------------------------- c if(acor.eq.c0) then do 320 k=1,km do 320 i=2,imum1 fx=c2*omega*sine(i,j) fV(i,k)= fx*v(i,k) fU(i,k)=-fx*u(i,k) 320 continue else do 330 k=1,km do 330 i=2,imum1 fx=c2*omega*sine(i,j) fV(i,k)= fx*vb(i,k) fU(i,k)=-fx*ub(i,k) 330 continue endif c #if defined bndy_rlx | defined botfrc | defined cstfrc c----------------------------------------------------------------------- c Compute boundary relaxation forcing. c----------------------------------------------------------------------- c do 335 k=1,km # if defined botfrc & !defined cstfrc fx = facbtfrc(k) # endif do 335 i=2,imum1 # if defined cstfrc & defined botfrc fx = faccsfrc(i,j)+facbtfrc(k) # elif defined cstfrc & !defined botfrc fx = faccsfrc(i,j) # endif # ifdef bndy_rlx if(itt.eq.1)then u_0(i,j,k)=ub(i,k) v_0(i,j,k)=vb(i,k) endif Ubrlx(i,k)=vfacbrlx(i,j)*(u_0(i,j,k)-ub(i,k)) Vbrlx(i,k)=vfacbrlx(i,j)*(v_0(i,j,k)-vb(i,k)) # if defined cstfrc | defined botfrc Ubrlx(i,k)=Ubrlx(i,k)-ub(i,k)*fx Vbrlx(i,k)=Vbrlx(i,k)-vb(i,k)*fx # endif # else Ubrlx(i,k)=-ub(i,k)*fx Vbrlx(i,k)=-vb(i,k)*fx # endif 335 continue c #endif c----------------------------------------------------------------------- c Construct the forcing terms for momentum and set land points to zero. c----------------------------------------------------------------------- c do 340 k=1,km do 340 i=2,imum1 ua(i,k)=gm(i,k)*( #if !defined notadvm & !defined linear_physics * UUx(i,k)+VUy(i,k)+WUz(i,k) # if defined ext_tide & defined tide_zero * +tUUx(i,k)+tVUy(i,k)+tWUz(i,k) # endif #endif * +Uxx(i,k)+Uyy(i,k)+Umet(i,k) #ifdef explicitvmix * +Uzz(i,k) #else * +Uzz(i,k)*(c1-aidif) #endif #if defined bndy_rlx | defined botfrc | defined cstfrc * +Ubrlx(i,k) #endif * +fV(i,k)-dpdx(i,k)) va(i,k)=gm(i,k)*( #if !defined notadvm & !defined linear_physics * UVx(i,k)+VVy(i,k)+WVz(i,k) # if defined ext_tide & defined tide_zero * +tUVx(i,k)+tVVy(i,k)+tWVz(i,k) # endif #endif * +Vxx(i,k)+Vyy(i,k)+Vmet(i,k) #ifdef explicitvmix * +Vzz(i,k) #else * +Vzz(i,k)*(c1-aidif) #endif #if defined bndy_rlx | defined botfrc | defined cstfrc * +Vbrlx(i,k) #endif * +fU(i,k)-dpdy(i,k)) 340 continue #ifndef explicitvmix c c----------------------------------------------------------------------- c Solve vertical diffusion implicitly. c----------------------------------------------------------------------- c c Update the velocity without implicit vertical diffusion. c do 350 k=1,km twodt(k)=c2dtuv do 350 i=2,imum1 ua(i,k)=gm(i,k)*ub(i,k)+c2dtuv*ua(i,k) va(i,k)=gm(i,k)*vb(i,k)+c2dtuv*va(i,k) 350 continue c c Store terms to compute implicit vertical diffusion on diagnostic c timesteps. c if(diagts.and.eots) then do 360 k=1,km do 360 i=2,imum1 tempa(i,k)=ua(i,k) tempb(i,k)=va(i,k) 360 continue endif c c Add in the implicit vertical diffusion. c #ifndef bottom call invtri(ua,smf(1,1),bmf(1,1),vvc,twodt,kmu,dzvqz(1,1,0), * dzvurq,dzvlrq,gm,2,imum1,aidif) call invtri(va,smf(1,2),bmf(1,2),vvc,twodt,kmu,dzvqz(1,1,0), * dzvurq,dzvlrq,gm,2,imum1,aidif) #else call invtri(ua,smf(1,1),cmf,vvc,twodt,kmu,dzvqz(1,1,0), * dzvurq,dzvlrq,gm,2,imum1,aidif,0) call invtri(va,smf(1,2),cmf,vvc,twodt,kmu,dzvqz(1,1,0), * dzvurq,dzvlrq,gm,2,imum1,aidif,0) #endif c c Compute residual implicit vertical diffusion. c r2dtuv=c1/c2dtuv if(diagts.and.eots) then do 370 k=1,km do 370 i=2,imum1 tempa(i,k)=r2dtuv*(ua(i,k)-tempa(i,k)) tempb(i,k)=r2dtuv*(va(i,k)-tempb(i,k)) 370 continue endif c c Convert back to time change of velocity. c do 380 k=1,km do 380 i=2,imum1 ua(i,k)=r2dtuv*(ua(i,k)-ub(i,k))*gm(i,k) va(i,k)=r2dtuv*(va(i,k)-vb(i,k))*gm(i,k) 380 continue #endif c c----------------------------------------------------------------------- c Form time change of vertically averaged forcing. c----------------------------------------------------------------------- c c 1st, integrate time change vertically c do 390 i=1,imu zun(i)=c0 zvn(i)=c0 390 continue fx=c2dtsf do 400 k=1,km do 400 i=2,imum1 zun(i)=zun(i)+ua(i,k)*fx*dzvqz(i,k,0) zvn(i)=zvn(i)+va(i,k)*fx*dzvqz(i,k,0) 400 continue c c 2nd, form average by dividing by depth c do 410 i=2,imum1 zun(i)=zun(i)*hv(i,j) zvn(i)=zvn(i)*hv(i,j) 410 continue c c Load values for reduced physics boundary conditions at the southern c boundary. c if((iopt(4).eq.6).and.(j.eq.2)) then do 415 i=2,imum1 zunz(i,2,south)=zun(i) zvnz(i,2,south)=zvn(i) 415 continue endif #ifndef cyclic c c Load values for reduced physics boundary conditions at the western c and eastern boundaries. c if(iopt(4).eq.6) then zunz(j,2,west)=zun(2) zunz(j,2,east)=zun(imum1) zvnz(j,2,west)=zvn(2) zvnz(j,2,west)=zvn(imum1) endif #else c c Set Cyclic boundary conditions. c zun(1 )=zun(imum1) zvn(1 )=zvn(imum1) zun(imu)=zun(2) zvn(imu)=zvn(2) #endif c c----------------------------------------------------------------------- c Do analysis of internal mode forcing on energy timestep. Also c form vertical average for use later in external mode analysis. c----------------------------------------------------------------------- c if(diagts.and.eots) then fx=cs(j)*dyu(j) c do 420 ll=1,11 do 420 i=1,imu zuneng(i,ll)=c0 zvneng(i,ll)=c0 420 continue c c Compute KE change due to pressure term globally. Also integrate c pressure term over the regional volume. c do 430 k=1,km do 430 i=2,imum1 ueng(i,k)=-gm(i,k)*dpdx(i,k) veng(i,k)=-gm(i,k)*dpdy(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,2)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,2)=engint(k,2)+eke(i,k,2)*boxvol zuneng(i,2)=zuneng(i,2)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,2)=zvneng(i,2)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,2,1)=termbm(k,2,1)+ueng(i,k)*boxvol termbm(k,2,2)=termbm(k,2,2)+veng(i,k)*boxvol 430 continue c c Compute KE change due to the zonal advection of momentum c globally. Also integrate zonal advection term over the c regional volume. c do 440 k=1,km do 440 i=2,imum1 ueng(i,k)=gm(i,k)*UUx(i,k) veng(i,k)=gm(i,k)*UVx(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,3)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,3)=engint(k,3)+eke(i,k,3)*boxvol zuneng(i,3)=zuneng(i,3)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,3)=zvneng(i,3)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,3,1)=termbm(k,3,1)+ueng(i,k)*boxvol termbm(k,3,2)=termbm(k,3,2)+veng(i,k)*boxvol 440 continue c c Compute KE change due to the meridional advection of momentum c globally. Also integrate meridional advection term over the c regional volume. c do 450 k=1,km do 450 i=2,imum1 ueng(i,k)=gm(i,k)*VUy(i,k) veng(i,k)=gm(i,k)*VVy(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,4)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,4)=engint(k,4)+eke(i,k,4)*boxvol zuneng(i,4)=zuneng(i,4)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,4)=zvneng(i,4)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,4,1)=termbm(k,4,1)+ueng(i,k)*boxvol termbm(k,4,2)=termbm(k,4,2)+veng(i,k)*boxvol 450 continue c c Compute KE change due to the vertical advection of momentum globally. c Also integrate vertical advection term over the regional volume. c do 460 k=1,km do 460 i=2,imum1 ueng(i,k)=gm(i,k)*WUz(i,k) veng(i,k)=gm(i,k)*WVz(i,k) eke(i,k,5)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) engint(k,5)=engint(k,5)+eke(i,k,5)*boxvol zuneng(i,5)=zuneng(i,5)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,5)=zvneng(i,5)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,5,1)=termbm(k,5,1)+ueng(i,k)*boxvol termbm(k,5,2)=termbm(k,5,2)+veng(i,k)*boxvol 460 continue c c Compute KE change due to zonal diffusion of momentum globally. Also c integrate zonal diffusion term over the regional volume. c if((mixvel.eq.2).or.(mixvel.eq.3)) then do 470 k=1,km do 470 i=2,imum1 ueng(i,k)=gm(i,k)*Uxx(i,k) veng(i,k)=gm(i,k)*Vxx(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,6)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,6)=engint(k,6)+eke(i,k,6)*boxvol zuneng(i,6)=zuneng(i,6)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,6)=zvneng(i,6)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,6,1)=termbm(k,6,1)+ueng(i,k)*boxvol termbm(k,6,2)=termbm(k,6,2)+veng(i,k)*boxvol 470 continue endif c c Compute KE change due to meridional diffusion of momentum globally. c Also integrate meridional diffusion term over the regional volume. c if((mixvel.eq.2).or.(mixvel.eq.3)) then do 480 k=1,km do 480 i=2,imum1 ueng(i,k)=gm(i,k)*Uyy(i,k) veng(i,k)=gm(i,k)*Vyy(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,7)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,7)=engint(k,7)+eke(i,k,7)*boxvol zuneng(i,7)=zuneng(i,7)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,7)=zvneng(i,7)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,7,1)=termbm(k,7,1)+ueng(i,k)*boxvol termbm(k,7,2)=termbm(k,7,2)+veng(i,k)*boxvol 480 continue endif c c Compute KE change due to vertical diffusion of momentum. c #ifndef barotropic do 490 k=2,kmm1 do 490 i=2,imum1 # ifdef explicitvmix ueng(i,k)=gm(i,k)*Uzz(i,k) veng(i,k)=gm(i,k)*Vzz(i,k) # else ueng(i,k)=gm(i,k)*(Uzz(i,k)*(c1-aidif)+tempa(i,k)) veng(i,k)=gm(i,k)*(Vzz(i,k)*(c1-aidif)+tempb(i,k)) # endif boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,8)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,8)=engint(k,8)+eke(i,k,8)*boxvol zuneng(i,8)=zuneng(i,8)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,8)=zvneng(i,8)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,8,1)=termbm(k,8,1)+ueng(i,k)*boxvol termbm(k,8,2)=termbm(k,8,2)+veng(i,k)*boxvol 490 continue #endif c c Compute KE change due to wind stress. c do 500 i=2,imum1 ueng(i,1)=gm(i,1)*smf(i,1) veng(i,1)=gm(i,1)*smf(i,2) boxvol=fx*dxuq(i,1)*dzvqz(i,1,0) eke(i,1,8)=(ucl(i,1)*ueng(i,1)+vcl(i,1)*veng(i,1)) engint(1,9)=engint(1,9)+eke(i,1,8)*boxvol zuneng(i,9)=zuneng(i,9)+ueng(i,1)*dzvqz(i,1,0)*hv(i,j) zvneng(i,9)=zvneng(i,9)+veng(i,1)*dzvqz(i,1,0)*hv(i,j) smflx(1)=smflx(1)+ueng(i,1)*fx*dxuq(i,1) smflx(2)=smflx(2)+ueng(i,1)*fx*dxuq(i,1) 500 continue c c Compute KE change due to bottom drag. c do 510 i=2,imum1 kz=km ueng(i,kz)=-gm(i,kz)*bmf(i,1) veng(i,kz)=-gm(i,kz)*bmf(i,2) boxvol=fx*dxuq(i,kz)*dzvqz(i,kz,0) eke(i,kz,8)=(ucl(i,kz)*ueng(i,kz)+vcl(i,kz)*veng(i,kz)) engint(kz,10)=engint(kz,10)+eke(i,kz,8)*boxvol zuneng(i,10)=zuneng(i,10)+ueng(i,kz)*dzvqz(i,kz,0)*hv(i,j) zvneng(i,10)=zvneng(i,10)+veng(i,kz)*dzvqz(i,kz,0)*hv(i,j) 510 continue c c Compute KE change due to metric diffusion of momentum globally. c Also integrate metric diffusion term over the regional volume. c if((mixvel.eq.2).or.(mixvel.eq.3)) then do 515 k=1,km do 515 i=2,imum1 ueng(i,k)=gm(i,k)*Umet(i,k) veng(i,k)=gm(i,k)*Vmet(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,9)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,11)=engint(k,11)+eke(i,k,9)*boxvol zuneng(i,11)=zuneng(i,11)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,11)=zvneng(i,11)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,11,1)=termbm(k,11,1)+ueng(i,k)*boxvol termbm(k,11,2)=termbm(k,11,2)+veng(i,k)*boxvol 515 continue endif c c Compute KE change due to Coriolis term. c do 520 k=1,km do 520 i=2,imum1 ueng(i,k)= gm(i,k)*fV(i,k) veng(i,k)=-gm(i,k)*fU(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,10)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,12)=engint(k,12)+eke(i,k,10)*boxvol zuneng(i,12)=zuneng(i,12)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,12)=zvneng(i,12)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,12,1)=termbm(k,12,1)+ueng(i,k)*boxvol termbm(k,12,2)=termbm(k,12,2)+veng(i,k)*boxvol 520 continue sofar=MTERMS00 # if defined bndy_rlx | defined botfrc | defined cstfrc c c Compute KE change due to boundary relaxation term. c sofar=sofar+1 do 523 k=1,km do 523 i=2,imum1 ueng(i,k)=gm(i,k)*Ubrlx(i,k) veng(i,k)=gm(i,k)*Vbrlx(i,k) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) eke(i,k,sofar-2)=(ucl(i,k)*ueng(i,k)+vcl(i,k)*veng(i,k)) engint(k,sofar)=engint(k,sofar)+eke(i,k,sofar-2)*boxvol zuneng(i,sofar)=zuneng(i,sofar)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,sofar)=zvneng(i,sofar)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,sofar,1)=termbm(k,sofar,1)+ueng(i,k)*boxvol termbm(k,sofar,2)=termbm(k,sofar,2)+veng(i,k)*boxvol 523 continue # endif c c Accumulate average velocity within volume. c do 530 k=1,km do 530 i=2,imum1 boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) avwv(1)=avwv(1)+u(i,k)*boxvol avwv(2)=avwv(2)+v(i,k)*boxvol avwv(3)=avwv(3)+p5*(wu(i,k)+wu(i,k))*gm(i,k)*boxvol 530 continue endif c c----------------------------------------------------------------------- c Compute new velocities (with incorrect vertical means); c also, add in remainder of coriolis term if treated implicitly. c----------------------------------------------------------------------- c if(acor.eq.c0) then do 540 k=1,km do 540 i=2,imum1 #if !defined bndy_rlx | !defined imp_bnd_rlx ua(i,k)=ub(i,k)+c2dtuv*ua(i,k) va(i,k)=vb(i,k)+c2dtuv*va(i,k) # else fx =c2dtuv/(c1+p5*c2dtuv*vfacbrlx(i,j)) ua(i,k)=ub(i,k)+fx*ua(i,k) va(i,k)=vb(i,k)+fx*va(i,k) #endif 540 continue else do 550 k=1,km do 550 i=2,imum1 #if !defined bndy_rlx | !defined imp_bnd_rlx fx=c2dtuv*acor*c2*omega*sine(i,j) detmr=c1/(c1+fx*fx) tempa(i,k)=(ua(i,k)+fx*va(i,k))*detmr tempb(i,k)=(va(i,k)-fx*ua(i,k))*detmr ua(i,k)=ub(i,k)+c2dtuv*tempa(i,k) va(i,k)=vb(i,k)+c2dtuv*tempb(i,k) # else diag1=c1+p5*c2dtuv*vfacbrlx(i,j) diag2=c2dtuv*acor*c2*omega*sine(i,j) detmr=c1/(diag1*diag1+diag2*diag2) tempa(i,k)=(diag1*ua(i,k)+diag2*va(i,k))*detmr tempb(i,k)=(diag1*va(i,k)-diag2*ua(i,k))*detmr ua(i,k)=ub(i,k)+c2dtuv*tempa(i,k) va(i,k)=vb(i,k)+c2dtuv*tempb(i,k) #endif 550 continue endif c c----------------------------------------------------------------------- c Determine the incorrect vertical means of the new velocities. c----------------------------------------------------------------------- c do 560 i=2,imum1 umean(i)=c0 vmean(i)=c0 560 continue do 570 k=1,km do 570 i=2,imum1 umean(i)=umean(i)+ua(i,k)*dzvqz(i,k,0) vmean(i)=vmean(i)+va(i,k)*dzvqz(i,k,0) 570 continue do 580 i=2,imum1 umean(i)=umean(i)*hvav(i,j) vmean(i)=vmean(i)*hvav(i,j) 580 continue c c----------------------------------------------------------------------- c Subtract incorrect vertical mean to get internal mode. c----------------------------------------------------------------------- c do 590 k=1,km do 590 i=2,imum1 ua(i,k)=gm(i,k)*(ua(i,k)-umean(i)) va(i,k)=gm(i,k)*(va(i,k)-vmean(i)) 590 continue #ifdef cyclic c c Set Cyclic boundary conditions. c do 595 k=1,km ua(1 ,k)=ua(imum1,k) va(1 ,k)=va(imum1,k) ua(imu,k)=ua(2 ,k) va(imu,k)=va(2 ,k) 595 continue #endif c c----------------------------------------------------------------------- c Compute total change of KE of internal mode on energy timestep. c----------------------------------------------------------------------- c if(diagts.and.eots) then fx=cs(j)*dyu(j)/c2dtuv do 600 k=1,km do 600 i=2,imum1 ueng(i,k)=gm(i,k)*(ua(i,k)-ub(i,k)) veng(i,k)=gm(i,k)*(va(i,k)-vb(i,k)) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) engint(k,1)=engint(k,1)+(ucl(i,k)*ueng(i,k)+ * vcl(i,k)*veng(i,k))*boxvol zuneng(i,1)=zuneng(i,1)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,1)=zvneng(i,1)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,1,1)=termbm(k,1,1)+ueng(i,k)*boxvol termbm(k,1,2)=termbm(k,1,2)+veng(i,k)*boxvol 600 continue # if defined bndy_rlx & defined imp_bnd_rlx c c Add implicit part of boundary relaxation term into momentum term balance. c do 605 k=1,km do 605 i=2,imum1 fx=p5*c2dtuv*vfacbrlx(i,j) ueng(i,k)=gm(i,k)*(ua(i,k)-ub(i,k)) veng(i,k)=gm(i,k)*(va(i,k)-vb(i,k)) boxvol=fx*dxuq(i,k)*dzvqz(i,k,0) zuneng(i,13)=zuneng(i,13)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j) zvneng(i,13)=zvneng(i,13)+veng(i,k)*dzvqz(i,k,0)*hv(i,j) termbm(k,13,1)=termbm(k,13,1)+ueng(i,k)*boxvol termbm(k,13,2)=termbm(k,13,2)+veng(i,k)*boxvol 605 continue # endif c c Add implicit part of the Coriolis term into momentum term balance. c if(acor.ne.c0) then do 610 k=1,km do 610 i=2,imum1 fx=acor*c2*omega*sine(i,j) ueng(i,k)=gm(i,k)*(ua(i,k)-ub(i,k)) veng(i,k)=gm(i,k)*(va(i,k)-vb(i,k)) boxfac=fx*dxuq(i,k)*cs(j)*dyu(j)*dzvqz(i,k,0) zuneng(i,12)=zuneng(i,12)+ueng(i,k)*dzvqz(i,k,0)*hv(i,j)*fx zvneng(i,12)=zvneng(i,12)+veng(i,k)*dzvqz(i,k,0)*hv(i,j)*fx termbm(k,12,1)=termbm(k,12,1)+ueng(i,k)*boxfac termbm(k,12,2)=termbm(k,12,2)+veng(i,k)*boxfac 610 continue endif #ifdef cyclic do 615 ll=2,8 zuneng(1,ll)=zuneng(imum1,ll) zvneng(1,ll)=zvneng(imum1,ll) 615 continue #endif endif c c======================================================================= c End computation of internal modes ================================== c======================================================================= c c======================================================================= c Begin computation of vorticity for input to "RELAX" ================ c======================================================================= c c----------------------------------------------------------------------- c Form curl of time change in vertically averaged equations. c----------------------------------------------------------------------- #ifndef islands c c Non-active vorticity points are kept at zero for convenience in c RELAX. To accomplish this, ZTD is computed only between given c starting and ending indices. c do 620 l=1,lseg istr=isz(j,l) if(istr.eq.0) go to 630 iend=iez(j,l) #else c c All vorticity points are computed to ensure the line integral of c hole relaxation is well defined. c istr=2 iend=imtm1 #endif do 620 i=istr,iend ztd(i,j)=((zun(i)*dxu(i)+zun(i-1)*dxu(i-1))*cs(j ) * -(zus(i)*dxu(i)+zus(i-1)*dxu(i-1))*cs(j-1)) ztd(i,j)=(((zvn(i)-zvn(i-1))*dyu(j ) * +(zvs(i)-zvs(i-1))*dyu(j-1) * -ztd(i,j))*dxt2r(i)*dytr(j))*cstr(j) 620 continue 630 continue #ifdef cyclic c c -- required in Shapiro filter [at least] c ztd(1,j)=ztd(imum1,j) ztd(imu,j)=ztd(2,j) #endif c c----------------------------------------------------------------------- c Do analysis of external mode forcing on energy timestep. c----------------------------------------------------------------------- c if(diagts.and.eots) then do 660 ll=1,mterms do 640 i=2,imtm1 tempa(i,1)=-p(i,j)*( * ((zvneng(i,ll)-zvneng(i-1,ll))*dyu(j) * +(zvseng(i,ll)-zvseng(i-1,ll))*dyu(j-1)) * *dxt2r(i)*cstr(j)*dytr(j) * -((zuneng(i,ll)*dxu(i)+zuneng(i-1,ll)*dxu(i-1)) * *cs(j) * -(zuseng(i,ll)*dxu(i)+zuseng(i-1,ll)*dxu(i-1)) * *cs(j-1))*dyt2r(j)*dxtr(i)*cstr(j)) * *dxt(i)*cst(j)*dyt(j) 640 continue do 650 i=2,imtm1 engext(ll)=engext(ll)+tempa(i,1) 650 continue 660 continue endif c c======================================================================= c End computation of vorticity ======================================= c======================================================================= c c----------------------------------------------------------------------- c Transfer quantities computed to the north of the present row c to be defined to the south in the computation of the next row. c----------------------------------------------------------------------- c fx=cs(j)*dyu(j)*csr(j+1)*dyur(j+1) do 670 k=1,km do 670 i=1,imt fvsu(i,k)=fvn(i,k)*fx 670 continue do 680 i=1,imt zus(i)=zun(i) zvs(i)=zvn(i) 680 continue if(diagts.and.eots) then do 690 ll=2,mterms do 690 i=1,imt zuseng(i,ll)=zuneng(i,ll) zvseng(i,ll)=zvneng(i,ll) 690 continue endif c c Set boundary condition arrays for reduced physics option. c if(iopt(4).eq.6) then #ifndef cyclic zusz(j,2,west)=zus(2) zusz(j,2,east)=zus(imum1) zvsz(j,2,west)=zvs(2) zvsz(j,2,east)=zvs(imum1) #endif if(j.eq.jmtm2) then do 700 i=2,imum1 zusz(i,2,north)=zus(i) zvsz(i,2,north)=zvs(i) 700 continue endif endif #endif #if defined usrdiagnostic & defined nesttime call dtime (ddco) dco(1) = dco(1) + ddco(1) dco(2) = dco(2) + ddco(2) #endif c return end