fortranを用いたせん断速度のプログラムの書き方

fortran90を使って2次元成分でx方向の速度をu(i,j)、y方向の速度をv(i,j)としたときの
(du/dy+dv/dx)というせん断速度を求めるプログラムの書き方を教えてください。

-### 試したこと
離散化はスタガード格子を用いて中心差分近似で行い、
以下のように定義されたu(i,j)点での不当間隔格子による差分をもちいて
[du/dy+dv/dx]=((u(i,j+1)+u(i,j))/dyh(j))-((u(i,j)-u(i,j-1)/dyh(j-1))
+((v(i+1,j)+v(i,j+1)-v(i,j)-v(i+1,j-1))/dxh(i))
というような式のプログラムを書いてみたのですが、うまくいきませんでした。

fortran
1    dudy_u   = ((u(i,j+1)-u(i,j))/dyh(j))-((u(i,j)-u(i,j-1))/dyh(j-1))
2    dvdx_u   = (v(i+1,j)+v(i,j-1)-v(i,j)-v(i+1,j-1))/dxh(i)
3    D(i,j) = dudy_u+dvdx_u

coco_bauer

2019/07/04 06:48

最初の２行は代入式(左辺が変数)なので実行可能ですが、３行目のコードは左辺が式になっているのでシンタックスエラー（文法エラー）になります（実行できません）「うまくいきませんでした」というのは、そういうエラーが出たという事ですか？

t.ewjspo

2019/07/04 06:52

すみませんそこは書き間違えてました。上手くいかなかったのは、結果が負の値になってしまいました。

coco_bauer

2019/07/04 07:04

境界条件や外部条件によって、注目している点の移動方向は変わりますから、それが負の方向（マイナスの速度）となる可能性はあります。負の値になったから、間違いとは言い切れません。（まぁ、どのような問題を解こうとしているのか、具体的な情報が何もないので、推測にすぎませんが）こういう問題は、質問に書かれているような図を書いて、そこにカラーマーカーで力の方向、点の移動方向などを書きこむとイメージがつかめます。

t.ewjspo

2019/07/04 07:26

そうなんですね。この値から粘度をpower lawモデルという式に代入して導出しているのですが、せん断速度が負の値の場合、粘度が負になってしまうのですが、この場合どう修正すればよいですか？

t.ewjspo

2019/07/04 07:29

粘度Dの式は D=D0*[dv/dx+du/dy]^n-1 という式です。