Changeset 15728
- Timestamp:
- 2022-02-23T16:55:36+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
r15727 r15728 2 2 3 3 !! v2: takes into account the stretching of the vertical grid in hpg_ffr 4 4 !! v3: takes into account the stretching of the vertical grid in the calculation of the reference profile 5 5 6 !!====================================================================== 6 7 !! *** MODULE dynhpg *** … … 94 95 REAL(wp) :: bco_bc_z_srf, aco_bc_z_bot, bco_bc_z_bot ! " 95 96 96 LOGICAL :: ln_hpg_djr_ref_ccs ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles 97 LOGICAL :: ln_hpg_ffr_ref ! T => reference profile is subtracted; F => no reference profile subtracted 98 LOGICAL :: ln_hpg_ffr_ref_ccs ! T => use constrained cubic spline to vertically interpolate reference to each profile 97 LOGICAL :: ln_hpg_ref ! T => reference profile is subtracted; F => no reference profile subtracted 98 LOGICAL :: ln_hpg_ref_str ! T => reference profile calculated taking into account vertical variation in grid spacing 99 LOGICAL :: ln_hpg_ref_ccs ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles 100 LOGICAL :: ln_hpg_ref_off ! only applies if ln_hpg_ref_ccs = T; T => use off-centred simple cubic at upper & lower boundaries 101 99 102 LOGICAL :: ln_hpg_ffr_hor_ccs ! T => use constrained cubic spline to interpolate z_rhd_pmr along upper faces of cells 100 103 LOGICAL :: ln_hpg_ffr_hor_cub ! T => use simple cubic polynomial to interpolate z_rhd_pmr along upper faces of cells … … 195 198 & ln_hpg_bcvN_rhd_hor, ln_hpg_bcvN_rhd_srf, & 196 199 & ln_hpg_bcvN_rhd_bot, ln_hpg_bcvN_z_hor, & 197 & ln_hpg_ djr_ref_ccs, ln_hpg_ffr_vrt_quad, &198 & ln_hpg_ ffr_ref, ln_hpg_ffr_ref_ccs, &200 & ln_hpg_ref, ln_hpg_ref_str, & 201 & ln_hpg_ref_ccs, ln_hpg_ref_off, & 199 202 & ln_hpg_ffr_hor_ccs, ln_hpg_ffr_hor_cub, & 203 & ln_hpg_ffr_vrt_quad, & 200 204 & ln_dbg_hpg, ln_dbg_ij, ln_dbg_ik, ln_dbg_jk, & 201 205 & ki_dbg_min, ki_dbg_max, ki_dbg_cen, & … … 322 326 END IF 323 327 324 IF ( lwp .AND. ln_hpg_djr ) THEN 325 IF ( ln_hpg_djr_ref_ccs ) THEN 326 WRITE(numout,*) ' interpolation of reference profile by constrained cubic spline' 327 ELSE 328 IF ( lwp .AND. ln_hpg_ref ) THEN 329 IF ( ln_hpg_ref_ccs .AND. ln_hpg_ref_off .AND. ln_hpg_ref_str ) THEN 330 WRITE(numout,*) ' interpolation of reference profile by constrained cubic spline with off-centring at boundaries' 331 ELSE IF ( ln_hpg_ref_ccs ) THEN 332 WRITE(numout,*) ' interpolation of reference profile by constrained cubic spline using boundary conditions' 333 ELSE 328 334 WRITE(numout,*) ' interpolation of reference profile by simple cubic spline with off-centring at boundaries' 329 335 END IF 336 IF ( ln_hpg_ref_str ) THEN 337 WRITE(numout,*) ' interpolation of reference profile takes account of variation in vertical grid spacing' 338 ELSE 339 WRITE(numout,*) ' interpolation of reference profile does NOT take account of variation in vertical grid spacing' 340 END IF 330 341 END IF 331 342 332 IF ( lwp .AND. ln_hpg_ffr ) THEN 333 IF ( ln_hpg_ffr_ref_ccs ) THEN 334 IF ( ln_hpg_ffr_ref_ccs ) THEN 335 WRITE(numout,*) ' interpolation of reference profile in vertical by constrained cubic spline ' 336 ELSE 337 WRITE(numout,*) ' interpolation of reference profile in vertical by pure cubic polynomial fit ' 338 END IF 339 ELSE 340 WRITE(numout,*) ' no reference profile subtracted ' 343 IF ( lwp .AND. ln_hpg_ffr ) THEN 344 IF ( .NOT. ln_hpg_ref ) THEN 345 WRITE(numout,*) ' no reference profile subtracted ' 341 346 END IF 342 347 IF ( ln_hpg_ffr_hor_ccs ) THEN … … 1751 1756 1752 1757 !-------------------------------------------------------------------------------------------------------- 1753 ! 2.2 IF ln_hpg_ djr_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom1758 ! 2.2 IF ln_hpg_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom 1754 1759 ! (bcs not needed for simple cubic off-centred at boundaries) 1755 1760 !-------------------------------------------------------------------------------------------------------- 1756 1761 1757 IF ( ln_hpg_ djr_ref_ccs ) THEN1762 IF ( ln_hpg_ref_ccs ) THEN 1758 1763 CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 1759 1764 IF ( ln_dbg_hpg ) CALL dbg_3dr( '2.3 zdrhd_k_ref', zdrhd_k_ref ) 1760 END IF ! ln_hpg_ djr_ref_ccs1765 END IF ! ln_hpg_ref_ccs 1761 1766 1762 1767 !---------------------------------------------------------------------------------------- … … 1773 1778 END IF 1774 1779 1775 IF ( ln_hpg_djr_ref_ccs ) THEN 1776 CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1777 ELSE 1778 CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1780 1781 IF ( ln_hpg_ref ) THEN 1782 IF ( ln_hpg_ref_str ) THEN 1783 IF ( ln_hpg_ref_ccs ) THEN 1784 IF ( ln_hpg_ref_off ) THEN 1785 CALL ref_to_tgt_ccs_str_off ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1786 ELSE 1787 CALL ref_to_tgt_ccs_str ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1788 END IF 1789 ELSE 1790 CALL ref_to_tgt_cub_str ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1791 END IF 1792 ELSE ! these calls are mainly retained to assist in testing 1793 IF ( ln_hpg_ref_ccs ) THEN 1794 CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1795 ELSE 1796 CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1797 END IF 1798 END IF 1779 1799 END IF 1780 1800 1781 1801 IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) ) 1782 1802 … … 2045 2065 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid 2046 2066 2067 !--------------------------------------------------------------------------------------------------------------- 2068 2047 2069 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid; target lies between jk_ref and jk_ref-1. 2070 LOGICAL, DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen ! T => off-centred interpolation (values not used in this routine) 2048 2071 2049 2072 !--------------------------------------------------------------------------------------------------------------- … … 2068 2091 2069 2092 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 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 )2093 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, ll_tgt_off_cen ) 2071 2094 2072 2095 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) … … 2106 2129 2107 2130 END SUBROUTINE ref_to_tgt_cub 2131 2132 !--------------------------------------------------------------------------------------------------------------- 2133 2134 SUBROUTINE ref_to_tgt_cub_str ( ki_off_tgt, kj_off_tgt, p_dep_tgt, p_fld_tgt, p_z_ref, p_fld_ref, kk_bot_ref, p_fld_tgt_ref) 2135 2136 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2137 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2138 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dep_tgt ! depths of target profiles 2139 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_tgt ! field values on the target grid 2140 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_z_ref ! heights of reference profiles 2141 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_ref ! field values to be interpolated (in the vertical) on reference grid 2142 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2143 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid 2144 2145 !--------------------------------------------------------------------------------------------------------------- 2146 2147 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid; target lies between jk_ref and jk_ref-1. 2148 LOGICAL, DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen ! T => off-centred interpolation (values not used in this routine) 2149 2150 !--------------------------------------------------------------------------------------------------------------- 2151 2152 INTEGER :: ji, jj, jk ! loop indices 2153 INTEGER :: jkr 2154 2155 REAL(wp) :: zf_p, zf_0, zf_m, zf_b ! plus, zero, minus, below ; note that zeta increases with depth ; zero is zeta = 1/2 2156 REAL(wp) :: zz_p, zz_0, zz_m, zz_b 2157 REAL(wp) :: zs_p, zs_0, zs_m, zs_b 2158 REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 2159 REAL(wp) :: za_1, za_2, za_3 2160 2161 REAL(wp) :: zz_tgt_lcl, zz_tgt_sq, zfld_ref_on_tgt 2162 LOGICAL :: ll_ccs 2163 2164 !----------------------------------------------------------------------------------------------------------------- 2165 2166 ll_ccs = .FALSE. ! set for the call to loc_ref_tgt 2167 2168 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 2169 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, ll_tgt_off_cen ) 2170 2171 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 2172 2173 ! it would probably be better computationally for fld_ref to have the jk index first. 2174 2175 !!! jkr >= 2 and p_fld_ref has jk = 0 available 2176 2177 jkr = jk_ref_for_tgt(ji,jj,jk) 2178 zf_b = p_fld_ref(ji,jj,jkr-2) 2179 zf_m = p_fld_ref(ji,jj,jkr-1) 2180 zf_0 = p_fld_ref(ji,jj,jkr ) 2181 zf_p = p_fld_ref(ji,jj,jkr+1) 2182 2183 zz_b = p_z_ref( ji, jj, jkr-2) 2184 zz_m = p_z_ref( ji, jj, jkr-1) 2185 zz_0 = p_z_ref( ji, jj, jkr ) 2186 zz_b = zz_b - zz_0 2187 zz_m = zz_m - zz_0 2188 zz_p = p_z_ref( ji, jj, jkr+1) - zz_0 2189 2190 zz_p_m = zz_p - zz_m 2191 zz_p_b = zz_p - zz_b 2192 zz_m_b = zz_m - zz_b 2193 2194 zs_0 = - zf_0 / ( zz_p * zz_m * zz_b ) 2195 2196 zs_p = zf_p / ( zz_p * zz_p_m * zz_p_b ) 2197 zs_b = zf_b / ( zz_b * zz_p_b * zz_m_b ) 2198 zs_m = - zf_m / ( zz_m * zz_p_m * zz_m_b ) 2199 2200 za_1 = zz_m*zz_b *zs_p + zz_p*zz_b *zs_m + zz_p*zz_m *zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0 2201 za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 2202 za_3 = zs_p + zs_m + zs_0 + zs_b 2203 2204 zz_tgt_lcl = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0 2205 zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl 2206 2207 zfld_ref_on_tgt = zf_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 2208 2209 ! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc. 2210 p_fld_tgt_ref(ji, jj, jk) = p_fld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt 2211 2212 END_3D 2213 2214 END SUBROUTINE ref_to_tgt_cub_str 2108 2215 2109 2216 !----------------------------------------------------------------------------------------------------------------- … … 2123 2230 !----------------------------------------------------------------------------------------------------------------- 2124 2231 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid 2232 LOGICAL, DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen ! T => off-centred interpolation (values not used in this routine) 2125 2233 2126 2234 INTEGER :: ji, jj, jk ! DO loop indices 2127 2235 INTEGER :: jkr 2128 REAL(wp) :: z_r_6 ! constant2129 2236 REAL(wp) :: zz_tgt_lcl, zz_ref_jkrm1, zz_ref_jkr, zeta, zetasq 2130 2237 REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave 2131 2238 REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3 2132 2239 REAL(wp) :: zfld_ref_on_tgt 2133 LOGICAL :: ll_ccs ! set to . FALSE. for the call to loc_ref_tgt2240 LOGICAL :: ll_ccs ! set to .TRUE. for the call to loc_ref_tgt 2134 2241 2135 2242 !----------------------------------------------------------------------------------------------------------------- 2136 2243 2137 2244 ll_ccs = .TRUE. ! set for the call to loc_ref_tgt 2138 z_r_6 = 1._wp / 6._wp2139 2245 2140 2246 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 2141 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt )2247 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen ) 2142 2248 2143 2249 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) … … 2182 2288 !----------------------------------------------------------------------------------------------------------------- 2183 2289 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 ) 2290 SUBROUTINE ref_to_tgt_ccs_str ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref) 2291 2292 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2293 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2294 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdep_tgt ! depths of target profiles 2295 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_tgt ! field values on the target grid 2296 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pz_ref ! heights of reference profiles 2297 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_ref ! reference field values 2298 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdfld_k_ref ! harmonic averages of vertical differences of reference field 2299 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2300 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: pfld_tgt_ref ! target minus reference on target grid 2301 2302 !----------------------------------------------------------------------------------------------------------------- 2303 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid 2304 LOGICAL, DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen ! T => off-centred interpolation (values not used in this routine) 2305 2306 INTEGER :: ji, jj, jk ! DO loop indices 2307 INTEGER :: jkr 2308 REAL(wp) :: zz_tgt_lcl, zz_tgt_sq 2309 2310 REAL(wp) :: zeta_p, zeta_0, zeta_m, zeta_b ! plus, zero, minus, below (note that zeta increases with depth) zero is zeta = 1/2 2311 REAL(wp) :: zz_p, zz_0, zz_m, zz_b 2312 REAL(wp) :: zs_p, zs_0, zs_m, zs_b 2313 REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 2314 REAL(wp) :: za_1, za_2, za_3 2315 REAL(wp) :: zeta, zetasq 2316 2317 REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave 2318 REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3 2319 REAL(wp) :: zfld_ref_on_tgt 2320 LOGICAL :: ll_ccs ! set to .FALSE. for the call to loc_ref_tgt 2321 2322 !----------------------------------------------------------------------------------------------------------------- 2323 2324 ll_ccs = .TRUE. ! set for the call to loc_ref_tgt 2325 2326 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point 2327 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen ) 2328 2329 zeta_p = 1.5_wp 2330 zeta_0 = 0.5_wp 2331 zeta_m = -0.5_wp 2332 zeta_b = -1.5_wp 2333 2334 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 2335 2336 jkr = jk_ref_for_tgt( ji, jj, jk ) 2337 2338 zz_0 = pz_ref( ji, jj, jkr ) 2339 zz_m = pz_ref( ji, jj, jkr-1) - zz_0 2340 2341 IF ( jkr .NE. 2) THEN 2342 zz_b = pz_ref( ji, jj, jkr-2) - zz_0 2343 ELSE 2344 zz_b = 2._wp * zz_m 2345 END IF 2346 2347 IF ( ll_tgt_off_cen(ji,jj,jk) ) THEN 2348 zz_p = - zz_m 2349 ELSE 2350 zz_p = pz_ref( ji, jj, jkr+1) - zz_0 2351 END IF 2352 2353 zz_p_m = zz_p - zz_m 2354 zz_p_b = zz_p - zz_b 2355 zz_m_b = zz_m - zz_b 2356 2357 zs_0 = - zeta_0 / ( zz_p * zz_m * zz_b ) 2358 2359 zs_p = zeta_p / ( zz_p * zz_p_m * zz_p_b ) 2360 zs_m = - zeta_m / ( zz_m * zz_p_m * zz_m_b ) 2361 zs_b = zeta_b / ( zz_b * zz_p_b * zz_m_b ) 2362 2363 za_1 = zz_m*zz_b*zs_p + zz_p*zz_b*zs_m + zz_p*zz_m*zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0 2364 za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 2365 za_3 = zs_p + zs_m + zs_0 + zs_b 2366 2367 zz_tgt_lcl = - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0 2368 zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl 2369 2370 zeta = zeta_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 2371 zetasq = zeta*zeta 2372 2373 zd_dif = pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1) 2374 zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) ) 2375 2376 zf_dif = pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1) 2377 zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) ) 2378 2379 zcoef_0 = zf_ave - 0.125_wp * zd_dif 2380 zcoef_1 = 1.5_wp * zf_dif - 0.5_wp * zd_ave 2381 zcoef_2 = 0.5_wp * zd_dif 2382 zcoef_3 = 2.0_wp * zd_ave - 2._wp * zf_dif 2383 2384 zfld_ref_on_tgt = zcoef_0 + zeta*zcoef_1 + zetasq * ( zcoef_2 + zeta*zcoef_3 ) 2385 2386 ! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc. 2387 2388 pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt 2389 2390 ! IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN 2391 ! WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave 2392 ! WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 2393 ! WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) 2394 ! WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt 2395 ! WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) 2396 ! END IF 2397 2398 END_3D 2399 2400 IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs_str: pfld_tgt_ref', pfld_tgt_ref ) 2401 2402 END SUBROUTINE ref_to_tgt_ccs_str 2403 2404 !----------------------------------------------------------------------------------------------------------------- 2405 2406 SUBROUTINE ref_to_tgt_ccs_str_off ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref) 2407 2408 ! Coded for a stretched grid and to perform off-centred pure cubic interpolation at the upper and lower boundaries 2409 2410 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2411 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2412 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdep_tgt ! depths of target profiles 2413 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_tgt ! field values on the target grid 2414 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pz_ref ! heights of reference profiles 2415 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_ref ! reference field values 2416 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdfld_k_ref ! harmonic averages of vertical differences of reference field 2417 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2418 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: pfld_tgt_ref ! target minus reference on target grid 2419 2420 !----------------------------------------------------------------------------------------------------------------- 2421 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid 2422 LOGICAL, DIMENSION(A2D(nn_hls),jpk) :: ll_tgt_off_cen ! T => off-centred interpolation 2423 2424 INTEGER :: ji, jj, jk ! DO loop indices 2425 INTEGER :: jkr 2426 REAL(wp) :: zz_tgt_lcl, zz_tgt_sq 2427 2428 REAL(wp) :: zz_p, zz_0, zz_m, zz_b 2429 REAL(wp) :: zeta_p, zeta_0, zeta_m, zeta_b 2430 REAL(wp) :: zf_p, zf_0, zf_m, zf_b 2431 REAL(wp) :: zs_p, zs_0, zs_m, zs_b 2432 REAL(wp) :: zz_p_m, zz_p_b, zz_m_b 2433 REAL(wp) :: za_1, za_2,za_3 2434 REAL(wp) :: zeta, zetasq 2435 2436 REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave 2437 REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3 2438 REAL(wp) :: zfld_ref_on_tgt 2439 LOGICAL :: ll_ccs ! set to .FALSE. for the call to loc_ref_tgt 2440 2441 !----------------------------------------------------------------------------------------------------------------- 2442 2443 ll_ccs = .FALSE. ! set for the call to loc_ref_tgt 2444 2445 ! find jk_ref_for_tgt (bounding levels on reference grid for each target point) and ll_tgt_off_cen 2446 CALL loc_ref_tgt ( ll_ccs, ki_off_tgt, kj_off_tgt, pdep_tgt, pz_ref, kk_bot_ref, jk_ref_for_tgt, ll_tgt_off_cen ) 2447 2448 zeta_p = 1.5_wp 2449 zeta_0 = 0.5_wp 2450 zeta_m = -0.5_wp 2451 zeta_b = -1.5_wp 2452 2453 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 2454 2455 jkr = jk_ref_for_tgt( ji, jj, jk ) 2456 2457 zz_b = pz_ref( ji, jj, jkr-2) 2458 zz_m = pz_ref( ji, jj, jkr-1) 2459 zz_0 = pz_ref( ji, jj, jkr ) 2460 zz_b = zz_b - zz_0 2461 zz_m = zz_m - zz_0 2462 zz_p = pz_ref( ji, jj, jkr+1) - zz_0 2463 2464 zz_tgt_lcl = - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) - zz_0 2465 zz_tgt_sq = zz_tgt_lcl*zz_tgt_lcl 2466 2467 IF ( ll_tgt_off_cen(ji,jj,jk) ) THEN 2468 2469 ! off centred pure cubic fit 2470 2471 zf_b = pfld_ref(ji,jj,jkr-2) 2472 zf_m = pfld_ref(ji,jj,jkr-1) 2473 zf_0 = pfld_ref(ji,jj,jkr ) 2474 zf_p = pfld_ref(ji,jj,jkr+1) 2475 2476 zz_p_m = zz_p - zz_m 2477 zz_p_b = zz_p - zz_b 2478 zz_m_b = zz_m - zz_b 2479 2480 zs_0 = - zf_0 / ( zz_p * zz_m * zz_b ) 2481 2482 zs_p = zf_p / ( zz_p * zz_p_m * zz_p_b ) 2483 zs_b = zf_b / ( zz_b * zz_p_b * zz_m_b ) 2484 zs_m = - zf_m / ( zz_m * zz_p_m * zz_m_b ) 2485 2486 za_1 = zz_m*zz_b *zs_p + zz_p*zz_b *zs_m + zz_p*zz_m *zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0 2487 za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 2488 za_3 = zs_p + zs_m + zs_0 + zs_b 2489 2490 zfld_ref_on_tgt = zf_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 2491 2492 ELSE 2493 2494 zz_p_m = zz_p - zz_m 2495 zz_p_b = zz_p - zz_b 2496 zz_m_b = zz_m - zz_b 2497 2498 zs_0 = - zeta_0 / ( zz_p * zz_m * zz_b ) 2499 2500 zs_p = zeta_p / ( zz_p * zz_p_m * zz_p_b ) 2501 zs_m = - zeta_m / ( zz_m * zz_p_m * zz_m_b ) 2502 zs_b = zeta_b / ( zz_b * zz_p_b * zz_m_b ) 2503 2504 za_1 = zz_m*zz_b*zs_p + zz_p*zz_b*zs_m + zz_p*zz_m*zs_b + (zz_p*zz_m+zz_m*zz_b+zz_p*zz_b)*zs_0 2505 za_2 = - (zz_m+zz_b)*zs_p - (zz_p+zz_b)*zs_m - (zz_p+zz_m)*zs_b - (zz_p+zz_m+zz_b)*zs_0 2506 za_3 = zs_p + zs_m + zs_0 + zs_b 2507 2508 zeta = zeta_0 + za_1*zz_tgt_lcl + za_2*zz_tgt_sq + za_3*zz_tgt_lcl*zz_tgt_sq 2509 zetasq = zeta*zeta 2510 2511 zd_dif = pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1) 2512 zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) ) 2513 2514 zf_dif = pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1) 2515 zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) ) 2516 2517 zcoef_0 = zf_ave - 0.125_wp * zd_dif 2518 zcoef_1 = 1.5_wp * zf_dif - 0.5_wp * zd_ave 2519 zcoef_2 = 0.5_wp * zd_dif 2520 zcoef_3 = 2.0_wp * zd_ave - 2._wp * zf_dif 2521 2522 zfld_ref_on_tgt = zcoef_0 + zeta*zcoef_1 + zetasq * ( zcoef_2 + zeta*zcoef_3 ) 2523 2524 END IF 2525 2526 ! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc. 2527 2528 pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt 2529 2530 ! IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN 2531 ! WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave 2532 ! WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 2533 ! WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) 2534 ! WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt 2535 ! WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) 2536 ! END IF 2537 2538 END_3D 2539 2540 IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs_str_off: pfld_tgt_ref', pfld_tgt_ref ) 2541 2542 END SUBROUTINE ref_to_tgt_ccs_str_off 2543 2544 !----------------------------------------------------------------------------------------------------------------- 2545 2546 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, ll_tgt_off_cen ) 2185 2547 2186 2548 LOGICAL, INTENT(in) :: ll_ccs ! true => constrained cubic spline; false => pure cubic fit … … 2191 2553 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2192 2554 INTEGER, DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: kk_ref_for_tgt ! reference index for interpolation to target grid (the lower point) 2193 ! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr 2555 ! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr 2556 LOGICAL, DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: ll_tgt_off_cen ! T => target point is off-centred (only applies when ll_ccs is False) 2557 2194 2558 !----------------------------------------------------------------------------------------------------------------- 2195 2559 … … 2211 2575 ! 1. Initialisation for search for points on reference grid bounding points on the target grid 2212 2576 ! the first point on the target grid is assumed to lie between the first two points on the reference grid 2213 IF ( ll_ccs ) THEN 2214 jk_init = 2 2215 ELSE 2216 jk_init = 3 ! for simple cubic use off-centred interpolation near the surface 2217 END IF 2577 2578 jk_init = 2 ! for all cases; enforcement of off-centred interpolation is now done at the end of the routine 2218 2579 2219 2580 jk_ref(:,:) = jk_init 2220 2581 kk_ref_for_tgt(:,:,1) = jk_init 2221 2582 jk_tgt(:,:) = 2 2583 2584 ll_tgt_off_cen(:,:,:) = .FALSE. 2222 2585 2223 2586 ! 2. Find kk_ref_for_tgt … … 2322 2685 ! This assumes that every sea point has at least 4 levels. 2323 2686 2324 IF ( .NOT.ll_ccs ) THEN2687 IF ( ll_ccs ) THEN 2325 2688 DO_3D(0, 0, 0, 0, 2, jpk-1) 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 2689 IF ( kk_ref_for_tgt(ji,jj,jk) == kk_bot_ref(ji,jj) ) THEN 2690 ll_tgt_off_cen(ji,jj,jk) = .TRUE. 2691 END IF 2327 2692 END_3D 2693 ELSE 2694 DO_3D(0, 0, 0, 0, 1, jpk-1) 2695 IF ( kk_ref_for_tgt(ji,jj,jk) == 2 ) THEN 2696 ll_tgt_off_cen(ji,jj,jk) = .TRUE. 2697 kk_ref_for_tgt(ji,jj,jk) = 3 2698 END IF 2699 IF ( kk_ref_for_tgt(ji,jj,jk) == kk_bot_ref(ji,jj) ) THEN 2700 ll_tgt_off_cen(ji,jj,jk) = .TRUE. 2701 kk_ref_for_tgt(ji,jj,jk) = kk_bot_ref(ji,jj) - 1 2702 END IF 2703 END_3D 2704 END IF 2705 2706 IF ( lwp .AND. ln_dbg_hpg ) THEN 2707 WRITE( numout,*) 'loc_ref_tgt at end of section 4 ', ki_off_tgt, kj_off_tgt 2708 CALL dbg_3di( 'kk_ref_for_tgt', kk_ref_for_tgt ) 2328 2709 END IF 2329 2710 … … 2411 2792 ! 1. Calculate the reference profile from all points in the stencil 2412 2793 2413 IF ( ln_hpg_ ffr_ref ) THEN2794 IF ( ln_hpg_ref ) THEN 2414 2795 CALL calc_rhd_ref(j_uv, jn_hor_pts, gdept(:,:,:,Kmm), zrhd_ref, zz_ref, jk_bot_ref) 2415 IF ( ln_hpg_ ffr_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref)2796 IF ( ln_hpg_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 2416 2797 END IF 2417 2798 … … 2434 2815 END IF 2435 2816 2436 IF ( ln_hpg_ffr_ref ) THEN 2437 2438 IF ( ln_hpg_ffr_ref_ccs ) THEN 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) ) 2440 ELSE 2441 CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2442 END IF 2817 IF ( ln_hpg_ref ) THEN 2818 IF ( ln_hpg_ref_str ) THEN 2819 IF ( ln_hpg_ref_ccs ) THEN 2820 IF ( ln_hpg_ref_off ) THEN 2821 CALL ref_to_tgt_ccs_str_off ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2822 ELSE 2823 CALL ref_to_tgt_ccs_str ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2824 END IF 2825 ELSE 2826 CALL ref_to_tgt_cub_str ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2827 END IF 2828 ELSE ! these calls are mainly retained to assist in testing 2829 IF ( ln_hpg_ref_ccs ) THEN 2830 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) ) 2831 ELSE 2832 CALL ref_to_tgt_cub ( ji_ro, jj_ro, gdept, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2833 END IF 2834 END IF 2443 2835 2444 2836 ELSE ! no subtraction of reference profile … … 2448 2840 END_3D 2449 2841 2450 END IF ! ln_hpg_ ffr_ref ) THEN2842 END IF ! ln_hpg_ref 2451 2843 2452 2844 DO_3D (0,0,0,0,1,jpk)
Note: See TracChangeset
for help on using the changeset viewer.
