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