forked from darwinproject/darwin3
-
Notifications
You must be signed in to change notification settings - Fork 1
/
ini_cartesian_grid.F
155 lines (136 loc) · 5.07 KB
/
ini_cartesian_grid.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
#include "CPP_OPTIONS.h"
CBOP
C !ROUTINE: INI_CARTESIAN_GRID
C !INTERFACE:
SUBROUTINE INI_CARTESIAN_GRID( myThid )
C !DESCRIPTION: \bv
C *==========================================================*
C | SUBROUTINE INI_CARTESIAN_GRID
C | o Initialise model coordinate system
C *==========================================================*
C | The grid arrays, initialised here, are used throughout
C | the code in evaluating gradients, integrals and spatial
C | avarages. This routine is called separately by each
C | thread and initialises only the region of the domain
C | it is "responsible" for.
C | Under the cartesian grid mode primitive distances
C | in X and Y are in metres. Distance in Z are in m or Pa
C | depending on the vertical gridding mode.
C *==========================================================*
C \ev
C !USES:
IMPLICIT NONE
C === Global variables ===
#include "SIZE.h"
#include "EEPARAMS.h"
#include "PARAMS.h"
#include "GRID.h"
C !INPUT/OUTPUT PARAMETERS:
C == Routine arguments ==
C myThid :: my Thread Id Number
INTEGER myThid
C !LOCAL VARIABLES:
C == Local variables ==
C bi,bj :: tile indices
C i, j :: loop counters
C delXloc :: mesh spacing in X direction
C delYloc :: mesh spacing in Y direction
C xGloc :: mesh corner-point location (local "Long" real array type)
C yGloc :: mesh corner-point location (local "Long" real array type)
INTEGER bi, bj
INTEGER i, j
INTEGER gridNx, gridNy
C NOTICE the extended range of indices!!
_RL delXloc(0-OLx:sNx+OLx)
_RL delYloc(0-OLy:sNy+OLy)
C NOTICE the extended range of indices!!
_RL xGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
_RL yGloc(1-OLx:sNx+OLx+1,1-OLy:sNy+OLy+1)
CEOP
C-- For each tile ...
DO bj = myByLo(myThid), myByHi(myThid)
DO bi = myBxLo(myThid), myBxHi(myThid)
C-- set tile local mesh (same units as delX,deY)
C corresponding to coordinates of cell corners for N+1 grid-lines
CALL INI_LOCAL_GRID(
O xGloc, yGloc,
O delXloc, delYloc,
O gridNx, gridNy,
I bi, bj, myThid )
C-- Make a permanent copy of [xGloc,yGloc] in [xG,yG]
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
xG(i,j,bi,bj) = xGloc(i,j)
yG(i,j,bi,bj) = yGloc(i,j)
ENDDO
ENDDO
C-- Calculate [xC,yC], coordinates of cell centers
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
C by averaging
xC(i,j,bi,bj) = 0.25 _d 0*(
& xGloc(i,j)+xGloc(i+1,j)+xGloc(i,j+1)+xGloc(i+1,j+1) )
yC(i,j,bi,bj) = 0.25 _d 0*(
& yGloc(i,j)+yGloc(i+1,j)+yGloc(i,j+1)+yGloc(i+1,j+1) )
ENDDO
ENDDO
C-- Calculate [dxF,dyF], lengths between cell faces (through center)
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dxF(i,j,bi,bj) = delXloc(i)
dyF(i,j,bi,bj) = delYloc(j)
ENDDO
ENDDO
C-- Calculate [dxG,dyG], lengths along cell boundaries
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
dxG(i,j,bi,bj) = delXloc(i)
dyG(i,j,bi,bj) = delYloc(j)
ENDDO
ENDDO
C-- The following arrays are not defined in some parts of the halo
C region. We set them to zero here for safety.
C Note: this is now done earlier in main S/R INI_GRID
C-- Calculate [dxC], zonal length between cell centers
DO j=1-OLy,sNy+OLy
DO i=1-OLx+1,sNx+OLx ! NOTE range
dxC(i,j,bi,bj) = 0.5 _d 0*(dxF(i,j,bi,bj)+dxF(i-1,j,bi,bj))
ENDDO
ENDDO
C-- Calculate [dyC], meridional length between cell centers
DO j=1-OLy+1,sNy+OLy ! NOTE range
DO i=1-OLx,sNx+OLx
dyC(i,j,bi,bj) = 0.5 _d 0*(dyF(i,j,bi,bj)+dyF(i,j-1,bi,bj))
ENDDO
ENDDO
C-- Calculate [dxV,dyU], length between velocity points (through corners)
DO j=1-OLy+1,sNy+OLy ! NOTE range
DO i=1-OLx+1,sNx+OLx ! NOTE range
C by averaging (method I)
dxV(i,j,bi,bj) = 0.5 _d 0*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
dyU(i,j,bi,bj) = 0.5 _d 0*(dyG(i,j,bi,bj)+dyG(i,j-1,bi,bj))
C by averaging (method II)
c dxV(i,j,bi,bj) = 0.5*(dxG(i,j,bi,bj)+dxG(i-1,j,bi,bj))
c dyU(i,j,bi,bj) = 0.5*(dyC(i,j,bi,bj)+dyC(i-1,j,bi,bj))
ENDDO
ENDDO
C-- Calculate vertical face area
DO j=1-OLy,sNy+OLy
DO i=1-OLx,sNx+OLx
rA (i,j,bi,bj) = dxF(i,j,bi,bj)*dyF(i,j,bi,bj)
rAw(i,j,bi,bj) = dxC(i,j,bi,bj)*dyG(i,j,bi,bj)
rAs(i,j,bi,bj) = dxG(i,j,bi,bj)*dyC(i,j,bi,bj)
rAz(i,j,bi,bj) = dxV(i,j,bi,bj)*dyU(i,j,bi,bj)
C-- Set trigonometric terms & grid orientation:
C Note: this is now done earlier in main S/R INI_GRID
c tanPhiAtU(i,j,bi,bj) = 0.
c tanPhiAtV(i,j,bi,bj) = 0.
c angleCosC(i,j,bi,bj) = 1.
c angleSinC(i,j,bi,bj) = 0.
ENDDO
ENDDO
C-- end bi,bj loops
ENDDO
ENDDO
RETURN
END