Calculates the Coriolis and momentum advection contributions to the acceleration.
112 type(ocean_grid_type),
intent(in) :: G
113 type(verticalGrid_type),
intent(in) :: GV
114 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u
115 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v
116 real,
dimension(SZI_(G),SZJ_(G),SZK_(G)),
intent(in) :: h
117 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: uh
119 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: vh
121 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(out) :: CAu
123 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(out) :: CAv
125 type(ocean_OBC_type),
pointer :: OBC
126 type(accel_diag_ptrs),
intent(inout) :: AD
127 type(unit_scale_type),
intent(in) :: US
128 type(CoriolisAdv_CS),
pointer :: CS
131 real,
dimension(SZIB_(G),SZJB_(G)) :: &
136 real,
dimension(SZIB_(G),SZJ_(G)) :: &
142 real,
dimension(SZI_(G),SZJ_(G)) :: &
146 real,
dimension(SZIB_(G),SZJ_(G)) :: &
152 real,
dimension(SZI_(G),SZJB_(G)) :: &
158 real,
dimension(SZI_(G),SZJ_(G)) :: &
164 real,
dimension(SZIB_(G),SZJB_(G)) :: &
172 real,
dimension(SZIB_(G),SZJB_(G),SZK_(G)) :: &
175 real :: fv1, fv2, fu1, fu2
176 real :: max_fv, max_fu
177 real :: min_fv, min_fu
179 real,
parameter :: C1_12=1.0/12.0
180 real,
parameter :: C1_24=1.0/24.0
181 real :: absolute_vorticity
182 real :: relative_vorticity
184 real :: max_Ihq, min_Ihq
190 real,
parameter :: eps_vel=1.0e-10
194 real :: c1, c2, c3, slope
210 real :: QUHeff,QVHeff
211 integer :: i, j, k, n, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
218 if (.not.
associated(cs))
call mom_error(fatal, &
219 "MOM_CoriolisAdv: Module must be initialized before it is used.")
220 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
221 isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB ; nz = g%ke
222 h_neglect = gv%H_subroundoff
223 h_tiny = gv%Angstrom_H
226 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+2
227 area_h(i,j) = g%mask2dT(i,j) * g%areaT(i,j)
229 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
230 if (.not. obc%segment(n)%on_pe) cycle
231 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
232 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1))
then
233 do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
234 if (obc%segment(n)%direction == obc_direction_n)
then
235 area_h(i,j+1) = area_h(i,j)
237 area_h(i,j) = area_h(i,j+1)
240 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1))
then
241 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
242 if (obc%segment(n)%direction == obc_direction_e)
then
243 area_h(i+1,j) = area_h(i,j)
245 area_h(i,j) = area_h(i+1,j)
251 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+1
252 area_q(i,j) = (area_h(i,j) + area_h(i+1,j+1)) + &
253 (area_h(i+1,j) + area_h(i,j+1))
264 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+1
265 dvdx(i,j) = v(i+1,j,k)*g%dyCv(i+1,j) - v(i,j,k)*g%dyCv(i,j)
266 dudy(i,j) = u(i,j+1,k)*g%dxCu(i,j+1) - u(i,j,k)*g%dxCu(i,j)
268 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+2
269 harea_v(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i,j+1) * h(i,j+1,k))
271 do j=jsq-1,jeq+2 ;
do i=isq-1,ieq+1
272 harea_u(i,j) = 0.5*(area_h(i,j) * h(i,j,k) + area_h(i+1,j) * h(i+1,j,k))
274 if (cs%Coriolis_En_Dis)
then
275 do j=jsq,jeq+1 ;
do i=is-1,ie
276 uh_center(i,j) = 0.5 * (g%dy_Cu(i,j) * u(i,j,k)) * (h(i,j,k) + h(i+1,j,k))
278 do j=js-1,je ;
do i=isq,ieq+1
279 vh_center(i,j) = 0.5 * (g%dx_Cv(i,j) * v(i,j,k)) * (h(i,j,k) + h(i,j+1,k))
285 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
286 if (.not. obc%segment(n)%on_pe) cycle
287 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
288 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1))
then
289 if (obc%zero_vorticity)
then ;
do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
290 dvdx(i,j) = 0. ; dudy(i,j) = 0.
292 if (obc%freeslip_vorticity)
then ;
do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
295 if (obc%computed_vorticity)
then ;
do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
296 if (obc%segment(n)%direction == obc_direction_n)
then
297 dudy(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - u(i,j,k))*g%dxCu(i,j)
299 dudy(i,j) = 2.0*(u(i,j+1,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dxCu(i,j+1)
302 if (obc%specified_vorticity)
then ;
do i=obc%segment(n)%HI%IsdB,obc%segment(n)%HI%IedB
303 if (obc%segment(n)%direction == obc_direction_n)
then
304 dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j)*g%dyBu(i,j)
306 dudy(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dxCu(i,j+1)*g%dyBu(i,j)
311 do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
312 if (obc%segment(n)%direction == obc_direction_n)
then
313 harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j,k)
315 harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
319 if (cs%Coriolis_En_Dis)
then
320 do i = max(isq-1,obc%segment(n)%HI%isd), min(ieq+2,obc%segment(n)%HI%ied)
321 if (obc%segment(n)%direction == obc_direction_n)
then
322 vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j,k)
324 vh_center(i,j) = g%dx_Cv(i,j) * v(i,j,k) * h(i,j+1,k)
328 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1))
then
329 if (obc%zero_vorticity)
then ;
do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
330 dvdx(i,j) = 0. ; dudy(i,j) = 0.
332 if (obc%freeslip_vorticity)
then ;
do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
335 if (obc%computed_vorticity)
then ;
do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
336 if (obc%segment(n)%direction == obc_direction_e)
then
337 dvdx(i,j) = 2.0*(obc%segment(n)%tangential_vel(i,j,k) - v(i,j,k))*g%dyCv(i,j)
339 dvdx(i,j) = 2.0*(v(i+1,j,k) - obc%segment(n)%tangential_vel(i,j,k))*g%dyCv(i+1,j)
342 if (obc%specified_vorticity)
then ;
do j=obc%segment(n)%HI%JsdB,obc%segment(n)%HI%JedB
343 if (obc%segment(n)%direction == obc_direction_e)
then
344 dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i,j)*g%dxBu(i,j)
346 dvdx(i,j) = obc%segment(n)%tangential_grad(i,j,k)*g%dyCv(i+1,j)*g%dxBu(i,j)
351 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
352 if (obc%segment(n)%direction == obc_direction_e)
then
353 harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i,j,k)
355 harea_u(i,j) = 0.5*(area_h(i,j) + area_h(i+1,j)) * h(i+1,j,k)
358 if (cs%Coriolis_En_Dis)
then
359 do j = max(jsq-1,obc%segment(n)%HI%jsd), min(jeq+2,obc%segment(n)%HI%jed)
360 if (obc%segment(n)%direction == obc_direction_e)
then
361 uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i,j,k)
363 uh_center(i,j) = g%dy_Cu(i,j) * u(i,j,k) * h(i+1,j,k)
370 if (
associated(obc))
then ;
do n=1,obc%number_of_segments
371 if (.not. obc%segment(n)%on_pe) cycle
374 i = obc%segment(n)%HI%IsdB ; j = obc%segment(n)%HI%JsdB
375 if (obc%segment(n)%is_N_or_S .and. (j >= jsq-1) .and. (j <= jeq+1))
then
376 do i = max(isq-1,obc%segment(n)%HI%IsdB), min(ieq+1,obc%segment(n)%HI%IedB)
377 if (obc%segment(n)%direction == obc_direction_n)
then
378 if (area_h(i,j) + area_h(i+1,j) > 0.0)
then
379 harea_u(i,j+1) = harea_u(i,j) * ((area_h(i,j+1) + area_h(i+1,j+1)) / &
380 (area_h(i,j) + area_h(i+1,j)))
381 else ; harea_u(i,j+1) = 0.0 ;
endif
383 if (area_h(i,j+1) + area_h(i+1,j+1) > 0.0)
then
384 harea_u(i,j) = harea_u(i,j+1) * ((area_h(i,j) + area_h(i+1,j)) / &
385 (area_h(i,j+1) + area_h(i+1,j+1)))
386 else ; harea_u(i,j) = 0.0 ;
endif
389 elseif (obc%segment(n)%is_E_or_W .and. (i >= isq-1) .and. (i <= ieq+1))
then
390 do j = max(jsq-1,obc%segment(n)%HI%JsdB), min(jeq+1,obc%segment(n)%HI%JedB)
391 if (obc%segment(n)%direction == obc_direction_e)
then
392 if (area_h(i,j) + area_h(i,j+1) > 0.0)
then
393 harea_v(i+1,j) = harea_v(i,j) * ((area_h(i+1,j) + area_h(i+1,j+1)) / &
394 (area_h(i,j) + area_h(i,j+1)))
395 else ; harea_v(i+1,j) = 0.0 ;
endif
397 harea_v(i,j) = 0.5 * (area_h(i,j) + area_h(i,j+1)) * h(i,j+1,k)
398 if (area_h(i+1,j) + area_h(i+1,j+1) > 0.0)
then
399 harea_v(i,j) = harea_v(i+1,j) * ((area_h(i,j) + area_h(i,j+1)) / &
400 (area_h(i+1,j) + area_h(i+1,j+1)))
401 else ; harea_v(i,j) = 0.0 ;
endif
407 do j=jsq-1,jeq+1 ;
do i=isq-1,ieq+1
408 if (cs%no_slip )
then
409 relative_vorticity = (2.0-g%mask2dBu(i,j)) * (dvdx(i,j) - dudy(i,j)) * &
412 relative_vorticity = g%mask2dBu(i,j) * (dvdx(i,j) - dudy(i,j)) * &
415 absolute_vorticity = us%s_to_T*g%CoriolisBu(i,j) + relative_vorticity
417 if (area_q(i,j) > 0.0)
then
418 harea_q = (harea_u(i,j) + harea_u(i,j+1)) + (harea_v(i,j) + harea_v(i+1,j))
419 ih = area_q(i,j) / (harea_q + h_neglect*area_q(i,j))
421 q(i,j) = absolute_vorticity * ih
422 abs_vort(i,j) = absolute_vorticity
425 if (cs%bound_Coriolis)
then
426 fv1 = absolute_vorticity*v(i+1,j,k)
427 fv2 = absolute_vorticity*v(i,j,k)
428 fu1 = -absolute_vorticity*u(i,j+1,k)
429 fu2 = -absolute_vorticity*u(i,j,k)
431 max_fvq(i,j) = fv1 ; min_fvq(i,j) = fv2
433 max_fvq(i,j) = fv2 ; min_fvq(i,j) = fv1
436 max_fuq(i,j) = fu1 ; min_fuq(i,j) = fu2
438 max_fuq(i,j) = fu2 ; min_fuq(i,j) = fu1
442 if (cs%id_rv > 0) rv(i,j,k) = relative_vorticity
443 if (cs%id_PV > 0) pv(i,j,k) = q(i,j)
444 if (
associated(ad%rv_x_v) .or.
associated(ad%rv_x_u)) &
445 q2(i,j) = relative_vorticity * ih
452 if (cs%Coriolis_Scheme == arakawa_hsu90)
then
455 a(i,j) = (q(i,j) + (q(i+1,j) + q(i,j-1))) * c1_12
456 d(i,j) = ((q(i,j) + q(i+1,j-1)) + q(i,j-1)) * c1_12
459 b(i,j) = (q(i,j) + (q(i-1,j) + q(i,j-1))) * c1_12
460 c(i,j) = ((q(i,j) + q(i-1,j-1)) + q(i,j-1)) * c1_12
463 elseif (cs%Coriolis_Scheme == arakawa_lamb81)
then
464 do j=jsq,jeq+1 ;
do i=isq,ieq+1
465 a(i-1,j) = (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
466 d(i-1,j) = ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
467 b(i,j) = ((q(i,j) + q(i-1,j-1)) + 2.0*(q(i-1,j) + q(i,j-1))) * c1_24
468 c(i,j) = (2.0*(q(i,j) + q(i-1,j-1)) + (q(i-1,j) + q(i,j-1))) * c1_24
469 ep_u(i,j) = ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
470 ep_v(i,j) = (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
472 elseif (cs%Coriolis_Scheme == al_blend)
then
473 fe_m2 = cs%F_eff_max_blend - 2.0
474 rat_lin = 1.5 * fe_m2 / max(cs%wt_lin_blend, 1.0e-16)
477 if (cs%F_eff_max_blend <= 2.0)
then ; fe_m2 = -1. ; rat_lin = -1.0 ;
endif
479 do j=jsq,jeq+1 ;
do i=isq,ieq+1
480 min_ihq = min(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
481 max_ihq = max(ih_q(i-1,j-1), ih_q(i,j-1), ih_q(i-1,j), ih_q(i,j))
483 if (max_ihq < 1.0e15*min_ihq) rat_m1 = max_ihq / min_ihq - 1.0
490 if (rat_m1 <= fe_m2)
then ; al_wt = 1.0
491 elseif (rat_m1 < 1.5*fe_m2)
then ; al_wt = 3.0*fe_m2 / rat_m1 - 2.0
492 else ; al_wt = 0.0 ;
endif
495 if (rat_m1 <= 1.5*fe_m2)
then ; sad_wt = 0.0
496 elseif (rat_m1 <= rat_lin)
then
497 sad_wt = 1.0 - (1.5*fe_m2) / rat_m1
498 elseif (rat_m1 < 2.0*rat_lin)
then
499 sad_wt = 1.0 - (cs%wt_lin_blend / rat_lin) * (rat_m1 - 2.0*rat_lin)
500 else ; sad_wt = 1.0 ;
endif
502 a(i-1,j) = sad_wt * 0.25 * q(i-1,j) + (1.0 - sad_wt) * &
503 ( ((2.0-al_wt)* q(i-1,j) + al_wt*q(i,j-1)) + &
504 2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
505 d(i-1,j) = sad_wt * 0.25 * q(i-1,j-1) + (1.0 - sad_wt) * &
506 ( ((2.0-al_wt)* q(i-1,j-1) + al_wt*q(i,j)) + &
507 2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
508 b(i,j) = sad_wt * 0.25 * q(i,j) + (1.0 - sad_wt) * &
509 ( ((2.0-al_wt)* q(i,j) + al_wt*q(i-1,j-1)) + &
510 2.0 * (q(i-1,j) + q(i,j-1)) ) * c1_24
511 c(i,j) = sad_wt * 0.25 * q(i,j-1) + (1.0 - sad_wt) * &
512 ( ((2.0-al_wt)* q(i,j-1) + al_wt*q(i-1,j)) + &
513 2.0 * (q(i,j) + q(i-1,j-1)) ) * c1_24
514 ep_u(i,j) = al_wt * ((q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
515 ep_v(i,j) = al_wt * (-(q(i,j) - q(i-1,j-1)) + (q(i-1,j) - q(i,j-1))) * c1_24
519 if (cs%Coriolis_En_Dis)
then
521 c1 = 1.0-1.5*0.5 ; c2 = 1.0-0.5 ; c3 = 2.0 ; slope = 0.5
523 do j=jsq,jeq+1 ;
do i=is-1,ie
527 if (g%dy_Cu(i,j) == 0.0) uhc = uhm
529 if (abs(uhc) < 0.1*abs(uhm))
then
531 elseif (abs(uhc) > c1*abs(uhm))
then
532 if (abs(uhc) < c2*abs(uhm))
then ; uhc = (3.0*uhc+(1.0-c2*3.0)*uhm)
533 elseif (abs(uhc) <= c3*abs(uhm))
then ; uhc = uhm
534 else ; uhc = slope*uhc+(1.0-c3*slope)*uhm
539 uh_min(i,j) = uhm ; uh_max(i,j) = uhc
541 uh_max(i,j) = uhm ; uh_min(i,j) = uhc
544 do j=js-1,je ;
do i=isq,ieq+1
548 if (g%dx_Cv(i,j) == 0.0) vhc = vhm
550 if (abs(vhc) < 0.1*abs(vhm))
then
552 elseif (abs(vhc) > c1*abs(vhm))
then
553 if (abs(vhc) < c2*abs(vhm))
then ; vhc = (3.0*vhc+(1.0-c2*3.0)*vhm)
554 elseif (abs(vhc) <= c3*abs(vhm))
then ; vhc = vhm
555 else ; vhc = slope*vhc+(1.0-c3*slope)*vhm
560 vh_min(i,j) = vhm ; vh_max(i,j) = vhc
562 vh_max(i,j) = vhm ; vh_min(i,j) = vhc
568 call gradke(u, v, h, ke, kex, key, k, obc, g, cs)
573 if (cs%Coriolis_Scheme == sadourny75_energy)
then
574 if (cs%Coriolis_En_Dis)
then
576 do j=js,je ;
do i=isq,ieq
577 if (q(i,j)*u(i,j,k) == 0.0)
then
578 temp1 = q(i,j) * ( (vh_max(i,j)+vh_max(i+1,j)) &
579 + (vh_min(i,j)+vh_min(i+1,j)) )*0.5
580 elseif (q(i,j)*u(i,j,k) < 0.0)
then
581 temp1 = q(i,j) * (vh_max(i,j)+vh_max(i+1,j))
583 temp1 = q(i,j) * (vh_min(i,j)+vh_min(i+1,j))
585 if (q(i,j-1)*u(i,j,k) == 0.0)
then
586 temp2 = q(i,j-1) * ( (vh_max(i,j-1)+vh_max(i+1,j-1)) &
587 + (vh_min(i,j-1)+vh_min(i+1,j-1)) )*0.5
588 elseif (q(i,j-1)*u(i,j,k) < 0.0)
then
589 temp2 = q(i,j-1) * (vh_max(i,j-1)+vh_max(i+1,j-1))
591 temp2 = q(i,j-1) * (vh_min(i,j-1)+vh_min(i+1,j-1))
593 cau(i,j,k) = 0.25 * g%IdxCu(i,j) * (temp1 + temp2)
597 do j=js,je ;
do i=isq,ieq
598 cau(i,j,k) = 0.25 * &
599 (q(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
600 q(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
603 elseif (cs%Coriolis_Scheme == sadourny75_enstro)
then
604 do j=js,je ;
do i=isq,ieq
605 cau(i,j,k) = 0.125 * (g%IdxCu(i,j) * (q(i,j) + q(i,j-1))) * &
606 ((vh(i+1,j,k) + vh(i,j,k)) + (vh(i,j-1,k) + vh(i+1,j-1,k)))
608 elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
609 (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
610 (cs%Coriolis_Scheme == al_blend))
then
612 do j=js,je ;
do i=isq,ieq
613 cau(i,j,k) = ((a(i,j) * vh(i+1,j,k) + &
614 c(i,j) * vh(i,j-1,k)) &
615 + (b(i,j) * vh(i,j,k) + &
616 d(i,j) * vh(i+1,j-1,k))) * g%IdxCu(i,j)
618 elseif (cs%Coriolis_Scheme == robust_enstro)
then
622 do j=js,je ;
do i=isq,ieq
623 heff1 = abs(vh(i,j,k)*g%IdxCv(i,j))/(eps_vel+abs(v(i,j,k)))
624 heff1 = max(heff1,min(h(i,j,k),h(i,j+1,k)))
625 heff1 = min(heff1,max(h(i,j,k),h(i,j+1,k)))
626 heff2 = abs(vh(i,j-1,k)*g%IdxCv(i,j-1))/(eps_vel+abs(v(i,j-1,k)))
627 heff2 = max(heff2,min(h(i,j-1,k),h(i,j,k)))
628 heff2 = min(heff2,max(h(i,j-1,k),h(i,j,k)))
629 heff3 = abs(vh(i+1,j,k)*g%IdxCv(i+1,j))/(eps_vel+abs(v(i+1,j,k)))
630 heff3 = max(heff3,min(h(i+1,j,k),h(i+1,j+1,k)))
631 heff3 = min(heff3,max(h(i+1,j,k),h(i+1,j+1,k)))
632 heff4 = abs(vh(i+1,j-1,k)*g%IdxCv(i+1,j-1))/(eps_vel+abs(v(i+1,j-1,k)))
633 heff4 = max(heff4,min(h(i+1,j-1,k),h(i+1,j,k)))
634 heff4 = min(heff4,max(h(i+1,j-1,k),h(i+1,j,k)))
635 if (cs%PV_Adv_Scheme == pv_adv_centered)
then
636 cau(i,j,k) = 0.5*(abs_vort(i,j)+abs_vort(i,j-1)) * &
637 ((vh(i ,j ,k)+vh(i+1,j-1,k)) + &
638 (vh(i ,j-1,k)+vh(i+1,j ,k)) ) / &
639 (h_tiny +((heff1+heff4) +(heff2+heff3)) ) * g%IdxCu(i,j)
640 elseif (cs%PV_Adv_Scheme == pv_adv_upwind1)
then
641 vheff = ((vh(i ,j ,k)+vh(i+1,j-1,k)) + &
642 (vh(i ,j-1,k)+vh(i+1,j ,k)) )
643 qvheff = 0.5*( (abs_vort(i,j)+abs_vort(i,j-1))*vheff &
644 -(abs_vort(i,j)-abs_vort(i,j-1))*abs(vheff) )
645 cau(i,j,k) = qvheff / &
646 (h_tiny +((heff1+heff4) +(heff2+heff3)) ) * g%IdxCu(i,j)
651 if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
652 (cs%Coriolis_Scheme == al_blend))
then ;
do j=js,je ;
do i=isq,ieq
653 cau(i,j,k) = cau(i,j,k) + &
654 (ep_u(i,j)*uh(i-1,j,k) - ep_u(i+1,j)*uh(i+1,j,k)) * g%IdxCu(i,j)
655 enddo ;
enddo ;
endif
658 if (cs%bound_Coriolis)
then
659 do j=js,je ;
do i=isq,ieq
660 max_fv = max(max_fvq(i,j),max_fvq(i,j-1))
661 min_fv = min(min_fvq(i,j),min_fvq(i,j-1))
664 if (cau(i,j,k) > max_fv)
then
667 if (cau(i,j,k) < min_fv) cau(i,j,k) = min_fv
673 do j=js,je ;
do i=isq,ieq
674 cau(i,j,k) = cau(i,j,k) - kex(i,j)
675 if (
associated(ad%gradKEu)) ad%gradKEu(i,j,k) = -kex(i,j)
682 if (cs%Coriolis_Scheme == sadourny75_energy)
then
683 if (cs%Coriolis_En_Dis)
then
685 do j=jsq,jeq ;
do i=is,ie
686 if (q(i-1,j)*v(i,j,k) == 0.0)
then
687 temp1 = q(i-1,j) * ( (uh_max(i-1,j)+uh_max(i-1,j+1)) &
688 + (uh_min(i-1,j)+uh_min(i-1,j+1)) )*0.5
689 elseif (q(i-1,j)*v(i,j,k) > 0.0)
then
690 temp1 = q(i-1,j) * (uh_max(i-1,j)+uh_max(i-1,j+1))
692 temp1 = q(i-1,j) * (uh_min(i-1,j)+uh_min(i-1,j+1))
694 if (q(i,j)*v(i,j,k) == 0.0)
then
695 temp2 = q(i,j) * ( (uh_max(i,j)+uh_max(i,j+1)) &
696 + (uh_min(i,j)+uh_min(i,j+1)) )*0.5
697 elseif (q(i,j)*v(i,j,k) > 0.0)
then
698 temp2 = q(i,j) * (uh_max(i,j)+uh_max(i,j+1))
700 temp2 = q(i,j) * (uh_min(i,j)+uh_min(i,j+1))
702 cav(i,j,k) = - 0.25 * g%IdyCv(i,j) * (temp1 + temp2)
706 do j=jsq,jeq ;
do i=is,ie
707 cav(i,j,k) = - 0.25* &
708 (q(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
709 q(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
712 elseif (cs%Coriolis_Scheme == sadourny75_enstro)
then
713 do j=jsq,jeq ;
do i=is,ie
714 cav(i,j,k) = -0.125 * (g%IdyCv(i,j) * (q(i-1,j) + q(i,j))) * &
715 ((uh(i-1,j,k) + uh(i-1,j+1,k)) + (uh(i,j,k) + uh(i,j+1,k)))
717 elseif ((cs%Coriolis_Scheme == arakawa_hsu90) .or. &
718 (cs%Coriolis_Scheme == arakawa_lamb81) .or. &
719 (cs%Coriolis_Scheme == al_blend))
then
721 do j=jsq,jeq ;
do i=is,ie
722 cav(i,j,k) = - ((a(i-1,j) * uh(i-1,j,k) + &
723 c(i,j+1) * uh(i,j+1,k)) &
724 + (b(i,j) * uh(i,j,k) + &
725 d(i-1,j+1) * uh(i-1,j+1,k))) * g%IdyCv(i,j)
727 elseif (cs%Coriolis_Scheme == robust_enstro)
then
731 do j=jsq,jeq ;
do i=is,ie
732 heff1 = abs(uh(i,j,k)*g%IdyCu(i,j))/(eps_vel+abs(u(i,j,k)))
733 heff1 = max(heff1,min(h(i,j,k),h(i+1,j,k)))
734 heff1 = min(heff1,max(h(i,j,k),h(i+1,j,k)))
735 heff2 = abs(uh(i-1,j,k)*g%IdyCu(i-1,j))/(eps_vel+abs(u(i-1,j,k)))
736 heff2 = max(heff2,min(h(i-1,j,k),h(i,j,k)))
737 heff2 = min(heff2,max(h(i-1,j,k),h(i,j,k)))
738 heff3 = abs(uh(i,j+1,k)*g%IdyCu(i,j+1))/(eps_vel+abs(u(i,j+1,k)))
739 heff3 = max(heff3,min(h(i,j+1,k),h(i+1,j+1,k)))
740 heff3 = min(heff3,max(h(i,j+1,k),h(i+1,j+1,k)))
741 heff4 = abs(uh(i-1,j+1,k)*g%IdyCu(i-1,j+1))/(eps_vel+abs(u(i-1,j+1,k)))
742 heff4 = max(heff4,min(h(i-1,j+1,k),h(i,j+1,k)))
743 heff4 = min(heff4,max(h(i-1,j+1,k),h(i,j+1,k)))
744 if (cs%PV_Adv_Scheme == pv_adv_centered)
then
745 cav(i,j,k) = - 0.5*(abs_vort(i,j)+abs_vort(i-1,j)) * &
746 ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
747 (uh(i-1,j ,k)+uh(i ,j+1,k)) ) / &
748 (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
749 elseif (cs%PV_Adv_Scheme == pv_adv_upwind1)
then
750 uheff = ((uh(i ,j ,k)+uh(i-1,j+1,k)) + &
751 (uh(i-1,j ,k)+uh(i ,j+1,k)) )
752 quheff = 0.5*( (abs_vort(i,j)+abs_vort(i-1,j))*uheff &
753 -(abs_vort(i,j)-abs_vort(i-1,j))*abs(uheff) )
754 cav(i,j,k) = - quheff / &
755 (h_tiny + ((heff1+heff4) +(heff2+heff3)) ) * g%IdyCv(i,j)
760 if ((cs%Coriolis_Scheme == arakawa_lamb81) .or. &
761 (cs%Coriolis_Scheme == al_blend))
then ;
do j=jsq,jeq ;
do i=is,ie
762 cav(i,j,k) = cav(i,j,k) + &
763 (ep_v(i,j)*vh(i,j-1,k) - ep_v(i,j+1)*vh(i,j+1,k)) * g%IdyCv(i,j)
764 enddo ;
enddo ;
endif
766 if (cs%bound_Coriolis)
then
767 do j=jsq,jeq ;
do i=is,ie
768 max_fu = max(max_fuq(i,j),max_fuq(i-1,j))
769 min_fu = min(min_fuq(i,j),min_fuq(i-1,j))
770 if (cav(i,j,k) > max_fu)
then
773 if (cav(i,j,k) < min_fu) cav(i,j,k) = min_fu
779 do j=jsq,jeq ;
do i=is,ie
780 cav(i,j,k) = cav(i,j,k) - key(i,j)
781 if (
associated(ad%gradKEv)) ad%gradKEv(i,j,k) = -key(i,j)
784 if (
associated(ad%rv_x_u) .or.
associated(ad%rv_x_v))
then
786 if (cs%Coriolis_Scheme == sadourny75_energy)
then
787 if (
associated(ad%rv_x_u))
then
788 do j=jsq,jeq ;
do i=is,ie
789 ad%rv_x_u(i,j,k) = - 0.25* &
790 (q2(i-1,j)*(uh(i-1,j,k) + uh(i-1,j+1,k)) + &
791 q2(i,j)*(uh(i,j,k) + uh(i,j+1,k))) * g%IdyCv(i,j)
795 if (
associated(ad%rv_x_v))
then
796 do j=js,je ;
do i=isq,ieq
797 ad%rv_x_v(i,j,k) = 0.25 * &
798 (q2(i,j) * (vh(i+1,j,k) + vh(i,j,k)) + &
799 q2(i,j-1) * (vh(i,j-1,k) + vh(i+1,j-1,k))) * g%IdxCu(i,j)
803 if (
associated(ad%rv_x_u))
then
804 do j=jsq,jeq ;
do i=is,ie
805 ad%rv_x_u(i,j,k) = -g%IdyCv(i,j) * c1_12 * &
806 ((q2(i,j) + q2(i-1,j) + q2(i-1,j-1)) * uh(i-1,j,k) + &
807 (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * uh(i,j,k) + &
808 (q2(i-1,j) + q2(i,j+1) + q2(i,j)) * uh(i,j+1,k) + &
809 (q2(i,j) + q2(i-1,j+1) + q2(i-1,j)) * uh(i-1,j+1,k))
813 if (
associated(ad%rv_x_v))
then
814 do j=js,je ;
do i=isq,ieq
815 ad%rv_x_v(i,j,k) = g%IdxCu(i,j) * c1_12 * &
816 ((q2(i+1,j) + q2(i,j) + q2(i,j-1)) * vh(i+1,j,k) + &
817 (q2(i-1,j) + q2(i,j) + q2(i,j-1)) * vh(i,j,k) + &
818 (q2(i-1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i,j-1,k) + &
819 (q2(i+1,j-1) + q2(i,j) + q2(i,j-1)) * vh(i+1,j-1,k))
828 if (query_averaging_enabled(cs%diag))
then
829 if (cs%id_rv > 0)
call post_data(cs%id_rv, rv, cs%diag)
830 if (cs%id_PV > 0)
call post_data(cs%id_PV, pv, cs%diag)
831 if (cs%id_gKEu>0)
call post_data(cs%id_gKEu, ad%gradKEu, cs%diag)
832 if (cs%id_gKEv>0)
call post_data(cs%id_gKEv, ad%gradKEv, cs%diag)
833 if (cs%id_rvxu > 0)
call post_data(cs%id_rvxu, ad%rv_x_u, cs%diag)
834 if (cs%id_rvxv > 0)
call post_data(cs%id_rvxv, ad%rv_x_v, cs%diag)