Line data Source code
1 : module dp_coupling
2 :
3 : !-------------------------------------------------------------------------------
4 : ! dynamics - physics coupling module
5 : !-------------------------------------------------------------------------------
6 :
7 : use shr_kind_mod, only: r8=>shr_kind_r8
8 : use ppgrid, only: begchunk, endchunk, pcols, pver, pverp
9 : use constituents, only: pcnst, cnst_type
10 :
11 : use spmd_dyn, only: local_dp_map
12 : use spmd_utils, only: iam
13 : use dyn_grid, only: TimeLevel, edgebuf
14 : use dyn_comp, only: dyn_export_t, dyn_import_t
15 :
16 : use physics_types, only: physics_state, physics_tend, physics_cnst_limit
17 : use phys_grid, only: get_ncols_p
18 : use phys_grid, only: get_dyn_col_p, columns_on_task, get_chunk_info_p, phys_columns_on_task
19 : use physics_buffer, only: physics_buffer_desc, pbuf_get_chunk, pbuf_get_field
20 :
21 : use dp_mapping, only: nphys_pts
22 :
23 : use perf_mod, only: t_startf, t_stopf
24 : use cam_abortutils, only: endrun
25 :
26 : use parallel_mod, only: par
27 : use thread_mod, only: horz_num_threads, max_num_threads
28 : use hybrid_mod, only: config_thread_region, get_loop_ranges, hybrid_t
29 : use dimensions_mod, only: np, nelemd, nlev, qsize, ntrac, fv_nphys
30 :
31 : use dof_mod, only: UniquePoints, PutUniquePoints
32 : use element_mod, only: element_t
33 :
34 : implicit none
35 : private
36 : save
37 :
38 : public :: d_p_coupling, p_d_coupling
39 :
40 : real (kind=r8), allocatable :: q_prev(:,:,:,:) ! Previous Q for computing tendencies
41 :
42 : !=========================================================================================
43 : CONTAINS
44 : !=========================================================================================
45 :
46 32256 : subroutine d_p_coupling(phys_state, phys_tend, pbuf2d, dyn_out)
47 :
48 : ! Convert the dynamics output state into the physics input state.
49 : ! Note that all pressures and tracer mixing ratios coming from the dycore are based on
50 : ! dry air mass.
51 :
52 : use gravity_waves_sources, only: gws_src_fnct,gws_src_vort
53 : use dyn_comp, only: frontgf_idx, frontga_idx, vort4gw_idx
54 : use phys_control, only: use_gw_front, use_gw_front_igw, use_gw_movmtn_pbl
55 : use hycoef, only: hyai, ps0
56 : use fvm_mapping, only: dyn2phys_vector, dyn2phys_all_vars
57 : use se_dyn_time_mod, only: timelevel_qdp
58 : use control_mod, only: qsplit
59 : use test_fvm_mapping, only: test_mapping_overwrite_dyn_state, test_mapping_output_phys_state
60 : use prim_advance_mod, only: tot_energy_dyn
61 : ! arguments
62 : type(dyn_export_t), intent(inout) :: dyn_out ! dynamics export
63 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
64 : type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
65 : type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
66 :
67 :
68 : ! LOCAL VARIABLES
69 10752 : type(element_t), pointer :: elem(:) ! pointer to dyn_out element array
70 : integer :: ie ! indices over elements
71 : integer :: lchnk, icol, ilyr ! indices over chunks, columns, layers
72 :
73 10752 : real (kind=r8), allocatable :: ps_tmp(:,:) ! temp array to hold ps
74 10752 : real (kind=r8), allocatable :: dp3d_tmp(:,:,:) ! temp array to hold dp3d
75 10752 : real (kind=r8), allocatable :: dp3d_tmp_tmp(:,:)
76 10752 : real (kind=r8), allocatable :: phis_tmp(:,:) ! temp array to hold phis
77 10752 : real (kind=r8), allocatable :: T_tmp(:,:,:) ! temp array to hold T
78 10752 : real (kind=r8), allocatable :: uv_tmp(:,:,:,:) ! temp array to hold u and v
79 10752 : real (kind=r8), allocatable :: q_tmp(:,:,:,:) ! temp to hold advected constituents
80 10752 : real (kind=r8), allocatable :: omega_tmp(:,:,:) ! temp array to hold omega
81 :
82 : ! Frontogenesis
83 10752 : real (kind=r8), allocatable :: frontgf(:,:,:) ! temp arrays to hold frontogenesis
84 10752 : real (kind=r8), allocatable :: frontga(:,:,:) ! function (frontgf) and angle (frontga)
85 10752 : real (kind=r8), allocatable :: frontgf_phys(:,:,:)
86 10752 : real (kind=r8), allocatable :: frontga_phys(:,:,:)
87 :
88 : ! Vorticity
89 10752 : real (kind=r8), allocatable :: vort4gw(:,:,:) ! temp arrays to hold vorticity
90 10752 : real (kind=r8), allocatable :: vort4gw_phys(:,:,:)
91 :
92 :
93 : ! Pointers to pbuf
94 10752 : real (kind=r8), pointer :: pbuf_frontgf(:,:)
95 10752 : real (kind=r8), pointer :: pbuf_frontga(:,:)
96 10752 : real (kind=r8), pointer :: pbuf_vort4gw(:,:)
97 :
98 : integer :: ncols, ierr
99 : integer :: col_ind, blk_ind(1), m
100 : integer :: nphys
101 :
102 10752 : real (kind=r8), allocatable :: qgll(:,:,:,:)
103 : real (kind=r8) :: inv_dp3d(np,np,nlev)
104 : integer :: tl_f, tl_qdp_np0, tl_qdp_np1
105 :
106 10752 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
107 : !----------------------------------------------------------------------------
108 :
109 10752 : if (.not. local_dp_map) then
110 0 : call endrun('d_p_coupling: Weak scaling does not support load balancing')
111 : end if
112 :
113 10752 : elem => dyn_out%elem
114 10752 : tl_f = TimeLevel%n0
115 10752 : call TimeLevel_Qdp(TimeLevel, qsplit, tl_qdp_np0,tl_qdp_np1)
116 :
117 10752 : nullify(pbuf_chnk)
118 10752 : nullify(pbuf_frontgf)
119 10752 : nullify(pbuf_frontga)
120 10752 : nullify(pbuf_vort4gw)
121 :
122 :
123 :
124 10752 : if (fv_nphys > 0) then
125 10752 : nphys = fv_nphys
126 : else
127 0 : allocate(qgll(np,np,nlev,pcnst))
128 0 : nphys = np
129 : end if
130 :
131 : ! Allocate temporary arrays to hold data for physics decomposition
132 10752 : allocate(ps_tmp(nphys_pts,nelemd))
133 10752 : allocate(dp3d_tmp(nphys_pts,pver,nelemd))
134 10752 : allocate(dp3d_tmp_tmp(nphys_pts,pver))
135 10752 : allocate(phis_tmp(nphys_pts,nelemd))
136 10752 : allocate(T_tmp(nphys_pts,pver,nelemd))
137 10752 : allocate(uv_tmp(nphys_pts,2,pver,nelemd))
138 10752 : allocate(q_tmp(nphys_pts,pver,pcnst,nelemd))
139 10752 : allocate(omega_tmp(nphys_pts,pver,nelemd))
140 :
141 10752 : call tot_energy_dyn(elem,dyn_out%fvm, 1, nelemd,tl_f , tl_qdp_np0,'dBF')
142 :
143 10752 : if (use_gw_front .or. use_gw_front_igw) then
144 0 : allocate(frontgf(nphys_pts,pver,nelemd), stat=ierr)
145 0 : if (ierr /= 0) call endrun("dp_coupling: Allocate of frontgf failed.")
146 0 : allocate(frontga(nphys_pts,pver,nelemd), stat=ierr)
147 0 : if (ierr /= 0) call endrun("dp_coupling: Allocate of frontga failed.")
148 : end if
149 10752 : if (use_gw_movmtn_pbl) then
150 0 : allocate(vort4gw(nphys_pts,pver,nelemd), stat=ierr)
151 0 : if (ierr /= 0) call endrun("dp_coupling: Allocate of vort4gw failed.")
152 : end if
153 :
154 10752 : if (iam < par%nprocs) then
155 10752 : if (use_gw_front .or. use_gw_front_igw ) then
156 0 : call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys)
157 : end if
158 10752 : if (use_gw_movmtn_pbl ) then
159 0 : call gws_src_vort(elem, tl_f, tl_qdp_np0, vort4gw, nphys)
160 : end if
161 :
162 10752 : if (fv_nphys > 0) then
163 10752 : call test_mapping_overwrite_dyn_state(elem,dyn_out%fvm)
164 : !******************************************************************
165 : ! physics runs on an FVM grid: map GLL vars to physics grid
166 : !******************************************************************
167 10752 : call t_startf('dyn2phys')
168 : ! note that the fvm halo has been filled in prim_run_subcycle
169 : ! if physics grid resolution is not equal to fvm resolution
170 : call dyn2phys_all_vars(1,nelemd,elem, dyn_out%fvm,&
171 : pcnst,hyai(1)*ps0,tl_f, &
172 : ! output
173 : dp3d_tmp, ps_tmp, q_tmp, T_tmp, &
174 : omega_tmp, phis_tmp &
175 10752 : )
176 124152 : do ie = 1, nelemd
177 0 : uv_tmp(:,:,:,ie) = &
178 124152 : dyn2phys_vector(elem(ie)%state%v(:,:,:,:,tl_f),elem(ie))
179 : end do
180 10752 : call t_stopf('dyn2phys')
181 : else
182 :
183 : !******************************************************************
184 : ! Physics runs on GLL grid: collect unique points before mapping to
185 : ! physics decomposition
186 : !******************************************************************
187 :
188 0 : if (qsize < pcnst) then
189 0 : call endrun('d_p_coupling: Fewer GLL tracers advected than required')
190 : end if
191 :
192 0 : call t_startf('UniquePoints')
193 0 : do ie = 1, nelemd
194 0 : inv_dp3d(:,:,:) = 1.0_r8/elem(ie)%state%dp3d(:,:,:,tl_f)
195 0 : do m=1,pcnst
196 0 : qgll(:,:,:,m) = elem(ie)%state%Qdp(:,:,:,m,tl_qdp_np0)*inv_dp3d(:,:,:)
197 : end do
198 0 : ncols = elem(ie)%idxP%NumUniquePts
199 0 : call UniquePoints(elem(ie)%idxP, elem(ie)%state%psdry(:,:), ps_tmp(1:ncols,ie))
200 0 : call UniquePoints(elem(ie)%idxP, nlev, elem(ie)%state%dp3d(:,:,:,tl_f), dp3d_tmp(1:ncols,:,ie))
201 0 : call UniquePoints(elem(ie)%idxP, nlev, elem(ie)%state%T(:,:,:,tl_f), T_tmp(1:ncols,:,ie))
202 0 : call UniquePoints(elem(ie)%idxV, 2, nlev, elem(ie)%state%V(:,:,:,:,tl_f), uv_tmp(1:ncols,:,:,ie))
203 0 : call UniquePoints(elem(ie)%idxV, nlev, elem(ie)%derived%omega, omega_tmp(1:ncols,:,ie))
204 :
205 0 : call UniquePoints(elem(ie)%idxP, elem(ie)%state%phis, phis_tmp(1:ncols,ie))
206 0 : call UniquePoints(elem(ie)%idxP, nlev, pcnst, qgll,Q_tmp(1:ncols,:,:,ie))
207 : end do
208 0 : call t_stopf('UniquePoints')
209 :
210 : end if ! if fv_nphys>0
211 :
212 : else
213 :
214 0 : ps_tmp(:,:) = 0._r8
215 0 : T_tmp(:,:,:) = 0._r8
216 0 : uv_tmp(:,:,:,:) = 0._r8
217 0 : omega_tmp(:,:,:) = 0._r8
218 0 : phis_tmp(:,:) = 0._r8
219 0 : Q_tmp(:,:,:,:) = 0._r8
220 :
221 0 : if (use_gw_front .or. use_gw_front_igw) then
222 0 : frontgf(:,:,:) = 0._r8
223 0 : frontga(:,:,:) = 0._r8
224 : end if
225 0 : if (use_gw_movmtn_pbl) then
226 0 : vort4gw(:,:,:) = 0._r8
227 : end if
228 :
229 : endif ! iam < par%nprocs
230 :
231 10752 : if (fv_nphys < 1) then
232 0 : deallocate(qgll)
233 : end if
234 :
235 : ! q_prev is for saving the tracer fields for calculating tendencies
236 10752 : if (.not. allocated(q_prev)) then
237 1024 : allocate(q_prev(pcols,pver,pcnst,begchunk:endchunk))
238 : end if
239 93632112 : q_prev = 0.0_R8
240 :
241 10752 : call t_startf('dpcopy')
242 10752 : if (use_gw_front .or. use_gw_front_igw) then
243 0 : allocate(frontgf_phys(pcols, pver, begchunk:endchunk))
244 0 : allocate(frontga_phys(pcols, pver, begchunk:endchunk))
245 : end if
246 10752 : if (use_gw_movmtn_pbl) then
247 0 : allocate(vort4gw_phys(pcols, pver, begchunk:endchunk))
248 : end if
249 : !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m)
250 1031352 : do col_ind = 1, phys_columns_on_task
251 1020600 : call get_dyn_col_p(col_ind, ie, blk_ind)
252 1020600 : call get_chunk_info_p(col_ind, lchnk, icol)
253 1020600 : phys_state(lchnk)%ps(icol) = ps_tmp(blk_ind(1), ie)
254 1020600 : phys_state(lchnk)%phis(icol) = phis_tmp(blk_ind(1), ie)
255 27556200 : do ilyr = 1, pver
256 26535600 : phys_state(lchnk)%pdel(icol, ilyr) = dp3d_tmp(blk_ind(1), ilyr, ie)
257 26535600 : phys_state(lchnk)%t(icol, ilyr) = T_tmp(blk_ind(1), ilyr, ie)
258 26535600 : phys_state(lchnk)%u(icol, ilyr) = uv_tmp(blk_ind(1), 1, ilyr, ie)
259 26535600 : phys_state(lchnk)%v(icol, ilyr) = uv_tmp(blk_ind(1), 2, ilyr, ie)
260 26535600 : phys_state(lchnk)%omega(icol, ilyr) = omega_tmp(blk_ind(1), ilyr, ie)
261 :
262 26535600 : if (use_gw_front .or. use_gw_front_igw) then
263 0 : frontgf_phys(icol, ilyr, lchnk) = frontgf(blk_ind(1), ilyr, ie)
264 0 : frontga_phys(icol, ilyr, lchnk) = frontga(blk_ind(1), ilyr, ie)
265 : end if
266 27556200 : if (use_gw_movmtn_pbl) then
267 0 : vort4gw_phys(icol, ilyr, lchnk) = vort4gw(blk_ind(1), ilyr, ie)
268 : end if
269 : end do
270 :
271 5113752 : do m = 1, pcnst
272 83689200 : do ilyr = 1, pver
273 82668600 : phys_state(lchnk)%q(icol, ilyr,m) = Q_tmp(blk_ind(1), ilyr,m, ie)
274 : end do
275 : end do
276 : end do
277 10752 : if (use_gw_front .or. use_gw_front_igw) then
278 : !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_frontgf, pbuf_frontga)
279 0 : do lchnk = begchunk, endchunk
280 0 : ncols = get_ncols_p(lchnk)
281 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
282 0 : call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf)
283 0 : call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga)
284 0 : do icol = 1, ncols
285 0 : do ilyr = 1, pver
286 0 : pbuf_frontgf(icol, ilyr) = frontgf_phys(icol, ilyr, lchnk)
287 0 : pbuf_frontga(icol, ilyr) = frontga_phys(icol, ilyr, lchnk)
288 : end do
289 : end do
290 : end do
291 0 : deallocate(frontgf_phys)
292 0 : deallocate(frontga_phys)
293 : end if
294 10752 : if (use_gw_movmtn_pbl) then
295 : !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, pbuf_chnk, pbuf_vort4gw)
296 0 : do lchnk = begchunk, endchunk
297 0 : ncols = get_ncols_p(lchnk)
298 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
299 0 : call pbuf_get_field(pbuf_chnk, vort4gw_idx, pbuf_vort4gw)
300 0 : do icol = 1, ncols
301 0 : do ilyr = 1, pver
302 0 : pbuf_vort4gw(icol, ilyr) = vort4gw_phys(icol, ilyr, lchnk)
303 : end do
304 : end do
305 : end do
306 0 : deallocate(vort4gw_phys)
307 : end if
308 :
309 10752 : call t_stopf('dpcopy')
310 :
311 : ! Save the tracer fields input to physics package for calculating tendencies
312 : ! The mixing ratios are all dry at this point.
313 81144 : do lchnk = begchunk, endchunk
314 70392 : ncols = phys_state(lchnk)%ncol
315 85389696 : q_prev(1:ncols,1:pver,1:pcnst,lchnk) = phys_state(lchnk)%q(1:ncols,1:pver,1:pcnst)
316 : end do
317 10752 : call test_mapping_output_phys_state(phys_state,dyn_out%fvm)
318 :
319 : ! Deallocate the temporary arrays
320 10752 : deallocate(ps_tmp)
321 10752 : deallocate(dp3d_tmp)
322 10752 : deallocate(phis_tmp)
323 10752 : deallocate(T_tmp)
324 10752 : deallocate(uv_tmp)
325 10752 : deallocate(q_tmp)
326 10752 : deallocate(omega_tmp)
327 :
328 : ! ps, pdel, and q in phys_state are all dry at this point. After return from derived_phys_dry
329 : ! ps and pdel include water vapor only, and the 'wet' constituents have been converted to wet mmr.
330 10752 : call t_startf('derived_phys')
331 10752 : call derived_phys_dry(phys_state, phys_tend, pbuf2d)
332 10752 : call t_stopf('derived_phys')
333 :
334 : !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, ilyr, icol)
335 81144 : do lchnk = begchunk, endchunk
336 70392 : ncols=get_ncols_p(lchnk)
337 81144 : if (pcols > ncols) then
338 116424 : phys_state(lchnk)%phis(ncols+1:) = 0.0_r8
339 : end if
340 : end do
341 21504 : end subroutine d_p_coupling
342 :
343 : !=========================================================================================
344 :
345 19456 : subroutine p_d_coupling(phys_state, phys_tend, dyn_in, tl_f, tl_qdp)
346 :
347 : ! Convert the physics output state into the dynamics input state.
348 :
349 10752 : use phys_grid, only: get_dyn_col_p, columns_on_task, get_chunk_info_p, phys_columns_on_task
350 : use bndry_mod, only: bndry_exchange
351 : use edge_mod, only: edgeVpack, edgeVunpack
352 : use fvm_mapping, only: phys2dyn_forcings_fvm
353 : use test_fvm_mapping, only: test_mapping_overwrite_tendencies
354 : use test_fvm_mapping, only: test_mapping_output_mapped_tendencies
355 : use dimensions_mod, only: use_cslam
356 : ! arguments
357 : type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
358 : type(physics_tend), intent(inout), dimension(begchunk:endchunk) :: phys_tend
359 : integer, intent(in) :: tl_qdp, tl_f
360 : type(dyn_import_t), intent(inout) :: dyn_in
361 : type(hybrid_t) :: hybrid
362 :
363 : ! LOCAL VARIABLES
364 : integer :: ncols ! index
365 9728 : type(element_t), pointer :: elem(:) ! pointer to dyn_in element array
366 : integer :: ie ! index for elements
367 : integer :: col_ind ! index over columns
368 : integer :: blk_ind(1) ! element offset
369 : integer :: lchnk, icol, ilyr ! indices for chunk, column, layer
370 :
371 9728 : real (kind=r8), allocatable :: dp_phys(:,:,:) ! temp array to hold dp on physics grid
372 9728 : real (kind=r8), allocatable :: T_tmp(:,:,:) ! temp array to hold T
373 9728 : real (kind=r8), allocatable :: dq_tmp(:,:,:,:) ! temp array to hold q
374 9728 : real (kind=r8), allocatable :: uv_tmp(:,:,:,:) ! temp array to hold uv
375 : integer :: m, i, j, k
376 :
377 : real (kind=r8) :: factor
378 : integer :: num_trac
379 : integer :: nets, nete
380 : integer :: kptr, ii
381 : !----------------------------------------------------------------------------
382 :
383 9728 : if (.not. local_dp_map) then
384 0 : call endrun('p_d_coupling: Weak scaling does not support load balancing')
385 : end if
386 :
387 9728 : if (iam < par%nprocs) then
388 9728 : elem => dyn_in%elem
389 : else
390 0 : nullify(elem)
391 : end if
392 :
393 9728 : allocate(T_tmp(nphys_pts,pver,nelemd))
394 9728 : allocate(uv_tmp(nphys_pts,2,pver,nelemd))
395 9728 : allocate(dq_tmp(nphys_pts,pver,pcnst,nelemd))
396 9728 : allocate(dp_phys(nphys_pts,pver,nelemd))
397 :
398 26788328 : T_tmp = 0.0_r8
399 56131928 : uv_tmp = 0.0_r8
400 80448128 : dq_tmp = 0.0_r8
401 :
402 9728 : if (.not. allocated(q_prev)) then
403 0 : call endrun('p_d_coupling: q_prev not allocated')
404 : end if
405 :
406 : ! Convert wet to dry mixing ratios and modify the physics temperature
407 : ! tendency to be thermodynamically consistent with the dycore.
408 : !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, icol, ilyr, m, factor)
409 73416 : do lchnk = begchunk, endchunk
410 63688 : ncols = get_ncols_p(lchnk)
411 996816 : do icol = 1, ncols
412 24995488 : do ilyr = 1, pver
413 : ! convert wet mixing ratios to dry
414 24008400 : factor = phys_state(lchnk)%pdel(icol,ilyr)/phys_state(lchnk)%pdeldry(icol,ilyr)
415 96957000 : do m = 1, pcnst
416 96033600 : if (cnst_type(m) == 'wet') then
417 72025200 : phys_state(lchnk)%q(icol,ilyr,m) = factor*phys_state(lchnk)%q(icol,ilyr,m)
418 : end if
419 : end do
420 : end do
421 : end do
422 : end do
423 :
424 9728 : call t_startf('pd_copy')
425 : !$omp parallel do num_threads(max_num_threads) private (col_ind, lchnk, icol, ie, blk_ind, ilyr, m)
426 933128 : do col_ind = 1, phys_columns_on_task
427 923400 : call get_dyn_col_p(col_ind, ie, blk_ind)
428 923400 : call get_chunk_info_p(col_ind, lchnk, icol)
429 :
430 : ! test code -- does nothing unless cpp macro debug_coupling is defined.
431 923400 : call test_mapping_overwrite_tendencies(phys_state(lchnk), &
432 923400 : phys_tend(lchnk), ncols, lchnk, q_prev(1:ncols,:,:,lchnk), &
433 2770200 : dyn_in%fvm)
434 :
435 25864928 : do ilyr = 1, pver
436 24008400 : dp_phys(blk_ind(1),ilyr,ie) = phys_state(lchnk)%pdeldry(icol,ilyr)
437 24008400 : T_tmp(blk_ind(1),ilyr,ie) = phys_tend(lchnk)%dtdt(icol,ilyr)
438 24008400 : uv_tmp(blk_ind(1),1,ilyr,ie) = phys_tend(lchnk)%dudt(icol,ilyr)
439 24008400 : uv_tmp(blk_ind(1),2,ilyr,ie) = phys_tend(lchnk)%dvdt(icol,ilyr)
440 96957000 : do m = 1, pcnst
441 0 : dq_tmp(blk_ind(1),ilyr,m,ie) = &
442 96033600 : (phys_state(lchnk)%q(icol,ilyr,m) - q_prev(icol,ilyr,m,lchnk))
443 : end do
444 : end do
445 : end do
446 9728 : call t_stopf('pd_copy')
447 :
448 9728 : if (iam < par%nprocs) then
449 :
450 9728 : if (fv_nphys > 0) then
451 :
452 : ! put forcings into fvm structure
453 9728 : num_trac = max(qsize,ntrac)
454 112328 : do ie = 1, nelemd
455 420128 : do j = 1, fv_nphys
456 1333800 : do i = 1, fv_nphys
457 923400 : ii = i + (j-1)*fv_nphys
458 24931800 : dyn_in%fvm(ie)%ft(i,j,1:pver) = T_tmp(ii,1:pver,ie)
459 72948600 : dyn_in%fvm(ie)%fm(i,j,1:2,1:pver) = uv_tmp(ii,1:2,1:pver,ie)
460 75718800 : dyn_in%fvm(ie)%fc_phys(i,j,1:pver,1:num_trac) = dq_tmp(ii,1:pver,1:num_trac,ie)
461 25239600 : dyn_in%fvm(ie)%dp_phys(i,j,1:pver) = dp_phys(ii,1:pver,ie)
462 : end do
463 : end do
464 : end do
465 :
466 : !JMD $OMP PARALLEL NUM_THREADS(horz_num_threads), DEFAULT(SHARED), PRIVATE(hybrid,nets,nete,n)
467 : !JMD hybrid = config_thread_region(par,'horizontal')
468 9728 : hybrid = config_thread_region(par,'serial')
469 9728 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
470 : !
471 : ! high-order mapping of ft and fm using fvm technology
472 : !
473 9728 : call t_startf('phys2dyn')
474 9728 : call phys2dyn_forcings_fvm(elem, dyn_in%fvm, hybrid,nets,nete,ntrac==0, tl_f, tl_qdp)
475 9728 : call t_stopf('phys2dyn')
476 : else
477 :
478 0 : call t_startf('putUniquePoints')
479 :
480 : !$omp parallel do num_threads(max_num_threads) private(ie,ncols)
481 0 : do ie = 1, nelemd
482 0 : ncols = elem(ie)%idxP%NumUniquePts
483 0 : call putUniquePoints(elem(ie)%idxP, nlev, T_tmp(1:ncols,:,ie), &
484 0 : elem(ie)%derived%fT(:,:,:))
485 0 : call putUniquePoints(elem(ie)%idxV, 2, nlev, uv_tmp(1:ncols,:,:,ie), &
486 0 : elem(ie)%derived%fM(:,:,:,:))
487 0 : call putUniquePoints(elem(ie)%idxV, nlev,pcnst, dq_tmp(1:ncols,:,:,ie), &
488 0 : elem(ie)%derived%fQ(:,:,:,:))
489 : end do
490 0 : call t_stopf('putUniquePoints')
491 : end if
492 : end if
493 :
494 9728 : deallocate(T_tmp)
495 9728 : deallocate(uv_tmp)
496 9728 : deallocate(dq_tmp)
497 :
498 : ! Boundary exchange for physics forcing terms.
499 : ! For physics on GLL grid, for points with duplicate degrees of freedom,
500 : ! putuniquepoints() set one of the element values and set the others to zero,
501 : ! so do a simple sum (boundary exchange with no weights).
502 : ! For physics grid, we interpolated into all points, so do weighted average.
503 :
504 9728 : call t_startf('p_d_coupling:bndry_exchange')
505 :
506 112328 : do ie = 1, nelemd
507 102600 : if (fv_nphys > 0) then
508 2770200 : do k = 1, nlev
509 2667600 : dyn_in%elem(ie)%derived%FM(:,:,1,k) = &
510 2667600 : dyn_in%elem(ie)%derived%FM(:,:,1,k) * &
511 61354800 : dyn_in%elem(ie)%spheremp(:,:)
512 2667600 : dyn_in%elem(ie)%derived%FM(:,:,2,k) = &
513 2667600 : dyn_in%elem(ie)%derived%FM(:,:,2,k) * &
514 61354800 : dyn_in%elem(ie)%spheremp(:,:)
515 2667600 : dyn_in%elem(ie)%derived%FT(:,:,k) = &
516 2667600 : dyn_in%elem(ie)%derived%FT(:,:,k) * &
517 61457400 : dyn_in%elem(ie)%spheremp(:,:)
518 : end do
519 : end if
520 102600 : kptr = 0
521 102600 : call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
522 102600 : kptr = kptr + 2*nlev
523 102600 : call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
524 112328 : if (.not. use_cslam) then
525 : !
526 : ! if using CSLAM qdp is being overwritten with CSLAM values in the dynamics
527 : ! so no need to do boundary exchange of tracer tendency on GLL grid here
528 : !
529 0 : kptr = kptr + nlev
530 0 : call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
531 : end if
532 : end do
533 :
534 9728 : if (iam < par%nprocs) then
535 9728 : call bndry_exchange(par, edgebuf, location='p_d_coupling')
536 : end if
537 :
538 112328 : do ie = 1, nelemd
539 102600 : kptr = 0
540 102600 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
541 102600 : kptr = kptr + 2*nlev
542 102600 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
543 102600 : kptr = kptr + nlev
544 102600 : if (.not. use_cslam) then
545 0 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
546 : end if
547 112328 : if (fv_nphys > 0) then
548 2770200 : do k = 1, nlev
549 2667600 : dyn_in%elem(ie)%derived%FM(:,:,1,k) = &
550 2667600 : dyn_in%elem(ie)%derived%FM(:,:,1,k) * &
551 61354800 : dyn_in%elem(ie)%rspheremp(:,:)
552 2667600 : dyn_in%elem(ie)%derived%FM(:,:,2,k) = &
553 2667600 : dyn_in%elem(ie)%derived%FM(:,:,2,k) * &
554 61354800 : dyn_in%elem(ie)%rspheremp(:,:)
555 2667600 : dyn_in%elem(ie)%derived%FT(:,:,k) = &
556 2667600 : dyn_in%elem(ie)%derived%FT(:,:,k) * &
557 61457400 : dyn_in%elem(ie)%rspheremp(:,:)
558 : end do
559 : end if
560 : end do
561 9728 : call t_stopf('p_d_coupling:bndry_exchange')
562 :
563 9728 : if (iam < par%nprocs .and. fv_nphys > 0) then
564 0 : call test_mapping_output_mapped_tendencies(dyn_in%fvm(1:nelemd), elem(1:nelemd), &
565 9728 : 1, nelemd, tl_f, tl_qdp)
566 : end if
567 19456 : end subroutine p_d_coupling
568 :
569 : !=========================================================================================
570 :
571 10752 : subroutine derived_phys_dry(phys_state, phys_tend, pbuf2d)
572 :
573 : ! The ps, pdel, and q components of phys_state are all dry on input.
574 : ! On output the psdry and pdeldry components are initialized; ps and pdel are
575 : ! updated to contain contribution from water vapor only; the 'wet' constituent
576 : ! mixing ratios are converted to a wet basis. Initialize geopotential heights.
577 : ! Finally compute energy and water column integrals of the physics input state.
578 :
579 9728 : use constituents, only: qmin
580 : use physconst, only: gravit, zvir
581 : use cam_thermo, only: cam_thermo_dry_air_update, cam_thermo_water_update
582 : use air_composition, only: thermodynamic_active_species_num
583 : use air_composition, only: thermodynamic_active_species_idx
584 : use air_composition, only: cpairv, rairv, cappav, dry_air_species_num
585 : use shr_const_mod, only: shr_const_rwv
586 : use phys_control, only: waccmx_is
587 : use geopotential, only: geopotential_t
588 : use static_energy, only: update_dry_static_energy_run
589 : use check_energy, only: check_energy_timestep_init
590 : use hycoef, only: hyai, ps0
591 : use shr_vmath_mod, only: shr_vmath_log
592 : use qneg_module, only: qneg3
593 : use dyn_tests_utils, only: vc_dry_pressure
594 : use shr_kind_mod, only: shr_kind_cx
595 : ! arguments
596 : type(physics_state), intent(inout), dimension(begchunk:endchunk) :: phys_state
597 : type(physics_tend ), intent(inout), dimension(begchunk:endchunk) :: phys_tend
598 : type(physics_buffer_desc), pointer :: pbuf2d(:,:)
599 :
600 : ! local variables
601 : integer :: lchnk
602 :
603 : real(r8) :: zvirv(pcols,pver) ! Local zvir array pointer
604 : real(r8) :: factor_array(pcols,nlev)
605 :
606 : integer :: m, i, k, ncol, m_cnst
607 10752 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
608 :
609 : !Needed for "update_dry_static_energy" CCPP scheme
610 : integer :: errflg
611 : character(len=shr_kind_cx) :: errmsg
612 : !----------------------------------------------------------------------------
613 :
614 : ! Evaluate derived quantities
615 :
616 : !!$omp parallel do num_threads(horz_num_threads) private (lchnk, ncol, k, i, m , zvirv, pbuf_chnk, factor_array)
617 81144 : do lchnk = begchunk,endchunk
618 :
619 70392 : ncol = get_ncols_p(lchnk)
620 :
621 : ! dry pressure variables
622 :
623 1090992 : do i = 1, ncol
624 27626592 : phys_state(lchnk)%psdry(i) = hyai(1)*ps0 + sum(phys_state(lchnk)%pdel(i,:))
625 : end do
626 1090992 : do i = 1, ncol
627 1090992 : phys_state(lchnk)%pintdry(i,1) = hyai(1)*ps0
628 : end do
629 140784 : call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,1), &
630 211176 : phys_state(lchnk)%lnpintdry(1:ncol,1),ncol)
631 1900584 : do k = 1, nlev
632 28365792 : do i = 1, ncol
633 53071200 : phys_state(lchnk)%pintdry(i,k+1) = phys_state(lchnk)%pintdry(i,k) + &
634 81436992 : phys_state(lchnk)%pdel(i,k)
635 : end do
636 3660384 : call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,k+1),&
637 5560968 : phys_state(lchnk)%lnpintdry(1:ncol,k+1),ncol)
638 : end do
639 :
640 1900584 : do k=1,nlev
641 28365792 : do i=1,ncol
642 26535600 : phys_state(lchnk)%pdeldry (i,k) = phys_state(lchnk)%pdel(i,k)
643 26535600 : phys_state(lchnk)%rpdeldry(i,k) = 1._r8/phys_state(lchnk)%pdeldry(i,k)
644 53071200 : phys_state(lchnk)%pmiddry (i,k) = 0.5D0*(phys_state(lchnk)%pintdry(i,k+1) + &
645 81436992 : phys_state(lchnk)%pintdry(i,k))
646 : end do
647 3660384 : call shr_vmath_log(phys_state(lchnk)%pmiddry(1:ncol,k), &
648 5560968 : phys_state(lchnk)%lnpmiddry(1:ncol,k),ncol)
649 : end do
650 :
651 : ! wet pressure variables (should be removed from physics!)
652 31183656 : factor_array(:,:) = 1.0_r8
653 281568 : do m_cnst=dry_air_species_num + 1,thermodynamic_active_species_num
654 211176 : m = thermodynamic_active_species_idx(m_cnst)
655 5772144 : do k=1,nlev
656 85308552 : do i=1,ncol
657 : ! at this point all q's are dry
658 85097376 : factor_array(i,k) = factor_array(i,k)+phys_state(lchnk)%q(i,k,m)
659 : end do
660 : end do
661 : end do
662 :
663 1900584 : do k=1,nlev
664 28436184 : do i=1,ncol
665 28365792 : phys_state(lchnk)%pdel (i,k) = phys_state(lchnk)%pdeldry(i,k)*factor_array(i,k)
666 : end do
667 : end do
668 :
669 : ! initialize vertical loop - model top pressure
670 :
671 1090992 : do i=1,ncol
672 1020600 : phys_state(lchnk)%ps(i) = phys_state(lchnk)%pintdry(i,1)
673 1090992 : phys_state(lchnk)%pint(i,1) = phys_state(lchnk)%pintdry(i,1)
674 : end do
675 1900584 : do k = 1, nlev
676 28365792 : do i=1,ncol
677 26535600 : phys_state(lchnk)%pint(i,k+1) = phys_state(lchnk)%pint(i,k)+phys_state(lchnk)%pdel(i,k)
678 26535600 : phys_state(lchnk)%pmid(i,k) = (phys_state(lchnk)%pint(i,k+1)+phys_state(lchnk)%pint(i,k))/2._r8
679 28365792 : phys_state(lchnk)%ps (i) = phys_state(lchnk)%ps(i) + phys_state(lchnk)%pdel(i,k)
680 : end do
681 1830192 : call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,k),phys_state(lchnk)%lnpint(1:ncol,k),ncol)
682 1900584 : call shr_vmath_log(phys_state(lchnk)%pmid(1:ncol,k),phys_state(lchnk)%lnpmid(1:ncol,k),ncol)
683 : end do
684 70392 : call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,pverp),phys_state(lchnk)%lnpint(1:ncol,pverp),ncol)
685 :
686 1900584 : do k = 1, nlev
687 28436184 : do i = 1, ncol
688 28365792 : phys_state(lchnk)%rpdel(i,k) = 1._r8/phys_state(lchnk)%pdel(i,k)
689 : end do
690 : end do
691 :
692 : !------------------------------------------------------------
693 : ! Apply limiters to mixing ratios of major species (waccmx)
694 : !------------------------------------------------------------
695 70392 : if (dry_air_species_num>0) then
696 0 : call physics_cnst_limit( phys_state(lchnk) )
697 : !-----------------------------------------------------------------------------
698 : ! Call cam_thermo_dry_air_update to compute cpairv, rairv, mbarv, and cappav as
699 : ! constituent dependent variables.
700 : ! Compute molecular viscosity(kmvis) and conductivity(kmcnd).
701 : ! Fill local zvirv variable; calculated for WACCM-X.
702 : !-----------------------------------------------------------------------------
703 0 : call cam_thermo_dry_air_update(phys_state(lchnk)%q, phys_state(lchnk)%t, lchnk, ncol)
704 0 : zvirv(:,:) = shr_const_rwv / rairv(:,:,lchnk) -1._r8
705 : else
706 31183656 : zvirv(:,:) = zvir
707 : end if
708 : !
709 : ! update cp_dycore in module air_composition.
710 : ! (note: at this point q is dry)
711 : !
712 70392 : call cam_thermo_water_update(phys_state(lchnk)%q(1:ncol,:,:), lchnk, ncol, vc_dry_pressure)
713 1900584 : do k = 1, nlev
714 28436184 : do i = 1, ncol
715 53071200 : phys_state(lchnk)%exner(i,k) = (phys_state(lchnk)%pint(i,pver+1) &
716 81436992 : / phys_state(lchnk)%pmid(i,k))**cappav(i,k,lchnk)
717 : end do
718 : end do
719 : !
720 : ! CAM physics: water tracers are moist; the rest dry
721 : !
722 28436184 : factor_array(1:ncol,1:nlev) = 1._r8/factor_array(1:ncol,1:nlev)
723 281568 : do m = 1,pcnst
724 281568 : if (cnst_type(m) == 'wet') then
725 5701752 : do k = 1, nlev
726 85308552 : do i = 1, ncol
727 85097376 : phys_state(lchnk)%q(i,k,m) = factor_array(i,k)*phys_state(lchnk)%q(i,k,m)
728 : end do
729 : end do
730 : end if
731 : end do
732 :
733 : ! Ensure tracers are all positive
734 : call qneg3('D_P_COUPLING',lchnk ,ncol ,pcols ,pver , &
735 70392 : 1, pcnst, qmin ,phys_state(lchnk)%q)
736 :
737 : ! Compute initial geopotential heights - based on full pressure
738 211176 : call geopotential_t(phys_state(lchnk)%lnpint, phys_state(lchnk)%lnpmid , phys_state(lchnk)%pint, &
739 211176 : phys_state(lchnk)%pmid , phys_state(lchnk)%pdel , phys_state(lchnk)%rpdel , &
740 211176 : phys_state(lchnk)%t , phys_state(lchnk)%q(:,:,:), rairv(:,:,lchnk), gravit, zvirv , &
741 703920 : phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol)
742 : ! Compute initial dry static energy, include surface geopotential
743 140784 : call update_dry_static_energy_run(pver, gravit, phys_state(lchnk)%t(1:ncol,:), &
744 140784 : phys_state(lchnk)%zm(1:ncol,:), &
745 140784 : phys_state(lchnk)%phis(1:ncol), &
746 140784 : phys_state(lchnk)%s(1:ncol,:), &
747 633528 : cpairv(1:ncol,:,lchnk), errflg, errmsg)
748 : ! Compute energy and water integrals of input state
749 70392 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
750 81144 : call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk)
751 :
752 :
753 : end do ! lchnk
754 :
755 21504 : end subroutine derived_phys_dry
756 : end module dp_coupling
|
