Skip to content

Commit cae15a1

Browse files
GeorgeGayno-NOAAGeorgeGayno-NOAA
committed
More cleanup to interp2 and output2.
Fixes ufs-community#709.
1 parent 46db839 commit cae15a1

File tree

2 files changed

+104
-165
lines changed

2 files changed

+104
-165
lines changed

‎sorc/sfc_climo_gen.fd/interp2.F90

+97-159
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,9 @@ subroutine interp2(localpet, input_file)
1313

1414
use esmf
1515
use netcdf
16-
use model_grid
16+
use model_grid, only : grid_mdl, i_mdl, j_mdl, &
17+
num_tiles, latitude_field_mdl, &
18+
longitude_field_mdl, mask_field_mdl
1719
use source_grid
1820
use utils
1921
use mpi
@@ -23,33 +25,46 @@ subroutine interp2(localpet, input_file)
2325
character(len=*), intent(in) :: input_file
2426

2527
integer :: rc, localpet
26-
integer :: i, j, tile, n, ncid, status
27-
integer :: t
28-
integer :: clb_mdl(3), cub_mdl(3)
29-
integer :: varid, record
28+
integer :: i, j, tile, ncid, status
29+
integer :: varid
3030
integer :: isrctermprocessing
3131
integer :: category, num_categories
3232

3333
integer(esmf_kind_i4), allocatable :: mask_mdl_one_tile(:,:)
3434
integer(esmf_kind_i4), pointer :: unmapped_ptr(:)
3535

36-
real(esmf_kind_r4), pointer :: data_mdl_ptr(:,:,:)
3736
real(esmf_kind_r4), allocatable :: data_src_global(:,:)
3837
real(esmf_kind_r4), allocatable :: data_src_global2(:,:,:)
3938
real(esmf_kind_r4), allocatable :: data_mdl_one_tile(:,:,:)
40-
real(esmf_kind_r4), allocatable :: vegt_mdl_one_tile(:,:)
4139
real(esmf_kind_r4), allocatable :: lat_mdl_one_tile(:,:)
4240
real(esmf_kind_r4), allocatable :: sum_mdl_one_tile(:,:)
4341
real(esmf_kind_r4), allocatable :: lon_mdl_one_tile(:,:)
4442

4543
type(esmf_regridmethod_flag) :: method
4644
type(esmf_field) :: data_field_src
47-
type(esmf_field) :: data_field_mdl2
45+
type(esmf_field) :: data_field_mdl
4846
type(esmf_routehandle) :: regrid_data
4947
type(esmf_polemethod_flag) :: pole
5048

51-
! get this from file.
52-
num_categories = 20
49+
if (localpet == 0) then
50+
allocate(data_src_global(i_src,j_src))
51+
else
52+
allocate(data_src_global(0,0))
53+
endif
54+
55+
if (localpet == 0) then
56+
print*,'- OPEN SOURCE FILE ', trim(input_file)
57+
status = nf90_open(input_file, nf90_nowrite, ncid)
58+
call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE")
59+
status = nf90_inq_varid(ncid, field_names(1), varid)
60+
call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID")
61+
status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,1/), count=(/i_src,j_src,1/))
62+
call netcdf_err(status, "IN ROUTINE INTERP READING FIELD")
63+
print*,'number of cats ',maxval(data_src_global)
64+
num_categories = nint(maxval(data_src_global))
65+
endif
66+
67+
call mpi_bcast(num_categories,1,MPI_INTEGER,0,MPI_COMM_WORLD,status)
5368

