Line data Source code
1 : module dp_coupling
2 : !BOP
3 : !
4 : ! !MODULE: dp_coupling --- dynamics-physics coupling module
5 : !
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use ppgrid, only: pcols, pver
8 : use phys_grid
9 :
10 : use physics_types, only: physics_state, physics_tend, physics_cnst_limit
11 : use constituents, only: pcnst
12 : use physconst, only: gravit, zvir
13 : use air_composition, only: cpairv, rairv
14 : use geopotential, only: geopotential_t
15 : use check_energy, only: check_energy_timestep_init
16 : use dynamics_vars, only: T_FVDYCORE_GRID, t_fvdycore_state
17 : use dyn_internal_state,only: get_dyn_state
18 : use dyn_comp, only: dyn_import_t, dyn_export_t, fv_print_dpcoup_warn
19 : use cam_abortutils, only: endrun
20 : #if defined ( SPMD )
21 : use spmd_dyn, only: local_dp_map, block_buf_nrecs, chunk_buf_nrecs
22 : #endif
23 : use perf_mod
24 : use cam_logfile, only: iulog
25 : !
26 : ! !PUBLIC MEMBER FUNCTIONS:
27 : PUBLIC d_p_coupling, p_d_coupling
28 :
29 : !
30 : ! !DESCRIPTION:
31 : !
32 : ! This module provides
33 : !
34 : ! \begin{tabular}{|l|l|} \hline \hline
35 : ! d\_p\_coupling & dynamics output to physics input \\ \hline
36 : ! p\_d\_coupling & physics output to dynamics input \\ \hline
37 : ! \hline
38 : ! \end{tabular}
39 : !
40 : ! !REVISION HISTORY:
41 : ! 00.06.01 Boville Creation
42 : ! 01.10.01 Lin Various revisions
43 : ! 01.03.26 Sawyer Added ProTeX documentation
44 : ! 01.06.27 Mirin Separate noncoupling coding into new routines
45 : ! 01.07.13 Mirin Some support for multi-2D decompositions
46 : ! 02.03.01 Worley Support for nontrivial physics remapping
47 : ! 03.03.28 Boville set all physics_state elements, add check_energy_timestep_init
48 : ! 03.08.13 Sawyer Removed ghost N1 region in u3sxy
49 : ! 05.06.28 Sawyer Simplified interfaces -- only XY decomposition
50 : ! 05.10.25 Sawyer Extensive refactoring, dyn_interface
51 : ! 05.11.10 Sawyer Now using dyn_import/export_t containers
52 : ! 06.07.01 Sawyer Transitioned constituents to T_TRACERS
53 : !
54 : !EOP
55 : !-----------------------------------------------------------------------
56 :
57 : private
58 : real(r8), parameter :: D0_5 = 0.5_r8
59 : real(r8), parameter :: D1_0 = 1.0_r8
60 :
61 : CONTAINS
62 :
63 : !-----------------------------------------------------------------------
64 : !BOP
65 : ! !IROUTINE: d_p_coupling --- convert dynamics output to physics input
66 : !
67 : ! !INTERFACE:
68 46848 : subroutine d_p_coupling(grid, phys_state, phys_tend, pbuf2d, dyn_out)
69 :
70 : ! !USES:
71 : use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, &
72 : pbuf_get_field
73 : use constituents, only: qmin
74 : use physics_types, only: set_state_pdry, set_wet_to_dry
75 :
76 : use pmgrid, only: plev
77 : use ctem, only: ctem_diags, do_circulation_diags
78 : use diag_module, only: fv_diag_am_calc
79 : use gravity_waves_sources, only: gws_src_fnct
80 : use cam_thermo, only: cam_thermo_dry_air_update
81 : use shr_const_mod, only: shr_const_rwv
82 : use dyn_comp, only: frontgf_idx, frontga_idx, uzm_idx
83 : use qbo, only: qbo_use_forcing
84 : use phys_control, only: use_gw_front, use_gw_front_igw, waccmx_is
85 : use zonal_mean, only: zonal_mean_3D
86 : use d2a3dikj_mod, only: d2a3dikj
87 : use qneg_module, only: qneg3
88 : use air_composition,only: dry_air_species_num
89 : !-----------------------------------------------------------------------
90 : implicit none
91 : !-----------------------------------------------------------------------
92 : ! !INPUT PARAMETERS:
93 : !
94 : type(t_fvdycore_grid), intent(in) :: grid
95 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
96 : type(dyn_export_t), intent(in) :: dyn_out ! dynamics export
97 :
98 : ! !OUTPUT PARAMETERS:
99 :
100 : type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
101 : type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
102 :
103 :
104 : ! !DESCRIPTION:
105 : !
106 : ! Coupler for converting dynamics output variables into physics
107 : ! input variables
108 : !
109 : ! !REVISION HISTORY:
110 : ! 00.06.01 Boville Creation
111 : ! 01.07.13 AAM Some support for multi-2D decompositions
112 : ! 02.03.01 Worley Support for nontrivial physics remapping
113 : ! 02.05.02 Sawyer u3s made inout due to ghosting in d2a3dikj
114 : ! 03.08.05 Sawyer Removed pe11k, pe11kln (for defunct Rayl fric)
115 : ! 04.08.29 Eaton Added lat, lon coords to physics_state type
116 : ! 05.06.28 Sawyer Simplified interface -- on XY decomp vars.
117 : ! 05.07.06 Sawyer Added dyn_state as argument
118 : ! 05.10.31 Sawyer Refactoring, replaced dyn_state by dyn_interface
119 : !
120 : !EOP
121 : !-----------------------------------------------------------------------
122 : !BOC
123 : ! !LOCAL VARIABLES:
124 :
125 : type(t_fvdycore_state), pointer :: dyn_state
126 :
127 : ! Variables from dynamics export container
128 16128 : real(r8), pointer :: phisxy(:,:) ! surface geopotential
129 16128 : real(r8), pointer :: psxy (:,:) ! surface pressure
130 16128 : real(r8), pointer :: u3sxy(:,:,:) ! u-wind on d-grid
131 16128 : real(r8), pointer :: v3sxy(:,:,:) ! v-wind on d-grid
132 16128 : real(r8), pointer :: du3sxy(:,:,:) ! u-wind increment on d-grid
133 16128 : real(r8), pointer :: dv3sxy(:,:,:) ! v-wind increment on d-grid
134 16128 : real(r8), pointer :: dua3sxy(:,:,:) ! u-wind adv. inc. on d-grid
135 16128 : real(r8), pointer :: dva3sxy(:,:,:) ! v-wind adv. inc. on d-grid
136 16128 : real(r8), pointer :: duf3sxy(:,:,:) ! u-wind fixer inc.on d-grid
137 16128 : real(r8), pointer :: ptxy (:,:,:) ! Virtual pot temp
138 16128 : real(r8), pointer :: tracer(:,:,:,:) ! constituents
139 16128 : real(r8), pointer :: omgaxy(:,:,:) ! vertical velocity
140 16128 : real(r8), pointer :: pexy (:,:,:) ! edge pressure
141 16128 : real(r8), pointer :: pelnxy(:,:,:) ! log(pe)
142 16128 : real(r8), pointer :: pkxy (:,:,:) ! pe**cappa
143 16128 : real(r8), pointer :: pkzxy (:,:,:) ! f-v mean of pk
144 :
145 : integer :: i,ib,j,k,m,lchnk ! indices
146 : integer :: ncol ! number of columns in current chunk
147 : integer :: lats(pcols) ! array of latitude indices
148 : integer :: lons(pcols) ! array of longitude indices
149 : integer :: blksiz ! number of columns in 2D block
150 : integer :: tsize ! amount of data per grid point passed to physics
151 16128 : integer, allocatable, dimension(:,:) :: bpter
152 : ! offsets into block buffer for packing data
153 : integer :: cpter(pcols,0:pver) ! offsets into chunk buffer for unpacking data
154 :
155 : real(r8) :: qmavl ! available q at level pver-1
156 : real(r8) :: dqreq ! q change at pver-1 required to remove q<qmin at pver
157 : real(r8) :: qbot ! bottom level q before change
158 : real(r8) :: qbotm1 ! bottom-1 level q before change
159 : real(r8) :: pic(pcols) ! ps**cappa
160 : real(r8) :: fraction
161 16128 : real(r8), allocatable :: u3(:, :, :) ! u-wind on a-grid
162 16128 : real(r8), allocatable :: v3(:, :, :) ! v-wind on a-grid
163 16128 : real(r8), allocatable :: du3 (:, :, :) ! u-wind increment on a-grid
164 16128 : real(r8), allocatable :: dv3 (:, :, :) ! v-wind increment on a-grid
165 16128 : real(r8), allocatable :: dua3(:, :, :) ! u-wind adv. inc. on a-grid
166 16128 : real(r8), allocatable :: dva3(:, :, :) ! v-wind adv. inc. on a-grid
167 16128 : real(r8), allocatable :: duf3(:, :, :) ! u-wind fixer inc.on a-grid
168 16128 : real(r8), allocatable :: dummy(:,:, :) ! dummy work array ~ v3sxy*0
169 :
170 16128 : real(r8), allocatable, dimension(:) :: bbuffer, cbuffer
171 : ! transpose buffers
172 :
173 : real(r8) :: zvirv(pcols,pver) ! Local zvir array pointer
174 :
175 16128 : real(r8) :: uzm(plev,grid%jfirstxy:grid%jlastxy) ! Zonal mean zonal wind
176 :
177 : integer :: km, kmp1, iam
178 : integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
179 : integer :: ic, jc
180 : integer :: astat
181 : integer :: boff
182 :
183 : !--------------------------------------------
184 : ! Variables needed for WACCM
185 : !--------------------------------------------
186 :
187 : ! frontogenesis function for gravity wave drag
188 16128 : real(r8), allocatable :: frontgf(:,:,:)
189 16128 : real(r8), pointer :: pbuf_frontgf(:,:)
190 : ! frontogenesis angle for gravity wave drag
191 16128 : real(r8), allocatable :: frontga(:,:,:)
192 16128 : real(r8), pointer :: pbuf_frontga(:,:)
193 : ! needed for qbo
194 16128 : real(r8), pointer :: pbuf_uzm(:,:)
195 : #if (! defined SPMD)
196 : integer :: block_buf_nrecs = 0
197 : integer :: chunk_buf_nrecs = 0
198 : logical :: local_dp_map=.true.
199 : #endif
200 16128 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
201 :
202 : !---------------------------End Local workspace-------------------------
203 :
204 32256 : dyn_state => get_dyn_state()
205 :
206 16128 : if (use_gw_front .or. use_gw_front_igw) then
207 :
208 80640 : allocate(frontgf(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy), stat=astat)
209 16128 : if( astat /= 0 ) then
210 0 : write(iulog,*) 'd_p_coupling: failed to allocate frontgf; error = ',astat
211 0 : call endrun
212 : end if
213 :
214 80640 : allocate(frontga(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy), stat=astat)
215 16128 : if( astat /= 0 ) then
216 0 : write(iulog,*) 'd_p_coupling: failed to allocate frontga; error = ',astat
217 0 : call endrun
218 : end if
219 :
220 : end if
221 :
222 16128 : nullify(pbuf_chnk)
223 16128 : nullify(pbuf_frontgf)
224 16128 : nullify(pbuf_frontga)
225 16128 : nullify(pbuf_uzm)
226 :
227 16128 : fraction = 0.1_r8
228 :
229 16128 : phisxy => dyn_out%phis
230 16128 : psxy => dyn_out%ps
231 16128 : u3sxy => dyn_out%u3s
232 16128 : v3sxy => dyn_out%v3s
233 16128 : ptxy => dyn_out%pt
234 16128 : tracer => dyn_out%tracer
235 :
236 16128 : omgaxy => dyn_out%omga
237 16128 : pexy => dyn_out%pe
238 16128 : pelnxy => dyn_out%peln
239 16128 : pkxy => dyn_out%pk
240 16128 : pkzxy => dyn_out%pkz
241 :
242 16128 : km = grid%km
243 16128 : kmp1 = km + 1
244 :
245 16128 : ifirstxy = grid%ifirstxy
246 16128 : ilastxy = grid%ilastxy
247 16128 : jfirstxy = grid%jfirstxy
248 16128 : jlastxy = grid%jlastxy
249 :
250 16128 : iam = grid%iam
251 : !-----------------------------------------------------------------------
252 : ! Transform dynamics staggered winds to physics grid (D=>A)
253 : !-----------------------------------------------------------------------
254 :
255 16128 : call t_startf ('d2a3dikj')
256 80640 : allocate (u3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))
257 64512 : allocate (v3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))
258 :
259 16128 : if (iam .lt. grid%npes_xy) then
260 16128 : call d2a3dikj(grid, dyn_state%am_geom_crrct, u3sxy, v3sxy, u3, v3)
261 : end if ! (iam .lt. grid%npes_xy)
262 :
263 16128 : call t_stopf ('d2a3dikj')
264 :
265 16128 : if ( do_circulation_diags ) then
266 0 : call t_startf('DP_CPLN_ctem')
267 0 : call ctem_diags( u3, v3, omgaxy, ptxy(:,jfirstxy:jlastxy,:), tracer(:,jfirstxy:jlastxy,:,1), &
268 0 : psxy, pexy, grid)
269 0 : call t_stopf('DP_CPLN_ctem')
270 : endif
271 :
272 16128 : if (dyn_state%am_diag) then
273 0 : du3sxy => dyn_out%du3s
274 0 : dv3sxy => dyn_out%dv3s
275 0 : dua3sxy => dyn_out%dua3s
276 0 : dva3sxy => dyn_out%dva3s
277 0 : duf3sxy => dyn_out%duf3s
278 0 : allocate (du3 (ifirstxy:ilastxy, km, jfirstxy:jlastxy))
279 0 : allocate (dv3 (ifirstxy:ilastxy, km, jfirstxy:jlastxy))
280 0 : allocate (dua3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))
281 0 : allocate (dva3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))
282 0 : allocate (duf3(ifirstxy:ilastxy, km, jfirstxy:jlastxy))
283 0 : allocate (dummy(ifirstxy:ilastxy,jfirstxy:jlastxy, km))
284 0 : du3(:,:,:) = 0._r8
285 0 : dv3(:,:,:) = 0._r8
286 0 : dua3(:,:,:) = 0._r8
287 0 : dva3(:,:,:) = 0._r8
288 0 : duf3(:,:,:) = 0._r8
289 0 : dummy(:,:,:) = 0._r8
290 :
291 0 : if (iam .lt. grid%npes_xy) then
292 : ! (note dummy use of dva3 hence call order matters)
293 0 : call d2a3dikj(grid, dyn_state%am_geom_crrct,duf3sxy, dummy, duf3 ,dva3)
294 0 : call d2a3dikj(grid, dyn_state%am_geom_crrct,dua3sxy, dva3sxy, dua3, dva3)
295 0 : call d2a3dikj(grid, dyn_state%am_geom_crrct, du3sxy, dv3sxy, du3 , dv3 )
296 : end if ! (iam .lt. grid%npes_xy)
297 :
298 0 : call t_startf('DP_CPLN_fv_am')
299 0 : call fv_diag_am_calc(grid, psxy, pexy, du3, dv3, dua3, dva3, duf3)
300 0 : call t_stopf('DP_CPLN_fv_am')
301 : endif
302 :
303 16128 : if ((iam .lt. grid%npes_xy) .and. (use_gw_front .or. use_gw_front_igw)) then
304 16128 : call t_startf('DP_CPLN_gw_sources')
305 16128 : call gws_src_fnct(grid, u3, v3, ptxy, tracer(:,jfirstxy:jlastxy,:,1), pexy, frontgf, frontga)
306 16128 : call t_stopf('DP_CPLN_gw_sources')
307 : end if
308 : if (qbo_use_forcing) then
309 : call zonal_mean_3D(grid, plev, u3, uzm)
310 : end if
311 :
312 : !-----------------------------------------------------------------------
313 : ! Copy data from dynamics data structure to physics data structure
314 : !-----------------------------------------------------------------------
315 : has_local_map : &
316 16128 : if (local_dp_map) then
317 :
318 : ! This declaration is too long; this parallel section needs some stuff
319 : ! pulled out into routines.
320 : !$omp parallel do private (lchnk, ncol, i, k, m, ic, jc, lons, lats, pic, pbuf_chnk, pbuf_uzm, pbuf_frontgf, pbuf_frontga)
321 : chnk_loop1 : &
322 0 : do lchnk = begchunk,endchunk
323 0 : ncol = phys_state(lchnk)%ncol
324 0 : call get_lon_all_p(lchnk, ncol, lons)
325 0 : call get_lat_all_p(lchnk, ncol, lats)
326 :
327 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
328 :
329 0 : if (use_gw_front .or. use_gw_front_igw) then
330 0 : call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf)
331 0 : call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga)
332 : end if
333 :
334 : if (qbo_use_forcing) then
335 : call pbuf_get_field(pbuf_chnk, uzm_idx, pbuf_uzm)
336 : end if
337 :
338 0 : do i=1,ncol
339 0 : ic = lons(i)
340 0 : jc = lats(i)
341 0 : phys_state(lchnk)%ps(i) = psxy(ic,jc)
342 0 : phys_state(lchnk)%phis(i) = phisxy(ic,jc)
343 0 : pic(i) = pkxy(ic,jc,pver+1)
344 : enddo
345 0 : do k=1,km
346 0 : do i=1,ncol
347 0 : ic = lons(i)
348 0 : jc = lats(i)
349 0 : phys_state(lchnk)%u (i,k) = u3(ic,k,jc)
350 0 : phys_state(lchnk)%v (i,k) = v3(ic,k,jc)
351 0 : phys_state(lchnk)%omega(i,k) = omgaxy(ic,k,jc)
352 0 : phys_state(lchnk)%t (i,k) = ptxy(ic,jc,k) / (D1_0 + zvir*tracer(ic,jc,k,1))
353 0 : phys_state(lchnk)%exner(i,k) = pic(i) / pkzxy(ic,jc,k)
354 :
355 0 : if (use_gw_front .or. use_gw_front_igw) then
356 0 : pbuf_frontgf(i,k) = frontgf(ic,k,jc)
357 0 : pbuf_frontga(i,k) = frontga(ic,k,jc)
358 : endif
359 :
360 0 : if (qbo_use_forcing) then
361 : pbuf_uzm(i,k) = uzm(k,jc)
362 : end if
363 :
364 : end do
365 : end do
366 :
367 0 : do k=1,kmp1
368 0 : do i=1,ncol
369 : !
370 : ! edge-level pressure arrays: copy from the arrays computed by dynpkg
371 : !
372 0 : ic = lons(i)
373 0 : jc = lats(i)
374 0 : phys_state(lchnk)%pint (i,k) = pexy (ic,k,jc)
375 0 : phys_state(lchnk)%lnpint(i,k) = pelnxy(ic,k,jc)
376 : end do
377 : end do
378 :
379 : !
380 : ! Copy constituents
381 : ! Dry types converted from moist to dry m.r. at bottom of this routine
382 : !
383 0 : do m=1,pcnst
384 0 : do k=1,km
385 0 : do i=1,ncol
386 0 : phys_state(lchnk)%q(i,k,m) = &
387 0 : tracer(lons(i),lats(i),k,m)
388 : end do
389 : end do
390 : end do
391 :
392 : end do chnk_loop1
393 :
394 : else has_local_map
395 :
396 16128 : boff = 6
397 16128 : if (use_gw_front .or. use_gw_front_igw) boff = boff+2
398 : if (qbo_use_forcing) boff = boff+1
399 :
400 16128 : tsize = boff + 1 + pcnst
401 :
402 16128 : blksiz = (jlastxy-jfirstxy+1)*(ilastxy-ifirstxy+1)
403 64512 : allocate( bpter(blksiz,0:km),stat=astat )
404 16128 : if( astat /= 0 ) then
405 0 : write(iulog,*) 'd_p_coupling: failed to allocate bpter; error = ',astat
406 0 : call endrun
407 : end if
408 48384 : allocate( bbuffer(tsize*block_buf_nrecs),stat=astat )
409 16128 : if( astat /= 0 ) then
410 0 : write(iulog,*) 'd_p_coupling: failed to allocate bbuffer; error = ',astat
411 0 : call endrun
412 : end if
413 48384 : allocate( cbuffer(tsize*chunk_buf_nrecs),stat=astat )
414 16128 : if( astat /= 0 ) then
415 0 : write(iulog,*) 'd_p_coupling: failed to allocate cbuffer; error = ',astat
416 0 : call endrun
417 : end if
418 :
419 16128 : if (iam .lt. grid%npes_xy) then
420 16128 : call block_to_chunk_send_pters(iam+1,blksiz,kmp1,tsize,bpter)
421 : endif
422 :
423 : !$omp parallel do private (j, i, ib, k, m)
424 64512 : do j=jfirstxy,jlastxy
425 1225728 : do i=ifirstxy,ilastxy
426 1161216 : ib = (j-jfirstxy)*(ilastxy-ifirstxy+1) + (i-ifirstxy+1)
427 :
428 285659136 : bbuffer(bpter(ib,0)+4:bpter(ib,0)+boff+pcnst) = 0.0_r8
429 :
430 1161216 : bbuffer(bpter(ib,0)) = pexy(i,kmp1,j)
431 1161216 : bbuffer(bpter(ib,0)+1) = pelnxy(i,kmp1,j)
432 1161216 : bbuffer(bpter(ib,0)+2) = psxy(i,j)
433 1161216 : bbuffer(bpter(ib,0)+3) = phisxy(i,j)
434 :
435 38368512 : do k=1,km
436 :
437 37158912 : bbuffer(bpter(ib,k)) = pexy(i,k,j)
438 37158912 : bbuffer(bpter(ib,k)+1) = pelnxy(i,k,j)
439 37158912 : bbuffer(bpter(ib,k)+2) = u3 (i,k,j)
440 37158912 : bbuffer(bpter(ib,k)+3) = v3 (i,k,j)
441 37158912 : bbuffer(bpter(ib,k)+4) = omgaxy(i,k,j)
442 37158912 : bbuffer(bpter(ib,k)+5) = ptxy(i,j,k) / (D1_0 + zvir*tracer(i,j,k,1))
443 37158912 : bbuffer(bpter(ib,k)+6) = pkxy(i,j,pver+1) / pkzxy(i,j,k)
444 :
445 37158912 : if (use_gw_front .or. use_gw_front_igw) then
446 37158912 : bbuffer(bpter(ib,k)+7) = frontgf(i,k,j)
447 37158912 : bbuffer(bpter(ib,k)+8) = frontga(i,k,j)
448 : end if
449 :
450 : if (qbo_use_forcing) then
451 : bbuffer(bpter(ib,k)+9) = uzm(k,j)
452 : end if
453 :
454 8956459008 : do m=1,pcnst
455 8955297792 : bbuffer(bpter(ib,k)+boff+m) = tracer(i,j,k,m)
456 : end do
457 :
458 : end do
459 : end do
460 : end do
461 :
462 16128 : call t_barrierf('sync_blk_to_chk', grid%commxy)
463 16128 : call t_startf ('block_to_chunk')
464 16128 : call transpose_block_to_chunk(tsize, bbuffer, cbuffer)
465 16128 : call t_stopf ('block_to_chunk')
466 :
467 : !$omp parallel do private (lchnk, ncol, i, k, m, cpter, pbuf_chnk, pbuf_uzm, pbuf_frontgf, pbuf_frontga)
468 : chnk_loop2 : &
469 96768 : do lchnk = begchunk,endchunk
470 80640 : ncol = phys_state(lchnk)%ncol
471 :
472 80640 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
473 :
474 80640 : if (use_gw_front .or. use_gw_front_igw) then
475 80640 : call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf)
476 80640 : call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga)
477 : end if
478 :
479 : if (qbo_use_forcing) then
480 : call pbuf_get_field(pbuf_chnk, uzm_idx, pbuf_uzm)
481 : end if
482 :
483 80640 : call block_to_chunk_recv_pters(lchnk,pcols,pver+1,tsize,cpter)
484 :
485 1257984 : do i=1,ncol
486 :
487 1161216 : phys_state(lchnk)%pint (i,pver+1) = cbuffer(cpter(i,0))
488 1161216 : phys_state(lchnk)%lnpint(i,pver+1) = cbuffer(cpter(i,0)+1)
489 1161216 : phys_state(lchnk)%ps(i) = cbuffer(cpter(i,0)+2)
490 1161216 : phys_state(lchnk)%phis(i) = cbuffer(cpter(i,0)+3)
491 :
492 38400768 : do k=1,km
493 :
494 37158912 : phys_state(lchnk)%pint (i,k) = cbuffer(cpter(i,k))
495 37158912 : phys_state(lchnk)%lnpint(i,k) = cbuffer(cpter(i,k)+1)
496 37158912 : phys_state(lchnk)%u (i,k) = cbuffer(cpter(i,k)+2)
497 37158912 : phys_state(lchnk)%v (i,k) = cbuffer(cpter(i,k)+3)
498 37158912 : phys_state(lchnk)%omega (i,k) = cbuffer(cpter(i,k)+4)
499 37158912 : phys_state(lchnk)%t (i,k) = cbuffer(cpter(i,k)+5)
500 37158912 : phys_state(lchnk)%exner (i,k) = cbuffer(cpter(i,k)+6)
501 :
502 37158912 : if (use_gw_front .or. use_gw_front_igw) then
503 37158912 : pbuf_frontgf(i,k) = cbuffer(cpter(i,k)+7)
504 37158912 : pbuf_frontga(i,k) = cbuffer(cpter(i,k)+8)
505 : end if
506 :
507 : if (qbo_use_forcing) then
508 : pbuf_uzm(i,k) = cbuffer(cpter(i,k)+9)
509 : end if
510 :
511 : ! dry type constituents converted from moist to dry at bottom of routine
512 8956459008 : do m=1,pcnst
513 8955297792 : phys_state(lchnk)%q(i,k,m) = cbuffer(cpter(i,k)+boff+m)
514 : end do
515 :
516 : end do
517 : end do
518 :
519 : end do chnk_loop2
520 :
521 16128 : deallocate(bpter)
522 16128 : deallocate(bbuffer)
523 16128 : deallocate(cbuffer)
524 :
525 : endif has_local_map
526 :
527 : !
528 : ! Evaluate derived quantities
529 : !
530 16128 : call t_startf ('derived_fields')
531 : !$omp parallel do private (lchnk, ncol, i, k, m, qmavl, dqreq, qbot, qbotm1, zvirv, pbuf_chnk)
532 96768 : do lchnk = begchunk,endchunk
533 80640 : ncol = phys_state(lchnk)%ncol
534 2661120 : do k=1,km
535 39820032 : do i=1,ncol
536 37158912 : phys_state(lchnk)%pdel (i,k) = phys_state(lchnk)%pint(i,k+1) - phys_state(lchnk)%pint(i,k)
537 37158912 : phys_state(lchnk)%rpdel(i,k) = D1_0/phys_state(lchnk)%pdel(i,k)
538 37158912 : phys_state(lchnk)%pmid (i,k) = D0_5*(phys_state(lchnk)%pint(i,k) + phys_state(lchnk)%pint(i,k+1))
539 39739392 : phys_state(lchnk)%lnpmid(i,k) = log(phys_state(lchnk)%pmid(i,k))
540 : end do
541 : end do
542 :
543 : ! Attempt to remove negative constituents in bottom layer only by moving from next level
544 : ! This is a BAB kludge to avoid masses of warning messages for cloud water and ice, since
545 : ! the vertical remapping operator currently being used for cam is not strictly monotonic
546 : ! at the endpoints.
547 19434240 : do m=1,pcnst
548 298126080 : do i=1,ncol
549 298045440 : if (phys_state(lchnk)%q(i,pver,m) < qmin(m)) then
550 : ! available q in 2nd level
551 5767995 : qmavl = phys_state(lchnk)%q (i,pver-1,m) - qmin(m)
552 : ! required q change in bottom level rescaled to mass fraction in 2nd level
553 : dqreq = (qmin(m) - phys_state(lchnk)%q(i,pver,m)) &
554 5767995 : * phys_state(lchnk)%pdel(i,pver) / phys_state(lchnk)%pdel(i,pver-1)
555 5767995 : qbot = phys_state(lchnk)%q(i,pver ,m)
556 5767995 : qbotm1 = phys_state(lchnk)%q(i,pver-1,m)
557 5767995 : if (dqreq < qmavl) then
558 199057 : phys_state(lchnk)%q(i,pver ,m) = qmin(m)
559 199057 : phys_state(lchnk)%q(i,pver-1,m) = phys_state(lchnk)%q(i,pver-1,m) - dqreq
560 : ! Comment out these log messages since they can make the log files so
561 : ! large that they're unusable.
562 199057 : if (dqreq>qmin(m) .and. dqreq>fraction*qbotm1 .and. (trim(fv_print_dpcoup_warn) == "full")) &
563 0 : write(iulog,*) 'dpcoup dqreq', m, lchnk, i, qbot, qbotm1, dqreq
564 : else
565 : ! Comment out these log messages since they can make the log files so
566 : ! large that they're unusable.
567 5568938 : if (dqreq>qmin(m) .and. (trim(fv_print_dpcoup_warn) == "full")) then
568 0 : write(iulog,*) 'dpcoup cant adjust', m, lchnk, i, qbot, qbotm1, dqreq
569 : end if
570 : end if
571 : end if
572 : end do
573 : end do
574 : !
575 : ! Convert dry type constituents from moist to dry mixing ratio
576 : ! (note: cam_thermo_dry_air_update assumes dry unless optional conversion factor provided)
577 : !
578 80640 : call set_state_pdry(phys_state(lchnk)) ! First get dry pressure to use for this timestep
579 80640 : call set_wet_to_dry(phys_state(lchnk), convert_cnst_type='dry') ! Dynamics had moist, physics wants dry
580 80640 : if (dry_air_species_num>0) then
581 : !------------------------------------------------------------
582 : ! Apply limiters to mixing ratios of major species
583 : !------------------------------------------------------------
584 0 : call physics_cnst_limit( phys_state(lchnk) )
585 : !-----------------------------------------------------------------------------
586 : ! Call cam_thermo_update to compute cpairv, rairv, mbarv, and cappav as constituent dependent variables
587 : ! and compute molecular viscosity(kmvis) and conductivity(kmcnd)
588 : !-----------------------------------------------------------------------------
589 0 : call cam_thermo_dry_air_update(phys_state(lchnk)%q, phys_state(lchnk)%t, lchnk, ncol)
590 : endif
591 :
592 : !------------------------------------------------------------------------
593 : ! Fill local zvirv variable; calculated for WACCM-X
594 : !------------------------------------------------------------------------
595 80640 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
596 0 : zvirv(:,:) = shr_const_rwv / rairv(:,:,lchnk) -1._r8
597 : else
598 43948800 : zvirv(:,:) = zvir
599 : endif
600 : !
601 : ! Compute initial geopotential heights
602 80640 : call geopotential_t (phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint , &
603 : phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , &
604 : phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv, &
605 161280 : phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol )
606 :
607 : ! Compute initial dry static energy, include surface geopotential
608 2661120 : do k = 1, pver
609 39820032 : do i=1,ncol
610 37158912 : phys_state(lchnk)%s(i,k) = cpairv(i,k,lchnk)*phys_state(lchnk)%t(i,k) &
611 76898304 : + gravit*phys_state(lchnk)%zm(i,k) + phys_state(lchnk)%phis(i)
612 : end do
613 : end do
614 :
615 : !
616 : ! Ensure tracers are all positive
617 : !
618 : call qneg3('D_P_COUPLING',lchnk ,ncol ,pcols ,pver , &
619 80640 : 1, pcnst, qmin ,phys_state(lchnk)%q)
620 :
621 : ! Compute energy and water integrals of input state
622 :
623 80640 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
624 258048 : call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk)
625 :
626 : end do
627 16128 : call t_stopf('derived_fields')
628 :
629 16128 : deallocate (u3)
630 16128 : deallocate (v3)
631 16128 : if (dyn_state%am_diag) then
632 0 : deallocate (du3)
633 0 : deallocate (dv3)
634 0 : deallocate (dua3)
635 0 : deallocate (dva3)
636 0 : deallocate (duf3)
637 0 : deallocate (dummy)
638 : end if
639 :
640 32256 : end subroutine d_p_coupling
641 : !-----------------------------------------------------------------------
642 :
643 : !-----------------------------------------------------------------------
644 : !BOP
645 : ! !IROUTINE: p_d_coupling --- convert physics output to dynamics input
646 : !
647 : ! !INTERFACE:
648 14592 : subroutine p_d_coupling(grid, phys_state, phys_tend, &
649 : dyn_in, dtime, zvir, cappa, ptop)
650 :
651 : ! !USES:
652 : #if ( defined OFFLINE_DYN )
653 : use metdata, only: get_met_fields
654 : #endif
655 16128 : use physics_buffer, only: physics_buffer_desc
656 : use cam_thermo, only: cam_thermo_calc_kappav
657 :
658 : !-----------------------------------------------------------------------
659 : implicit none
660 :
661 : ! Variables ending in xy are xy-decomposition instanciations.
662 :
663 : type(T_FVDYCORE_GRID), intent(in) :: grid ! FV Dynamics grid
664 :
665 : ! !INPUT PARAMETERS:
666 : type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
667 : type(physics_tend), intent(inout), dimension(begchunk:endchunk) :: phys_tend
668 : type(dyn_import_t), intent(inout) :: dyn_in
669 :
670 : real(r8), intent(in) :: dtime
671 : real(r8), intent(in) :: zvir
672 : real(r8), intent(in) :: cappa
673 : real(r8), intent(in) :: ptop
674 :
675 : ! !DESCRIPTION:
676 : !
677 : ! Coupler for converting physics output variables into dynamics input variables
678 : !
679 : ! !REVISION HISTORY:
680 : ! 00.06.01 Boville Creation
681 : ! 01.06.08 AAM Compactified
682 : ! 01.07.13 AAM Some support for multi-2D decompositions
683 : ! 02.03.01 Worley Support for nontrivial physics remapping
684 : ! 02.08.06 Sawyer T3 added -- updated to current temperature
685 : ! 05.07.12 Sawyer Added dyn_state as argument
686 : ! 05.09.23 Sawyer Transitioned to XY decomposition vars. only
687 : ! 05.10.31 Sawyer Replaced dyn_state with dyn_interface
688 : !
689 : !EOP
690 : !-----------------------------------------------------------------------
691 : !BOC
692 : ! !LOCAL VARIABLES:
693 :
694 : type(t_fvdycore_state), pointer :: dyn_state
695 :
696 : ! Variables from the dynamics import container
697 :
698 14592 : real(r8), pointer :: psxy(:,:)
699 14592 : real(r8), pointer :: u3sxy(:,:,:)
700 14592 : real(r8), pointer :: v3sxy(:,:,:)
701 14592 : real(r8), pointer :: t3xy(:,:,:) ! Temperature
702 14592 : real(r8), pointer :: ptxy(:,:,:) ! Virt. pot. temp.
703 14592 : real(r8), pointer :: tracer(:,:,:,:) ! Constituents
704 :
705 14592 : real(r8), pointer :: pexy(:,:,:)
706 14592 : real(r8), pointer :: delpxy(:,:,:)
707 14592 : real(r8), pointer :: pkxy(:,:,:)
708 14592 : real(r8), pointer :: pkzxy(:,:,:)
709 :
710 : ! Local workspace
711 :
712 : real(r8):: dudtxy(grid%ifirstxy:grid%ilastxy,&
713 29184 : grid%km,grid%jfirstxy:grid%jlastxy)
714 : real(r8):: dvdtxy(grid%ifirstxy:grid%ilastxy,&
715 29184 : grid%km,grid%jfirstxy:grid%jlastxy)
716 : real(r8):: dummy_pelnxy(grid%ifirstxy:grid%ilastxy,grid%km+1, &
717 29184 : grid%jfirstxy:grid%jlastxy)
718 :
719 : integer :: i, ib, k, m, j, lchnk ! indices
720 : integer :: ncol ! number of columns in current chunk
721 : integer :: lats(pcols) ! array of latitude indices
722 : integer :: lons(pcols) ! array of longitude indices
723 : integer :: blksiz ! number of columns in 2D block
724 : integer :: tsize ! amount of data per grid point passed to physics
725 14592 : integer, allocatable, dimension(:,:) :: bpter
726 : ! offsets into block buffer for unpacking data
727 : integer :: cpter(pcols,0:pver) ! offsets into chunk buffer for packing data
728 :
729 : real(r8) :: dt5
730 : real(r8), allocatable, dimension(:) :: &
731 14592 : bbuffer, cbuffer ! transpose buffers
732 : #if (! defined SPMD)
733 : integer :: block_buf_nrecs = 0
734 : integer :: chunk_buf_nrecs = 0
735 : logical :: local_dp_map=.true.
736 : #endif
737 : integer :: km, iam
738 : integer :: ifirstxy, ilastxy, jfirstxy, jlastxy
739 :
740 : real(r8) :: cappa3v( grid%ifirstxy:grid%ilastxy,&
741 29184 : grid%jfirstxy:grid%jlastxy, grid%km )
742 :
743 29184 : dyn_state => get_dyn_state()
744 :
745 : ! Pull the variables out of the dynamics export container
746 :
747 14592 : psxy => dyn_in%ps
748 14592 : u3sxy => dyn_in%u3s
749 14592 : v3sxy => dyn_in%v3s
750 14592 : t3xy => dyn_in%t3
751 14592 : ptxy => dyn_in%pt
752 14592 : tracer => dyn_in%tracer
753 :
754 14592 : pexy => dyn_in%pe
755 14592 : delpxy => dyn_in%delp
756 14592 : pkxy => dyn_in%pk
757 14592 : pkzxy => dyn_in%pkz
758 :
759 14592 : km = grid%km
760 :
761 14592 : ifirstxy = grid%ifirstxy
762 14592 : ilastxy = grid%ilastxy
763 14592 : jfirstxy = grid%jfirstxy
764 14592 : jlastxy = grid%jlastxy
765 :
766 14592 : iam = grid%iam
767 :
768 : !---------------------------End Local workspace-------------------------
769 :
770 : #if ( defined OFFLINE_DYN )
771 : !
772 : ! set the dyn flds to offline meteorological data
773 : !
774 : call get_met_fields( phys_state, phys_tend, dtime )
775 : #endif
776 : ! -------------------------------------------------------------------------
777 : ! Copy temperature, tendencies and constituents to dynamics data structures
778 : ! -------------------------------------------------------------------------
779 :
780 : ! -------------------------------------------------------------------------
781 : ! Copy onto xy decomposition, then transpose to yz decomposition
782 : ! -------------------------------------------------------------------------
783 :
784 14592 : if (local_dp_map) then
785 :
786 : !$omp parallel do private(lchnk, i, k, ncol, m, lons, lats)
787 :
788 0 : do lchnk = begchunk,endchunk
789 0 : ncol = get_ncols_p(lchnk)
790 0 : call get_lon_all_p(lchnk, ncol, lons)
791 0 : call get_lat_all_p(lchnk, ncol, lats)
792 :
793 0 : do k = 1, km
794 0 : do i = 1, ncol
795 0 : dvdtxy(lons(i),k,lats(i)) = phys_tend(lchnk)%dvdt(i,k)
796 0 : dudtxy(lons(i),k,lats(i)) = phys_tend(lchnk)%dudt(i,k)
797 0 : ptxy (lons(i),lats(i),k) = phys_state(lchnk)%t(i,k)
798 0 : delpxy(lons(i),lats(i),k) = phys_state(lchnk)%pdel(i,k)
799 : enddo
800 : enddo
801 :
802 0 : do m=1,pcnst
803 0 : do k=1,km
804 0 : do i=1,ncol
805 0 : tracer(lons(i),lats(i),k,m) = &
806 0 : phys_state(lchnk)%q(i,k,m)
807 : end do
808 : end do
809 : end do
810 :
811 : enddo
812 :
813 : else
814 :
815 14592 : tsize = 4 + pcnst
816 :
817 14592 : blksiz = (jlastxy-jfirstxy+1)*(ilastxy-ifirstxy+1)
818 58368 : allocate(bpter(blksiz,0:km))
819 43776 : allocate(bbuffer(tsize*block_buf_nrecs))
820 43776 : allocate(cbuffer(tsize*chunk_buf_nrecs))
821 :
822 : !$omp parallel do private (lchnk, ncol, i, k, m, cpter)
823 87552 : do lchnk = begchunk,endchunk
824 72960 : ncol = get_ncols_p(lchnk)
825 :
826 72960 : call chunk_to_block_send_pters(lchnk,pcols,km+1,tsize,cpter)
827 :
828 1123584 : do i=1,ncol
829 257475840 : cbuffer(cpter(i,0):cpter(i,0)+3+pcnst) = 0.0_r8
830 : end do
831 :
832 2422272 : do k=1,km
833 36027648 : do i=1,ncol
834 :
835 33619968 : cbuffer(cpter(i,k)) = phys_tend(lchnk)%dvdt(i,k)
836 33619968 : cbuffer(cpter(i,k)+1) = phys_tend(lchnk)%dudt(i,k)
837 33619968 : cbuffer(cpter(i,k)+2) = phys_state(lchnk)%t(i,k)
838 33619968 : cbuffer(cpter(i,k)+3) = phys_state(lchnk)%pdel(i,k)
839 :
840 8104747008 : do m=1,pcnst
841 8102412288 : cbuffer(cpter(i,k)+3+m) = phys_state(lchnk)%q(i,k,m)
842 : end do
843 :
844 : end do
845 :
846 : end do
847 :
848 : end do
849 :
850 14592 : call t_barrierf('sync_chk_to_blk', grid%commxy)
851 14592 : call t_startf ('chunk_to_block')
852 14592 : call transpose_chunk_to_block(tsize, cbuffer, bbuffer)
853 14592 : call t_stopf ('chunk_to_block')
854 :
855 14592 : if (iam .lt. grid%npes_xy) then
856 14592 : call chunk_to_block_recv_pters(iam+1,blksiz,km+1,tsize,bpter)
857 : endif
858 :
859 : !$omp parallel do private (j, i, ib, k, m)
860 58368 : do j=jfirstxy,jlastxy
861 1459200 : do k=1,km
862 35064576 : do i=ifirstxy,ilastxy
863 33619968 : ib = (j-jfirstxy)*(ilastxy-ifirstxy+1) + (i-ifirstxy+1)
864 :
865 33619968 : dvdtxy(i,k,j) = bbuffer(bpter(ib,k))
866 33619968 : dudtxy(i,k,j) = bbuffer(bpter(ib,k)+1)
867 33619968 : ptxy (i,j,k) = bbuffer(bpter(ib,k)+2)
868 33619968 : delpxy(i,j,k) = bbuffer(bpter(ib,k)+3)
869 :
870 8103813120 : do m=1,pcnst
871 8102412288 : tracer(i,j,k,m) = bbuffer(bpter(ib,k)+3+m)
872 : end do
873 :
874 : enddo
875 : enddo
876 : enddo
877 :
878 14592 : deallocate(bpter)
879 14592 : deallocate(bbuffer)
880 14592 : deallocate(cbuffer)
881 :
882 : endif
883 :
884 : ! WS: 02.08.06: Update t3 to temperature
885 : !$omp parallel do private(i,j,k)
886 481536 : do k=1,km
887 1882368 : do j = jfirstxy,jlastxy
888 35487744 : do i = ifirstxy,ilastxy
889 35020800 : t3xy(i,j,k) = ptxy(i,j,k)
890 : enddo
891 : enddo
892 : enddo
893 :
894 : ! -------------------------------------------------------------------------
895 : ! Update u3s and v3s from tendencies dudt and dvdt.
896 : ! -------------------------------------------------------------------------
897 14592 : dt5 = D0_5*dtime
898 :
899 14592 : call t_barrierf('sync_uv3s_update', grid%commxy)
900 14592 : call t_startf('uv3s_update')
901 14592 : if (iam .lt. grid%npes_xy) then
902 : call uv3s_update(grid, dudtxy, u3sxy, dvdtxy, v3sxy, dt5, &
903 14592 : dyn_state%am_geom_crrct)
904 : end if ! (iam .lt. grid%npes_xy)
905 14592 : call t_stopf('uv3s_update')
906 :
907 : ! -------------------------------------------------------------------------
908 : ! Compute pt, q3, pe, delp, ps, peln, pkz and pk.
909 : ! For 2-D decomposition, delp is transposed to delpxy, pexy is computed
910 : ! from delpxy (and ptop), and pexy is transposed back to pe.
911 : ! Note that pt, q3, delp and pe are input parameters as well.
912 : ! -------------------------------------------------------------------------
913 14592 : call t_barrierf('sync_p_d_adjust', grid%commxy)
914 14592 : call t_startf ('p_d_adjust')
915 14592 : if (iam .lt. grid%npes_xy) then
916 14592 : if (grid%high_alt) then
917 0 : call cam_thermo_calc_kappav(tracer, cappa3v )
918 : else
919 35502336 : cappa3v = cappa
920 : endif
921 : call p_d_adjust(grid, tracer, dummy_pelnxy, pkxy, pkzxy, zvir, cappa3v, &
922 14592 : delpxy, ptxy, pexy, psxy, ptop)
923 : end if ! (iam .lt. grid%npes_xy)
924 14592 : call t_stopf ('p_d_adjust')
925 :
926 : !EOC
927 29184 : end subroutine p_d_coupling
928 : !-----------------------------------------------------------------------
929 : end module dp_coupling
|
