Line data Source code
1 : !----------------------------------------------------------------------------------
2 : ! circulation diagnostics -- terms of the Transformed Eulerian Mean (TEM) equation
3 : !
4 : !----------------------------------------------------------------------------------
5 : module phys_grid_ctem
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use ppgrid, only: begchunk, endchunk, pcols, pver
8 : use physics_types, only: physics_state
9 : use cam_history, only: addfld, outfld
10 : use zonal_mean_mod,only: ZonalAverage_t, ZonalMean_t
11 : use physconst, only: pi
12 : use cam_logfile, only: iulog
13 : use cam_abortutils,only: endrun, handle_allocate_error
14 : use namelist_utils,only: find_group_name
15 : use spmd_utils, only: masterproc, mpi_integer, masterprocid, mpicom
16 : use time_manager, only: get_step_size, get_nstep
17 :
18 : use shr_const_mod, only: rgas => shr_const_rgas ! J/K/kmole
19 : use shr_const_mod, only: grav => shr_const_g ! m/s2
20 : use air_composition, only: mbarv ! g/mole
21 : use string_utils, only: int2str
22 :
23 : implicit none
24 :
25 : private
26 : public :: phys_grid_ctem_readnl
27 : public :: phys_grid_ctem_reg
28 : public :: phys_grid_ctem_init
29 : public :: phys_grid_ctem_diags
30 : public :: phys_grid_ctem_final
31 :
32 : type(ZonalMean_t) :: ZMobj
33 : type(ZonalAverage_t) :: ZAobj
34 :
35 : integer :: nzalat = -huge(1)
36 : integer :: nzmbas = -huge(1)
37 :
38 : integer :: ntimesteps = -huge(1) ! number of time steps bewteen TEM calculations
39 :
40 : logical :: do_tem_diags = .false.
41 :
42 : contains
43 :
44 : !------------------------------------------------------------------------------
45 : !------------------------------------------------------------------------------
46 1536 : subroutine phys_grid_ctem_readnl(nlfile)
47 : character(len=*), intent(in) :: nlfile
48 : integer :: ierr, unitn
49 :
50 : character(len=*), parameter :: prefix = 'phys_grid_ctem_readnl: '
51 : real(r8) :: dtime
52 :
53 : integer :: phys_grid_ctem_zm_nbas
54 : integer :: phys_grid_ctem_za_nlat
55 : integer :: phys_grid_ctem_nfreq
56 :
57 : namelist /phys_grid_ctem_opts/ phys_grid_ctem_zm_nbas, phys_grid_ctem_za_nlat, phys_grid_ctem_nfreq
58 :
59 1536 : phys_grid_ctem_zm_nbas = 0
60 1536 : phys_grid_ctem_za_nlat = 0
61 1536 : phys_grid_ctem_nfreq = 0
62 :
63 : ! Read in namelist values
64 : !------------------------
65 1536 : if(masterproc) then
66 2 : open(newunit=unitn, file=trim(nlfile), status='old')
67 2 : call find_group_name(unitn, 'phys_grid_ctem_opts', status=ierr)
68 2 : if(ierr == 0) then
69 0 : read(unitn,phys_grid_ctem_opts,iostat=ierr)
70 0 : if(ierr /= 0) then
71 0 : call endrun(prefix//'ERROR reading namelist')
72 : end if
73 : end if
74 2 : close(unitn)
75 : end if
76 :
77 1536 : call MPI_bcast(phys_grid_ctem_zm_nbas, 1, mpi_integer, masterprocid, mpicom, ierr)
78 1536 : if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_zm_nbas')
79 1536 : call MPI_bcast(phys_grid_ctem_za_nlat, 1, mpi_integer, masterprocid, mpicom, ierr)
80 1536 : if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_za_nlat')
81 1536 : call MPI_bcast(phys_grid_ctem_nfreq, 1, mpi_integer, masterprocid, mpicom, ierr)
82 1536 : if (ierr /= 0) call endrun(prefix//'FATAL: mpi_bcast: phys_grid_ctem_nfreq')
83 :
84 1536 : do_tem_diags = .false.
85 1536 : if (phys_grid_ctem_nfreq/=0) then
86 0 : if (.not.(phys_grid_ctem_zm_nbas>0 .and. phys_grid_ctem_za_nlat>0)) then
87 : call endrun(prefix//'inconsistent phys_grid_ctem namelist settings -- phys_grid_ctem_zm_nbas=' &
88 0 : //int2str(phys_grid_ctem_zm_nbas)//', phys_grid_ctem_za_nlat='//int2str(phys_grid_ctem_za_nlat))
89 : end if
90 0 : if (phys_grid_ctem_nfreq>0) then
91 0 : ntimesteps = phys_grid_ctem_nfreq
92 : else
93 0 : dtime = get_step_size()
94 0 : ntimesteps = nint( -phys_grid_ctem_nfreq*3600._r8/dtime )
95 : end if
96 0 : if (ntimesteps<1) then
97 : call endrun(prefix//'invalid ntimesteps -- phys_grid_ctem_nfreq needs to be a larger negative value ' &
98 0 : //'or the model time step needs to be shorter')
99 : end if
100 0 : do_tem_diags = .true.
101 : end if
102 :
103 1536 : if (masterproc) then
104 2 : if (do_tem_diags) then
105 0 : write(iulog,*) 'TEM diagnostics will be calculated every ',ntimesteps,' time steps'
106 0 : write(iulog,*) ' phys_grid_ctem_zm_nbas = ', phys_grid_ctem_zm_nbas
107 0 : write(iulog,*) ' phys_grid_ctem_za_nlat = ', phys_grid_ctem_za_nlat
108 0 : write(iulog,*) ' phys_grid_ctem_nfreq = ', phys_grid_ctem_nfreq
109 : else
110 2 : write(iulog,*) 'TEM diagnostics will not be performed'
111 : end if
112 : endif
113 :
114 1536 : if (do_tem_diags) then
115 0 : nzalat = phys_grid_ctem_za_nlat
116 0 : nzmbas = phys_grid_ctem_zm_nbas
117 : end if
118 :
119 1536 : end subroutine phys_grid_ctem_readnl
120 :
121 : !-----------------------------------------------------------------------------
122 : !-----------------------------------------------------------------------------
123 1536 : subroutine phys_grid_ctem_reg
124 :
125 : use cam_grid_support, only: horiz_coord_t, horiz_coord_create, iMap, cam_grid_register
126 :
127 : type(horiz_coord_t), pointer :: zalon_coord
128 : type(horiz_coord_t), pointer :: zalat_coord
129 1536 : integer(iMap), pointer :: grid_map(:,:)
130 :
131 3072 : real(r8) :: zalats(nzalat)
132 3072 : real(r8) :: area(nzalat)
133 : real(r8) :: zalons(1)
134 : real(r8) :: dlatrad, dlatdeg, lat1, lat2
135 : real(r8) :: total_area
136 : real(r8) :: total_wght
137 : integer :: j, astat
138 :
139 : real(r8), parameter :: latdeg0 = -90._r8
140 : real(r8), parameter :: latrad0 = -pi*0.5_r8
141 : real(r8), parameter :: fourpi = pi*4._r8
142 :
143 : integer, parameter :: ctem_zavg_phys_decomp = 333 ! Must be unique within CAM
144 :
145 1536 : if (.not.do_tem_diags) return
146 :
147 0 : nullify(zalat_coord)
148 0 : nullify(zalon_coord)
149 0 : nullify(grid_map)
150 :
151 0 : zalons(1) = 0._r8
152 :
153 0 : dlatrad = pi/real(nzalat,kind=r8)
154 0 : dlatdeg = 180._r8/real(nzalat,kind=r8)
155 0 : total_area = 0._r8
156 0 : total_wght = 0._r8
157 :
158 : ! calculate latitudes and areas of zonal average grid boxes
159 0 : do j = 1,nzalat
160 0 : zalats(j) = latdeg0 + (real(j,kind=r8)-0.5_r8)*dlatdeg
161 0 : lat1 = latrad0 + real(j-1,kind=r8)*dlatrad
162 0 : lat2 = latrad0 + real(j ,kind=r8)*dlatrad
163 0 : area(j) = 2._r8*pi*(sin(lat2)-sin(lat1))
164 0 : total_area = total_area + area(j)
165 0 : total_wght = total_wght + 0.5_r8*(sin(lat2)-sin(lat1))
166 : end do
167 :
168 : ! sanity check
169 0 : if ( abs(1._r8-total_wght)>1.e-12_r8 .or. abs(fourpi-total_area)>1.e-12_r8 ) then
170 0 : call endrun('phys_grid_ctem_reg: problem with area/wght calc')
171 : end if
172 :
173 : ! initialize zonal-average and zonal-mean utility objects
174 0 : call ZAobj%init(zalats,area,nzalat,GEN_GAUSSLATS=.false.)
175 0 : call ZMobj%init(nzmbas)
176 :
177 : ! Zonal average grid for history fields
178 :
179 0 : zalat_coord => horiz_coord_create('zalat', '', nzalat, 'latitude', 'degrees_north', 1, nzalat, zalats)
180 0 : zalon_coord => horiz_coord_create('zalon', '', 1, 'longitude', 'degrees_east', 1, 1, zalons)
181 :
182 : ! grid decomposition map
183 0 : allocate(grid_map(4,nzalat), stat=astat)
184 0 : call handle_allocate_error(astat, 'phys_grid_ctem_reg', 'grid_map')
185 :
186 0 : do j = 1,nzalat
187 0 : grid_map(1,j) = 1
188 0 : grid_map(2,j) = j
189 0 : if (masterproc) then
190 0 : grid_map(3,j) = 1
191 0 : grid_map(4,j) = j
192 : else
193 0 : grid_map(3,j) = 0
194 0 : grid_map(4,j) = 0
195 : end if
196 : end do
197 :
198 : ! register the zonal average grid
199 : call cam_grid_register('ctem_zavg_phys', ctem_zavg_phys_decomp, zalat_coord, zalon_coord, grid_map, &
200 0 : unstruct=.false., zonal_grid=.true.)
201 :
202 1536 : end subroutine phys_grid_ctem_reg
203 :
204 : !-----------------------------------------------------------------------------
205 : !-----------------------------------------------------------------------------
206 1536 : subroutine phys_grid_ctem_init
207 :
208 1536 : if (.not.do_tem_diags) return
209 :
210 0 : call addfld ('Uzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean zonal wind', gridname='ctem_zavg_phys' )
211 0 : call addfld ('Vzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean meridional wind', gridname='ctem_zavg_phys' )
212 0 : call addfld ('Wzm', (/'lev'/), 'A','m s-1', 'Zonal-Mean vertical wind', gridname='ctem_zavg_phys' )
213 0 : call addfld ('THzm', (/'lev'/), 'A','K', 'Zonal-Mean potential temp', gridname='ctem_zavg_phys' )
214 0 : call addfld ('VTHzm',(/'lev'/), 'A','K m s-1','Meridional Heat Flux:', gridname='ctem_zavg_phys')
215 0 : call addfld ('WTHzm',(/'lev'/), 'A','K m s-1','Vertical Heat Flux:', gridname='ctem_zavg_phys')
216 0 : call addfld ('UVzm', (/'lev'/), 'A','m2 s-2', 'Meridional Flux of Zonal Momentum', gridname='ctem_zavg_phys')
217 0 : call addfld ('UWzm', (/'lev'/), 'A','m2 s-2', 'Vertical Flux of Zonal Momentum', gridname='ctem_zavg_phys')
218 0 : call addfld ('THphys',(/'lev'/), 'A', 'K', 'Potential temp', gridname='physgrid' )
219 :
220 1536 : end subroutine phys_grid_ctem_init
221 :
222 : !-----------------------------------------------------------------------------
223 : !-----------------------------------------------------------------------------
224 370944 : subroutine phys_grid_ctem_diags(phys_state)
225 : type(physics_state), intent(in) :: phys_state(begchunk:endchunk)
226 :
227 : character(len=*), parameter :: prefix = 'phys_grid_ctem_diags: '
228 :
229 741888 : real(r8) :: u(pcols,pver,begchunk:endchunk)
230 741888 : real(r8) :: v(pcols,pver,begchunk:endchunk)
231 741888 : real(r8) :: w(pcols,pver,begchunk:endchunk)
232 :
233 741888 : real(r8) :: uzm(pcols,pver,begchunk:endchunk)
234 741888 : real(r8) :: vzm(pcols,pver,begchunk:endchunk)
235 741888 : real(r8) :: wzm(pcols,pver,begchunk:endchunk)
236 :
237 741888 : real(r8) :: ud(pcols,pver,begchunk:endchunk)
238 741888 : real(r8) :: vd(pcols,pver,begchunk:endchunk)
239 741888 : real(r8) :: wd(pcols,pver,begchunk:endchunk)
240 741888 : real(r8) :: thd(pcols,pver,begchunk:endchunk)
241 :
242 741888 : real(r8) :: uvp(pcols,pver,begchunk:endchunk)
243 741888 : real(r8) :: uwp(pcols,pver,begchunk:endchunk)
244 741888 : real(r8) :: vthp(pcols,pver,begchunk:endchunk)
245 741888 : real(r8) :: wthp(pcols,pver,begchunk:endchunk)
246 :
247 : integer :: lchnk, ncol, j, k
248 :
249 : ! potential temperature
250 741888 : real(r8) :: theta(pcols,pver,begchunk:endchunk)
251 741888 : real(r8) :: thzm(pcols,pver,begchunk:endchunk)
252 :
253 741888 : real(r8) :: uvza(nzalat,pver)
254 741888 : real(r8) :: uwza(nzalat,pver)
255 741888 : real(r8) :: vthza(nzalat,pver)
256 741888 : real(r8) :: wthza(nzalat,pver)
257 :
258 741888 : real(r8) :: uza(nzalat,pver)
259 741888 : real(r8) :: vza(nzalat,pver)
260 741888 : real(r8) :: wza(nzalat,pver)
261 370944 : real(r8) :: thza(nzalat,pver)
262 :
263 : real(r8) :: sheight(pcols,pver) ! pressure scale height (m)
264 :
265 370944 : if (.not.do_calc()) return
266 :
267 0 : do lchnk = begchunk,endchunk
268 :
269 0 : ncol = phys_state(lchnk)%ncol
270 :
271 : ! scale height
272 0 : sheight(:ncol,:) = phys_state(lchnk)%t(:ncol,:) * rgas / ( mbarv(:ncol,:,lchnk) * grav ) ! meters
273 :
274 : ! potential temperature
275 0 : theta(:ncol,:,lchnk) = phys_state(lchnk)%t(:ncol,:) * phys_state(lchnk)%exner(:ncol,:)
276 :
277 : ! vertical velocity
278 0 : w(:ncol,:,lchnk) = -sheight(:ncol,:) * phys_state(lchnk)%omega(:ncol,:) / phys_state(lchnk)%pmid(:ncol,:)
279 :
280 0 : u(:ncol,:,lchnk) = phys_state(lchnk)%u(:ncol,:)
281 0 : v(:ncol,:,lchnk) = phys_state(lchnk)%v(:ncol,:)
282 :
283 : end do
284 :
285 : ! zonal means evaluated on the physics grid (3D) to be used in the deviations calculation below
286 0 : uzm(:,:,:) = zmean_fld(u(:,:,:))
287 0 : vzm(:,:,:) = zmean_fld(v(:,:,:))
288 0 : wzm(:,:,:) = zmean_fld(w(:,:,:))
289 0 : thzm(:,:,:) = zmean_fld(theta(:,:,:))
290 :
291 : ! diagnostic output
292 0 : do lchnk = begchunk, endchunk
293 0 : call outfld( 'THphys', theta(:,:,lchnk), pcols, lchnk)
294 : end do
295 :
296 0 : do lchnk = begchunk,endchunk
297 0 : ncol = phys_state(lchnk)%ncol
298 0 : do k = 1,pver
299 : ! zonal deviations
300 0 : thd(:ncol,k,lchnk) = theta(:ncol,k,lchnk) - thzm(:ncol,k,lchnk)
301 0 : ud(:ncol,k,lchnk) = u(:ncol,k,lchnk) - uzm(:ncol,k,lchnk)
302 0 : vd(:ncol,k,lchnk) = v(:ncol,k,lchnk) - vzm(:ncol,k,lchnk)
303 0 : wd(:ncol,k,lchnk) = w(:ncol,k,lchnk) - wzm(:ncol,k,lchnk)
304 : ! fluxes
305 0 : uvp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * vd(:ncol,k,lchnk)
306 0 : uwp(:ncol,k,lchnk) = ud(:ncol,k,lchnk) * wd(:ncol,k,lchnk)
307 0 : vthp(:ncol,k,lchnk) = vd(:ncol,k,lchnk) * thd(:ncol,k,lchnk)
308 0 : wthp(:ncol,k,lchnk) = wd(:ncol,k,lchnk) * thd(:ncol,k,lchnk)
309 : end do
310 : end do
311 :
312 : ! evaluate and output fluxes on the zonal-average grid
313 0 : call ZAobj%binAvg(uvp, uvza)
314 0 : call ZAobj%binAvg(uwp, uwza)
315 0 : call ZAobj%binAvg(vthp, vthza)
316 0 : call ZAobj%binAvg(wthp, wthza)
317 :
318 0 : if (any(abs(uvza)>1.e20_r8)) call endrun(prefix//'bad values in uvza')
319 0 : if (any(abs(uwza)>1.e20_r8)) call endrun(prefix//'bad values in uwza')
320 0 : if (any(abs(vthza)>1.e20_r8)) call endrun(prefix//'bad values in vthza')
321 0 : if (any(abs(wthza)>1.e20_r8)) call endrun(prefix//'bad values in wthza')
322 :
323 0 : call ZAobj%binAvg(uzm, uza)
324 0 : call ZAobj%binAvg(vzm, vza)
325 0 : call ZAobj%binAvg(wzm, wza)
326 0 : call ZAobj%binAvg(thzm, thza)
327 :
328 0 : if (any(abs(uza)>1.e20_r8)) call endrun(prefix//'bad values in uza')
329 0 : if (any(abs(vza)>1.e20_r8)) call endrun(prefix//'bad values in vza')
330 0 : if (any(abs(wza)>1.e20_r8)) call endrun(prefix//'bad values in wza')
331 0 : if (any(abs(thza)>1.e20_r8)) call endrun(prefix//'bad values in thza')
332 :
333 : ! diagnostic output
334 0 : do j = 1,nzalat
335 0 : call outfld('Uzm',uza(j,:),1,j)
336 0 : call outfld('Vzm',vza(j,:),1,j)
337 0 : call outfld('Wzm',wza(j,:),1,j)
338 0 : call outfld('THzm',thza(j,:),1,j)
339 0 : call outfld('UVzm',uvza(j,:),1,j)
340 0 : call outfld('UWzm',uwza(j,:),1,j)
341 0 : call outfld('VTHzm',vthza(j,:),1,j)
342 0 : call outfld('WTHzm',wthza(j,:),1,j)
343 : end do
344 :
345 : contains
346 :
347 : !------------------------------------------------------------------------------
348 : ! utility function for evaluating 3D zonal mean fields
349 : !------------------------------------------------------------------------------
350 0 : function zmean_fld( fld ) result(fldzm)
351 :
352 : real(r8), intent(in) :: fld(pcols,pver,begchunk:endchunk)
353 :
354 : real(r8) :: fldzm(pcols,pver,begchunk:endchunk)
355 :
356 0 : real(r8) :: Zonal_Bamp3d(nzmbas,pver)
357 :
358 0 : call ZMobj%calc_amps(fld,Zonal_Bamp3d)
359 0 : call ZMobj%eval_grid(Zonal_Bamp3d,fldzm)
360 :
361 0 : end function zmean_fld
362 :
363 : !------------------------------------------------------------------------------
364 : ! utility function returns TRUE when time to update TEM diags
365 : !------------------------------------------------------------------------------
366 370944 : logical function do_calc()
367 :
368 : integer :: nstep
369 370944 : nstep = get_nstep()
370 370944 : do_calc = do_tem_diags .and. mod(nstep,ntimesteps) == 0
371 :
372 370944 : end function do_calc
373 :
374 : end subroutine phys_grid_ctem_diags
375 :
376 : !-----------------------------------------------------------------------------
377 : !-----------------------------------------------------------------------------
378 1536 : subroutine phys_grid_ctem_final
379 1536 : call ZAobj%final()
380 1536 : call ZMobj%final()
381 1536 : end subroutine phys_grid_ctem_final
382 :
383 : end module phys_grid_ctem
|
