Ocean Modeling Discussion

First, I run the standard riverplume2 case without any change.
Then I add some print command in step2d_LF_AM3.h and step3d_uv.F to print out the value of point sources input (like: Qbar(is), cff) for checking.

################################################################################################
In the step2d_LF_AM3.h the print statement are add like:

!
!-----------------------------------------------------------------------
! Apply mass point sources.
!-----------------------------------------------------------------------
!
DO j=Jstr-1,Jend+1
DO i=Istr-1,Iend+1
Dnew(i,j)=zeta(i,j,knew)+h(i,j)
END DO
END DO
DO is=1,Nsrc
i=Isrc(is)
j=Jsrc(is)
IF (((IstrR.le.i).and.(i.le.IendR)).and. &
& ((JstrR.le.j).and.(j.le.JendR))) THEN
IF (INT(Dsrc(is)).eq.0) THEN
cff=1.0_r8/(on_u(i,j)*0.5_r8*(Dnew(i-1,j)+Dnew(i,j)))
ubar(i,j,knew)=Qbar(is)*cff
if (is.eq.1)then !! rex
print *,"rex ubar(i,j,knew)", is," ", i," ", j," ",knew ," ",Qbar(is), " ", cff !! rex
end if !! rex

# ifdef SOLVE3D
DU_avg1(i,j)=Qbar(is)
if (is.eq.1)then !! rex
print *,"rex DU_avg1(i,j)", is," ", i," ", j," ",Qbar(is), " ", cff !! rex
end if !! rex

# endif
ELSE
cff=1.0_r8/(om_v(i,j)*0.5_r8*(Dnew(i,j-1)+Dnew(i,j)))
vbar(i,j,knew)=Qbar(is)*cff
if (is.eq.1)then !! rex
print *,"rex vbar(i,j,knew)", is," ", i," ", j," ",knew ," ",Qbar(is), " ", cff !! rex
end if !! rex

# ifdef SOLVE3D
DV_avg1(i,j)=Qbar(is)
if (is.eq.1)then !! rex
print *,"rex DV_avg1(i,j)", is," ", i," ", j," ",Qbar(is), " ", cff !! rex
end if !! rex

# endif
END IF
END IF
END DO
# endif
############################################################
############################################################
And in the step3d_uv.F the print command add in:
!
!-----------------------------------------------------------------------
! Apply mass point sources.
!-----------------------------------------------------------------------
!
DO is=1,Nsrc
i=Isrc(is)
j=Jsrc(is)
IF (((IstrR.le.i).and.(i.le.IendR)).and. &
& ((JstrR.le.j).and.(j.le.JendR))) THEN
IF (INT(Dsrc(is)).eq.0) THEN
DO k=1,N(ng)
cff1=1.0_r8/(on_u(i,j)* &
& 0.5_r8*(z_w(i-1,j,k)-z_w(i-1,j,k-1)+ &
& z_w(i ,j,k)-z_w(i ,j,k-1)))
u(i,j,k,nnew)=Qsrc(is,k)*cff1
if (is.eq.1)then !! rex
print *,"rex u(i,j,k,nnew)", is," ", i," ", j," ", k," ",nnew ," ",Qsrc(is,k), " ", cff1 !! rex
end if !! rex

END DO
ELSE
DO k=1,N(ng)
cff1=1.0_r8/(om_v(i,j)* &
& 0.5_r8*(z_w(i,j-1,k)-z_w(i,j-1,k-1)+ &
& z_w(i,j ,k)-z_w(i,j ,k-1)))
v(i,j,k,nnew)=Qsrc(is,k)*cff1
if (is.eq.1)then !! rex
print *,"rex v(i,j,k,nnew)", is," ", i," ", j," ", k," ",nnew ," ",Qsrc(is,k), " ", cff1 !! rex
end if !! rex

END DO
END IF
END IF
END DO
# endif

#####################################################################################################


And after adding print, the case result of u,v, S are changed as the figures attached.
Furthermore I test to add a print command in the upwelling case, it wouldn't change the result.
I have some questions:
1, why adding the print commands will change result in the riverplume2 case?
2, how to solve this problem?