Line data Source code
1 : module getinterpnetcdfdata
2 :
3 : ! Description:
4 : ! Routines for extracting a column from a netcdf file
5 : !
6 : ! Author:
7 : !
8 : ! Modules Used:
9 : !
10 : use cam_abortutils, only: endrun
11 : use pmgrid, only: plev
12 : use cam_logfile, only: iulog
13 :
14 : implicit none
15 : private
16 : !
17 : ! Public Methods:
18 : !
19 : public getinterpncdata
20 :
21 : contains
22 :
23 0 : subroutine getinterpncdata( NCID, camlat, camlon, TimeIdx, &
24 : varName, have_surfdat, surfdat, fill_ends, scm_crm_mode, &
25 0 : press, npress, ps, hyam, hybm, outData, STATUS )
26 :
27 : ! getinterpncdata: extracts the entire level dimension for a
28 : ! particular lat,lon,time from a netCDF file
29 : ! and interpolates it onto the input pressure levels, placing
30 : ! result in outData, and the error status inx STATUS
31 :
32 : use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
33 : use shr_scam_mod, only: shr_scam_GetCloseLatLon
34 : use netcdf
35 : implicit none
36 : !-----------------------------------------------------------------------
37 :
38 :
39 : ! ---------- inputs ------------
40 :
41 : integer, intent(in) :: NCID ! NetCDF ID
42 : integer, intent(in) :: TimeIdx ! time index
43 : real(r8), intent(in) :: camlat,camlon ! target lat and lon to be extracted
44 : logical, intent(in) :: have_surfdat ! is surfdat provided
45 : logical, intent(in) :: fill_ends ! extrapolate the end values
46 : logical, intent(in) :: scm_crm_mode ! scam column radiation mode
47 : integer, intent(in) :: npress ! number of dataset pressure levels
48 : real(r8), intent(in) :: press(npress) ! dataset pressure levels
49 : real(r8), intent(in) :: ps ! surface pressure
50 : real(r8), intent(in) :: hyam(:) ! dataset hybrid midpoint pressure levels
51 : real(r8), intent(in) :: hybm(:) ! dataset hybrid midpoint pressure levels
52 :
53 : ! ---------- outputs ----------
54 :
55 : real(r8), intent(inout) :: outData(:) ! interpolated output data
56 : integer, intent(out) :: STATUS ! return status of netcdf calls
57 :
58 : ! ------- locals ---------
59 :
60 : real(r8) surfdat ! surface value to be added before interpolation
61 : integer nlev ! number of levels in dataset
62 : integer latIdx ! latitude index
63 : integer lonIdx ! longitude index
64 0 : real(r8),allocatable :: tmp(:)
65 : real(r8) closelat,closelon
66 : real(r8) dx, dy, m ! slope for interpolation of surfdat
67 : integer varID
68 : integer var_ndims
69 : integer dims_set
70 : integer i
71 : integer var_dimIDs( NF90_MAX_VAR_DIMS )
72 : integer start( NF90_MAX_VAR_DIMS )
73 : integer count( NF90_MAX_VAR_DIMS )
74 :
75 : character varName*(*)
76 : character dim_name*( 256 )
77 : real(r8) missing_val
78 : logical usable_var
79 :
80 : ! ------- code ---------
81 :
82 0 : call shr_scam_GetCloseLatLon(ncid,camlat,camlon,closelat,closelon,latidx,lonidx)
83 :
84 : !
85 : ! Check mode: double or single precision
86 : !
87 :
88 : !
89 : ! Get var ID. Check to make sure it's there.
90 : !
91 0 : STATUS = NF90_INQ_VARID( NCID, varName, varID )
92 :
93 0 : if ( STATUS .NE. NF90_NOERR ) return
94 :
95 : !
96 : ! Check the var variable's information with what we are expecting
97 : ! it to be.
98 : !
99 :
100 0 : STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, ndims=var_ndims )
101 0 : if ( var_ndims .GT. 4 ) then
102 0 : write(iulog,* ) 'ERROR - extractdata.F: The input var',varName, &
103 0 : 'has', var_ndims, 'dimensions'
104 0 : STATUS = -1
105 : endif
106 :
107 : !
108 : ! surface variables
109 : !
110 0 : if ( var_ndims .EQ. 0 ) then
111 0 : STATUS = NF90_GET_VAR( NCID, varID, outData )
112 0 : return
113 : endif
114 :
115 0 : STATUS = NF90_INQUIRE_VARIABLE( NCID, varID, dimids=var_dimIDs )
116 0 : if ( STATUS .NE. NF90_NOERR ) then
117 0 : write(iulog,* ) 'ERROR - extractdata.F:Cant get dimension IDs for', varName
118 0 : return
119 : endif
120 : !
121 : ! Initialize the start and count arrays
122 : !
123 0 : dims_set = 0
124 0 : nlev = 1
125 0 : do i = var_ndims, 1, -1
126 :
127 0 : usable_var = .false.
128 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), dim_name )
129 :
130 0 : if ( dim_name .EQ. 'lat' ) then
131 0 : start( i ) = latIdx
132 0 : count( i ) = 1 ! Extract a single value
133 0 : dims_set = dims_set + 1
134 0 : usable_var = .true.
135 : endif
136 :
137 0 : if ( dim_name .EQ. 'lon' .or. dim_name .EQ. 'ncol' .or. dim_name .EQ. 'ncol_d' ) then
138 0 : start( i ) = lonIdx
139 0 : count( i ) = 1 ! Extract a single value
140 0 : dims_set = dims_set + 1
141 0 : usable_var = .true.
142 : endif
143 :
144 0 : if ( dim_name .EQ. 'lev' ) then
145 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
146 0 : start( i ) = 1
147 0 : count( i ) = nlev ! Extract all levels
148 0 : dims_set = dims_set + 1
149 0 : usable_var = .true.
150 : endif
151 :
152 0 : if ( dim_name .EQ. 'ilev' ) then
153 0 : STATUS = NF90_INQUIRE_DIMENSION( NCID, var_dimIDs( i ), len=nlev )
154 0 : start( i ) = 1
155 0 : count( i ) = nlev ! Extract all levels
156 0 : dims_set = dims_set + 1
157 0 : usable_var = .true.
158 : endif
159 :
160 0 : if ( dim_name .EQ. 'time' .OR. dim_name .EQ. 'tsec' ) then
161 0 : start( i ) = TimeIdx
162 0 : count( i ) = 1 ! Extract a single value
163 0 : dims_set = dims_set + 1
164 : usable_var = .true.
165 : endif
166 :
167 0 : if ( usable_var .EQV. .false. ) then
168 0 : write(iulog,* )'ERROR - extractdata.F: The input var ',varName, &
169 0 : ' has an unusable dimension ', dim_name
170 0 : STATUS = 1
171 : endif
172 : end do
173 :
174 0 : if ( dims_set .NE. var_ndims ) then
175 0 : write(iulog,* )'ERROR - extractdata.F: Could not find all the', &
176 0 : ' dimensions for input var ', varName
177 0 : write(iulog,* )'Found ',dims_set, ' of ',var_ndims
178 0 : STATUS = 1
179 : endif
180 :
181 0 : allocate(tmp(nlev+1))
182 :
183 0 : STATUS = NF90_GET_VAR( NCID, varID, tmp, start, count )
184 :
185 0 : if ( STATUS .NE. NF90_NOERR ) then
186 0 : write(iulog,* )'ERROR - extractdata.F: Could not get data for input var ', varName
187 0 : return
188 : endif
189 :
190 0 : if ( nlev .eq. 1 ) then
191 0 : outdata(1) = tmp(1)
192 0 : return ! no need to do interpolation
193 : endif
194 : ! if ( use_camiop .and. nlev.eq.plev) then
195 0 : if ( nlev.eq.plev .or. nlev.eq.plev+1) then
196 0 : outData(:nlev)= tmp(:nlev)! no need to do interpolation
197 : else
198 : !
199 : ! add the surface data if available, else
200 : ! fill in missing surface data by extrapolation
201 : !
202 0 : if(.not.scm_crm_mode) then
203 0 : if ( have_surfdat ) then
204 0 : tmp(npress) = surfdat
205 : else
206 0 : dy = press(npress-1) - press(npress-2)
207 0 : dx = tmp(npress-1) - tmp(npress-2)
208 0 : if ( dx .ne. 0.0_r8 ) then
209 0 : m = dy/dx
210 0 : tmp(npress) = ((press(npress) - press(npress-1)) / m ) + tmp(npress-1)
211 : else
212 0 : tmp(npress) = tmp(npress-1)
213 : endif
214 0 : surfdat = tmp(npress)
215 : endif
216 : endif
217 :
218 : #if DEBUG > 1
219 : !
220 : ! check data for missing values
221 : !
222 :
223 : STATUS = NF90_GET_ATT( NCID, varID, 'missing_value', missing_val )
224 : if ( STATUS .NE. NF90_NOERR ) then
225 : missing_val = -9999999.0_r8
226 : endif
227 : !
228 : ! reset status to zero
229 : !
230 : STATUS = 0
231 : !
232 : do i=1, npress
233 : if ( tmp(i) .eq. missing_val ) then
234 : write(iulog,*) 'ERROR - missing value found in ', varname
235 : write(iulog,*) 'time,lat,lon,lev = ' ,timeidx, latidx, lonidx, i
236 : stop
237 : endif
238 : enddo
239 : #endif
240 : !
241 0 : call interplevs( tmp(:npress), press, npress, ps, fill_ends, hyam, hybm, outdata )
242 :
243 : endif
244 :
245 0 : deallocate(tmp)
246 0 : return
247 0 : end subroutine getinterpncdata
248 :
249 0 : subroutine interplevs( inputdata, dplevs, nlev, &
250 0 : ps, fill_ends, hyam, hybm, outdata)
251 :
252 : use shr_kind_mod, only: r8 => shr_kind_r8, i8 => shr_kind_i8
253 : use interpolate_data, only: lininterp
254 : implicit none
255 :
256 : !
257 : ! WARNING: ps, siga and sigb must be initialized before calling this routine
258 : !
259 :
260 : !------------------------------Commons----------------------------------
261 :
262 : !-----------------------------------------------------------------------
263 :
264 :
265 : ! ------- inputs -----------
266 : integer, intent(in) :: nlev ! num press levels in dataset
267 :
268 : real(r8), intent(in) :: ps ! surface pressure
269 : real(r8), intent(in) :: hyam(:) ! a midpoint pressure
270 : real(r8), intent(in) :: hybm(:) ! b midpoint pressure
271 : real(r8), intent(in) :: inputdata(nlev) ! data from netcdf dataset
272 : real(r8), intent(in) :: dplevs(nlev) ! input data pressure levels
273 :
274 : logical, intent(in) :: fill_ends ! fill in missing end values(used for
275 : ! global model datasets)
276 :
277 :
278 : ! ------- outputs ----------
279 : real(r8), intent(inout) :: outdata(:) ! interpolated column data
280 :
281 : ! ------- locals -----------
282 :
283 : real(r8) mplevs( PLEV )
284 : real(r8) interpdata( PLEV )
285 :
286 :
287 : integer dstart_lev, dend_lev
288 : integer mstart_lev, mend_lev
289 : integer data_nlevs, model_nlevs, i
290 : integer STATUS
291 :
292 : !
293 : ! Initialize model_pressure_levels. ps should be set in the calling
294 : ! routine to the value in the dataset
295 : !
296 0 : do i = 1, plev
297 0 : mplevs( i ) = 1000.0_r8 * hyam( i ) + ps * hybm( i ) / 100.0_r8
298 : end do
299 : !
300 : ! the following algorithm assumes that pressures are increasing in the
301 : ! arrays
302 : !
303 : !
304 : ! Find the data pressure levels that are just outside the range
305 : ! of the model pressure levels, and that contain valid values
306 : !
307 : dstart_lev = 1
308 0 : do i= 1, nlev
309 0 : if ( dplevs(i) .LE. mplevs(1) ) dstart_lev = i
310 : end do
311 :
312 : dend_lev = nlev
313 0 : do i= nlev, 1, -1
314 0 : if ( dplevs(i) .GE. mplevs(plev) ) then
315 0 : dend_lev = i
316 : endif
317 : end do
318 : !
319 : ! Find the model pressure levels that are just inside the range
320 : ! of the data pressure levels
321 : !
322 : mstart_lev = 1
323 0 : do i=plev, 1, -1
324 0 : if ( mplevs( i ) .GE. dplevs( dstart_lev ) ) mstart_lev = i
325 : end do
326 :
327 : mend_lev = plev
328 0 : do i=1,plev
329 0 : if ( mplevs( i ) .LE. dplevs( dend_lev ) ) mend_lev = i
330 : end do
331 :
332 0 : data_nlevs = dend_lev - dstart_lev +1
333 0 : model_nlevs = mend_lev - mstart_lev +1
334 :
335 0 : call lininterp (inputdata(dstart_lev:dend_lev),dplevs(dstart_lev:dend_lev),data_nlevs, &
336 0 : interpdata,mplevs(mstart_lev:mend_lev),model_nlevs)
337 : !
338 : ! interpolate data onto the model pressure levels
339 : !
340 : !!$ call lininterp (inputdata,dplevs,nlev, &
341 : !!$ outdata(:plev),mplevs,plev)
342 0 : do i=1 , model_nlevs
343 0 : outdata( i+mstart_lev-1 ) = interpdata( i )
344 : end do
345 : !
346 : ! fill in the missing end values
347 : ! (usually done if this is global model dataset)
348 : !
349 0 : if ( fill_ends ) then
350 0 : do i=1, mstart_lev
351 0 : outdata(i) = inputdata(1)
352 : end do
353 0 : do i= mend_lev, plev
354 0 : outdata(i) = inputdata(nlev)
355 : end do
356 : end if
357 :
358 0 : return
359 : end subroutine interplevs
360 : end module getinterpnetcdfdata
|
