chgres_cube  1.12.0
 All Data Structures Files Functions Variables
thompson_mp_climo_data.F90
Go to the documentation of this file.
1 
4 
10 
11  use esmf
12  use netcdf
13  use program_setup, only : cycle_mon, cycle_day, cycle_hour, &
14  thomp_mp_climo_file
16 
17  implicit none
18 
19  private
20 
21  integer :: i_thomp_mp_climo
22 
23  integer :: j_thomp_mp_climo
24 
25  integer, public :: lev_thomp_mp_climo
26 
27 
28  type(esmf_grid) :: thomp_mp_climo_grid
29 
30 
31  type(esmf_field), public :: qnifa_climo_input_grid
32 
34  type(esmf_field), public :: qnwfa_climo_input_grid
35 
37  type(esmf_field), public :: thomp_pres_climo_input_grid
38 
40 
41  public :: read_thomp_mp_climo_data
43 
44  contains
45 
51 
52  implicit none
53 
54  integer :: error, ncid, rc, clb(2), cub(2)
55  integer :: i, j, localpet, npets, id_var
56  integer :: jda(8), jdow, jdoy, jday, id_dim
57  integer :: mm, mmm, mmp, mon1, mon2
58 
59  real(esmf_kind_r8), allocatable :: dummy3d(:,:,:)
60  real(esmf_kind_r8), allocatable :: dummy3d_mon1(:,:,:)
61  real(esmf_kind_r8), allocatable :: dummy3d_mon2(:,:,:)
62  real(esmf_kind_r8), pointer :: lat_ptr(:,:), lon_ptr(:,:)
63  real(esmf_kind_r8), allocatable :: lons(:), lats(:)
64  real :: rjday, dayhf(13), wei1m, wei2m
65 
66  type(esmf_vm) :: vm
67 
68  type(esmf_polekind_flag) :: polekindflag(2)
69 
70  data dayhf/ 15.5, 45.0, 74.5,105.0,135.5,166.0, &
71  196.5,227.5,258.0,288.5,319.0,349.5,380.5/
72 
73 !-----------------------------------------------------------------------------------
74 ! Open the file and read the grid dimensions and latitude/longitude.
75 !-----------------------------------------------------------------------------------
76 
77  print*,"- READ THOMP_MP_CLIMO_FILE: ", trim(thomp_mp_climo_file)
78  error=nf90_open(trim(thomp_mp_climo_file),nf90_nowrite,ncid)
79  call netcdf_err(error, 'opening: '//trim(thomp_mp_climo_file) )
80 
81  error=nf90_inq_dimid(ncid, 'lat', id_dim)
82  call netcdf_err(error, 'reading lat id')
83  error=nf90_inquire_dimension(ncid,id_dim,len=j_thomp_mp_climo)
84  call netcdf_err(error, 'reading lat')
85 
86  error=nf90_inq_dimid(ncid, 'lon', id_dim)
87  call netcdf_err(error, 'reading lon id')
88  error=nf90_inquire_dimension(ncid,id_dim,len=i_thomp_mp_climo)
89  call netcdf_err(error, 'reading lon')
90 
91  error=nf90_inq_dimid(ncid, 'plev', id_dim)
92  call netcdf_err(error, 'reading plev id')
93  error=nf90_inquire_dimension(ncid,id_dim,len=lev_thomp_mp_climo)
94  call netcdf_err(error, 'reading plev')
95 
96  allocate(lons(i_thomp_mp_climo))
97  allocate(lats(j_thomp_mp_climo))
98  error=nf90_inq_varid(ncid, 'lon', id_var)
99  call netcdf_err(error, 'reading lon field id' )
100  error=nf90_get_var(ncid, id_var, lons)
101  call netcdf_err(error, 'reading grid longitude' )
102  error=nf90_inq_varid(ncid, 'lat', id_var)
103  call netcdf_err(error, 'reading lat field id' )
104  error=nf90_get_var(ncid, id_var, lats)
105  call netcdf_err(error, 'reading grid latitude' )
106 
107 !-----------------------------------------------------------------------------------
108 ! Now that we have the grid information, create the esmf grid object.
109 !-----------------------------------------------------------------------------------
110 
111  print*,"- CALL VMGetGlobal"
112  call esmf_vmgetglobal(vm, rc=rc)
113  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
114  call error_handler("IN VMGetGlobal", rc)
115 
116  print*,"- CALL VMGet"
117  call esmf_vmget(vm, localpet=localpet, petcount=npets, rc=rc)
118  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
119  call error_handler("IN VMGet", rc)
120 
121  polekindflag(1:2) = esmf_polekind_monopole
122 
123  print*,"- CALL GridCreate1PeriDim FOR THOMP MP CLIMO GRID."
124  thomp_mp_climo_grid = esmf_gridcreate1peridim(minindex=(/1,1/), &
125  maxindex=(/i_thomp_mp_climo,j_thomp_mp_climo/), &
126  polekindflag=polekindflag, &
127  periodicdim=1, &
128  poledim=2, &
129  coordsys=esmf_coordsys_sph_deg, &
130  regdecomp=(/1,npets/), &
131  indexflag=esmf_index_global, rc=rc)
132  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
133  call error_handler("IN GridCreate1PeriDim", rc)
134 
135  print*,"- CALL GridAddCoord FOR THOMP MP CLIMO GRID."
136  call esmf_gridaddcoord(thomp_mp_climo_grid, &
137  staggerloc=esmf_staggerloc_center, rc=rc)
138  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
139  call error_handler("IN GridAddCoord", rc)
140 
141 !-----------------------------------------------------------------------------------
142 ! Set the grid object lat/lon.
143 !-----------------------------------------------------------------------------------
144 
145  print*,"- CALL GridGetCoord FOR INPUT GRID X-COORD."
146  nullify(lon_ptr)
147  call esmf_gridgetcoord(thomp_mp_climo_grid, &
148  staggerloc=esmf_staggerloc_center, &
149  coorddim=1, &
150  farrayptr=lon_ptr, rc=rc)
151  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
152  call error_handler("IN GridGetCoord", rc)
153 
154  print*,"- CALL GridGetCoord FOR INPUT GRID Y-COORD."
155  nullify(lat_ptr)
156  call esmf_gridgetcoord(thomp_mp_climo_grid, &
157  staggerloc=esmf_staggerloc_center, &
158  coorddim=2, &
159  computationallbound=clb, &
160  computationalubound=cub, &
161  farrayptr=lat_ptr, rc=rc)
162  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
163  call error_handler("IN GridGetCoord", rc)
164 
165  do i = clb(1), cub(1)
166  lon_ptr(i,:) = lons(i)
167  enddo
168 
169  do j = clb(2), cub(2)
170  lat_ptr(:,j) = lats(j)
171  enddo
172 
173 !-----------------------------------------------------------------------------------
174 ! Create esmf fields for the two tracers and 3-d pressure.
175 !-----------------------------------------------------------------------------------
176 
177  print*,"- CALL FieldCreate FOR QNIFA INPUT CLIMO."
178  qnifa_climo_input_grid = esmf_fieldcreate(thomp_mp_climo_grid, &
179  typekind=esmf_typekind_r8, &
180  staggerloc=esmf_staggerloc_center, &
181  ungriddedlbound=(/1/), &
182  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
183  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
184  call error_handler("IN FieldCreate", rc)
185 
186  print*,"- CALL FieldCreate FOR QNWFA INPUT CLIMO."
187  qnwfa_climo_input_grid = esmf_fieldcreate(thomp_mp_climo_grid, &
188  typekind=esmf_typekind_r8, &
189  staggerloc=esmf_staggerloc_center, &
190  ungriddedlbound=(/1/), &
191  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
192  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
193  call error_handler("IN FieldCreate", rc)
194 
195  print*,"- CALL FieldCreate FOR THOMP PRESS CLIMO."
196  thomp_pres_climo_input_grid = esmf_fieldcreate(thomp_mp_climo_grid, &
197  typekind=esmf_typekind_r8, &
198  staggerloc=esmf_staggerloc_center, &
199  ungriddedlbound=(/1/), &
200  ungriddedubound=(/lev_thomp_mp_climo/), rc=rc)
201  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
202  call error_handler("IN FieldCreate", rc)
203 
204 !-----------------------------------------------------------------------------------
205 ! Data are monthly and valid at the 15th of the month. Compute time interpolation
206 ! weights for the current cycle.
207 !-----------------------------------------------------------------------------------
208 
209  jda=0
210  jda(1) = 2007
211  if (cycle_mon == 2 .and. cycle_day == 29) then ! leap year
212  jda(2) = 3
213  jda(3) = 1
214  else
215  jda(2) = cycle_mon
216  jda(3) = cycle_day
217  endif
218 
219  jda(5) = cycle_hour
220  jdow = 0
221  jdoy = 0
222  jday = 0
223  call w3doxdat(jda,jdow,jdoy,jday)
224  rjday = float(jdoy) + float(jda(5)) / 24.
225  if(rjday < dayhf(1)) rjday = rjday + 365.
226 
227  do mm=1,12
228  mmm = mm
229  mmp = mm + 1
230  if(rjday >= dayhf(mmm) .and. rjday < dayhf(mmp)) then
231  mon1 = mmm
232  mon2 = mmp
233  exit
234  endif
235  enddo
236 
237  wei1m = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
238  wei2m = 1.0 - wei1m
239 
240  if (mon2==13) mon2=1
241 
242  print*,"- BOUNDING MONTHS AND INTERPOLATION WEIGHTS: ", mon1, wei1m, mon2, wei2m
243 
244 !-----------------------------------------------------------------------------------
245 ! Read tracers and 3-d pressure for each bounding month. Then linearly
246 ! interpolate in time.
247 !-----------------------------------------------------------------------------------
248 
249  if (localpet == 0) then
250  allocate(dummy3d(i_thomp_mp_climo, j_thomp_mp_climo, lev_thomp_mp_climo))
251  dummy3d = 0.0
252  allocate(dummy3d_mon1(i_thomp_mp_climo, j_thomp_mp_climo, lev_thomp_mp_climo))
253  dummy3d_mon1 = 0.0
254  allocate(dummy3d_mon2(i_thomp_mp_climo, j_thomp_mp_climo, lev_thomp_mp_climo))
255  dummy3d_mon2 = 0.0
256  else
257  allocate(dummy3d(0,0,0))
258  allocate(dummy3d_mon1(0,0,0))
259  allocate(dummy3d_mon2(0,0,0))
260  endif
261 
262  if (localpet == 0) then
263  print*,"- READ QNIFA FOR BOUNDING MONTH 1"
264  error=nf90_inq_varid(ncid, 'nifa', id_var)
265  call netcdf_err(error, 'reading nifa field id' )
266  error=nf90_get_var(ncid, id_var, dummy3d_mon1, start=(/1,1,1,mon1/), &
267  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
268  call netcdf_err(error, 'reading nifa month1 field' )
269  print*,"- READ QNIFA FOR BOUNDING MONTH 2"
270  error=nf90_get_var(ncid, id_var, dummy3d_mon2, start=(/1,1,1,mon2/), &
271  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
272  call netcdf_err(error, 'reading nifa month2 field' )
273  dummy3d(:,:,:) = wei1m * dummy3d_mon1 + wei2m * dummy3d_mon2
274  endif
275 
276  print*,"- CALL FieldScatter FOR qnifa input climo."
277  call esmf_fieldscatter(qnifa_climo_input_grid, dummy3d, rootpet=0, rc=rc)
278  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
279  call error_handler("IN FieldScatter", rc)
280 
281  if (localpet == 0) then
282  print*,"- READ QNWFA FOR BOUNDING MONTH 1"
283  error=nf90_inq_varid(ncid, 'nwfa', id_var)
284  call netcdf_err(error, 'reading nwfa field id' )
285  error=nf90_get_var(ncid, id_var, dummy3d_mon1, start=(/1,1,1,mon1/), &
286  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
287  call netcdf_err(error, 'reading nwfa month1 field' )
288  print*,"- READ QNWFA FOR BOUNDING MONTH 2"
289  error=nf90_get_var(ncid, id_var, dummy3d_mon2, start=(/1,1,1,mon2/), &
290  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
291  call netcdf_err(error, 'reading nwfa month2 field' )
292  dummy3d(:,:,:) = wei1m * dummy3d_mon1 + wei2m * dummy3d_mon2
293  endif
294 
295  print*,"- CALL FieldScatter FOR qnwfa input climo."
296  call esmf_fieldscatter(qnwfa_climo_input_grid, dummy3d, rootpet=0, rc=rc)
297  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
298  call error_handler("IN FieldScatter", rc)
299 
300  if (localpet == 0) then
301  print*,"- READ PRESSURE FOR BOUNDING MONTH 1"
302  error=nf90_inq_varid(ncid, 'prs', id_var)
303  call netcdf_err(error, 'reading prs field id' )
304  error=nf90_get_var(ncid, id_var, dummy3d_mon1, start=(/1,1,1,mon1/), &
305  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
306  call netcdf_err(error, 'reading prs month1 field' )
307  print*,"- READ PRESSURE FOR BOUNDING MONTH 2"
308  error=nf90_get_var(ncid, id_var, dummy3d_mon2, start=(/1,1,1,mon2/), &
309  count=(/i_thomp_mp_climo,j_thomp_mp_climo,lev_thomp_mp_climo,1/) )
310  call netcdf_err(error, 'reading prs month2 field' )
311  dummy3d(:,:,:) = wei1m * dummy3d_mon1 + wei2m * dummy3d_mon2
312  endif
313 
314  print*,"- CALL FieldScatter FOR thomp press."
315  call esmf_fieldscatter(thomp_pres_climo_input_grid, dummy3d, rootpet=0, rc=rc)
316  if(esmf_logfounderror(rctocheck=rc,msg=esmf_logerr_passthru,line=__line__,file=__file__)) &
317  call error_handler("IN FieldScatter", rc)
318 
319  error=nf90_close(ncid)
320 
321  deallocate(lons, lats, dummy3d, dummy3d_mon1, dummy3d_mon2)
322 
323  end subroutine read_thomp_mp_climo_data
324 
329 
330  implicit none
331 
332  integer :: rc
333 
334  call esmf_griddestroy(thomp_mp_climo_grid, rc=rc)
335  call esmf_fielddestroy(thomp_pres_climo_input_grid, rc=rc)
336  call esmf_fielddestroy(qnifa_climo_input_grid, rc=rc)
337  call esmf_fielddestroy(qnwfa_climo_input_grid, rc=rc)
338 
339  end subroutine cleanup_thomp_mp_climo_input_data
340 
341  end module thompson_mp_climo_data
subroutine netcdf_err(err, string)
Error handler for netcdf.
Definition: utils.F90:34
subroutine, public cleanup_thomp_mp_climo_input_data
Free up memory associated with this module.
Module to read the Thompson climatological MP data file and set up the associated esmf field and grid...
subroutine error_handler(string, rc)
General error handler.
Definition: utils.F90:12
subroutine, public read_thomp_mp_climo_data
Read Thompson climatological MP data file and time interpolate data to current cycle time...
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...