cpld_gridgen 1.14.0
Loading...
Searching...
No Matches
angles.F90
Go to the documentation of this file.
1
9
10module angles
11
12 use gengrid_kinds, only : dbl_kind, int_kind
13 use grdvars, only : debug
14
15 implicit none
16
17contains
26
27 subroutine find_angq(iind,jind,xangCt,anglet,angle)
28
29 integer , intent(in) :: iind(2),jind(2)
30 real(dbl_kind), intent(in) :: xangCt(:)
31 real(dbl_kind), intent(in) :: anglet(:,:)
32 real(dbl_kind), intent(out) :: angle(:,:)
33
34 ! local variables
35 integer :: i,j
36
37 real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
38 real(dbl_kind) :: p25 = 0.25
39
40 !---------------------------------------------------------------------
41 ! find the angle on corners using the same relationship CICE uses
42 ! internally to calculate angles on Ct using angles on Bu
43 !
44 ! w-----------------0 Ct(i+1,j+1)
45 ! | |
46 ! ----------Bu(i,j)---------- Bu lies on seam at j=nj
47 ! | |
48 ! Ct(i,j) sw----------------s
49 !
50 !---------------------------------------------------------------------
51
52 angle = 0.0
53 do j = jind(1)+1,jind(2)
54 do i = iind(1),iind(2)-1
55 if (j .lt. jind(2)) then
56 angle_0 = anglet(i+1,j+1)
57 angle_w = anglet(i, j+1)
58 angle_s = anglet(i+1,j )
59 angle_sw = anglet(i ,j )
60 else
61 angle_0 = xangct(i+1 )
62 angle_w = xangct(i )
63 angle_s = anglet(i+1,j)
64 angle_sw = anglet(i, j)
65 end if
66 angle(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
67 p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))
68
69 if (abs(angle(i,j)) .le. 1.0e-10)angle(i,j) = 0.0
70 enddo
71 enddo
72
73 end subroutine find_angq
74
83 subroutine find_angchk(iind,jind,angle,angchk)
84
85 integer , intent(in) :: iind(2),jind(2)
86 real(dbl_kind), intent(in) :: angle(:,:)
87 real(dbl_kind), intent(out) :: angchk(:,:)
88
89 ! local variables
90 integer :: i,j
91 real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
92 real(dbl_kind) :: p25 = 0.25
93
94 !---------------------------------------------------------------------
95 ! check: calculate anglet from angle on corners as CICE does internally.
96 ! since angle changes sign between CICE and MOM6, (-1)*angchk ~ anglet
97 !
98 ! w-----------------0 Bu(i,j)
99 ! | |
100 ! | Ct(i,j) |
101 ! | |
102 ! Bu(i-1,j-1) sw----------------s
103 !
104 !---------------------------------------------------------------------
105
106 angchk = 0.0
107 do j = jind(1)+1,jind(2)
108 do i = iind(1)+1,iind(2)
109 angle_0 = angle(i ,j )
110 angle_w = angle(i-1,j )
111 angle_s = angle(i, j-1)
112 angle_sw = angle(i-1,j-1)
113 angchk(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
114 p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))
115 enddo
116 enddo
117
118 end subroutine find_angchk
119
129
130 subroutine find_ang(iind,jind,lonBu,latBu,lonCt,anglet)
131
132 integer , intent(in) :: iind(2),jind(2)
133 real(dbl_kind), intent(in) :: lonBu(:,:)
134 real(dbl_kind), intent(in) :: latBu(:,:)
135 real(dbl_kind), intent(in) :: lonCt(:,:)
136 real(dbl_kind), intent(out) :: anglet(:,:)
137
138 ! local variables
139 integer :: i,j,m,n
140 integer :: ii,jj
141
142 ! from geolonB fix in MOM6
143 real(dbl_kind) :: len_lon ! The periodic range of longitudes, usually 360 degrees.
144 real(dbl_kind) :: pi_720deg ! One quarter the conversion factor from degrees to radians.
145 real(dbl_kind) :: lonB(2,2) ! The longitude of a point, shifted to have about the same value.
146 real(dbl_kind) :: lon_scale = 0.0
147
148 !---------------------------------------------------------------------
149 ! rotation angle for "use_bugs" = false case from MOM6
150 ! src/initialization/MOM_shared_initialization.F90 but allow for not
151 ! having halo values
152 ! note this does not reproduce sin_rot,cos_rot found in MOM6 output
153 ! differences are ~O 10-6
154 !---------------------------------------------------------------------
155
156 anglet = 0.0
157 pi_720deg = atan(1.0) / 180.0
158 len_lon = 360.0
159 do j=jind(1),jind(2); do i = iind(1),iind(2)
160 do n=1,2 ; do m=1,2
161 jj = j+n-2; ii = i+m-2
162 if(jj .eq. 0)jj = 1
163 if(ii .eq. 0)ii = iind(2)
164 lonb(m,n) = modulo_around_point(lonbu(ii,jj), lonct(i,j), len_lon)
165 ! lonB(m,n) = modulo_around_point(LonBu(I+m-2,J+n-2), LonCt(i,j), len_lon)
166 enddo; enddo
167 jj = j-1; ii = i-1
168 if(jj .eq. 0)jj = 1
169 if(ii .eq. 0)ii = iind(2)
170 lon_scale = cos(pi_720deg*((latbu(ii,jj) + latbu(i,j)) + &
171 (latbu(i,jj) + latbu(ii,j)) ) )
172 anglet(i,j) = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
173 (latbu(ii,j) - latbu(i,jj)) + &
174 (latbu(i,j) - latbu(ii,jj)) )
175
176 !lon_scale = cos(pi_720deg*((LatBu(I-1,J-1) + LatBu(I,J)) + &
177 ! (LatBu(I,J-1) + LatBu(I-1,J)) ) )
178 !anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
179 ! (LatBu(I-1,J) - LatBu(I,J-1)) + &
180 ! (LatBu(I,J) - LatBu(I-1,J-1)) )
181 enddo; enddo
182
183 end subroutine find_ang
184 ! -----------------------------------------------------------------------------
193 function modulo_around_point(x, xc, Lx) result(x_mod)
194 use gengrid_kinds, only : dbl_kind
195
196 real(dbl_kind), intent(in) :: x
197 real(dbl_kind), intent(in) :: xc
198 real(dbl_kind), intent(in) :: lx
199 real(dbl_kind) :: x_mod
200
201 if (lx > 0.0) then
202 x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
203 else
204 x_mod = x
205 endif
206 end function modulo_around_point
207end module angles