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