-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathangles.F90
209 lines (180 loc) · 7.78 KB
/
angles.F90
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
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
!> @file
!! @brief Determine the rotation angle on center and
!! corner points
!! @author [email protected]
!!
!> This module finds the rotation angle for at both center and corner points
!! It utilizes the MOM6 function modulo_around_point
!! @author [email protected]
module angles
use gengrid_kinds, only : dbl_kind, int_kind
use grdvars, only : debug
implicit none
contains
!> Find the rotation angle on center (Bu) grid points
!!
!! @param[in] ni the i-dimension of the grid
!! @param[in] nj the j-dimension of the grid
!! @param[in] xangCt the angle across the tripole seam
!! @param[in] anglet the rotation angle on Ct points
!! @param[out] angle the rotation angle on Bu points
!! @author [email protected]
subroutine find_angq(iind,jind,xangCt,anglet,angle)
integer , intent(in) :: iind(2),jind(2)
real(dbl_kind), intent(in) :: xangCt(:)
real(dbl_kind), intent(in) :: anglet(:,:)
real(dbl_kind), intent(out) :: angle(:,:)
! local variables
integer :: i,j
real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
real(dbl_kind) :: p25 = 0.25
!---------------------------------------------------------------------
! find the angle on corners using the same relationship CICE uses
! internally to calculate angles on Ct using angles on Bu
!
! w-----------------0 Ct(i+1,j+1)
! | |
! ----------Bu(i,j)---------- Bu lies on seam at j=nj
! | |
! Ct(i,j) sw----------------s
!
!---------------------------------------------------------------------
angle = 0.0
do j = jind(1)+1,jind(2)
do i = iind(1),iind(2)-1
if (j .lt. jind(2)) then
angle_0 = anglet(i+1,j+1)
angle_w = anglet(i, j+1)
angle_s = anglet(i+1,j )
angle_sw = anglet(i ,j )
else
angle_0 = xangCt(i+1 )
angle_w = xangCt(i )
angle_s = anglet(i+1,j)
angle_sw = anglet(i, j)
end if
angle(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))
if (abs(angle(i,j)) .le. 1.0e-10)angle(i,j) = 0.0
enddo
enddo
!angle(iind(2),:) = -angle(1,:)
end subroutine find_angq
!> Verify the rotation angle on center (Ct) grid points using angle on corner
!! (Bu) grid points
!!
!! @param[in] ni the i-dimension of the grid
!! @param[in] nj the j-dimension of the grid
!! @param[in] angle the rotation angle on Bu points
!! @param[out] angchk the rotation angle on Ct points
!! @author [email protected]
subroutine find_angchk(iind,jind,angle,angchk)
integer , intent(in) :: iind(2),jind(2)
real(dbl_kind), intent(in) :: angle(:,:)
real(dbl_kind), intent(out) :: angchk(:,:)
! local variables
integer :: i,j
real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
real(dbl_kind) :: p25 = 0.25
!---------------------------------------------------------------------
! check: calculate anglet from angle on corners as CICE does internally.
! since angle changes sign between CICE and MOM6, (-1)*angchk ~ anglet
!
! w-----------------0 Bu(i,j)
! | |
! | Ct(i,j) |
! | |
! Bu(i-1,j-1) sw----------------s
!
!---------------------------------------------------------------------
angchk = 0.0
do j = jind(1)+1,jind(2)
do i = iind(1)+1,iind(2)
angle_0 = angle(i ,j )
angle_w = angle(i-1,j )
angle_s = angle(i, j-1)
angle_sw = angle(i-1,j-1)
angchk(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))
enddo
enddo
!angchk(1,:) = -angchk(iind(2),:)
end subroutine find_angchk
!> Find the rotation angle on center (Ct) grid points
!!
!! @param[in] ni the i-dimension of the grid
!! @param[in] nj the j-dimension of the grid
!! @param[in] lonBu the longitudes of the corner grid points
!! @param[in] latBu the latitudes of the corner grid points
!! @param[in] lonCt the longitudes of the center grid points
!! @param[out] anglet the rotation angle on Ct points
!! @author [email protected]
subroutine find_ang(iind,jind,lonBu,latBu,lonCt,anglet)
integer , intent(in) :: iind(2),jind(2)
real(dbl_kind), intent(in) :: lonBu(:,:)
real(dbl_kind), intent(in) :: latBu(:,:)
real(dbl_kind), intent(in) :: lonCt(:,:)
real(dbl_kind), intent(out) :: anglet(:,:)
! local variables
integer :: i,j,m,n
integer :: ii,jj
! from geolonB fix in MOM6
real(dbl_kind) :: len_lon ! The periodic range of longitudes, usually 360 degrees.
real(dbl_kind) :: pi_720deg ! One quarter the conversion factor from degrees to radians.
real(dbl_kind) :: lonB(2,2) ! The longitude of a point, shifted to have about the same value.
real(dbl_kind) :: lon_scale = 0.0
!---------------------------------------------------------------------
! rotation angle for "use_bugs" = false case from MOM6
! src/initialization/MOM_shared_initialization.F90 but allow for not
! having halo values
! note this does not reproduce sin_rot,cos_rot found in MOM6 output
! differences are ~O 10-6
!---------------------------------------------------------------------
anglet = 0.0
pi_720deg = atan(1.0) / 180.0
len_lon = 360.0
do j=jind(1),jind(2); do i = iind(1),iind(2)
do n=1,2 ; do m=1,2
jj = J+n-2; ii = I+m-2
if(jj .eq. 0)jj = 1
if(ii .eq. 0)ii = iind(2)
lonB(m,n) = modulo_around_point(LonBu(ii,jj), LonCt(i,j), len_lon)
! lonB(m,n) = modulo_around_point(LonBu(I+m-2,J+n-2), LonCt(i,j), len_lon)
enddo; enddo
jj = j-1; ii = i-1
if(jj .eq. 0)jj = 1
if(ii .eq. 0)ii = iind(2)
lon_scale = cos(pi_720deg*((LatBu(ii,jj) + LatBu(I,J)) + &
(LatBu(I,jj) + LatBu(ii,J)) ) )
anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
(LatBu(ii,J) - LatBu(I,jj)) + &
(LatBu(I,J) - LatBu(ii,jj)) )
!lon_scale = cos(pi_720deg*((LatBu(I-1,J-1) + LatBu(I,J)) + &
! (LatBu(I,J-1) + LatBu(I-1,J)) ) )
!anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
! (LatBu(I-1,J) - LatBu(I,J-1)) + &
! (LatBu(I,J) - LatBu(I-1,J-1)) )
enddo; enddo
end subroutine find_ang
! -----------------------------------------------------------------------------
!> Return the modulo value of x in an interval [xc-(Lx/2) xc+(Lx/2)]
!! If Lx<=0, then it returns x without applying modulo arithmetic.
!!
!! From src/initialization/MOM_shared_initialization.F90:
!! @param[in] x Value to which to apply modulo arithmetic
!! @param[in] xc Center of modulo range
!! @param[in] Lx Modulo range width
!! @return x_mod Value x shifted by an integer multiple of Lx to be close to xc
function modulo_around_point(x, xc, Lx) result(x_mod)
use gengrid_kinds, only : dbl_kind
real(dbl_kind), intent(in) :: x
real(dbl_kind), intent(in) :: xc
real(dbl_kind), intent(in) :: Lx
real(dbl_kind) :: x_mod
if (Lx > 0.0) then
x_mod = modulo(x - (xc - 0.5*Lx), Lx) + (xc - 0.5*Lx)
else
x_mod = x
endif
end function modulo_around_point
end module angles