Changeset 15727
- Timestamp:
- 2022-02-23T16:54:50+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0-TRUNK_r14960_HPG/src/OCE/DYN/dynhpg.F90
r15726 r15727 1 1 MODULE dynhpg 2 3 !! v2: takes into account the stretching of the vertical grid in hpg_ffr 4 2 5 !!====================================================================== 3 6 !! *** MODULE dynhpg *** … … 1738 1741 1739 1742 IF ( ln_dbg_hpg ) CALL dbg_2dr( '2.1 ht_0', ht_0 ) 1740 CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref) ! Uses rhd (IN) to calculate all other fields (OUT)1743 CALL calc_rhd_ref(j_uv, jn_hor_pts, gde3w, zrhd_ref, zz_ref, jk_bot_ref) ! Uses rhd (IN) to calculate all other fields (OUT) 1741 1744 1742 1745 IF ( ln_dbg_hpg ) THEN … … 2056 2059 2057 2060 REAL(wp) :: zz_tgt_lcl, zfld_ref_on_tgt 2061 LOGICAL :: ll_ccs 2058 2062 2059 2063 !----------------------------------------------------------------------------------------------------------------- 2060 2064 2065 ll_ccs = .FALSE. ! set for the call to loc_ref_tgt 2061 2066 z_r_6 = 1._wp / 6._wp 2062 2067 z_r_24 = 1._wp / 24._wp 2063 2068 2064 2069 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 2065 CALL loc_ref_tgt ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt )2070 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, jk_ref_for_tgt ) 2066 2071 2067 2072 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) … … 2099 2104 2100 2105 END_3D 2101 2102 RETURN2103 2106 2104 2107 END SUBROUTINE ref_to_tgt_cub … … 2128 2131 REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3 2129 2132 REAL(wp) :: zfld_ref_on_tgt 2133 LOGICAL :: ll_ccs ! set to .FALSE. for the call to loc_ref_tgt 2134 2130 2135 !----------------------------------------------------------------------------------------------------------------- 2131 2136 2137 ll_ccs = .TRUE. ! set for the call to loc_ref_tgt 2132 2138 z_r_6 = 1._wp / 6._wp 2133 2139 2134 2140 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 2135 CALL loc_ref_tgt ( ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt )2141 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt ) 2136 2142 2137 2143 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) … … 2172 2178 IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs: pfld_tgt_ref', pfld_tgt_ref ) 2173 2179 2174 RETURN2175 2176 2180 END SUBROUTINE ref_to_tgt_ccs 2177 2181 2178 2182 !----------------------------------------------------------------------------------------------------------------- 2179 2183 2180 SUBROUTINE loc_ref_tgt ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt ) 2181 2184 SUBROUTINE loc_ref_tgt (ll_ccs, ki_off_tgt, kj_off_tgt, p_dep_tgt, p_z_ref, kk_bot_ref, kk_ref_for_tgt ) 2185 2186 LOGICAL, INTENT(in) :: ll_ccs ! true => constrained cubic spline; false => pure cubic fit 2182 2187 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2183 2188 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction … … 2206 2211 ! 1. Initialisation for search for points on reference grid bounding points on the target grid 2207 2212 ! the first point on the target grid is assumed to lie between the first two points on the reference grid 2208 IF ( l n_hpg_djr_ref_ccs ) THEN2213 IF ( ll_ccs ) THEN 2209 2214 jk_init = 2 2210 2215 ELSE … … 2288 2293 IF ( ( z_above < z_tgt .AND. jkr > jk_init ) .OR. z_below > (z_tgt + tol_wp) ) THEN 2289 2294 IF ( jcount < 10 ) THEN 2290 WRITE(numout,*) 'loc_ref_tgt: ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above '2291 WRITE(numout,*) ji, jj, jk, ki_off_tgt, kj_off_tgt, z_tgt, z_below, z_above2295 WRITE(numout,*) 'loc_ref_tgt: ji, jj, jk, ki_off_tgt, kj_off_tgt, jkr, z_tgt, z_below, z_above ' 2296 WRITE(numout,*) ji, jj, jk, ki_off_tgt, kj_off_tgt, jkr, z_tgt, z_below, z_above 2292 2297 END IF 2293 2298 jcount = jcount + 1 … … 2317 2322 ! This assumes that every sea point has at least 4 levels. 2318 2323 2319 IF ( .NOT. l n_hpg_djr_ref_ccs ) THEN2324 IF ( .NOT. ll_ccs ) THEN 2320 2325 DO_3D(0, 0, 0, 0, 2, jpk-1) 2321 2326 IF ( kk_ref_for_tgt(ji,jj, jk) == kk_bot_ref(ji,jj) ) kk_ref_for_tgt(ji,jj, jk) = kk_ref_for_tgt(ji,jj, jk) - 1 2322 2327 END_3D 2323 2328 END IF 2324 2325 RETURN2326 2329 2327 2330 END SUBROUTINE loc_ref_tgt … … 2350 2353 2351 2354 ! The following fields could probably be made single level or at most 2 level fields 2352 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_e3t_pmr, z_depw_pmr ! corresponding e3t and gdepw pmr profiles 2353 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rhd_mid, z_e3t_mid ! profiles and e3t interpolated to middle of cell (using ccs) 2355 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_dept_pmr, z_depw_pmr ! corresponding gdept and gdepw pmr profiles 2356 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rhd_mid ! rhd profiles interpolated to middle of cell (using ccs or cub) 2357 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_dept_mid, z_depw_mid ! dept and depw similarly interpolated to middle of cell (using ccs or cub) 2358 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_e3t_mid ! e3t profiles interpolated to middle of cell; only printed as a diagnostic 2354 2359 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_ddepw_ij ! constrained spline horizontal differences in gdepw 2355 2360 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_p_pmr, z_F_pmr ! pressures and forces on vertical faces … … 2365 2370 INTEGER :: ji_ro, jj_ro, jr_offset ! offsets for the reference profile 2366 2371 INTEGER :: jn_hor_pts ! number of profiles required in horizontal (4 if cubic interpolarion or 2 if not) 2367 2372 INTEGER :: j_t_levs, j_w_levs ! indicators that field passed in is valid on t-levels or w-levels 2373 2368 2374 REAL(wp) :: z_r_6 ! 1/6 2369 2375 REAL(wp) :: zr_a, zr_b, zr_c ! rhd evaluated at i, i+1/2 and i+1 … … 2379 2385 2380 2386 z_r_6 = 1.0_wp/6.0_wp 2387 2388 j_t_levs = 1 ! this is not very neat - it is duplicated in the calling routines 2389 j_w_levs = 2 2381 2390 2382 2391 IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN … … 2400 2409 END IF 2401 2410 2402 ! 1. Calculate the referen e profile from all points in the stencil2411 ! 1. Calculate the reference profile from all points in the stencil 2403 2412 2404 2413 IF ( ln_hpg_ffr_ref ) THEN 2405 CALL calc_rhd_ref(j_uv, jn_hor_pts, zrhd_ref, zz_ref, jk_bot_ref) ! Uses rhd (IN) to calculate all other fields (OUT)2414 CALL calc_rhd_ref(j_uv, jn_hor_pts, gdept(:,:,:,Kmm), zrhd_ref, zz_ref, jk_bot_ref) 2406 2415 IF ( ln_hpg_ffr_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 2407 2416 END IF … … 2425 2434 END IF 2426 2435 2427 IF ( ln_hpg_ffr_ref ) THEN !!!! MJB should gde3w below now be gdept ?? !!!!2428 2436 IF ( ln_hpg_ffr_ref ) THEN 2437 2429 2438 IF ( ln_hpg_ffr_ref_ccs ) THEN 2430 CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gde 3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )2439 CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2431 2440 ELSE 2432 CALL ref_to_tgt_cub ( ji_ro, jj_ro, gde 3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) )2441 CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2433 2442 END IF 2434 2443 2435 2444 ELSE ! no subtraction of reference profile 2436 2445 2437 DO_3D (0,0,0,0,1,jpk -1)2446 DO_3D (0,0,0,0,1,jpk) 2438 2447 z_rhd_pmr(ji,jj,jk,jr) = rhd(ji+ji_ro,jj+jj_ro,jk) 2439 2448 END_3D … … 2442 2451 2443 2452 DO_3D (0,0,0,0,1,jpk) 2444 z_ e3t_pmr (ji,jj,jk,jr) = e3t(ji+ji_ro,jj+jj_ro,jk,Kmm)2453 z_dept_pmr(ji,jj,jk,jr) = gdept(ji+ji_ro,jj+jj_ro,jk,Kmm) 2445 2454 z_depw_pmr(ji,jj,jk,jr) = gdepw(ji+ji_ro,jj+jj_ro,jk,Kmm) 2446 2455 END_3D … … 2458 2467 2459 2468 IF ( ln_hpg_ffr_hor_ccs ) THEN 2460 CALL calc_mid_ccs(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid) 2461 CALL calc_mid_ccs(j_uv, aco_bc_z_hor, bco_bc_z_hor, z_e3t_pmr, z_e3t_mid) 2469 CALL calc_mid_ccs(j_uv, j_t_levs, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid) 2470 CALL calc_mid_ccs(j_uv, j_t_levs, aco_bc_z_hor, bco_bc_z_hor, z_dept_pmr, z_dept_mid) 2471 CALL calc_mid_ccs(j_uv, j_w_levs, aco_bc_z_hor, bco_bc_z_hor, z_depw_pmr, z_depw_mid) 2462 2472 CALL calc_dz_dij_ccs( j_uv, z_depw_pmr, z_ddepw_ij ) 2463 2473 ELSE ! ln_hpg_ffr_hor_cub 2464 CALL calc_mid_cub(j_uv, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid) 2465 CALL calc_mid_cub(j_uv, aco_bc_z_hor, bco_bc_z_hor, z_e3t_pmr, z_e3t_mid) 2474 CALL calc_mid_cub(j_uv, j_t_levs, aco_bc_rhd_hor, bco_bc_rhd_hor, z_rhd_pmr, z_rhd_mid) 2475 CALL calc_mid_cub(j_uv, j_t_levs, aco_bc_z_hor, bco_bc_z_hor, z_dept_pmr, z_dept_mid) 2476 CALL calc_mid_cub(j_uv, j_w_levs, aco_bc_z_hor, bco_bc_z_hor, z_depw_pmr, z_depw_mid) 2466 2477 CALL calc_dz_dij_cub( j_uv, z_depw_pmr, z_ddepw_ij ) 2467 2478 END IF 2468 2479 2469 DO_3D (0,0,0,0,1,jpk -1)2480 DO_3D (0,0,0,0,1,jpk) 2470 2481 z_rhd_pmr(ji,jj,jk,1) = z_rhd_pmr(ji,jj,jk,2) 2471 2482 z_rhd_pmr(ji,jj,jk,2) = z_rhd_mid(ji,jj,jk) 2472 z_e3t_pmr(ji,jj,jk,1) = z_e3t_pmr(ji,jj,jk,2) 2473 z_e3t_pmr(ji,jj,jk,2) = z_e3t_mid(ji,jj,jk) 2483 z_dept_pmr(ji,jj,jk,1) = z_dept_pmr(ji,jj,jk,2) 2484 z_dept_pmr(ji,jj,jk,2) = z_dept_mid(ji,jj,jk) 2485 z_depw_pmr(ji,jj,jk,1) = z_depw_pmr(ji,jj,jk,2) 2486 z_depw_pmr(ji,jj,jk,2) = z_depw_mid(ji,jj,jk) 2474 2487 END_3D 2475 2488 2476 2489 ELSE ! simple linear interpolation 2477 2490 2478 DO_3D (0,0,0,0,1,jpk-1) 2479 z_rhd_pmr(ji,jj,jk,3) = z_rhd_pmr(ji,jj,jk,2) 2480 z_rhd_pmr(ji,jj,jk,2) = 0.5_wp*( z_rhd_pmr(ji,jj,jk,1) + z_rhd_pmr(ji,jj,jk,2) ) 2481 END_3D 2482 DO_3D (0,0,0,0,1,jpk-1) 2483 z_e3t_pmr(ji,jj,jk,1) = e3t(ji, jj, jk, Kmm) 2484 z_e3t_pmr(ji,jj,jk,2) = 0.5_wp*( e3t(ji, jj, jk, Kmm) + e3t(ji+jio, jj+jjo, jk, Kmm) ) 2485 z_e3t_pmr(ji,jj,jk,3) = e3t(ji+jio, jj+jjo, jk, Kmm) 2491 DO_3D (0,0,0,0,1,jpk) 2492 z_rhd_pmr(ji,jj,jk,3) = z_rhd_pmr(ji,jj,jk,2) 2493 z_rhd_pmr(ji,jj,jk,2) = 0.5_wp*( z_rhd_pmr(ji,jj,jk,1) + z_rhd_pmr(ji,jj,jk,2) ) 2494 z_dept_pmr(ji,jj,jk,3) = z_dept_pmr(ji,jj,jk,2) 2495 z_dept_pmr(ji,jj,jk,2) = 0.5_wp*( z_dept_pmr(ji,jj,jk,1) + z_dept_pmr(ji,jj,jk,2) ) 2496 z_depw_pmr(ji,jj,jk,3) = z_depw_pmr(ji,jj,jk,2) 2497 z_depw_pmr(ji,jj,jk,2) = 0.5_wp*( z_depw_pmr(ji,jj,jk,1) + z_depw_pmr(ji,jj,jk,2) ) 2486 2498 END_3D 2487 2499 END IF … … 2494 2506 2495 2507 IF ( ln_dbg_hpg ) THEN 2496 CALL dbg_3dr( '3. e3t ', e3t )2497 2508 CALL dbg_3dr( '3. z_rhd_pmr: 1', z_rhd_pmr(:,:,:,1) ) 2498 2509 CALL dbg_3dr( '3. z_rhd_pmr: 2', z_rhd_pmr(:,:,:,2) ) 2499 2510 CALL dbg_3dr( '3. z_rhd_pmr: 3', z_rhd_pmr(:,:,:,3) ) 2511 CALL dbg_3dr( '3. z_depw_pmr: 1', z_depw_pmr(:,:,:,1) ) 2512 CALL dbg_3dr( '3. z_depw_pmr: 2', z_depw_pmr(:,:,:,2) ) 2513 CALL dbg_3dr( '3. z_depw_pmr: 3', z_depw_pmr(:,:,:,3) ) 2514 DO jr = 1, 3 2515 DO_3D (0,0,0,0,1,jpk-1) 2516 z_e3t_mid(ji,jj,jk) = z_depw_pmr(ji,jj,jk,jr) - z_depw_pmr(ji,jj,jk+1,jr) 2517 END_3D 2518 CALL dbg_3dr( '3. z_e3t_mid jr : ', z_e3t_mid(:,:,:) ) 2519 END DO 2500 2520 END IF 2501 2521 … … 2504 2524 IF ( ln_hpg_ffr_vrt_quad ) THEN 2505 2525 DO jr = 1, 3 2506 CALL vrt_int_quad(j_mbkt(:,:,jr), z_rhd_pmr(:,:,:,jr), z_ e3t_pmr(:,:,:,jr), z_p_pmr(:,:,:,jr), z_F_pmr(:,:,:,jr))2526 CALL vrt_int_quad(j_mbkt(:,:,jr), z_rhd_pmr(:,:,:,jr), z_dept_pmr(:,:,:,jr), z_depw_pmr(:,:,:,jr), z_p_pmr(:,:,:,jr), z_F_pmr(:,:,:,jr)) 2507 2527 END DO ! jr 2508 2528 … … 2515 2535 DO jk = 1, jpk - 1 2516 2536 DO_2D(0,0,0,0) 2517 z_p_pmr(ji,jj,jk+1,jr) = z_p_pmr(ji,jj,jk,jr) + grav* z_e3t_pmr(ji,jj,jk,jr)*z_rhd_pmr(ji,jj,jk,jr)2537 z_p_pmr(ji,jj,jk+1,jr) = z_p_pmr(ji,jj,jk,jr) + grav*(z_depw_pmr(ji,jj,jk+1,jr) - z_depw_pmr(ji,jj,jk,jr)) *z_rhd_pmr(ji,jj,jk,jr) 2518 2538 END_2D 2519 2539 DO_2D(0,0,0,0) 2520 z_F_pmr(ji,jj,jk,jr) = 0.5_wp* z_e3t_pmr(ji,jj,jk,jr)*( z_p_pmr(ji,jj,jk,jr) + z_p_pmr(ji,jj,jk+1,jr) )2540 z_F_pmr(ji,jj,jk,jr) = 0.5_wp*(z_depw_pmr(ji,jj,jk+1,jr) - z_depw_pmr(ji,jj,jk,jr))*( z_p_pmr(ji,jj,jk,jr) + z_p_pmr(ji,jj,jk+1,jr) ) 2521 2541 END_2D 2522 2542 END DO ! jk … … 2531 2551 CALL dbg_3dr( '4. z_F_pmr: 1', z_F_pmr(:,:,:,1) ) 2532 2552 CALL dbg_3dr( '4. z_F_pmr: 3', z_F_pmr(:,:,:,3) ) 2533 CALL dbg_3dr( ' z_e3t_pmr: 1', z_e3t_pmr(:,:,:,1) )2534 CALL dbg_3dr( ' z_e3t_pmr: 2', z_e3t_pmr(:,:,:,2) )2535 CALL dbg_3dr( ' z_e3t_pmr: 3', z_e3t_pmr(:,:,:,3) )2536 2553 END IF 2537 2554 … … 2587 2604 CALL iom_put( "u_force_west", z_F_pmr(:,:,:,1) ) 2588 2605 CALL iom_put( "u_force_upper", z_F_upp ) 2606 ! CALL iom_put( "mid_x_press", z_p_pmr(:,:,:,2) ) 2607 ! CALL iom_put( "x_mid_slope", z_ddepw_ij(:,:,:,2) ) 2608 ! CALL iom_put( "lhs_x_press", z_p_pmr(:,:,:,1) ) 2609 ! CALL iom_put( "x_lhs_slope", z_ddepw_ij(:,:,:,1) ) 2610 ! CALL iom_put( "rhs_x_press", z_p_pmr(:,:,:,3) ) 2611 ! CALL iom_put( "x_rhs_slope", z_ddepw_ij(:,:,:,3) ) 2589 2612 ELSE 2590 2613 CALL iom_put( "v_force_south", z_F_pmr(:,:,:,1) ) 2591 2614 CALL iom_put( "v_force_upper", z_F_upp ) 2615 ! CALL iom_put( "mid_y_press", z_p_pmr(:,:,:,2) ) 2616 ! CALL iom_put( "y_mid_slope", z_ddepw_ij(:,:,:,2) ) 2592 2617 END IF 2593 2618 … … 2598 2623 !------------------------------------------------------------------------------------------ 2599 2624 2600 SUBROUTINE calc_rhd_ref (k_uv, kn_hor, p rhd_ref, pz_ref, kk_bot_ref)2625 SUBROUTINE calc_rhd_ref (k_uv, kn_hor, pdep, prhd_ref, pz_ref, kk_bot_ref) 2601 2626 2602 2627 !!--------------------------------------------------------------------- … … 2612 2637 INTEGER , INTENT( in ) :: k_uv ! 1 for u-vel; 2 for v_vel 2613 2638 INTEGER , INTENT( in ) :: kn_hor ! 4 for cubic; 2 for linear 2614 2639 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: pdep ! depths (either gdept or gde3w) 2615 2640 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: prhd_ref ! densities of reference profile 2616 2641 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pz_ref ! heights of " " … … 2667 2692 jir = ji_ref(ji,jj) 2668 2693 jjr = jj_ref(ji,jj) 2669 prhd_ref (ji,jj,jk) = 2670 pz_ref (ji,jj,jk) = - gde3w(jir, jjr, jk)2694 prhd_ref (ji,jj,jk) = rhd(jir, jjr, jk) 2695 pz_ref (ji,jj,jk) = - pdep(jir, jjr, jk) 2671 2696 END_3D 2672 2697 … … 2725 2750 !---------------------------------------------------------------------------- 2726 2751 2727 SUBROUTINE calc_mid_ccs( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid )2752 SUBROUTINE calc_mid_ccs( k_uv, k_lev_type, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 2728 2753 2729 2754 !!--------------------------------------------------------------------- … … 2731 2756 !! 2732 2757 !! ** Method : Use constrained cubic spline to interpolate 4 profiles to their central point 2733 !! This version can only be used to interpolate fields on tracer levels (not w-levels)2734 2758 !! 2735 2759 !! ** Action : - set p_fld_mid … … 2737 2761 2738 2762 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2763 INTEGER , INTENT(in) :: k_lev_type ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask 2739 2764 REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation) 2740 2765 REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation) … … 2745 2770 2746 2771 INTEGER ji, jj, jk ! standard loop indices 2747 INTEGER :: j_t_levs ! indicator that field passed is on tracer levels 2748 2772 2773 INTEGER :: j_t_levs, jk_msk ! indicators that field passed in is valid on t-levels or w-levels, jk level to use with masks 2774 2749 2775 j_t_levs = 1 2750 2776 2751 2777 DO jk = 1, jpk 2752 2778 2753 CALL calc_dfld_pmr_ij( k_uv, j_t_levs, jk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, zz_dfld_ij ) 2779 IF ( k_lev_type == j_t_levs ) THEN 2780 jk_msk = jk 2781 ELSE 2782 IF ( jk == 1 ) THEN 2783 jk_msk = 1 2784 ELSE 2785 jk_msk = jk - 1 2786 END IF 2787 END IF 2788 2789 CALL calc_dfld_pmr_ij( k_uv, k_lev_type, jk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, zz_dfld_ij ) 2790 2754 2791 2755 2792 ! first contribution is simple centred average; p_rhd_pmr has 4 points; 2 and 3 are the central ones … … 2760 2797 IF ( k_uv == 1 ) THEN 2761 2798 DO_2D( 0, 0, 0, 0) 2762 IF ( umask(ji-1, jj, jk ) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN2799 IF ( umask(ji-1, jj, jk_msk) > 0.5 .OR. umask(ji+1, jj, jk_msk) > 0.5 ) THEN 2763 2800 p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) ) 2764 2801 END IF … … 2766 2803 ELSE ! k_uv == 2 2767 2804 DO_2D( 0, 0, 0, 0) 2768 IF ( vmask(ji, jj-1, jk ) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN2805 IF ( vmask(ji, jj-1, jk_msk) > 0.5 .OR. vmask(ji, jj+1, jk_msk) > 0.5 ) THEN 2769 2806 p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - 0.125_wp * ( zz_dfld_ij(ji,jj,jk,2) - zz_dfld_ij(ji,jj,jk,1) ) 2770 2807 END IF … … 2774 2811 END DO ! jk 2775 2812 2813 IF ( ln_dbg_hpg ) THEN 2814 CALL dbg_3dr( 'calc_mid_ccs: zz_dfld_ij: 1', zz_dfld_ij(:,:,:,1) ) 2815 CALL dbg_3dr( 'calc_mid_ccs: zz_dfld_ij: 2', zz_dfld_ij(:,:,:,2) ) 2816 END IF 2817 2818 2776 2819 END SUBROUTINE calc_mid_ccs 2777 2820 2778 2821 !---------------------------------------------------------------------------- 2779 2822 2780 SUBROUTINE calc_mid_cub( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid )2823 SUBROUTINE calc_mid_cub( k_uv, k_lev_type, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 2781 2824 2782 2825 !!--------------------------------------------------------------------- … … 2790 2833 2791 2834 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2835 INTEGER , INTENT(in) :: k_lev_type ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask 2792 2836 REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation) 2793 2837 REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation) … … 2802 2846 REAL(wp) z_cor ! correction to the central value 2803 2847 2848 INTEGER :: j_t_levs, jk_msk ! indicators that field passed in is valid on t-levels or w-levels, jk level to use with masks 2849 2850 j_t_levs = 1 2851 2804 2852 z_1_16 = 1._wp / 16._wp ; z_9_16 = 9._wp / 16._wp 2805 2853 2806 2854 DO jk = 1, jpk 2855 2856 IF ( k_lev_type == j_t_levs ) THEN 2857 jk_msk = jk 2858 ELSE 2859 IF ( jk == 1 ) THEN 2860 jk_msk = 1 2861 ELSE 2862 jk_msk = jk - 1 2863 END IF 2864 END IF 2807 2865 2808 2866 ! first contribution is simple centred average plus 1/16 of it; p_rhd_pmr has 4 points; 2 and 3 are the central ones … … 2813 2871 IF ( k_uv == 1 ) THEN 2814 2872 DO_2D( 0, 0, 0, 0) 2815 IF ( umask(ji-1, jj, jk ) > 0.5 ) THEN2873 IF ( umask(ji-1, jj, jk_msk) > 0.5 ) THEN 2816 2874 z_cor = p_fld_pmr(ji,jj,jk,1) 2817 2875 ELSE 2818 2876 z_cor = p_fld_pmr(ji,jj,jk,2) 2819 2877 END IF 2820 IF ( umask(ji+1, jj, jk ) > 0.5 ) THEN2878 IF ( umask(ji+1, jj, jk_msk) > 0.5 ) THEN 2821 2879 z_cor = z_cor + p_fld_pmr(ji,jj,jk,4) 2822 2880 ELSE … … 2827 2885 ELSE ! k_uv == 2 2828 2886 DO_2D( 0, 0, 0, 0) 2829 IF ( vmask(ji, jj-1, jk ) > 0.5 ) THEN2887 IF ( vmask(ji, jj-1, jk_msk) > 0.5 ) THEN 2830 2888 z_cor = p_fld_pmr(ji,jj,jk,1) 2831 2889 ELSE 2832 2890 z_cor = p_fld_pmr(ji,jj,jk,2) 2833 2891 END IF 2834 IF ( vmask(ji, jj+1, jk ) > 0.5 ) THEN2892 IF ( vmask(ji, jj+1, jk_msk) > 0.5 ) THEN 2835 2893 z_cor = z_cor + p_fld_pmr(ji,jj,jk,4) 2836 2894 ELSE … … 2864 2922 2865 2923 INTEGER ji, jj, jk ! standard loop indices 2866 INTEGER jk_msk 2924 INTEGER jk_msk ! jk level to use for umask and vmask (we need z on the upper and lower faces) 2867 2925 INTEGER j_w_levs ! indicator for data on w-levels 2868 2926 2869 2927 j_w_levs = 2 2870 2928 2871 2929 DO jk = 1, jpk 2872 2930 … … 2878 2936 END_2D 2879 2937 2880 ! set the central element if it is not an isolated velocity cell ; this evaluates dz_dij at zeta = 0 ; that is f^{(1)}; and that is given by SMcW(5.8)2938 ! set the central element if it is not an isolated velocity cell 2881 2939 IF ( jk == 1 ) THEN 2882 2940 jk_msk = 1 … … 3120 3178 !------------------------------------------------------------------------------------------ 3121 3179 3122 SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_ e3t, p_p, p_F)3180 SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_dept, p_depw, p_p, p_F) 3123 3181 3124 3182 !!--------------------------------------------------------------------- … … 3126 3184 !! 3127 3185 !! ** Method : Use quadratic reconstruction of p_rhd (density profile) to calculate pressure profile and 3128 !! forces on vertical faces between w levels 3186 !! forces on vertical faces between w levels. Formulated for variable vertical grid spacing. 3129 3187 !! 3130 3188 !! ** Action : - set p_p and p_F … … 3133 3191 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: k_mbkt ! bottom level 3134 3192 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_rhd ! densities on tracer levels 3135 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_e3t ! layer thickness between w levels 3193 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dept ! depths of T points 3194 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_depw ! depths of w points (top & bottom of T cells) 3136 3195 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_p ! pressures on w levels 3137 3196 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_F ! force on the vertical face between w levels 3138 3197 3139 3198 INTEGER :: ji, jj, jk ! dummy loop indices 3140 REAL(wp) :: z_2_3, z_4_3, z_8_3 ! 2/3 and 4/3 and 8/3 3141 REAL(wp) :: z_22_3, z_28_3 ! 22/3 and 28/3 3142 3143 REAL(wp) :: zr_a, zr_b, zr_c ! rhd evaluated at i, i+1/2 and i+1 3144 REAL(wp) :: za_0, za_1, za_2 ! polynomial coefficients 3145 REAL(wp) :: ze3t ! 3146 3147 z_2_3 = 2.0_wp/3.0_wp 3148 z_4_3 = 4.0_wp/3.0_wp 3149 z_8_3 = 8.0_wp/3.0_wp 3150 z_22_3 = 22.0_wp/3.0_wp 3151 z_28_3 = 28.0_wp/3.0_wp 3199 REAL(wp) :: z_1_3 ! 1/3 3200 3201 INTEGER, DIMENSION(A2D(nn_hls)) :: jklev 3202 REAL(wp), DIMENSION(A2D(nn_hls)) :: zr_a, zr_0, zr_b ! rhd values (r) at points a (above), 0 (central), b (below) 3203 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_a, zz_0, zz_b ! heights at points a (above), 0 (central), b (below) 3204 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_upp, zz_low ! upper & lower w levels 3205 INTEGER :: jkt ! level being worked on 3206 REAL(wp) :: zet_a, zet_0, zet_b ! zeta values for above (a), central (0) and below (b) 3207 REAL(wp) :: zr_d_a, zr_d_0, zr_d_b ! reciprocals of denominators for terms involving r_a, r_0 and r_b 3208 REAL(wp) :: za_0, za_1, za_2 ! quadratic polynomial coefficients for r (rhd) 3209 REAL(wp) :: zet_upp, zet_low ! zeta for upper and lower w levels 3210 REAL(wp) :: zet_upp_sq, zet_low_sq ! squares of above values 3211 REAL(wp) :: zp_upp, zp_low ! contributions to the pressure 3212 REAL(wp) :: zf_uml ! force "upper minus lower" 3213 3214 z_1_3 = 1.0_wp/3.0_wp 3152 3215 3153 3216 ! Integrate densities in the vertical using quadratic fits to the densities 3154 3217 ! zr_a is r_{-1} zr_b is r_1 ; zr_c is r_3 3155 3218 3219 DO_2D(0,0,0,0) 3220 p_p(ji,jj,1) = 0.0_wp 3221 END_2D 3222 3156 3223 ! one-sided quadratic at the upper boundary. pressure is zero at the upper surface 3157 DO_2D(0,0,0,0) 3158 zr_a = p_rhd(ji,jj,1) 3159 zr_b = p_rhd(ji,jj,2) 3160 zr_c = p_rhd(ji,jj,3) 3161 ze3t = p_e3t(ji,jj,1) 3162 za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 3163 za_1 = 0.5_wp * ( zr_b - zr_a ) 3164 za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a ) 3165 p_p(ji,jj,1) = 0.0_wp 3166 p_p(ji,jj,2) = grav*ze3t* ( za_0 - za_1 + z_4_3*za_2 ) 3167 p_F(ji,jj,1) = 0.5_wp*grav*ze3t*ze3t*( za_0 - z_4_3*za_1 + 2._wp*za_2 ) 3168 END_2D 3169 3170 ! centred quadratic in the interior 3171 DO_3D(0,0,0,0,2,jpk-1) 3172 zr_a = p_rhd(ji,jj,jk-1) 3173 zr_b = p_rhd(ji,jj,jk) 3174 zr_c = p_rhd(ji,jj,jk+1) 3175 ze3t = p_e3t(ji,jj,jk) 3176 za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 3177 za_1 = 0.5_wp * ( zr_b - zr_a ) 3178 za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a ) 3179 p_p(ji,jj,jk+1) = p_p(ji,jj,jk) + grav*ze3t* ( za_0 + za_1 + z_4_3*za_2 ) 3180 p_F(ji,jj,jk ) = ze3t*p_p(ji,jj,jk) + 0.5*grav*ze3t*ze3t*( za_0 + z_2_3*za_1 + z_2_3*za_2 ) 3181 END_3D 3182 3183 ! one-sided quadratic at the lower boundary 3184 DO_2D(0,0,0,0) 3185 jk = k_mbkt(ji,jj) 3186 zr_a = p_rhd(ji,jj,jk-2) 3187 zr_b = p_rhd(ji,jj,jk-1) 3188 zr_c = p_rhd(ji,jj,jk) 3189 ze3t = p_e3t(ji,jj,jk) 3190 za_0 = 0.125_wp * ( - zr_c + 6._wp*zr_b + 3._wp*zr_a ) 3191 za_1 = 0.5_wp * ( zr_b - zr_a ) 3192 za_2 = 0.125_wp * ( zr_c - 2._wp*zr_b + zr_a ) 3193 p_p(ji,jj,jk+1) = p_p(ji,jj,jk) + grav*ze3t* ( za_0 + 3.*za_1 + z_28_3*za_2 ) 3194 p_F(ji,jj,jk ) = ze3t*p_p(ji,jj,jk) + 0.5*grav*ze3t*ze3t*( za_0 + z_8_3*za_1 + z_22_3*za_2 ) 3195 END_2D 3196 3197 RETURN 3224 DO jk = 1, jpk 3225 3226 IF (jk == 1) THEN ! perform off-centred calculations for the first level 3227 DO_2D(0,0,0,0) 3228 jklev(ji,jj) = jk 3229 zr_a(ji,jj) = p_rhd(ji,jj,1) 3230 zr_0(ji,jj) = p_rhd(ji,jj,2) 3231 zr_b(ji,jj) = p_rhd(ji,jj,3) 3232 zz_a(ji,jj) = - p_dept(ji,jj,1) 3233 zz_0(ji,jj) = - p_dept(ji,jj,2) 3234 zz_b(ji,jj) = - p_dept(ji,jj,3) 3235 zz_upp(ji,jj) = - p_depw(ji,jj,1) 3236 zz_low(ji,jj) = - p_depw(ji,jj,2) 3237 END_2D 3238 3239 ELSE IF ( jk < jpk ) THEN 3240 DO_2D(0,0,0,0) 3241 jklev(ji,jj) = jk 3242 zr_a(ji,jj) = p_rhd(ji,jj,jk-1) 3243 zr_0(ji,jj) = p_rhd(ji,jj,jk) 3244 zr_b(ji,jj) = p_rhd(ji,jj,jk+1) 3245 zz_a(ji,jj) = - p_dept(ji,jj,jk-1) 3246 zz_0(ji,jj) = - p_dept(ji,jj,jk) 3247 zz_b(ji,jj) = - p_dept(ji,jj,jk+1) 3248 zz_upp(ji,jj) = - p_depw(ji,jj,jk) 3249 zz_low(ji,jj) = - p_depw(ji,jj,jk+1) 3250 END_2D 3251 3252 ELSE IF ( jk == jpk ) THEN ! repeat calculations at bottom level using off-centred quadratic 3253 3254 DO_2D(0,0,0,0) 3255 jkt = k_mbkt(ji,jj) 3256 jklev(ji,jj) = jkt 3257 zr_a(ji,jj) = p_rhd(ji,jj,jkt-2) 3258 zr_0(ji,jj) = p_rhd(ji,jj,jkt-1) 3259 zr_b(ji,jj) = p_rhd(ji,jj,jkt) 3260 zz_a(ji,jj) = - p_dept(ji,jj,jkt-2) 3261 zz_0(ji,jj) = - p_dept(ji,jj,jkt-1) 3262 zz_b(ji,jj) = - p_dept(ji,jj,jkt) 3263 zz_upp(ji,jj) = - p_depw(ji,jj,jkt) 3264 zz_low(ji,jj) = - p_depw(ji,jj,jkt+1) 3265 END_2D 3266 3267 END IF ! jk 3268 3269 DO_2D(0,0,0,0) 3270 jkt = jklev(ji,jj) 3271 3272 zet_a = - ( zz_a(ji,jj) - zz_0(ji,jj) ) 3273 zet_0 = 0.0_wp 3274 zet_b = - ( zz_b(ji,jj) - zz_0(ji,jj) ) 3275 3276 zr_d_a = 1.0_wp / ( zet_a*(zet_a - zet_b) ) 3277 zr_d_0 = 1.0_wp / ( zet_a*zet_b ) 3278 zr_d_b = 1.0_wp / ( zet_b*(zet_b - zet_a) ) 3279 3280 za_0 = zr_0(ji,jj) 3281 za_1 = - zr_d_0*zr_0(ji,jj)*(zet_a+zet_b) - zr_d_a*zr_a(ji,jj)*zet_b - zr_d_b*zr_b(ji,jj)*zet_a 3282 za_2 = zr_d_0*zr_0(ji,jj) + zr_d_a*zr_a(ji,jj) + zr_d_b*zr_b(ji,jj) 3283 3284 zet_upp = - ( zz_upp(ji,jj) - zz_0(ji,jj) ) 3285 zet_low = - ( zz_low(ji,jj) - zz_0(ji,jj) ) 3286 3287 zet_upp_sq = zet_upp*zet_upp ; zet_low_sq = zet_low*zet_low 3288 3289 zp_upp = za_0*zet_upp + 0.5_wp*za_1*zet_upp_sq + z_1_3*za_2*zet_upp_sq*zet_upp 3290 zp_low = za_0*zet_low + 0.5_wp*za_1*zet_low_sq + z_1_3*za_2*zet_low_sq*zet_low 3291 3292 p_p(ji,jj,jkt+1) = p_p(ji,jj,jkt) + grav * ( zp_low - zp_upp ) 3293 3294 zf_uml = za_0 * ( 0.5_wp *(zet_low_sq - zet_upp_sq) - zet_upp *(zet_low-zet_upp) ) & 3295 & + 0.5_wp*za_1 * ( z_1_3 *(zet_low_sq*zet_low - zet_upp_sq*zet_upp) - zet_upp_sq *(zet_low-zet_upp) ) & 3296 & + z_1_3 *za_2 * ( 0.25_wp*(zet_low_sq*zet_low_sq - zet_upp_sq*zet_upp_sq) - zet_upp_sq*zet_upp*(zet_low-zet_upp) ) 3297 3298 ! the force on the face is positive (the opposite sign from my notes) 3299 p_F(ji,jj,jkt) = p_p(ji,jj,jkt)*( zz_upp(ji,jj) - zz_low(ji,jj) ) + grav * zf_uml 3300 3301 END_2D 3302 3303 END DO ! jk 3304 3198 3305 END SUBROUTINE vrt_int_quad 3199 3306
Note: See TracChangeset
for help on using the changeset viewer.
