Line data Source code
1 : module hycoef
2 :
3 : use shr_kind_mod, only: r8 => shr_kind_r8
4 : use spmd_utils, only: masterproc
5 : use pmgrid, only: plev, plevp
6 : use cam_logfile, only: iulog
7 : use cam_abortutils, only: endrun
8 : use pio, only: file_desc_t, var_desc_t, &
9 : pio_inq_dimid, pio_inq_dimlen, pio_inq_varid, &
10 : pio_double, pio_def_dim, pio_def_var, &
11 : pio_put_var, pio_get_var, &
12 : pio_seterrorhandling, PIO_BCAST_ERROR, PIO_NOERR
13 :
14 : implicit none
15 : private
16 : save
17 :
18 : !-----------------------------------------------------------------------
19 : !
20 : ! Purpose: Hybrid level definitions: p = a*p0 + b*ps
21 : ! interfaces p(k) = hyai(k)*ps0 + hybi(k)*ps
22 : ! midpoints p(k) = hyam(k)*ps0 + hybm(k)*ps
23 : !
24 : ! Note: Module data with a target attribute are targets of pointers in hist_coord_t
25 : ! objects in the cam_history_support module. They are associated by the calls
26 : ! to add_hist_coord and add_vert_coord
27 : !
28 : !-----------------------------------------------------------------------
29 :
30 : real(r8), public, target :: hyai(plevp) ! ps0 component of hybrid coordinate - interfaces
31 : real(r8), public, target :: hyam(plev) ! ps0 component of hybrid coordinate - midpoints
32 : real(r8), public, target :: hybi(plevp) ! ps component of hybrid coordinate - interfaces
33 : real(r8), public, target :: hybm(plev) ! ps component of hybrid coordinate - midpoints
34 :
35 : real(r8), public :: etamid(plev) ! hybrid coordinate - midpoints
36 :
37 : real(r8), public :: hybd(plev) ! difference in b (hybi) across layers
38 : real(r8), public :: hypi(plevp) ! reference pressures at interfaces
39 : real(r8), public :: hypm(plev) ! reference pressures at midpoints
40 : real(r8), public :: hypd(plev) ! reference pressure layer thickness
41 : #ifdef planet_mars
42 : real(r8), public, protected :: ps0 = 6.0e1_r8 ! Base state surface pressure (pascals)
43 : real(r8), public, protected :: psr = 6.0e1_r8 ! Reference surface pressure (pascals)
44 : #else
45 : real(r8), public, protected :: ps0 = 1.0e5_r8 ! Base state surface pressure (pascals)
46 : real(r8), public, protected :: psr = 1.0e5_r8 ! Reference surface pressure (pascals)
47 : #endif
48 : real(r8), target :: alev(plev) ! level values (hPa) for 'lev' coord
49 : real(r8), target :: ailev(plevp) ! interface level values for 'ilev' coord
50 :
51 : integer, public :: nprlev ! number of pure pressure levels at top
52 :
53 : public hycoef_init
54 :
55 : type(var_desc_t) :: hyam_desc, hyai_desc, hybm_desc, hybi_desc, p0_desc
56 : public init_restart_hycoef, write_restart_hycoef
57 :
58 : !=======================================================================
59 : contains
60 : !=======================================================================
61 :
62 1536 : subroutine hycoef_init(file, psdry)
63 :
64 : use cam_history_support, only: add_hist_coord, add_vert_coord, formula_terms_t
65 :
66 : !-----------------------------------------------------------------------
67 : !
68 : ! Purpose:
69 : ! Defines the locations of model interfaces from input data in the
70 : ! hybrid coordinate scheme. Actual pressure values of model level
71 : ! interfaces are determined elsewhere from the fields set here.
72 : !
73 : ! Method:
74 : ! the following fields are set:
75 : ! hyai fraction of reference pressure used for interface pressures
76 : ! hyam fraction of reference pressure used for midpoint pressures
77 : ! hybi fraction of surface pressure used for interface pressures
78 : ! hybm fraction of surface pressure used for midpoint pressures
79 : ! hybd difference of hybi's
80 : ! hypi reference state interface pressures
81 : ! hypm reference state midpoint pressures
82 : ! hypd reference state layer thicknesses
83 : ! hypdln reference state layer thicknesses (log p)
84 : ! hyalph distance from interface to level (used in integrals)
85 : ! prsfac log pressure extrapolation factor (used to compute psl)
86 : !
87 : ! Author: B. Boville
88 : !
89 : !-----------------------------------------------------------------------
90 :
91 : ! arguments
92 : type(file_desc_t), intent(inout) :: file
93 : logical, optional, intent(in) :: psdry ! set true when coordinate is based
94 : ! on dry surface pressure
95 :
96 : ! local variables
97 : integer :: k ! Level index
98 : logical :: dry_coord
99 : real(r8) :: amean, bmean, atest, btest, eps
100 : type(formula_terms_t) :: formula_terms ! For the 'lev' and 'ilev' coords
101 : !-----------------------------------------------------------------------
102 :
103 : ! check for dry pressure coordinate (default is moist)
104 1536 : dry_coord = .false.
105 1536 : if (present(psdry)) dry_coord = psdry
106 :
107 : ! read hybrid coeficients
108 1536 : call hycoef_read(file)
109 :
110 : ! Set layer locations
111 1536 : nprlev = 0
112 144384 : do k=1,plev
113 :
114 : ! Interfaces. Set nprlev to the interface above, the first time a
115 : ! nonzero surface pressure contribution is found. "nprlev"
116 : ! identifies the lowest pure pressure interface.
117 :
118 144384 : if (nprlev==0 .and. hybi(k).ne.0.0_r8) nprlev = k - 1
119 : end do
120 :
121 : ! Set nprlev if no nonzero b's have been found. All interfaces are
122 : ! pure pressure. A pure pressure model requires other changes as well.
123 1536 : if (nprlev==0) nprlev = plev + 2
124 :
125 : ! Set delta sigma part of layer thickness and reference state midpoint
126 : ! pressures
127 144384 : do k=1,plev
128 142848 : hybd(k) = hybi(k+1) - hybi(k)
129 142848 : hypm(k) = hyam(k)*ps0 + hybm(k)*psr
130 144384 : etamid(k) = hyam(k) + hybm(k)
131 : end do
132 :
133 : ! Reference state interface pressures
134 145920 : do k=1,plevp
135 145920 : hypi(k) = hyai(k)*ps0 + hybi(k)*psr
136 : end do
137 :
138 : ! Reference state layer thicknesses
139 144384 : do k=1,plev
140 144384 : hypd(k) = hypi(k+1) - hypi(k)
141 : end do
142 :
143 : ! Test that A's and B's at full levels are arithmetic means of A's and
144 : ! B's at interfaces
145 1536 : eps = 1.e-05_r8
146 144384 : do k = 1,plev
147 142848 : amean = ( hyai(k+1) + hyai(k) )*0.5_r8
148 142848 : bmean = ( hybi(k+1) + hybi(k) )*0.5_r8
149 142848 : if(amean == 0._r8 .and. hyam(k) == 0._r8) then
150 1536 : atest = 0._r8
151 : else
152 141312 : atest = abs( amean - hyam(k) )/ ( 0.5_r8*( abs(amean + hyam(k)) ) )
153 : endif
154 142848 : if(bmean == 0._r8 .and. hybm(k) == 0._r8) then
155 66048 : btest = 0._r8
156 : else
157 76800 : btest = abs( bmean - hybm(k) )/ ( 0.5_r8*( abs(bmean + hybm(k)) ) )
158 : endif
159 142848 : if (atest > eps) then
160 0 : if (masterproc) then
161 0 : write(iulog,9850)
162 0 : write(iulog,*)'k,atest,eps=',k,atest,eps
163 : end if
164 : end if
165 :
166 144384 : if (btest > eps) then
167 0 : if (masterproc) then
168 0 : write(iulog,9850)
169 0 : write(iulog,*)'k,btest,eps=',k,btest,eps
170 : end if
171 : end if
172 : end do
173 :
174 : ! Add the information for the 'lev' and 'ilev' mdim history coordinates
175 : !
176 : ! The hybrid coordinate used by the SE dycore is based on a dry surface
177 : ! pressure. Hence it is the dry pressure rather than actual pressure
178 : ! that is computed by the formula_terms attribute. This coordinate is
179 : ! not described by the formula
180 : ! atmosphere_hybrid_sigma_pressure_coordinate since the formula
181 : ! associated with that name uses actual pressure values. Furthermore,
182 : ! the actual pressure field cannot be reconstructed from the hybrid
183 : ! coefficients and the surface pressure field. Hence in the case of a
184 : ! dry coordinate we add neither the standard_name nor the formula_terms
185 : ! attributes to the lev and ilev coordinates.
186 :
187 : ! 0.01 converts Pascals to millibars
188 144384 : alev(:plev) = 0.01_r8*ps0*(hyam(:plev) + hybm(:plev))
189 145920 : ailev(:plevp) = 0.01_r8*ps0*(hyai(:plevp) + hybi(:plevp))
190 :
191 1536 : if (dry_coord) then
192 : call add_vert_coord('lev', plev, &
193 : 'hybrid level at midpoints (1000*(A+B))', 'hPa', alev, &
194 1536 : positive='down')
195 : call add_hist_coord('hyam', plev, &
196 1536 : 'hybrid A coefficient at layer midpoints', '1', hyam, dimname='lev')
197 : call add_hist_coord('hybm', plev, &
198 1536 : 'hybrid B coefficient at layer midpoints', '1', hybm, dimname='lev')
199 : else
200 :
201 0 : formula_terms%a_name = 'hyam'
202 0 : formula_terms%a_long_name = 'hybrid A coefficient at layer midpoints'
203 0 : formula_terms%a_values => hyam
204 0 : formula_terms%b_name = 'hybm'
205 0 : formula_terms%b_long_name = 'hybrid B coefficient at layer midpoints'
206 0 : formula_terms%b_values => hybm
207 0 : formula_terms%p0_name = 'P0'
208 0 : formula_terms%p0_long_name = 'reference pressure'
209 0 : formula_terms%p0_units = 'Pa'
210 0 : formula_terms%p0_value = ps0
211 0 : formula_terms%ps_name = 'PS'
212 :
213 : call add_vert_coord('lev', plev, &
214 : 'hybrid level at midpoints (1000*(A+B))', 'hPa', alev, &
215 : positive='down', &
216 : standard_name='atmosphere_hybrid_sigma_pressure_coordinate', &
217 0 : formula_terms=formula_terms)
218 : end if
219 :
220 1536 : if (dry_coord) then
221 : call add_vert_coord('ilev', plevp, &
222 : 'hybrid level at interfaces (1000*(A+B))', 'hPa', ailev, &
223 1536 : positive='down')
224 : call add_hist_coord('hyai', plevp, &
225 1536 : 'hybrid A coefficient at layer interfaces', '1', hyai, dimname='ilev')
226 : call add_hist_coord('hybi', plevp, &
227 1536 : 'hybrid B coefficient at layer interfaces', '1', hybi, dimname='ilev')
228 : else
229 0 : formula_terms%a_name = 'hyai'
230 0 : formula_terms%a_long_name = 'hybrid A coefficient at layer interfaces'
231 0 : formula_terms%a_values => hyai
232 0 : formula_terms%b_name = 'hybi'
233 0 : formula_terms%b_long_name = 'hybrid B coefficient at layer interfaces'
234 0 : formula_terms%b_values => hybi
235 0 : formula_terms%p0_name = 'P0'
236 0 : formula_terms%p0_long_name = 'reference pressure'
237 0 : formula_terms%p0_units = 'Pa'
238 0 : formula_terms%p0_value = ps0
239 0 : formula_terms%ps_name = 'PS'
240 :
241 : call add_vert_coord('ilev', plevp, &
242 : 'hybrid level at interfaces (1000*(A+B))', 'hPa', ailev, &
243 : positive='down', &
244 : standard_name='atmosphere_hybrid_sigma_pressure_coordinate', &
245 0 : formula_terms=formula_terms)
246 : end if
247 :
248 1536 : if (masterproc) then
249 2 : write(iulog,'(a)')' Layer Locations (*1000) '
250 188 : do k=1,plev
251 186 : write(iulog,9800)k,hyai(k),hybi(k),hyai(k)+hybi(k)
252 188 : write(iulog,9810) hyam(k), hybm(k), hyam(k)+hybm(k)
253 : end do
254 :
255 2 : write(iulog,9800)plevp,hyai(plevp),hybi(plevp),hyai(plevp)+hybi(plevp)
256 2 : write(iulog,9820)
257 188 : do k=1,plev
258 186 : write(iulog,9830) k, hypi(k)
259 188 : write(iulog,9840) hypm(k), hypd(k)
260 : end do
261 2 : write(iulog,9830) plevp, hypi(plevp)
262 : end if
263 :
264 : 9800 format( 1x, i3, 3p, 3(f10.4,10x) )
265 : 9810 format( 1x, 3x, 3p, 3(10x,f10.4) )
266 : 9820 format(1x,'reference pressures (Pa)')
267 : 9830 format(1x,i3,f15.4)
268 : 9840 format(1x,3x,15x,2f15.4)
269 : 9850 format('HYCOEF: A and/or B vertical level coefficients at full',/, &
270 : ' levels are not the arithmetic mean of half-level values')
271 :
272 3072 : end subroutine hycoef_init
273 :
274 : !=======================================================================
275 :
276 1536 : subroutine init_restart_hycoef(File, vdimids)
277 :
278 : type(file_desc_t), intent(inout) :: File
279 : integer, intent(out) :: vdimids(:)
280 :
281 : ! PIO traps errors internally, no need to check ierr
282 :
283 : integer :: ierr
284 :
285 1536 : ierr = PIO_Def_Dim(File, 'lev', plev, vdimids(1))
286 1536 : ierr = PIO_Def_Dim(File, 'ilev', plevp, vdimids(2))
287 :
288 1536 : ierr = pio_def_var(File, 'hyai', pio_double, vdimids(2:2), hyai_desc)
289 1536 : ierr = pio_def_var(File, 'hyam', pio_double, vdimids(1:1), hyam_desc)
290 1536 : ierr = pio_def_var(File, 'hybi', pio_double, vdimids(2:2), hybi_desc)
291 1536 : ierr = pio_def_var(File, 'hybm', pio_double, vdimids(1:1), hybm_desc)
292 :
293 1536 : ierr = pio_def_var(File, 'P0', pio_double, p0_desc)
294 :
295 1536 : end subroutine init_restart_hycoef
296 :
297 : !=======================================================================
298 :
299 1536 : subroutine write_restart_hycoef(file)
300 :
301 : type(file_desc_t), intent(inout) :: File
302 :
303 : ! PIO traps errors internally, no need to check ierr
304 :
305 : integer :: ierr
306 :
307 1536 : ierr = pio_put_var(File, hyai_desc, hyai)
308 1536 : ierr = pio_put_var(File, hyam_desc, hyam)
309 1536 : ierr = pio_put_var(File, hybi_desc, hybi)
310 1536 : ierr = pio_put_var(File, hybm_desc, hybm)
311 :
312 1536 : ierr = pio_put_var(File, p0_desc, ps0)
313 :
314 1536 : end subroutine write_restart_hycoef
315 :
316 : !=======================================================================
317 :
318 1536 : subroutine hycoef_read(File)
319 :
320 : ! This code is used both for initial and restart reading.
321 :
322 : type(file_desc_t), intent(inout) :: File
323 :
324 : integer :: flev, filev, lev_dimid, ierr
325 : integer :: pio_errtype
326 :
327 : type(var_desc_t) :: p0_desc
328 :
329 : character(len=*), parameter :: routine = 'hycoef_read'
330 : !----------------------------------------------------------------------------
331 :
332 : ! PIO traps errors internally, no need to check ierr
333 :
334 1536 : ierr = PIO_Inq_DimID(File, 'lev', lev_dimid)
335 1536 : ierr = PIO_Inq_dimlen(File, lev_dimid, flev)
336 1536 : if (plev /= flev) then
337 0 : write(iulog,*) routine//': ERROR: file lev does not match model. lev (file, model):',flev, plev
338 0 : call endrun(routine//': ERROR: file lev does not match model.')
339 : end if
340 :
341 1536 : ierr = PIO_Inq_DimID(File, 'ilev', lev_dimid)
342 1536 : ierr = PIO_Inq_dimlen(File, lev_dimid, filev)
343 1536 : if (plevp /= filev) then
344 0 : write(iulog,*) routine//':ERROR: file ilev does not match model plevp (file, model):',filev, plevp
345 0 : call endrun(routine//':ERROR: file ilev does not match model.')
346 : end if
347 :
348 1536 : ierr = pio_inq_varid(File, 'hyai', hyai_desc)
349 1536 : ierr = pio_inq_varid(File, 'hyam', hyam_desc)
350 1536 : ierr = pio_inq_varid(File, 'hybi', hybi_desc)
351 1536 : ierr = pio_inq_varid(File, 'hybm', hybm_desc)
352 :
353 1536 : ierr = pio_get_var(File, hyai_desc, hyai)
354 1536 : ierr = pio_get_var(File, hybi_desc, hybi)
355 1536 : ierr = pio_get_var(File, hyam_desc, hyam)
356 1536 : ierr = pio_get_var(File, hybm_desc, hybm)
357 :
358 1536 : if (masterproc) then
359 2 : write(iulog,*) routine//': read hyai, hybi, hyam, hybm'
360 : end if
361 :
362 : ! Check whether file contains value for P0. If it does then use it
363 :
364 : ! Set PIO to return error codes.
365 1536 : call pio_seterrorhandling(file, PIO_BCAST_ERROR, pio_errtype)
366 :
367 1536 : ierr = pio_inq_varid(file, 'P0', p0_desc)
368 1536 : if (ierr == PIO_NOERR) then
369 768 : ierr = pio_get_var(file, p0_desc, ps0)
370 768 : if (ierr /= PIO_NOERR) then
371 0 : call endrun(routine//': reading P0.')
372 : end if
373 768 : psr = ps0
374 :
375 768 : if (masterproc) then
376 1 : write(iulog,*) routine//': read P0 value: ', ps0
377 : end if
378 :
379 : end if
380 :
381 : ! Put the error handling back the way it was
382 1536 : call pio_seterrorhandling(file, pio_errtype)
383 :
384 : #if ( defined OFFLINE_DYN )
385 : ! make sure top interface is non zero for fv dycore
386 : if (hyai(1) .eq. 0._r8) then
387 : if (hybm(1) .ne. 0.0_r8) then
388 : hyai(1) = hybm(1)*1.e-2_r8
389 : else if (hyam(1) .ne. 0.0_r8) then
390 : hyai(1) = hyam(1)*1.e-2_r8
391 : else
392 : call endrun('Not able to set hyai(1) to non-zero.')
393 : end if
394 : end if
395 : #endif
396 :
397 1536 : end subroutine hycoef_read
398 :
399 : !=======================================================================
400 :
401 : end module hycoef
|
