Line data Source code
1 : module test_fvm_mapping
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 : use fvm_control_volume_mod, only: fvm_struct
4 : use cam_history, only: outfld
5 : use physconst, only: pi
6 : use dimensions_mod, only: np, nelemd, nlev, npsq, ntrac, use_cslam
7 : use element_mod, only: element_t
8 : implicit none
9 : private
10 :
11 : real(r8), parameter, private :: deg2rad = pi/180.0_r8
12 : real(r8), parameter, private :: psurf_moist = 100000.0_r8 !moist surface pressure
13 : integer, parameter, private :: cl_idx = 12
14 : integer, parameter, private :: cl2_idx = 13
15 : real(r8), parameter, private :: cly_constant = 4.e-6_r8
16 : integer, parameter, private :: num_fnc = 26
17 : integer, parameter, private :: offset = 15
18 :
19 : public :: test_mapping_overwrite_tendencies, test_mapping_addfld
20 : public :: test_mapping_output_mapped_tendencies, test_mapping_overwrite_dyn_state
21 : public :: test_mapping_output_phys_state
22 : contains
23 :
24 372480 : subroutine test_mapping_addfld
25 : #ifdef debug_coupling
26 : use cam_history, only: addfld, add_default, horiz_only, register_vector_field
27 : use constituents, only: cnst_get_ind,cnst_name
28 : character(LEN=128) :: name
29 : integer :: nq,m_cnst
30 :
31 : name = 'd2p_u_gll'
32 : call addfld(trim(name), (/ 'lev' /), 'I','m/2','Exact zonal wind on GLL grid',gridname='GLL')
33 : !call add_default (trim(name), 1, ' ')
34 :
35 : name = 'd2p_v_gll'
36 : call addfld(trim(name), (/ 'lev' /), 'I','m/2','Exact meridional wind on GLL grid',gridname='GLL')
37 : !call add_default (trim(name), 1, ' ')
38 :
39 : name = 'd2p_scalar_gll'
40 : call addfld(trim(name), (/ 'lev' /), 'I','','Exact scalar on GLL grid',gridname='GLL')
41 : !call add_default (trim(name), 1, ' ')
42 :
43 : name = 'd2p_u'
44 : call addfld(trim(name), (/ 'lev' /), 'I','m/2','Zonal wind mapped to physics grid')
45 : !call add_default (trim(name), 1, ' ')
46 :
47 : name = 'd2p_u_err'
48 : call addfld(trim(name), (/ 'lev' /), 'I','m/2','Error in zonal wind mapped to physics grid')
49 : !call add_default (trim(name), 1, ' ')
50 :
51 : name = 'd2p_v_err'
52 : call addfld(trim(name), (/ 'lev' /), 'I','m/2','Error in meridional wind mapped to physics grid')
53 : !call add_default (trim(name), 1, ' ')
54 :
55 : name = 'd2p_v'
56 : call addfld(trim(name), (/ 'lev' /), 'I','m/s','Meridional wind mapped to physics grid')
57 : !call add_default (trim(name), 1, ' ')
58 :
59 : name = 'd2p_scalar'
60 : call addfld(trim(name), (/ 'lev' /), 'I','','Scalar mapped to physics grid')
61 : call add_default (trim(name), 1, ' ')
62 :
63 : name = 'd2p_scalar_err'
64 : call addfld(trim(name), (/ 'lev' /), 'I','','Error in scalar mapped to physics grid')
65 : call add_default (trim(name), 1, ' ')
66 :
67 : do nq=ntrac,ntrac
68 : m_cnst = nq
69 : name = 'f2p_'//trim(cnst_name(m_cnst))//'_fvm'
70 : call addfld(trim(name), (/ 'lev' /), 'I','','Exact water tracer on fvm grid',gridname='FVM')
71 : call add_default (trim(name), 1, ' ')
72 : name = 'f2p_'//trim(cnst_name(m_cnst))//'_err'
73 : call addfld(trim(name), (/ 'lev' /), 'I','','Error in water tracer on physics grid (mapped from fvm grid)')
74 : call add_default (trim(name), 1, ' ')
75 : name = 'f2p_'//trim(cnst_name(m_cnst))//''
76 : call addfld(trim(name), (/ 'lev' /), 'I','','Water tracer on physics grid (mapped from fvm grid')
77 : call add_default (trim(name), 1, ' ')
78 : !
79 : ! physgrid to gll (condensate loading tracers)
80 : !
81 : name = 'p2d_'//trim(cnst_name(m_cnst))//''
82 : call addfld(trim(name), (/ 'lev' /), 'I','','Water tracer on physics grid')
83 : !call add_default (trim(name), 1, ' ')
84 : name = 'p2d_'//trim(cnst_name(m_cnst))//'_gll'
85 : call addfld(trim(name), (/ 'lev' /), 'I','','Water tracer on GLL grid',gridname='GLL')
86 : !call add_default (trim(name), 1, ' ')
87 : name = 'p2d_'//trim(cnst_name(m_cnst))//'_err_gll'
88 : call addfld(trim(name), (/ 'lev' /), 'I','','Error in water tracer mapped to GLL grid',gridname='GLL')
89 : !call add_default (trim(name), 1, ' ')
90 : !
91 : ! physgrid to fvm (condensate loading tracers)
92 : !
93 : name = 'p2f_'//trim(cnst_name(m_cnst))//''
94 : call addfld(trim(name), (/ 'lev' /), 'I','','Water tracer on physics grid')
95 : call add_default (trim(name), 1, ' ')
96 : name = 'p2f_'//trim(cnst_name(m_cnst))//'_fvm'
97 : call addfld(trim(name), (/ 'lev' /), 'I','','Water tracer on FVM grid',gridname='FVM')
98 : call add_default (trim(name), 1, ' ')
99 : name = 'p2f_'//trim(cnst_name(m_cnst))//'_err_fvm'
100 : call addfld(trim(name), (/ 'lev' /), 'I','','Error in water tracer mapped to FVM grid',gridname='FVM')
101 : call add_default (trim(name), 1, ' ')
102 : end do
103 : !
104 : ! temperature tendency
105 : !
106 : name = 'p2d_ptend'
107 : call addfld(trim(name), (/ 'lev' /), 'I','','T tendency on physics grid')
108 : !call add_default (trim(name), 1, ' ')
109 : name = 'p2d_ptend_gll'
110 : call addfld(trim(name), (/ 'lev' /), 'I','','T tendency on GLL grid',gridname='GLL')
111 : !call add_default (trim(name), 1, ' ')
112 : name = 'p2d_ptend_err_gll'
113 : call addfld(trim(name), (/ 'lev' /), 'I','','Error in T tendency mapped to GLL grid',gridname='GLL')
114 : !call add_default (trim(name), 1, ' ')
115 :
116 : call addfld('p2d_u', (/ 'lev' /), 'I','m/2','Zonal wind on physics grid')
117 : !call add_default ('p2d_u', 1, ' ')
118 : call addfld('p2d_v', (/ 'lev' /), 'I','m/2','Meridional wind on physics grid')
119 : !call add_default ('p2d_v', 1, ' ')
120 : call addfld('p2d_u_gll', (/ 'lev' /), 'I','m/2','Zonal wind on physics grid',gridname='GLL')
121 : !call add_default ('p2d_u_gll', 1, ' ')
122 : call addfld('p2d_v_gll', (/ 'lev' /), 'I','m/2','Meridional wind on physics grid',gridname='GLL')
123 : !call add_default ('p2d_v_gll', 1, ' ')
124 : call addfld('p2d_u_gll_err', (/ 'lev' /), 'I','m/2','Error in zonal wind interpolation to GLL grid',gridname='GLL')
125 : !call add_default ('p2d_u_gll_err', 1, ' ')
126 : call addfld('p2d_v_gll_err', (/ 'lev' /), 'I','m/2','Error in meridional wind interpolation to GLL grid',&
127 : gridname='GLL')
128 : !call add_default ('p2d_v_gll_err', 1, ' ')
129 :
130 : ! name = 'phys2dyn_'//trim(cnst_name(m_cnst))//'_physgrid'
131 : ! call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
132 : #endif
133 1536 : end subroutine test_mapping_addfld
134 :
135 23376600 : subroutine test_mapping_overwrite_tendencies(phys_state,phys_tend,ncols,lchnk,q_prev,fvm)
136 : use dimensions_mod, only: fv_nphys
137 : use constituents, only: cnst_get_ind,pcnst,cnst_name
138 : use physics_types, only: physics_state, physics_tend
139 : type(physics_state), intent(inout) :: phys_state
140 : type(physics_tend), intent(inout) :: phys_tend
141 : real(r8), dimension(:,:,:), intent(inout) :: q_prev
142 : type(fvm_struct), intent(inout):: fvm(:)
143 : integer, intent(in) :: ncols,lchnk
144 : #ifdef debug_coupling
145 : integer :: icol,k
146 : character(LEN=128) :: name
147 : integer :: m_cnst, nq, ie
148 :
149 : q_prev(:,:,ntrac) = 0.0_r8
150 : phys_state%pdel(1:ncols,:) = phys_state%pdeldry(1:ncols,:) !make sure there is no conversion from wet to dry
151 : do nq=ntrac,ntrac
152 : m_cnst = nq
153 : do icol=1,ncols
154 : do k=1,num_fnc
155 : phys_state%q(icol,k,m_cnst) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
156 : end do
157 : enddo
158 : name = 'p2f_'//trim(cnst_name(m_cnst))//''
159 : call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
160 : name = 'p2d_'//trim(cnst_name(m_cnst))//''
161 : call outfld(trim(name),phys_state%q(:ncols,:,m_cnst),ncols,lchnk)
162 : end do
163 :
164 : do icol=1,ncols
165 : do k=ntrac,ntrac
166 : phys_tend%dudt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
167 : phys_tend%dvdt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
168 : phys_tend%dtdt(icol,k) = test_func(phys_state%lat(icol), phys_state%lon(icol), k, k)
169 : end do
170 : enddo
171 : name = 'p2d_u'
172 : call outfld(trim(name),phys_tend%dudt(:ncols,:),ncols,lchnk)
173 : name = 'p2d_v'
174 : call outfld(trim(name),phys_tend%dvdt(:ncols,:),ncols,lchnk)
175 : name = 'p2d_ptend'
176 : call outfld(trim(name),phys_tend%dtdt(:ncols,:),ncols,lchnk)
177 :
178 :
179 : do icol=1,ncols
180 : do k=1,nlev
181 : ! phys_tend%dudt(icol,k) = 0.0_r8
182 : ! phys_tend%dvdt(icol,k) = 0.0_r8
183 : ! phys_tend%dtdt(icol,k) = 0.0_r8
184 : end do
185 : enddo
186 : #endif
187 23376600 : end subroutine test_mapping_overwrite_tendencies
188 :
189 369408 : subroutine test_mapping_output_mapped_tendencies(fvm,elem,nets,nete,tl_f,tl_qdp)
190 23376600 : use dimensions_mod, only: fv_nphys,nlev,nc
191 : use constituents, only: cnst_get_ind,cnst_name
192 : integer, intent(in) :: nets,nete,tl_f,tl_qdp
193 : type(fvm_struct), intent(inout):: fvm(nets:nete)
194 : type(element_t), intent(inout):: elem(nets:nete) ! pointer to dyn_out element array
195 : #ifdef debug_coupling
196 : integer :: ie,i,j,k
197 : character(LEN=128) :: name
198 : integer :: nq,m_cnst
199 : real(r8) :: diff(nc,nc,nlev,ntrac)
200 : diff = 0.0_r8
201 : do ie = nets,nete
202 : call outfld('p2d_u_gll', RESHAPE(elem(ie)%derived%fm(:,:,1,:),(/npsq,nlev/)), npsq, ie)
203 : call outfld('p2d_v_gll', RESHAPE(elem(ie)%derived%fm(:,:,2,:),(/npsq,nlev/)), npsq, ie)
204 : call outfld('p2d_ptend_gll', RESHAPE(elem(ie)%derived%ft(:,:,:),(/npsq,nlev/)), npsq, ie)
205 : do k=ntrac,ntrac
206 : do j=1,np
207 : do i=1,np
208 : elem(ie)%derived%fm(i,j,1,k) = elem(ie)%derived%fm(i,j,1,k) -&
209 : test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
210 : elem(ie)%derived%fm(i,j,2,k) = elem(ie)%derived%fm(i,j,2,k) - &
211 : test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
212 : elem(ie)%derived%ft(i,j,k) = elem(ie)%derived%ft(i,j,k) - &
213 : test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k,k)
214 : end do
215 : end do
216 : end do
217 : call outfld('p2d_u_gll_err' , RESHAPE(elem(ie)%derived%fm(:,:,1,:),(/npsq,nlev/)), npsq, ie)
218 : call outfld('p2d_v_gll_err' , RESHAPE(elem(ie)%derived%fm(:,:,2,:),(/npsq,nlev/)), npsq, ie)
219 : call outfld('p2d_ptend_err_gll', RESHAPE(elem(ie)%derived%ft(:,:,:),(/npsq,nlev/)), npsq, ie)
220 : elem(ie)%derived%ft(:,:,:) = 0.0_r8
221 : end do
222 :
223 : do ie = nets,nete
224 : do nq=ntrac,ntrac
225 : m_cnst = nq
226 : name = 'p2d_'//trim(cnst_name(m_cnst))//'_gll'
227 : call outfld(TRIM(name), RESHAPE(elem(ie)%derived%fq(:,:,:,nq),(/npsq,nlev/)), npsq, ie)
228 : ! call outfld(trim(name),&
229 : ! RESHAPE(fvm(ie)%fc(1:nc,1:nc,:,m_cnst),&
230 : ! (/nc*nc,nlev/)),nc*nc,ie)
231 : do k=1,num_fnc
232 : do j=1,np
233 : do i=1,np
234 : elem(ie)%derived%fq(i,j,k,nq) = elem(ie)%derived%fq(i,j,k,nq)-&
235 : test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
236 : end do
237 : end do
238 : end do
239 : name = 'p2d_'//trim(cnst_name(m_cnst))//'_err_gll'
240 : call outfld(TRIM(name), RESHAPE(elem(ie)%derived%fq(:,:,:,nq),(/npsq,nlev/)), npsq, ie)
241 : end do
242 : if (use_cslam) then
243 : do nq=ntrac,ntrac
244 : m_cnst = nq
245 : name = 'p2f_'//trim(cnst_name(m_cnst))//'_fvm'
246 : !
247 : ! cly
248 : !
249 : ! k=num_tracer+1
250 : ! fvm(ie)%fc(1:nc,1:nc,k,:) = fvm(ie)%fc(1:nc,1:nc,cl_idx,:)+&
251 : ! 2.0_r8*fvm(ie)%fc(1:nc,1:nc,cl2_idx,:)
252 : call outfld(trim(name),&
253 : RESHAPE(fvm(ie)%fc(1:nc,1:nc,:,m_cnst)/fvm(ie)%dp_fvm(1:nc,1:nc,:),&
254 : (/nc*nc,nlev/)),nc*nc,ie)
255 : do k=1,num_fnc
256 : do j=1,nc
257 : do i=1,nc
258 : diff(i,j,k,m_cnst) = fvm(ie)%fc(i,j,k,m_cnst)/fvm(ie)%dp_fvm(i,j,k)-&
259 : test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k, k)
260 : end do
261 : end do
262 : end do
263 : name = 'p2f_'//trim(cnst_name(m_cnst))//'_err_fvm'
264 : call outfld(TRIM(name), RESHAPE(diff(:,:,:,m_cnst),(/nc*nc,nlev/)), nc*nc, ie)
265 :
266 : end do
267 : endif
268 : end do
269 : #endif
270 369408 : end subroutine test_mapping_output_mapped_tendencies
271 :
272 370944 : subroutine test_mapping_overwrite_dyn_state(elem,fvm)
273 : use fvm_control_volume_mod, only: fvm_struct
274 : use constituents, only: cnst_name
275 : use dimensions_mod, only: nc,nhc
276 : use hybrid_mod, only: get_loop_ranges, hybrid_t,config_thread_region
277 : use control_mod, only: north, south, east, west, neast, nwest, seast, swest
278 : ! use fvm_mod, only: fill_halo_fvm_noprealloc
279 : use fvm_mod, only: fill_halo_fvm,ghostBufQnhc_h
280 : use parallel_mod, only: par
281 : type (fvm_struct), intent(inout) :: fvm(:)
282 : type(element_t), intent(inout) :: elem(:) ! pointer to dyn_out element array
283 : #ifdef debug_coupling
284 : integer :: i,j,k,ie,nq,m_cnst
285 : character(LEN=128) :: name
286 : integer :: nets,nete
287 : type(hybrid_t) :: hybrid
288 :
289 : hybrid = config_thread_region(par,'serial')
290 : call get_loop_ranges(hybrid,ibeg=nets,iend=nete)
291 : do ie=nets,nete
292 : do nq=ntrac,ntrac
293 : m_cnst = nq
294 : name = 'f2p_'//trim(cnst_name(m_cnst))//'_fvm'
295 : do k=1,num_fnc
296 : do j=1,nc
297 : do i=1,nc
298 : fvm(ie)%c(i,j,k,m_cnst) = test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k, k)
299 : end do
300 : end do
301 : end do
302 : !
303 : ! cly
304 : !
305 : ! k=num_tracer+1
306 : ! do j=1,nc
307 : ! do i=1,nc
308 : ! fvm(ie)%c(i,j,k,m_cnst) = test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k,cl_idx)+&
309 : ! 2.0_r8*test_func(fvm(ie)%center_cart(i,j)%lat,fvm(ie)%center_cart(i,j)%lon, k,cl2_idx)
310 : ! end do
311 : ! end do
312 : call outfld(TRIM(name), RESHAPE(fvm(ie)%c(1:nc,1:nc,:,m_cnst),(/nc*nc,nlev/)), nc*nc, ie)
313 : end do
314 :
315 : elem(ie)%state%Qdp(:,:,:,:,:) = 0.0_r8 !for testing the p2d map
316 : do k=1,num_fnc
317 : do j=1,np
318 : do i=1,np
319 : elem(ie)%state%v(i,j,1,k,:) = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k )
320 : elem(ie)%state%v(i,j,2,k,:) = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
321 : end do
322 : end do
323 : end do
324 : do k=1,num_fnc
325 : do j=1,np
326 : do i=1,np
327 : elem(ie)%derived%omega(i,j,k) = test_func(elem(ie)%spherep(i,j)%lat, elem(ie)%spherep(i,j)%lon, k, k)
328 : end do
329 : end do
330 : end do
331 : call outfld('d2p_scalar_gll', RESHAPE(elem(ie)%derived%omega(:,:,:) ,(/npsq,nlev/)), npsq, ie)
332 : call outfld('d2p_u_gll', RESHAPE(elem(ie)%state%v(:,:,1,:,1),(/npsq,nlev/)), npsq, ie)
333 : call outfld('d2p_v_gll', RESHAPE(elem(ie)%state%v(:,:,2,:,1),(/npsq,nlev/)), npsq, ie)
334 : end do
335 : !
336 : ! do boundary exchange (this call should be indentical to call in prim_driver)
337 : !
338 : call fill_halo_fvm(ghostBufQnhc_h,elem,fvm,hybrid,nets,nete,nhc,1,nlev,nlev)
339 : do ie=nets,nete
340 : if (fvm(ie)%cubeboundary>4) then
341 : do k=ntrac,ntrac
342 : select case(fvm(ie)%cubeboundary)
343 : case (nwest)
344 : fvm(ie)%c(0,nc+1,:,k) = fvm(ie)%c(1,nc+1,:,k)
345 : case (swest)
346 : fvm(ie)%c(0,0,:,k) = fvm(ie)%c(0,1,:,k)
347 : case (seast)
348 : fvm(ie)%c(nc+1,0,:,k) = fvm(ie)%c(0,nc,:,k)
349 : case (neast)
350 : fvm(ie)%c(nc+1,nc+1,:,k) = fvm(ie)%c(nc,nc+1,:,k)
351 : end select
352 : end do
353 : end if
354 : end do
355 : #endif
356 370944 : end subroutine test_mapping_overwrite_dyn_state
357 :
358 370944 : subroutine test_mapping_output_phys_state(phys_state,fvm)
359 370944 : use physics_types, only: physics_state
360 : use ppgrid, only: begchunk, endchunk, pver, pcols
361 : use constituents, only: cnst_get_ind,cnst_name
362 : type(physics_state), intent(inout) :: phys_state(begchunk:endchunk)
363 : type(fvm_struct), pointer:: fvm(:)
364 : #ifdef debug_coupling
365 : integer :: lchnk, ncol,k,icol,m_cnst,nq,ie
366 : character(LEN=128) :: name
367 :
368 : do lchnk = begchunk, endchunk
369 : call outfld('d2p_scalar', phys_state(lchnk)%omega(1:pcols,1:pver), pcols, lchnk)
370 : call outfld('d2p_u', phys_state(lchnk)%U(1:pcols,1:pver), pcols, lchnk)
371 : call outfld('d2p_v', phys_state(lchnk)%V(1:pcols,1:pver), pcols, lchnk)
372 : if (use_cslam) then
373 : do nq=ntrac,ntrac
374 : m_cnst = nq
375 : name = 'f2p_'//trim(cnst_name(m_cnst))
376 : !
377 : ! cly
378 : !
379 : !phys_state(lchnk)%q(1:pcols,num_tracer+1,m_cnst)=phys_state(lchnk)%q(1:pcols,cl_idx,m_cnst)+&
380 : ! 2.0_r8*phys_state(lchnk)%q(1:pcols,12,m_cnst)
381 : call outfld(TRIM(name), phys_state(lchnk)%q(1:pcols,1:pver,m_cnst), pcols, lchnk)
382 : ! k=num_tracer+1
383 : ! do icol=1,phys_state(lchnk)%ncol
384 : ! phys_state(lchnk)%q(icol,k,m_cnst) = phys_state(lchnk)%q(icol,cl_idx,m_cnst)+&
385 : ! 2.0_r8*phys_state(lchnk)%q(icol,cl2_idx,m_cnst)-&
386 : ! cly_constant
387 : ! end do
388 : do k=1,num_fnc
389 : do icol=1,phys_state(lchnk)%ncol
390 : phys_state(lchnk)%q(icol,k,m_cnst) = phys_state(lchnk)%q(icol,k,m_cnst)&
391 : -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,k)
392 : end do
393 : enddo
394 : name = 'f2p_'//trim(cnst_name(m_cnst))//'_err'
395 : call outfld(TRIM(name), phys_state(lchnk)%q(1:pcols,1:pver,m_cnst), pcols, lchnk)
396 : phys_state(lchnk)%q(1:pcols,1:pver,m_cnst) = 0.0_r8
397 : end do
398 : end if
399 : end do
400 :
401 :
402 : do lchnk = begchunk, endchunk
403 : do k=1,nlev
404 : do icol=1,phys_state(lchnk)%ncol
405 : phys_state(lchnk)%U(icol,k) = phys_state(lchnk)%U(icol,k)&
406 : -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k, 9)
407 : phys_state(lchnk)%V(icol,k) = phys_state(lchnk)%V(icol,k)&
408 : -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,10)
409 : end do
410 : enddo
411 : name = 'd2p_u_err'
412 : call outfld(trim(name),phys_state(lchnk)%U(:pcols,:),pcols,lchnk)
413 : name = 'd2p_v_err'
414 : call outfld(trim(name),phys_state(lchnk)%V(:pcols,:),pcols,lchnk)
415 : do k=1,num_fnc
416 : do icol=1,phys_state(lchnk)%ncol
417 : phys_state(lchnk)%omega(icol,k) = phys_state(lchnk)%omega(icol,k)&
418 : -test_func(phys_state(lchnk)%lat(icol), phys_state(lchnk)%lon(icol), k,k)
419 : end do
420 : end do
421 : name = 'd2p_scalar_err'
422 : call outfld(trim(name),phys_state(lchnk)%omega(:pcols,:),pcols,lchnk)
423 : end do
424 : #endif
425 370944 : end subroutine test_mapping_output_phys_state
426 :
427 : ! subroutine test_mapping_overwrite_state(phys_tend,nets,nete)
428 : !#ifdef debug_coupling
429 : ! phys_tend(lchnk)%dtdt(icol,ilyr) = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,9)
430 : ! phys_tend(lchnk)%dudt(icol,ilyr) = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,12)
431 : ! phys_tend(lchnk)%dvdt(icol,ilyr) = 0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,13)
432 : ! q_prev(icol, ilyr, 2:pcnst, lchnk) = 0.0D0
433 : ! do m=2,pcnst
434 : ! phys_state(lchnk)%Q(icol,ilyr,m)=0!test_func(phys_state(lchnk)%lat(icol),phys_state(lchnk)%lon(icol),ilyr,m)
435 : ! end do
436 : !#endif
437 : ! end subroutine test_mapping_overwrite_state
438 :
439 : #ifdef debug_coupling
440 : function test_func(lat_in, lon_in, k, funcnum) result(fout)
441 : use hycoef, only: hyai, hybi, hyam, hybm, ps0
442 : use shr_sys_mod, only: shr_sys_flush
443 : use cam_abortutils, only: endrun
444 : real(r8), intent(in) :: lon_in
445 : real(r8), intent(in) :: lat_in
446 : integer, intent(in) :: k
447 : integer, intent(in) :: funcnum
448 : real(r8) :: fout
449 : real(r8) :: lon1,lat1,R0,Rg1,Rg2,lon2,lat2,cl,cl2
450 : real(r8) :: eta_c
451 :
452 : real(r8) :: radius = 10.0_r8 ! radius of the perturbation
453 : real(r8) :: perturb_lon = 20.0_r8 ! longitudinal position, 20E
454 : real(r8) :: perturb_lat = 40.0_r8 ! latitudinal position, 40N
455 : real(r8) :: cos_tmp, sin_tmp, eta
456 : real(r8) :: u_wind, v_wind, lat, lon, u_tmp, v_tmp
457 : real(r8) :: rotation_angle
458 : real(r8) :: det,r,k1,k2
459 : real(r8), parameter :: pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164_r8
460 : real(r8), parameter :: half_pi = pi*0.5_r8
461 : real(r8), parameter :: degrees_to_radians = pi/180.0_r8
462 : real(r8), parameter :: k1_lat_center = 20.d0*degrees_to_radians
463 : real(r8), parameter :: k1_lon_center = 300.d0*degrees_to_radians
464 :
465 : lon = lon_in
466 : lat = lat_in
467 :
468 :
469 : select case(MOD(funcnum,8)+1)
470 : case(1)
471 : !
472 : ! Non-smooth scalar field (slotted cylinder)
473 : !
474 : R0 = 0.5_r8
475 : lon1 = 5.0_r8 * PI / 6.0_r8
476 : lat1 = 0.0_r8
477 : Rg1 = acos(sin(lat1)*sin(lat)+cos(lat1)*cos(lat)*cos(lon-lon1))
478 : lon2 = 7.0_r8 * PI / 6.0_r8
479 : lat2 = 0.0_r8
480 : Rg2 = acos(sin(lat2)*sin(lat)+cos(lat2)*cos(lat)*cos(lon-lon2))
481 :
482 : if ((Rg1 <= R0) .AND. (abs(lon-lon1) >= R0/6)) then
483 : fout = 2.0_r8
484 : elseif ((Rg2 <= R0) .AND. (abs(lon-lon2) >= R0/6)) then
485 : fout = 2.0_r8
486 : elseif ((Rg1 <= R0) .AND. (abs(lon-lon1) < R0/6) &
487 : .AND. (lat-lat1 < -5.0_r8*R0/12.0_r8)) then
488 : fout = 2.0_r8
489 : elseif ((Rg2 <= R0) .AND. (abs(lon-lon2) < R0/6) &
490 : .AND. (lat-lat2 > 5.0_r8*R0/12.0_r8)) then
491 : fout = 2.0_r8
492 : else
493 : fout = 1.0_r8
494 : endif
495 : case(2)
496 : !
497 : ! Smooth Gaussian "ball"
498 : !
499 : R0 = 10.0_r8 ! radius of the perturbation
500 : lon1 = 20.0_r8*deg2rad ! longitudinal position, 20E
501 : lat1 = 40.0_r8 *deg2rad ! latitudinal position, 40N
502 : eta_c = 0.6_r8
503 : sin_tmp = SIN(lat1)*SIN(lat)
504 : cos_tmp = COS(lat1)*COS(lat)
505 : Rg1 = ACOS( sin_tmp + cos_tmp*COS(lon-lon1) ) ! great circle distance
506 : eta = (hyam(k)*ps0 + hybm(k)*psurf_moist)/psurf_moist
507 : fout = EXP(- ((Rg1*R0)**2 + ((eta-eta_c)/0.1_r8)**2))
508 : if (ABS(fout) < 1.0E-8_r8) then
509 : fout = 0.0_r8
510 : end IF
511 : case(3)
512 : !
513 : !
514 : !
515 : fout = 0.5_r8 * ( tanh( 3.0_r8*abs(lat)-pi ) + 1.0_r8)
516 : case(4)
517 : fout = 2.0_r8+cos(5.0_r8+40*lon)!1.0e-8_r8
518 : fout = -0.5_r8-0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
519 : case(5)
520 : !
521 : ! approximately Y^2_2 spherical harmonic
522 : !
523 : fout = sin(lon)*cos(40*lat)!1.0e-8_r8
524 : fout = 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
525 : case(6)
526 : !
527 : ! approximately Y32_16 spherical harmonic
528 : !
529 : fout = 0.5_r8 + 0.5_r8*(cos(16*lon)*(sin(2_r8*lat)**16))
530 : case(7)
531 : fout = 2.0_r8 + lat
532 : case(8)
533 : fout = 2.0_r8 + cos(lon)
534 : case(9)
535 : rotation_angle = 45.0_r8*pi/180.0_r8
536 : CALL regrot(lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,1)
537 : call Rossby_Haurwitz (lon, lat,u_wind, v_wind)
538 : CALL turnwi(u_wind,v_wind,u_tmp,v_tmp,lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,-1)
539 : fout = u_tmp
540 : case(10)
541 : rotation_angle = 45.0_r8*pi/180.0_r8
542 : CALL regrot(lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,1)
543 : call Rossby_Haurwitz (lon, lat,u_wind, v_wind)
544 : CALL turnwi(u_wind,v_wind,u_tmp,v_tmp,lon_in,lat_in,lon,lat,0.0_r8,-0.5_r8*pi+rotation_angle,-1)
545 : fout = v_tmp
546 : case(11)
547 : fout = 1.0E-8_r8
548 : case(12)
549 : !
550 : ! Terminator chemistry initial condition
551 : !
552 : k1 = 1.0_r8*max(0.d0,sin(lat)*sin(k1_lat_center) + cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
553 : k2 = 1._r8
554 :
555 : r = k1 / (4._r8*k2)
556 : det = sqrt(r*r + 2._r8*cly_constant*r)
557 :
558 : fout = (det-r)
559 : ! fout = cly_constant/2._r8 - (det-r)/2._r8
560 : case(13)
561 : !
562 : ! Terminator chemistry initial condition
563 : !
564 : k1 = 1.0_r8*max(0.d0,sin(lat)*sin(k1_lat_center) + cos(lat)*cos(k1_lat_center)*cos(lon-k1_lon_center))
565 : k2 = 1._r8
566 :
567 : r = k1 / (4._r8*k2)
568 : det = sqrt(r*r + 2._r8*cly_constant*r)
569 :
570 : ! fout = (det-r)
571 : fout = cly_constant/2._r8 - (det-r)/2._r8
572 : case default
573 : ! call endrun("Illegal funcnum_arg in test_func")
574 : fout = 1.0_r8
575 : end select
576 : end function test_func
577 :
578 : function test_wind(lat, lon, iwind) result(fout)
579 : use cam_abortutils, only: endrun
580 : real(r8), intent(in) :: lon
581 : real(r8), intent(in) :: lat
582 : integer, intent(in) :: iwind
583 :
584 : real(r8) :: fout
585 :
586 :
587 : fout = 0
588 : end function test_wind
589 :
590 :
591 : SUBROUTINE regrot(pxreg,pyreg,pxrot,pyrot,pxcen,pycen,kcall)
592 : use physconst, only: pi
593 : !
594 : !----------------------------------------------------------------------
595 : !
596 : !* conversion between regular and rotated spherical coordinates.
597 : !*
598 : !* pxreg longitudes of the regular coordinates
599 : !* pyreg latitudes of the regular coordinates
600 : !* pxrot longitudes of the rotated coordinates
601 : !* pyrot latitudes of the rotated coordinates
602 : !* all coordinates given in degrees n (negative for s)
603 : !* and degrees e (negative values for w)
604 : !* pxcen regular longitude of the south pole of the rotated grid
605 : !* pycen regular latitude of the south pole of the rotated grid
606 : !*
607 : !* kcall=-1: find regular as functions of rotated coordinates.
608 : !* kcall= 1: find rotated as functions of regular coordinates.
609 : !
610 : !-----------------------------------------------------------------------
611 : !
612 : integer kxdim,kydim,kx,ky,kcall
613 : real(r8) :: pxreg,pyreg,&
614 : pxrot,pyrot,&
615 : pxcen,pycen
616 : !
617 : !-----------------------------------------------------------------------
618 : !
619 : real(r8) zsycen,zcycen,zxmxc,zsxmxc,zcxmxc,zsyreg,zcyreg, &
620 : zsyrot,zcyrot,zcxrot,zsxrot,zpi,zpih
621 : integer jy,jx
622 :
623 : zpih = pi*0.5_r8
624 : !
625 : !----------------------------------------------------------------------
626 : !
627 : zsycen = SIN((pycen+zpih))
628 : zcycen = COS((pycen+zpih))
629 : !
630 : IF (kcall.eq.1) then
631 : !
632 : zxmxc = pxreg - pxcen
633 : zsxmxc = SIN(zxmxc)
634 : zcxmxc = COS(zxmxc)
635 : zsyreg = SIN(pyreg)
636 : zcyreg = COS(pyreg)
637 : zsyrot = zcycen*zsyreg - zsycen*zcyreg*zcxmxc
638 : zsyrot = max(zsyrot,-1.0_r8)
639 : zsyrot = min(zsyrot,+1.0_r8)
640 : !
641 : pyrot = ASIN(zsyrot)
642 : !
643 : zcyrot = COS(pyrot)
644 : zcxrot = (zcycen*zcyreg*zcxmxc +zsycen*zsyreg)/zcyrot
645 : zcxrot = max(zcxrot,-1.0_r8)
646 : zcxrot = min(zcxrot,+1.0_r8)
647 : zsxrot = zcyreg*zsxmxc/zcyrot
648 : !
649 : pxrot = ACOS(zcxrot)
650 : !
651 : IF (zsxrot < 0.0_r8) then
652 : pxrot = -pxrot
653 : end IF
654 : !
655 : ELSEIF (kcall.eq.-1) then
656 : !
657 : zsxrot = SIN(pxrot)
658 : zcxrot = COS(pxrot)
659 : zsyrot = SIN(pyrot)
660 : zcyrot = COS(pyrot)
661 : zsyreg = zcycen*zsyrot + zsycen*zcyrot*zcxrot
662 : zsyreg = max(zsyreg,-1.0_r8)
663 : zsyreg = min(zsyreg,+1.0_r8)
664 : !
665 : pyreg = ASIN(zsyreg)
666 : !
667 : zcyreg = COS(pyreg)
668 : zcxmxc = (zcycen*zcyrot*zcxrot -&
669 : zsycen*zsyrot)/zcyreg
670 : zcxmxc = max(zcxmxc,-1.0_r8)
671 : zcxmxc = min(zcxmxc,+1.0_r8)
672 : zsxmxc = zcyrot*zsxrot/zcyreg
673 : zxmxc = ACOS(zcxmxc)
674 : IF (zsxmxc < 0.0_r8) then
675 : zxmxc = -zxmxc
676 : end IF
677 : !
678 : pxreg = zxmxc + pxcen
679 : !
680 : ELSE
681 : WRITE(6,'(1x,''invalid kcall in regrot'')')
682 : STOP
683 : ENDIF
684 : END SUBROUTINE regrot
685 :
686 : SUBROUTINE turnwi(puarg,pvarg,pures,pvres,pxreg,pyreg,pxrot,pyrot,pxcen,pycen,kcall)
687 : use physconst, only: pi
688 : !
689 : !-----------------------------------------------------------------------
690 : !
691 : !* turn horizontal velocity components between regular and
692 : !* rotated spherical coordinates.
693 : !
694 : !* puarg : input u components
695 : !* pvarg : input v components
696 : !* pures : output u components
697 : !* pvres : output v components
698 : !* pa : transformation coefficients
699 : !* pb : -"-
700 : !* pc : -"-
701 : !* pd : -"-
702 : !* pxreg : regular longitudes
703 : !* pyreg : regular latitudes
704 : !* pxrot : rotated longitudes
705 : !* pyrot : rotated latitudes
706 : !* kxdim : dimension in the x (longitude) direction
707 : !* kydim : dimension in the y (latitude) direction
708 : !* kx : number of gridpoints in the x direction
709 : !* ky : number of gridpoints in the y direction
710 : !* pxcen : regular longitude of the south pole of the
711 : !* transformed grid
712 : !* pycen : regular latitude of the south pole of the
713 : !* transformed grid
714 : !*
715 : !* kcall < 0 : find wind components in regular coordinates
716 : !* from wind components in rotated coordinates
717 : !* kcall > 0 : find wind components in rotated coordinates
718 : !* from wind components in regular coordinates
719 : !* note that all coordinates are given in degrees n and degrees e.
720 : !* (negative values for s and w)
721 : !
722 : !-----------------------------------------------------------------------
723 :
724 : integer kxdim,kydim,kx,ky,kcall
725 : real(r8) puarg,pvarg, &
726 : pures,pvres, &
727 : pa, pb, &
728 : pc, pd, &
729 : pxreg,pyreg, &
730 : pxrot,pyrot
731 : real(r8) pxcen,pycen
732 : !
733 : !-----------------------------------------------------------------------
734 : !
735 : integer jy,jx
736 : real(r8) zpih,zsyc,zcyc,zsxreg,zcxreg,zsyreg,zcyreg,zxmxc,&
737 : zsxmxc,zcxmxc,zsxrot,zcxrot,zsyrot,zcyrot
738 : !
739 : !-----------------------------------------------------------------------
740 : !
741 : IF (kcall.eq.1) then
742 : zpih = pi*0.5_r8
743 : zsyc = SIN(pycen+zpih)
744 : zcyc = COS(pycen+zpih)
745 : !
746 : zsxreg = SIN(pxreg)
747 : zcxreg = COS(pxreg)
748 : zsyreg = SIN(pyreg)
749 : zcyreg = COS(pyreg)
750 : !
751 : zxmxc = pxreg - pxcen
752 : zsxmxc = SIN(zxmxc)
753 : zcxmxc = COS(zxmxc)
754 : !
755 : zsxrot = SIN(pxrot)
756 : zcxrot = COS(pxrot)
757 : zsyrot = SIN(pyrot)
758 : zcyrot = COS(pyrot)
759 : !
760 : pa = zcyc*zsxmxc*zsxrot + zcxmxc*zcxrot
761 : pb = zcyc*zcxmxc*zsyreg*zsxrot - zsyc*zcyreg*zsxrot - &
762 : zsxmxc*zsyreg*zcxrot
763 : pc = zsyc*zsxmxc/zcyrot
764 : pd = (zsyc*zcxmxc*zsyreg + zcyc*zcyreg)/zcyrot
765 : !
766 : pures = pa*puarg + pb*pvarg
767 : pvres = pc*puarg + pd*pvarg
768 : ELSEIF (kcall.eq.-1) then
769 : zpih = pi*0.5_r8
770 : zsyc = SIN(pycen+zpih)
771 : zcyc = COS(pycen+zpih)
772 : !
773 : zsxreg = SIN(pxreg)
774 : zcxreg = COS(pxreg)
775 : zsyreg = SIN(pyreg)
776 : zcyreg = COS(pyreg)
777 : !
778 : zxmxc = pxreg - pxcen
779 : zsxmxc = SIN(zxmxc)
780 : zcxmxc = COS(zxmxc)
781 : !
782 : zsxrot = SIN(pxrot)
783 : zcxrot = COS(pxrot)
784 : zsyrot = SIN(pyrot)
785 : zcyrot = COS(pyrot)
786 : !
787 : pa = zcxmxc*zcxrot + zcyc*zsxmxc*zsxrot
788 : pb = zcyc*zsxmxc*zcxrot*zsyrot + zsyc*zsxmxc*zcyrot -&
789 : zcxmxc*zsxrot*zsyrot
790 : pc =-zsyc*zsxrot/zcyreg
791 : pd = (zcyc*zcyrot - zsyc*zcxrot*zsyrot)/zcyreg
792 : !
793 : pures = pa*puarg + pb*pvarg
794 : pvres = pc*puarg + pd*pvarg
795 : ELSE
796 : write(6,'(1x,''invalid kcall in turnwi'')')
797 : STOP
798 : ENDIF
799 : END SUBROUTINE turnwi
800 :
801 : SUBROUTINE Rossby_Haurwitz (lon, lat,u_wind, v_wind)
802 : use physconst, only: rearth
803 : !-----------------------------------------------------------------------
804 : ! input parameters
805 : !-----------------------------------------------------------------------
806 : real(r8), intent(in) :: lon, & ! longitude in radians
807 : lat ! latitude in radians
808 : ! both coefficients 'a' and 'b' are needed at the full model level
809 : !-----------------------------------------------------------------------
810 : ! input parameters
811 : !-----------------------------------------------------------------------
812 : real(r8), intent(out) :: u_wind, & ! zonal wind in m/s
813 : v_wind ! meridional wind in m/s
814 :
815 : !-----------------------------------------------------------------------
816 : ! test case parameters
817 : !-----------------------------------------------------------------------
818 : real(r8),parameter :: u0 = 50._r8, & ! reference wind
819 : n = 4._r8 ! wavenumber
820 :
821 : !-----------------------------------------------------------------------
822 : ! local
823 : !-----------------------------------------------------------------------
824 : real(r8) :: tmp1, tmp2, tmp3, KK, MM
825 : real(r8) :: sin_lat, cos_lat, sin_slat, cos_slat
826 :
827 : !-----------------------------------------------------------------------
828 : ! initialize the wind components
829 : !-----------------------------------------------------------------------
830 : MM = u0/(n*rearth) ! parameter M
831 : KK = u0/(n*rearth) ! parameter K
832 :
833 :
834 : cos_lat = cos(lat)
835 : sin_lat = sin(lat)
836 : tmp1 = rearth * MM * cos_lat
837 : tmp2 = rearth * KK * cos_lat**(n-1._r8)*(n*sin_lat**2 - cos_lat**2)
838 : tmp3 = -rearth * KK * n * cos_lat**(n-1._r8) * sin_lat
839 : u_wind = tmp1 + tmp2 * cos(n*lon)
840 : v_wind = tmp3 * sin(n*lon)
841 : end subroutine Rossby_Haurwitz
842 :
843 : #endif
844 : end module test_fvm_mapping
|
