-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathsmooth_topo_cube.F90
342 lines (264 loc) · 10.6 KB
/
smooth_topo_cube.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
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
!-----------------------------------------------------------------------------
! MODULE subgrid_topo_ana
!
! Purpose:
!
! Date Programmer Affiliation Description of change
! ==== ========== =========== =====================
!
!-----------------------------------------------------------------------------
#define WRITERESULTS
MODULE smooth_topo_cube_sph
USE reconstruct
!USE ridge_ana
IMPLICIT NONE
PRIVATE
PUBLIC smooth_intermediate_topo
CONTAINS
!=============================================================================
SUBROUTINE smooth_intermediate_topo(terr, da, ncube,nhalo, NSCL_f,NSCL_c, SMITER &
, terr_sm,terr_dev, lread_smooth_topofile )
REAL (KIND=dbl_kind), PARAMETER :: pi = 3.14159265358979323846264338327
INTEGER (KIND=int_kind), INTENT(IN) :: ncube, nhalo,NSCL_f,NSCL_c,SMITER
REAL (KIND=dbl_kind), &
DIMENSION(ncube,ncube,6), INTENT(INOUT) :: terr
REAL (KIND=dbl_kind), &
DIMENSION(ncube,ncube), INTENT(INOUT) :: da
#if 1
REAL (KIND=dbl_kind), &
DIMENSION(ncube,ncube,6), INTENT(INOUT) :: terr_dev
REAL (KIND=dbl_kind), &
DIMENSION(ncube,ncube,6), INTENT(INOUT) :: terr_sm
#endif
LOGICAL, INTENT(IN) :: lread_smooth_topofile
!-----------------------------------------------------------------
!PRIMARY Outputs
!-----------------------------------------------------------------
!------------------------------
REAL (KIND=dbl_kind), &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: terr_halo
REAL (KIND=dbl_kind), &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: terr_halo_sm
REAL (KIND=dbl_kind), &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: da_halo
REAL (KIND=dbl_kind) , &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: terr_halo_r4
REAL (KIND=dbl_kind) , &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: terr_halo_rw
REAL (KIND=dbl_kind) , &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo, 6) :: terr_halo_fx,terr_halo_fx_sv
REAL (KIND=dbl_kind) , &
DIMENSION(1-nhalo:ncube+nhalo, 1-nhalo:ncube+nhalo) :: smwt,ggaa,ggbb,ggab
REAL (KIND=dbl_kind) , &
DIMENSION(1-nhalo:ncube+nhalo ) :: xv,yv,alph,beta
INTEGER (KIND=int_kind):: np,i,j, ncube_halo,norx,nory,ipanel,x0,&
x1,y0,y1,initd,ii0,ii1,jj0,jj1,nctest,NSM,NS2,ismi,NSB
REAL (KIND=dbl_kind), allocatable :: daxx(:,:,:)
REAL (KIND=dbl_kind), allocatable :: wt1p(:,:),terr_patch(:,:)
REAL(KIND=dbl_kind) :: cosll, dx, dy ,dbet,dalp,diss,diss00,lon_ij,lat_ij,latfactor
INTEGER :: NOCTV , isx0, isx1, jsy0, jsy1,i2,j2,iix,jjx,i00,ncube_in_file
REAL(KIND=dbl_kind) :: RSM_scl, smoo,irho,volt0,volt1,volume_after,volume_before
CHARACTER(len=1024) :: ofile$
logical :: read_in_precomputed
!read_in_precomputed = .FALSE.
read_in_precomputed = lread_smooth_topofile !.TRUE.
IF (read_in_precomputed) then
write( ofile$ , &
"('./output/topo_smooth_nc',i0.4,'_Co',i0.3,'_Fi',i0.3)" ) &
ncube, NSCL_c/2, NSCL_f/2
ofile$= trim(ofile$)//'.dat'
OPEN (unit = 711, file= trim(ofile$) ,form="UNFORMATTED" )
READ(711) ncube_in_file
READ(711) terr
READ(711) terr_sm
READ(711) terr_dev
close(711)
write(*,*) " Read precomputed filtered topography from "
write(*,*) ofile$
RETURN
ENDIF
write(*,*) " NCUBE !!! " , ncube
write( ofile$ , &
"('./output/topo_smooth_nc',i0.4,'_Co',i0.3,'_Fi',i0.3)" ) &
ncube, NSCL_c/2, NSCL_f/2
#ifdef USELATFAC
ofile$= trim(ofile$)//'_LATFACTOR.dat'
#else
ofile$= trim(ofile$)//'.dat'
#endif
write(*,*) " Will do smoothing of topo on cubed sphere "
write(*,*) " Output will go to:"
write(*,*) ofile$
allocate( daxx(ncube,ncube,6) )
DO np = 1, 6
daxx(:,:,np) = da
end do
DO np = 1, 6
!CALL CubedSphereFillHalo_Linear_extended(terr, terr_halo(:,:,np), np, ncube+1,nhalo)
!CALL CubedSphereFillHalo_Linear_extended(daxx, da_halo(:,:,np), np, ncube+1,nhalo)
CALL CubedSphereFillHalo(terr, terr_halo, np, ncube+1,nhalo)
CALL CubedSphereFillHalo(daxx, da_halo, np, ncube+1,nhalo)
END DO
deallocate( daxx )
ncube_halo = size( terr_halo(:,:,1), 1 )
!terr_halo(1-nhalo:1,1-nhalo:1,:)=0.
!terr_halo(ncube:ncube+nhalo,ncube:ncube+nhalo , :)=0.
DO i=1-nhalo,ncube+nhalo
xv(i)=1.*i
yv(i)=1.*i
END DO
DO i=1-nhalo,ncube+nhalo
alph(i)=(pi/4.)*(1.*i - 0.5 + nhalo - (ncube+2.*nhalo)/2.) / ((ncube+2.*nhalo)/2.)
beta(i)=(pi/4.)*(1.*i - 0.5 + nhalo - (ncube+2.*nhalo)/2.) / ((ncube+2.*nhalo)/2.)
END DO
DO j=1-nhalo,ncube+nhalo
DO i=1-nhalo,ncube+nhalo
irho = ( 1. + (tan(alph(i))**2) + (tan(beta(j))**2 ) )**2
irho = 1. / ( ( cos(alph(i))**2 ) * (cos(beta(j))**2) * irho )
!irho = 1./ ( ( cos(alph(i))**2)*(cos(beta(j))**2)* ( ( 1. + (tan(alph(i))**2) + (tan(beta(j))**2 ) )**2 ))
ggaa(i,j) = irho * ( 1. + ( tan( alph(i) ) )**2 )
ggbb(i,j) = irho * ( 1. + ( tan( beta(j) ) )**2 )
ggab(i,j) = -irho *( tan( beta(j) ) ) * ( tan( alph(i) ) )
END DO
END DO
terr_halo_sm = terr_halo
terr_halo_fx = terr_halo
terr_halo_rw = terr_halo
if (NSCL_f > 1 ) then
NSM=NSCL_f
NS2=NSM/2
write(*,*)" smoothing fine scle w/" ,NSCL_f
allocate( wt1p(-ns2:ns2, -ns2:ns2 ) )
allocate( terr_patch(-ns2:ns2, -ns2:ns2 ) )
i00 = ncube/2
dalp = alph(i00+ns2 )-alph(i00)
diss00 = 1./sqrt( ggaa(i00,i00)*dalp*dalp )
terr_halo_fx = 0.0
#ifdef DEBUGRIDGE
DO np=4,4
DO j=2500,3000
DO i=550,1350
#else
DO np=1,6
DO j=1-nhalo+ns2,ncube+nhalo-ns2
DO i=1-nhalo+ns2,ncube+nhalo-ns2
#endif
volt0 = terr_halo(i,j,np)*da_halo(i,j,np)
volt1 = 0.
do j2=-ns2,ns2
do i2=-ns2,ns2
jjx = j+j2
iix = i+i2
dalp = alph(iix)-alph(i)
dbet = beta(jjx)-beta(j)
diss = ggaa(i,j)*dalp*dalp + ggbb(i,j)*dbet*dbet + 2.*ggab(i,j)*dalp*dbet
wt1p(i2,j2) = da_halo(iix,jjx,np)
terr_patch(i2,j2) = terr_halo(i,j,np)*( 1. - diss00 * sqrt( diss ) ) !*da_halo(iix,jjx,np)
if ((volt0*terr_patch(i2,j2)<=0.).or.(wt1p(i2,j2)<=0.) ) then
terr_patch(i2,j2)=0.
wt1p(i2,j2) =0.
end if
volt1 = volt1 + terr_patch(i2,j2)*wt1p(i2,j2)
end do
end do
if ( abs(volt1) > 0.) terr_patch = (volt0 / volt1) * terr_patch
do j2=-ns2,ns2
do i2=-ns2,ns2
jjx = j+j2
iix = i+i2
terr_halo_fx(iix,jjx,np) = terr_halo_fx(iix,jjx,np) + terr_patch(i2,j2)
end do
end do
END DO
END DO
END DO
deallocate( wt1p )
deallocate( terr_patch )
else
write(*,*)" No fine scale smoother "
terr_halo_fx = terr_halo
endif
NSM=NSCL_c
NS2=NSM/2
allocate( wt1p(-ns2:ns2, -ns2:ns2 ) )
allocate( terr_patch(-ns2:ns2, -ns2:ns2 ) )
i00 = ncube/2
dalp = alph(i00+ns2 )-alph(i00)
diss00 = 1./sqrt( ggaa(i00,i00)*dalp*dalp )
!terr_halo_sm = 0.0
terr_halo_fx_sv = terr_halo_fx
write(*,*) "LIMITS in smoother "
write(*,*) 1-nhalo+ns2,ncube+nhalo-ns2
do ismi = 1,SMITER
terr_halo_sm = 0.0
#ifdef DEBUGRIDGE
DO np=4,4
DO j=2500,3000
DO i=300,1100
#else
DO np=1,6
DO j=1-nhalo+ns2,ncube+nhalo-ns2
DO i=1-nhalo+ns2,ncube+nhalo-ns2
#endif
volt0 = terr_halo_fx(i,j,np)*da_halo(i,j,np)
volt1 = 0.
#ifdef USELATFAC
call CubedSphereRLLFromABP(alph(i), beta(j) , np , lon_ij, lat_ij ) ! Results in radians
latfactor = 1. / cos( lat_ij )
latfactor = min(latfactor, 3.0 )
#endif
do j2=-ns2,ns2
do i2=-ns2,ns2
jjx = j+j2
iix = i+i2
dalp = alph(iix)-alph(i)
dbet = beta(jjx)-beta(j)
diss = ggaa(i,j)*dalp*dalp + ggbb(i,j)*dbet*dbet + 2.*ggab(i,j)*dalp*dbet
wt1p(i2,j2) = da_halo(iix,jjx,np)
#ifdef USELATFAC
terr_patch(i2,j2) = terr_halo_fx(i,j,np)*( 1. - latfactor * diss00 * sqrt( diss ) ) !*da_halo(iix,jjx,np)
#else
terr_patch(i2,j2) = terr_halo_fx(i,j,np)*( 1. - diss00 * sqrt( diss ) ) !*da_halo(iix,jjx,np)
#endif
if ((volt0*terr_patch(i2,j2)<=0.).or.(wt1p(i2,j2)<=0.) ) then
terr_patch(i2,j2)=0.
wt1p(i2,j2) =0.
end if
volt1 = volt1 + terr_patch(i2,j2)*wt1p(i2,j2)
end do
end do
if ( abs(volt1) > 0.) terr_patch = (volt0 / volt1) * terr_patch
do j2=-ns2,ns2
do i2=-ns2,ns2
jjx = j+j2
iix = i+i2
terr_halo_sm(iix,jjx,np) = terr_halo_sm(iix,jjx,np) + terr_patch(i2,j2)
end do
end do
END DO
if (mod(j,1) ==0 ) write(*,*) "Crs Sm J = ",J, " Panel=",np," iter=",ismi
END DO
END DO
!terr_halo_fx = terr_halo_sm
do np=1,6
terr_sm (1:ncube,1:ncube,np) = terr_halo_sm(1:ncube,1:ncube,np )
end do
do np=1,6
CALL CubedSphereFillHalo(terr_sm, terr_halo_fx, np, ncube+1,nhalo)
end do
end do
deallocate( wt1p )
deallocate( terr_patch )
do np=1,6
terr_dev(1:ncube,1:ncube,np) = terr_halo_fx_sv(1:ncube,1:ncube,np ) - terr_halo_sm(1:ncube,1:ncube,np )
terr_sm (1:ncube,1:ncube,np) = terr_halo_sm(1:ncube,1:ncube,np )
end do
OPEN (unit = 711, file= trim(ofile$) ,form="UNFORMATTED" )
write(711) ncube
WRITE(711) terr
WRITE(711) terr_sm
WRITE(711) terr_dev
close(711)
STOP
END SUBROUTINE smooth_intermediate_topo
END MODULE smooth_topo_cube_sph