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