FV3-JEDI
wind_variables_mod.f90
Go to the documentation of this file.
1 ! (C) Copyright 2018-2019 UCAR
2 !
3 ! This software is licensed under the terms of the Apache Licence Version 2.0
4 ! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
5 
6 module wind_vt_mod
7 
8 use fv3jedi_constants_mod, only: pi, rad2deg
9 use fv3jedi_geom_mod, only: fv3jedi_geom
10 use fv3jedi_kinds_mod, only: kind_real
11 use fv3jedi_communication_mod, only: gather_field, scatter_field
12 
13 use fv_mp_adm_mod, only: mpp_update_domains_adm
14 use mpp_domains_mod, only: mpp_update_domains, dgrid_ne
15 use mpp_domains_mod, only: mpp_get_boundary, mpp_get_boundary_ad
16 use mpp_parameter_mod, only: cgrid_ne
17 
18 use femps_fv3_mod, only: fv3field_to_ufield, ufield_to_fv3field
19 use femps_grid_mod, only: fempsgrid
20 use femps_operators_mod, only: fempsoprs
21 use femps_solve_mod, only: laplace, inverselaplace
22 
23 implicit none
24 private
25 
26 public sfc_10m_winds
27 
28 public psichi_to_uava
29 public psichi_to_uava_adm
30 
31 public psichi_to_udvd
32 public udvd_to_psichi
33 
34 public psichi_to_vortdivg
35 
36 public a2d
37 public a2d_ad
38 public d2a
39 public d2a_ad
40 
41 contains
42 
43 !----------------------------------------------------------------------------
44 ! Lowest model level winds to 10m speed and direction -----------------------
45 !----------------------------------------------------------------------------
46 
47 subroutine sfc_10m_winds(geom,usrf,vsrf,f10r,spd10m,dir10m)
48 
49  implicit none
50 
51  !Arguments
52  type(fv3jedi_geom) , intent(in ) :: geom !Geometry for the model
53  real(kind=kind_real), intent(in ) :: usrf(geom%isc:geom%iec,geom%jsc:geom%jec) !Lowest model level u m/s
54  real(kind=kind_real), intent(in ) :: vsrf(geom%isc:geom%iec,geom%jsc:geom%jec) !Lowest model level v m/s
55  real(kind=kind_real), intent(in ) :: f10r(geom%isc:geom%iec,geom%jsc:geom%jec) !Ratio of lowest level to 10m
56  real(kind=kind_real), intent(out) :: spd10m(geom%isc:geom%iec,geom%jsc:geom%jec) !10m wind speed u m/s
57  real(kind=kind_real), intent(out) :: dir10m(geom%isc:geom%iec,geom%jsc:geom%jec) !10m model wind direction
58 
59  !Locals
60  integer :: isc,iec,jsc,jec,i,j
61  integer :: iquad
62  real(kind=kind_real) :: windangle, windratio
63  real(kind=kind_real), parameter :: windscale = 999999.0_kind_real
64  real(kind=kind_real), parameter :: windlimit = 0.0001_kind_real
65  real(kind=kind_real), parameter :: quadcof(4,2) = reshape((/ 0.0_kind_real, 1.0_kind_real, &
66  1.0_kind_real, 2.0_kind_real, &
67  1.0_kind_real, -1.0_kind_real, &
68  1.0_kind_real, -1.0_kind_real /), (/4, 2/))
69 
70  !In GSI these calculations are done after interpolation to obs location
71 
72  isc = geom%isc
73  iec = geom%iec
74  jsc = geom%jsc
75  jec = geom%jec
76 
77  !10m wind speed
78  spd10m(isc:iec,jsc:jec) = f10r(isc:iec,jsc:jec)*sqrt( usrf(isc:iec,jsc:jec)*usrf(isc:iec,jsc:jec) + &
79  vsrf(isc:iec,jsc:jec)*vsrf(isc:iec,jsc:jec) )
80 
81  !10m wind direction
82  do j = jsc,jec
83  do i = isc,iec
84 
85  if (usrf(i,j)*f10r(i,j) >= 0.0_kind_real .and. vsrf(i,j)*f10r(i,j) >= 0.0_kind_real) iquad = 1
86  if (usrf(i,j)*f10r(i,j) >= 0.0_kind_real .and. vsrf(i,j)*f10r(i,j) < 0.0_kind_real) iquad = 2
87  if (usrf(i,j)*f10r(i,j) < 0.0_kind_real .and. vsrf(i,j)*f10r(i,j) >= 0.0_kind_real) iquad = 3
88  if (usrf(i,j)*f10r(i,j) < 0.0_kind_real .and. vsrf(i,j)*f10r(i,j) < 0.0_kind_real) iquad = 4
89 
90  if (abs(vsrf(i,j)*f10r(i,j)) >= windlimit) then
91  windratio = (usrf(i,j)*f10r(i,j)) / (vsrf(i,j)*f10r(i,j))
92  else
93  windratio = 0.0_kind_real
94  if (abs(usrf(i,j)*f10r(i,j)) > windlimit) then
95  windratio = windscale * usrf(i,j)*f10r(i,j)
96  endif
97  endif
98 
99  windangle = atan(abs(windratio))
100 
101  dir10m(i,j) = rad2deg*(quadcof(iquad,1) * pi + windangle * quadcof(iquad, 2))
102 
103  enddo
104  enddo
105 
106 end subroutine sfc_10m_winds
107 
108 !----------------------------------------------------------------------------
109 
110 subroutine psichi_to_uava(geom,psi,chi,ua,va)
111 
112  implicit none
113  type(fv3jedi_geom), intent(inout) :: geom
114  real(kind=kind_real), intent(inout) :: psi(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz) !Stream function
115  real(kind=kind_real), intent(inout) :: chi(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz) !Velocity potential
116  real(kind=kind_real), intent(out) :: ua(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Agrid winds (u)
117  real(kind=kind_real), intent(out) :: va(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Agrid winds (v)
118 
119  integer :: i,j,k
120 
121  ! x-----------------x
122  ! | |
123  ! | |
124  ! | |
125  ! | x |
126  ! | psi,chi |
127  ! | ua,va |
128  ! | |
129  ! | |
130  ! x-----------------x
131 
132  call mpp_update_domains(psi, geom%domain, complete=.true.)
133  call mpp_update_domains(chi, geom%domain, complete=.true.)
134 
135  do k=1,geom%npz
136  do j=geom%jsc,geom%jec
137  do i=geom%isc,geom%iec
138 
139  ua(i,j,k) = (psi(i,j+1,k) - psi(i,j-1,k))/(geom%dyc(i,j) + geom%dyc(i,j+1)) + &
140  (chi(i+1,j,k) - chi(i-1,j,k))/(geom%dxc(i,j) + geom%dxc(i+1,j))
141  va(i,j,k) = -(psi(i+1,j,k) - psi(i-1,j,k))/(geom%dxc(i,j) + geom%dxc(i+1,j)) + &
142  (chi(i,j+1,k) - chi(i,j-1,k))/(geom%dyc(i,j) + geom%dyc(i,j+1))
143 
144  enddo
145  enddo
146  enddo
147 
148 end subroutine psichi_to_uava
149 
150 !----------------------------------------------------------------------------
151 
152 subroutine psichi_to_udvd(geom,psi,chi,u,v)
153 
154  implicit none
155  type(fv3jedi_geom), intent(inout) :: geom
156  real(kind=kind_real), intent(in) :: psi(geom%isc:geom%iec ,geom%jsc:geom%jec ,1:geom%npz) !Stream function
157  real(kind=kind_real), intent(in) :: chi(geom%isc:geom%iec ,geom%jsc:geom%jec ,1:geom%npz) !Velocity potential
158  real(kind=kind_real), intent(out) :: u(geom%isc:geom%iec ,geom%jsc:geom%jec+1,1:geom%npz) !Dgrid winds (u)
159  real(kind=kind_real), intent(out) :: v(geom%isc:geom%iec+1,geom%jsc:geom%jec ,1:geom%npz) !Dgrid winds (v)
160 
161  integer :: i,j,k
162  real(kind=kind_real), allocatable, dimension(:,:,:) :: psid, chid
163  real(kind=kind_real) :: u_rot(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz) !Rotational winds (u)
164  real(kind=kind_real) :: v_rot(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz) !Rotational winds (v)
165  real(kind=kind_real) :: u_div(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz) !Divergent winds (u)
166  real(kind=kind_real) :: v_div(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz) !Divergent winds (v)
167 
168  real(kind=kind_real) :: uc_div(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz) !Divergent winds (u)
169  real(kind=kind_real) :: vc_div(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz) !Divergent winds (v)
170  real(kind=kind_real) :: ua_div(geom%isd:geom%ied ,geom%jsd:geom%jed ,1:geom%npz) !Divergent winds (u)
171  real(kind=kind_real) :: va_div(geom%isd:geom%ied ,geom%jsd:geom%jed ,1:geom%npz) !Divergent winds (v)
172 
173  ! Fill halos on psi and chi
174  !--------------------------
175  allocate(psid(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz))
176  allocate(chid(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz))
177  psid = 0.0_kind_real
178  chid = 0.0_kind_real
179 
180  psid(geom%isc:geom%iec,geom%jsc:geom%jec,:) = psi(geom%isc:geom%iec,geom%jsc:geom%jec,:)
181  chid(geom%isc:geom%iec,geom%jsc:geom%jec,:) = chi(geom%isc:geom%iec,geom%jsc:geom%jec,:)
182 
183  call mpp_update_domains(psid, geom%domain, complete=.true.)
184  call mpp_update_domains(chid, geom%domain, complete=.true.)
185 
186 
187  ! Compute rotational and divergent winds
188  ! --------------------------------------
189 
190  do k=1,geom%npz
191  do j=geom%jsc,geom%jec+1
192  do i=geom%isc,geom%iec
193  u_rot(i,j,k) = -(psid(i,j,k) - psid(i,j-1,k))/geom%dyc(i,j)
194  enddo
195  enddo
196  enddo
197 
198  do k=1,geom%npz
199  do j=geom%jsc,geom%jec
200  do i=geom%isc,geom%iec+1
201  v_rot(i,j,k) = (psid(i,j,k) - psid(i-1,j,k))/geom%dxc(i,j)
202  enddo
203  enddo
204  enddo
205 
206  do k=1,geom%npz
207  do j=geom%jsc,geom%jec
208  do i=geom%isc,geom%iec+1
209  uc_div(i,j,k) = (chid(i,j,k) - chid(i-1,j,k))/geom%dxc(i,j)
210  enddo
211  enddo
212  enddo
213 
214  do k=1,geom%npz
215  do j=geom%jsc,geom%jec+1
216  do i=geom%isc,geom%iec
217  vc_div(i,j,k) = (chid(i,j,k) - chid(i,j-1,k))/geom%dyc(i,j)
218  enddo
219  enddo
220  enddo
221 
222  ! Convert C-grid divergent winds to A-grid divergent winds
223  ! --------------------------------------------------------
224  call fill_cgrid_winds(geom,uc_div,vc_div,fillhalo=.true.)
225 
226  do k = 1,geom%npz
227  call ctoa(geom, uc_div(:,:,k), vc_div(:,:,k), ua_div(:,:,k), va_div(:,:,k))
228  enddo
229 
230  ! Convert A-grid divergent winds to D-grid divergent winds
231  ! --------------------------------------------------------
232 
233  do k = 1,geom%npz
234  call atod(geom, ua_div(:,:,k), va_div(:,:,k), u_div(:,:,k), v_div(:,:,k))
235  enddo
236 
237 
238  ! Sum rotational and divergent parts
239  ! ----------------------------------
240 
241  call fill_dgrid_winds(geom,u_div,v_div,fillhalo=.true.)
242 
243  do k=1,geom%npz
244  do j=geom%jsc,geom%jec+1
245  do i=geom%isc,geom%iec
246  u(i,j,k) = u_rot(i,j,k) + u_div(i,j,k)
247  enddo
248  enddo
249  enddo
250 
251  do k=1,geom%npz
252  do j=geom%jsc,geom%jec
253  do i=geom%isc,geom%iec+1
254  v(i,j,k) = v_rot(i,j,k) + v_div(i,j,k)
255  enddo
256  enddo
257  enddo
258 
259 end subroutine psichi_to_udvd
260 
261 !----------------------------------------------------------------------------
262 
263 subroutine psichi_to_uava_adm(geom,psi_ad,chi_ad,ua_ad,va_ad)
264 
265  implicit none
266  type(fv3jedi_geom), intent(inout) :: geom
267  real(kind=kind_real), intent(inout) :: psi_ad(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz) !Stream function
268  real(kind=kind_real), intent(inout) :: chi_ad(geom%isd:geom%ied,geom%jsd:geom%jed,1:geom%npz) !Velocity potential
269  real(kind=kind_real), intent(inout) :: ua_ad(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Agrid winds (u)
270  real(kind=kind_real), intent(inout) :: va_ad(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Agrid winds (v)
271 
272  integer :: i,j,k
273  real(kind=kind_real) :: temp1
274  real(kind=kind_real) :: temp2
275  real(kind=kind_real) :: temp3
276  real(kind=kind_real) :: temp4
277 
278  ! x-----------------x
279  ! | |
280  ! | |
281  ! | |
282  ! | x |
283  ! | psi,chi |
284  ! | ua,va |
285  ! | |
286  ! | |
287  ! x-----------------x
288 
289  do k=geom%npz,1,-1
290  do j=geom%jec,geom%jsc,-1
291  do i=geom%iec,geom%isc,-1
292 
293  temp1 = va_ad(i,j,k)/(geom%dyc(i,j)+geom%dyc(i,j+1))
294  temp2 = -(va_ad(i,j,k)/(geom%dxc(i,j)+geom%dxc(i+1,j)))
295  chi_ad(i,j+1,k) = chi_ad(i,j+1,k) + temp1
296  chi_ad(i,j-1,k) = chi_ad(i,j-1,k) - temp1
297  psi_ad(i+1,j,k) = psi_ad(i+1,j,k) + temp2
298  psi_ad(i-1,j,k) = psi_ad(i-1,j,k) - temp2
299  va_ad(i,j,k) = 0.0_8
300 
301  temp3 = ua_ad(i,j,k)/(geom%dyc(i,j)+geom%dyc(i,j+1))
302  temp4 = ua_ad(i,j,k)/(geom%dxc(i,j)+geom%dxc(i+1,j))
303  psi_ad(i,j+1,k) = psi_ad(i,j+1,k) + temp3
304  psi_ad(i,j-1,k) = psi_ad(i,j-1,k) - temp3
305  chi_ad(i+1,j,k) = chi_ad(i+1,j,k) + temp4
306  chi_ad(i-1,j,k) = chi_ad(i-1,j,k) - temp4
307  ua_ad(i,j,k) = 0.0_8
308 
309  end do
310  end do
311  end do
312 
313  call mpp_update_domains_adm(ua_ad, chi_ad, geom%domain, complete=.true.)
314  call mpp_update_domains_adm(va_ad, psi_ad, geom%domain, complete=.true.)
315 
316 end subroutine psichi_to_uava_adm
317 
318 !----------------------------------------------------------------------------
319 
320 subroutine udvd_to_psichi(geom,grid,oprs,u_in,v_in,psi,chi,lsize,lev_start,lev_final,vor_out,div_out)
321 
322  implicit none
323  type(fv3jedi_geom), intent(inout) :: geom
324  type(fempsgrid), intent(in) :: grid
325  type(fempsoprs), intent(in) :: oprs
326  real(kind=kind_real), intent(inout) :: u_in(geom%isc:geom%iec ,geom%jsc:geom%jec+1,1:geom%npz) !Dgrid winds (u)
327  real(kind=kind_real), intent(inout) :: v_in(geom%isc:geom%iec+1,geom%jsc:geom%jec ,1:geom%npz) !Dgrid winds (v)
328  real(kind=kind_real), intent(out) :: psi(geom%isc:geom%iec ,geom%jsc:geom%jec ,1:geom%npz) !Stream function
329  real(kind=kind_real), intent(out) :: chi(geom%isc:geom%iec ,geom%jsc:geom%jec ,1:geom%npz) !Velocity potential
330  integer, intent(in) :: lsize
331  integer, intent(in) :: lev_start(lsize),lev_final(lsize)
332  real(kind=kind_real), optional, intent(out) :: vor_out(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Stream function
333  real(kind=kind_real), optional, intent(out) :: div_out(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Velocity potential
334 
335  integer :: i,j,k
336  real(kind=kind_real), allocatable, dimension(:,:,:) :: u, v
337  real(kind=kind_real), allocatable, dimension(:,:,:) :: uc, vc
338  real(kind=kind_real), allocatable, dimension(:,:,:) :: ut, vt
339 
340  real(kind=kind_real), allocatable, dimension(:,:,:) :: vor, div
341 
342  real(kind=kind_real), allocatable, dimension(:,:,:) :: vorgcomm, divgcomm
343  real(kind=kind_real), allocatable, dimension(:,:,:) :: psigcomm, chigcomm
344  real(kind=kind_real), allocatable, dimension(:,:,:,:) :: vorg, divg !Global level of vor and div
345  real(kind=kind_real), allocatable, dimension(:,:,:,:) :: psig, chig !Global level of psi and chi
346 
347  real(kind=kind_real), allocatable, dimension(:) :: voru, divu !Unstructured vor and div
348  real(kind=kind_real), allocatable, dimension(:) :: psiu, chiu !Unstructured psi and chi
349 
350  integer :: ranki, n
351  integer, allocatable :: lev_proc(:)
352 
353  ! ------------------------------------------ !
354  ! Convert D-grid winds to A-grid psi and chi !
355  ! ------------------------------------------ !
356 
357 
358  ! Fill edge of D grid winds
359  ! -------------------------
360  allocate(u(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz))
361  allocate(v(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz))
362 
363  ! Copy internal part
364  u(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) = u_in(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)
365  v(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) = v_in(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz)
366 
367  call fill_dgrid_winds(geom,u,v,fillhalo=.true.)
368 
369 
370  ! Get C grid winds
371  ! ----------------
372  allocate(uc(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz))
373  allocate(vc(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz))
374 
375  do k = 1,geom%npz
376  call udvd_to_ucvc(geom,u(:,:,k),v(:,:,k),uc(:,:,k),vc(:,:,k))
377  enddo
378 
379  call fill_cgrid_winds(geom,uc,vc,fillhalo=.true.)
380 
381 
382  ! Compute ut,vt
383  ! -------------
384  allocate(ut(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz))
385  allocate(vt(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz))
386 
387  do k = 1,geom%npz
388 
389  call compute_utvt(geom, uc(:,:,k), vc(:,:,k), ut(:,:,k), vt(:,:,k), 1.0_kind_real)
390 
391  do j=geom%jsc,geom%jec
392  do i=geom%isc,geom%iec+1
393  if ( ut(i,j,k) > 0. ) then
394  ut(i,j,k) = geom%dy(i,j)*ut(i,j,k)*geom%sin_sg(i-1,j,3)
395  else
396  ut(i,j,k) = geom%dy(i,j)*ut(i,j,k)*geom%sin_sg(i,j,1)
397  endif
398  enddo
399  enddo
400  do j=geom%jsc,geom%jec+1
401  do i=geom%isc,geom%iec
402  if ( vt(i,j,k) > 0. ) then
403  vt(i,j,k) = geom%dx(i,j)*vt(i,j,k)*geom%sin_sg(i,j-1,4)
404  else
405  vt(i,j,k) = geom%dx(i,j)*vt(i,j,k)*geom%sin_sg(i,j,2)
406  endif
407  enddo
408  enddo
409 
410  enddo
411 
412  !D/C grid u and v to A grid vorticity and divergence
413  !---------------------------------------------------
414  allocate(vor(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz))
415  allocate(div(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz))
416  vor = 0.0_kind_real
417  div = 0.0_kind_real
418 
419  do k=1,geom%npz
420  do j=geom%jsc,geom%jec
421  do i=geom%isc,geom%iec
422  vor(i,j,k) = geom%rarea(i,j)*( u(i,j,k)*geom%dx(i,j)-u(i,j+1,k)*geom%dx(i,j+1) - &
423  v(i,j,k)*geom%dy(i,j)+v(i+1,j,k)*geom%dy(i+1,j))
424  div(i,j,k) = geom%rarea(i,j)*( ut(i+1,j,k)-ut(i,j,k) + vt(i,j+1,k)-vt(i,j,k) )
425  !geom%rarea(i,j)*( uc(i+1,j,k)*geom%dy(i+1,j)-uc(i,j,k)*geom%dy(i,j) + &
426  ! vc(i,j+1,k)*geom%dx(i,j+1)-vc(i,j,k)*geom%dx(i,j+1) )
427  enddo
428  enddo
429  enddo
430 
431  if (present(vor_out)) vor_out = vor
432  if (present(div_out)) div_out = div
433 
434  ! Gather voricity and divergence to one processor and compute psi and chi
435  ! -----------------------------------------------------------------------
436 
437  ranki = geom%f_comm%rank() + 1
438 
439  allocate(vorgcomm(1:geom%npx-1,1:geom%npy-1,6))
440  allocate(divgcomm(1:geom%npx-1,1:geom%npy-1,6))
441 
442  if (ranki <= lsize) then
443  allocate(vorg(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
444  allocate(divg(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
445  allocate(psig(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
446  allocate(chig(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
447 
448  allocate(voru(grid%nface(grid%ngrids)))
449  allocate(divu(grid%nface(grid%ngrids)))
450  allocate(psiu(grid%nface(grid%ngrids)))
451  allocate(chiu(grid%nface(grid%ngrids)))
452  endif
453 
454  ! Processor that each level is associated with
455  allocate(lev_proc(geom%npz))
456  do k = 1,geom%npz
457  do n = 1,lsize
458  if (k >= lev_start(n) .and. k <= lev_final(n)) then
459  lev_proc(k) = n
460  endif
461  enddo
462  enddo
463 
464  ! Gather field to respective processor
465  do k=1,geom%npz
466  call gather_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,vor(:,:,k),vorgcomm)
467  call gather_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,div(:,:,k),divgcomm)
468  if (ranki == lev_proc(k)) then
469  vorg(:,:,:,k) = vorgcomm
470  divg(:,:,:,k) = divgcomm
471  endif
472  enddo
473 
474  deallocate(vorgcomm,divgcomm)
475 
476  ! Loop over level and compute psi/chi
477  if (ranki <= lsize) then ! Only processors with a level
478  do k = lev_start(ranki),lev_final(ranki)
479 
480  if (geom%f_comm%rank()==0) &
481  print*, "Doing femps level ", k, ". Proc range:",lev_start(ranki),"-",lev_final(ranki)
482 
483  call fv3field_to_ufield(grid,geom%npx-1,vorg(:,:,:,k),voru)
484  call fv3field_to_ufield(grid,geom%npx-1,divg(:,:,:,k),divu)
485 
486  ! Convert to area integrals, required by femps
487  voru = voru*grid%farea(:,grid%ngrids)
488  divu = divu*grid%farea(:,grid%ngrids)
489 
490  ! Solve poisson equation (\psi=\nabla^{-2}\zeta, \chi=\nabla^{-2}D)
491  call inverselaplace(grid,oprs,grid%ngrids,voru,psiu,level=k)
492  call inverselaplace(grid,oprs,grid%ngrids,divu,chiu,level=k)
493 
494  ! Convert from area integrals
495  psiu = psiu/grid%farea(:,grid%ngrids)
496  chiu = chiu/grid%farea(:,grid%ngrids)
497 
498  call ufield_to_fv3field(grid,geom%npx-1,psiu,psig(:,:,:,k))
499  call ufield_to_fv3field(grid,geom%npx-1,chiu,chig(:,:,:,k))
500 
501  enddo
502  endif
503 
504  ! Scatter field from respective processor
505  allocate(psigcomm(1:geom%npx-1,1:geom%npy-1,6))
506  allocate(chigcomm(1:geom%npx-1,1:geom%npy-1,6))
507 
508  do k=1,geom%npz
509  if (ranki == lev_proc(k)) then
510  psigcomm = psig(:,:,:,k)
511  chigcomm = chig(:,:,:,k)
512  endif
513  call scatter_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,psigcomm,psi(:,:,k))
514  call scatter_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,chigcomm,chi(:,:,k))
515  enddo
516 
517 end subroutine udvd_to_psichi
518 
519 ! ------------------------------------------------------------------------------
520 
521 subroutine psichi_to_vortdivg(geom,grid,oprs,psi,chi,lsize,lev_start,lev_final,vor,div)
522 
523  implicit none
524  type(fv3jedi_geom), intent(inout) :: geom
525  type(fempsgrid), intent(in) :: grid
526  type(fempsoprs), intent(in) :: oprs
527  real(kind=kind_real), intent(in) :: psi(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Stream function
528  real(kind=kind_real), intent(in) :: chi(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Velocity potential
529  integer, intent(in) :: lsize
530  integer, intent(in) :: lev_start(lsize),lev_final(lsize)
531  real(kind=kind_real), optional, intent(out) :: vor(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Vorticity
532  real(kind=kind_real), optional, intent(out) :: div(geom%isc:geom%iec,geom%jsc:geom%jec,1:geom%npz) !Divergence
533 
534  integer :: k
535  real(kind=kind_real), allocatable, dimension(:,:,:) :: vorgcomm, divgcomm
536  real(kind=kind_real), allocatable, dimension(:,:,:) :: psigcomm, chigcomm
537  real(kind=kind_real), allocatable, dimension(:,:,:,:) :: vorg, divg !Global level of vor and div
538  real(kind=kind_real), allocatable, dimension(:,:,:,:) :: psig, chig !Global level of psi and chi
539  real(kind=kind_real), allocatable, dimension(:) :: voru, divu !Unstructured vor and div
540  real(kind=kind_real), allocatable, dimension(:) :: psiu, chiu !Unstructured psi and chi
541 
542  integer :: ranki, n
543  integer, allocatable :: lev_proc(:)
544 
545  ! Gather voricity and divergence to one processor and compute psi and chi
546  ! -----------------------------------------------------------------------
547 
548  ranki = geom%f_comm%rank() + 1
549 
550  if (ranki <= lsize) then
551  allocate(vorg(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
552  allocate(divg(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
553  allocate(psig(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
554  allocate(chig(1:geom%npx-1,1:geom%npy-1,6,lev_start(ranki):lev_final(ranki)))
555  allocate(voru(grid%nface(grid%ngrids)))
556  allocate(divu(grid%nface(grid%ngrids)))
557  allocate(psiu(grid%nface(grid%ngrids)))
558  allocate(chiu(grid%nface(grid%ngrids)))
559  endif
560 
561  ! Processor that each level is associated with
562  allocate(lev_proc(geom%npz))
563  do k = 1,geom%npz
564  do n = 1,lsize
565  if (k >= lev_start(n) .and. k <= lev_final(n)) then
566  lev_proc(k) = n
567  endif
568  enddo
569  enddo
570 
571  allocate(psigcomm(1:geom%npx-1,1:geom%npy-1,6))
572  allocate(chigcomm(1:geom%npx-1,1:geom%npy-1,6))
573  do k=1,geom%npz
574  call gather_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,psi(:,:,k),psigcomm)
575  call gather_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,chi(:,:,k),chigcomm)
576  if (ranki == lev_proc(k)) then
577  psig(:,:,:,k) = psigcomm
578  chig(:,:,:,k) = chigcomm
579  endif
580  enddo
581  deallocate(psigcomm,chigcomm)
582 
583  ! Loop over level and compute vort/divg
584  if (ranki <= lsize) then ! Only processors with a level
585  do k = lev_start(ranki),lev_final(ranki)
586 
587  if (geom%f_comm%rank()==0) &
588  print*, "Doing femps level ", k, ". Proc range:",lev_start(ranki),"-",lev_final(ranki)
589 
590  call fv3field_to_ufield(grid,geom%npx-1,psig(:,:,:,k),psiu)
591  call fv3field_to_ufield(grid,geom%npx-1,chig(:,:,:,k),chiu)
592 
593  ! Convert to area integrals, required by femps
594  psiu = psiu*grid%farea(:,grid%ngrids)
595  chiu = chiu*grid%farea(:,grid%ngrids)
596 
597  ! Solve poisson equation (\psi=\nabla^{-2}\zeta, \chi=\nabla^{-2}D)
598  call laplace(grid,oprs,grid%ngrids,psiu,voru)
599  call laplace(grid,oprs,grid%ngrids,chiu,divu)
600 
601  ! Convert from area integrals
602  voru = voru/grid%farea(:,grid%ngrids)
603  divu = divu/grid%farea(:,grid%ngrids)
604 
605  call ufield_to_fv3field(grid,geom%npx-1,voru,vorg(:,:,:,k))
606  call ufield_to_fv3field(grid,geom%npx-1,divu,divg(:,:,:,k))
607 
608  enddo
609  endif
610 
611  ! Scatter field from root to all
612  allocate(vorgcomm(1:geom%npx-1,1:geom%npy-1,6))
613  allocate(divgcomm(1:geom%npx-1,1:geom%npy-1,6))
614  do k=1,geom%npz
615  if (ranki == lev_proc(k)) then
616  vorgcomm = vorg(:,:,:,k)
617  divgcomm = divg(:,:,:,k)
618  endif
619  call scatter_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,vorgcomm,vor(:,:,k))
620  call scatter_field(geom,geom%f_comm%communicator(),lev_proc(k)-1,divgcomm,div(:,:,k))
621  enddo
622  deallocate(vorgcomm,divgcomm)
623 
624 end subroutine psichi_to_vortdivg
625 
626 ! ------------------------------------------------------------------------------
627 
628 subroutine udvd_to_ucvc(geom,u,v,uc,vc,ua_out,va_out)
629 
630 implicit none
631 type(fv3jedi_geom), target, intent(inout) :: geom
632 real(kind=kind_real), intent(in) :: u(geom%isd:geom%ied ,geom%jsd:geom%jed+1) !Dgrid winds
633 real(kind=kind_real), intent(in) :: v(geom%isd:geom%ied+1,geom%jsd:geom%jed ) !Dgrid winds
634 real(kind=kind_real), intent(out) :: uc(geom%isd:geom%ied+1,geom%jsd:geom%jed ) !Cgrid winds
635 real(kind=kind_real), intent(out) :: vc(geom%isd:geom%ied ,geom%jsd:geom%jed+1) !Cgrid winds
636 real(kind=kind_real), optional, intent(out) :: ua_out(geom%isd:geom%ied ,geom%jsd:geom%jed) !Agrid winds
637 real(kind=kind_real), optional, intent(out) :: va_out(geom%isd:geom%ied ,geom%jsd:geom%jed) !Agrid winds
638 
639 ! Normally input/output
640 real(kind=kind_real) :: ut(geom%isd:geom%ied ,geom%jsd:geom%jed )
641 real(kind=kind_real) :: vt(geom%isd:geom%ied ,geom%jsd:geom%jed )
642 real(kind=kind_real) :: utmp(geom%isd:geom%ied ,geom%jsd:geom%jed )
643 real(kind=kind_real) :: vtmp(geom%isd:geom%ied ,geom%jsd:geom%jed )
644 real(kind=kind_real) :: ua(geom%isd:geom%ied ,geom%jsd:geom%jed )
645 real(kind=kind_real) :: va(geom%isd:geom%ied ,geom%jsd:geom%jed )
646 logical :: dord4
647 logical :: nested
648 integer :: grid_type
649 integer :: npx, npy
650 
651 ! Locals
652 real(kind=kind_real), parameter:: a1 = 0.5625_kind_real
653 real(kind=kind_real), parameter:: a2 = -0.0625_kind_real
654 real(kind=kind_real), parameter:: c1 = -2.0_kind_real/14.0_kind_real
655 real(kind=kind_real), parameter:: c2 = 11.0_kind_real/14.0_kind_real
656 real(kind=kind_real), parameter:: c3 = 5.0_kind_real/14.0_kind_real
657 
658 integer npt, i, j, ifirst, ilast, id
659 integer :: is, ie, js, je
660 integer :: isd, ied, jsd, jed
661 real(kind=kind_real), pointer, dimension(:,:,:) :: sin_sg
662 real(kind=kind_real), pointer, dimension(:,:) :: cosa_u, cosa_v, cosa_s
663 real(kind=kind_real), pointer, dimension(:,:) :: rsin_u, rsin_v, rsin2
664 real(kind=kind_real), pointer, dimension(:,:) :: dxa,dya
665 
666 is = geom%isc
667 ie = geom%iec
668 js = geom%jsc
669 je = geom%jec
670 isd = geom%isd
671 ied = geom%ied
672 jsd = geom%jsd
673 jed = geom%jed
674 npx = geom%npx
675 npy = geom%npy
676 
677 dord4 = geom%dord4
678 nested = geom%nested
679 grid_type = geom%grid_type
680 
681 sin_sg => geom%sin_sg
682 cosa_u => geom%cosa_u
683 cosa_v => geom%cosa_v
684 cosa_s => geom%cosa_s
685 rsin_u => geom%rsin_u
686 rsin_v => geom%rsin_v
687 rsin2 => geom%rsin2
688 dxa => geom%dxa
689 dya => geom%dya
690 
691 if ( geom%dord4 ) then
692  id = 1
693 else
694  id = 0
695 endif
696 
697 if (grid_type < 3 .and. .not. nested) then
698  npt = 4
699 else
700  npt = -2
701 endif
702 
703 ! Initialize the non-existing corner regions
704 utmp(:,:) = 1.0e30_kind_real
705 vtmp(:,:) = 1.0e30_kind_real
706 
707 if ( nested) then
708 
709  do j=jsd+1,jed-1
710  do i=isd,ied
711  utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
712  enddo
713  enddo
714 
715  do i=isd,ied
716  j = jsd
717  utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
718  j = jed
719  utmp(i,j) = 0.5*(u(i,j)+u(i,j+1))
720  end do
721 
722  do j=jsd,jed
723  do i=isd+1,ied-1
724  vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
725  enddo
726  i = isd
727  vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
728  i = ied
729  vtmp(i,j) = 0.5*(v(i,j)+v(i+1,j))
730  enddo
731 
732  do j=jsd,jed
733  do i=isd,ied
734  ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
735  va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
736  enddo
737  enddo
738 
739 else
740 
741  !----------
742  ! Interior:
743  !----------
744  do j=max(npt,js-1),min(npy-npt,je+1)
745  do i=max(npt,isd),min(npx-npt,ied)
746  utmp(i,j) = a2*(u(i,j-1)+u(i,j+2)) + a1*(u(i,j)+u(i,j+1))
747  enddo
748  enddo
749  do j=max(npt,jsd),min(npy-npt,jed)
750  do i=max(npt,is-1),min(npx-npt,ie+1)
751  vtmp(i,j) = a2*(v(i-1,j)+v(i+2,j)) + a1*(v(i,j)+v(i+1,j))
752  enddo
753  enddo
754 
755  !----------
756  ! edges:
757  !----------
758  if (grid_type < 3) then
759 
760  if ( js==1 .or. jsd<npt) then
761  do j=jsd,npt-1
762  do i=isd,ied
763  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
764  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
765  enddo
766  enddo
767  endif
768 
769  if ( (je+1)==npy .or. jed>=(npy-npt)) then
770  do j=npy-npt+1,jed
771  do i=isd,ied
772  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
773  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
774  enddo
775  enddo
776  endif
777 
778  if ( is==1 .or. isd<npt ) then
779  do j=max(npt,jsd),min(npy-npt,jed)
780  do i=isd,npt-1
781  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
782  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
783  enddo
784  enddo
785  endif
786  if ( (ie+1)==npx .or. ied>=(npx-npt)) then
787  do j=max(npt,jsd),min(npy-npt,jed)
788  do i=npx-npt+1,ied
789  utmp(i,j) = 0.5*(u(i,j) + u(i,j+1))
790  vtmp(i,j) = 0.5*(v(i,j) + v(i+1,j))
791  enddo
792  enddo
793  endif
794 
795  endif
796 
797  ! Contra-variant components at cell center:
798  do j=js-1-id,je+1+id
799  do i=is-1-id,ie+1+id
800  ua(i,j) = (utmp(i,j)-vtmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
801  va(i,j) = (vtmp(i,j)-utmp(i,j)*cosa_s(i,j)) * rsin2(i,j)
802  enddo
803  enddo
804 
805 end if
806 
807 ! A -> C
808 !--------------
809 ! Fix the edges
810 !--------------
811 ! Xdir:
812 if( geom%sw_corner ) then
813  do i=-2,0
814  utmp(i,0) = -vtmp(0,1-i)
815  enddo
816 endif
817 if( geom%se_corner ) then
818  do i=0,2
819  utmp(npx+i,0) = vtmp(npx,i+1)
820  enddo
821 endif
822 if( geom%ne_corner ) then
823  do i=0,2
824  utmp(npx+i,npy) = -vtmp(npx,je-i)
825  enddo
826 endif
827 if( geom%nw_corner ) then
828  do i=-2,0
829  utmp(i,npy) = vtmp(0,je+i)
830  enddo
831 endif
832 
833 if (grid_type < 3 .and. .not. nested) then
834  ifirst = max(3, is-1)
835  ilast = min(npx-2,ie+2)
836 else
837  ifirst = is-1
838  ilast = ie+2
839 endif
840 
841 !---------------------------------------------
842 ! 4th order interpolation for interior points:
843 !---------------------------------------------
844 do j=js-1,je+1
845  do i=ifirst,ilast
846  uc(i,j) = a2*(utmp(i-2,j)+utmp(i+1,j)) + a1*(utmp(i-1,j)+utmp(i,j))
847  ut(i,j) = (uc(i,j) - v(i,j)*cosa_u(i,j))*rsin_u(i,j)
848  enddo
849 enddo
850 
851 if (grid_type < 3) then
852  ! Xdir:
853  if( geom%sw_corner ) then
854  ua(-1,0) = -va(0,2)
855  ua( 0,0) = -va(0,1)
856  endif
857  if( geom%se_corner ) then
858  ua(npx, 0) = va(npx,1)
859  ua(npx+1,0) = va(npx,2)
860  endif
861  if( geom%ne_corner ) then
862  ua(npx, npy) = -va(npx,npy-1)
863  ua(npx+1,npy) = -va(npx,npy-2)
864  endif
865  if( geom%nw_corner ) then
866  ua(-1,npy) = va(0,npy-2)
867  ua( 0,npy) = va(0,npy-1)
868  endif
869 
870  if( is==1 .and. .not. nested ) then
871  do j=js-1,je+1
872  uc(0,j) = c1*utmp(-2,j) + c2*utmp(-1,j) + c3*utmp(0,j)
873  ut(1,j) = edge_interpolate4(ua(-1:2,j), dxa(-1:2,j))
874  !Want to use the UPSTREAM value
875  if (ut(1,j) > 0.) then
876  uc(1,j) = ut(1,j)*sin_sg(0,j,3)
877  else
878  uc(1,j) = ut(1,j)*sin_sg(1,j,1)
879  end if
880  uc(2,j) = c1*utmp(3,j) + c2*utmp(2,j) + c3*utmp(1,j)
881  ut(0,j) = (uc(0,j) - v(0,j)*cosa_u(0,j))*rsin_u(0,j)
882  ut(2,j) = (uc(2,j) - v(2,j)*cosa_u(2,j))*rsin_u(2,j)
883  enddo
884  endif
885 
886  if( (ie+1)==npx .and. .not. nested ) then
887  do j=js-1,je+1
888  uc(npx-1,j) = c1*utmp(npx-3,j)+c2*utmp(npx-2,j)+c3*utmp(npx-1,j)
889  ut(npx, j) = edge_interpolate4(ua(npx-2:npx+1,j), dxa(npx-2:npx+1,j))
890  if (ut(npx,j) > 0.) then
891  uc(npx,j) = ut(npx,j)*sin_sg(npx-1,j,3)
892  else
893  uc(npx,j) = ut(npx,j)*sin_sg(npx,j,1)
894  end if
895  uc(npx+1,j) = c3*utmp(npx,j) + c2*utmp(npx+1,j) + c1*utmp(npx+2,j)
896  ut(npx-1,j) = (uc(npx-1,j)-v(npx-1,j)*cosa_u(npx-1,j))*rsin_u(npx-1,j)
897  ut(npx+1,j) = (uc(npx+1,j)-v(npx+1,j)*cosa_u(npx+1,j))*rsin_u(npx+1,j)
898  enddo
899  endif
900 
901 endif
902 
903 !------
904 ! Ydir:
905 !------
906 if( geom%sw_corner ) then
907  do j=-2,0
908  vtmp(0,j) = -utmp(1-j,0)
909  enddo
910  endif
911  if( geom%nw_corner ) then
912  do j=0,2
913  vtmp(0,npy+j) = utmp(j+1,npy)
914  enddo
915  endif
916  if( geom%se_corner ) then
917  do j=-2,0
918  vtmp(npx,j) = utmp(ie+j,0)
919  enddo
920  endif
921  if( geom%ne_corner ) then
922  do j=0,2
923  vtmp(npx,npy+j) = -utmp(ie-j,npy)
924  enddo
925  endif
926  if( geom%sw_corner ) then
927  va(0,-1) = -ua(2,0)
928  va(0, 0) = -ua(1,0)
929  endif
930  if( geom%se_corner ) then
931  va(npx, 0) = ua(npx-1,0)
932  va(npx,-1) = ua(npx-2,0)
933  endif
934  if( geom%ne_corner ) then
935  va(npx,npy ) = -ua(npx-1,npy)
936  va(npx,npy+1) = -ua(npx-2,npy)
937  endif
938  if( geom%nw_corner ) then
939  va(0,npy) = ua(1,npy)
940  va(0,npy+1) = ua(2,npy)
941  endif
942 
943 if (grid_type < 3) then
944 
945  do j=js-1,je+2
946  if ( j==1 .and. .not. nested ) then
947  do i=is-1,ie+1
948  vt(i,j) = edge_interpolate4(va(i,-1:2), dya(i,-1:2))
949  if (vt(i,j) > 0.) then
950  vc(i,j) = vt(i,j)*sin_sg(i,j-1,4)
951  else
952  vc(i,j) = vt(i,j)*sin_sg(i,j,2)
953  end if
954  enddo
955  elseif ( j==0 .or. j==(npy-1) .and. .not. nested ) then
956  do i=is-1,ie+1
957  vc(i,j) = c1*vtmp(i,j-2) + c2*vtmp(i,j-1) + c3*vtmp(i,j)
958  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
959  enddo
960  elseif ( j==2 .or. j==(npy+1) .and. .not. nested ) then
961  do i=is-1,ie+1
962  vc(i,j) = c1*vtmp(i,j+1) + c2*vtmp(i,j) + c3*vtmp(i,j-1)
963  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
964  enddo
965  elseif ( j==npy .and. .not. nested ) then
966  do i=is-1,ie+1
967  vt(i,j) = edge_interpolate4(va(i,j-2:j+1), dya(i,j-2:j+1))
968  if (vt(i,j) > 0.) then
969  vc(i,j) = vt(i,j)*sin_sg(i,j-1,4)
970  else
971  vc(i,j) = vt(i,j)*sin_sg(i,j,2)
972  end if
973  enddo
974  else
975 
976  ! 4th order interpolation for interior points:
977  do i=is-1,ie+1
978  vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1)) + a1*(vtmp(i,j-1)+vtmp(i,j))
979  vt(i,j) = (vc(i,j) - u(i,j)*cosa_v(i,j))*rsin_v(i,j)
980  enddo
981  endif
982  enddo
983 
984 else
985 
986  ! 4th order interpolation:
987  do j=js-1,je+2
988  do i=is-1,ie+1
989  vc(i,j) = a2*(vtmp(i,j-2)+vtmp(i,j+1)) + a1*(vtmp(i,j-1)+vtmp(i,j))
990  vt(i,j) = vc(i,j)
991  enddo
992  enddo
993 
994 endif
995 
996 if (present(ua_out) .and. present(va_out)) then
997  ua_out(is:ie,js:je) = ua(is:ie,js:je)
998  va_out(is:ie,js:je) = va(is:ie,js:je)
999 endif
1000 
1001 end subroutine udvd_to_ucvc
1002 
1003 ! ------------------------------------------------------------------------------
1004 
1005 real(kind=kind_real) function edge_interpolate4(ua, dxa)
1006 
1007 implicit none
1008 real(kind=kind_real), intent(in) :: ua(4)
1009 real(kind=kind_real), intent(in) :: dxa(4)
1010 real(kind=kind_real):: t1, t2
1011 
1012 t1 = dxa(1) + dxa(2)
1013 t2 = dxa(3) + dxa(4)
1014 
1015 edge_interpolate4 = 0.5*( ((t1+dxa(2))*ua(2)-dxa(2)*ua(1)) / t1 + &
1016  ((t2+dxa(3))*ua(3)-dxa(3)*ua(4)) / t2 )
1017 
1018 end function edge_interpolate4
1019 
1020 ! ------------------------------------------------------------------------------
1021 
1022 subroutine a2d(geom, ua, va, ud, vd)
1023 
1024 implicit none
1025 
1026 type(fv3jedi_geom), intent(inout) :: geom
1027 real(kind=kind_real), intent(in) :: ua(geom%isc:geom%iec, geom%jsc:geom%jec, geom%npz)
1028 real(kind=kind_real), intent(in) :: va(geom%isc:geom%iec, geom%jsc:geom%jec, geom%npz)
1029 real(kind=kind_real), intent(inout) :: ud(geom%isc:geom%iec, geom%jsc:geom%jec+1,geom%npz)
1030 real(kind=kind_real), intent(inout) :: vd(geom%isc:geom%iec+1,geom%jsc:geom%jec, geom%npz)
1031 
1032 integer :: is ,ie , js ,je
1033 integer :: npx, npy, npz
1034 integer :: i,j,k, im2,jm2
1035 
1036 real(kind=kind_real) :: uatemp(geom%isd:geom%ied,geom%jsd:geom%jed,geom%npz)
1037 real(kind=kind_real) :: vatemp(geom%isd:geom%ied,geom%jsd:geom%jed,geom%npz)
1038 
1039 real(kind=kind_real) :: v3(geom%isc-1:geom%iec+1,geom%jsc-1:geom%jec+1,3)
1040 real(kind=kind_real) :: ue(geom%isc-1:geom%iec+1,geom%jsc :geom%jec+1,3) ! 3D winds at edges
1041 real(kind=kind_real) :: ve(geom%isc :geom%iec+1,geom%jsc-1:geom%jec+1,3) ! 3D winds at edges
1042 real(kind=kind_real), dimension(geom%isc:geom%iec):: ut1, ut2, ut3
1043 real(kind=kind_real), dimension(geom%jsc:geom%jec):: vt1, vt2, vt3
1044 
1045 npx = geom%npx
1046 npy = geom%npy
1047 npz = geom%npz
1048 is = geom%isc
1049 ie = geom%iec
1050 js = geom%jsc
1051 je = geom%jec
1052 
1053 im2 = (npx-1)/2
1054 jm2 = (npy-1)/2
1055 
1056 uatemp(:,:,:) = 0.0
1057 vatemp(:,:,:) = 0.0
1058 
1059 uatemp(is:ie,js:je,:) = ua
1060 vatemp(is:ie,js:je,:) = va
1061 
1062 call mpp_update_domains(uatemp, geom%domain, complete=.true.)
1063 call mpp_update_domains(vatemp, geom%domain, complete=.true.)
1064 
1065 do k=1, npz
1066 
1067  do j=js-1,je+1
1068  do i=is-1,ie+1
1069  v3(i,j,1) = uatemp(i,j,k)*geom%vlon(i,j,1) + vatemp(i,j,k)*geom%vlat(i,j,1)
1070  v3(i,j,2) = uatemp(i,j,k)*geom%vlon(i,j,2) + vatemp(i,j,k)*geom%vlat(i,j,2)
1071  v3(i,j,3) = uatemp(i,j,k)*geom%vlon(i,j,3) + vatemp(i,j,k)*geom%vlat(i,j,3)
1072  enddo
1073  enddo
1074 
1075  do j=js,je+1
1076  do i=is-1,ie+1
1077  ue(i,j,1) = v3(i,j-1,1) + v3(i,j,1)
1078  ue(i,j,2) = v3(i,j-1,2) + v3(i,j,2)
1079  ue(i,j,3) = v3(i,j-1,3) + v3(i,j,3)
1080  enddo
1081  enddo
1082 
1083  do j=js-1,je+1
1084  do i=is,ie+1
1085  ve(i,j,1) = v3(i-1,j,1) + v3(i,j,1)
1086  ve(i,j,2) = v3(i-1,j,2) + v3(i,j,2)
1087  ve(i,j,3) = v3(i-1,j,3) + v3(i,j,3)
1088  enddo
1089  enddo
1090 
1091  if ( is==1 ) then
1092  i = 1
1093  do j=js,je
1094  if ( j>jm2 ) then
1095  vt1(j) = geom%edge_vect_w(j)*ve(i,j-1,1)+(1.-geom%edge_vect_w(j))*ve(i,j,1)
1096  vt2(j) = geom%edge_vect_w(j)*ve(i,j-1,2)+(1.-geom%edge_vect_w(j))*ve(i,j,2)
1097  vt3(j) = geom%edge_vect_w(j)*ve(i,j-1,3)+(1.-geom%edge_vect_w(j))*ve(i,j,3)
1098  else
1099  vt1(j) = geom%edge_vect_w(j)*ve(i,j+1,1)+(1.-geom%edge_vect_w(j))*ve(i,j,1)
1100  vt2(j) = geom%edge_vect_w(j)*ve(i,j+1,2)+(1.-geom%edge_vect_w(j))*ve(i,j,2)
1101  vt3(j) = geom%edge_vect_w(j)*ve(i,j+1,3)+(1.-geom%edge_vect_w(j))*ve(i,j,3)
1102  endif
1103  enddo
1104  do j=js,je
1105  ve(i,j,1) = vt1(j)
1106  ve(i,j,2) = vt2(j)
1107  ve(i,j,3) = vt3(j)
1108  enddo
1109  endif
1110 
1111  if ( (ie+1)==npx ) then
1112  i = npx
1113  do j=js,je
1114  if ( j>jm2 ) then
1115  vt1(j) = geom%edge_vect_e(j)*ve(i,j-1,1)+(1.-geom%edge_vect_e(j))*ve(i,j,1)
1116  vt2(j) = geom%edge_vect_e(j)*ve(i,j-1,2)+(1.-geom%edge_vect_e(j))*ve(i,j,2)
1117  vt3(j) = geom%edge_vect_e(j)*ve(i,j-1,3)+(1.-geom%edge_vect_e(j))*ve(i,j,3)
1118  else
1119  vt1(j) = geom%edge_vect_e(j)*ve(i,j+1,1)+(1.-geom%edge_vect_e(j))*ve(i,j,1)
1120  vt2(j) = geom%edge_vect_e(j)*ve(i,j+1,2)+(1.-geom%edge_vect_e(j))*ve(i,j,2)
1121  vt3(j) = geom%edge_vect_e(j)*ve(i,j+1,3)+(1.-geom%edge_vect_e(j))*ve(i,j,3)
1122  endif
1123  enddo
1124  do j=js,je
1125  ve(i,j,1) = vt1(j)
1126  ve(i,j,2) = vt2(j)
1127  ve(i,j,3) = vt3(j)
1128  enddo
1129  endif
1130 
1131  if ( js==1 ) then
1132  j = 1
1133  do i=is,ie
1134  if ( i>im2 ) then
1135  ut1(i) = geom%edge_vect_s(i)*ue(i-1,j,1)+(1.-geom%edge_vect_s(i))*ue(i,j,1)
1136  ut2(i) = geom%edge_vect_s(i)*ue(i-1,j,2)+(1.-geom%edge_vect_s(i))*ue(i,j,2)
1137  ut3(i) = geom%edge_vect_s(i)*ue(i-1,j,3)+(1.-geom%edge_vect_s(i))*ue(i,j,3)
1138  else
1139  ut1(i) = geom%edge_vect_s(i)*ue(i+1,j,1)+(1.-geom%edge_vect_s(i))*ue(i,j,1)
1140  ut2(i) = geom%edge_vect_s(i)*ue(i+1,j,2)+(1.-geom%edge_vect_s(i))*ue(i,j,2)
1141  ut3(i) = geom%edge_vect_s(i)*ue(i+1,j,3)+(1.-geom%edge_vect_s(i))*ue(i,j,3)
1142  endif
1143  enddo
1144  do i=is,ie
1145  ue(i,j,1) = ut1(i)
1146  ue(i,j,2) = ut2(i)
1147  ue(i,j,3) = ut3(i)
1148  enddo
1149  endif
1150 
1151  if ( (je+1)==npy ) then
1152  j = npy
1153  do i=is,ie
1154  if ( i>im2 ) then
1155  ut1(i) = geom%edge_vect_n(i)*ue(i-1,j,1)+(1.-geom%edge_vect_n(i))*ue(i,j,1)
1156  ut2(i) = geom%edge_vect_n(i)*ue(i-1,j,2)+(1.-geom%edge_vect_n(i))*ue(i,j,2)
1157  ut3(i) = geom%edge_vect_n(i)*ue(i-1,j,3)+(1.-geom%edge_vect_n(i))*ue(i,j,3)
1158  else
1159  ut1(i) = geom%edge_vect_n(i)*ue(i+1,j,1)+(1.-geom%edge_vect_n(i))*ue(i,j,1)
1160  ut2(i) = geom%edge_vect_n(i)*ue(i+1,j,2)+(1.-geom%edge_vect_n(i))*ue(i,j,2)
1161  ut3(i) = geom%edge_vect_n(i)*ue(i+1,j,3)+(1.-geom%edge_vect_n(i))*ue(i,j,3)
1162  endif
1163  enddo
1164  do i=is,ie
1165  ue(i,j,1) = ut1(i)
1166  ue(i,j,2) = ut2(i)
1167  ue(i,j,3) = ut3(i)
1168  enddo
1169  endif
1170 
1171  do j=js,je+1
1172  do i=is,ie
1173  ud(i,j,k) = 0.5*( ue(i,j,1)*geom%es(1,i,j,1) + &
1174  ue(i,j,2)*geom%es(2,i,j,1) + &
1175  ue(i,j,3)*geom%es(3,i,j,1) )
1176  enddo
1177  enddo
1178 
1179  do j=js,je
1180  do i=is,ie+1
1181  vd(i,j,k) = 0.5*( ve(i,j,1)*geom%ew(1,i,j,2) + &
1182  ve(i,j,2)*geom%ew(2,i,j,2) + &
1183  ve(i,j,3)*geom%ew(3,i,j,2) )
1184  enddo
1185  enddo
1186 
1187 enddo
1188 
1189 end subroutine a2d
1190 
1191 ! ------------------------------------------------------------------------------
1192 
1193 subroutine a2d_ad(geom, ua_ad, va_ad, ud_ad, vd_ad)
1194 
1195 implicit none
1196 type(fv3jedi_geom), intent(inout) :: geom
1197 real(kind=kind_real), intent(inout) :: ua_ad(geom%isc:geom%iec, geom%jsc:geom%jec, geom%npz)
1198 real(kind=kind_real), intent(inout) :: va_ad(geom%isc:geom%iec, geom%jsc:geom%jec, geom%npz)
1199 real(kind=kind_real), intent(inout) :: ud_ad(geom%isc:geom%iec, geom%jsc:geom%jec+1,geom%npz)
1200 real(kind=kind_real), intent(inout) :: vd_ad(geom%isc:geom%iec+1,geom%jsc:geom%jec ,geom%npz)
1201 
1202 integer :: is, ie, js, je
1203 integer :: npx, npy, npz
1204 integer :: i, j, k, im2, jm2
1205 
1206 real(kind=kind_real) :: uatemp(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1207 real(kind=kind_real) :: vatemp(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1208 
1209 real(kind=kind_real) :: uatemp_ad(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1210 real(kind=kind_real) :: vatemp_ad(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1211 real(kind=kind_real) :: v3_ad(geom%isc-1:geom%iec+1, geom%jsc-1:geom%jec+1, 3)
1212 real(kind=kind_real) :: ue_ad(geom%isc-1:geom%iec+1, geom%jsc:geom%jec+1, 3)
1213 real(kind=kind_real) :: ve_ad(geom%isc:geom%iec+1, geom%jsc-1:geom%jec+1, 3)
1214 real(kind=kind_real), dimension(geom%isc:geom%iec) :: ut1_ad, ut2_ad, ut3_ad
1215 real(kind=kind_real), dimension(geom%jsc:geom%jec) :: vt1_ad, vt2_ad, vt3_ad
1216 real(kind=kind_real) :: temp_ad
1217 real(kind=kind_real) :: temp_ad0
1218 
1219 npx = geom%npx
1220 npy = geom%npy
1221 npz = geom%npz
1222 is = geom%isc
1223 ie = geom%iec
1224 js = geom%jsc
1225 je = geom%jec
1226 im2 = (npx-1)/2
1227 jm2 = (npy-1)/2
1228 
1229 uatemp(:, :, :) = 0.0
1230 vatemp(:, :, :) = 0.0
1231 
1232 v3_ad = 0.0_8
1233 uatemp_ad = 0.0_8
1234 ue_ad = 0.0_8
1235 vt1_ad = 0.0_8
1236 vt2_ad = 0.0_8
1237 vt3_ad = 0.0_8
1238 ve_ad = 0.0_8
1239 vatemp_ad = 0.0_8
1240 ut1_ad = 0.0_8
1241 ut2_ad = 0.0_8
1242 ut3_ad = 0.0_8
1243 
1244 do k=npz,1,-1
1245 
1246  do j=je,js,-1
1247  do i=ie+1,is,-1
1248  temp_ad0 = 0.5*vd_ad(i, j, k)
1249  ve_ad(i, j, 1) = ve_ad(i, j, 1) + geom%ew(1, i, j, 2)*temp_ad0
1250  ve_ad(i, j, 2) = ve_ad(i, j, 2) + geom%ew(2, i, j, 2)*temp_ad0
1251  ve_ad(i, j, 3) = ve_ad(i, j, 3) + geom%ew(3, i, j, 2)*temp_ad0
1252  vd_ad(i, j, k) = 0.0_8
1253  end do
1254  end do
1255 
1256  do j=je+1,js,-1
1257  do i=ie,is,-1
1258  temp_ad = 0.5*ud_ad(i, j, k)
1259  ue_ad(i, j, 1) = ue_ad(i, j, 1) + geom%es(1, i, j, 1)*temp_ad
1260  ue_ad(i, j, 2) = ue_ad(i, j, 2) + geom%es(2, i, j, 1)*temp_ad
1261  ue_ad(i, j, 3) = ue_ad(i, j, 3) + geom%es(3, i, j, 1)*temp_ad
1262  ud_ad(i, j, k) = 0.0_8
1263  end do
1264  end do
1265 
1266  if (je + 1 .eq. npy) then
1267  do i=ie,is,-1
1268  ut3_ad(i) = ut3_ad(i) + ue_ad(i, npy, 3)
1269  ue_ad(i, npy, 3) = 0.0_8
1270  ut2_ad(i) = ut2_ad(i) + ue_ad(i, npy, 2)
1271  ue_ad(i, npy, 2) = 0.0_8
1272  ut1_ad(i) = ut1_ad(i) + ue_ad(i, npy, 1)
1273  ue_ad(i, npy, 1) = 0.0_8
1274  end do
1275  do i=ie,is,-1
1276  if (i .le. im2) then
1277  ue_ad(i+1, npy, 3) = ue_ad(i+1, npy, 3) + geom%edge_vect_n(i)*ut3_ad(i)
1278  ue_ad(i, npy, 3) = ue_ad(i, npy, 3) + (1.-geom%edge_vect_n(i))*ut3_ad(i)
1279  ut3_ad(i) = 0.0_8
1280  ue_ad(i+1, npy, 2) = ue_ad(i+1, npy, 2) + geom%edge_vect_n(i)*ut2_ad(i)
1281  ue_ad(i, npy, 2) = ue_ad(i, npy, 2) + (1.-geom%edge_vect_n(i))*ut2_ad(i)
1282  ut2_ad(i) = 0.0_8
1283  ue_ad(i+1, npy, 1) = ue_ad(i+1, npy, 1) + geom%edge_vect_n(i)*ut1_ad(i)
1284  ue_ad(i, npy, 1) = ue_ad(i, npy, 1) + (1.-geom%edge_vect_n(i))*ut1_ad(i)
1285  ut1_ad(i) = 0.0_8
1286  else
1287  ue_ad(i-1, npy, 3) = ue_ad(i-1, npy, 3) + geom%edge_vect_n(i)*ut3_ad(i)
1288  ue_ad(i, npy, 3) = ue_ad(i, npy, 3) + (1.-geom%edge_vect_n(i))*ut3_ad(i)
1289  ut3_ad(i) = 0.0_8
1290  ue_ad(i-1, npy, 2) = ue_ad(i-1, npy, 2) + geom%edge_vect_n(i)*ut2_ad(i)
1291  ue_ad(i, npy, 2) = ue_ad(i, npy, 2) + (1.-geom%edge_vect_n(i))*ut2_ad(i)
1292  ut2_ad(i) = 0.0_8
1293  ue_ad(i-1, npy, 1) = ue_ad(i-1, npy, 1) + geom%edge_vect_n(i)*ut1_ad(i)
1294  ue_ad(i, npy, 1) = ue_ad(i, npy, 1) + (1.-geom%edge_vect_n(i))*ut1_ad(i)
1295  ut1_ad(i) = 0.0_8
1296  end if
1297  end do
1298  end if
1299 
1300  if (js .eq. 1) then
1301  do i=ie,is,-1
1302  ut3_ad(i) = ut3_ad(i) + ue_ad(i, 1, 3)
1303  ue_ad(i, 1, 3) = 0.0_8
1304  ut2_ad(i) = ut2_ad(i) + ue_ad(i, 1, 2)
1305  ue_ad(i, 1, 2) = 0.0_8
1306  ut1_ad(i) = ut1_ad(i) + ue_ad(i, 1, 1)
1307  ue_ad(i, 1, 1) = 0.0_8
1308  end do
1309  do i=ie,is,-1
1310  if (i .le. im2) then
1311  ue_ad(i+1, 1, 3) = ue_ad(i+1, 1, 3) + geom%edge_vect_s(i)*ut3_ad(i)
1312  ue_ad(i, 1, 3) = ue_ad(i, 1, 3) + (1.-geom%edge_vect_s(i))*ut3_ad(i)
1313  ut3_ad(i) = 0.0_8
1314  ue_ad(i+1, 1, 2) = ue_ad(i+1, 1, 2) + geom%edge_vect_s(i)*ut2_ad(i)
1315  ue_ad(i, 1, 2) = ue_ad(i, 1, 2) + (1.-geom%edge_vect_s(i))*ut2_ad(i)
1316  ut2_ad(i) = 0.0_8
1317  ue_ad(i+1, 1, 1) = ue_ad(i+1, 1, 1) + geom%edge_vect_s(i)*ut1_ad(i)
1318  ue_ad(i, 1, 1) = ue_ad(i, 1, 1) + (1.-geom%edge_vect_s(i))*ut1_ad(i)
1319  ut1_ad(i) = 0.0_8
1320  else
1321  ue_ad(i-1, 1, 3) = ue_ad(i-1, 1, 3) + geom%edge_vect_s(i)*ut3_ad(i)
1322  ue_ad(i, 1, 3) = ue_ad(i, 1, 3) + (1.-geom%edge_vect_s(i))*ut3_ad(i)
1323  ut3_ad(i) = 0.0_8
1324  ue_ad(i-1, 1, 2) = ue_ad(i-1, 1, 2) + geom%edge_vect_s(i)*ut2_ad(i)
1325  ue_ad(i, 1, 2) = ue_ad(i, 1, 2) + (1.-geom%edge_vect_s(i))*ut2_ad(i)
1326  ut2_ad(i) = 0.0_8
1327  ue_ad(i-1, 1, 1) = ue_ad(i-1, 1, 1) + geom%edge_vect_s(i)*ut1_ad(i)
1328  ue_ad(i, 1, 1) = ue_ad(i, 1, 1) + (1.-geom%edge_vect_s(i))*ut1_ad(i)
1329  ut1_ad(i) = 0.0_8
1330  end if
1331  end do
1332  end if
1333 
1334  if (ie + 1 .eq. npx) then
1335  do j=je,js,-1
1336  vt3_ad(j) = vt3_ad(j) + ve_ad(npx, j, 3)
1337  ve_ad(npx, j, 3) = 0.0_8
1338  vt2_ad(j) = vt2_ad(j) + ve_ad(npx, j, 2)
1339  ve_ad(npx, j, 2) = 0.0_8
1340  vt1_ad(j) = vt1_ad(j) + ve_ad(npx, j, 1)
1341  ve_ad(npx, j, 1) = 0.0_8
1342  end do
1343  do j=je,js,-1
1344  if (j .le. jm2) then
1345  ve_ad(npx, j+1, 3) = ve_ad(npx, j+1, 3) + geom%edge_vect_e(j)*vt3_ad(j)
1346  ve_ad(npx, j, 3) = ve_ad(npx, j, 3) + (1.-geom%edge_vect_e(j))*vt3_ad(j)
1347  vt3_ad(j) = 0.0_8
1348  ve_ad(npx, j+1, 2) = ve_ad(npx, j+1, 2) + geom%edge_vect_e(j)*vt2_ad(j)
1349  ve_ad(npx, j, 2) = ve_ad(npx, j, 2) + (1.-geom%edge_vect_e(j))*vt2_ad(j)
1350  vt2_ad(j) = 0.0_8
1351  ve_ad(npx, j+1, 1) = ve_ad(npx, j+1, 1) + geom%edge_vect_e(j)*vt1_ad(j)
1352  ve_ad(npx, j, 1) = ve_ad(npx, j, 1) + (1.-geom%edge_vect_e(j))*vt1_ad(j)
1353  vt1_ad(j) = 0.0_8
1354  else
1355  ve_ad(npx, j-1, 3) = ve_ad(npx, j-1, 3) + geom%edge_vect_e(j)*vt3_ad(j)
1356  ve_ad(npx, j, 3) = ve_ad(npx, j, 3) + (1.-geom%edge_vect_e(j))*vt3_ad(j)
1357  vt3_ad(j) = 0.0_8
1358  ve_ad(npx, j-1, 2) = ve_ad(npx, j-1, 2) + geom%edge_vect_e(j)*vt2_ad(j)
1359  ve_ad(npx, j, 2) = ve_ad(npx, j, 2) + (1.-geom%edge_vect_e(j))*vt2_ad(j)
1360  vt2_ad(j) = 0.0_8
1361  ve_ad(npx, j-1, 1) = ve_ad(npx, j-1, 1) + geom%edge_vect_e(j)*vt1_ad(j)
1362  ve_ad(npx, j, 1) = ve_ad(npx, j, 1) + (1.-geom%edge_vect_e(j))*vt1_ad(j)
1363  vt1_ad(j) = 0.0_8
1364  end if
1365  end do
1366  end if
1367 
1368  if (is .eq. 1) then
1369  do j=je,js,-1
1370  vt3_ad(j) = vt3_ad(j) + ve_ad(1, j, 3)
1371  ve_ad(1, j, 3) = 0.0_8
1372  vt2_ad(j) = vt2_ad(j) + ve_ad(1, j, 2)
1373  ve_ad(1, j, 2) = 0.0_8
1374  vt1_ad(j) = vt1_ad(j) + ve_ad(1, j, 1)
1375  ve_ad(1, j, 1) = 0.0_8
1376  end do
1377  do j=je,js,-1
1378  if (j .le. jm2) then
1379  ve_ad(1, j+1, 3) = ve_ad(1, j+1, 3) + geom%edge_vect_w(j)*vt3_ad(j)
1380  ve_ad(1, j, 3) = ve_ad(1, j, 3) + (1.-geom%edge_vect_w(j))*vt3_ad(j)
1381  vt3_ad(j) = 0.0_8
1382  ve_ad(1, j+1, 2) = ve_ad(1, j+1, 2) + geom%edge_vect_w(j)*vt2_ad(j)
1383  ve_ad(1, j, 2) = ve_ad(1, j, 2) + (1.-geom%edge_vect_w(j))*vt2_ad(j)
1384  vt2_ad(j) = 0.0_8
1385  ve_ad(1, j+1, 1) = ve_ad(1, j+1, 1) + geom%edge_vect_w(j)*vt1_ad(j)
1386  ve_ad(1, j, 1) = ve_ad(1, j, 1) + (1.-geom%edge_vect_w(j))*vt1_ad(j)
1387  vt1_ad(j) = 0.0_8
1388  else
1389  ve_ad(1, j-1, 3) = ve_ad(1, j-1, 3) + geom%edge_vect_w(j)*vt3_ad(j)
1390  ve_ad(1, j, 3) = ve_ad(1, j, 3) + (1.-geom%edge_vect_w(j))*vt3_ad(j)
1391  vt3_ad(j) = 0.0_8
1392  ve_ad(1, j-1, 2) = ve_ad(1, j-1, 2) + geom%edge_vect_w(j)*vt2_ad(j)
1393  ve_ad(1, j, 2) = ve_ad(1, j, 2) + (1.-geom%edge_vect_w(j))*vt2_ad(j)
1394  vt2_ad(j) = 0.0_8
1395  ve_ad(1, j-1, 1) = ve_ad(1, j-1, 1) + geom%edge_vect_w(j)*vt1_ad(j)
1396  ve_ad(1, j, 1) = ve_ad(1, j, 1) + (1.-geom%edge_vect_w(j))*vt1_ad(j)
1397  vt1_ad(j) = 0.0_8
1398  end if
1399  end do
1400  end if
1401 
1402  do j=je+1,js-1,-1
1403  do i=ie+1,is,-1
1404  v3_ad(i-1, j, 3) = v3_ad(i-1, j, 3) + ve_ad(i, j, 3)
1405  v3_ad(i, j, 3) = v3_ad(i, j, 3) + ve_ad(i, j, 3)
1406  ve_ad(i, j, 3) = 0.0_8
1407  v3_ad(i-1, j, 2) = v3_ad(i-1, j, 2) + ve_ad(i, j, 2)
1408  v3_ad(i, j, 2) = v3_ad(i, j, 2) + ve_ad(i, j, 2)
1409  ve_ad(i, j, 2) = 0.0_8
1410  v3_ad(i-1, j, 1) = v3_ad(i-1, j, 1) + ve_ad(i, j, 1)
1411  v3_ad(i, j, 1) = v3_ad(i, j, 1) + ve_ad(i, j, 1)
1412  ve_ad(i, j, 1) = 0.0_8
1413  end do
1414  end do
1415 
1416  do j=je+1,js,-1
1417  do i=ie+1,is-1,-1
1418  v3_ad(i, j-1, 3) = v3_ad(i, j-1, 3) + ue_ad(i, j, 3)
1419  v3_ad(i, j, 3) = v3_ad(i, j, 3) + ue_ad(i, j, 3)
1420  ue_ad(i, j, 3) = 0.0_8
1421  v3_ad(i, j-1, 2) = v3_ad(i, j-1, 2) + ue_ad(i, j, 2)
1422  v3_ad(i, j, 2) = v3_ad(i, j, 2) + ue_ad(i, j, 2)
1423  ue_ad(i, j, 2) = 0.0_8
1424  v3_ad(i, j-1, 1) = v3_ad(i, j-1, 1) + ue_ad(i, j, 1)
1425  v3_ad(i, j, 1) = v3_ad(i, j, 1) + ue_ad(i, j, 1)
1426  ue_ad(i, j, 1) = 0.0_8
1427  end do
1428  end do
1429 
1430  do j=je+1,js-1,-1
1431  do i=ie+1,is-1,-1
1432  uatemp_ad(i, j, k) = uatemp_ad(i, j, k) + geom%vlon(i, j, 3)*v3_ad(i, j, 3)
1433  vatemp_ad(i, j, k) = vatemp_ad(i, j, k) + geom%vlat(i, j, 3)*v3_ad(i, j, 3)
1434  v3_ad(i, j, 3) = 0.0_8
1435  uatemp_ad(i, j, k) = uatemp_ad(i, j, k) + geom%vlon(i, j, 2)*v3_ad(i, j, 2)
1436  vatemp_ad(i, j, k) = vatemp_ad(i, j, k) + geom%vlat(i, j, 2)*v3_ad(i, j, 2)
1437  v3_ad(i, j, 2) = 0.0_8
1438  uatemp_ad(i, j, k) = uatemp_ad(i, j, k) + geom%vlon(i, j, 1)*v3_ad(i, j, 1)
1439  vatemp_ad(i, j, k) = vatemp_ad(i, j, k) + geom%vlat(i, j, 1)*v3_ad(i, j, 1)
1440  v3_ad(i, j, 1) = 0.0_8
1441  end do
1442  end do
1443 end do
1444 
1445 call mpp_update_domains_adm(vatemp, vatemp_ad, geom%domain, complete=.true.)
1446 call mpp_update_domains_adm(uatemp, uatemp_ad, geom%domain, complete=.true.)
1447 
1448 va_ad = va_ad + vatemp_ad(is:ie, js:je, :)
1449 ua_ad = ua_ad + uatemp_ad(is:ie, js:je, :)
1450 
1451 end subroutine a2d_ad
1452 
1453 ! ------------------------------------------------------------------------------
1454 
1455 subroutine d2a(geom, u_comp, v_comp, ua_comp, va_comp)
1456 
1457 implicit none
1458 
1459 type(fv3jedi_geom), intent(inout) :: geom
1460 real(kind=kind_real), intent(inout) :: u_comp(geom%isc:geom%iec ,geom%jsc:geom%jec+1,geom%npz)
1461 real(kind=kind_real), intent(inout) :: v_comp(geom%isc:geom%iec+1,geom%jsc:geom%jec ,geom%npz)
1462 real(kind=kind_real), intent(inout) :: ua_comp(geom%isc:geom%iec ,geom%jsc:geom%jec ,geom%npz)
1463 real(kind=kind_real), intent(inout) :: va_comp(geom%isc:geom%iec ,geom%jsc:geom%jec ,geom%npz)
1464 
1465 real(kind=kind_real) :: c1 = 1.125
1466 real(kind=kind_real) :: c2 = -0.125
1467 real(kind=kind_real) :: utmp(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1468 real(kind=kind_real) :: vtmp(geom%isc:geom%iec+1,geom%jsc:geom%jec)
1469 real(kind=kind_real) :: wu(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1470 real(kind=kind_real) :: wv(geom%isc:geom%iec+1,geom%jsc:geom%jec)
1471 
1472 real(kind=kind_real) :: u(geom%isd:geom%ied ,geom%jsd:geom%jed+1,geom%npz)
1473 real(kind=kind_real) :: v(geom%isd:geom%ied+1,geom%jsd:geom%jed ,geom%npz)
1474 real(kind=kind_real) :: ua(geom%isd:geom%ied ,geom%jsd:geom%jed ,geom%npz)
1475 real(kind=kind_real) :: va(geom%isd:geom%ied ,geom%jsd:geom%jed ,geom%npz)
1476 
1477 real(kind=kind_real) :: ebuffery(geom%jsc:geom%jec,1:geom%npz)
1478 real(kind=kind_real) :: nbufferx(geom%isc:geom%iec,1:geom%npz)
1479 real(kind=kind_real) :: wbuffery(geom%jsc:geom%jec,1:geom%npz)
1480 real(kind=kind_real) :: sbufferx(geom%isc:geom%iec,1:geom%npz)
1481 
1482 integer i, j, k
1483 integer :: is, ie, js, je, npx, npy, npz
1484 
1485 !Shorten notation
1486 is = geom%isc
1487 ie = geom%iec
1488 js = geom%jsc
1489 je = geom%jec
1490 npx = geom%npx
1491 npy = geom%npy
1492 npz = geom%npz
1493 
1494 !Fill compute part from input
1495 u = 0.0_kind_real
1496 v = 0.0_kind_real
1497 u(is:ie ,js:je+1,:) = u_comp
1498 v(is:ie+1,js:je ,:) = v_comp
1499 
1500 !Fill edges
1501 ebuffery = 0.0_kind_real
1502 nbufferx = 0.0_kind_real
1503 wbuffery = 0.0_kind_real
1504 sbufferx = 0.0_kind_real
1505 call mpp_get_boundary( u, v, geom%domain, &
1506  wbuffery=wbuffery, ebuffery=ebuffery, &
1507  sbufferx=sbufferx, nbufferx=nbufferx, &
1508  gridtype=dgrid_ne, complete=.true. )
1509 do k=1,npz
1510  do i=is,ie
1511  u(i,je+1,k) = nbufferx(i,k)
1512  enddo
1513 enddo
1514 do k=1,npz
1515  do j=js,je
1516  v(ie+1,j,k) = ebuffery(j,k)
1517  enddo
1518 enddo
1519 
1520 !Fill halos
1521 call mpp_update_domains(u, v, geom%domain, gridtype=dgrid_ne)
1522 
1523 ! TODO: OpenMP Re-enable this parallel section by fixing variable shared/private assignment
1524 ! !$OMP parallel do default(none) shared(is,ie,js,je,npz,npx,npy,c2,c1, &
1525 ! !$OMP u,v,ua,va) &
1526 ! !$OMP private(utmp, vtmp, wu, wv)
1527 do k=1,npz
1528 
1529  do j=max(2,js),min(npy-2,je)
1530  do i=max(2,is),min(npx-2,ie)
1531  utmp(i,j) = c2*(u(i,j-1,k)+u(i,j+2,k)) + c1*(u(i,j,k)+u(i,j+1,k))
1532  vtmp(i,j) = c2*(v(i-1,j,k)+v(i+2,j,k)) + c1*(v(i,j,k)+v(i+1,j,k))
1533  enddo
1534  enddo
1535 
1536  if ( js==1 ) then
1537  do i=is,ie+1
1538  wv(i,1) = v(i,1,k)*geom%dy(i,1)
1539  enddo
1540  do i=is,ie
1541  vtmp(i,1) = 2.*(wv(i,1) + wv(i+1,1)) / (geom%dy(i,1)+geom%dy(i+1,1))
1542  utmp(i,1) = 2.*(u(i,1,k)*geom%dx(i,1) + u(i,2,k)*geom%dx(i,2)) &
1543  / ( geom%dx(i,1) + geom%dx(i,2))
1544  enddo
1545  endif
1546 
1547  if ( (je+1)==npy ) then
1548  j = npy-1
1549  do i=is,ie+1
1550  wv(i,j) = v(i,j,k)*geom%dy(i,j)
1551  enddo
1552  do i=is,ie
1553  vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j)) / (geom%dy(i,j)+geom%dy(i+1,j))
1554  utmp(i,j) = 2.*(u(i,j,k)*geom%dx(i,j) + u(i,j+1,k)*geom%dx(i,j+1)) &
1555  / ( geom%dx(i,j) + geom%dx(i,j+1))
1556  enddo
1557  endif
1558 
1559  if ( is==1 ) then
1560  i = 1
1561  do j=js,je
1562  wv(1,j) = v(1,j,k)*geom%dy(1,j)
1563  wv(2,j) = v(2,j,k)*geom%dy(2,j)
1564  enddo
1565  do j=js,je+1
1566  wu(i,j) = u(i,j,k)*geom%dx(i,j)
1567  enddo
1568  do j=js,je
1569  utmp(i,j) = 2.*(wu(i,j) + wu(i,j+1))/(geom%dx(i,j)+geom%dx(i,j+1))
1570  vtmp(i,j) = 2.*(wv(1,j) + wv(2,j ))/(geom%dy(1,j)+geom%dy(2,j))
1571  enddo
1572  endif
1573 
1574  if ( (ie+1)==npx) then
1575  i = npx-1
1576  do j=js,je
1577  wv(i, j) = v(i, j,k)*geom%dy(i, j)
1578  wv(i+1,j) = v(i+1,j,k)*geom%dy(i+1,j)
1579  enddo
1580  do j=js,je+1
1581  wu(i,j) = u(i,j,k)*geom%dx(i,j)
1582  enddo
1583  do j=js,je
1584  utmp(i,j) = 2.*(wu(i,j) + wu(i, j+1))/(geom%dx(i,j)+geom%dx(i,j+1))
1585  vtmp(i,j) = 2.*(wv(i,j) + wv(i+1,j ))/(geom%dy(i,j)+geom%dy(i+1,j))
1586  enddo
1587  endif
1588 
1589  !Transform local a-grid winds into latitude-longitude coordinates
1590  do j=js,je
1591  do i=is,ie
1592  ua(i,j,k) = geom%a11(i,j)*utmp(i,j) + geom%a12(i,j)*vtmp(i,j)
1593  va(i,j,k) = geom%a21(i,j)*utmp(i,j) + geom%a22(i,j)*vtmp(i,j)
1594  enddo
1595  enddo
1596 
1597 enddo
1598 
1599 ua_comp = ua(is:ie,js:je,:)
1600 va_comp = va(is:ie,js:je,:)
1601 
1602 end subroutine d2a
1603 
1604 ! ------------------------------------------------------------------------------
1605 
1606 subroutine d2a_ad(geom, u_ad_comp, v_ad_comp, ua_ad_comp, va_ad_comp)
1607 
1608 implicit none
1609 type(fv3jedi_geom), intent(inout) :: geom
1610 real(kind=kind_real), intent(inout) :: u_ad_comp(geom%isc:geom%iec, geom%jsc:geom%jec+1,1:geom%npz)
1611 real(kind=kind_real), intent(inout) :: v_ad_comp(geom%isc:geom%iec+1,geom%jsc:geom%jec, 1:geom%npz)
1612 real(kind=kind_real), intent(inout) :: ua_ad_comp(geom%isc:geom%iec, geom%jsc:geom%jec, 1:geom%npz)
1613 real(kind=kind_real), intent(inout) :: va_ad_comp(geom%isc:geom%iec, geom%jsc:geom%jec, 1:geom%npz)
1614 
1615 real(kind=kind_real) :: u(geom%isd:geom%ied, geom%jsd:geom%jed+1,geom%npz)
1616 real(kind=kind_real) :: v(geom%isd:geom%ied+1,geom%jsd:geom%jed, geom%npz)
1617 real(kind=kind_real) :: ua(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1618 real(kind=kind_real) :: va(geom%isd:geom%ied, geom%jsd:geom%jed, geom%npz)
1619 
1620 real(kind=kind_real), save :: c1=1.125
1621 real(kind=kind_real), save :: c2=-0.125
1622 real(kind=kind_real) :: utmp(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1623 real(kind=kind_real) :: utmp_ad(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1624 real(kind=kind_real) :: vtmp(geom%isc:geom%iec+1, geom%jsc:geom%jec)
1625 real(kind=kind_real) :: vtmp_ad(geom%isc:geom%iec+1, geom%jsc:geom%jec)
1626 real(kind=kind_real) :: wu(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1627 real(kind=kind_real) :: wu_ad(geom%isc:geom%iec, geom%jsc:geom%jec+1)
1628 real(kind=kind_real) :: wv(geom%isc:geom%iec+1, geom%jsc:geom%jec)
1629 real(kind=kind_real) :: wv_ad(geom%isc:geom%iec+1, geom%jsc:geom%jec)
1630 integer :: i, j, k
1631 integer :: is, ie, js, je, npx, npy, npz
1632 intrinsic max
1633 intrinsic min
1634 integer :: max1
1635 integer :: max2
1636 integer :: min1
1637 integer :: min2
1638 real(kind=kind_real) :: temp_ad
1639 real(kind=kind_real) :: temp_ad0
1640 real(kind=kind_real) :: temp_ad1
1641 real(kind=kind_real) :: temp_ad2
1642 real(kind=kind_real) :: temp_ad3
1643 real(kind=kind_real) :: temp_ad4
1644 real(kind=kind_real) :: temp_ad5
1645 real(kind=kind_real) :: temp_ad6
1646 integer :: ad_from
1647 integer :: ad_to
1648 integer :: ad_from0
1649 integer :: ad_to0
1650 integer :: branch
1651 
1652 real(kind=kind_real) :: u_ad(geom%isd:geom%ied, geom%jsd:geom%jed+1,1:geom%npz)
1653 real(kind=kind_real) :: v_ad(geom%isd:geom%ied+1,geom%jsd:geom%jed, 1:geom%npz)
1654 real(kind=kind_real) :: ua_ad(geom%isd:geom%ied, geom%jsd:geom%jed, 1:geom%npz)
1655 real(kind=kind_real) :: va_ad(geom%isd:geom%ied, geom%jsd:geom%jed, 1:geom%npz)
1656 
1657 real(kind=kind_real) :: ebuffery(geom%jsc:geom%jec,1:geom%npz)
1658 real(kind=kind_real) :: nbufferx(geom%isc:geom%iec,1:geom%npz)
1659 real(kind=kind_real) :: wbuffery(geom%jsc:geom%jec,1:geom%npz)
1660 real(kind=kind_real) :: sbufferx(geom%isc:geom%iec,1:geom%npz)
1661 
1662 is = geom%isc
1663 ie = geom%iec
1664 js = geom%jsc
1665 je = geom%jec
1666 npx = geom%npx
1667 npy = geom%npy
1668 npz = geom%npz
1669 
1670 ua_ad = 0.0_kind_real
1671 va_ad = 0.0_kind_real
1672 u_ad = 0.0_kind_real
1673 v_ad = 0.0_kind_real
1674 
1675 ua_ad(is:ie,js:je,:) = ua_ad_comp
1676 va_ad(is:ie,js:je,:) = va_ad_comp
1677 
1678 ! TODO: OpenMP Re-enable this parallel section by fixing variable shared/private assignment
1679 ! !$omp parallel do default(none) shared(is,ie,js,je,npz,npx,npy,c2,c1, &
1680 ! !$omp u,v,ua,va) &
1681 ! !$omp private(utmp, vtmp, wu, wv)
1682 do k=1,npz
1683  if (2 .lt. js) then
1684  max1 = js
1685  else
1686  max1 = 2
1687  end if
1688  if (npy - 2 .gt. je) then
1689  min1 = je
1690  else
1691  min1 = npy - 2
1692  end if
1693  ad_from0 = max1
1694  do j=ad_from0,min1
1695  if (2 .lt. is) then
1696  max2 = is
1697  else
1698  max2 = 2
1699  end if
1700  if (npx - 2 .gt. ie) then
1701  min2 = ie
1702  else
1703  min2 = npx - 2
1704  end if
1705  ad_from = max2
1706  call pushinteger4(i)
1707  i = min2 + 1
1708  call pushinteger4(i - 1)
1709  call pushinteger4(ad_from)
1710  end do
1711  call pushinteger4(j - 1)
1712  call pushinteger4(ad_from0)
1713  if (js .eq. 1) then
1714  call pushinteger4(i)
1715  call pushcontrol1b(0)
1716  else
1717  call pushcontrol1b(1)
1718  end if
1719  if (je + 1 .eq. npy) then
1720  j = npy - 1
1721  call pushinteger4(i)
1722  call pushcontrol1b(0)
1723  else
1724  call pushcontrol1b(1)
1725  end if
1726  if (is .eq. 1) then
1727  call pushinteger4(i)
1728  i = 1
1729  call pushinteger4(j)
1730  call pushcontrol1b(0)
1731  else
1732  call pushcontrol1b(1)
1733  end if
1734  if (ie + 1 .eq. npx) then
1735  call pushinteger4(i)
1736  i = npx - 1
1737  call pushinteger4(j)
1738  call pushcontrol1b(1)
1739  else
1740  call pushcontrol1b(0)
1741  end if
1742  call pushinteger4(j)
1743  do j=js,je
1744  call pushinteger4(i)
1745  end do
1746 end do
1747 vtmp_ad = 0.0_8
1748 wu_ad = 0.0_8
1749 wv_ad = 0.0_8
1750 utmp_ad = 0.0_8
1751 do k=npz,1,-1
1752  do j=je,js,-1
1753  do i=ie,is,-1
1754  utmp_ad(i, j) = utmp_ad(i, j) + geom%a11(i, j)*ua_ad(i, j, k) + geom%a21(i, j)*va_ad(i, j, k)
1755  vtmp_ad(i, j) = vtmp_ad(i, j) + geom%a12(i, j)*ua_ad(i, j, k) + geom%a22(i, j)*va_ad(i, j, k)
1756  va_ad(i, j, k) = 0.0_8
1757  ua_ad(i, j, k) = 0.0_8
1758  end do
1759  call popinteger4(i)
1760  end do
1761  call popinteger4(j)
1762  call popcontrol1b(branch)
1763  if (branch .ne. 0) then
1764  do j=je,js,-1
1765  temp_ad5 = 2.*vtmp_ad(i, j)/(geom%dy(i, j)+geom%dy(i+1, j))
1766  wv_ad(i, j) = wv_ad(i, j) + temp_ad5
1767  wv_ad(i+1, j) = wv_ad(i+1, j) + temp_ad5
1768  vtmp_ad(i, j) = 0.0_8
1769  temp_ad6 = 2.*utmp_ad(i, j)/(geom%dx(i, j)+geom%dx(i, j+1))
1770  wu_ad(i, j) = wu_ad(i, j) + temp_ad6
1771  wu_ad(i, j+1) = wu_ad(i, j+1) + temp_ad6
1772  utmp_ad(i, j) = 0.0_8
1773  end do
1774  do j=je+1,js,-1
1775  u_ad(i, j, k) = u_ad(i, j, k) + geom%dx(i, j)*wu_ad(i, j)
1776  wu_ad(i, j) = 0.0_8
1777  end do
1778  do j=je,js,-1
1779  v_ad(i+1, j, k) = v_ad(i+1, j, k) + geom%dy(i+1, j)*wv_ad(i+1, j)
1780  wv_ad(i+1, j) = 0.0_8
1781  v_ad(i, j, k) = v_ad(i, j, k) + geom%dy(i, j)*wv_ad(i, j)
1782  wv_ad(i, j) = 0.0_8
1783  end do
1784  call popinteger4(j)
1785  call popinteger4(i)
1786  end if
1787  call popcontrol1b(branch)
1788  if (branch .eq. 0) then
1789  do j=je,js,-1
1790  temp_ad3 = 2.*vtmp_ad(i, j)/(geom%dy(1, j)+geom%dy(2, j))
1791  wv_ad(1, j) = wv_ad(1, j) + temp_ad3
1792  wv_ad(2, j) = wv_ad(2, j) + temp_ad3
1793  vtmp_ad(i, j) = 0.0_8
1794  temp_ad4 = 2.*utmp_ad(i, j)/(geom%dx(i, j)+geom%dx(i, j+1))
1795  wu_ad(i, j) = wu_ad(i, j) + temp_ad4
1796  wu_ad(i, j+1) = wu_ad(i, j+1) + temp_ad4
1797  utmp_ad(i, j) = 0.0_8
1798  end do
1799  do j=je+1,js,-1
1800  u_ad(i, j, k) = u_ad(i, j, k) + geom%dx(i, j)*wu_ad(i, j)
1801  wu_ad(i, j) = 0.0_8
1802  end do
1803  do j=je,js,-1
1804  v_ad(2, j, k) = v_ad(2, j, k) + geom%dy(2, j)*wv_ad(2, j)
1805  wv_ad(2, j) = 0.0_8
1806  v_ad(1, j, k) = v_ad(1, j, k) + geom%dy(1, j)*wv_ad(1, j)
1807  wv_ad(1, j) = 0.0_8
1808  end do
1809  call popinteger4(j)
1810  call popinteger4(i)
1811  end if
1812  call popcontrol1b(branch)
1813  if (branch .eq. 0) then
1814  do i=ie,is,-1
1815  temp_ad1 = 2.*utmp_ad(i, j)/(geom%dx(i, j)+geom%dx(i, j+1))
1816  u_ad(i, j, k) = u_ad(i, j, k) + geom%dx(i, j)*temp_ad1
1817  u_ad(i, j+1, k) = u_ad(i, j+1, k) + geom%dx(i, j+1)*temp_ad1
1818  utmp_ad(i, j) = 0.0_8
1819  temp_ad2 = 2.*vtmp_ad(i, j)/(geom%dy(i, j)+geom%dy(i+1, j))
1820  wv_ad(i, j) = wv_ad(i, j) + temp_ad2
1821  wv_ad(i+1, j) = wv_ad(i+1, j) + temp_ad2
1822  vtmp_ad(i, j) = 0.0_8
1823  end do
1824  do i=ie+1,is,-1
1825  v_ad(i, j, k) = v_ad(i, j, k) + geom%dy(i, j)*wv_ad(i, j)
1826  wv_ad(i, j) = 0.0_8
1827  end do
1828  call popinteger4(i)
1829  end if
1830  call popcontrol1b(branch)
1831  if (branch .eq. 0) then
1832  do i=ie,is,-1
1833  temp_ad = 2.*utmp_ad(i, 1)/(geom%dx(i, 1)+geom%dx(i, 2))
1834  u_ad(i, 1, k) = u_ad(i, 1, k) + geom%dx(i, 1)*temp_ad
1835  u_ad(i, 2, k) = u_ad(i, 2, k) + geom%dx(i, 2)*temp_ad
1836  utmp_ad(i, 1) = 0.0_8
1837  temp_ad0 = 2.*vtmp_ad(i, 1)/(geom%dy(i, 1)+geom%dy(i+1, 1))
1838  wv_ad(i, 1) = wv_ad(i, 1) + temp_ad0
1839  wv_ad(i+1, 1) = wv_ad(i+1, 1) + temp_ad0
1840  vtmp_ad(i, 1) = 0.0_8
1841  end do
1842  do i=ie+1,is,-1
1843  v_ad(i, 1, k) = v_ad(i, 1, k) + geom%dy(i, 1)*wv_ad(i, 1)
1844  wv_ad(i, 1) = 0.0_8
1845  end do
1846  call popinteger4(i)
1847  end if
1848  call popinteger4(ad_from0)
1849  call popinteger4(ad_to0)
1850  do j=ad_to0,ad_from0,-1
1851  call popinteger4(ad_from)
1852  call popinteger4(ad_to)
1853  do i=ad_to,ad_from,-1
1854  v_ad(i-1, j, k) = v_ad(i-1, j, k) + c2*vtmp_ad(i, j)
1855  v_ad(i+2, j, k) = v_ad(i+2, j, k) + c2*vtmp_ad(i, j)
1856  v_ad(i, j, k) = v_ad(i, j, k) + c1*vtmp_ad(i, j)
1857  v_ad(i+1, j, k) = v_ad(i+1, j, k) + c1*vtmp_ad(i, j)
1858  vtmp_ad(i, j) = 0.0_8
1859  u_ad(i, j-1, k) = u_ad(i, j-1, k) + c2*utmp_ad(i, j)
1860  u_ad(i, j+2, k) = u_ad(i, j+2, k) + c2*utmp_ad(i, j)
1861  u_ad(i, j, k) = u_ad(i, j, k) + c1*utmp_ad(i, j)
1862  u_ad(i, j+1, k) = u_ad(i, j+1, k) + c1*utmp_ad(i, j)
1863  utmp_ad(i, j) = 0.0_8
1864  end do
1865  call popinteger4(i)
1866  end do
1867 end do
1868 
1869 !Adjoint fill halos
1870 call mpp_update_domains_adm(u, u_ad, v, v_ad, geom%domain, gridtype=dgrid_ne)
1871 
1872 !Adjoint fill edges
1873 ebuffery = 0.0_kind_real
1874 nbufferx = 0.0_kind_real
1875 wbuffery = 0.0_kind_real
1876 sbufferx = 0.0_kind_real
1877 
1878 do k=1,npz
1879  do i=is,ie
1880  nbufferx(i,k) = u_ad(i,je+1,k)
1881  enddo
1882 enddo
1883 do k=1,npz
1884  do j=js,je
1885  ebuffery(j,k) = v_ad(ie+1,j,k)
1886  enddo
1887 enddo
1888 
1889 call mpp_get_boundary_ad( u_ad, v_ad, geom%domain, &
1890  wbuffery=wbuffery, ebuffery=ebuffery, sbufferx=sbufferx, nbufferx=nbufferx, &
1891  gridtype=dgrid_ne, complete=.true. )
1892 
1893 !Output compute part
1894 u_ad_comp = u_ad(is:ie ,js:je+1,:)
1895 v_ad_comp = v_ad(is:ie+1,js:je ,:)
1896 
1897 end subroutine d2a_ad
1898 
1899 ! ------------------------------------------------------------------------------
1900 
1901 subroutine fill_dgrid_winds(geom,u,v,fillhalo)
1902 
1903 implicit none
1904 type(fv3jedi_geom), intent(inout) :: geom
1905 real(kind=kind_real), intent(inout) :: u(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz)
1906 real(kind=kind_real), intent(inout) :: v(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz)
1907 logical, optional, intent(in) :: fillhalo
1908 
1909 integer :: isc, iec, jsc, jec, npz, i, j, k
1910 real(kind=kind_real) :: ebuffery(geom%jsc:geom%jec,1:geom%npz)
1911 real(kind=kind_real) :: nbufferx(geom%isc:geom%iec,1:geom%npz)
1912 real(kind=kind_real) :: wbuffery(geom%jsc:geom%jec,1:geom%npz)
1913 real(kind=kind_real) :: sbufferx(geom%isc:geom%iec,1:geom%npz)
1914 
1915 ! ---------------------------------------- !
1916 ! Fill edge and then halo of D-grid winds !
1917 ! ---------------------------------------- !
1918 
1919 ! Shortcuts
1920 ! ---------
1921 isc = geom%isc
1922 iec = geom%iec
1923 jsc = geom%jsc
1924 jec = geom%jec
1925 npz = geom%npz
1926 
1927 ! Fill north/east edges
1928 ! ---------------------
1929 ebuffery = 0.0_kind_real
1930 nbufferx = 0.0_kind_real
1931 wbuffery = 0.0_kind_real
1932 sbufferx = 0.0_kind_real
1933 call mpp_get_boundary( u, v, geom%domain, &
1934  wbuffery=wbuffery, ebuffery=ebuffery, &
1935  sbufferx=sbufferx, nbufferx=nbufferx, &
1936  gridtype=dgrid_ne, complete=.true. )
1937 do k=1,npz
1938  do i=isc,iec
1939  u(i,jec+1,k) = nbufferx(i,k)
1940  enddo
1941 enddo
1942 do k=1,npz
1943  do j=jsc,jec
1944  v(iec+1,j,k) = ebuffery(j,k)
1945  enddo
1946 enddo
1947 
1948 ! Fill halos
1949 ! ----------
1950 if (present(fillhalo)) then
1951  if (fillhalo) then
1952  call mpp_update_domains(u, v, geom%domain, gridtype=dgrid_ne)
1953  endif
1954 endif
1955 
1956 end subroutine fill_dgrid_winds
1957 
1958 ! ------------------------------------------------------------------------------
1959 
1960 subroutine fill_cgrid_winds(geom,uc,vc,fillhalo)
1961 
1962 implicit none
1963 type(fv3jedi_geom), intent(inout) :: geom
1964 real(kind=kind_real), intent(inout) :: uc(geom%isd:geom%ied+1,geom%jsd:geom%jed ,1:geom%npz)
1965 real(kind=kind_real), intent(inout) :: vc(geom%isd:geom%ied ,geom%jsd:geom%jed+1,1:geom%npz)
1966 logical, optional, intent(in) :: fillhalo
1967 
1968 integer :: isc, iec, jsc, jec, npz, i, j, k
1969 real(kind=kind_real) :: wbufferx(geom%jsc:geom%jec,1:geom%npz)
1970 real(kind=kind_real) :: sbuffery(geom%isc:geom%iec,1:geom%npz)
1971 real(kind=kind_real) :: ebufferx(geom%jsc:geom%jec,1:geom%npz)
1972 real(kind=kind_real) :: nbuffery(geom%isc:geom%iec,1:geom%npz)
1973 
1974 ! ---------------------------------------- !
1975 ! Fill edge and then halo of C-grid winds !
1976 ! ---------------------------------------- !
1977 
1978 ! Shortcuts
1979 ! ---------
1980 isc = geom%isc
1981 iec = geom%iec
1982 jsc = geom%jsc
1983 jec = geom%jec
1984 npz = geom%npz
1985 
1986 ! Fill north/east edges
1987 ! ---------------------
1988 ebufferx = 0.0_kind_real
1989 nbuffery = 0.0_kind_real
1990 wbufferx = 0.0_kind_real
1991 sbuffery = 0.0_kind_real
1992 
1993 call mpp_get_boundary(uc, vc, geom%domain, &
1994  wbufferx=wbufferx, ebufferx=ebufferx, &
1995  sbuffery=sbuffery, nbuffery=nbuffery, &
1996  gridtype=cgrid_ne, complete=.true.)
1997 do k=1,npz
1998  do j=jsc,jec
1999  uc(iec+1,j,k) = ebufferx(j,k)
2000  enddo
2001  do i=isc,iec
2002  vc(i,jec+1,k) = nbuffery(i,k)
2003  enddo
2004 enddo
2005 
2006 ! Fill halos
2007 ! ----------
2008 if (present(fillhalo)) then
2009  if (fillhalo) then
2010  call mpp_update_domains(uc, vc, geom%domain, gridtype=cgrid_ne)
2011  endif
2012 endif
2013 
2014 end subroutine fill_cgrid_winds
2015 
2016 ! ------------------------------------------------------------------------------
2017 
2018 subroutine compute_utvt(geom, uc, vc, ut, vt, dt)
2019 
2020  implicit none
2021  type(fv3jedi_geom), intent(inout) :: geom
2022  real(kind=kind_real), intent(in ) :: uc(geom%isd:geom%ied+1,geom%jsd:geom%jed )
2023  real(kind=kind_real), intent(in ) :: vc(geom%isd:geom%ied ,geom%jsd:geom%jed+1)
2024  real(kind=kind_real), intent(inout) :: ut(geom%isd:geom%ied+1,geom%jsd:geom%jed )
2025  real(kind=kind_real), intent(inout) :: vt(geom%isd:geom%ied ,geom%jsd:geom%jed+1)
2026  real(kind=kind_real), intent(in ) :: dt
2027 
2028 ! Local vars
2029  integer :: is,js,ie,je,isd,jsd,ied,jed,npx,npy
2030  real(kind=kind_real) :: damp
2031  integer i,j
2032 
2033  is = geom%isc
2034  js = geom%jsc
2035  ie = geom%iec
2036  je = geom%jec
2037  isd = geom%isd
2038  jsd = geom%jsd
2039  ied = geom%ied
2040  jed = geom%jed
2041  npx = geom%npx
2042  npy = geom%npy
2043 
2044  if ( geom%grid_type < 3 ) then
2045 
2046 ! Center of Cube Faces
2047  do j=jsd,jed
2048  if(j/=0 .and. j/=1 .and. j/=(npy-1) .and. j/=npy) then
2049  do i=is-1,ie+2
2050  ut(i,j) = ( uc(i,j) - 0.25 * geom%cosa_u(i,j) * &
2051  (vc(i-1,j)+vc(i,j)+vc(i-1,j+1)+vc(i,j+1)))*geom%rsin_u(i,j)
2052  enddo
2053  endif
2054  enddo
2055  do j=js-1,je+2
2056  if( j/=1 .and. j/=npy ) then
2057  do i=isd,ied
2058  vt(i,j) = ( vc(i,j) - 0.25 * geom%cosa_v(i,j) * &
2059  (uc(i,j-1)+uc(i+1,j-1)+uc(i,j)+uc(i+1,j)))*geom%rsin_v(i,j)
2060  enddo
2061  endif
2062  enddo
2063 
2064 ! West edge:
2065  if ( is==1 ) then
2066  do j=jsd,jed
2067  if ( uc(1,j)*dt > 0. ) then
2068  ut(1,j) = uc(1,j) / geom%sin_sg(0,j,3)
2069  else
2070  ut(1,j) = uc(1,j) / geom%sin_sg(1,j,1)
2071  endif
2072  enddo
2073  do j=max(3,js), min(npy-2,je+1)
2074  vt(0,j) = vc(0,j) - 0.25*geom%cosa_v(0,j)* &
2075  (ut(0,j-1)+ut(1,j-1)+ut(0,j)+ut(1,j))
2076  vt(1,j) = vc(1,j) - 0.25*geom%cosa_v(1,j)* &
2077  (ut(1,j-1)+ut(2,j-1)+ut(1,j)+ut(2,j))
2078  enddo
2079  endif ! West face
2080 
2081 ! East edge:
2082  if ( (ie+1)==npx ) then
2083  do j=jsd,jed
2084  if ( uc(npx,j)*dt > 0. ) then
2085  ut(npx,j) = uc(npx,j) / geom%sin_sg(npx-1,j,3)
2086  else
2087  ut(npx,j) = uc(npx,j) / geom%sin_sg(npx,j,1)
2088  endif
2089  enddo
2090 
2091  do j=max(3,js), min(npy-2,je+1)
2092  vt(npx-1,j) = vc(npx-1,j) - 0.25*geom%cosa_v(npx-1,j)* &
2093  (ut(npx-1,j-1)+ut(npx,j-1)+ut(npx-1,j)+ut(npx,j))
2094  vt(npx,j) = vc(npx,j) - 0.25*geom%cosa_v(npx,j)* &
2095  (ut(npx,j-1)+ut(npx+1,j-1)+ut(npx,j)+ut(npx+1,j))
2096  enddo
2097  endif
2098 
2099 ! South (Bottom) edge:
2100  if ( js==1 ) then
2101 
2102  do i=isd,ied
2103  if ( vc(i,1)*dt > 0. ) then
2104  vt(i,1) = vc(i,1) / geom%sin_sg(i,0,4)
2105  else
2106  vt(i,1) = vc(i,1) / geom%sin_sg(i,1,2)
2107  endif
2108  enddo
2109 
2110  do i=max(3,is),min(npx-2,ie+1)
2111  ut(i,0) = uc(i,0) - 0.25*geom%cosa_u(i,0)* &
2112  (vt(i-1,0)+vt(i,0)+vt(i-1,1)+vt(i,1))
2113  ut(i,1) = uc(i,1) - 0.25*geom%cosa_u(i,1)* &
2114  (vt(i-1,1)+vt(i,1)+vt(i-1,2)+vt(i,2))
2115  enddo
2116  endif
2117 
2118 ! North edge:
2119  if ( (je+1)==npy ) then
2120  do i=isd,ied
2121  if ( vc(i,npy)*dt > 0. ) then
2122  vt(i,npy) = vc(i,npy) / geom%sin_sg(i,npy-1,4)
2123  else
2124  vt(i,npy) = vc(i,npy) / geom%sin_sg(i,npy,2)
2125  endif
2126  enddo
2127  do i=max(3,is),min(npx-2,ie+1)
2128  ut(i,npy-1) = uc(i,npy-1) - 0.25*geom%cosa_u(i,npy-1)* &
2129  (vt(i-1,npy-1)+vt(i,npy-1)+vt(i-1,npy)+vt(i,npy))
2130  ut(i,npy) = uc(i,npy) - 0.25*geom%cosa_u(i,npy)* &
2131  (vt(i-1,npy)+vt(i,npy)+vt(i-1,npy+1)+vt(i,npy+1))
2132  enddo
2133  endif
2134 
2135 !The following code solves a 2x2 system to get the interior parallel-to-edge
2136 ! uc,vc values near the corners (ex: for the sw corner ut(2,1) and vt(1,2) are solved for simultaneously).
2137 ! It then computes the halo uc, vc values so as to be consistent with the computations on the facing panel.
2138 
2139  !The system solved is:
2140  ! ut(2,1) = uc(2,1) - avg(vt)*geom%cosa_u(2,1)
2141  ! vt(1,2) = vc(1,2) - avg(ut)*geom%cosa_v(1,2)
2142  ! in which avg(vt) includes vt(1,2) and avg(ut) includes ut(2,1)
2143 
2144  if( geom%sw_corner ) then
2145  damp = 1. / (1.-0.0625*geom%cosa_u(2,0)*geom%cosa_v(1,0))
2146  ut(2,0) = (uc(2,0)-0.25*geom%cosa_u(2,0)*(vt(1,1)+vt(2,1)+vt(2,0) +vc(1,0) - &
2147  0.25*geom%cosa_v(1,0)*(ut(1,0)+ut(1,-1)+ut(2,-1))) ) * damp
2148  damp = 1. / (1.-0.0625*geom%cosa_u(0,1)*geom%cosa_v(0,2))
2149  vt(0,2) = (vc(0,2)-0.25*geom%cosa_v(0,2)*(ut(1,1)+ut(1,2)+ut(0,2)+uc(0,1) - &
2150  0.25*geom%cosa_u(0,1)*(vt(0,1)+vt(-1,1)+vt(-1,2))) ) * damp
2151 
2152  damp = 1. / (1.-0.0625*geom%cosa_u(2,1)*geom%cosa_v(1,2))
2153  ut(2,1) = (uc(2,1)-0.25*geom%cosa_u(2,1)*(vt(1,1)+vt(2,1)+vt(2,2)+vc(1,2) - &
2154  0.25*geom%cosa_v(1,2)*(ut(1,1)+ut(1,2)+ut(2,2))) ) * damp
2155 
2156  vt(1,2) = (vc(1,2)-0.25*geom%cosa_v(1,2)*(ut(1,1)+ut(1,2)+ut(2,2)+uc(2,1) - &
2157  0.25*geom%cosa_u(2,1)*(vt(1,1)+vt(2,1)+vt(2,2))) ) * damp
2158  endif
2159 
2160  if( geom%se_corner ) then
2161  damp = 1. / (1. - 0.0625*geom%cosa_u(npx-1,0)*geom%cosa_v(npx-1,0))
2162  ut(npx-1,0) = ( uc(npx-1,0)-0.25*geom%cosa_u(npx-1,0)*( &
2163  vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,0)+vc(npx-1,0) - &
2164  0.25*geom%cosa_v(npx-1,0)*(ut(npx,0)+ut(npx,-1)+ut(npx-1,-1))) ) * damp
2165  damp = 1. / (1. - 0.0625*geom%cosa_u(npx+1,1)*geom%cosa_v(npx,2))
2166  vt(npx, 2) = ( vc(npx,2)-0.25*geom%cosa_v(npx,2)*( &
2167  ut(npx,1)+ut(npx,2)+ut(npx+1,2)+uc(npx+1,1) - &
2168  0.25*geom%cosa_u(npx+1,1)*(vt(npx,1)+vt(npx+1,1)+vt(npx+1,2))) ) * damp
2169 
2170  damp = 1. / (1. - 0.0625*geom%cosa_u(npx-1,1)*geom%cosa_v(npx-1,2))
2171  ut(npx-1,1) = ( uc(npx-1,1)-0.25*geom%cosa_u(npx-1,1)*( &
2172  vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,2)+vc(npx-1,2) - &
2173  0.25*geom%cosa_v(npx-1,2)*(ut(npx,1)+ut(npx,2)+ut(npx-1,2))) ) * damp
2174  vt(npx-1,2) = ( vc(npx-1,2)-0.25*geom%cosa_v(npx-1,2)*( &
2175  ut(npx,1)+ut(npx,2)+ut(npx-1,2)+uc(npx-1,1) - &
2176  0.25*geom%cosa_u(npx-1,1)*(vt(npx-1,1)+vt(npx-2,1)+vt(npx-2,2))) ) * damp
2177  endif
2178 
2179  if( geom%ne_corner ) then
2180  damp = 1. / (1. - 0.0625*geom%cosa_u(npx-1,npy)*geom%cosa_v(npx-1,npy+1))
2181  ut(npx-1,npy) = ( uc(npx-1,npy)-0.25*geom%cosa_u(npx-1,npy)*( &
2182  vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy+1)+vc(npx-1,npy+1) - &
2183  0.25*geom%cosa_v(npx-1,npy+1)*(ut(npx,npy)+ut(npx,npy+1)+ut(npx-1,npy+1))) ) * damp
2184  damp = 1. / (1. - 0.0625*geom%cosa_u(npx+1,npy-1)*geom%cosa_v(npx,npy-1))
2185  vt(npx, npy-1) = ( vc(npx,npy-1)-0.25*geom%cosa_v(npx,npy-1)*( &
2186  ut(npx,npy-1)+ut(npx,npy-2)+ut(npx+1,npy-2)+uc(npx+1,npy-1) - &
2187  0.25*geom%cosa_u(npx+1,npy-1)*(vt(npx,npy)+vt(npx+1,npy)+vt(npx+1,npy-1))) ) * damp
2188 
2189  damp = 1. / (1. - 0.0625*geom%cosa_u(npx-1,npy-1)*geom%cosa_v(npx-1,npy-1))
2190  ut(npx-1,npy-1) = ( uc(npx-1,npy-1)-0.25*geom%cosa_u(npx-1,npy-1)*( &
2191  vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy-1)+vc(npx-1,npy-1) - &
2192  0.25*geom%cosa_v(npx-1,npy-1)*(ut(npx,npy-1)+ut(npx,npy-2)+ut(npx-1,npy-2))) ) * damp
2193  vt(npx-1,npy-1) = ( vc(npx-1,npy-1)-0.25*geom%cosa_v(npx-1,npy-1)*( &
2194  ut(npx,npy-1)+ut(npx,npy-2)+ut(npx-1,npy-2)+uc(npx-1,npy-1) - &
2195  0.25*geom%cosa_u(npx-1,npy-1)*(vt(npx-1,npy)+vt(npx-2,npy)+vt(npx-2,npy-1))) ) * damp
2196  endif
2197 
2198  if( geom%nw_corner ) then
2199  damp = 1. / (1. - 0.0625*geom%cosa_u(2,npy)*geom%cosa_v(1,npy+1))
2200  ut(2,npy) = ( uc(2,npy)-0.25*geom%cosa_u(2,npy)*( &
2201  vt(1,npy)+vt(2,npy)+vt(2,npy+1)+vc(1,npy+1) - &
2202  0.25*geom%cosa_v(1,npy+1)*(ut(1,npy)+ut(1,npy+1)+ut(2,npy+1))) ) * damp
2203  damp = 1. / (1. - 0.0625*geom%cosa_u(0,npy-1)*geom%cosa_v(0,npy-1))
2204  vt(0,npy-1) = ( vc(0,npy-1)-0.25*geom%cosa_v(0,npy-1)*( &
2205  ut(1,npy-1)+ut(1,npy-2)+ut(0,npy-2)+uc(0,npy-1) - &
2206  0.25*geom%cosa_u(0,npy-1)*(vt(0,npy)+vt(-1,npy)+vt(-1,npy-1))) ) * damp
2207 
2208  damp = 1. / (1. - 0.0625*geom%cosa_u(2,npy-1)*geom%cosa_v(1,npy-1))
2209  ut(2,npy-1) = ( uc(2,npy-1)-0.25*geom%cosa_u(2,npy-1)*( &
2210  vt(1,npy)+vt(2,npy)+vt(2,npy-1)+vc(1,npy-1) - &
2211  0.25*geom%cosa_v(1,npy-1)*(ut(1,npy-1)+ut(1,npy-2)+ut(2,npy-2))) ) * damp
2212 
2213  vt(1,npy-1) = ( vc(1,npy-1)-0.25*geom%cosa_v(1,npy-1)*( &
2214  ut(1,npy-1)+ut(1,npy-2)+ut(2,npy-2)+uc(2,npy-1) - &
2215  0.25*geom%cosa_u(2,npy-1)*(vt(1,npy)+vt(2,npy)+vt(2,npy-1))) ) * damp
2216  endif
2217 
2218  else
2219 ! grid_type >= 3
2220 
2221  do j=jsd,jed
2222  do i=is,ie+1
2223  ut(i,j) = uc(i,j)
2224  enddo
2225  enddo
2226 
2227  do j=js,je+1
2228  do i=isd,ied
2229  vt(i,j) = vc(i,j)
2230  enddo
2231  enddo
2232  endif
2233 
2234 end subroutine compute_utvt
2235 
2236 ! ------------------------------------------------------------------------------
2237 
2238 subroutine ctoa(geom, uc, vc, ua, va)
2239 
2240  implicit none
2241  type(fv3jedi_geom), intent(in) :: geom
2242  real(kind=kind_real), intent(in) :: uc(geom%isd:geom%ied+1,geom%jsd:geom%jed ) ! c-grid u-wind field
2243  real(kind=kind_real), intent(in) :: vc(geom%isd:geom%ied ,geom%jsd:geom%jed+1) ! c-grid v-wind field
2244  real(kind=kind_real), intent(out) :: ua(geom%isd:geom%ied ,geom%jsd:geom%jed ) ! a-grid u-wind field
2245  real(kind=kind_real), intent(out) :: va(geom%isd:geom%ied ,geom%jsd:geom%jed ) ! a-grid v-wind field
2246 
2247  integer :: i,j
2248  real(kind=kind_real) :: tmp1i(geom%isd:geom%ied+1)
2249  real(kind=kind_real) :: tmp2i(geom%isd:geom%ied+1)
2250  real(kind=kind_real) :: tmp1j(geom%jsd:geom%jed+1)
2251  real(kind=kind_real) :: tmp2j(geom%jsd:geom%jed+1)
2252 
2253  do i=geom%isd,geom%ied
2254  tmp1j(:) = 0.0_kind_real
2255  tmp2j(:) = vc(i,:)*geom%dx(i,:)
2256  call interp_left_edge_1d(geom%jsd, geom%jed+1, tmp2j, tmp1j)
2257  va(i,geom%jsd:geom%jed) = tmp1j(geom%jsd+1:geom%jed+1)/geom%dxa(i,geom%jsd:geom%jed)
2258  enddo
2259 
2260  do j=geom%jsd,geom%jed
2261  tmp1i(:) = 0.0_kind_real
2262  tmp2i(:) = uc(:,j)*geom%dy(:,j)
2263  call interp_left_edge_1d(geom%isd, geom%ied+1, tmp2i, tmp1i)
2264  ua(geom%isd:geom%ied,j) = tmp1i(geom%isd+1:geom%ied+1)/geom%dya(geom%isd:geom%ied,j)
2265  enddo
2266 
2267 end subroutine ctoa
2268 
2269 ! ------------------------------------------------------------------------------
2270 
2271 subroutine atod(geom, ua, va, u, v)
2272 
2273  implicit none
2274  type(fv3jedi_geom), intent(in) :: geom
2275  real(kind=kind_real), intent(in) :: ua(geom%isd:geom%ied ,geom%jsd:geom%jed ) ! a-grid u-wind field
2276  real(kind=kind_real), intent(in) :: va(geom%isd:geom%ied ,geom%jsd:geom%jed ) ! a-grid v-wind field
2277  real(kind=kind_real), intent(out) :: u(geom%isd:geom%ied ,geom%jsd:geom%jed+1) ! d-grid u-wind field
2278  real(kind=kind_real), intent(out) :: v(geom%isd:geom%ied+1,geom%jsd:geom%jed ) ! d-grid v-wind field
2279 
2280  integer :: i,j
2281  real(kind=kind_real) :: tmp1i(geom%isd:geom%ied+1)
2282  real(kind=kind_real) :: tmp2i(geom%isd:geom%ied)
2283  real(kind=kind_real) :: tmp1j(geom%jsd:geom%jed+1)
2284  real(kind=kind_real) :: tmp2j(geom%jsd:geom%jed)
2285 
2286  do j=geom%jsd+1,geom%jed
2287  tmp1i(:) = 0.0
2288  tmp2i(:) = va(:,j)*geom%dxa(:,j)
2289  call interp_left_edge_1d(geom%isd, geom%ied, tmp1i, tmp2i)
2290  v(:,j) = tmp1i(:)/geom%dxc(:,j)
2291  enddo
2292  do i=geom%isd+1,geom%ied
2293  tmp1j(:) = 0.0
2294  tmp2j(:) = ua(i,:)*geom%dya(i,:)
2295  call interp_left_edge_1d(geom%jsd, geom%jed, tmp1j, tmp2j)
2296  u(i,:) = tmp1j(:)/geom%dyc(i,:)
2297  enddo
2298 
2299 end subroutine atod
2300 
2301 ! ------------------------------------------------------------------------------
2302 
2303 subroutine interp_left_edge_1d(ifirst, ilast, qin, qout)
2304 
2305  implicit none
2306  integer, intent(in) :: ifirst,ilast
2307  real(kind=kind_real), intent(in) :: qin(ifirst:ilast)
2308  real(kind=kind_real), intent(out) :: qout(ifirst:ilast)
2309 
2310  integer :: i
2311 
2312  qout(:) = 0.0_kind_real
2313  do i=ifirst+1,ilast
2314  qout(i) = 0.5_kind_real * (qin(i-1) + qin(i))
2315  enddo
2316 
2317 end subroutine interp_left_edge_1d
2318 
2319 ! ------------------------------------------------------------------------------
2320 
2321 end module wind_vt_mod
wind_vt_mod::fill_dgrid_winds
subroutine fill_dgrid_winds(geom, u, v, fillhalo)
Definition: wind_variables_mod.f90:1902
wind_vt_mod::interp_left_edge_1d
subroutine interp_left_edge_1d(ifirst, ilast, qin, qout)
Definition: wind_variables_mod.f90:2304
wind_vt_mod::edge_interpolate4
real(kind=kind_real) function edge_interpolate4(ua, dxa)
Definition: wind_variables_mod.f90:1006
fv3jedi_constants_mod::rad2deg
real(kind=kind_real), parameter, public rad2deg
Definition: fv3jedi_constants_mod.f90:13
wind_vt_mod::compute_utvt
subroutine compute_utvt(geom, uc, vc, ut, vt, dt)
Definition: wind_variables_mod.f90:2019
fv3jedi_geom_mod
Fortran module handling geometry for the FV3 model.
Definition: fv3jedi_geom_mod.f90:8
wind_vt_mod::sfc_10m_winds
subroutine, public sfc_10m_winds(geom, usrf, vsrf, f10r, spd10m, dir10m)
Definition: wind_variables_mod.f90:48
wind_vt_mod::psichi_to_vortdivg
subroutine, public psichi_to_vortdivg(geom, grid, oprs, psi, chi, lsize, lev_start, lev_final, vor, div)
Definition: wind_variables_mod.f90:522
fv3jedi_communication_mod
Definition: fv3jedi_communication_mod.f90:1
fv3jedi_geom_mod::fv3jedi_geom
Fortran derived type to hold geometry data for the FV3JEDI model.
Definition: fv3jedi_geom_mod.f90:46
fv3jedi_communication_mod::scatter_field
subroutine, public scatter_field(geom, comm, gproc, field_in, field_out)
Definition: fv3jedi_communication_mod.f90:116
wind_vt_mod::d2a_ad
subroutine, public d2a_ad(geom, u_ad_comp, v_ad_comp, ua_ad_comp, va_ad_comp)
Definition: wind_variables_mod.f90:1607
wind_vt_mod
Definition: wind_variables_mod.f90:6
wind_vt_mod::udvd_to_psichi
subroutine, public udvd_to_psichi(geom, grid, oprs, u_in, v_in, psi, chi, lsize, lev_start, lev_final, vor_out, div_out)
Definition: wind_variables_mod.f90:321
fv3jedi_communication_mod::gather_field
subroutine, public gather_field(geom, comm, gproc, field_in, field_out)
Definition: fv3jedi_communication_mod.f90:19
wind_vt_mod::atod
subroutine atod(geom, ua, va, u, v)
Definition: wind_variables_mod.f90:2272
fv3jedi_constants_mod
Definition: fv3jedi_constants_mod.f90:6
wind_vt_mod::ctoa
subroutine ctoa(geom, uc, vc, ua, va)
Definition: wind_variables_mod.f90:2239
wind_vt_mod::psichi_to_uava
subroutine, public psichi_to_uava(geom, psi, chi, ua, va)
Definition: wind_variables_mod.f90:111
wind_vt_mod::psichi_to_uava_adm
subroutine, public psichi_to_uava_adm(geom, psi_ad, chi_ad, ua_ad, va_ad)
Definition: wind_variables_mod.f90:264
fv3jedi_kinds_mod::kind_real
integer, parameter, public kind_real
Definition: fv3jedi_kinds_mod.f90:14
wind_vt_mod::d2a
subroutine, public d2a(geom, u_comp, v_comp, ua_comp, va_comp)
Definition: wind_variables_mod.f90:1456
wind_vt_mod::a2d
subroutine, public a2d(geom, ua, va, ud, vd)
Definition: wind_variables_mod.f90:1023
fv3jedi_kinds_mod
Definition: fv3jedi_kinds_mod.f90:6
wind_vt_mod::psichi_to_udvd
subroutine, public psichi_to_udvd(geom, psi, chi, u, v)
Definition: wind_variables_mod.f90:153
fv3jedi_constants_mod::pi
real(kind=kind_real), parameter, public pi
Definition: fv3jedi_constants_mod.f90:16
wind_vt_mod::udvd_to_ucvc
subroutine udvd_to_ucvc(geom, u, v, uc, vc, ua_out, va_out)
Definition: wind_variables_mod.f90:629
wind_vt_mod::a2d_ad
subroutine, public a2d_ad(geom, ua_ad, va_ad, ud_ad, vd_ad)
Definition: wind_variables_mod.f90:1194
wind_vt_mod::fill_cgrid_winds
subroutine fill_cgrid_winds(geom, uc, vc, fillhalo)
Definition: wind_variables_mod.f90:1961