module_cnvdiff.f90
Go to the documentation of this file.
1 !*************************************************************************
2 ! This file is part of HERCULES.
3 !
4 ! HERCULES is free software; you can redistribute it and/or modify
5 ! it under the terms of the GNU General Public License as published by
6 ! the Free Software Foundation; either version 3 of the License, or
7 ! (at your option) any later version.
8 !
9 ! HERCULES is distributed in the hope that it will be useful,
10 ! but WITHOUT ANY WARRANTY; without even the implied warranty of
11 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 ! GNU General Public License for more details.
13 !
14 ! You should have received a copy of the GNU General Public License
15 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
16 !
17 ! Copyright 2015, Ping He, MEAS, NCSU
18 !*************************************************************************
19 
20 
21 
23 
25 
26  implicit none
27 
28 contains
29 
30 
31 
32  ! this function calculates the spatial derivatives and
33  ! the right hand side terms for the momentum equations
34  ! i.e., f,g,h,qu,qv,qw etc
35  subroutine get_momentum_cnvdiff &
36  (u,v,w,p,t,uhx,uhy,uhz,vhx,vhy,vhz,whx,why,ff,ff0,gg,gg0,hh,hh0,qu,qv,qw)
37 
38  implicit none
39 
40  real(mytype),dimension(:,:,:),intent(in) ::u,v,w,p,t,ff0,gg0, &
41  hh0,uhx,uhy,uhz,vhx,vhy,vhz,whx,why
42  real(mytype),dimension(:,:,:),intent(out)::ff,gg,hh,qu,qv,qw
43 
44  integer :: i,j,k
45  real(mytype) :: a1,a2
46  real(mytype) :: uux, uvy, uwz, uxx, uyy, uzz, dpdx, ux
47  real(mytype) :: vux, vvy, vwz, vxx, vyy, vzz, dpdy, vx
48  real(mytype) :: wux, wvy, wwz, wxx, wyy, wzz, dpdz, wx
49 
50  ! first calculate dpx_drive, it is value is dependent on dp_opt
51  call update_dpx_drive
52 
53 
54  !!!!!!!!!! U equation !!!!!!!!!!!
55  do k=kstr3,kend3,1
56  do j=jstr3,jend3,1
57  do i=istr3,iend3,1
58 
59  ! horizontal derivatives, check to use 2nd or 4th order schemes
60  if ( cds .eq. 1) then
61 
62  uux= ( uhx(i+1,j,k)**2.d0 - &
63  uhx(i,j,k)**2.d0 ) / dx
64 
65  uvy= ( uhy(i,j,k)*vhx(i,j,k) - &
66  uhy(i,j-1,k)*vhx(i,j-1,k) ) / dy
67 
68  ux = ( uhx(i+1,j,k)-uhx(i,j,k) ) / dx
69 
70  elseif ( cds .eq. 2) then
71 
72  uux= ( -uhx(i+2,j,k)**2.d0 + &
73  27.d0*uhx(i+1,j,k)**2.d0 - &
74  27.d0*uhx(i,j,k)**2.d0 + &
75  uhx(i-1,j,k)**2.d0 ) / 24.d0 / dx
76 
77  uvy= ( -uhy(i,j+1,k)*vhx(i,j+1,k) + &
78  27.d0*uhy(i,j,k) *vhx(i,j,k) - &
79  27.d0*uhy(i,j-1,k)*vhx(i,j-1,k) + &
80  uhy(i,j-2,k)*vhx(i,j-2,k) ) / 24.d0 /dy
81 
82  ux = ( -uhx(i+2,j,k) + &
83  27.d0*uhx(i+1,j,k) - &
84  27.d0*uhx(i,j,k) + &
85  uhx(i-1,j,k) ) / 24.d0 / dx
86 
87  endif
88 
89  ! we always use the 2nd order scheme for vertical derivatives
90  uwz= ( uhz(i,j,k) *whx(i,j,k) - &
91  uhz(i,j,k-1)*whx(i,j,k-1) ) / dz(k)
92 
93 
94  ! horizontal derivatives, check to use 2nd or 4th schemes
95  if ( cds .eq. 1) then
96 
97  uxx= ( u(i+1,j,k) - &
98  2.d0*u(i,j,k) + &
99  u(i-1,j,k) ) / dx2
100 
101  uyy= ( u(i,j+1,k) - &
102  2.d0*u(i,j,k) + &
103  u(i,j-1,k) ) / dy2
104 
105  elseif ( cds .eq. 2) then
106 
107  uxx= ( -u(i+2,j,k) + &
108  16.d0*u(i+1,j,k) - &
109  30.d0*u(i,j,k) + &
110  16.d0*u(i-1,j,k) - &
111  u(i-2,j,k) ) / 12.d0 / dx2
112 
113  uyy= ( -u(i,j+2,k) + &
114  16.d0*u(i,j+1,k) - &
115  30.d0*u(i,j,k) + &
116  16.d0*u(i,j-1,k) - &
117  u(i,j-2,k) ) / 12.d0 / dy2
118 
119  endif
120 
121  ! we always use the 2nd order scheme for vertical derivatives
122  uzz= ( (u(i,j,k+1)-u(i,j,k))/dz_t(k) - &
123  (u(i,j,k)-u(i,j,k-1))/dz_b(k) ) / dz(k)
124 
125  ! calculate the nonlinear terms for U
126  ! check if it is Ekman layer, if yes, add the Coriolis term
127  if (dp_opt .eq. 3) then
128 
129  ! constant geostrophic wind
130  ff(i,j,k) = -uux-uvy-uwz - u_mrf*ux + &
131  fc*( 0.5d0*(vhx(i,j,k)+vhx(i,j-1,k)) - vg )
132 
133  else
134 
135  ! + dpx_drive to drive the flow for dp_opt=1 or 2
136  ff(i,j,k) = -uux-uvy-uwz + dpx_drive - u_mrf*ux
137 
138  endif
139 
140  ! dpdx, check to use 2nd or 4th schemes
141  if ( cds .eq. 1) then
142 
143  dpdx= ( p(i+1,j,k)-p(i,j,k) ) / dx
144 
145  elseif ( cds .eq. 2) then
146 
147  dpdx= ( -p(i+2,j,k) + &
148  27.d0*p(i+1,j,k) - &
149  27.d0*p(i,j,k) + &
150  p(i-1,j,k) ) / 24.d0 / dx
151 
152  endif
153 
154  ! the right hand side term for U eqn
155  qu(i,j,k)= u(i,j,k) + dt*( rk_a*ff(i,j,k)+rk_b*ff0(i,j,k) ) + &
156  rk_c*0.5d0*dt*nu*(uxx+uyy+uzz) - rk_c*dt*dpdx
157 
158 
159  enddo
160  enddo
161  enddo
162 
163 
164 
165  !!!!!!!!!! V equation !!!!!!!!!!!
166  do k=kstr3,kend3,1
167  do j=jstr3,jend3,1
168  do i=istr3,iend3,1
169 
170  ! horizontal derivatives, check to use 2nd or 4th schemes
171  if ( cds .eq. 1) then
172 
173  vux= ( vhx(i,j,k) *uhy(i,j,k) - &
174  vhx(i-1,j,k)*uhy(i-1,j,k) ) / dx
175 
176  vvy= ( vhy(i,j+1,k)**2.d0 - &
177  vhy(i,j,k)**2.d0 ) / dy
178 
179  vx = ( vhx(i,j,k)-vhx(i-1,j,k) ) / dx
180 
181  elseif ( cds .eq. 2) then
182 
183  vux= ( -vhx(i+1,j,k)*uhy(i+1,j,k) + &
184  27.d0*vhx(i,j,k) *uhy(i,j,k) - &
185  27.d0*vhx(i-1,j,k)*uhy(i-1,j,k) + &
186  vhx(i-2,j,k)*uhy(i-2,j,k) ) / 24.d0 / dx
187 
188  vvy= ( -vhy(i,j+2,k)**2.d0 + &
189  27.d0*vhy(i,j+1,k)**2.d0 - &
190  27.d0*vhy(i,j,k)**2.d0 + &
191  vhy(i,j-1,k)**2.d0 ) / 24.d0 / dy
192 
193  vx = ( -vhx(i+1,j,k) + &
194  27.d0*vhx(i,j,k) - &
195  27.d0*vhx(i-1,j,k) + &
196  vhx(i-2,j,k) ) / 24.d0 / dx
197 
198  endif
199 
200  ! we always use the 2nd order scheme for vertical derivatives
201  vwz= ( vhz(i,j,k) *why(i,j,k) - &
202  vhz(i,j,k-1)*why(i,j,k-1) ) / dz(k)
203 
204  ! horizontal derivatives, check to use 2nd or 4th schemes
205  if ( cds .eq. 1) then
206 
207  vxx= ( v(i+1,j,k) - &
208  2.d0*v(i,j,k) + &
209  v(i-1,j,k) ) / dx2
210 
211  vyy= ( v(i,j+1,k) - &
212  2.d0*v(i,j,k) + &
213  v(i,j-1,k) ) / dy2
214 
215  elseif ( cds .eq. 2) then
216 
217  vxx= ( -v(i+2,j,k) + &
218  16.d0*v(i+1,j,k) - &
219  30.d0*v(i,j,k) + &
220  16.d0*v(i-1,j,k) - &
221  v(i-2,j,k) ) / 12.d0 / dx2
222 
223  vyy= ( -v(i,j+2,k) + &
224  16.d0*v(i,j+1,k) - &
225  30.d0*v(i,j,k) + &
226  16.d0*v(i,j-1,k) - &
227  v(i,j-2,k) ) / 12.d0 / dy2
228 
229  endif
230 
231  ! we always use the 2nd order scheme for vertical derivatives
232  vzz= ( (v(i,j,k+1)-v(i,j,k))/dz_t(k) - &
233  (v(i,j,k)-v(i,j,k-1))/dz_b(k) ) / dz(k)
234 
235  ! calculate the nonlinear term for V
236  ! check if it is Ekman layer
237  if (dp_opt .eq. 3) then
238 
239  ! constant geostrophci wind
240  gg(i,j,k) = -vux-vvy-vwz - u_mrf*vx + &
241  fc*( ug - 0.5d0*(uhy(i,j,k)+uhy(i-1,j,k)) - u_mrf )
242 
243  else
244 
245  gg(i,j,k) = -vux-vvy-vwz - u_mrf*vx
246 
247  endif
248 
249  ! dpdy, check to use 2nd or 4th schemes
250  if ( cds .eq. 1) then
251 
252  dpdy = ( p(i,j+1,k)-p(i,j,k) ) / dy
253 
254  elseif ( cds .eq. 2) then
255 
256  dpdy= ( -p(i,j+2,k) + &
257  27.d0*p(i,j+1,k) - &
258  27.d0*p(i,j,k) + &
259  p(i,j-1,k) ) / 24.d0 / dy
260 
261  endif
262 
263  ! the right hand side term for V eqn
264  qv(i,j,k)=v(i,j,k) + dt*( rk_a*gg(i,j,k)+rk_b*gg0(i,j,k) ) + &
265  rk_c*0.5d0*dt*nu*(vxx+vyy+vzz) - rk_c*dt*dpdy
266 
267 
268  enddo
269  enddo
270  enddo
271 
272 
273  !!!!!!!!!!!!!!!!!!!!!!!! W equation !!!!!!!!!!!!!!!!!
274  do k=kstr3,kend3-1,1 ! note that here the z range is different
275  do j=jstr3,jend3,1
276  do i=istr3,iend3,1
277 
278  ! horizontal derivatives, check to use 2nd or 4th schemes
279  if ( cds .eq. 1) then
280 
281  wux= ( whx(i,j,k) *uhz(i,j,k) - &
282  whx(i-1,j,k)*uhz(i-1,j,k) ) / dx
283 
284  wvy= ( why(i,j,k) *vhz(i,j,k) - &
285  why(i,j-1,k)*vhz(i,j-1,k) ) / dy
286 
287  wx = ( whx(i,j,k)-whx(i-1,j,k) ) / dx
288 
289  elseif ( cds .eq. 2) then
290 
291  wux= ( -whx(i+1,j,k)*uhz(i+1,j,k) + &
292  27.d0*whx(i,j,k) *uhz(i,j,k) - &
293  27.d0*whx(i-1,j,k)*uhz(i-1,j,k) + &
294  whx(i-2,j,k)*uhz(i-2,j,k) ) / 24.d0 /dx
295 
296  wvy= ( -why(i,j+1,k)*vhz(i,j+1,k) + &
297  27.d0*why(i,j,k) *vhz(i,j,k) - &
298  27.d0*why(i,j-1,k)*vhz(i,j-1,k) + &
299  why(i,j-2,k)*vhz(i,j-2,k) ) / 24.d0 /dy
300 
301  wx = ( -whx(i+1,j,k) + &
302  27.d0*whx(i,j,k) - &
303  27.d0*whx(i-1,j,k) + &
304  whx(i-2,j,k) ) / 24.d0 /dx
305 
306  endif
307 
308  ! we always use the 2nd order scheme for vertical derivatives
309  a1 = ( 0.5d0*(w(i,j,k)+w(i,j,k+1)) )**2.d0
310  a2 = ( 0.5d0*(w(i,j,k)+w(i,j,k-1)) )**2.d0
311  wwz= (a1-a2) / dz_t(k)
312 
313  ! horizontal derivatives, check to use 2nd or 4th schemes
314  if ( cds .eq. 1) then
315 
316  wxx= ( w(i+1,j,k) - &
317  2.d0*w(i,j,k) + &
318  w(i-1,j,k) ) / dx2
319 
320  wyy= ( w(i,j+1,k) - &
321  2.d0*w(i,j,k) + &
322  w(i,j-1,k) ) / dy2
323 
324  elseif ( cds .eq. 2) then
325 
326  wxx= ( -w(i+2,j,k) + &
327  16.d0*w(i+1,j,k) - &
328  30.d0*w(i,j,k) + &
329  16.d0*w(i-1,j,k) - &
330  w(i-2,j,k) ) / 12.d0 / dx2
331 
332  wyy= ( -w(i,j+2,k) + &
333  16.d0*w(i,j+1,k) - &
334  30.d0*w(i,j,k) + &
335  16.d0*w(i,j-1,k) - &
336  w(i,j-2,k) ) / 12.d0 / dy2
337 
338  endif
339 
340  ! we always use the 2nd order scheme for vertical derivatives
341  wzz= ( (w(i,j,k+1)-w(i,j,k))/dz(k+1) - &
342  (w(i,j,k)-w(i,j,k-1))/dz(k) ) / dz_t(k)
343 
344  ! calculate the nonlinear term for W
345  ! check whether to add the buoyancy term
346  if (isscalar .eq. 1) then
347 
348  hh(i,j,k) = -wux-wvy-wwz -u_mrf*wx + &
349  got*( (1.d0-l_t(k))*t(i,j,k)+l_t(k)*t(i,j,k+1) -t_ref )
350 
351  elseif (isscalar .eq. 0) then
352 
353  hh(i,j,k) = -wux-wvy-wwz -u_mrf*wx
354 
355  else
356  print *, 'isscalar error!'
357  stop
358  endif
359 
360  ! dpdz term
361  dpdz = ( p(i,j,k+1)-p(i,j,k) ) / dz_t(k)
362 
363  ! the right hand side term for W eqn
364  qw(i,j,k)=w(i,j,k) + dt*( rk_a*hh(i,j,k)+rk_b*hh0(i,j,k) ) + &
365  rk_c*0.5d0*dt*nu*(wxx+wyy+wzz) - rk_c*dt*dpdz
366 
367 
368  enddo
369  enddo
370  enddo
371 
372 
373  return
374  end subroutine
375 
376 
377 
378  ! this function calculates the spatial derivatives and
379  ! the right hand side terms for the temperature equation
380  ! i.e., ss, qt etc
381  subroutine get_temperature_cnvdiff(u,v,w,t,thx,thy,thz,ss,ss0,qt)
383  implicit none
384 
385  real(mytype),dimension(:,:,:),intent(in) :: u,v,w,t,ss0,thx,thy,thz
386  real(mytype),dimension(:,:,:),intent(out) :: ss,qt
387 
388  integer :: i,j,k
389  real(mytype) :: tux, tvy, twz, txx, tyy, tzz, tx
390 
391  do k=kstr3,kend3,1
392  do j=jstr3,jend3,1
393  do i=istr3,iend3,1
394 
395  ! horizontal derivatives, check to use 2nd or 4th schemes
396  if ( cds .eq. 1) then
397 
398  tux= ( thx(i,j,k) *u(i,j,k) - &
399  thx(i-1,j,k)*u(i-1,j,k) ) / dx
400 
401  tvy= ( thy(i,j,k) *v(i,j,k) - &
402  thy(i,j-1,k)*v(i,j-1,k) ) / dy
403 
404  tx = ( thx(i,j,k)-thx(i-1,j,k) ) / dx
405 
406  elseif ( cds .eq. 2) then
407 
408  tux= ( -thx(i+1,j,k)*u(i+1,j,k) + &
409  27.d0*thx(i,j,k) *u(i,j,k) - &
410  27.d0*thx(i-1,j,k)*u(i-1,j,k) + &
411  thx(i-2,j,k)*u(i-2,j,k) ) / 24.d0 / dx
412 
413  tvy= ( -thy(i,j+1,k)*v(i,j+1,k) + &
414  27.d0*thy(i,j,k) *v(i,j,k) - &
415  27.d0*thy(i,j-1,k)*v(i,j-1,k) + &
416  thy(i,j-2,k)*v(i,j-2,k) ) / 24.d0 / dy
417 
418  tx = ( -thx(i+1,j,k) + &
419  27.d0*thx(i,j,k) - &
420  27.d0*thx(i-1,j,k) + &
421  thx(i-2,j,k) ) / 24.d0 / dx
422 
423  endif
424 
425  ! we always use the 2nd order scheme for vertical derivatives
426  twz= ( thz(i,j,k)*w(i,j,k) - thz(i,j,k-1)*w(i,j,k-1) ) / dz(k)
427 
428  ! horizontal derivative, check to use 2nd or 4th schemes
429  if ( cds .eq. 1) then
430 
431  txx= ( t(i+1,j,k) - &
432  2.d0*t(i,j,k) + &
433  t(i-1,j,k) ) / dx2
434 
435  tyy= ( t(i,j+1,k) - &
436  2.d0*t(i,j,k) + &
437  t(i,j-1,k) ) / dy2
438 
439  elseif ( cds .eq. 2) then
440 
441  txx= ( -t(i+2,j,k) + &
442  16.d0*t(i+1,j,k) - &
443  30.d0*t(i,j,k) + &
444  16.d0*t(i-1,j,k) - &
445  t(i-2,j,k) ) / 12.d0 / dx2
446 
447  tyy= ( -t(i,j+2,k) + &
448  16.d0*t(i,j+1,k) - &
449  30.d0*t(i,j,k) + &
450  16.d0*t(i,j-1,k) - &
451  t(i,j-2,k) ) / 12.d0 / dy2
452 
453  endif
454 
455  ! we always use the 2nd order scheme for vertical derivatives
456  tzz= ( (t(i,j,k+1)-t(i,j,k))/dz_t(k) - &
457  (t(i,j,k)-t(i,j,k-1))/dz_b(k) ) / dz(k)
458 
459  ! nonlinear term for T
460  ss(i,j,k) = -tux-tvy-twz - u_mrf*tx
461 
462  ! the right hand side term for T
463  qt(i,j,k)= t(i,j,k) + dt*( rk_a*ss(i,j,k)+rk_b*ss0(i,j,k) ) + &
464  rk_c*0.5d0*dt * alpha * (txx+tyy+tzz)
465 
466 
467  enddo
468  enddo
469  enddo
470 
471  return
472 
473  end subroutine
474 
475 
476 
477 
478 
479 
480 
481 end module
real(mytype), save, protected rk_b
real(mytype), save, protected fc
real(mytype), save, protected dt
real(mytype), save, protected rk_a
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ug
real(mytype), dimension(:), allocatable, save, protected dz
real(mytype), save, protected dy2
integer, save, protected isscalar
real(mytype), save, protected alpha
integer, save, protected istr3
real(mytype), save, protected vg
real(mytype), save, protected t_ref
real(mytype), save, protected dx
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
real(mytype), save, protected dx2
subroutine get_temperature_cnvdiff(u, v, w, t, thx, thy, thz, ss, ss0, qt)
real(mytype), save, protected u_mrf
integer, save, protected cds
real(mytype), save, protected dpx_drive
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), save, protected dy
real(mytype), save, protected got
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), dimension(:), allocatable, save, protected l_t
subroutine get_momentum_cnvdiff(u, v, w, p, t, uhx, uhy, uhz, vhx, vhy, vhz, whx, why, ff, ff0, gg, gg0, hh, hh0, qu, qv, qw)
integer, save, protected dp_opt
integer, save, protected kstr3
integer, save, protected kend3