Line data Source code
1 : !-----------------------------------------------------------------------------
2 : ! circulation diagnostics -- terms of the Transformed Eulerian Mean (TEM) equation
3 : !-----------------------------------------------------------------------------
4 : module ctem
5 :
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use spmd_utils, only: masterproc
8 : use pmgrid, only: plon, plev, plevp
9 : use cam_logfile, only: iulog
10 : use cam_history, only: addfld, outfld, add_default, horiz_only
11 : use cam_abortutils, only: endrun
12 :
13 : implicit none
14 : private
15 :
16 : public :: ctem_readnl
17 : public :: ctem_init
18 : public :: ctem_diags
19 : public :: do_circulation_diags
20 :
21 : real(r8) :: rplon
22 : real(r8) :: iref_p(plevp) ! interface reference pressure for vertical interpolation
23 : integer :: ip_b ! level index where hybrid levels become purely pressure
24 : integer :: zm_limit
25 :
26 : logical :: do_circulation_diags = .false.
27 :
28 : contains
29 :
30 : !================================================================================
31 :
32 0 : subroutine ctem_diags( u3, v3, omga, pt, h2o, ps, pe, grid)
33 :
34 : use physconst, only : zvir, cappa
35 : use dynamics_vars, only : T_FVDYCORE_GRID
36 : use hycoef, only : ps0
37 : use interpolate_data, only : vertinterp
38 : #ifdef SPMD
39 : use mpishorthand, only : mpilog, mpiint
40 : use parutilitiesmodule, only : pargatherint
41 : #endif
42 :
43 : !-------------------------------------------------------------
44 : ! ... dummy arguments
45 : !-------------------------------------------------------------
46 : type(T_FVDYCORE_GRID), intent(in) :: grid ! FV Dynamics grid
47 :
48 : real(r8), intent(in) :: ps(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! surface pressure (pa)
49 : real(r8), intent(in) :: u3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! zonal velocity (m/s)
50 : real(r8), intent(in) :: v3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! meridional velocity (m/s)
51 : real(r8), intent(in) :: omga(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! pressure velocity
52 : real(r8), intent(in) :: pe(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy) ! interface pressure (pa)
53 : real(r8), intent(in) :: pt(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,plev) ! virtual temperature
54 : real(r8), intent(in) :: h2o(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,plev) ! water constituent (kg/kg)
55 :
56 : !-------------------------------------------------------------
57 : ! ... local variables
58 : !-------------------------------------------------------------
59 : real(r8), parameter :: hscale = 7000._r8 ! pressure scale height
60 : real(r8), parameter :: navp = 1.e35_r8
61 :
62 : real(r8) :: pinterp
63 0 : real(r8) :: w(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! vertical velocity
64 0 : real(r8) :: th(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! pot. temperature
65 :
66 0 : real(r8) :: pm(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! mid-point pressure
67 :
68 : real(r8) :: pexf ! Exner function
69 : real(r8) :: psurf
70 :
71 0 : real(r8) :: ui(grid%ifirstxy:grid%ilastxy,plevp) ! interpolated zonal velocity
72 0 : real(r8) :: vi(grid%ifirstxy:grid%ilastxy,plevp) ! interpolated meridional velocity
73 0 : real(r8) :: wi(grid%ifirstxy:grid%ilastxy,plevp) ! interpolated vertical velocity
74 0 : real(r8) :: thi(grid%ifirstxy:grid%ilastxy,plevp) ! interpolated pot. temperature
75 :
76 : real(r8) :: um(plevp) ! zonal mean zonal velocity
77 : real(r8) :: vm(plevp) ! zonal mean meridional velocity
78 : real(r8) :: wm(plevp) ! zonal mean vertical velocity
79 : real(r8) :: thm(plevp) ! zonal mean pot. temperature
80 :
81 0 : real(r8) :: ud(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of zonal velocity
82 0 : real(r8) :: vd(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of meridional velocity
83 0 : real(r8) :: wd(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of vertical velocity
84 0 : real(r8) :: thd(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of pot. temperature
85 :
86 0 : real(r8) :: vthp(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of zonal velocity
87 0 : real(r8) :: wthp(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of meridional velocity
88 0 : real(r8) :: uvp(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of vertical velocity
89 0 : real(r8) :: uwp(grid%ifirstxy:grid%ilastxy,plevp) ! zonal deviation of pot. temperature
90 :
91 : real(r8) :: rdiv(plevp)
92 :
93 0 : integer :: ip_gm1g(plon,grid%jfirstxy:grid%jlastxy) ! contains level index-1 where blocked points begin
94 : integer :: zm_cnt(plevp) ! counter
95 : integer :: i,j,k
96 : integer :: nlons
97 :
98 0 : logical :: has_zm(plevp,grid%jfirstxy:grid%jlastxy) ! .true. the (z,y) point is a valid zonal mean
99 0 : integer :: ip_gm1(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! contains level index-1 where blocked points begin
100 :
101 0 : real(r8) :: vth(plevp,grid%jfirstxy:grid%jlastxy) ! VTH flux
102 0 : real(r8) :: uv(plevp,grid%jfirstxy:grid%jlastxy) ! UV flux
103 0 : real(r8) :: wth(plevp,grid%jfirstxy:grid%jlastxy) ! WTH flux
104 0 : real(r8) :: uw(plevp,grid%jfirstxy:grid%jlastxy) ! UW flux
105 0 : real(r8) :: u2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged U
106 0 : real(r8) :: v2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged V
107 0 : real(r8) :: th2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged TH
108 0 : real(r8) :: w2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged W
109 0 : real(r8) :: thig(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy) ! interpolated pot. temperature
110 :
111 0 : real(r8) :: tmp2(grid%ifirstxy:grid%ilastxy)
112 0 : real(r8) :: tmp3(grid%ifirstxy:grid%ilastxy,plevp)
113 :
114 : integer :: beglat, endlat ! begin,end latitude indicies
115 : integer :: beglon, endlon ! begin,end longitude indicies
116 :
117 0 : beglon = grid%ifirstxy
118 0 : endlon = grid%ilastxy
119 0 : beglat = grid%jfirstxy
120 0 : endlat = grid%jlastxy
121 :
122 : !omp parallel do private (i,j,k,pexf,psurf)
123 : lat_loop1 : &
124 0 : do j = beglat, endlat
125 0 : do k = 1, plev
126 0 : do i = beglon, endlon
127 : !-------------------------------------------------------------
128 : ! Calculate pressure and Exner function
129 : !-------------------------------------------------------------
130 0 : pm(i,k,j) = 0.5_r8 * ( pe(i,k,j) + pe(i,k+1,j) )
131 0 : pexf = (ps0 / pm(i,k,j))**cappa
132 : !-------------------------------------------------------------
133 : ! Convert virtual temperature to temperature and calculate potential temperature
134 : !-------------------------------------------------------------
135 0 : th(i,k,j) = pt(i,j,k) / (1._r8 + zvir*h2o(i,j,k))
136 0 : th(i,k,j) = th(i,k,j) * pexf
137 : !-------------------------------------------------------------
138 : ! Calculate vertical velocity
139 : !-------------------------------------------------------------
140 0 : w(i,k,j) = - hscale * omga(i,k,j) / pm(i,k,j)
141 : end do
142 : end do
143 : !-------------------------------------------------------------
144 : ! Keep track of where the bottom is in each column
145 : ! (i.e., largest index for which P(k) <= PS)
146 : !-------------------------------------------------------------
147 0 : ip_gm1(:,j) = plevp
148 0 : do i = beglon, endlon
149 0 : psurf = ps(i,j)
150 0 : do k = ip_b+1, plevp
151 0 : if( iref_p(k) <= psurf ) then
152 0 : ip_gm1(i,j) = k
153 : end if
154 : end do
155 : end do
156 : end do lat_loop1
157 :
158 0 : nlons = endlon - beglon + 1
159 :
160 : #ifdef SPMD
161 0 : if( grid%twod_decomp == 1 ) then
162 0 : if (grid%iam .lt. grid%npes_xy) then
163 0 : call pargatherint( grid%commxy_x, 0, ip_gm1, grid%strip2dx, ip_gm1g )
164 : endif
165 : else
166 0 : ip_gm1g(:,:) = ip_gm1(:,:)
167 : end if
168 : #else
169 : ip_gm1g(:,:) = ip_gm1(:,:)
170 : #endif
171 : #ifdef CTEM_DIAGS
172 : write(iulog,*) '===================================================='
173 : do j = beglat,endlat
174 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
175 : write(iulog,'(20i3)') ip_gm1(:,j)
176 : end do
177 : if( grid%myidxy_x == 0 ) then
178 : do j = beglat,endlat
179 : write(iulog,*) '===================================================='
180 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
181 : write(iulog,'(20i3)') ip_gm1g(:,j)
182 : end do
183 : write(iulog,*) '===================================================='
184 : #else
185 : #ifdef SPMD
186 0 : if( grid%myidxy_x == 0 ) then
187 : #endif
188 : #endif
189 : lat_loop2 : &
190 0 : do j = beglat, endlat
191 0 : zm_cnt(:ip_b) = plon
192 0 : do k = ip_b+1, plevp
193 0 : zm_cnt(k) = count( ip_gm1g(:,j) >= k )
194 : end do
195 0 : has_zm(:ip_b,j) = .true.
196 0 : do k = ip_b+1, plevp
197 0 : has_zm(k,j) = zm_cnt(k) >= zm_limit
198 : end do
199 : end do lat_loop2
200 : #ifdef SPMD
201 : end if
202 0 : if( grid%twod_decomp == 1 ) then
203 0 : call mpibcast( has_zm, plevp*(endlat-beglat+1), mpilog, 0, grid%commxy_x )
204 0 : call mpibcast( ip_gm1g, plon*(endlat-beglat+1), mpiint, 0, grid%commxy_x )
205 : end if
206 : #endif
207 :
208 : #ifdef CTEM_DIAGS
209 : if( grid%myidxy_y == 12 ) then
210 : write(iulog,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
211 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,beglat
212 : write(iulog,*) 'has_zm'
213 : write(iulog,'(20l2)') has_zm(:,beglat)
214 : write(iulog,*) 'ip_gm1g'
215 : write(iulog,'(20i4)') ip_gm1g(:,beglat)
216 : write(iulog,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
217 : end if
218 : #endif
219 :
220 : lat_loop3 : &
221 0 : do j = beglat, endlat
222 : !-------------------------------------------------------------
223 : ! Vertical interpolation
224 : !-------------------------------------------------------------
225 0 : do k = 1, plevp
226 0 : pinterp = iref_p(k)
227 : !-------------------------------------------------------------
228 : ! Zonal velocity
229 : !-------------------------------------------------------------
230 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
231 0 : u3(beglon,1,j), ui(beglon,k) )
232 : !-------------------------------------------------------------
233 : ! Meridional velocity
234 : !-------------------------------------------------------------
235 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
236 0 : v3(beglon,1,j), vi(beglon,k) )
237 : !-------------------------------------------------------------
238 : ! Vertical velocity
239 : !-------------------------------------------------------------
240 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
241 0 : w(beglon,1,j), wi(beglon,k) )
242 : !-------------------------------------------------------------
243 : ! Pot. Temperature
244 : !-------------------------------------------------------------
245 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
246 0 : th(beglon,1,j), thi(beglon,k) )
247 : end do
248 : #ifdef CTEM_DIAGS
249 : if( j == endlat ) then
250 : write(iulog,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
251 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
252 : write(iulog,*) 'iref_p'
253 : write(iulog,'(5g15.7)') iref_p(:)
254 : write(iulog,'(''pm(endlon,:,'',i2,'')'')') j
255 : write(iulog,'(5g15.7)') pm(endlon,:,j)
256 : write(iulog,'(''u3(endlon,:,'',i2,'')'')') j
257 : write(iulog,'(5g15.7)') u3(endlon,:,j)
258 : write(iulog,*) 'ui(endlon,:)'
259 : write(iulog,'(5g15.7)') ui(endlon,:)
260 : write(iulog,*) '^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^'
261 : end if
262 : #endif
263 :
264 : !-------------------------------------------------------------
265 : ! Calculate zonal averages
266 : !-------------------------------------------------------------
267 0 : do k = ip_b+1, plevp
268 0 : if( has_zm(k,j) ) then
269 0 : where( ip_gm1(beglon:endlon,j) < k )
270 : ui(beglon:endlon,k) = 0._r8
271 : vi(beglon:endlon,k) = 0._r8
272 : wi(beglon:endlon,k) = 0._r8
273 : thi(beglon:endlon,k) = 0._r8
274 : endwhere
275 : end if
276 : end do
277 :
278 0 : call par_xsum( grid, ui, plevp, um )
279 0 : call par_xsum( grid, vi, plevp, vm )
280 0 : call par_xsum( grid, wi, plevp, wm )
281 0 : call par_xsum( grid, thi, plevp, thm )
282 : #ifdef CTEM_DIAGS
283 : if( j == endlat .and. grid%myidxy_y == 12 ) then
284 : write(iulog,*) '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
285 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
286 : write(iulog,*) 'um after par_xsum'
287 : write(iulog,'(5g15.7)') um(:)
288 : write(iulog,*) '$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$'
289 : end if
290 : #endif
291 0 : do k = 1, ip_b
292 0 : um(k) = um(k) * rplon
293 0 : vm(k) = vm(k) * rplon
294 0 : wm(k) = wm(k) * rplon
295 0 : thm(k) = thm(k) * rplon
296 0 : u2d(k,j) = um(k)
297 0 : v2d(k,j) = vm(k)
298 0 : th2d(k,j) = thm(k)
299 0 : w2d(k,j) = wm(k)
300 : end do
301 0 : do k = ip_b+1, plevp
302 0 : if( has_zm(k,j) ) then
303 0 : rdiv(k) = 1._r8/count( ip_gm1g(:,j) >= k )
304 0 : um(k) = um(k) * rdiv(k)
305 0 : vm(k) = vm(k) * rdiv(k)
306 0 : wm(k) = wm(k) * rdiv(k)
307 0 : thm(k) = thm(k) * rdiv(k)
308 0 : u2d(k,j) = um(k)
309 0 : v2d(k,j) = vm(k)
310 0 : th2d(k,j) = thm(k)
311 0 : w2d(k,j) = wm(k)
312 : else
313 0 : u2d(k,j) = navp
314 0 : v2d(k,j) = navp
315 0 : th2d(k,j) = navp
316 0 : w2d(k,j) = navp
317 : end if
318 : end do
319 :
320 : !-------------------------------------------------------------
321 : ! Calculate zonal deviations
322 : !-------------------------------------------------------------
323 0 : do k = 1, ip_b
324 0 : ud(beglon:endlon,k) = ui(beglon:endlon,k) - um(k)
325 0 : vd(beglon:endlon,k) = vi(beglon:endlon,k) - vm(k)
326 0 : wd(beglon:endlon,k) = wi(beglon:endlon,k) - wm(k)
327 0 : thd(beglon:endlon,k) = thi(beglon:endlon,k) - thm(k)
328 : end do
329 :
330 0 : do k = ip_b+1, plevp
331 0 : if( has_zm(k,j) ) then
332 0 : where( ip_gm1g(beglon:endlon,j) >= k )
333 : ud(beglon:endlon,k) = ui(beglon:endlon,k) - um(k)
334 : vd(beglon:endlon,k) = vi(beglon:endlon,k) - vm(k)
335 : wd(beglon:endlon,k) = wi(beglon:endlon,k) - wm(k)
336 : thd(beglon:endlon,k) = thi(beglon:endlon,k) - thm(k)
337 : elsewhere
338 : ud(beglon:endlon,k) = 0._r8
339 : vd(beglon:endlon,k) = 0._r8
340 : wd(beglon:endlon,k) = 0._r8
341 : thd(beglon:endlon,k) = 0._r8
342 : endwhere
343 : end if
344 : end do
345 :
346 : !-------------------------------------------------------------
347 : ! Calculate fluxes
348 : !-------------------------------------------------------------
349 0 : do k = 1, ip_b
350 0 : vthp(:,k) = vd(:,k) * thd(:,k)
351 0 : wthp(:,k) = wd(:,k) * thd(:,k)
352 0 : uwp(:,k) = wd(:,k) * ud(:,k)
353 0 : uvp(:,k) = vd(:,k) * ud(:,k)
354 : end do
355 :
356 0 : do k = ip_b+1, plevp
357 0 : if( has_zm(k,j) ) then
358 0 : vthp(:,k) = vd(:,k) * thd(:,k)
359 0 : wthp(:,k) = wd(:,k) * thd(:,k)
360 0 : uwp(:,k) = wd(:,k) * ud(:,k)
361 0 : uvp(:,k) = vd(:,k) * ud(:,k)
362 : else
363 0 : vthp(:,k) = 0._r8
364 0 : wthp(:,k) = 0._r8
365 0 : uwp(:,k) = 0._r8
366 0 : uvp(:,k) = 0._r8
367 : end if
368 : end do
369 :
370 : #ifdef CTEM_DIAGS
371 : if( j == endlat .and. grid%myidxy_y == 12 ) then
372 : write(iulog,*) '#################################################'
373 : write(iulog,*) 'DIAGNOSTICS before par_xsum'
374 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
375 : write(iulog,*) 'has_zm'
376 : write(iulog,*) has_zm(:,j)
377 : write(iulog,*) 'rdiv'
378 : write(iulog,'(5g15.7)') rdiv(:)
379 : write(iulog,*) 'wm'
380 : write(iulog,'(5g15.7)') wm(:)
381 : write(iulog,*) 'um'
382 : write(iulog,'(5g15.7)') um(:)
383 : write(iulog,*) 'uw'
384 : write(iulog,'(5g15.7)') uw(:)
385 : write(iulog,*) '#################################################'
386 : end if
387 : #endif
388 0 : call par_xsum( grid, vthp, plevp, vth(1,j) )
389 0 : call par_xsum( grid, wthp, plevp, wth(1,j) )
390 0 : call par_xsum( grid, uvp, plevp, uv(1,j) )
391 0 : call par_xsum( grid, uwp, plevp, uw(1,j) )
392 : #ifdef CTEM_DIAGS
393 : if( j == endlat .and. grid%myidxy_y == 12 ) then
394 : write(iulog,*) '#################################################'
395 : write(iulog,'(''iam,myidxy_x,myidxy_y,j = '',4i4)') grid%iam,grid%myidxy_x,grid%myidxy_y,j
396 : write(iulog,*) 'uw after par_xsum'
397 : write(iulog,'(5g15.7)') uw(:,j)
398 : write(iulog,*) '#################################################'
399 : end if
400 : #endif
401 0 : do k = 1, ip_b
402 0 : vth(k,j) = vth(k,j) * rplon
403 0 : wth(k,j) = wth(k,j) * rplon
404 0 : uw(k,j) = uw(k,j) * rplon
405 0 : uv(k,j) = uv(k,j) * rplon
406 : end do
407 0 : do k = ip_b+1, plevp
408 0 : if( has_zm(k,j) ) then
409 0 : vth(k,j) = vth(k,j) * rdiv(k)
410 0 : wth(k,j) = wth(k,j) * rdiv(k)
411 0 : uw(k,j) = uw(k,j) * rdiv(k)
412 0 : uv(k,j) = uv(k,j) * rdiv(k)
413 : else
414 0 : vth(k,j) = navp
415 0 : wth(k,j) = navp
416 0 : uw(k,j) = navp
417 0 : uv(k,j) = navp
418 : end if
419 : end do
420 :
421 0 : thig(:,:,j) = thi(:,:)
422 : end do lat_loop3
423 :
424 : !-------------------------------------------------------------
425 : ! Do the output
426 : !-------------------------------------------------------------
427 0 : latloop: do j = beglat,endlat
428 :
429 : !-------------------------------------------------------------
430 : ! zonal-mean output
431 : !-------------------------------------------------------------
432 0 : do k = 1,plevp
433 0 : tmp3(grid%ifirstxy,k) = vth(k,j)
434 : enddo
435 0 : call outfld( 'VTHzm', tmp3(grid%ifirstxy,:), 1, j )
436 :
437 0 : do k = 1,plevp
438 0 : tmp3(grid%ifirstxy,k) = wth(k,j)
439 : enddo
440 0 : call outfld( 'WTHzm', tmp3(grid%ifirstxy,:), 1, j )
441 :
442 0 : do k = 1,plevp
443 0 : tmp3(grid%ifirstxy,k) = uv(k,j)
444 : enddo
445 0 : call outfld( 'UVzm', tmp3(grid%ifirstxy,:), 1, j )
446 :
447 0 : do k = 1,plevp
448 0 : tmp3(grid%ifirstxy,k) = uw(k,j)
449 : enddo
450 0 : call outfld( 'UWzm', tmp3(grid%ifirstxy,:), 1, j )
451 0 : do k = 1,plevp
452 0 : tmp3(grid%ifirstxy,k) = u2d(k,j)
453 : enddo
454 0 : call outfld( 'Uzm', tmp3(grid%ifirstxy,:), 1, j )
455 0 : do k = 1,plevp
456 0 : tmp3(grid%ifirstxy,k) = v2d(k,j)
457 : enddo
458 0 : call outfld( 'Vzm', tmp3(grid%ifirstxy,:), 1, j )
459 0 : do k = 1,plevp
460 0 : tmp3(grid%ifirstxy,k) = w2d(k,j)
461 : enddo
462 0 : call outfld( 'Wzm', tmp3(grid%ifirstxy,:), 1, j )
463 0 : do k = 1,plevp
464 0 : tmp3(grid%ifirstxy,k) = th2d(k,j)
465 : enddo
466 0 : call outfld( 'THzm', tmp3(grid%ifirstxy,:), 1, j )
467 :
468 : !-------------------------------------------------------------
469 : ! 3D output
470 : !-------------------------------------------------------------
471 0 : do k = 1,plevp
472 0 : do i = beglon,endlon
473 0 : tmp3(i,k) = thig(i,k,j)
474 : enddo
475 : enddo
476 0 : call outfld( 'TH', tmp3, nlons, j )
477 :
478 : !-------------------------------------------------------------
479 : ! horizontal output
480 : !-------------------------------------------------------------
481 0 : tmp2(beglon:endlon) = ip_gm1(beglon:endlon,j)
482 0 : call outfld( 'MSKtem', tmp2, nlons, j )
483 :
484 : enddo latloop
485 :
486 0 : end subroutine ctem_diags
487 :
488 : !=================================================================================
489 :
490 1536 : subroutine ctem_init()
491 :
492 : use hycoef, only : hyai, hybi, ps0
493 : use phys_control, only : phys_getopts
494 :
495 : !-------------------------------------------------------------
496 : ! ... local variables
497 : !-------------------------------------------------------------
498 : integer :: k
499 : logical :: history_waccm
500 :
501 1536 : if (.not.do_circulation_diags) return
502 :
503 0 : rplon = 1._r8/plon
504 0 : zm_limit = plon/3
505 :
506 : !-------------------------------------------------------------
507 : ! Calculate reference pressure
508 : !-------------------------------------------------------------
509 0 : do k = 1, plevp
510 0 : iref_p(k) = (hyai(k) + hybi(k)) * ps0
511 : end do
512 0 : if( masterproc ) then
513 0 : write(iulog,*) 'ctem_inti: iref_p'
514 0 : write(iulog,'(1p5g15.7)') iref_p(:)
515 : end if
516 :
517 : !-------------------------------------------------------------
518 : ! Find level where hybrid levels become purely pressure
519 : !-------------------------------------------------------------
520 0 : ip_b = -1
521 0 : do k = 1,plev
522 0 : if( hybi(k) == 0._r8 ) ip_b = k
523 : end do
524 :
525 0 : call phys_getopts( history_waccm_out = history_waccm )
526 :
527 : !-------------------------------------------------------------
528 : ! Initialize output buffer
529 : !-------------------------------------------------------------
530 0 : call addfld ('VTHzm',(/ 'ilev' /),'A','MK/S','Meridional Heat Flux: 3D zon. mean', gridname='fv_centers_zonal' )
531 0 : call addfld ('WTHzm',(/ 'ilev' /),'A','MK/S','Vertical Heat Flux: 3D zon. mean', gridname='fv_centers_zonal' )
532 0 : call addfld ('UVzm', (/ 'ilev' /),'A','M2/S2','Meridional Flux of Zonal Momentum: 3D zon. mean', gridname='fv_centers_zonal' )
533 0 : call addfld ('UWzm', (/ 'ilev' /),'A','M2/S2','Vertical Flux of Zonal Momentum: 3D zon. mean', gridname='fv_centers_zonal' )
534 :
535 0 : call addfld ('Uzm', (/ 'ilev' /),'A','M/S','Zonal-Mean zonal wind - defined on ilev', gridname='fv_centers_zonal' )
536 0 : call addfld ('Vzm', (/ 'ilev' /),'A','M/S','Zonal-Mean meridional wind - defined on ilev', gridname='fv_centers_zonal' )
537 0 : call addfld ('Wzm', (/ 'ilev' /),'A','M/S','Zonal-Mean vertical wind - defined on ilev', gridname='fv_centers_zonal' )
538 0 : call addfld ('THzm', (/ 'ilev' /),'A', 'K','Zonal-Mean potential temp - defined on ilev', gridname='fv_centers_zonal' )
539 :
540 0 : call addfld ('TH', (/ 'ilev' /),'A','K', 'Potential Temperature', gridname='fv_centers' )
541 0 : call addfld ('MSKtem',horiz_only, 'A','1', 'TEM mask', gridname='fv_centers' )
542 :
543 : !-------------------------------------------------------------
544 : ! primary tapes: 3D fields
545 : !-------------------------------------------------------------
546 0 : call add_default ('VTHzm', 1, ' ')
547 0 : call add_default ('WTHzm', 1, ' ')
548 0 : call add_default ('UVzm' , 1, ' ')
549 0 : call add_default ('UWzm' , 1, ' ')
550 0 : call add_default ('TH' , 1, ' ')
551 0 : call add_default ('MSKtem',1, ' ')
552 :
553 0 : if (history_waccm) then
554 0 : call add_default ('MSKtem',7, ' ')
555 0 : call add_default ('VTHzm', 7, ' ')
556 0 : call add_default ('UVzm', 7, ' ')
557 0 : call add_default ('UWzm', 7, ' ')
558 0 : call add_default ('Uzm', 7, ' ')
559 0 : call add_default ('Vzm', 7, ' ')
560 0 : call add_default ('Wzm', 7, ' ')
561 0 : call add_default ('THzm', 7, ' ')
562 : end if
563 :
564 0 : if (masterproc) then
565 0 : write(iulog,*) 'ctem_inti: do_circulation_diags = ',do_circulation_diags
566 : endif
567 :
568 : end subroutine ctem_init
569 :
570 : !================================================================================
571 :
572 1536 : subroutine ctem_readnl(nlfile)
573 :
574 : use namelist_utils, only: find_group_name
575 : use units, only: getunit, freeunit
576 : use spmd_utils, only: mpicom, mstrid=>masterprocid, mpi_logical
577 :
578 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
579 :
580 : ! Local variables
581 : integer :: unitn, ierr
582 : character(len=*), parameter :: subname = 'ctem_readnl'
583 :
584 : namelist /circ_diag_nl/ do_circulation_diags
585 : !-----------------------------------------------------------------------------
586 :
587 1536 : if (masterproc) then
588 2 : unitn = getunit()
589 2 : open( unitn, file=trim(nlfile), status='old' )
590 2 : call find_group_name(unitn, 'circ_diag_nl', status=ierr)
591 2 : if (ierr == 0) then
592 0 : read(unitn, circ_diag_nl, iostat=ierr)
593 0 : if (ierr /= 0) then
594 0 : call endrun(subname // ':: ERROR reading namelist')
595 : end if
596 : end if
597 2 : close(unitn)
598 2 : call freeunit(unitn)
599 : end if
600 :
601 1536 : call mpi_bcast(do_circulation_diags, 1, mpi_logical, mstrid, mpicom, ierr)
602 1536 : if (ierr /= 0) call endrun(subname//": FATAL: mpi_bcast: do_circulation_diags")
603 :
604 1536 : end subroutine ctem_readnl
605 :
606 : end module ctem
|
