OOPS
qg_advect_q_mod.F90
Go to the documentation of this file.
1 ! (C) Copyright 2009-2016 ECMWF.
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 ! In applying this licence, ECMWF does not waive the privileges and immunities
6 ! granted to it by virtue of its status as an intergovernmental organisation nor
7 ! does it submit to any jurisdiction.
8 
9 module qg_advect_q_mod
10 
11 use kinds
12 !$ use omp_lib
13 use qg_constants_mod
14 use qg_geom_mod
15 use qg_interp_mod
16 
17 implicit none
18 
19 private
20 public :: advect_q,advect_q_tl,advect_q_ad
21 ! ------------------------------------------------------------------------------
22 contains
23 ! ------------------------------------------------------------------------------
24 !> Advect potential vorticity
25 subroutine advect_q(geom,dt,u,v,q,q_north,q_south,qnew)
26 
27 implicit none
28 
29 ! Passed variables
30 type(qg_geom),intent(in) :: geom !< Geometry
31 real(kind_real),intent(in) :: dt !< Time step
32 real(kind_real),intent(in) :: u(geom%nx,geom%ny,geom%nz) !< Zonal wind
33 real(kind_real),intent(in) :: v(geom%nx,geom%ny,geom%nz) !< Meridional wind
34 real(kind_real),intent(in) :: q(geom%nx,geom%ny,geom%nz) !< Input potential vorticity
35 real(kind_real),intent(in) :: q_north(geom%nx,geom%nz) !< Potential vorticity on northern wall
36 real(kind_real),intent(in) :: q_south(geom%nx,geom%nz) !< Potential vorticity on southern wall
37 real(kind_real),intent(out) :: qnew(geom%nx,geom%ny,geom%nz) !< Output potential vorticity
38 
39 ! Local variables
40 integer :: ix,iy,iz
41 real(kind_real) :: x,y
42 real(kind_real) :: qext(geom%nx,-1:geom%ny+2,geom%nz)
43 
44 ! Extend q field
45 qext(:,1:geom%ny,:) = q
46 qext(:,-1,:) = q_south
47 qext(:,0,:) = q_south
48 qext(:,geom%ny+1,:) = q_north
49 qext(:,geom%ny+2,:) = q_north
50 
51 ! Advect q
52 !$omp parallel do schedule(static) private(iz,iy,ix,x,y)
53 do iz=1,geom%nz
54  do iy=1,geom%ny
55  do ix=1,geom%nx
56  ! Find the interpolation point
57  x = geom%x(ix)-u(ix,iy,iz)*dt
58  y = geom%y(iy)-v(ix,iy,iz)*dt
59 
60  ! Interpolate
61  call qg_interp_bicubic(geom,x,y,qext(:,:,iz),qnew(ix,iy,iz))
62  enddo
63  enddo
64 enddo
65 !$omp end parallel do
66 
67 end subroutine advect_q
68 ! ------------------------------------------------------------------------------
69 !> Advect potential vorticity - tangent linear
70 subroutine advect_q_tl(geom,dt,u_traj,v_traj,q_traj,q_traj_north,q_traj_south,u,v,q,qnew)
71 
72 implicit none
73 
74 ! Passed variables
75 type(qg_geom),intent(in) :: geom !< Geometry
76 real(kind_real),intent(in) :: dt !< Time step
77 real(kind_real),intent(in) :: u_traj(geom%nx,geom%ny,geom%nz) !< Zonal wind (trajectory)
78 real(kind_real),intent(in) :: v_traj(geom%nx,geom%ny,geom%nz) !< Meridional wind (trajectory)
79 real(kind_real),intent(in) :: q_traj(geom%nx,geom%ny,geom%nz) !< Potential vorticity (trajectory)
80 real(kind_real),intent(in) :: q_traj_north(geom%nx,geom%nz) !< Potential vorticity on northern wall (trajectory)
81 real(kind_real),intent(in) :: q_traj_south(geom%nx,geom%nz) !< Potential vorticity on southern wall (trajectory)
82 real(kind_real),intent(in) :: u(geom%nx,geom%ny,geom%nz) !< Zonal wind (perturbation)
83 real(kind_real),intent(in) :: v(geom%nx,geom%ny,geom%nz) !< Meridional wind (perturbation)
84 real(kind_real),intent(in) :: q(geom%nx,geom%ny,geom%nz) !< Input potential vorticity (perturbation)
85 real(kind_real),intent(out) :: qnew(geom%nx,geom%ny,geom%nz) !< Output potential vorticity (perturbation)
86 
87 ! Local variables
88 integer :: ix,iy,iz
89 real(kind_real) :: x,y,dx,dy
90 real(kind_real) :: qext_traj(geom%nx,-1:geom%ny+2,geom%nz),qext(geom%nx,-1:geom%ny+2,geom%nz)
91 
92 ! Extend q field (trajectory)
93 qext_traj(:,1:geom%ny,:) = q_traj
94 qext_traj(:,-1,:) = q_traj_south
95 qext_traj(:,0,:) = q_traj_south
96 qext_traj(:,geom%ny+1,:) = q_traj_north
97 qext_traj(:,geom%ny+2,:) = q_traj_north
98 
99 ! Extend q field (perturbation)
100 qext(:,1:geom%ny,:) = q
101 qext(:,-1,:) = 0.0
102 qext(:,0,:) = 0.0
103 qext(:,geom%ny+1,:) = 0.0
104 qext(:,geom%ny+2,:) = 0.0
105 
106 ! Advect q
107 !$omp parallel do schedule(static) private(iz,iy,ix,x,y,dx,dy)
108 do iz=1,geom%nz
109  do iy=1,geom%ny
110  do ix=1,geom%nx
111  ! Find the interpolation point
112  x = geom%x(ix)-u_traj(ix,iy,iz)*dt
113  y = geom%y(iy)-v_traj(ix,iy,iz)*dt
114 
115  ! Find interpolation point perturbation
116  dx = -u(ix,iy,iz)*dt
117  dy = -v(ix,iy,iz)*dt
118 
119  ! Interpolate
120  call qg_interp_bicubic_tl(geom,x,y,qext_traj(:,:,iz),dx,dy,qext(:,:,iz),qnew(ix,iy,iz))
121  enddo
122  enddo
123 enddo
124 !$omp end parallel do
125 
126 end subroutine advect_q_tl
127 ! ------------------------------------------------------------------------------
128 !> Advect potential vorticity - adjoint
129 subroutine advect_q_ad(geom,dt,u_traj,v_traj,q_traj,q_traj_north,q_traj_south,qnew,u,v,q)
130 
131 implicit none
132 
133 ! Passed variables
134 type(qg_geom),intent(in) :: geom !< Geometry
135 real(kind_real),intent(in) :: dt !< Time step
136 real(kind_real),intent(in) :: u_traj(geom%nx,geom%ny,geom%nz) !< Zonal wind (trajectory)
137 real(kind_real),intent(in) :: v_traj(geom%nx,geom%ny,geom%nz) !< Meridional wind (trajectory)
138 real(kind_real),intent(in) :: q_traj(geom%nx,geom%ny,geom%nz) !< Potential vorticity (trajectory)
139 real(kind_real),intent(in) :: q_traj_north(geom%nx,geom%nz) !< Potential vorticity on northern wall (trajectory)
140 real(kind_real),intent(in) :: q_traj_south(geom%nx,geom%nz) !< Potential vorticity on southern wall (trajectory)
141 real(kind_real),intent(in) :: qnew(geom%nx,geom%ny,geom%nz) !< Output potential vorticity (perturbation)
142 real(kind_real),intent(inout) :: u(geom%nx,geom%ny,geom%nz) !< Zonal wind (perturbation)
143 real(kind_real),intent(inout) :: v(geom%nx,geom%ny,geom%nz) !< Meridional wind (perturbation)
144 real(kind_real),intent(inout) :: q(geom%nx,geom%ny,geom%nz) !< Input potential vorticity (perturbation)
145 
146 ! Local variables
147 integer :: ix,iy,iz
148 real(kind_real) :: x,y,dx,dy
149 real(kind_real) :: qext_traj(geom%nx,-1:geom%ny+2,geom%nz),qext(geom%nx,-1:geom%ny+2,geom%nz)
150 
151 ! Extend q field (trajectory)
152 qext_traj(:,1:geom%ny,:) = q_traj
153 qext_traj(:,-1,:) = q_traj_south
154 qext_traj(:,0,:) = q_traj_south
155 qext_traj(:,geom%ny+1,:) = q_traj_north
156 qext_traj(:,geom%ny+2,:) = q_traj_north
157 
158 ! Initialization
159 qext = 0.0
160 
161 ! Advect q
162 do iz=geom%nz,1,-1
163  do iy=geom%ny,1,-1
164  do ix=geom%nx,1,-1
165  ! Find the interpolation point
166  x = geom%x(ix)-u_traj(ix,iy,iz)*dt
167  y = geom%y(iy)-v_traj(ix,iy,iz)*dt
168 
169  ! Initialization
170  dx = 0.0
171  dy = 0.0
172 
173  ! Interpolate,adjoint
174  call qg_interp_bicubic_ad(geom,x,y,qext_traj(:,:,iz),qnew(ix,iy,iz),dx,dy,qext(:,:,iz))
175 
176  ! Find interpolation point perturbation,adjoint
177  u(ix,iy,iz) = u(ix,iy,iz)-dx*dt
178  v(ix,iy,iz) = v(ix,iy,iz)-dy*dt
179  enddo
180  enddo
181 enddo
182 
183 ! Extend q field (perturbation)
184 q = q+qext(:,1:geom%ny,:)
185 
186 end subroutine advect_q_ad
187 ! ------------------------------------------------------------------------------
188 end module qg_advect_q_mod
qg_advect_q_mod
Definition: qg_advect_q_mod.F90:9
qg_advect_q_mod::advect_q_tl
subroutine, public advect_q_tl(geom, dt, u_traj, v_traj, q_traj, q_traj_north, q_traj_south, u, v, q, qnew)
Advect potential vorticity - tangent linear.
Definition: qg_advect_q_mod.F90:71
qg_advect_q_mod::advect_q_ad
subroutine, public advect_q_ad(geom, dt, u_traj, v_traj, q_traj, q_traj_north, q_traj_south, qnew, u, v, q)
Advect potential vorticity - adjoint.
Definition: qg_advect_q_mod.F90:130
qg_advect_q_mod::advect_q
subroutine, public advect_q(geom, dt, u, v, q, q_north, q_south, qnew)
Advect potential vorticity.
Definition: qg_advect_q_mod.F90:26
qg_constants_mod
Definition: qg_constants_mod.F90:9
qg_geom_mod
Definition: qg_geom_mod.F90:9
qg_interp_mod
Definition: qg_interp_mod.F90:9
qg_interp_mod::qg_interp_bicubic_ad
subroutine, public qg_interp_bicubic_ad(geom, x, y, gfld2dext_traj, val, dx, dy, gfld2dext)
Bicubic horizontal interpolation - adjoint.
Definition: qg_interp_mod.F90:221
qg_interp_mod::qg_interp_bicubic_tl
subroutine, public qg_interp_bicubic_tl(geom, x, y, gfld2dext_traj, dx, dy, gfld2dext, val)
Bicubic horizontal interpolation - tangent linear.
Definition: qg_interp_mod.F90:174
qg_interp_mod::qg_interp_bicubic
subroutine, public qg_interp_bicubic(geom, x, y, gfld2dext, val)
Bicubic horizontal interpolation.
Definition: qg_interp_mod.F90:144
qg_geom_mod::qg_geom
Definition: qg_geom_mod.F90:26