5469
print*,"- CALL FieldCreate FOR SOURCE GRID DATA."
5570
data_field_src = ESMF_FieldCreate(grid_src, &
@@ -64,7 +79,7 @@ subroutine interp2(localpet, input_file)
6479
call error_handler("IN FieldCreate", rc)
6580

6681
print*,"- CALL FieldCreate FOR model GRID veg DATA."
67-
data_field_mdl2 = ESMF_FieldCreate(grid_mdl, &
82+
data_field_mdl = ESMF_FieldCreate(grid_mdl, &
6883
typekind=ESMF_TYPEKIND_R4, &
6984
indexflag=ESMF_INDEX_GLOBAL, &
7085
staggerloc=ESMF_STAGGERLOC_CENTER, &
@@ -75,85 +90,49 @@ subroutine interp2(localpet, input_file)
7590
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
7691
call error_handler("IN FieldCreate", rc)
7792

78-
print*,"- CALL FieldGet FOR MODEL GRID DATA."
79-
nullify(data_mdl_ptr)
80-
call ESMF_FieldGet(data_field_mdl2, &
81-
farrayPtr=data_mdl_ptr, &
82-
computationalLBound=clb_mdl, &
83-
computationalUBound=cub_mdl, &
84-
rc=rc)
85-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
86-
call error_handler("IN FieldGet", rc)
87-
88-
89-
print*,'got here2 ',localpet,clb_mdl,cub_mdl
90-
9193
if (localpet == 0) then
92-
allocate(data_src_global(i_src,j_src))
9394
allocate(data_src_global2(i_src,j_src,num_categories))
94-
else
95-
allocate(data_src_global(0,0))
96-
allocate(data_src_global2(0,0,0))
97-
endif
98-
99-
print*,'- OPEN SOURCE FILE ', trim(input_file)
100-
status = nf90_open(input_file, nf90_nowrite, ncid)
101-
call netcdf_err(status, "IN ROUTINE INTERP OPENING SOURCE FILE")
102-
103-
if (localpet == 0) then
10495
allocate(data_mdl_one_tile(i_mdl,j_mdl,num_categories))
10596
allocate(mask_mdl_one_tile(i_mdl,j_mdl))
10697
allocate(lat_mdl_one_tile(i_mdl,j_mdl))
10798
allocate(sum_mdl_one_tile(i_mdl,j_mdl))
10899
allocate(lon_mdl_one_tile(i_mdl,j_mdl))
109100
else
101+
allocate(data_src_global2(0,0,0))
110102
allocate(data_mdl_one_tile(0,0,0))
111103
allocate(mask_mdl_one_tile(0,0))
112104
allocate(lat_mdl_one_tile(0,0))
113105
allocate(sum_mdl_one_tile(0,0))
114106
allocate(lon_mdl_one_tile(0,0))
115107
endif
116108

117-
record = 0
118-
119-
TIME_LOOP : do t = 1, num_time_recs ! loop over each time period
120-
121-
FIELD_LOOP : do n = 1, num_fields ! loop over each surface field.
122-
123-
record = record + 1
124-
125-
if (localpet == 0) then
126-
status = nf90_inq_varid(ncid, field_names(n), varid)
127-
call netcdf_err(status, "IN ROUTINE INTERP READING FIELD ID")
128-
status = nf90_get_var(ncid, varid, data_src_global, start=(/1,1,t/), count=(/i_src,j_src,1/))
129-
call netcdf_err(status, "IN ROUTINE INTERP READING FIELD")
130-
data_src_global2 = 0.0
131-
do j = 1, j_src
132-
do i = 1, i_src
133-
category = nint(data_src_global(i,j))
134-
! if (category < 1) category = 17
135-
if (category < 1) cycle
136-
data_src_global2(i,j,category) = 1.0
137-
enddo
138-
enddo
139-
endif
109+
if (localpet == 0) then
110+
data_src_global2 = 0.0
111+
do j = 1, j_src
112+
do i = 1, i_src
113+
category = nint(data_src_global(i,j))
114+
if (category < 1) cycle
115+
data_src_global2(i,j,category) = 1.0
116+
enddo
117+
enddo
118+
endif
140119

141-
print*,"- CALL FieldScatter FOR SOURCE GRID DATA."
142-
call ESMF_FieldScatter(data_field_src, data_src_global2, rootPet=0, rc=rc)
143-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
144-
call error_handler("IN FieldScatter.", rc)
120+
deallocate(data_src_global)
145121

146-
isrctermprocessing = 1
122+
print*,"- CALL FieldScatter FOR SOURCE GRID DATA."
123+
call ESMF_FieldScatter(data_field_src, data_src_global2, rootPet=0, rc=rc)
124+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
125+
call error_handler("IN FieldScatter.", rc)
147126

148-
if (record == 1) then
127+
isrctermprocessing = 1
149128

150-
method = ESMF_REGRIDMETHOD_CONSERVE
151-
pole = ESMF_POLEMETHOD_NONE
129+
method = ESMF_REGRIDMETHOD_CONSERVE
130+
pole = ESMF_POLEMETHOD_NONE
152131

153-
print*,"- CALL FieldRegridStore."
154-
nullify(unmapped_ptr)
155-
call ESMF_FieldRegridStore(data_field_src, &
156-
data_field_mdl2, &
132+
print*,"- CALL FieldRegridStore."
133+
nullify(unmapped_ptr)
134+
call ESMF_FieldRegridStore(data_field_src, &
135+
data_field_mdl, &
157136
srcmaskvalues=(/0/), &
158137
dstmaskvalues=(/0/), &
159138
polemethod=pole, &
@@ -164,111 +143,70 @@ subroutine interp2(localpet, input_file)
164143
regridmethod=method, &
165144
unmappedDstList=unmapped_ptr, &
166145
rc=rc)
167-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
168-
call error_handler("IN FieldRegridStore.", rc)
146+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
147+
call error_handler("IN FieldRegridStore.", rc)
148+
149+
print*,"- CALL Field_Regrid."
150+
call ESMF_FieldRegrid(data_field_src, &
151+
data_field_mdl, &
152+
routehandle=regrid_data, &
153+
termorderflag=ESMF_TERMORDER_SRCSEQ, &
154+
rc=rc)
155+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
156+
call error_handler("IN FieldRegrid.", rc)
169157

170-
endif
158+
OUTPUT_LOOP : do tile = 1, num_tiles
171159

172-
print*,"- CALL Field_Regrid."
173-
call ESMF_FieldRegrid(data_field_src, &
174-
data_field_mdl2, &
175-
routehandle=regrid_data, &
176-
termorderflag=ESMF_TERMORDER_SRCSEQ, &
177-
rc=rc)
160+
print*,"- CALL FieldGather FOR MODEL LATITUDE."
161+
call ESMF_FieldGather(latitude_field_mdl, lat_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
178162
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
179-
call error_handler("IN FieldRegrid.", rc)
180-
181-
!-----------------------------------------------------------------------
182-
! Unmapped points are stored in "unmapped_ptr". The pointer contains
183-
! "ij" global indices as follows:
184-
!
185-
! tile 1: 1 thru (itile*jtile)
186-
! tile n: (n-1)*(itile*jtile) thru n*(itile*jtile)
187-
!
188-
! This "ij" index is converted to the tile number and i/j index of the
189-
! field object. This logic assumes the model grid object was
190-
! created using "GLOBAL" indices.
191-
!
192-
! Unmapped data points are given the flag value of -9999.9, which
193-
! will be replaced in routine "search".
194-
!-----------------------------------------------------------------------
195-
196-
OUTPUT_LOOP : do tile = 1, num_tiles
197-
198-
print*,"- CALL FieldGather FOR MODEL LATITUDE."
199-
call ESMF_FieldGather(latitude_field_mdl, lat_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
200-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
201-
call error_handler("IN FieldGather.", rc)
202-
203-
print*,"- CALL FieldGather FOR MODEL LONGITUDE."
204-
call ESMF_FieldGather(longitude_field_mdl, lon_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
205-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
206-
call error_handler("IN FieldGather.", rc)
207-
208-
print*,"- CALL FieldGather FOR MODEL GRID DATA."
209-
call ESMF_FieldGather(data_field_mdl2, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
210-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
211-
call error_handler("IN FieldGather.", rc)
212-
213-
print*,"- CALL FieldGather FOR MODEL GRID LAND MASK."
214-
call ESMF_FieldGather(mask_field_mdl, mask_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
215-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
216-
call error_handler("IN FieldGather.", rc)
217-
218-
select case (trim(field_names(n)))
219-
case ('substrate_temperature','vegetation_greenness','leaf_area_index','slope_type','soil_type')
220-
print*,"- CALL FieldGather FOR MODEL GRID VEG TYPE."
221-
call ESMF_FieldGather(vegt_field_mdl, vegt_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
222-
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
223-
call error_handler("IN FieldGather.", rc)
224-
end select
225-
226-
if (localpet == 0) then
227-
print*,'- CALL SEARCH FOR TILE ',tile
228-
sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3) ! use unused variable to now.
229-
do j = 1, j_mdl
230-
do i = 1, i_mdl
231-
232-
if (mask_mdl_one_tile(i,j) == 1 .and. sum_mdl_one_tile(i,j) == 0.0) then
233-
data_mdl_one_tile(i,j,:) = -9999.9
234-
endif
235-
236-
enddo
237-
enddo
163+
call error_handler("IN FieldGather.", rc)
238164

239-
call search2 (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, field_names(n))
240-
! where(mask_mdl_one_tile == 0) data_mdl_one_tile = missing
241-
print*,'after regrid ',data_mdl_one_tile(i_mdl/2,j_mdl/2,:)
242-
call output2 (data_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, t, n)
243-
endif
244-
245-
print*,'after output ', localpet
246-
call mpi_barrier(mpi_comm_world, rc)
247-
stop
165+
print*,"- CALL FieldGather FOR MODEL LONGITUDE."
166+
call ESMF_FieldGather(longitude_field_mdl, lon_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
167+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
168+
call error_handler("IN FieldGather.", rc)
248169

249-
! if (field_names(n) == 'vegetation_type') then
250-
! print*,"- CALL FieldScatter FOR MODEL GRID VEGETATION TYPE."
251-
! call ESMF_FieldScatter(vegt_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
252-
! if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
253-
! call error_handler("IN FieldScatter.", rc)
254-
! endif
170+
print*,"- CALL FieldGather FOR MODEL GRID DATA."
171+
call ESMF_FieldGather(data_field_mdl, data_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
172+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
173+
call error_handler("IN FieldGather.", rc)
255174

256-
enddo OUTPUT_LOOP
175+
print*,"- CALL FieldGather FOR MODEL GRID LAND MASK."
176+
call ESMF_FieldGather(mask_field_mdl, mask_mdl_one_tile, rootPet=0, tile=tile, rc=rc)
177+
if(ESMF_logFoundError(rcToCheck=rc,msg=ESMF_LOGERR_PASSTHRU,line=__LINE__,file=__FILE__)) &
178+
call error_handler("IN FieldGather.", rc)
257179

258-
if (allocated(vegt_mdl_one_tile)) deallocate(vegt_mdl_one_tile)
180+
if (localpet == 0) then
181+
print*,'- CALL SEARCH FOR TILE ',tile
182+
sum_mdl_one_tile = sum(data_mdl_one_tile, dim=3) ! use unused variable to now.
183+
do j = 1, j_mdl
184+
do i = 1, i_mdl
185+
if (mask_mdl_one_tile(i,j) == 1 .and. sum_mdl_one_tile(i,j) == 0.0) then
186+
data_mdl_one_tile(i,j,:) = -9999.9
187+
endif
188+
enddo
189+
enddo
190+
call search2 (data_mdl_one_tile, mask_mdl_one_tile, i_mdl, j_mdl, num_categories, tile, field_names(1))
191+
print*,'after regrid ',data_mdl_one_tile(i_mdl/2,j_mdl/2,:)
192+
call output2 (data_mdl_one_tile, lat_mdl_one_tile, lon_mdl_one_tile, i_mdl, j_mdl, num_categories, tile)
193+
endif
259194

260-
enddo FIELD_LOOP
261-
enddo TIME_LOOP
195+
enddo OUTPUT_LOOP
262196

263197
status=nf90_close(ncid)
264198

265-
deallocate(data_mdl_one_tile, mask_mdl_one_tile)
266-
deallocate(data_src_global, lat_mdl_one_tile, lon_mdl_one_tile)
199+
deallocate(data_mdl_one_tile, mask_mdl_one_tile, data_src_global2)
200+
deallocate(lat_mdl_one_tile, lon_mdl_one_tile, sum_mdl_one_tile)
267201

268202
print*,"- CALL FieldRegridRelease."
269203
call ESMF_FieldRegridRelease(routehandle=regrid_data, rc=rc)
270204

271205
print*,"- CALL FieldDestroy."
272206
call ESMF_FieldDestroy(data_field_src, rc=rc)
273207

208+
print*,'after output ', localpet
209+
call mpi_barrier(mpi_comm_world, rc)
210+
stop
211+
274212
end subroutine interp2

0 commit comments

Comments
 (0)