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 1112832 : 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 370944 : 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 370944 : real (kind=r8), allocatable :: ps_tmp(:,:) ! temp array to hold ps
74 370944 : real (kind=r8), allocatable :: dp3d_tmp(:,:,:) ! temp array to hold dp3d
75 370944 : real (kind=r8), allocatable :: dp3d_tmp_tmp(:,:)
76 370944 : real (kind=r8), allocatable :: phis_tmp(:,:) ! temp array to hold phis
77 370944 : real (kind=r8), allocatable :: T_tmp(:,:,:) ! temp array to hold T
78 370944 : real (kind=r8), allocatable :: uv_tmp(:,:,:,:) ! temp array to hold u and v
79 370944 : real (kind=r8), allocatable :: q_tmp(:,:,:,:) ! temp to hold advected constituents
80 370944 : real (kind=r8), allocatable :: omega_tmp(:,:,:) ! temp array to hold omega
81 :
82 : ! Frontogenesis
83 370944 : real (kind=r8), allocatable :: frontgf(:,:,:) ! temp arrays to hold frontogenesis
84 370944 : real (kind=r8), allocatable :: frontga(:,:,:) ! function (frontgf) and angle (frontga)
85 370944 : real (kind=r8), allocatable :: frontgf_phys(:,:,:)
86 370944 : real (kind=r8), allocatable :: frontga_phys(:,:,:)
87 : ! Pointers to pbuf
88 370944 : real (kind=r8), pointer :: pbuf_frontgf(:,:)
89 370944 : real (kind=r8), pointer :: pbuf_frontga(:,:)
90 :
91 : integer :: ncols, ierr
92 : integer :: col_ind, blk_ind(1), m
93 : integer :: nphys
94 :
95 370944 : 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 370944 : type(physics_buffer_desc), pointer :: pbuf_chnk(:)
100 : !----------------------------------------------------------------------------
101 :
102 370944 : 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 370944 : elem => dyn_out%elem
107 370944 : tl_f = TimeLevel%n0
108 370944 : call TimeLevel_Qdp(TimeLevel, qsplit, tl_qdp_np0,tl_qdp_np1)
109 :
110 370944 : nullify(pbuf_chnk)
111 370944 : nullify(pbuf_frontgf)
112 370944 : nullify(pbuf_frontga)
113 :
114 370944 : if (fv_nphys > 0) then
115 370944 : 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 1483776 : allocate(ps_tmp(nphys_pts,nelemd))
123 1483776 : allocate(dp3d_tmp(nphys_pts,pver,nelemd))
124 1112832 : allocate(dp3d_tmp_tmp(nphys_pts,pver))
125 1112832 : allocate(phis_tmp(nphys_pts,nelemd))
126 1112832 : allocate(T_tmp(nphys_pts,pver,nelemd))
127 1854720 : allocate(uv_tmp(nphys_pts,2,pver,nelemd))
128 1854720 : allocate(q_tmp(nphys_pts,pver,pcnst,nelemd))
129 1112832 : allocate(omega_tmp(nphys_pts,pver,nelemd))
130 :
131 370944 : call tot_energy_dyn(elem,dyn_out%fvm, 1, nelemd,tl_f , tl_qdp_np0,'dBF')
132 :
133 370944 : if (use_gw_front .or. use_gw_front_igw) then
134 0 : allocate(frontgf(nphys_pts,pver,nelemd), stat=ierr)
135 0 : if (ierr /= 0) call endrun("dp_coupling: Allocate of frontgf failed.")
136 0 : allocate(frontga(nphys_pts,pver,nelemd), stat=ierr)
137 0 : if (ierr /= 0) call endrun("dp_coupling: Allocate of frontga failed.")
138 : end if
139 :
140 370944 : if (iam < par%nprocs) then
141 370944 : if (use_gw_front .or. use_gw_front_igw) then
142 0 : call gws_src_fnct(elem, tl_f, tl_qdp_np0, frontgf, frontga, nphys)
143 : end if
144 :
145 370944 : if (fv_nphys > 0) then
146 370944 : 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 370944 : 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 370944 : )
159 2979144 : do ie = 1, nelemd
160 : uv_tmp(:,:,:,ie) = &
161 2979144 : dyn2phys_vector(elem(ie)%state%v(:,:,:,:,tl_f),elem(ie))
162 : end do
163 370944 : 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 370944 : 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 370944 : if (.not. allocated(q_prev)) then
217 4608 : allocate(q_prev(pcols,pver,pcnst,begchunk:endchunk))
218 : end if
219 1989210384 : q_prev = 0.0_R8
220 :
221 370944 : call t_startf('dpcopy')
222 370944 : if (use_gw_front .or. use_gw_front_igw) then
223 0 : allocate(frontgf_phys(pcols, pver, begchunk:endchunk))
224 0 : 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 23844744 : do col_ind = 1, phys_columns_on_task
228 23473800 : call get_dyn_col_p(col_ind, ie, blk_ind)
229 23473800 : call get_chunk_info_p(col_ind, lchnk, icol)
230 23473800 : phys_state(lchnk)%ps(icol) = ps_tmp(blk_ind(1), ie)
231 23473800 : phys_state(lchnk)%phis(icol) = phis_tmp(blk_ind(1), ie)
232 633792600 : do ilyr = 1, pver
233 610318800 : phys_state(lchnk)%pdel(icol, ilyr) = dp3d_tmp(blk_ind(1), ilyr, ie)
234 610318800 : phys_state(lchnk)%t(icol, ilyr) = T_tmp(blk_ind(1), ilyr, ie)
235 610318800 : phys_state(lchnk)%u(icol, ilyr) = uv_tmp(blk_ind(1), 1, ilyr, ie)
236 610318800 : phys_state(lchnk)%v(icol, ilyr) = uv_tmp(blk_ind(1), 2, ilyr, ie)
237 610318800 : phys_state(lchnk)%omega(icol, ilyr) = omega_tmp(blk_ind(1), ilyr, ie)
238 :
239 633792600 : if (use_gw_front .or. use_gw_front_igw) then
240 0 : frontgf_phys(icol, ilyr, lchnk) = frontgf(blk_ind(1), ilyr, ie)
241 0 : frontga_phys(icol, ilyr, lchnk) = frontga(blk_ind(1), ilyr, ie)
242 : end if
243 : end do
244 :
245 117739944 : do m = 1, pcnst
246 1924851600 : do ilyr = 1, pver
247 1901377800 : 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 370944 : 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 0 : do lchnk = begchunk, endchunk
254 0 : ncols = get_ncols_p(lchnk)
255 0 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
256 0 : call pbuf_get_field(pbuf_chnk, frontgf_idx, pbuf_frontgf)
257 0 : call pbuf_get_field(pbuf_chnk, frontga_idx, pbuf_frontga)
258 0 : do icol = 1, ncols
259 0 : do ilyr = 1, pver
260 0 : pbuf_frontgf(icol, ilyr) = frontgf_phys(icol, ilyr, lchnk)
261 0 : pbuf_frontga(icol, ilyr) = frontga_phys(icol, ilyr, lchnk)
262 : end do
263 : end do
264 : end do
265 0 : deallocate(frontgf_phys)
266 0 : deallocate(frontga_phys)
267 : end if
268 :
269 370944 : 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 1866312 : do lchnk = begchunk, endchunk
274 1495368 : ncols = phys_state(lchnk)%ncol
275 1953947520 : q_prev(1:ncols,1:pver,1:pcnst,lchnk) = phys_state(lchnk)%q(1:ncols,1:pver,1:pcnst)
276 : end do
277 370944 : call test_mapping_output_phys_state(phys_state,dyn_out%fvm)
278 :
279 : ! Deallocate the temporary arrays
280 370944 : deallocate(ps_tmp)
281 370944 : deallocate(dp3d_tmp)
282 370944 : deallocate(phis_tmp)
283 370944 : deallocate(T_tmp)
284 370944 : deallocate(uv_tmp)
285 370944 : deallocate(q_tmp)
286 370944 : 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 370944 : call t_startf('derived_phys')
291 370944 : call derived_phys_dry(phys_state, phys_tend, pbuf2d)
292 370944 : call t_stopf('derived_phys')
293 :
294 : !$omp parallel do num_threads(max_num_threads) private (lchnk, ncols, ilyr, icol)
295 1866312 : do lchnk = begchunk, endchunk
296 1495368 : ncols=get_ncols_p(lchnk)
297 1866312 : if (pcols > ncols) then
298 1193976 : phys_state(lchnk)%phis(ncols+1:) = 0.0_r8
299 : end if
300 : end do
301 741888 : end subroutine d_p_coupling
302 :
303 : !=========================================================================================
304 :
305 738816 : 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 370944 : use phys_grid, only: get_dyn_col_p, columns_on_task, get_chunk_info_p, phys_columns_on_task
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 369408 : 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 369408 : real (kind=r8), allocatable :: dp_phys(:,:,:) ! temp array to hold dp on physics grid
332 369408 : real (kind=r8), allocatable :: T_tmp(:,:,:) ! temp array to hold T
333 369408 : real (kind=r8), allocatable :: dq_tmp(:,:,:,:) ! temp array to hold q
334 369408 : 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 369408 : 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 369408 : if (iam < par%nprocs) then
348 369408 : elem => dyn_in%elem
349 : else
350 0 : nullify(elem)
351 : end if
352 :
353 1477632 : allocate(T_tmp(nphys_pts,pver,nelemd))
354 1847040 : allocate(uv_tmp(nphys_pts,2,pver,nelemd))
355 1847040 : allocate(dq_tmp(nphys_pts,pver,pcnst,nelemd))
356 1108224 : allocate(dp_phys(nphys_pts,pver,nelemd))
357 :
358 678290808 : T_tmp = 0.0_r8
359 1421147208 : uv_tmp = 0.0_r8
360 2036731008 : dq_tmp = 0.0_r8
361 :
362 369408 : 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 1858584 : do lchnk = begchunk, endchunk
370 1489176 : ncols = get_ncols_p(lchnk)
371 25235184 : do icol = 1, ncols
372 632657376 : do ilyr = 1, pver
373 : ! convert wet mixing ratios to dry
374 607791600 : factor = phys_state(lchnk)%pdel(icol,ilyr)/phys_state(lchnk)%pdeldry(icol,ilyr)
375 2454543000 : do m = 1, pcnst
376 2431166400 : if (cnst_type(m) == 'wet') then
377 1823374800 : 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 369408 : 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 23746008 : do col_ind = 1, phys_columns_on_task
387 23376600 : call get_dyn_col_p(col_ind, ie, blk_ind)
388 23376600 : call get_chunk_info_p(col_ind, lchnk, icol)
389 :
390 : ! test code -- does nothing unless cpp macro debug_coupling is defined.
391 23376600 : call test_mapping_overwrite_tendencies(phys_state(lchnk), &
392 23376600 : phys_tend(lchnk), ncols, lchnk, q_prev(1:ncols,:,:,lchnk), &
393 70129800 : dyn_in%fvm)
394 :
395 654914208 : do ilyr = 1, pver
396 607791600 : dp_phys(blk_ind(1),ilyr,ie) = phys_state(lchnk)%pdeldry(icol,ilyr)
397 607791600 : T_tmp(blk_ind(1),ilyr,ie) = phys_tend(lchnk)%dtdt(icol,ilyr)
398 607791600 : uv_tmp(blk_ind(1),1,ilyr,ie) = phys_tend(lchnk)%dudt(icol,ilyr)
399 607791600 : uv_tmp(blk_ind(1),2,ilyr,ie) = phys_tend(lchnk)%dvdt(icol,ilyr)
400 2454543000 : do m = 1, pcnst
401 1823374800 : dq_tmp(blk_ind(1),ilyr,m,ie) = &
402 4254541200 : (phys_state(lchnk)%q(icol,ilyr,m) - q_prev(icol,ilyr,m,lchnk))
403 : end do
404 : end do
405 : end do
406 369408 : call t_stopf('pd_copy')
407 :
408 369408 : if (iam < par%nprocs) then
409 :
410 369408 : if (fv_nphys > 0) then
411 :
412 : ! put forcings into fvm structure
413 369408 : num_trac = max(qsize,ntrac)
414 2966808 : do ie = 1, nelemd
415 10759008 : do j = 1, fv_nphys
416 33766200 : do i = 1, fv_nphys
417 23376600 : ii = i + (j-1)*fv_nphys
418 631168200 : dyn_in%fvm(ie)%ft(i,j,1:pver) = T_tmp(ii,1:pver,ie)
419 1846751400 : dyn_in%fvm(ie)%fm(i,j,1:2,1:pver) = uv_tmp(ii,1:2,1:pver,ie)
420 1916881200 : dyn_in%fvm(ie)%fc_phys(i,j,1:pver,1:num_trac) = dq_tmp(ii,1:pver,1:num_trac,ie)
421 638960400 : 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 369408 : hybrid = config_thread_region(par,'serial')
429 369408 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
430 : !
431 : ! high-order mapping of ft and fm using fvm technology
432 : !
433 369408 : call t_startf('phys2dyn')
434 369408 : call phys2dyn_forcings_fvm(elem, dyn_in%fvm, hybrid,nets,nete,ntrac==0, tl_f, tl_qdp)
435 369408 : 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 369408 : deallocate(T_tmp)
455 369408 : deallocate(uv_tmp)
456 369408 : 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 369408 : call t_startf('p_d_coupling:bndry_exchange')
465 :
466 2966808 : do ie = 1, nelemd
467 2597400 : if (fv_nphys > 0) then
468 70129800 : do k = 1, nlev
469 : dyn_in%elem(ie)%derived%FM(:,:,1,k) = &
470 67532400 : dyn_in%elem(ie)%derived%FM(:,:,1,k) * &
471 1485712800 : 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 1418180400 : dyn_in%elem(ie)%spheremp(:,:)
475 : dyn_in%elem(ie)%derived%FT(:,:,k) = &
476 0 : dyn_in%elem(ie)%derived%FT(:,:,k) * &
477 1420777800 : dyn_in%elem(ie)%spheremp(:,:)
478 : end do
479 : end if
480 2597400 : kptr = 0
481 2597400 : call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
482 2597400 : kptr = kptr + 2*nlev
483 2597400 : call edgeVpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
484 2966808 : 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 369408 : if (iam < par%nprocs) then
495 369408 : call bndry_exchange(par, edgebuf, location='p_d_coupling')
496 : end if
497 :
498 2966808 : do ie = 1, nelemd
499 2597400 : kptr = 0
500 2597400 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FM(:,:,:,:), 2*nlev, kptr, ie)
501 2597400 : kptr = kptr + 2*nlev
502 2597400 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FT(:,:,:), nlev, kptr, ie)
503 2597400 : kptr = kptr + nlev
504 2597400 : if (.not. use_cslam) then
505 0 : call edgeVunpack(edgebuf, dyn_in%elem(ie)%derived%FQ(:,:,:,:), nlev*qsize, kptr, ie)
506 : end if
507 2966808 : if (fv_nphys > 0) then
508 70129800 : do k = 1, nlev
509 : dyn_in%elem(ie)%derived%FM(:,:,1,k) = &
510 67532400 : dyn_in%elem(ie)%derived%FM(:,:,1,k) * &
511 1485712800 : 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 1418180400 : dyn_in%elem(ie)%rspheremp(:,:)
515 : dyn_in%elem(ie)%derived%FT(:,:,k) = &
516 0 : dyn_in%elem(ie)%derived%FT(:,:,k) * &
517 1420777800 : dyn_in%elem(ie)%rspheremp(:,:)
518 : end do
519 : end if
520 : end do
521 369408 : call t_stopf('p_d_coupling:bndry_exchange')
522 :
523 369408 : if (iam < par%nprocs .and. fv_nphys > 0) then
524 738816 : call test_mapping_output_mapped_tendencies(dyn_in%fvm(1:nelemd), elem(1:nelemd), &
525 1108224 : 1, nelemd, tl_f, tl_qdp)
526 : end if
527 738816 : end subroutine p_d_coupling
528 :
529 : !=========================================================================================
530 :
531 370944 : 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 369408 : 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 370944 : 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 1866312 : do lchnk = begchunk,endchunk
578 :
579 1495368 : ncol = get_ncols_p(lchnk)
580 :
581 : ! dry pressure variables
582 :
583 24969168 : do i = 1, ncol
584 635287968 : phys_state(lchnk)%psdry(i) = hyai(1)*ps0 + sum(phys_state(lchnk)%pdel(i,:))
585 : end do
586 24969168 : do i = 1, ncol
587 24969168 : phys_state(lchnk)%pintdry(i,1) = hyai(1)*ps0
588 : end do
589 2990736 : call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,1), &
590 4486104 : phys_state(lchnk)%lnpintdry(1:ncol,1),ncol)
591 40374936 : do k = 1, nlev
592 649198368 : do i = 1, ncol
593 1220637600 : phys_state(lchnk)%pintdry(i,k+1) = phys_state(lchnk)%pintdry(i,k) + &
594 1869835968 : phys_state(lchnk)%pdel(i,k)
595 : end do
596 77759136 : call shr_vmath_log(phys_state(lchnk)%pintdry(1:ncol,k+1),&
597 118134072 : phys_state(lchnk)%lnpintdry(1:ncol,k+1),ncol)
598 : end do
599 :
600 40374936 : do k=1,nlev
601 649198368 : do i=1,ncol
602 610318800 : phys_state(lchnk)%pdeldry (i,k) = phys_state(lchnk)%pdel(i,k)
603 610318800 : 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 649198368 : phys_state(lchnk)%pintdry(i,k))
606 : end do
607 77759136 : call shr_vmath_log(phys_state(lchnk)%pmiddry(1:ncol,k), &
608 118134072 : phys_state(lchnk)%lnpmiddry(1:ncol,k),ncol)
609 : end do
610 :
611 : ! wet pressure variables (should be removed from physics!)
612 662448024 : factor_array(:,:) = 1.0_r8
613 5981472 : do m_cnst=dry_air_species_num + 1,thermodynamic_active_species_num
614 4486104 : m = thermodynamic_active_species_idx(m_cnst)
615 122620176 : do k=1,nlev
616 1952081208 : do i=1,ncol
617 : ! at this point all q's are dry
618 1947595104 : 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 40374936 : do k=1,nlev
624 650693736 : do i=1,ncol
625 649198368 : 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 24969168 : do i=1,ncol
632 23473800 : phys_state(lchnk)%ps(i) = phys_state(lchnk)%pintdry(i,1)
633 24969168 : phys_state(lchnk)%pint(i,1) = phys_state(lchnk)%pintdry(i,1)
634 : end do
635 40374936 : do k = 1, nlev
636 649198368 : do i=1,ncol
637 610318800 : phys_state(lchnk)%pint(i,k+1) = phys_state(lchnk)%pint(i,k)+phys_state(lchnk)%pdel(i,k)
638 610318800 : phys_state(lchnk)%pmid(i,k) = (phys_state(lchnk)%pint(i,k+1)+phys_state(lchnk)%pint(i,k))/2._r8
639 649198368 : phys_state(lchnk)%ps (i) = phys_state(lchnk)%ps(i) + phys_state(lchnk)%pdel(i,k)
640 : end do
641 38879568 : call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,k),phys_state(lchnk)%lnpint(1:ncol,k),ncol)
642 40374936 : call shr_vmath_log(phys_state(lchnk)%pmid(1:ncol,k),phys_state(lchnk)%lnpmid(1:ncol,k),ncol)
643 : end do
644 1495368 : call shr_vmath_log(phys_state(lchnk)%pint(1:ncol,pverp),phys_state(lchnk)%lnpint(1:ncol,pverp),ncol)
645 :
646 40374936 : do k = 1, nlev
647 650693736 : do i = 1, ncol
648 649198368 : 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 1495368 : 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 662448024 : zvirv(:,:) = zvir
667 : end if
668 : !
669 : ! update cp_dycore in module air_composition.
670 : ! (note: at this point q is dry)
671 : !
672 1495368 : call cam_thermo_water_update(phys_state(lchnk)%q(1:ncol,:,:), lchnk, ncol, vc_dry_pressure)
673 40374936 : do k = 1, nlev
674 650693736 : do i = 1, ncol
675 610318800 : phys_state(lchnk)%exner(i,k) = (phys_state(lchnk)%pint(i,pver+1) &
676 1259517168 : / 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 650693736 : factor_array(1:ncol,1:nlev) = 1._r8/factor_array(1:ncol,1:nlev)
683 5981472 : do m = 1,pcnst
684 5981472 : if (cnst_type(m) == 'wet') then
685 121124808 : do k = 1, nlev
686 1952081208 : do i = 1, ncol
687 1947595104 : 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 1495368 : 1, pcnst, qmin ,phys_state(lchnk)%q)
696 :
697 : ! Compute initial geopotential heights - based on full pressure
698 1495368 : 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 2990736 : phys_state(lchnk)%zi , phys_state(lchnk)%zm , ncol)
702 : ! Compute initial dry static energy, include surface geopotential
703 2990736 : call update_dry_static_energy_run(pver, gravit, phys_state(lchnk)%t(1:ncol,:), &
704 1495368 : phys_state(lchnk)%zm(1:ncol,:), &
705 1495368 : phys_state(lchnk)%phis(1:ncol), &
706 1495368 : phys_state(lchnk)%s(1:ncol,:), &
707 8972208 : cpairv(1:ncol,:,lchnk), errflg, errmsg)
708 : ! Compute energy and water integrals of input state
709 1495368 : pbuf_chnk => pbuf_get_chunk(pbuf2d, lchnk)
710 1866312 : call check_energy_timestep_init(phys_state(lchnk), phys_tend(lchnk), pbuf_chnk)
711 :
712 :
713 : end do ! lchnk
714 :
715 741888 : end subroutine derived_phys_dry
716 : end module dp_coupling
|
