Line data Source code
1 : !
2 : ! pg3->GLL and GLL->pg3 mapping algorithm described in:
3 : !
4 : ! Adam R. Herrington, Peter H. Lauritzen, Mark A. Taylor, Steve Goldhaber, Brian Eaton, Kevin A Reed and Paul A. Ullrich, 2018:
5 : ! Physics-dynamics coupling with element-based high-order Galerkin methods: quasi equal-area physics grid:
6 : ! Mon. Wea. Rev., DOI:MWR-D-18-0136.1
7 : !
8 : ! pg2->pg3 mapping algorithm described in:
9 : !
10 : ! Adam R. Herrington, Peter H. Lauritzen, Kevin A Reed, Steve Goldhaber, and Brian Eaton, 2019:
11 : ! Exploring a lower resolution physics grid in CAM-SE-CSLAM. J. Adv. Model. Earth Syst.
12 : !
13 : !#define PCoM !replace PPM with PCoM for mass variables for fvm2phys and phys2fvm
14 : !#define skip_high_order_fq_map !do mass and correlation preserving phys2fvm mapping but no high-order pre-mapping of fq
15 : #define mass_fix
16 : module fvm_mapping
17 : use shr_kind_mod, only: r8=>shr_kind_r8
18 : use dimensions_mod, only: irecons_tracer
19 : use element_mod, only: element_t
20 : use fvm_control_volume_mod, only: fvm_struct
21 : use perf_mod, only: t_startf, t_stopf
22 : use cam_abortutils, only: endrun
23 : use cam_logfile, only: iulog
24 : implicit none
25 : private
26 :
27 : public :: phys2dyn_forcings_fvm, dyn2phys, dyn2phys_vector, dyn2phys_all_vars,dyn2fvm_mass_vars
28 : public :: phys2dyn,fvm2dyn,dyn2fvm,cslam2gll
29 : save
30 : integer :: save_max_overlap
31 : real(kind=r8), allocatable, dimension(:,:,:,:,:) :: save_air_mass_overlap
32 : real(kind=r8), allocatable, dimension(:,:,:,:,:,:) :: save_q_overlap
33 : real(kind=r8), allocatable, dimension(:,:,:,:,:) :: save_q_phys
34 : real(kind=r8), allocatable, dimension(:,:,:,:) :: save_dp_phys
35 : real(kind=r8), allocatable, dimension(:,:,:,:) :: save_overlap_area
36 : integer , allocatable, dimension(:,:,:,:,:) :: save_overlap_idx
37 : integer , allocatable, dimension(:,:,:,:) :: save_num_overlap
38 :
39 : interface fvm2dyn
40 : module procedure fvm2dynt1
41 : module procedure fvm2dyntn
42 : end interface fvm2dyn
43 :
44 : contains
45 : !
46 : ! map all mass variables from gll to fvm
47 : !
48 369408 : subroutine phys2dyn_forcings_fvm(elem, fvm, hybrid,nets,nete,no_cslam, tl_f, tl_qdp)
49 : use dimensions_mod, only: np, nc,nlev
50 : use dimensions_mod, only: fv_nphys, nhc_phys,ntrac,nhc,ksponge_end, nu_scale_top
51 : use hybrid_mod, only: hybrid_t
52 : use air_composition, only: thermodynamic_active_species_num, thermodynamic_active_species_idx
53 : type (element_t), intent(inout):: elem(:)
54 : type(fvm_struct), intent(inout):: fvm(:)
55 :
56 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
57 : logical, intent(in) :: no_cslam
58 : integer, intent(in) :: nets, nete, tl_f, tl_qdp
59 :
60 : integer :: ie,i,j,k,m_cnst,nq
61 369408 : real (kind=r8), dimension(:,:,:,:,:) , allocatable :: fld_phys, fld_gll
62 : real (kind=r8) :: element_ave
63 : !
64 : ! for tensor product Lagrange interpolation
65 : !
66 : integer :: nflds
67 369408 : logical, allocatable :: llimiter(:)
68 :
69 : integer :: ierr
70 :
71 369408 : if (no_cslam) then
72 0 : call endrun("phys2dyn_forcings_fvm: no cslam case: NOT SUPPORTED")
73 369408 : else if (nc.ne.fv_nphys) then
74 : !
75 : !***********************************************************
76 : !
77 : ! using cslam and different resolution physics grid
78 : !
79 : !***********************************************************
80 : !
81 0 : call t_startf('p2d-pg2:copying')
82 0 : nflds = 4+ntrac
83 0 : allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,nflds,nets:nete), stat=ierr)
84 0 : if( ierr /= 0 ) then
85 0 : write(iulog,*) 'phys2dyn_forcings_fvm: fld_phys allocation error = ',ierr
86 0 : call endrun('phys2dyn_forcings_fvm: failed to allocate fld_phys array')
87 : end if
88 0 : allocate(fld_gll(np,np,nlev,3,nets:nete), stat=ierr)
89 0 : if( ierr /= 0 ) then
90 0 : write(iulog,*) 'phys2dyn_forcings_fvm: fld_gll allocation error = ',ierr
91 0 : call endrun('phys2dyn_forcings_fvm: failed to allocate fld_gll array')
92 : end if
93 0 : allocate(llimiter(3), stat=ierr)
94 0 : if( ierr /= 0 ) then
95 0 : write(iulog,*) 'phys2dyn_forcings_fvm: llimiter allocation error = ',ierr
96 0 : call endrun('phys2dyn_forcings_fvm: failed to allocate llimiter array')
97 : end if
98 0 : fld_phys = -9.99E99_r8!xxx necessary?
99 :
100 0 : llimiter = .false.
101 :
102 0 : do ie=nets,nete
103 : !
104 : ! pack fields that need to be interpolated
105 : !
106 0 : fld_phys(1:fv_nphys,1:fv_nphys,:,1,ie) = fvm(ie)%ft(1:fv_nphys,1:fv_nphys,:)
107 0 : fld_phys(1:fv_nphys,1:fv_nphys,:,2,ie) = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,1,:)
108 0 : fld_phys(1:fv_nphys,1:fv_nphys,:,3,ie) = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,2,:)
109 0 : fld_phys(1:fv_nphys,1:fv_nphys,:,4,ie) = fvm(ie)%dp_phys(1:fv_nphys,1:fv_nphys,:)
110 0 : do m_cnst=1,ntrac
111 0 : fld_phys(1:fv_nphys,1:fv_nphys,:,4+m_cnst,ie) = &
112 0 : fvm(ie)%fc_phys(1:fv_nphys,1:fv_nphys,:,m_cnst)
113 : end do
114 : end do
115 0 : call t_stopf('p2d-pg2:copying')
116 0 : call t_startf('p2d-pg2:fill_halo_phys')
117 0 : call fill_halo_phys(fld_phys,hybrid,nets,nete,nlev,nflds)
118 : !
119 : ! do mapping of fu,fv,ft
120 : !
121 0 : call phys2dyn(hybrid,elem,fld_phys(:,:,:,1:3,:),fld_gll,nets,nete,nlev,3,fvm,llimiter, &
122 0 : istart_vector=2,halo_filled=.true.)
123 0 : do ie=nets,nete
124 0 : elem(ie)%derived%fT(:,:,:) = fld_gll(:,:,:,1,ie)
125 0 : elem(ie)%derived%fM(:,:,1,:) = fld_gll(:,:,:,2,ie)
126 0 : elem(ie)%derived%fM(:,:,2,:) = fld_gll(:,:,:,3,ie)
127 : end do
128 0 : call t_stopf('p2d-pg2:fill_halo_phys')
129 :
130 0 : deallocate(fld_gll)
131 : !
132 : ! map fq from phys to fvm
133 : !
134 0 : call t_startf('p2d-pg2:phys2fvm')
135 :
136 0 : do ie=nets,nete
137 0 : do k=1,nlev
138 0 : call phys2fvm(ie,k,fvm(ie),&
139 0 : fld_phys(:,:,k,5:4+ntrac,ie),fvm(ie)%fc(:,:,k,1:ntrac),ntrac)
140 : end do
141 : end do
142 0 : call t_stopf('p2d-pg2:phys2fvm')
143 : !deallocate arrays allocated in dyn2phys_all_vars
144 0 : deallocate(save_air_mass_overlap,save_q_phys,save_q_overlap,&
145 0 : save_overlap_area,save_num_overlap,save_overlap_idx,save_dp_phys)
146 : else
147 : !
148 : !
149 : !*****************************************************************************************
150 : !
151 : ! using cslam with same physics grid resolution as cslam resolution
152 : !
153 : !*****************************************************************************************
154 : !
155 : ! nflds is ft, fu, fv, + thermo species
156 369408 : nflds = 3
157 2585856 : allocate(fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,nlev,nflds,nets:nete))
158 1108224 : allocate(fld_gll(np,np,nlev,nflds,nets:nete))
159 369408 : allocate(llimiter(nflds))
160 1477632 : llimiter(1:nflds) = .false.
161 2966808 : do ie=nets,nete
162 : !
163 : ! pack fields that need to be interpolated
164 : !
165 880518600 : fld_phys(1:fv_nphys,1:fv_nphys,:,1,ie) = fvm(ie)%ft(1:fv_nphys,1:fv_nphys,:)
166 880518600 : fld_phys(1:fv_nphys,1:fv_nphys,:,2,ie) = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,1,:)
167 880888008 : fld_phys(1:fv_nphys,1:fv_nphys,:,3,ie) = fvm(ie)%fm(1:fv_nphys,1:fv_nphys,2,:)
168 : end do
169 : !
170 : ! do mapping
171 : !
172 369408 : call phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,nlev,nflds,fvm,llimiter,2)
173 2966808 : do ie=nets,nete
174 1420777800 : elem(ie)%derived%fT(:,:,:) = fld_gll(:,:,:,1,ie)
175 1420777800 : elem(ie)%derived%fM(:,:,1,:) = fld_gll(:,:,:,2,ie)
176 1421147208 : elem(ie)%derived%fM(:,:,2,:) = fld_gll(:,:,:,3,ie)
177 : end do
178 369408 : deallocate(fld_gll)
179 2966808 : do ie=nets,nete
180 10759008 : do m_cnst = 1,ntrac
181 2644153200 : fvm(ie)%fc(1:nc,1:nc,:,m_cnst) = fvm(ie)%fc_phys(1:nc,1:nc,:,m_cnst)*fvm(ie)%dp_fvm(1:nc,1:nc,:)
182 : end do
183 : end do
184 : end if
185 369408 : deallocate(fld_phys,llimiter)
186 369408 : end subroutine phys2dyn_forcings_fvm
187 :
188 : ! for multiple fields
189 1477632 : subroutine fvm2dyntn(fld_fvm,fld_gll,hybrid,nets,nete,numlev,num_flds,fvm,llimiter,halo_filled)
190 : use dimensions_mod, only: np, nhc, nc
191 : use hybrid_mod , only: hybrid_t
192 : use bndry_mod , only: ghost_exchange
193 : use edge_mod , only: ghostpack,ghostunpack
194 : use fvm_mod , only: ghostBufQnhc_s
195 : !
196 : integer , intent(in) :: nets,nete,num_flds,numlev
197 : real (kind=r8), intent(inout) :: fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,numlev,num_flds,nets:nete)
198 : real (kind=r8), intent(out) :: fld_gll(np,np,numlev,num_flds,nets:nete)
199 : type (hybrid_t) , intent(in) :: hybrid
200 : type(fvm_struct) , intent(in) :: fvm(nets:nete)
201 : logical , intent(in) :: llimiter(num_flds)
202 : logical, optional , intent(in) :: halo_filled !optional if boundary exchange for fld_fvm has already been called
203 :
204 : integer :: ie, iwidth
205 : logical :: fill_halo
206 : !
207 : !*********************************************
208 : !
209 : ! halo exchange
210 : !
211 : !*********************************************
212 : !
213 738816 : fill_halo = .true.
214 738816 : if (present(halo_filled)) then
215 738816 : fill_halo = .not. halo_filled
216 : end if
217 :
218 738816 : if (fill_halo) then
219 0 : do ie=nets,nete
220 0 : call ghostpack(ghostBufQnhc_s, fld_fvm(:,:,:,:,ie),numlev*num_flds,0,ie)
221 : end do
222 0 : call ghost_exchange(hybrid,ghostbufQnhc_s,location='fvm2dyntn')
223 0 : do ie=nets,nete
224 0 : call ghostunpack(ghostbufQnhc_s, fld_fvm(:,:,:,:,ie),numlev*num_flds,0,ie)
225 : end do
226 : end if
227 : !
228 : ! mapping
229 : !
230 738816 : iwidth=2
231 : ! iwidth=1 !low-order mapping
232 5933616 : do ie=nets,nete
233 5194800 : call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,nc,nhc,numlev,num_flds,fld_fvm(:,:,:,:,ie),&
234 11128416 : fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord)
235 : end do
236 738816 : end subroutine fvm2dyntn
237 :
238 : ! for single field
239 0 : subroutine fvm2dynt1(fld_fvm,fld_gll,hybrid,nets,nete,numlev,fvm,llimiter,halo_filled)
240 738816 : use dimensions_mod, only: np, nhc, nc
241 : use hybrid_mod , only: hybrid_t
242 : use bndry_mod , only: ghost_exchange
243 : use edge_mod , only: ghostpack,ghostunpack
244 : use fvm_mod , only: ghostBufQnhc_t1
245 : !
246 : integer , intent(in) :: nets,nete,numlev
247 : real (kind=r8), intent(inout) :: fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,numlev,1,nets:nete)
248 : real (kind=r8), intent(out) :: fld_gll(np,np,numlev,1,nets:nete)
249 : type (hybrid_t) , intent(in) :: hybrid
250 : type(fvm_struct) , intent(in) :: fvm(nets:nete)
251 : logical , intent(in) :: llimiter(1)
252 : logical, optional , intent(in) :: halo_filled!optional if boundary exchange for fld_fvm has already been called
253 :
254 : integer :: ie, iwidth
255 : logical :: fill_halo
256 : !
257 : !*********************************************
258 : !
259 : ! halo exchange
260 : !
261 : !*********************************************
262 : !
263 0 : fill_halo = .true.
264 0 : if (present(halo_filled)) then
265 0 : fill_halo = .not. halo_filled
266 : end if
267 :
268 0 : if (fill_halo) then
269 0 : do ie=nets,nete
270 0 : call ghostpack(ghostBufQnhc_t1, fld_fvm(:,:,:,1,ie),numlev,0,ie)
271 : end do
272 0 : call ghost_exchange(hybrid,ghostbufQnhc_t1,location='fvm2dynt1')
273 0 : do ie=nets,nete
274 0 : call ghostunpack(ghostbufQnhc_t1, fld_fvm(:,:,:,1,ie),numlev,0,ie)
275 : end do
276 : end if
277 : !
278 : ! mapping
279 : !
280 0 : iwidth=2
281 0 : do ie=nets,nete
282 0 : call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,nc,nhc,numlev,1,fld_fvm(:,:,:,:,ie),&
283 0 : fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord)
284 : end do
285 0 : end subroutine fvm2dynt1
286 :
287 738816 : subroutine fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
288 0 : use dimensions_mod, only: nhc_phys, fv_nphys
289 : use hybrid_mod , only: hybrid_t
290 : use bndry_mod , only: ghost_exchange
291 : use edge_mod , only: ghostpack, ghostunpack
292 : use fvm_mod , only: ghostBufPG_s
293 :
294 : integer , intent(in) :: nets,nete,num_lev,num_flds
295 : real (kind=r8), intent(inout) :: fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_lev,num_flds, &
296 : nets:nete)
297 : type (hybrid_t) , intent(in) :: hybrid ! distributed parallel structure (shared)
298 :
299 : integer :: ie
300 : !
301 : !*********************************************
302 : !
303 : ! halo exchange
304 : !
305 : !*********************************************
306 : !
307 369408 : call t_startf('fvm:fill_halo_phys')
308 2966808 : do ie=nets,nete
309 2966808 : call ghostpack(ghostBufPG_s, fld_phys(:,:,:,:,ie),num_lev*num_flds,0,ie)
310 : end do
311 :
312 369408 : call ghost_exchange(hybrid,ghostBufPG_s,location='fill_halo_phys')
313 :
314 2966808 : do ie=nets,nete
315 2966808 : call ghostunpack(ghostBufPG_s, fld_phys(:,:,:,:,ie),num_lev*num_flds,0,ie)
316 : end do
317 : !
318 369408 : call t_stopf('fvm:fill_halo_phys')
319 369408 : end subroutine fill_halo_phys
320 : !
321 : ! must call fill_halo_phys before calling this subroutine
322 : !
323 369408 : subroutine phys2dyn(hybrid,elem,fld_phys,fld_gll,nets,nete,num_lev,num_flds,fvm,llimiter,istart_vector,halo_filled)
324 369408 : use dimensions_mod, only: np, nhc_phys, fv_nphys
325 : use hybrid_mod, only : hybrid_t
326 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
327 : integer , intent(in) :: nets,nete,num_flds,num_lev
328 : real (kind=r8), intent(inout) :: fld_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_lev,num_flds, &
329 : nets:nete)
330 : real (kind=r8), intent(out) :: fld_gll(np,np,num_lev,num_flds,nets:nete)
331 : type (element_t) , intent(inout) :: elem(:)
332 : type(fvm_struct) , intent(in) :: fvm(:)
333 : integer, optional , intent(in) :: istart_vector
334 : logical , intent(in) :: llimiter(num_flds)
335 : logical, optional , intent(in) :: halo_filled!optional if boundary exchange for fld_fvm has already been called
336 :
337 : integer :: i, j, ie, k, iwidth
338 : real (kind=r8) :: v1,v2
339 :
340 369408 : if (present(halo_filled)) then
341 0 : if (.not.halo_filled) call fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
342 : else
343 369408 : call fill_halo_phys(fld_phys,hybrid,nets,nete,num_lev,num_flds)
344 : end if
345 369408 : if (present(istart_vector)) then
346 2966808 : do ie=nets,nete
347 70499208 : do k=1,num_lev
348 677921400 : do j=1-nhc_phys,fv_nphys+nhc_phys
349 6145448400 : do i=1-nhc_phys,fv_nphys+nhc_phys
350 : !
351 : ! convert lat-lon vectors to contra-variant gnomonic
352 : !
353 5470124400 : v1 = fld_phys(i,j,k,istart_vector ,ie)
354 5470124400 : v2 = fld_phys(i,j,k,istart_vector+1,ie)
355 5470124400 : fld_phys(i,j,k,istart_vector ,ie)=fvm(ie)%Dinv_physgrid(i,j,1,1)*v1 + fvm(ie)%Dinv_physgrid(i,j,1,2)*v2
356 6077916000 : fld_phys(i,j,k,istart_vector+1,ie)=fvm(ie)%Dinv_physgrid(i,j,2,1)*v1 + fvm(ie)%Dinv_physgrid(i,j,2,2)*v2
357 : end do
358 : end do
359 : end do
360 : end do
361 : end if
362 : !
363 : ! mapping
364 : !
365 369408 : iwidth=2
366 : ! iwidth=1
367 369408 : if (fv_nphys==1) iwidth=1
368 2966808 : do ie=nets,nete
369 2597400 : call tensor_lagrange_interp(fvm(ie)%cubeboundary,np,fv_nphys,nhc_phys,num_lev,num_flds,fld_phys(:,:,:,:,ie),&
370 5564208 : fld_gll(:,:,:,:,ie),llimiter,iwidth,fvm(ie)%norm_elem_coord_physgrid)
371 : end do
372 :
373 369408 : if (present(istart_vector)) then
374 : !
375 : ! convert contra-variant to lat-lon
376 : !
377 2966808 : do ie=nets,nete
378 70499208 : do k=1,num_lev
379 340259400 : do j=1,np
380 1418180400 : do i=1,np
381 1080518400 : v1 = fld_gll(i,j,k,istart_vector ,ie)
382 1080518400 : v2 = fld_gll(i,j,k,istart_vector+1,ie)
383 1080518400 : fld_gll(i,j,k,istart_vector ,ie) = elem(ie)%D(i,j,1,1)*v1 + elem(ie)%D(i,j,1,2)*v2
384 1350648000 : fld_gll(i,j,k,istart_vector+1,ie) = elem(ie)%D(i,j,2,1)*v1 + elem(ie)%D(i,j,2,2)*v2
385 : end do
386 : end do
387 : end do
388 : end do
389 : end if
390 369408 : end subroutine phys2dyn
391 : !
392 : ! map all mass variables from gll to fvm
393 : !
394 5400 : subroutine dyn2fvm_mass_vars(dp_gll,ps_gll,q_gll,&
395 5400 : dp_fvm,ps_fvm,q_fvm,num_trac,metdet,inv_area)
396 : use dimensions_mod, only: np, nc,nlev
397 : integer, intent(in) :: num_trac
398 : real (kind=r8), dimension(np,np,nlev) , intent(in) :: dp_gll
399 : real (kind=r8), dimension(np,np,nlev,num_trac), intent(in) :: q_gll
400 : real (kind=r8), dimension(np,np) , intent(in) :: ps_gll
401 :
402 :
403 : real (kind=r8), dimension(nc,nc,nlev) , intent(inout) :: dp_fvm
404 : real (kind=r8), dimension(nc,nc,nlev,num_trac), intent(inout) :: q_fvm
405 : real (kind=r8), dimension(nc,nc) , intent(inout) :: ps_fvm
406 : real (kind=r8), dimension(nc,nc) , intent(out) :: inv_area
407 :
408 : real (kind=r8), intent(in) :: metdet(np,np)
409 :
410 : real (kind=r8) :: se_area_sphere(nc,nc), tmp(np,np)
411 : real (kind=r8) :: inv_darea_dp_fvm(nc,nc)
412 : integer :: k,m_cnst
413 :
414 113400 : tmp = 1.0_r8
415 5400 : se_area_sphere = dyn2fvm(tmp,metdet)
416 70200 : inv_area = 1.0_r8/se_area_sphere
417 :
418 70200 : ps_fvm(:,:) = dyn2fvm(ps_gll,metdet,inv_area)
419 145800 : do k=1,nlev
420 1825200 : dp_fvm(:,:,k) = dyn2fvm(dp_gll(:,:,k),metdet,inv_area)
421 1825200 : inv_darea_dp_fvm = inv_area/dp_fvm(:,:,k)
422 567000 : do m_cnst=1,num_trac
423 : q_fvm(:,:,k,m_cnst) = &
424 421200 : dyn2fvm(q_gll(:,:,k,m_cnst)*dp_gll(:,:,k),metdet,&
425 14461200 : inv_darea_dp_fvm,q_gll(:,:,k,m_cnst))
426 : end do
427 : end do
428 5400 : end subroutine dyn2fvm_mass_vars
429 :
430 : !
431 : ! this subroutine assumes that the fvm halo has already been filled
432 : ! (if nc/=fv_nphys)
433 : !
434 :
435 370944 : subroutine dyn2phys_all_vars(nets,nete,elem,fvm,&
436 : num_trac,ptop,tl,&
437 370944 : dp3d_phys,ps_phys,q_phys,T_phys,omega_phys,phis_phys)
438 : use dimensions_mod, only: np, nc,nlev,fv_nphys
439 : use dp_mapping, only: nphys_pts
440 : use element_mod, only: element_t
441 : integer, intent(in) :: nets,nete,num_trac,tl
442 :
443 : type(fvm_struct), dimension(nets:nete), intent(inout):: fvm
444 : type(element_t), dimension(nets:nete), intent(in) :: elem
445 :
446 : real (kind=r8), intent(in) :: ptop
447 :
448 : real (kind=r8), dimension(nphys_pts,nets:nete) , intent(out) :: ps_phys,phis_phys
449 : real (kind=r8), dimension(nphys_pts,nlev,nets:nete) , intent(out) :: dp3d_phys,T_phys,omega_phys
450 : real (kind=r8), dimension(nphys_pts,nlev,num_trac,nets:nete) , intent(out) :: q_phys
451 :
452 :
453 : real (kind=r8) :: tmp(np,np)
454 741888 : real (kind=r8), dimension(fv_nphys,fv_nphys) :: inv_area,inv_darea_dp_phys,dp3d_tmp
455 741888 : real (kind=r8), dimension(fv_nphys,fv_nphys,num_trac) :: q_phys_tmp
456 : real (kind=r8), dimension(nc,nc) :: inv_darea_dp_fvm
457 : integer :: k,m_cnst,ie
458 :
459 :
460 :
461 : !OMP BARRIER OMP MASTER needed
462 370944 : if (nc.ne.fv_nphys) then
463 0 : save_max_overlap = 4 !max number of mass overlap areas between phys and fvm grids
464 0 : allocate(save_air_mass_overlap(save_max_overlap,fv_nphys,fv_nphys,nlev,nets:nete))
465 0 : allocate(save_q_overlap(save_max_overlap,fv_nphys,fv_nphys,nlev,num_trac,nets:nete))
466 0 : allocate(save_q_phys(fv_nphys,fv_nphys,nlev,num_trac,nets:nete))
467 0 : allocate(save_dp_phys(fv_nphys,fv_nphys,nlev,nets:nete))
468 0 : allocate(save_overlap_area(save_max_overlap,fv_nphys,fv_nphys,nets:nete))
469 0 : allocate(save_num_overlap(fv_nphys,fv_nphys,nlev,nets:nete))
470 0 : save_num_overlap = 0
471 0 : allocate(save_overlap_idx(2,save_max_overlap,fv_nphys,fv_nphys,nets:nete))
472 : end if
473 :
474 2979144 : do ie=nets,nete
475 54772200 : tmp = 1.0_r8
476 33906600 : inv_area = 1.0_r8/dyn2phys(tmp,elem(ie)%metdet(:,:))
477 5216400 : phis_phys(:,ie) = RESHAPE(dyn2phys(elem(ie)%state%phis(:,:),elem(ie)%metdet(:,:),inv_area),SHAPE(phis_phys(:,ie)))
478 26082000 : ps_phys(:,ie) = ptop
479 2979144 : if (nc.ne.fv_nphys) then
480 0 : tmp = 1.0_r8
481 0 : do k=1,nlev
482 0 : inv_darea_dp_fvm = dyn2fvm(elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:))
483 0 : inv_darea_dp_fvm = 1.0_r8/inv_darea_dp_fvm
484 0 : T_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%state%T(:,:,k,tl),elem(ie)%metdet(:,:),inv_area),SHAPE(T_phys(:,k,ie)))
485 0 : Omega_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%derived%omega(:,:,k),elem(ie)%metdet(:,:),inv_area), &
486 0 : SHAPE(Omega_phys(:,k,ie)))
487 0 : call fvm2phys(ie,k,fvm(ie),fvm(ie)%c(:,:,k,:),q_phys_tmp,num_trac)
488 0 : dp3d_phys(:,k,ie) = RESHAPE(save_dp_phys(:,:,k,ie),SHAPE(dp3d_phys(:,k,ie)))
489 0 : ps_phys(:,ie) = ps_phys(:,ie)+RESHAPE(save_dp_phys(:,:,k,ie),SHAPE(ps_phys(:,ie)))
490 0 : do m_cnst=1,num_trac
491 0 : q_phys(:,k,m_cnst,ie) = RESHAPE(q_phys_tmp(:,:,m_cnst),SHAPE(q_phys(:,k,m_cnst,ie)))
492 : end do
493 : end do
494 : else
495 70421400 : do k=1,nlev
496 67813200 : dp3d_tmp = dyn2phys(elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:),inv_area)
497 881571600 : inv_darea_dp_phys = inv_area/dp3d_tmp
498 203439600 : T_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%state%T(:,:,k,tl)*elem(ie)%state%dp3d(:,:,k,tl),elem(ie)%metdet(:,:),&
499 1695330000 : inv_darea_dp_phys),SHAPE(T_phys(:,k,ie)))
500 67813200 : Omega_phys(:,k,ie) = RESHAPE(dyn2phys(elem(ie)%derived%OMEGA(:,:,k),elem(ie)%metdet(:,:),inv_area), &
501 203439600 : SHAPE(Omega_phys(:,k,ie)))
502 : !
503 : ! no mapping needed - just copy fields into physics structure
504 : !
505 135626400 : dp3d_phys(:,k,ie) = RESHAPE(fvm(ie)%dp_fvm(1:nc,1:nc,k),SHAPE(dp3d_phys(:,k,ie)))
506 745945200 : ps_phys(:,ie) = ps_phys(:,ie)+RESHAPE(fvm(ie)%dp_fvm(1:nc,1:nc,k),SHAPE(ps_phys(:,ie)))
507 273861000 : do m_cnst=1,num_trac
508 474692400 : q_phys(:,k,m_cnst,ie) = RESHAPE(fvm(ie)%c(1:nc,1:nc,k,m_cnst),SHAPE(q_phys(:,k,m_cnst,ie)))
509 : end do
510 : end do
511 : end if
512 : end do
513 370944 : end subroutine dyn2phys_all_vars
514 :
515 :
516 208656000 : function dyn2phys(qdp_gll,metdet,inv_dp_darea_phys) result(qdp_phys)
517 : use dimensions_mod, only: np, nc, fv_nphys
518 : use derivative_mod, only: subcell_integration
519 : real (kind=r8), intent(in) :: qdp_gll(np,np)
520 : real (kind=r8) :: qdp_phys(fv_nphys,fv_nphys)
521 : real (kind=r8), intent(in) :: metdet(np,np)
522 : real (kind=r8), intent(in), optional :: inv_dp_darea_phys(fv_nphys,fv_nphys)
523 :
524 208656000 : call subcell_integration(qdp_gll(:,:), np, fv_nphys, metdet,qdp_phys,nc.ne.fv_nphys)
525 208656000 : if (present(inv_dp_darea_phys)) &
526 2884669200 : qdp_phys = qdp_phys*inv_dp_darea_phys ! convert qdp to q
527 625968000 : end function dyn2phys
528 :
529 :
530 572400 : function dyn2fvm(qdp_gll,metdet,inv_dp_darea_phys,q_gll) result(qdp_phys)
531 208656000 : use dimensions_mod, only: np, nc
532 : use derivative_mod, only: subcell_integration
533 : real (kind=r8), intent(in) :: qdp_gll(np,np)
534 : real (kind=r8), intent(in) :: metdet(np,np)
535 : real (kind=r8), intent(in), optional :: inv_dp_darea_phys(nc,nc)
536 : real (kind=r8), intent(in), optional :: q_gll(np,np)
537 : real (kind=r8) :: qdp_phys(nc,nc), min_val, max_val
538 : integer :: i,j
539 :
540 572400 : call subcell_integration(qdp_gll(:,:), np, nc, metdet,qdp_phys)
541 572400 : if (present(inv_dp_darea_phys)) then
542 : !
543 : ! convert qdp to q
544 : !
545 7371000 : qdp_phys = qdp_phys*inv_dp_darea_phys
546 : !
547 : ! simple limiter
548 : !
549 567000 : if (present(q_gll)) then
550 8845200 : min_val = minval(q_gll)
551 8845200 : max_val = maxval(q_gll)
552 1684800 : do j = 1, nc
553 5475600 : do i = 1, nc
554 : !
555 : ! simple limiter: only coded for nc=3 and np4
556 : !
557 5054400 : qdp_phys(i,j) = max(min_val,min(max_val,qdp_phys(i,j)))
558 : end do
559 : end do
560 : end if
561 : end if
562 1144800 : end function dyn2fvm
563 :
564 13041000 : function dyn2phys_vector(v_gll,elem) result(v_phys)
565 572400 : use dimensions_mod, only: np, nlev, fv_nphys
566 : use interpolate_mod,only: interpdata_t,interpolate_2d,interpolate_t
567 : use cube_mod ,only: dmap
568 : use control_mod ,only: cubed_sphere_map
569 :
570 : type (interpdata_t):: interpdata
571 : type (element_t), intent(in) :: elem
572 : type (interpolate_t) , target :: interp_p
573 : real (kind=r8), intent(in) :: v_gll(np,np,2,nlev)
574 : real (kind=r8) :: v_phys(fv_nphys*fv_nphys,2,nlev)
575 :
576 : integer :: i,j,k
577 :
578 : ! Local variables
579 : real (kind=r8) :: fld_contra(np,np,2,nlev) ! vector field
580 :
581 : real (kind=r8) :: v1,v2
582 5216400 : real (kind=r8) :: D(2,2,fv_nphys*fv_nphys) ! derivative of gnomonic mapping
583 : !
584 : ! this could be done at initialization and does not need to be repeated
585 : !
586 2608200 : call setup_interpdata_for_gll_to_phys_vec_mapping(interpdata, interp_p)
587 : ! convert to contra
588 70421400 : do k=1,nlev
589 341674200 : do j=1,np
590 1424077200 : do i=1,np
591 : ! latlon->contra
592 1085011200 : fld_contra(i,j,1,k) = elem%Dinv(i,j,1,1)*v_gll(i,j,1,k) + elem%Dinv(i,j,1,2)*v_gll(i,j,2,k)
593 1356264000 : fld_contra(i,j,2,k) = elem%Dinv(i,j,2,1)*v_gll(i,j,1,k) + elem%Dinv(i,j,2,2)*v_gll(i,j,2,k)
594 : enddo
595 : enddo
596 : end do
597 :
598 70421400 : do k=1,nlev
599 680740200 : do i=1,interpdata%n_interp
600 610318800 : v_phys(i,1,k)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,1,k),interp_p,np)
601 678132000 : v_phys(i,2,k)=interpolate_2d(interpdata%interp_xy(i),fld_contra(:,:,2,k),interp_p,np)
602 : end do
603 : end do
604 26082000 : do i=1,interpdata%n_interp
605 : ! convert fld from contra->latlon
606 0 : call dmap(D(:,:,i),interpdata%interp_xy(i)%x,interpdata%interp_xy(i)%y,&
607 26082000 : elem%corners3D,cubed_sphere_map,elem%corners,elem%u2qmap,elem%facenum)
608 : end do
609 70421400 : do k=1,nlev
610 680740200 : do i=1,interpdata%n_interp
611 : ! convert fld from contra->latlon
612 610318800 : v1 = v_phys(i,1,k)
613 610318800 : v2 = v_phys(i,2,k)
614 :
615 610318800 : v_phys(i,1,k)=D(1,1,i)*v1 + D(1,2,i)*v2
616 678132000 : v_phys(i,2,k)=D(2,1,i)*v1 + D(2,2,i)*v2
617 : end do
618 : end do
619 2608200 : end function dyn2phys_vector
620 :
621 2608200 : subroutine setup_interpdata_for_gll_to_phys_vec_mapping(interpdata,interp_p)
622 : !
623 : ! initialize interpolation data structures to interpolate to phys grid
624 : ! using interpolate_mod subroutines
625 : !
626 : use interpolate_mod, only: interpolate_t, interpdata_t, interpolate_create
627 : use dimensions_mod, only : np
628 : use quadrature_mod, only : quadrature_t, gausslobatto
629 : use dimensions_mod, only : fv_nphys
630 : type (interpdata_t) , intent(out) :: interpdata
631 : type (interpolate_t) , intent(out), target :: interp_p
632 :
633 : ! local
634 : type (quadrature_t) :: gp_quadrature
635 : integer i,j,ioff,ngrid
636 : real (kind=r8) :: dx
637 :
638 2608200 : ngrid = fv_nphys*fv_nphys
639 2608200 : interpdata%n_interp=ngrid
640 : !
641 : ! initialize interpolation stuff related to basis functions
642 : !
643 2608200 : gp_quadrature = gausslobatto(np)
644 2608200 : call interpolate_create(gp_quadrature,interp_p)
645 7824600 : allocate(interpdata%interp_xy(ngrid))
646 7824600 : allocate(interpdata%ilat(ngrid) )
647 5216400 : allocate(interpdata%ilon(ngrid) )
648 : !
649 : !WARNING: THIS CODE INTERFERES WITH LAT-LON OUTPUT
650 : ! OF REGULAR SE IF nc>0
651 : !
652 2608200 : ioff=1
653 2608200 : dx = 2.0_r8/dble(fv_nphys)
654 10432800 : do j=1,fv_nphys
655 33906600 : do i=1,fv_nphys
656 23473800 : interpdata%interp_xy(ioff)%x = -1_r8+(i-0.5_r8)*dx
657 23473800 : interpdata%interp_xy(ioff)%y = -1_r8+(j-0.5_r8)*dx
658 23473800 : interpdata%ilon(ioff) = i
659 23473800 : interpdata%ilat(ioff) = j
660 31298400 : ioff=ioff+1
661 : enddo
662 : enddo
663 2608200 : end subroutine setup_interpdata_for_gll_to_phys_vec_mapping
664 :
665 :
666 26769843360 : function lagrange_1d(src_grid,src_val,ngrid,dst_point,iwidth) result(val)
667 : integer , intent(in) :: ngrid,iwidth
668 : real (kind=r8), intent(in) :: src_grid(ngrid), src_val(ngrid)
669 : real (kind=r8) :: val
670 :
671 : real (kind=r8), intent(in) :: dst_point
672 :
673 : integer :: iref, j,k
674 53539686720 : real (kind=r8) :: w(ngrid)
675 :
676 26769843360 : if (dst_point.LE.src_grid(1)) then
677 : iref=1
678 : else
679 : iref=1
680 >13743*10^7 : do while (dst_point>src_grid(iref))
681 >11066*10^7 : iref = iref + 1
682 >13743*10^7 : if (iref>ngrid) then
683 : exit
684 : end if
685 : end do
686 26767929942 : iref=iref-1
687 : end if
688 :
689 26769843360 : iref=MIN(MAX(iref,iwidth),ngrid-iwidth)
690 :
691 >24810*10^7 : w = 1.0_r8
692 >13384*10^7 : do j=iref-(iwidth-1),iref+iwidth
693 >56216*10^7 : do k=iref-(iwidth-1),iref+iwidth
694 >53539*10^7 : if (k.ne.j) then
695 >32123*10^7 : w(j)=w(j)*(dst_point-src_grid(k))/(src_grid(j)-src_grid(k))
696 : end if
697 : end do
698 : end do
699 :
700 : val=0.0_r8
701 >13384*10^7 : do j=iref-(iwidth-1),iref+iwidth
702 >13384*10^7 : val=val+w(j)*src_val(j)
703 : end do
704 26769843360 : end function lagrange_1d
705 :
706 7792200 : subroutine tensor_lagrange_interp(cubeboundary,np,nc,nhc,num_lev,nflds,psi,interp_value,llimiter,iwidth,norm_elem_coord)
707 : use control_mod, only : north, south, east, west, neast, nwest, seast, swest
708 : implicit none
709 :
710 : integer , intent(in) :: cubeboundary,nc, np, iwidth,nhc,num_lev,nflds
711 : logical , intent(in) :: llimiter(nflds) !apply limiter
712 : real (kind=r8), intent(inout) :: psi(1-nhc:nc+nhc,1-nhc:nc+nhc,num_lev,nflds) !fvm grid values with filled halo
713 : real (kind=r8), intent(out) :: interp_value(np,np,num_lev,nflds) !interpolated field
714 : real (kind=r8), intent(in) :: norm_elem_coord(2,1-nhc:nc+nhc,1-nhc:nc+nhc)
715 15584400 : integer :: which_nc_cell(np)
716 :
717 15584400 : real (kind=r8):: dx,gll_points(np)
718 15584400 : real (kind=r8):: nc_points(1-nc:nc+nc)
719 :
720 15584400 : real (kind=r8):: value(1-iwidth:nc+iwidth)
721 15584400 : real (kind=r8):: val_tmp(1-nhc:nc+nhc,1-nhc:nc+nhc)
722 :
723 15584400 : real (kind=r8):: min_value(np,np,num_lev,nflds), max_value(np,np,num_lev,nflds)
724 :
725 15584400 : integer :: imin(1-nhc:nc+nhc), imax(1-nhc:nc+nhc)
726 : integer :: k,i,j,isearch,igll,jgll,jrow,h,irow,itr
727 :
728 7792200 : gll_points(1) = -1.0_r8
729 7792200 : gll_points(2) = -sqrt(1.0_r8/5.0_r8)
730 7792200 : gll_points(3) = sqrt(1.0_r8/5.0_r8)
731 7792200 : gll_points(4) = 1.0_r8
732 :
733 7792200 : dx = 2_r8/dble(nc)
734 77922000 : do k=1-nc,2*nc
735 77922000 : nc_points(k) = -1.0_r8+dx*0.5_r8+dble(k-1)*dx
736 : end do
737 : !
738 : ! find fvm point surrounding gll points for simple limiter
739 : !
740 38961000 : do k=1,np
741 77922000 : do isearch=0,nc+1
742 77922000 : if (nc_points(isearch)<gll_points(k).and.nc_points(isearch+1).ge.gll_points(k)) exit
743 : end do
744 38961000 : which_nc_cell(k)=isearch
745 : end do
746 31168800 : do itr=1,nflds
747 31168800 : if (llimiter(itr)) then
748 : !
749 : ! fill non-existent halo cells for limiter
750 : !
751 0 : if (cubeboundary>4) then
752 0 : h=1
753 0 : select case(cubeboundary)
754 : case (nwest)
755 0 : psi(0,nc+h ,:,itr) = psi(1-h,nc ,:,itr)
756 0 : psi(1-h,nc+1,:,itr) = psi(1 ,nc+h,:,itr)
757 : case (swest)
758 0 : psi(1-h,0,:,itr) = psi(1,1-h,:,itr)
759 0 : psi(0,1-h,:,itr) = psi(1-h,1,:,itr)
760 : case (seast)
761 0 : psi(nc+h,0,:,itr) = psi(nc,1-h,:,itr)
762 0 : psi(nc+1,1-h,:,itr) = psi(nc+h,1,:,itr)
763 : case (neast)
764 0 : psi(nc+h,nc+1,:,itr) = psi(nc,nc+h,:,itr)
765 0 : psi(nc+1,nc+h,:,itr) = psi(nc+h,nc,:,itr)
766 : end select
767 : end if
768 0 : do k=1,num_lev
769 0 : do j=1,np
770 0 : do i=1,np
771 0 : max_value(i,j,k,itr) = max(&
772 0 : psi(which_nc_cell(i) ,which_nc_cell(j) ,k,itr),&
773 0 : psi(which_nc_cell(i)+1,which_nc_cell(j) ,k,itr),&
774 0 : psi(which_nc_cell(i) ,which_nc_cell(j)+1,k,itr),&
775 : psi(which_nc_cell(i)+1,which_nc_cell(j)+1,k,itr) &
776 0 : )
777 : min_value(i,j,k,itr) = min(&
778 0 : psi(which_nc_cell(i) ,which_nc_cell(j) ,k,itr),&
779 0 : psi(which_nc_cell(i)+1,which_nc_cell(j) ,k,itr),&
780 0 : psi(which_nc_cell(i) ,which_nc_cell(j)+1,k,itr),&
781 : psi(which_nc_cell(i)+1,which_nc_cell(j)+1,k,itr) &
782 0 : )
783 : end do
784 : end do
785 : end do
786 : end if
787 : end do
788 :
789 77922000 : imin=1-nhc
790 77922000 : imax=nc+nhc
791 : !
792 : ! special corner treatment
793 : !
794 7792200 : if (cubeboundary==swest) then
795 34632 : do itr=1,nflds
796 709956 : do k=1,num_lev
797 4051944 : do jrow=1,nc+iwidth
798 : !
799 : ! cubic along constant x (i=irow) in west halo to fvm points in halo
800 : !
801 10805184 : do irow=1-iwidth,0
802 27012960 : val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1:nc+nhc),psi(irow,1:nc+nhc,k,itr),nc+nhc,&
803 118181700 : norm_elem_coord(2,1,jrow),iwidth)
804 : end do
805 : end do
806 10831158 : psi(1-iwidth:0,1:nc+iwidth,k,itr) = val_tmp(1-iwidth:0,1:nc+iwidth)
807 : enddo
808 : end do
809 34632 : imin(1-nhc:0) = 1
810 : end if
811 7792200 : if (cubeboundary==nwest) then
812 34632 : do itr=1,nflds
813 709956 : do k=1,num_lev
814 4051944 : do jrow=1-iwidth,nc
815 : !
816 : ! cubic along constant x (i=irow) in west halo to fvm points in halo
817 : !
818 10805184 : do irow=1-iwidth,0
819 27012960 : val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc),psi(irow,1-nhc:nc,k,itr),nc+nhc,&
820 118181700 : norm_elem_coord(2,1,jrow),iwidth)
821 : end do
822 : end do
823 10831158 : psi(1-iwidth:0,1-iwidth:nc,k,itr) = val_tmp(1-iwidth:0,1-iwidth:nc)
824 : end do
825 : end do
826 34632 : imin(nc+1:nc+nhc) = 1
827 : end if
828 :
829 7792200 : if (cubeboundary==seast) then
830 34632 : do itr=1,nflds
831 709956 : do k=1,num_lev
832 4051944 : do jrow=1,nc+iwidth
833 27012960 : value=0.0_r8
834 : !
835 : ! cubic along constant y in ease halo to fvm points in halo
836 : !
837 10805184 : do irow=nc+1,nc+iwidth
838 27012960 : val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1:nc+nhc),psi(irow,1:nc+nhc,k,itr),nc+nhc,&
839 118181700 : norm_elem_coord(2,1,jrow),iwidth)
840 : end do
841 : end do
842 10831158 : psi(nc+1:nc+iwidth,1:nc+iwidth,k,itr) = val_tmp(nc+1:nc+iwidth,1:nc+iwidth)
843 : end do
844 : end do
845 34632 : imax(1-nhc:0) = nc
846 : end if
847 :
848 7792200 : if (cubeboundary==neast) then
849 34632 : do itr=1,nflds
850 709956 : do k=1,num_lev
851 4051944 : do jrow=1-iwidth,nc
852 : !
853 : ! cubic along constant y in ease halo to fvm points in halo
854 : !
855 10805184 : do irow=nc+1,nc+iwidth
856 27012960 : val_tmp(irow,jrow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc),psi(irow,1-nhc:nc,k,itr),nc+nhc,&
857 118181700 : norm_elem_coord(2,1,jrow),iwidth)
858 : end do
859 : end do
860 10831158 : psi(nc+1:nc+iwidth,1-iwidth:nc,k,itr) = val_tmp(nc+1:nc+iwidth,1-iwidth:nc)
861 : end do
862 : end do
863 34632 : imax(nc+1:nc+nhc) = nc
864 : end if
865 : !
866 : ! mapping
867 : !
868 : !
869 : if (cubeboundary==0.or.cubeboundary==north.or.cubeboundary==south.or.&
870 : cubeboundary==swest.or.cubeboundary==nwest.or.&
871 7792200 : cubeboundary==seast.or.cubeboundary==neast) then
872 29229408 : do itr=1,nflds
873 599202864 : do k=1,num_lev
874 2871789336 : do igll=1,np
875 : !
876 : ! cubic along constant y (j=jrow)
877 : !
878 18239150592 : do jrow=1-iwidth,nc+iwidth
879 15959256768 : value(jrow) = lagrange_1d(norm_elem_coord(1,imin(jrow):imax(jrow),jrow),psi(imin(jrow):imax(jrow),jrow,k,itr),&
880 >17776*10^7 : imax(jrow)-imin(jrow)+1,gll_points(igll),iwidth)
881 : end do
882 11969442576 : do jgll=1,np
883 45597876480 : interp_value(igll,jgll,k,itr) = lagrange_1d(norm_elem_coord(2,1,1-iwidth:nc+iwidth),value,nc+2*iwidth,&
884 >12083*10^7 : gll_points(jgll),iwidth)
885 : end do
886 : end do
887 : end do
888 : end do
889 484848 : else if (cubeboundary==east.or.cubeboundary==west) then
890 1939392 : do itr=1,nflds
891 39757536 : do k=1,num_lev
892 190545264 : do jgll=1,np
893 : !
894 : ! cubic along constant x (i=irow)
895 : !
896 1210180608 : do irow=1-iwidth,nc+iwidth
897 5294540160 : value(irow) = lagrange_1d(norm_elem_coord(2,irow,1-nhc:nc+nhc),psi(irow,1-nhc:nc+nhc,k,itr),nc+2*nhc,&
898 25565065344 : gll_points(jgll),iwidth)
899 : end do
900 794181024 : do igll=1,np
901 3025451520 : interp_value(igll,jgll,k,itr) = lagrange_1d(norm_elem_coord(1,1-iwidth:nc+iwidth,1),value,nc+2*iwidth,&
902 8017446528 : gll_points(igll),iwidth)
903 : end do
904 : end do
905 : end do
906 : end do
907 : end if
908 31168800 : do itr=1,nflds
909 31168800 : if (llimiter(itr)) then
910 0 : do k=1,num_lev
911 0 : do j=1,np
912 0 : do i=1,np
913 0 : interp_value(i,j,k,itr)=max(min_value(i,j,k,itr),min(max_value(i,j,k,itr),interp_value(i,j,k,itr)))
914 : end do
915 : enddo
916 : end do
917 : end if
918 : end do
919 7792200 : end subroutine tensor_lagrange_interp
920 :
921 0 : subroutine fvm2phys(ie,k,fvm,q_fvm,q_phys,num_trac)
922 : use dimensions_mod, only: nc,nhc,fv_nphys
923 : !
924 : ! weights must be initialized in fvm2phys_init before using these functions
925 : !
926 : type(fvm_struct) , intent(inout) :: fvm
927 : integer , intent(in) :: ie,k
928 : integer , intent(in) :: num_trac
929 :
930 : real (kind=r8), intent(inout) :: q_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,num_trac)
931 : real (kind=r8), intent(out) :: q_phys(fv_nphys,fv_nphys,num_trac)
932 :
933 : real (kind=r8) :: recons (irecons_tracer,1:nc,1:nc,1)
934 :
935 : integer :: jx,jy
936 :
937 0 : call get_dp_overlap_save(ie,k,fvm,recons)
938 0 : do jy=1,fv_nphys
939 0 : do jx=1,fv_nphys
940 0 : save_dp_phys(jx,jy,k,ie) = SUM(save_air_mass_overlap(1:save_num_overlap(jx,jy,k,ie),jx,jy,k,ie))
941 : end do
942 : end do
943 0 : call get_q_overlap_save(ie,k,fvm,q_fvm,num_trac,q_phys)
944 0 : save_dp_phys(:,:,k,ie) = save_dp_phys(:,:,k,ie)/fvm%area_sphere_physgrid
945 0 : end subroutine fvm2phys
946 :
947 :
948 0 : subroutine phys2fvm(ie,k,fvm,fq_phys,fqdp_fvm,num_trac)
949 : use dimensions_mod, only: nhc_phys,fv_nphys,nc
950 : integer , intent(in) :: ie,k
951 : type(fvm_struct) , intent(inout) :: fvm
952 : integer , intent(in) :: num_trac
953 : real (kind=r8), intent(inout) :: fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_trac)
954 : real (kind=r8), intent(out) :: fqdp_fvm (nc,nc,num_trac)
955 :
956 :
957 : integer :: h,jx,jy,jdx,jdy,m_cnst
958 :
959 0 : real(kind=r8), dimension(fv_nphys,fv_nphys) :: phys_cdp_max, phys_cdp_min
960 :
961 : integer :: num
962 : real(kind=r8) :: tmp,sum_dq_min,sum_dq_max,fq
963 :
964 0 : real(kind=r8) :: mass_phys(fv_nphys,fv_nphys)
965 : real(kind=r8) :: min_patch,max_patch,gamma
966 : real (kind=r8):: q_prev,mass_forcing,mass_forcing_phys
967 :
968 0 : real(kind=r8), allocatable, dimension(:,:,:) :: dq_min_overlap,dq_max_overlap
969 0 : real(kind=r8), allocatable, dimension(:,:,:) :: dq_overlap
970 0 : real(kind=r8), allocatable, dimension(:,:,:) :: fq_phys_overlap
971 :
972 0 : allocate(dq_min_overlap (save_max_overlap,fv_nphys,fv_nphys))
973 0 : allocate(dq_max_overlap (save_max_overlap,fv_nphys,fv_nphys))
974 0 : allocate(dq_overlap (save_max_overlap,fv_nphys,fv_nphys))
975 0 : allocate(fq_phys_overlap (save_max_overlap,fv_nphys,fv_nphys))
976 :
977 0 : do m_cnst=1,num_trac
978 0 : fqdp_fvm(:,:,m_cnst) = 0.0_r8
979 : call get_fq_overlap(ie,k,fvm,&
980 0 : fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,m_cnst),save_max_overlap,&
981 0 : fq_phys_overlap,1)
982 0 : mass_phys(1:fv_nphys,1:fv_nphys) = fq_phys(1:fv_nphys,1:fv_nphys,m_cnst)*&
983 0 : (save_dp_phys(1:fv_nphys,1:fv_nphys,k,ie)*fvm%area_sphere_physgrid)
984 :
985 0 : min_patch = MINVAL(fvm%c(0:nc+1,0:nc+1,k,m_cnst))
986 0 : max_patch = MAXVAL(fvm%c(0:nc+1,0:nc+1,k,m_cnst))
987 0 : do jy=1,fv_nphys
988 0 : do jx=1,fv_nphys
989 0 : num = save_num_overlap(jx,jy,k,ie)
990 : #ifdef debug_coupling
991 : save_q_overlap(:,jx,jy,k,m_cnst,ie) = 0.0_r8
992 : save_q_phys(jx,jy,k,m_cnst,ie) = 0.0_r8
993 : tmp = save_q_phys(jx,jy,k,m_cnst,ie)+fq_phys(jx,jy,m_cnst) !updated physics grid mixing ratio
994 : phys_cdp_max(jx,jy)= MAX(max_patch,tmp)
995 : phys_cdp_min(jx,jy)= MIN(min_patch,tmp)
996 : #else
997 0 : tmp = save_q_phys(jx,jy,k,m_cnst,ie)+fq_phys(jx,jy,m_cnst) !updated physics grid mixing ratio
998 0 : phys_cdp_max(jx,jy)= MAX(MAX(MAXVAL(save_q_overlap(1:num,jx,jy,k,m_cnst,ie)),max_patch),tmp)
999 0 : phys_cdp_min(jx,jy)= MIN(MIN(MINVAL(save_q_overlap(1:num,jx,jy,k,m_cnst,ie)),min_patch),tmp)
1000 : #endif
1001 : !
1002 : ! add high-order fq change when it does not violate monotonicity
1003 : !
1004 0 : mass_forcing_phys = 0.0_r8
1005 0 : do h=1,num
1006 0 : jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
1007 0 : q_prev = save_q_overlap(h,jx,jy,k,m_cnst,ie)
1008 : #ifndef skip_high_order_fq_map
1009 0 : save_q_overlap(h,jx,jy,k,m_cnst,ie) = save_q_overlap(h,jx,jy,k,m_cnst,ie)+fq_phys_overlap(h,jx,jy)
1010 0 : save_q_overlap(h,jx,jy,k,m_cnst,ie) = MIN(save_q_overlap(h,jx,jy,k,m_cnst,ie),phys_cdp_max(jx,jy))
1011 0 : save_q_overlap(h,jx,jy,k,m_cnst,ie) = MAX(save_q_overlap(h,jx,jy,k,m_cnst,ie),phys_cdp_min(jx,jy))
1012 0 : mass_forcing = (save_q_overlap(h,jx,jy,k,m_cnst,ie)-q_prev)*save_air_mass_overlap(h,jx,jy,k,ie)
1013 0 : mass_forcing_phys = mass_forcing_phys + mass_forcing
1014 0 : fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)+mass_forcing
1015 : #endif
1016 : !
1017 : ! prepare for mass fixing algorithm
1018 : !
1019 0 : dq_min_overlap(h,jx,jy) = save_q_overlap(h,jx,jy,k,m_cnst,ie)-phys_cdp_min(jx,jy)
1020 0 : dq_max_overlap (h,jx,jy) = save_q_overlap(h,jx,jy,k,m_cnst,ie)-phys_cdp_max(jx,jy)
1021 : end do
1022 0 : mass_phys(jx,jy) = mass_phys(jx,jy) -mass_forcing_phys
1023 : end do
1024 : end do
1025 : !
1026 : ! let physics mass tendency remove excess mass (as defined above) first proportional to how much is availabe
1027 : !
1028 : #ifdef mass_fix
1029 0 : do jy=1,fv_nphys
1030 0 : do jx=1,fv_nphys
1031 : !
1032 : ! total mass change from physics on physics grid
1033 : !
1034 0 : num = save_num_overlap(jx,jy,k,ie)
1035 0 : fq = mass_phys(jx,jy)/(fvm%area_sphere_physgrid(jx,jy)*save_dp_phys(jx,jy,k,ie))
1036 0 : if (fq<0.0_r8) then
1037 0 : sum_dq_min = SUM(dq_min_overlap(1:num,jx,jy)*save_air_mass_overlap(1:num,jx,jy,k,ie))
1038 0 : if (sum_dq_min>1.0E-14_r8) then
1039 0 : gamma=mass_phys(jx,jy)/sum_dq_min
1040 0 : do h=1,num
1041 0 : jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
1042 0 : fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)&
1043 0 : +gamma*dq_min_overlap(h,jx,jy)*save_air_mass_overlap(h,jx,jy,k,ie)
1044 : end do
1045 : end if
1046 : end if
1047 :
1048 0 : if (fq>0.0_r8) then
1049 0 : sum_dq_max = SUM(dq_max_overlap(1:num,jx,jy)*save_air_mass_overlap(1:num,jx,jy,k,ie))
1050 0 : if (sum_dq_max<-1.0E-14_r8) then
1051 0 : gamma=mass_phys(jx,jy)/sum_dq_max
1052 0 : do h=1,num
1053 0 : jdx = save_overlap_idx(1,h,jx,jy,ie); jdy = save_overlap_idx(2,h,jx,jy,ie)
1054 0 : fqdp_fvm(jdx,jdy,m_cnst) = fqdp_fvm(jdx,jdy,m_cnst)&
1055 0 : +gamma*dq_max_overlap(h,jx,jy)*save_air_mass_overlap(h,jx,jy,k,ie)
1056 : end do
1057 : end if
1058 : end if
1059 : end do
1060 : end do
1061 : #endif
1062 : !
1063 : ! convert to mass per unit area
1064 : !
1065 0 : fqdp_fvm(:,:,m_cnst) = fqdp_fvm(:,:,m_cnst)*fvm%inv_area_sphere(:,:)
1066 : end do
1067 0 : deallocate(dq_min_overlap)
1068 0 : deallocate(dq_max_overlap)
1069 0 : deallocate(fq_phys_overlap)
1070 0 : end subroutine phys2fvm
1071 :
1072 :
1073 0 : subroutine get_dp_overlap_save(ie,k,fvm,recons)
1074 : use dimensions_mod, only: nc,nhr,nhc
1075 : !
1076 : ! weights must be initialized in fvm2phys_init before using these functions
1077 : !
1078 : use dp_mapping, only: weights_all_fvm2phys, weights_eul_index_all_fvm2phys
1079 : use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
1080 : !
1081 : ! setting nhe=0 because we do not need reconstruction outside of element
1082 : !
1083 : integer, parameter :: nh = nhr!+(nhe-1) ! = 2 (nhr=2; nhe_local=1),! = 3 (nhr=2; nhe_local=2)
1084 :
1085 : type(fvm_struct) , intent(inout):: fvm
1086 : integer , intent(in) :: ie, k
1087 :
1088 : real (kind=r8) , intent(out) :: recons (irecons_tracer,nc,nc)
1089 : logical :: llimiter(1)
1090 : integer :: h,jx,jy,jdx,jdy,idx
1091 0 : llimiter=.false.
1092 0 : call get_fvm_recons(fvm,fvm%dp_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,k),recons,1,llimiter)
1093 :
1094 0 : do h=1,jall_fvm2phys(ie)
1095 0 : jx = weights_lgr_index_all_fvm2phys(h,1,ie); jy = weights_lgr_index_all_fvm2phys(h,2,ie)
1096 0 : jdx = weights_eul_index_all_fvm2phys(h,1,ie); jdy = weights_eul_index_all_fvm2phys(h,2,ie)
1097 0 : save_num_overlap(jx,jy,k,ie) = save_num_overlap(jx,jy,k,ie)+1!could be pre-computed
1098 0 : idx = save_num_overlap(jx,jy,k,ie)
1099 0 : save_overlap_idx(1,idx,jx,jy,ie) = jdx; save_overlap_idx(2,idx,jx,jy,ie) = jdy;
1100 0 : save_overlap_area(idx,jx,jy,ie) = weights_all_fvm2phys(h,1,ie)
1101 0 : save_air_mass_overlap(idx,jx,jy,k,ie) = SUM(weights_all_fvm2phys(h,:,ie)*recons(:,jdx,jdy))
1102 : #ifdef PCoM
1103 : save_air_mass_overlap(idx,jx,jy,k,ie) = fvm%dp_fvm(jdx,jdy,k)*weights_all_fvm2phys(h,1,ie)!PCoM
1104 : #endif
1105 : end do
1106 :
1107 0 : end subroutine get_dp_overlap_save
1108 :
1109 :
1110 0 : subroutine get_fq_overlap(ie,k,fvm,fq_phys,max_overlap,fq_phys_overlap,num_trac)
1111 : use dimensions_mod, only: fv_nphys, nhc_phys, nc
1112 : use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
1113 : use dp_mapping, only: weights_eul_index_all_fvm2phys
1114 : use dp_mapping, only: weights_lgr_index_all_phys2fvm, weights_eul_index_all_phys2fvm,jall_phys2fvm
1115 : use dp_mapping, only: weights_all_phys2fvm
1116 :
1117 : integer , intent(in) :: ie,k
1118 : type(fvm_struct) , intent(in) :: fvm
1119 : integer , intent(in) :: num_trac, max_overlap
1120 : real(kind=r8), dimension(max_overlap,fv_nphys,fv_nphys,num_trac),intent(out) :: fq_phys_overlap
1121 :
1122 : real (kind=r8), intent(inout) :: fq_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,num_trac)
1123 0 : real (kind=r8) :: recons_q (irecons_tracer,fv_nphys,fv_nphys,num_trac)
1124 0 : integer :: num_overlap(fv_nphys,fv_nphys)
1125 0 : logical :: llimiter_q(num_trac)
1126 : integer :: h,jx,jy,m_cnst,jdx,jdy
1127 : real (kind=r8) :: dp_tmp
1128 : integer :: idx,jxx,jyy,jdxx,jdyy,hh
1129 0 : real(kind=r8) :: weights_all_phys2fvm_local((nc+fv_nphys)**2,irecons_tracer)
1130 : !
1131 : ! could be pre-computed
1132 : !
1133 0 : do h=1,jall_fvm2phys(ie)
1134 0 : jx = weights_lgr_index_all_fvm2phys(h,1,ie)
1135 0 : jy = weights_lgr_index_all_fvm2phys(h,2,ie)
1136 0 : jdx = weights_eul_index_all_fvm2phys(h,1,ie)
1137 0 : jdy = weights_eul_index_all_fvm2phys(h,2,ie)
1138 0 : do hh=1,jall_phys2fvm(ie)
1139 0 : jxx = weights_lgr_index_all_phys2fvm(hh,1,ie)
1140 0 : jyy = weights_lgr_index_all_phys2fvm(hh,2,ie)
1141 0 : jdxx = weights_eul_index_all_phys2fvm(hh,1,ie)
1142 0 : jdyy = weights_eul_index_all_phys2fvm(hh,2,ie)
1143 0 : if (jx==jdxx.and.jy==jdyy.and.jdx==jxx.and.jdy==jyy) then
1144 0 : weights_all_phys2fvm_local(h,:) = weights_all_phys2fvm(hh,:,ie)
1145 : exit
1146 : end if
1147 : end do
1148 : end do
1149 :
1150 0 : llimiter_q=.false.
1151 0 : call get_physgrid_recons(fvm,fq_phys,recons_q,num_trac,llimiter_q)
1152 : !
1153 : ! q-dp coupling as described in equation (55) in Appendinx B of
1154 : ! Nair and Lauritzen, 2010: A Class of Deformational Flow Test Cases for Linear Transport Problems on the Sphere.
1155 : ! J. Comput. Phys.: Vol. 229, Issue 23, pp. 8868-8887, DOI:10.1016/j.jcp.2010.08.014.
1156 : !
1157 0 : num_overlap = 0
1158 0 : do h=1,jall_fvm2phys(ie)
1159 0 : jx = weights_lgr_index_all_fvm2phys(h,1,ie)
1160 0 : jy = weights_lgr_index_all_fvm2phys(h,2,ie)
1161 0 : jdx = weights_eul_index_all_fvm2phys(h,1,ie)
1162 0 : jdy = weights_eul_index_all_fvm2phys(h,2,ie)
1163 0 : num_overlap(jx,jy) = num_overlap(jx,jy)+1
1164 0 : idx = num_overlap(jx,jy)
1165 0 : dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)-fvm%dp_fvm(jdx,jdy,k)*weights_all_phys2fvm_local(h,1)
1166 0 : do m_cnst=1,num_trac
1167 0 : fq_phys_overlap(idx,jx,jy,m_cnst) = &
1168 : (fvm%dp_fvm(jdx,jdy,k)*SUM(weights_all_phys2fvm_local(h,:)*recons_q(:,jx,jy,m_cnst))+&
1169 0 : fq_phys(jx,jy,m_cnst)*dp_tmp)/save_air_mass_overlap(idx,jx,jy,k,ie)
1170 : end do
1171 : end do
1172 0 : end subroutine get_fq_overlap
1173 :
1174 0 : subroutine get_physgrid_recons(fvm,field_phys,recons_phys,num_trac,llimiter)
1175 : use dimensions_mod, only: fv_nphys,nhr_phys,nhc_phys,ns_phys
1176 : use fvm_reconstruction_mod, only: reconstruction
1177 : type(fvm_struct), intent(in) :: fvm
1178 : integer, intent(in) :: num_trac
1179 : real (kind=r8), intent(inout) :: field_phys(1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys,1,num_trac)
1180 : real (kind=r8), intent(out) :: recons_phys(irecons_tracer,1:fv_nphys,1:fv_nphys,num_trac)
1181 : logical, intent(in) :: llimiter(num_trac)
1182 :
1183 : integer, dimension(3) :: jx_min_local, jx_max_local, jy_min_local, jy_max_local
1184 :
1185 0 : jx_min_local(1) = 1 ; jx_max_local(1) = fv_nphys+1
1186 0 : jy_min_local(1) = 1 ; jy_max_local(1) = fv_nphys+1
1187 0 : jx_min_local(2) = 0 ; jx_max_local(2) = -1
1188 0 : jy_min_local(2) = 0 ; jy_max_local(2) = -1
1189 0 : jx_min_local(3) = 0 ; jx_max_local(3) = -1
1190 0 : jy_min_local(3) = 0 ; jy_max_local(3) = -1
1191 :
1192 : call reconstruction(field_phys,1,1,recons_phys,irecons_tracer,llimiter,num_trac,&
1193 : fv_nphys,0,nhr_phys,nhc_phys,nhr_phys,ns_phys,nhr_phys,&
1194 : jx_min_local,jx_max_local,jy_min_local,jy_max_local,&
1195 0 : fvm%cubeboundary,fvm%halo_interp_weight_physgrid(1:ns_phys,1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,:),&
1196 0 : fvm%ibase_physgrid(1-nhr_phys:fv_nphys+nhr_phys,1:nhr_phys,:),&
1197 0 : fvm%spherecentroid_physgrid(:,1:fv_nphys,1:fv_nphys),&
1198 0 : fvm%recons_metrics_physgrid(:,1:fv_nphys,1:fv_nphys),&
1199 0 : fvm%recons_metrics_integral_physgrid(:,1:fv_nphys,1:fv_nphys) ,&
1200 : fvm%rot_matrix_physgrid,&
1201 0 : fvm%centroid_stretch_physgrid(1:7,1:fv_nphys,1:fv_nphys),&
1202 0 : fvm%vertex_recons_weights_physgrid(:,1:irecons_tracer-1,1:fv_nphys,1:fv_nphys),&
1203 0 : fvm%vtx_cart_physgrid(:,:,1-nhc_phys:fv_nphys+nhc_phys,1-nhc_phys:fv_nphys+nhc_phys))
1204 0 : end subroutine get_physgrid_recons
1205 :
1206 0 : subroutine get_fvm_recons(fvm,field_fvm,recons_fvm,num_trac,llimiter)
1207 0 : use dimensions_mod, only: nc,nhr,nhc,ns
1208 : use fvm_reconstruction_mod, only: reconstruction
1209 :
1210 : type(fvm_struct), intent(in) :: fvm
1211 : integer, intent(in) :: num_trac
1212 : logical, intent(in) :: llimiter(num_trac)
1213 : real (kind=r8), intent(inout):: field_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,1,num_trac)
1214 : real (kind=r8), intent(out) :: recons_fvm(irecons_tracer,nc,nc,num_trac)
1215 : integer :: jx_min_local(3), jx_max_local(3), jy_min_local(3), jy_max_local(3)
1216 :
1217 0 : jx_min_local(1) = 1 ; jx_max_local(1) = nc+1
1218 0 : jy_min_local(1) = 1 ; jy_max_local(1) = nc+1
1219 0 : jx_min_local(2) = 0 ; jx_max_local(2) = -1
1220 0 : jy_min_local(2) = 0 ; jy_max_local(2) = -1
1221 0 : jx_min_local(3) = 0 ; jx_max_local(3) = -1
1222 0 : jy_min_local(3) = 0 ; jy_max_local(3) = -1
1223 :
1224 : call reconstruction(field_fvm,1,1,recons_fvm,irecons_tracer,&
1225 : llimiter,num_trac,nc,0,nhr,nhc,nhr,ns,nhr,&
1226 : jx_min_local,jx_max_local,jy_min_local,jy_max_local,&
1227 : fvm%cubeboundary,fvm%halo_interp_weight(1:ns,1-nhr:nc+nhr,1:nhr,:),fvm%ibase(1-nhr:nc+nhr,1:nhr,:),&
1228 : fvm%spherecentroid(:,1:nc,1:nc),&
1229 : fvm%recons_metrics(:,1:nc,1:nc),&
1230 : fvm%recons_metrics_integral(:,1:nc,1:nc) ,&
1231 : fvm%rot_matrix,fvm%centroid_stretch(1:7,1:nc,1:nc),&
1232 : fvm%vertex_recons_weights(:,1:irecons_tracer-1,1:nc,1:nc),&
1233 0 : fvm%vtx_cart(:,:,1-nhc:nc+nhc,1-nhc:nc+nhc))
1234 0 : end subroutine get_fvm_recons
1235 :
1236 0 : subroutine get_q_overlap_save(ie,k,fvm,q_fvm,num_trac,q_phys)
1237 0 : use dimensions_mod, only: nc,nhr,nhc,fv_nphys
1238 : !
1239 : ! weights must be initialized in fvm2phys_init before using these functions
1240 : !
1241 : use dp_mapping, only: weights_all_fvm2phys, weights_eul_index_all_fvm2phys
1242 : use dp_mapping, only: weights_lgr_index_all_fvm2phys, jall_fvm2phys
1243 : !
1244 : ! setting nhe=0 because we do not need reconstruction outside of element
1245 : !
1246 : integer, parameter :: nhe_local=0
1247 : integer, parameter :: nh = nhr!+(nhe-1) ! = 2 (nhr=2; nhe_local=1),! = 3 (nhr=2; nhe_local=2)
1248 :
1249 : type(fvm_struct) , intent(inout) :: fvm
1250 : integer , intent(in) :: ie, k
1251 : integer , intent(in) :: num_trac
1252 :
1253 : real(kind=r8), dimension(1-nhc:nc+nhc,1-nhc:nc+nhc,1:num_trac) :: q_fvm
1254 : real(kind=r8), dimension(fv_nphys,fv_nphys, num_trac), intent(out):: q_phys
1255 :
1256 0 : real (kind=r8) :: recons_q (irecons_tracer,1:nc,1:nc,num_trac)
1257 0 : logical :: llimiter_q(num_trac)
1258 : integer :: h,jx,jy,jdx,jdy,m_cnst,idx
1259 : real (kind=r8) :: dp_tmp, dp_fvm_tmp, tmp
1260 0 : real (kind=r8) :: dp_phys_inv(fv_nphys,fv_nphys)
1261 0 : integer, dimension(fv_nphys,fv_nphys):: num_overlap
1262 :
1263 0 : llimiter_q=.true.
1264 0 : call get_fvm_recons(fvm,q_fvm,recons_q,num_trac,llimiter_q)
1265 0 : num_overlap(:,:) = 0
1266 0 : q_phys = 0.0_r8
1267 0 : do h=1,jall_fvm2phys(ie)
1268 0 : jx = weights_lgr_index_all_fvm2phys(h,1,ie); jy = weights_lgr_index_all_fvm2phys(h,2,ie)
1269 0 : jdx = weights_eul_index_all_fvm2phys(h,1,ie); jdy = weights_eul_index_all_fvm2phys(h,2,ie)
1270 :
1271 0 : num_overlap(jx,jy) = num_overlap(jx,jy)+1
1272 0 : idx = num_overlap(jx,jy)
1273 :
1274 0 : dp_fvm_tmp = fvm%dp_fvm(jdx,jdy,k)
1275 0 : dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)-dp_fvm_tmp*weights_all_fvm2phys(h,1,ie)
1276 : #ifdef PCoM
1277 : dp_tmp = save_air_mass_overlap(idx,jx,jy,k,ie)
1278 : #endif
1279 0 : do m_cnst=1,num_trac
1280 0 : tmp = dp_fvm_tmp*SUM(weights_all_fvm2phys(h,:,ie)*recons_q(:,jdx,jdy,m_cnst))+q_fvm(jdx,jdy,m_cnst)*dp_tmp
1281 : #ifdef PCoM
1282 : tmp = dp_fvm_tmp*weights_all_fvm2phys(h,1,ie)*q_fvm(jdx,jdy,m_cnst)
1283 : #endif
1284 0 : save_q_overlap(idx,jx,jy,k,m_cnst,ie) = tmp/save_air_mass_overlap(idx,jx,jy,k,ie)
1285 0 : q_phys(jx,jy,m_cnst) = q_phys(jx,jy,m_cnst)+tmp
1286 : end do
1287 : end do
1288 : !
1289 : ! q_phys holds mass - convert to mixing ratio
1290 : !
1291 0 : dp_phys_inv = 1.0_r8/save_dp_phys(:,:,k,ie)!*fvm%area_sphere_physgrid)
1292 0 : do m_cnst=1,num_trac
1293 0 : q_phys(:,:,m_cnst) = q_phys(:,:,m_cnst)*dp_phys_inv
1294 0 : save_q_phys(:,:,k,m_cnst,ie) = q_phys(:,:,m_cnst)
1295 : end do
1296 0 : end subroutine get_q_overlap_save
1297 : !
1298 : ! Routine to overwrite thermodynamic active tracers on the GLL grid with CSLAM values
1299 : ! by Lagrange interpolation from 3x3 CSLAM grid to GLL grid.
1300 : !
1301 738816 : subroutine cslam2gll(elem, fvm, hybrid,nets,nete, tl_f, tl_qdp)
1302 : use dimensions_mod, only: nc,nlev,np,nhc
1303 : use hybrid_mod, only: hybrid_t
1304 : use air_composition, only: thermodynamic_active_species_num, thermodynamic_active_species_idx
1305 : use fvm_mod, only: ghostBuf_cslam2gll
1306 : use bndry_mod, only: ghost_exchange
1307 : use edge_mod, only: ghostpack,ghostunpack
1308 :
1309 : type (element_t), intent(inout):: elem(:)
1310 : type(fvm_struct), intent(inout):: fvm(:)
1311 :
1312 : type (hybrid_t), intent(in) :: hybrid ! distributed parallel structure (shared)
1313 : integer, intent(in) :: nets, nete, tl_f, tl_qdp
1314 :
1315 : integer :: ie,i,j,k,m_cnst,nq,ierr
1316 738816 : real (kind=r8), dimension(:,:,:,:,:) , allocatable :: fld_fvm, fld_gll
1317 : !
1318 : ! for tensor product Lagrange interpolation
1319 : !
1320 : integer :: nflds
1321 738816 : logical, allocatable :: llimiter(:)
1322 738816 : call t_startf('cslam2gll')
1323 738816 : nflds = thermodynamic_active_species_num
1324 :
1325 : !Allocate variables
1326 : !------------------
1327 2955264 : allocate(fld_fvm(1-nhc:nc+nhc,1-nhc:nc+nhc,nlev,nflds,nets:nete), stat=ierr)
1328 738816 : if( ierr /= 0 ) then
1329 0 : write(iulog,*) 'cslam2gll: fld_fvm allocation error = ', ierr
1330 0 : call endrun('cslam2gll: failed to allocate fld_fvm array')
1331 : end if
1332 :
1333 2955264 : allocate(fld_gll(np,np,nlev,thermodynamic_active_species_num,nets:nete),stat=ierr)
1334 738816 : if( ierr /= 0 ) then
1335 0 : write(iulog,*) 'cslam2gll: fld_gll allocation error = ', ierr
1336 0 : call endrun('cslam2gll: failed to allocate fld_gll array')
1337 : end if
1338 :
1339 2216448 : allocate(llimiter(nflds), stat=ierr)
1340 738816 : if( ierr /= 0 ) then
1341 0 : write(iulog,*) 'cslam2gll: llimiter allocation error = ', ierr
1342 0 : call endrun('cslam2gll: failed to allocate llimiter array')
1343 : end if
1344 : !------------------
1345 :
1346 3694080 : llimiter(1:nflds) = .false.
1347 5933616 : do ie=nets,nete
1348 21518016 : do m_cnst=1,thermodynamic_active_species_num
1349 425973600 : do k=1,nlev
1350 810388800 : fld_fvm(1:nc,1:nc,k,m_cnst,ie) = &
1351 6093500400 : fvm(ie)%c(1:nc,1:nc,k,thermodynamic_active_species_idx(m_cnst))
1352 : end do
1353 : end do
1354 : end do
1355 738816 : call t_startf('fvm:fill_halo_cslam2gll')
1356 5933616 : do ie=nets,nete
1357 5933616 : call ghostpack(ghostBuf_cslam2gll, fld_fvm(:,:,:,:,ie),nlev*nflds,0,ie)
1358 : end do
1359 :
1360 738816 : call ghost_exchange(hybrid,ghostBuf_cslam2gll,location='cslam2gll')
1361 :
1362 5933616 : do ie=nets,nete
1363 5933616 : call ghostunpack(ghostBuf_cslam2gll, fld_fvm(:,:,:,:,ie),nlev*nflds,0,ie)
1364 : end do
1365 738816 : call t_stopf('fvm:fill_halo_cslam2gll')
1366 : !
1367 : ! do mapping
1368 : !
1369 738816 : call fvm2dyn(fld_fvm,fld_gll,hybrid,nets,nete,nlev,nflds,fvm,llimiter,halo_filled=.true.)
1370 :
1371 5933616 : do ie=nets,nete
1372 21518016 : do m_cnst=1,thermodynamic_active_species_num
1373 31168800 : elem(ie)%state%qdp(:,:,:,m_cnst,tl_qdp) = fld_gll(:,:,:,m_cnst,ie)*&
1374 8561030400 : elem(ie)%state%dp3d(:,:,:,tl_f)
1375 : end do
1376 : end do
1377 738816 : deallocate(fld_fvm, fld_gll, llimiter)
1378 738816 : call t_stopf('cslam2gll')
1379 738816 : end subroutine cslam2gll
1380 : end module fvm_mapping
|
