chgres_cube  1.13.0
 All Data Structures Files Functions Variables
thompson_mp_climo_data.F90
Go to the documentation of this file.
1 
4 
9  module thompson_mp_climo_data
10 
11  use esmf
12  use netcdf
13  use program_setup, only : cycle_mon, cycle_day, cycle_hour, &
14  thomp_mp_climo_file
15  use utilities, only : error_handler, netcdf_err
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
42  public :: cleanup_thomp_mp_climo_input_data
43 
44  contains
45 
50  subroutine read_thomp_mp_climo_data
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 
328  subroutine cleanup_thomp_mp_climo_input_data
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
utilities::netcdf_err
subroutine netcdf_err(err, string)
Error handler for netcdf.
Definition: utils.F90:34
thompson_mp_climo_data::cleanup_thomp_mp_climo_input_data
subroutine, public cleanup_thomp_mp_climo_input_data
Free up memory associated with this module.
Definition: thompson_mp_climo_data.F90:328
thompson_mp_climo_data
Module to read the Thompson climatological MP data file and set up the associated esmf field and grid...
Definition: thompson_mp_climo_data.F90:9
utilities::error_handler
subroutine error_handler(string, rc)
General error handler.
Definition: utils.F90:12
thompson_mp_climo_data::read_thomp_mp_climo_data
subroutine, public read_thomp_mp_climo_data
Read Thompson climatological MP data file and time interpolate data to current cycle time...
Definition: thompson_mp_climo_data.F90:50
program_setup
This module contains code to read the setup namelist file, handle the varmap file for GRIB2 data...
Definition: program_setup.F90:9