Changeset 15726 for NEMO/branches/UKMO
- Timestamp:
- 2022-02-23T16:54:01+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
r14834 r15726 28 28 !! hpg_sco : s-coordinate (standard jacobian formulation) 29 29 !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf 30 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial)30 !! hpg_djc : s-coordinate (Density Jacobian with constrained cubic splines (ccs)) 31 31 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 32 !! hpg_djr : s-coordinate (Density Jacobian with ccs subtracting a reference) 33 !! hpg_ffr : s-coordinate (Forces on faces subtracting a reference profile) 32 34 !!---------------------------------------------------------------------- 33 35 USE oce ! ocean dynamics and tracers … … 63 65 LOGICAL, PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 64 66 LOGICAL, PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 67 LOGICAL, PUBLIC :: ln_hpg_djr !: s-coordinate (density Jacobian with cubic polynomial, subtracting a local reference profile) 68 LOGICAL, PUBLIC :: ln_hpg_ffr !: s-coordinate (forces on faces with subtraction of a local reference profile) 65 69 66 70 ! !! Flag to control the type of hydrostatic pressure gradient … … 72 76 INTEGER, PARAMETER :: np_prj = 4 ! s-coordinate (Pressure Jacobian scheme) 73 77 INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf 78 INTEGER, PARAMETER :: np_djr = 6 ! s-coordinate (density Jacobian with cubic polynomial, subtracting a local reference profile) 79 INTEGER, PARAMETER :: np_ffr = 7 ! s-coordinate (forces on faces with subtraction of a local reference profile) 80 74 81 ! 75 82 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 76 83 ! 77 LOGICAL :: ln_hpg_djc_vnh, ln_hpg_djc_vnv ! flag to specify hpg_djc boundary condition type 78 REAL(wp), PUBLIC :: aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions 84 85 LOGICAL :: ln_hpg_bcvN_rhd_hor, ln_hpg_bcvN_rhd_srf ! flags to specify constrained cubic spline (ccs) bdy conditions for rhd & z 86 LOGICAL :: ln_hpg_bcvN_rhd_bot, ln_hpg_bcvN_z_hor ! True implies von Neumann bcs; False implies linear extrapolation 87 88 REAL(wp) :: aco_bc_rhd_hor, bco_bc_rhd_hor, aco_bc_rhd_srf ! coefficients for hpg_djc & hpg_djr boundary conditions 89 REAL(wp) :: bco_bc_rhd_srf, aco_bc_rhd_bot, bco_bc_rhd_bot ! " 90 REAL(wp) :: aco_bc_z_hor, bco_bc_z_hor, aco_bc_z_srf ! " 91 REAL(wp) :: bco_bc_z_srf, aco_bc_z_bot, bco_bc_z_bot ! " 92 93 LOGICAL :: ln_hpg_djr_ref_ccs ! T => constrained cubic, F => simple cubic used for interpolation of reference to target profiles 94 LOGICAL :: ln_hpg_ffr_ref ! T => reference profile is subtracted; F => no reference profile subtracted 95 LOGICAL :: ln_hpg_ffr_ref_ccs ! T => use constrained cubic spline to vertically interpolate reference to each profile 96 LOGICAL :: ln_hpg_ffr_hor_ccs ! T => use constrained cubic spline to interpolate z_rhd_pmr along upper faces of cells 97 LOGICAL :: ln_hpg_ffr_hor_cub ! T => use simple cubic polynomial to interpolate z_rhd_pmr along upper faces of cells 98 LOGICAL :: ln_hpg_ffr_vrt_quad ! T => use quadratic fit to z_rhd_pmr in vertical interpoln & integration; else standard 2nd order scheme 99 100 LOGICAL :: ln_dbg_hpg ! T => debug outputs generated 101 LOGICAL :: ln_dbg_ij, ln_dbg_ik, ln_dbg_jk 102 INTEGER :: ki_dbg_min, ki_dbg_max, ki_dbg_cen 103 INTEGER :: kj_dbg_min, kj_dbg_max, kj_dbg_cen 104 INTEGER :: kk_dbg_min, kk_dbg_max, kk_dbg_cen 105 79 106 80 107 !! * Substitutions … … 121 148 CASE ( np_prj ) ; CALL hpg_prj ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Pressure Jacobian scheme) 122 149 CASE ( np_isf ) ; CALL hpg_isf ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate similar to sco modify for ice shelf 150 CASE ( np_djr ) ; CALL hpg_djr ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (Density Jacobian with Cubic polynomial subtracting a local reference profile) 151 CASE ( np_ffr ) ; CALL hpg_ffr ( kt, Kmm, puu, pvv, Krhs ) ! s-coordinate (forces on faces with subtraction of a local reference profile) 123 152 END SELECT 124 153 ! … … 158 187 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 159 188 !! 160 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 161 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, & 162 & ln_hpg_djc_vnh, ln_hpg_djc_vnv 189 190 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_isf, & 191 & ln_hpg_djc, ln_hpg_prj, ln_hpg_djr, ln_hpg_ffr, & 192 & ln_hpg_bcvN_rhd_hor, ln_hpg_bcvN_rhd_srf, & 193 & ln_hpg_bcvN_rhd_bot, ln_hpg_bcvN_z_hor, & 194 & ln_hpg_djr_ref_ccs, ln_hpg_ffr_vrt_quad, & 195 & ln_hpg_ffr_ref, ln_hpg_ffr_ref_ccs, & 196 & ln_hpg_ffr_hor_ccs, ln_hpg_ffr_hor_cub, & 197 & ln_dbg_hpg, ln_dbg_ij, ln_dbg_ik, ln_dbg_jk, & 198 & ki_dbg_min, ki_dbg_max, ki_dbg_cen, & 199 & kj_dbg_min, kj_dbg_max, kj_dbg_cen, & 200 & kk_dbg_min, kk_dbg_max, kk_dbg_cen 201 163 202 !!---------------------------------------------------------------------- 164 203 ! … … 175 214 WRITE(numout,*) '~~~~~~~~~~~~' 176 215 WRITE(numout,*) ' Namelist namdyn_hpg : choice of hpg scheme' 177 WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco 178 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 179 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 180 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 181 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 182 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 216 WRITE(numout,*) ' z-coord. - full steps ln_hpg_zco = ', ln_hpg_zco 217 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 218 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 219 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 220 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 221 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj 222 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic minus reference) ln_hpg_djr = ', ln_hpg_djr 223 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic minus reference) ln_hpg_ffr = ', ln_hpg_djr 224 WRITE(numout,*) ' s-coord. (forces on faces minus local reference) ln_hpg_ffr = ', ln_hpg_ffr 183 225 ENDIF 184 226 ! … … 188 230 IF( (.NOT.ln_hpg_isf .AND. ln_isfcav) .OR. (ln_hpg_isf .AND. .NOT.ln_isfcav) ) & 189 231 & CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' ) 232 190 233 ! 191 234 #if defined key_qco … … 194 237 ENDIF 195 238 #endif 239 240 IF( ln_hpg_djr .AND. nn_hls < 2) THEN 241 CALL ctl_stop( 'dyn_hpg_init : nn_hls < 2 and ln_hpg_djr not yet compatible' ) 242 ENDIF 243 196 244 ! 197 245 ! ! Set nhpg from ln_hpg_... flags & consistency check … … 204 252 IF( ln_hpg_prj ) THEN ; nhpg = np_prj ; ioptio = ioptio +1 ; ENDIF 205 253 IF( ln_hpg_isf ) THEN ; nhpg = np_isf ; ioptio = ioptio +1 ; ENDIF 254 IF( ln_hpg_djr ) THEN ; nhpg = np_djr ; ioptio = ioptio +1 ; ENDIF 255 IF( ln_hpg_ffr ) THEN ; nhpg = np_ffr ; ioptio = ioptio +1 ; ENDIF 206 256 ! 207 257 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) … … 216 266 CASE( np_prj ) ; WRITE(numout,*) ' ==>>> s-coord. (Pressure Jacobian: Cubic polynomial)' 217 267 CASE( np_isf ) ; WRITE(numout,*) ' ==>>> s-coord. (standard jacobian formulation) for isf' 268 CASE( np_djr ) ; WRITE(numout,*) ' ==>>> s-coord. (Density Jacobian: Cubic polynomial subtracting a local reference profile)' 269 CASE( np_ffr ) ; WRITE(numout,*) ' ==>>> s-coord. (forces on faces subtracting a local reference profile)' 270 218 271 END SELECT 219 272 WRITE(numout,*) 220 273 ENDIF 221 274 ! 222 IF ( ln_hpg_djc ) THEN223 IF (ln_hpg_ djc_vnh) THEN ! Von Neumann boundary condition224 IF(lwp) WRITE(numout,*) ' horizontal bc: von Neumann '225 aco_bc_ hor = 6.0_wp/5.0_wp226 bco_bc_ hor = 7.0_wp/15.0_wp275 IF ( ln_hpg_djc .OR. ln_hpg_djr .OR. ln_hpg_ffr ) THEN 276 IF (ln_hpg_bcvN_rhd_hor) THEN ! Von Neumann boundary condition 277 IF(lwp) WRITE(numout,*) ' ccs rhd horizontal bc: von Neumann ' 278 aco_bc_rhd_hor = 6.0_wp/5.0_wp 279 bco_bc_rhd_hor = 7.0_wp/15.0_wp 227 280 ELSE ! Linear extrapolation 228 IF(lwp) WRITE(numout,*) ' horizontal bc: linear extrapolation'229 aco_bc_ hor = 3.0_wp/2.0_wp230 bco_bc_ hor = 1.0_wp/2.0_wp281 IF(lwp) WRITE(numout,*) ' ccs rhd horizontal bc: linear extrapolation' 282 aco_bc_rhd_hor = 3.0_wp/2.0_wp 283 bco_bc_rhd_hor = 1.0_wp/2.0_wp 231 284 END IF 232 IF (ln_hpg_ djc_vnv) THEN ! Von Neumann boundary condition233 IF(lwp) WRITE(numout,*) ' verticalbc: von Neumann '234 aco_bc_ vrt= 6.0_wp/5.0_wp235 bco_bc_ vrt= 7.0_wp/15.0_wp285 IF (ln_hpg_bcvN_rhd_srf) THEN ! Von Neumann boundary condition 286 IF(lwp) WRITE(numout,*) ' ccs rhd surface bc: von Neumann ' 287 aco_bc_rhd_srf = 6.0_wp/5.0_wp 288 bco_bc_rhd_srf = 7.0_wp/15.0_wp 236 289 ELSE ! Linear extrapolation 237 IF(lwp) WRITE(numout,*) ' verticalbc: linear extrapolation'238 aco_bc_ vrt= 3.0_wp/2.0_wp239 bco_bc_ vrt= 1.0_wp/2.0_wp290 IF(lwp) WRITE(numout,*) ' ccs rhd surface bc: linear extrapolation' 291 aco_bc_rhd_srf = 3.0_wp/2.0_wp 292 bco_bc_rhd_srf = 1.0_wp/2.0_wp 240 293 END IF 294 IF (ln_hpg_bcvN_rhd_bot) THEN ! Von Neumann boundary condition 295 IF(lwp) WRITE(numout,*) ' ccs rhd bottom bc: von Neumann ' 296 aco_bc_rhd_bot = 6.0_wp/5.0_wp 297 bco_bc_rhd_bot = 7.0_wp/15.0_wp 298 ELSE ! Linear extrapolation 299 IF(lwp) WRITE(numout,*) ' ccs rhd bottom bc: linear extrapolation' 300 aco_bc_rhd_bot = 3.0_wp/2.0_wp 301 bco_bc_rhd_bot = 1.0_wp/2.0_wp 302 END IF 303 IF (ln_hpg_bcvN_z_hor) THEN ! Von Neumann boundary condition 304 IF(lwp) WRITE(numout,*) ' ccs z horizontal bc: von Neumann ' 305 aco_bc_z_hor = 6.0_wp/5.0_wp 306 bco_bc_z_hor = 7.0_wp/15.0_wp 307 ELSE ! Linear extrapolation 308 IF(lwp) WRITE(numout,*) ' ccs z horizontal bc: linear extrapolation' 309 aco_bc_z_hor = 3.0_wp/2.0_wp 310 bco_bc_z_hor = 1.0_wp/2.0_wp 311 END IF 312 313 ! linear extrapolation used for z in the vertical (surface and bottom) 314 aco_bc_z_srf = 3.0_wp/2.0_wp 315 bco_bc_z_srf = 1.0_wp/2.0_wp 316 aco_bc_z_bot = 3.0_wp/2.0_wp 317 bco_bc_z_bot = 1.0_wp/2.0_wp 318 241 319 END IF 242 ! 320 321 IF ( lwp .AND. ln_hpg_djr ) THEN 322 IF ( ln_hpg_djr_ref_ccs ) THEN 323 WRITE(numout,*) ' interpolation of reference profile by constrained cubic spline' 324 ELSE 325 WRITE(numout,*) ' interpolation of reference profile by simple cubic spline with off-centring at boundaries' 326 END IF 327 END IF 328 329 IF ( lwp .AND. ln_hpg_ffr ) THEN 330 IF ( ln_hpg_ffr_ref_ccs ) THEN 331 IF ( ln_hpg_ffr_ref_ccs ) THEN 332 WRITE(numout,*) ' interpolation of reference profile in vertical by constrained cubic spline ' 333 ELSE 334 WRITE(numout,*) ' interpolation of reference profile in vertical by pure cubic polynomial fit ' 335 END IF 336 ELSE 337 WRITE(numout,*) ' no reference profile subtracted ' 338 END IF 339 IF ( ln_hpg_ffr_hor_ccs ) THEN 340 WRITE(numout,*) ' interpolation of z_rhd_pmr profile in horizontal by constrained cubic spline ' 341 ELSE if ( ln_hpg_ffr_hor_cub ) THEN 342 WRITE(numout,*) ' interpolation of z_rhd_pmr profile in horizontal by simple cubic polynomial ' 343 ELSE 344 WRITE(numout,*) ' simple linear interpolation in horizontal of z_rhd_pmr profile' 345 END IF 346 IF ( ln_hpg_ffr_vrt_quad ) THEN 347 WRITE(numout,*) ' quadratic fit to z_rhd_pmr profile in vertical for interpolation and integration ' 348 ELSE 349 WRITE(numout,*) ' simple second order accurate vertical integration of z_rhd_pmr profile in vertical' 350 END IF 351 352 END IF 353 354 ! 355 IF ( ln_dbg_hpg .AND. lwp ) THEN 356 WRITE(numout,*) ' dyn_hpg diagnostic settings' 357 WRITE(numout,*) ' ki_dbg_min = ', ki_dbg_min, '; ki_dbg_max = ', ki_dbg_max 358 WRITE(numout,*) ' kj_dbg_min = ', kj_dbg_min, '; kj_dbg_max = ', kj_dbg_max 359 WRITE(numout,*) ' kk_dbg_min = ', kk_dbg_min, '; kk_dbg_max = ', kk_dbg_max 360 WRITE(numout,*) ' ki_dbg_cen = ', ki_dbg_cen, '; kj_dbg_cen = ', kj_dbg_cen 361 WRITE(numout,*) ' kk_dbg_cen = ', kk_dbg_cen 362 END IF 363 243 364 END SUBROUTINE dyn_hpg_init 244 365 … … 741 862 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 742 863 DO_2D( 1, 1, 1, 1 ) 743 zdrho_k(ji,jj,1) = aco_bc_ vrt * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_vrt* zdrho_k(ji,jj,2)744 zdz_k (ji,jj,1) = aco_bc_ vrt * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_vrt* zdz_k (ji,jj,2)864 zdrho_k(ji,jj,1) = aco_bc_rhd_srf * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_rhd_srf * zdrho_k(ji,jj,2) 865 zdz_k (ji,jj,1) = aco_bc_z_srf * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_z_srf * zdz_k (ji,jj,2) 745 866 END_2D 746 867 … … 748 869 IF ( mbkt(ji,jj)>1 ) THEN 749 870 iktb = mbkt(ji,jj) 750 zdrho_k(ji,jj,iktb) = aco_bc_ vrt * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1)751 zdz_k (ji,jj,iktb) = aco_bc_ vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt* zdz_k (ji,jj,iktb-1)871 zdrho_k(ji,jj,iktb) = aco_bc_rhd_bot * ( rhd(ji,jj,iktb) - rhd(ji,jj,iktb-1) ) - bco_bc_rhd_bot * zdrho_k(ji,jj,iktb-1) 872 zdz_k (ji,jj,iktb) = aco_bc_z_bot * ( -gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_z_bot * zdz_k (ji,jj,iktb-1) 752 873 END IF 753 874 END_2D 875 876 IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. zdz_k', zdz_k ) 877 IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. zdrho_k', zdrho_k ) 754 878 755 879 !-------------------------------------------------------------- … … 787 911 END_3D 788 912 913 IF ( ln_dbg_hpg ) CALL dbg_3dr( '4. z_rho_k', z_rho_k ) 914 789 915 !---------------------------------------------------------------------------------------- 790 916 ! 5. compute and store elementary horizontal differences in provisional arrays … … 840 966 DO_2D( 0, 0, 0, 1 ) 841 967 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 842 zz_drho_i(ji,jj) = aco_bc_ hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk)843 zz_dz_i (ji,jj) = aco_bc_ hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor* zdz_i (ji+1,jj,jk)968 zz_drho_i(ji,jj) = aco_bc_rhd_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_rhd_hor * zdrho_i(ji+1,jj,jk) 969 zz_dz_i (ji,jj) = aco_bc_z_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zdz_i (ji+1,jj,jk) 844 970 END IF 845 971 END_2D … … 847 973 DO_2D( -1, 1, 0, 1 ) 848 974 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 849 zz_drho_i(ji,jj) = aco_bc_ hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)850 zz_dz_i (ji,jj) = aco_bc_ hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor* zdz_i (ji-1,jj,jk)975 zz_drho_i(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_rhd_hor * zdrho_i(ji-1,jj,jk) 976 zz_dz_i (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_z_hor * zdz_i (ji-1,jj,jk) 851 977 END IF 852 978 END_2D … … 854 980 DO_2D( 0, 1, 0, 0 ) 855 981 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 856 zz_drho_j(ji,jj) = aco_bc_ hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)857 zz_dz_j (ji,jj) = aco_bc_ hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor* zdz_j (ji,jj+1,jk)982 zz_drho_j(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_rhd_hor * zdrho_j(ji,jj+1,jk) 983 zz_dz_j (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zdz_j (ji,jj+1,jk) 858 984 END IF 859 985 END_2D … … 861 987 DO_2D( 0, 1, -1, 1 ) 862 988 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 863 zz_drho_j(ji,jj) = aco_bc_ hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)864 zz_dz_j (ji,jj) = aco_bc_ hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor* zdz_j (ji,jj-1,jk)989 zz_drho_j(ji,jj) = aco_bc_rhd_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_rhd_hor * zdrho_j(ji,jj-1,jk) 990 zz_dz_j (ji,jj) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_z_hor * zdz_j (ji,jj-1,jk) 865 991 END IF 866 992 END_2D … … 869 995 zdrho_j(:,:,jk) = zz_drho_j(:,:) 870 996 zdz_j (:,:,jk) = zz_dz_j (:,:) 871 END DO 997 END DO ! jk 998 999 IF ( ln_dbg_hpg ) THEN 1000 CALL dbg_3dr( '6. zdrho_i', zdrho_i ) 1001 CALL dbg_3dr( '6. zdrho_j', zdrho_j ) 1002 END IF 872 1003 873 1004 !-------------------------------------------------------------- 874 ! 7. Calculate integrals on sidefaces1005 ! 7. Calculate integrals on upper/lower faces 875 1006 !------------------------------------------------------------- 876 1007 … … 900 1031 END_3D 901 1032 1033 IF ( ln_dbg_hpg ) THEN 1034 CALL dbg_3dr( '7. z_rho_i', z_rho_i ) 1035 CALL dbg_3dr( '7. z_rho_j', z_rho_j ) 1036 END IF 1037 902 1038 !-------------------------------------------------------------- 903 1039 ! 8. Integrate in the vertical … … 939 1075 END_3D 940 1076 ! 1077 1078 IF ( ln_dbg_hpg ) THEN 1079 CALL dbg_3dr( '8. zhpi', zhpi ) 1080 CALL dbg_3dr( '8. zhpj', zhpj ) 1081 END IF 1082 941 1083 IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy ) 942 1084 ! … … 1435 1577 END FUNCTION integ_spline 1436 1578 1579 SUBROUTINE hpg_djr( kt, Kmm, puu, pvv, Krhs ) 1580 !!--------------------------------------------------------------------- 1581 !! *** ROUTINE hpg_djr *** 1582 !! 1583 !! ** Method : Density Jacobian with Cubic polynomial scheme subtracting a local reference profile (pmr is profile minus reference) 1584 !! This code assumes a 2-point halo 1585 !! 1586 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 1587 !!---------------------------------------------------------------------- 1588 INTEGER , INTENT( in ) :: kt ! ocean time-step index 1589 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 1590 REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 1591 !! 1592 INTEGER :: ji, jj, jk, jr ! loop indices 1593 INTEGER :: jn_hor_pts ! number of points in the horizontal stencil 1594 INTEGER :: j_uv ! 1 for u-cell; 2 for v-cell 1595 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 1596 INTEGER :: jia, jib, jja, jjb ! 1597 INTEGER :: jir, jjr ! reference (expand) 1598 INTEGER :: jio, jjo ! offset (expand) 1599 1600 REAL(wp) :: z_grav_10, z1_12 ! constants 1601 REAL(wp) :: zhta, zhtb ! temporary scalars 1602 REAL(wp) :: zcoef0, zep, cffw ! " " 1603 REAL(wp) :: aco, bco ! " " 1604 REAL(wp) :: cffu, cffx, z1_cff ! " " 1605 REAL(wp) :: cffv, cffy ! " " 1606 REAL(wp) :: cff_31, cff_42 ! " " 1607 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 1608 1609 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztmp, zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences 1610 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrhd_ref ! Reference density 1611 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zz_ref ! Reference heights 1612 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhd_k_ref ! Harmonic average of primitive differences for reference field 1613 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_rhd_pmr ! rhd_prm = rhd - rhd_ref (values on the original grid) 1614 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhd_k_pmr ! 1615 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_lmr_k ! left minus right density integrals on vertical faces 1616 1617 INTEGER, DIMENSION(A2D(nn_hls)) :: jk_bot_ref ! bottom levels in the reference profile 1618 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdzx, zdzy ! primitive differences in x and y 1619 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j 1620 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdrhd_21, zdrhd_32, zdrhd_43 1621 REAL(wp), DIMENSION(A2D(nn_hls),2) :: zz_drhd_ij, zdrhd_ij 1622 REAL(wp), DIMENSION(A2D(nn_hls)) :: z_low_ij, z_upp_ij 1623 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj 1624 1625 !!---------------------------------------------------------------------- 1626 1627 IF( kt == nit000 ) THEN 1628 IF(lwp) WRITE(numout,*) 1629 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 1630 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 1631 ENDIF 1632 1633 ! Local constant initialization 1634 zcoef0 = - grav * 0.5_wp 1635 z_grav_10 = grav / 10._wp 1636 z1_12 = 1.0_wp / 12._wp 1637 zep = 1.e-15 1638 1639 jn_hor_pts = 4 ! 4 points in the horizontal stencil 1640 1641 !------------------------------------------------------------------------------------------------------ 1642 ! 1. calculate harmonic averages of differences for z (grid heights) in i, j and k directions 1643 !------------------------------------------------------------------------------------------------------ 1644 1645 !------------------------------------------------------------------------------------------------------ 1646 ! 1.1 compute and store elementary vertical differences then harmonic averages for z using eqn 5.18 1647 ! Full domain covered so that _ref profiles can be taken from zdz_k 1648 !------------------------------------------------------------------------------------------------------ 1649 1650 DO_3D( 2, 2, 2, 2, 2, jpk ) 1651 ztmp (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1) 1652 END_3D 1653 1654 zdz_k (:,:,1) = 0._wp ! jk index changed from : to 1 to make computationally less wasteful 1655 1656 DO_3D( 2, 2, 2, 2, 2, jpk-1 ) 1657 zdz_k(ji,jj,jk) = 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1) & 1658 & / ( ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1) ) 1659 END_3D 1660 1661 !---------------------------------------------------------------------------------- 1662 ! 1.2 apply boundary conditions at top and bottom using 5.36-5.37 1663 !---------------------------------------------------------------------------------- 1664 1665 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 1666 zdz_k (:,:,1) = aco_bc_z_srf * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_z_srf * zdz_k (:,:,2) 1667 1668 DO_2D( 2, 2, 2, 2 ) 1669 iktb = mbkt(ji,jj) 1670 IF ( iktb > 1 ) THEN 1671 zdz_k (ji,jj,iktb) = aco_bc_z_bot * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_z_bot * zdz_k (ji,jj,iktb-1) 1672 END IF 1673 END_2D 1674 1675 !---------------------------------------------------------------------------------------- 1676 ! 1.3 compute and store elementary horizontal differences then harmonic averages for z using eqn 5.18 1677 !---------------------------------------------------------------------------------------- 1678 1679 DO jk = 1, jpkm1 1680 DO_2D( 1, 1, 1, 1 ) 1681 zdzx (ji,jj) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 1682 zdzy (ji,jj) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 1683 END_2D 1684 1685 ! this subroutine requires a 2-point halo CALL lbc_lnk_multi( 'dynhpg', zdzx, 'U', 1., zdzy, 'V', 1. ) 1686 1687 DO_2D( 0, 1, 0, 1 ) 1688 cffx = MAX( 2._wp * zdzx (ji-1,jj) * zdzx (ji,jj), 0._wp) 1689 z1_cff = zdzx(ji-1,jj) + zdzx(ji,jj) 1690 zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 1691 1692 cffy = 2._wp * zdzy (ji ,jj-1) * zdzy (ji,jj ) 1693 z1_cff = zdzy(ji,jj-1) + zdzy(ji,jj) 1694 zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 1695 END_2D 1696 1697 !---------------------------------------------------------------------------------- 1698 ! 1.4 apply boundary conditions at sides using 5.36-5.37 1699 !---------------------------------------------------------------------------------- 1700 1701 DO_2D( 0, 1, 0, 1) 1702 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 1703 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 1704 zdz_i(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zz_dz_i(ji+1,jj) 1705 END IF 1706 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 1707 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 1708 zdz_i(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_z_hor * zz_dz_i(ji-1,jj) 1709 END IF 1710 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 1711 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 1712 zdz_j(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_z_hor * zz_dz_j(ji,jj+1) 1713 END IF 1714 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 1715 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 1716 zdz_j(ji,jj,jk) = aco_bc_z_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_z_hor * zz_dz_j(ji,jj-1) 1717 END IF 1718 END_2D 1719 END DO ! k 1720 1721 IF ( ln_dbg_hpg ) THEN 1722 CALL dbg_3dr( '1.4 gde3w', gde3w ) 1723 CALL dbg_3dr( '1.4 zdz_i', zdz_i ) 1724 CALL dbg_3dr( '1.4 zdz_j', zdz_j ) 1725 CALL dbg_3dr( '1.4 zdz_k', zdz_k ) 1726 CALL dbg_2di( 'mbkt', mbkt) 1727 END IF 1728 !---------------------------------------------------------------------------------------- 1729 ! 2. Start loop over the u and v components and find the reference profile 1730 ! The loop ends in section 5.4 1731 !---------------------------------------------------------------------------------------- 1732 1733 DO j_uv = 1, 2 ! j_uv = 1 is for u-cell ; j_uv = 2 for v-cell 1734 1735 !---------------------------------------------------------------------------------------- 1736 ! 2.1 find reference profiles zrhd_ref and zz_ref and the bottom level of the reference profile 1737 !---------------------------------------------------------------------------------------- 1738 1739 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) 1741 1742 IF ( ln_dbg_hpg ) THEN 1743 CALL dbg_3dr( '2.1 rhd', rhd ) 1744 CALL dbg_3dr( '2.1 zrhd_ref', zrhd_ref ) 1745 CALL dbg_3dr( '2.1 zz_ref', zz_ref ) 1746 CALL dbg_2di( '2.1 jk_bot_ref', jk_bot_ref ) 1747 END IF 1748 1749 !-------------------------------------------------------------------------------------------------------- 1750 ! 2.2 IF ln_hpg_djr_ref_ccs compute zdrhd_k_ref then set bcs at top & bottom 1751 ! (bcs not needed for simple cubic off-centred at boundaries) 1752 !-------------------------------------------------------------------------------------------------------- 1753 1754 IF ( ln_hpg_djr_ref_ccs ) THEN 1755 CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 1756 IF ( ln_dbg_hpg ) CALL dbg_3dr( '2.3 zdrhd_k_ref', zdrhd_k_ref ) 1757 END IF ! ln_hpg_djr_ref_ccs 1758 1759 !---------------------------------------------------------------------------------------- 1760 ! 3. interpolate reference profiles to target profiles and form difference profiles z_rhd_pmr 1761 !---------------------------------------------------------------------------------------- 1762 1763 DO jr = 1, 4 1764 IF ( j_uv == 1 ) THEN 1765 jio = jr - 2 ! range of jio is -1 to 2 1766 jjo = 0 1767 ELSE 1768 jio = 0 1769 jjo = jr - 2 1770 END IF 1771 1772 IF ( ln_hpg_djr_ref_ccs ) THEN 1773 CALL ref_to_tgt_ccs ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1774 ELSE 1775 CALL ref_to_tgt_cub ( jio, jjo, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 1776 END IF 1777 1778 IF ( ln_dbg_hpg ) CALL dbg_3dr( '3. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) ) 1779 1780 END DO 1781 1782 !---------------------------------------------------------------------------------------- 1783 ! 4. Calculations for side-face integrals 1784 !---------------------------------------------------------------------------------------- 1785 1786 !---------------------------------------------------------------------------------------- 1787 ! 4.1 compute and store elementary vertical differences then harmonic averages 1788 ! based on z_rhd_pmr arrays (zdz_k has already been calculated) 1789 !---------------------------------------------------------------------------------------- 1790 1791 ! start loop over the two side faces jr = 2 "left" face; jr = 3 "right" face 1792 DO jr = 2, 3 1793 1794 IF ( j_uv == 1 ) THEN 1795 jio = jr - 2 1796 jjo = 0 1797 ELSE 1798 jio = 0 1799 jjo = jr - 2 1800 END IF 1801 1802 DO_3D( 0, 0, 0, 0, 2, jpk ) 1803 ztmp(ji,jj,jk) = z_rhd_pmr (ji,jj,jk,jr) - z_rhd_pmr(ji,jj,jk-1,jr) 1804 END_3D 1805 1806 zdrhd_k_pmr(:,:,:) = 0._wp ! should be unnecessary 1807 1808 DO_3D( 0, 0, 0, 0, 2, jpk-1 ) 1809 cffw = MAX( 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1), 0._wp ) 1810 z1_cff = ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1) 1811 zdrhd_k_pmr(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 1812 END_3D 1813 1814 ! apply boundary conditions at top and bottom 1815 DO_2D( 0, 0, 0, 0 ) 1816 zdrhd_k_pmr(ji,jj,1) = aco_bc_rhd_srf * ( z_rhd_pmr(ji,jj,2,jr) - z_rhd_pmr(ji,jj,1,jr) ) - bco_bc_rhd_srf * zdrhd_k_pmr(ji,jj,2) 1817 iktb = mbkt(ji+jio,jj+jjo) 1818 IF ( iktb > 1 ) THEN 1819 zdrhd_k_pmr(ji,jj,iktb) = aco_bc_rhd_bot * (z_rhd_pmr(ji,jj,iktb,jr) - z_rhd_pmr(ji,jj,iktb-1,jr) ) - bco_bc_rhd_bot * zdrhd_k_pmr(ji,jj,iktb-1) 1820 END IF 1821 END_2D 1822 1823 1824 !-------------------------------------------------------------- 1825 ! 4.2 Upper half of top-most grid box, compute and store 1826 !------------------------------------------------------------- 1827 !! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights 1828 !! rho_k stores grav * FX / rho_0 1829 !! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1,Kmm) 1830 DO_2D( 0, 0, 0, 0) 1831 ztmp(ji,jj,1) = grav * ( ssh(ji+jio,jj+jjo,Kmm) + gde3w(ji+jio,jj+jjo,1) ) & 1832 & * ( z_rhd_pmr(ji,jj,1,jr) & 1833 & + 0.5_wp * ( z_rhd_pmr(ji,jj,2,jr) - z_rhd_pmr(ji,jj,1,jr) ) & 1834 & * ( ssh (ji+jio,jj+jjo,Kmm) + gde3w(ji+jio,jj+jjo,1) ) & 1835 & / ( - gde3w(ji+jio,jj+jjo,2) + gde3w(ji+jio,jj+jjo,1) ) ) 1836 END_2D 1837 1838 !-------------------------------------------------------------- 1839 ! 4.3 Interior faces, compute and store 1840 !------------------------------------------------------------- 1841 1842 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1843 ztmp(ji,jj,jk) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,jr) + z_rhd_pmr(ji,jj,jk-1,jr) ) & 1844 & * ( - gde3w(ji+jio,jj+jjo,jk) + gde3w(ji+jio,jj+jjo,jk-1) ) & 1845 & + z_grav_10 * ( & 1846 & ( zdrhd_k_pmr(ji,jj,jk) - zdrhd_k_pmr(ji,jj,jk-1) ) & 1847 & * ( - gde3w(ji+jio,jj+jjo,jk) + gde3w(ji+jio,jj+jjo,jk-1) - z1_12 * ( zdz_k(ji+jio,jj+jjo,jk) + zdz_k (ji+jio,jj+jjo,jk-1) ) ) & 1848 & - ( zdz_k(ji+jio,jj+jjo,jk) - zdz_k(ji+jio,jj+jjo,jk-1) ) & 1849 & * ( z_rhd_pmr(ji,jj,jk,jr) - z_rhd_pmr(ji,jj,jk-1,jr) - z1_12 * ( zdrhd_k_pmr(ji,jj,jk) + zdrhd_k_pmr(ji,jj,jk-1) ) ) & 1850 & ) 1851 END_3D 1852 1853 ! the force on the right face could be set equal to the average of the right face for this cell and the left face for the cell to the right 1854 ! this would require an lbc_lnk call 1855 1856 ! lmr stands for left minus right 1857 1858 IF ( jr == 2 ) THEN 1859 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1860 z_lmr_k(ji,jj,jk) = ztmp(ji,jj,jk) ! values on left face; 1861 END_3D 1862 ELSE 1863 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1864 z_lmr_k(ji,jj,jk) = z_lmr_k(ji,jj,jk) - ztmp(ji,jj,jk) ! subtract the values on the right face 1865 END_3D 1866 END IF 1867 1868 IF ( ln_dbg_hpg ) CALL dbg_3dr( '4. z_lmr_k', z_lmr_k ) 1869 1870 END DO ! jr 1871 1872 1873 !---------------------------------------------------------------------------------------- 1874 ! 5. Calculations for upper and lower faces and the vertical integration 1875 !---------------------------------------------------------------------------------------- 1876 1877 z_upp_ij(:,:) = 0._wp 1878 zhpi(:,:) = 0._wp 1879 zhpj(:,:) = 0._wp 1880 1881 DO jk = 1, jpk -1 1882 1883 IF ( ln_dbg_hpg .AND. lwp ) THEN 1884 WRITE(numout,*) 1885 WRITE(numout,*) ' jk = ', jk 1886 END IF 1887 1888 !---------------------------------------------------------------------------------------- 1889 ! 5.1 compute and store elementary horizontal differences zfor z_rhd_pmr arrays 1890 !---------------------------------------------------------------------------------------- 1891 1892 DO_2D( 0, 0, 0, 0 ) 1893 zdrhd_21(ji,jj) = z_rhd_pmr(ji,jj,jk,2) - z_rhd_pmr(ji,jj,jk,1) 1894 zdrhd_32(ji,jj) = z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) 1895 zdrhd_43(ji,jj) = z_rhd_pmr(ji,jj,jk,4) - z_rhd_pmr(ji,jj,jk,3) 1896 END_2D 1897 1898 DO_2D( 0, 0, 0, 0 ) 1899 cff_31 = MAX( 2._wp * zdrhd_21(ji,jj) * zdrhd_32(ji,jj), 0._wp ) 1900 z1_cff = zdrhd_21(ji,jj) + zdrhd_32(ji,jj) 1901 zz_drhd_ij(ji,jj,1) = cff_31 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 1902 1903 cff_42 = MAX( 2._wp * zdrhd_32(ji,jj) * zdrhd_43(ji,jj), 0._wp ) 1904 z1_cff = zdrhd_32(ji,jj) + zdrhd_43(ji,jj) 1905 zz_drhd_ij(ji,jj,2) = cff_42 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 1906 END_2D 1907 1908 !---------------------------------------------------------------------------------- 1909 ! 5.2 apply boundary conditions at side boundaries using 5.36-5.37 1910 !---------------------------------------------------------------------------------- 1911 1912 1913 ! need to check this sub-section more carefully 1914 1915 zdrhd_ij(:,:,:) = zz_drhd_ij(:,:,:) 1916 1917 IF ( j_uv == 1 ) THEN 1918 1919 DO_2D( 0, 0, 0, 0) 1920 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 1921 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 1922 zdrhd_ij(ji,jj,1) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,2) 1923 END IF 1924 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 1925 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 1926 zdrhd_ij(ji,jj,2) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,1) 1927 END IF 1928 END_2D 1929 1930 ELSE ! j_uv == 2 1931 1932 DO_2D( 0, 0, 0, 0) 1933 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 1934 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 1935 zdrhd_ij(ji,jj,1) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,2) 1936 END IF 1937 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 1938 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 1939 zdrhd_ij(ji,jj,2) = aco_bc_rhd_hor * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) ) - bco_bc_rhd_hor * zz_drhd_ij(ji,jj,1) 1940 END IF 1941 END_2D 1942 1943 END IF ! j_uv == 2 1944 1945 IF ( ln_dbg_hpg ) THEN 1946 CALL dbg_2dr( '5.2 zdrhd_ij(:,:,1)', zdrhd_ij(:,:,1) ) 1947 CALL dbg_2dr( '5.2 zdrhd_ij(:,:,2)', zdrhd_ij(:,:,2) ) 1948 END IF 1949 1950 !-------------------------------------------------------------- 1951 ! 5.3 Calculate integrals on lower faces 1952 !------------------------------------------------------------- 1953 1954 IF ( j_uv == 1 ) THEN 1955 1956 DO_2D( 0, 0, 0, 0 ) 1957 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 1958 z_low_ij(ji,jj) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,3) + z_rhd_pmr(ji,jj,jk,2) ) & 1959 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) 1960 1961 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 1962 z_low_ij(ji,jj) = z_low_ij(ji,jj) - z_grav_10 * ( & 1963 & ( zdrhd_ij(ji,jj,2) - zdrhd_ij(ji,jj,1) ) & 1964 & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i(ji+1,jj,jk) + zdz_i(ji,jj,jk) ) ) & 1965 & - ( zdz_i(ji+1,jj,jk) - zdz_i(ji,jj,jk) ) & 1966 & * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) - z1_12 * ( zdrhd_ij(ji,jj,2) + zdrhd_ij(ji,jj,1) ) ) & 1967 & ) 1968 END IF 1969 END_2D 1970 1971 IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.3 z_low_ij 1', z_low_ij ) 1972 1973 ELSE ! j_uv == 2 1974 1975 DO_2D( 0, 0, 0, 0 ) 1976 z_low_ij(ji,jj) = zcoef0 * ( z_rhd_pmr(ji,jj,jk,3) + z_rhd_pmr(ji,jj,jk,2) ) & 1977 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) 1978 1979 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 1980 z_low_ij(ji,jj) = z_low_ij(ji,jj) - z_grav_10 * ( & 1981 & ( zdrhd_ij(ji,jj,2) - zdrhd_ij(ji,jj,1) ) & 1982 & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j(ji,jj+1,jk) + zdz_j(ji,jj,jk) ) ) & 1983 & - ( zdz_j(ji,jj+1,jk) - zdz_j(ji,jj,jk) ) & 1984 & * ( z_rhd_pmr(ji,jj,jk,3) - z_rhd_pmr(ji,jj,jk,2) - z1_12 * ( zdrhd_ij(ji,jj,2) + zdrhd_ij(ji,jj,1) ) ) & 1985 & ) 1986 END IF 1987 END_2D 1988 1989 IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.3 z_low_ij 2', z_low_ij ) 1990 1991 END IF ! j_uv 1992 1993 !-------------------------------------------------------------- 1994 ! 5.4 Integrate in the vertical (including contributions from both upper and lower and side-faces) 1995 !------------------------------------------------------------- 1996 ! 1997 IF ( j_uv == 1 ) THEN 1998 1999 DO_2D( 0, 0, 0, 0 ) 2000 zhpi(ji,jj) = zhpi(ji,jj) + & 2001 & ( z_lmr_k(ji,jj,jk) - ( z_low_ij(ji,jj) - z_upp_ij(ji,jj) ) ) * r1_e1u(ji,jj) 2002 ! add to the general momentum trend 2003 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zhpi(ji,jj) 2004 END_2D 2005 2006 IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.4 zhpi', zhpi ) 2007 2008 ELSE ! j_uv == 2 2009 2010 DO_2D( 0, 0, 0, 0 ) 2011 zhpj(ji,jj) = zhpj(ji,jj) + & 2012 & ( z_lmr_k(ji,jj,jk) - ( z_low_ij(ji,jj) - z_upp_ij(ji,jj) ) ) * r1_e2v(ji,jj) 2013 ! add to the general momentum trend 2014 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zhpj(ji,jj) 2015 END_2D 2016 2017 IF ( ln_dbg_hpg ) CALL dbg_2dr( '5.4 zhpj', zhpj ) 2018 2019 END IF ! j_uv 2020 2021 DO_2D( 0, 0, 0, 0 ) 2022 z_upp_ij(ji,jj) = z_low_ij(ji,jj) 2023 END_2D 2024 2025 END DO ! k 2026 2027 END DO ! j_uv 2028 2029 END SUBROUTINE hpg_djr 2030 2031 !----------------------------------------------------------------------------------------------------------------- 2032 2033 SUBROUTINE ref_to_tgt_cub ( 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) 2034 2035 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2036 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2037 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dep_tgt ! depths of target profiles 2038 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_tgt ! field values on the target grid 2039 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_z_ref ! heights of reference profiles 2040 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_fld_ref ! field values to be interpolated (in the vertical) on reference grid 2041 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2042 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: p_fld_tgt_ref ! target minus reference on target grid 2043 2044 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. 2045 2046 !--------------------------------------------------------------------------------------------------------------- 2047 2048 INTEGER :: ji, jj, jk ! loop indices 2049 INTEGER :: jkr 2050 REAL(wp) :: z_r_6, z_r_24 ! constants 2051 2052 REAL(wp) :: zf_a, zf_b, zf_c, zf_d 2053 REAL(wp) :: zz_ref_jkr, zz_ref_jkrm1, zeta, zetasq 2054 REAL(wp) :: zf_0, zf_1, zf_2, zf_3 2055 REAL(wp) :: zave_bc, zave_ad, zdif_cb, zdif_da 2056 2057 REAL(wp) :: zz_tgt_lcl, zfld_ref_on_tgt 2058 2059 !----------------------------------------------------------------------------------------------------------------- 2060 2061 z_r_6 = 1._wp / 6._wp 2062 z_r_24 = 1._wp / 24._wp 2063 2064 ! 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 ) 2066 2067 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 2068 zz_tgt_lcl = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) 2069 2070 ! it would probably be better computationally for fld_ref to have the jk index first. 2071 2072 !!! jkr >= 2 and p_fld_ref has jk = 0 available 2073 2074 jkr = jk_ref_for_tgt(ji,jj,jk) 2075 zf_a = p_fld_ref(ji,jj,jkr-2) 2076 zf_b = p_fld_ref(ji,jj,jkr-1) 2077 zf_c = p_fld_ref(ji,jj,jkr ) 2078 zf_d = p_fld_ref(ji,jj,jkr+1) 2079 2080 zz_ref_jkrm1 = p_z_ref( ji, jj, jkr - 1 ) 2081 zz_ref_jkr = p_z_ref( ji, jj, jkr ) 2082 zeta = ( zz_tgt_lcl - 0.5_wp*(zz_ref_jkr+zz_ref_jkrm1) ) / ( zz_ref_jkr - zz_ref_jkrm1 ) 2083 zetasq = zeta*zeta 2084 2085 zave_bc = 0.5_wp*(zf_b+zf_c) 2086 zave_ad = 0.5_wp*(zf_a+zf_d) 2087 zdif_cb = zf_c - zf_b 2088 zdif_da = zf_d - zf_a 2089 2090 zf_0 = 1.125_wp*zave_bc - 0.125_wp*zave_ad 2091 zf_1 = 1.125_wp*zdif_cb - z_r_24*zdif_da 2092 zf_2 = 0.5_wp*(zave_ad - zave_bc) ! corrected 12/09/2021 2093 zf_3 = z_r_6 * zdif_da - 0.5_wp*zdif_cb ! corrected 12/09/2021 2094 2095 zfld_ref_on_tgt = zf_0 + zeta*zf_1 + zetasq*(zf_2 + zeta*zf_3) 2096 2097 ! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc. 2098 p_fld_tgt_ref(ji, jj, jk) = p_fld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt 2099 2100 END_3D 2101 2102 RETURN 2103 2104 END SUBROUTINE ref_to_tgt_cub 2105 2106 !----------------------------------------------------------------------------------------------------------------- 2107 2108 SUBROUTINE ref_to_tgt_ccs ( ki_off_tgt, kj_off_tgt, pdep_tgt, pfld_tgt, pz_ref, pfld_ref, pdfld_k_ref, kk_bot_ref, pfld_tgt_ref) 2109 2110 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2111 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2112 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdep_tgt ! depths of target profiles 2113 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_tgt ! field values on the target grid 2114 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pz_ref ! heights of reference profiles 2115 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pfld_ref ! reference field values 2116 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: pdfld_k_ref ! harmonic averages of vertical differences of reference field 2117 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2118 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: pfld_tgt_ref ! target minus reference on target grid 2119 2120 !----------------------------------------------------------------------------------------------------------------- 2121 INTEGER, DIMENSION(A2D(nn_hls),jpk) :: jk_ref_for_tgt ! reference index for interpolation to target grid 2122 2123 INTEGER :: ji, jj, jk ! DO loop indices 2124 INTEGER :: jkr 2125 REAL(wp) :: z_r_6 ! constant 2126 REAL(wp) :: zz_tgt_lcl, zz_ref_jkrm1, zz_ref_jkr, zeta, zetasq 2127 REAL(wp) :: zd_dif, zd_ave, zf_dif, zf_ave 2128 REAL(wp) :: zcoef_0, zcoef_1, zcoef_2, zcoef_3 2129 REAL(wp) :: zfld_ref_on_tgt 2130 !----------------------------------------------------------------------------------------------------------------- 2131 2132 z_r_6 = 1._wp / 6._wp 2133 2134 ! 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 ) 2136 2137 DO_3D( 0, 0, 0, 0, 1, jpk-1 ) 2138 zz_tgt_lcl = - pdep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk ) 2139 jkr = jk_ref_for_tgt( ji, jj, jk ) 2140 zz_ref_jkrm1 = pz_ref( ji, jj, jkr - 1 ) 2141 zz_ref_jkr = pz_ref( ji, jj, jkr ) 2142 zeta = ( zz_tgt_lcl - 0.5_wp*(zz_ref_jkr+zz_ref_jkrm1) ) / ( zz_ref_jkr - zz_ref_jkrm1 ) 2143 zetasq = zeta*zeta 2144 2145 zd_dif = pdfld_k_ref(ji,jj,jkr) - pdfld_k_ref(ji,jj,jkr-1) 2146 zd_ave = 0.5_wp * ( pdfld_k_ref(ji,jj,jkr) + pdfld_k_ref(ji,jj,jkr-1) ) 2147 2148 zf_dif = pfld_ref(ji,jj,jkr) - pfld_ref(ji,jj,jkr-1) 2149 zf_ave = 0.5_wp * ( pfld_ref(ji,jj,jkr) + pfld_ref(ji,jj,jkr-1) ) 2150 2151 zcoef_0 = zf_ave - 0.125_wp * zd_dif 2152 zcoef_1 = 1.5_wp * zf_dif - 0.5_wp * zd_ave 2153 zcoef_2 = 0.5_wp * zd_dif 2154 zcoef_3 = 2.0_wp * zd_ave - 2._wp * zf_dif 2155 2156 zfld_ref_on_tgt = zcoef_0 + zeta*zcoef_1 + zetasq * ( zcoef_2 + zeta*zcoef_3 ) 2157 2158 ! when zfld_ref_on_tgt is commented out in the next line, the results for hpg_djr should agree with those for hpg_djc. 2159 2160 pfld_tgt_ref(ji, jj, jk) = pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) - zfld_ref_on_tgt 2161 2162 ! IF ( ln_dbg_hpg .AND. lwp .AND. ji == ki_dbg_cen .AND. jj == kj_dbg_cen .AND. jk == kk_dbg_cen ) THEN 2163 ! WRITE(numout,*) ' zeta, zd_dif, zd_ave, zf_dif, zf_ave = ', zeta, zd_dif, zd_ave, zf_dif, zf_ave 2164 ! WRITE(numout,*) ' zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 =', zz_tgt_lcl, jkr, zz_ref_jkr, zz_ref_jkrm1 2165 ! WRITE(numout,*) ' pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) = ', pfld_ref(ji,jj,jkr), pfld_ref(ji,jj,jkr-1) 2166 ! WRITE(numout,*) ' zfld_ref_on_tgt = ', zfld_ref_on_tgt 2167 ! WRITE(numout,*) ' pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) = ', pfld_tgt(ji+ki_off_tgt, jj+kj_off_tgt, jk) 2168 ! END IF 2169 2170 END_3D 2171 2172 IF ( ln_dbg_hpg ) CALL dbg_3dr( 'ref_to_tgt_ccs: pfld_tgt_ref', pfld_tgt_ref ) 2173 2174 RETURN 2175 2176 END SUBROUTINE ref_to_tgt_ccs 2177 2178 !----------------------------------------------------------------------------------------------------------------- 2179 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 2182 INTEGER, INTENT(in) :: ki_off_tgt ! offset of points in target array in i-direction 2183 INTEGER, INTENT(in) :: kj_off_tgt ! offset of points in target array in j-direction 2184 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_dep_tgt ! depths of target profiles 2185 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_z_ref ! heights of reference profiles 2186 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot_ref ! bottom levels in the reference profile 2187 INTEGER, DIMENSION(A2D(nn_hls),jpk), INTENT(OUT) :: kk_ref_for_tgt ! reference index for interpolation to target grid (the lower point) 2188 ! with jkr = kk_ref_for_tgt(ji,jj,jk) the reference levels are jkr-1 and jkr 2189 !----------------------------------------------------------------------------------------------------------------- 2190 2191 INTEGER, DIMENSION(A2D(nn_hls)) :: jk_tgt, jk_ref ! vertical level being processed on target and reference grids respectively 2192 2193 INTEGER :: ji, jj, jk_comb, jk_bot 2194 2195 INTEGER :: jk, jkr, jcount ! for debugging only 2196 REAL(wp):: z_tgt, z_below, z_above ! for debugging only 2197 2198 INTEGER :: jk_init ! initial jk value for search 2199 2200 REAL(wp):: tol_wp ! this should be the working precision tolerance 2201 2202 tol_wp = 1.E-4 ! the inferred bottom depth seems to be accurate to about 1.E-6. 2203 2204 !----------------------------------------------------------------------------------------------------------------- 2205 2206 ! 1. Initialisation for search for points on reference grid bounding points on the target grid 2207 ! the first point on the target grid is assumed to lie between the first two points on the reference grid 2208 IF ( ln_hpg_djr_ref_ccs ) THEN 2209 jk_init = 2 2210 ELSE 2211 jk_init = 3 ! for simple cubic use off-centred interpolation near the surface 2212 END IF 2213 2214 jk_ref(:,:) = jk_init 2215 kk_ref_for_tgt(:,:,1) = jk_init 2216 jk_tgt(:,:) = 2 2217 2218 ! 2. Find kk_ref_for_tgt 2219 2220 DO jk_comb = 1, 2*(jpk-1) 2221 2222 ! if level jk_tgt in target profile is lower than jk_ref in reference profile add one to jk_ref 2223 DO_2D( 0, 0, 0, 0 ) 2224 IF ( - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk_tgt(ji,jj) ) < p_z_ref( ji, jj, jk_ref(ji,jj) ) ) THEN 2225 IF ( jk_ref(ji,jj) < jpk-1 ) jk_ref(ji,jj) = jk_ref(ji,jj) + 1 2226 END IF 2227 END_2D 2228 2229 ! if level jk_tgt lies above or at same level as jk_ref, store jk_ref for this level and add one to jk_tgt 2230 DO_2D( 0, 0, 0, 0 ) 2231 IF ( - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk_tgt(ji,jj) ) + tol_wp > p_z_ref( ji, jj, jk_ref(ji,jj) ) ) THEN 2232 IF ( jk_tgt(ji,jj) < jpk ) THEN 2233 kk_ref_for_tgt( ji, jj, jk_tgt(ji,jj) ) = jk_ref(ji,jj) 2234 jk_tgt(ji,jj) = jk_tgt(ji,jj) + 1 2235 END IF 2236 END IF 2237 END_2D 2238 2239 IF ( lwp .AND. ln_dbg_hpg ) THEN 2240 CALL dbg_2di_k( 'jk_ref', jk_ref, jk_comb ) 2241 CALL dbg_2di_k( 'jk_tgt', jk_tgt, jk_comb ) 2242 END IF 2243 2244 END DO ! jk_comb 2245 2246 2247 ! 3. Checks to make sure we do not read or write out of bounds 2248 2249 ! 3.1 First test on jk_tgt to check that it reaches the bottom level mbkt 2250 2251 jcount = 0 2252 DO_2D(0,0,0,0) 2253 IF ( jk_tgt(ji,jj) < (mbkt(ji+ki_off_tgt, jj+kj_off_tgt) + 1) ) jcount = jcount + 1 2254 END_2D 2255 2256 IF ( jcount > 0 ) THEN 2257 WRITE( numout,*) 'loc_ref_tgt: stopping because kk_ref_for_tgt will not cover all sea-points; jcount = ', jcount 2258 CALL ctl_stop( 'dyn_hpg_djr : kk_ref_for_tgt will not cover all sea-points' ) 2259 END IF 2260 2261 ! 3.2 kk_ref_for_tgt pointing to invalid level in reference profile 2262 2263 jcount = 0 2264 DO_2D(0,0,0,0) 2265 jk_bot = mbkt(ji+ki_off_tgt, jj+kj_off_tgt) 2266 IF ( kk_ref_for_tgt (ji,jj,jk_bot) > kk_bot_ref(ji,jj) ) jcount = jcount + 1 2267 END_2D 2268 2269 IF ( jcount > 0 ) THEN 2270 WRITE( numout,*) 'loc_ref_tgt: stopping because kk_ref_for_tgt points to a level below the bottom of the reference profile; jcount = ', jcount 2271 CALL ctl_stop( 'dyn_hpg_djr : kk_ref_for_tgt points to a level below the bottom of the reference profile' ) 2272 END IF 2273 2274 2275 ! 3.3 diagnostic check that the right reference levels are chosen for the target profile 2276 IF ( lwp .AND. ln_dbg_hpg ) THEN 2277 WRITE( numout,*) 'loc_ref_tgt: ki_off_tgt, kj_off_tgt = ', ki_off_tgt, kj_off_tgt 2278 CALL dbg_3dr( '-p_dep_tgt', -p_dep_tgt ) 2279 CALL dbg_3dr( 'p_z_ref', p_z_ref ) 2280 CALL dbg_3di( 'kk_ref_for_tgt', kk_ref_for_tgt ) 2281 2282 jcount = 0 2283 DO_3D(0,0,0,0,2,jpk-1) 2284 z_tgt = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk) 2285 jkr = kk_ref_for_tgt(ji,jj,jk) 2286 z_below = p_z_ref( ji, jj, jkr ) 2287 z_above = p_z_ref( ji, jj, jkr-1 ) 2288 IF ( ( z_above < z_tgt .AND. jkr > jk_init ) .OR. z_below > (z_tgt + tol_wp) ) THEN 2289 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_above 2292 END IF 2293 jcount = jcount + 1 2294 END IF 2295 END_3D 2296 WRITE(numout,*) 'loc_ref_tgt: jcount = ', jcount 2297 END IF 2298 2299 ! 3.4 Same check as in 4.2 but detailed diagnostics not written out. 2300 jcount = 0 2301 DO_3D(0,0,0,0,2,jpk-1) 2302 z_tgt = - p_dep_tgt( ji+ki_off_tgt, jj+kj_off_tgt, jk) 2303 jkr = kk_ref_for_tgt(ji,jj,jk) 2304 z_below = p_z_ref( ji, jj, jkr ) 2305 z_above = p_z_ref( ji, jj, jkr-1 ) 2306 IF ( ( z_above < z_tgt .AND. jkr > jk_init) .OR. z_below > (z_tgt + tol_wp) ) THEN 2307 jcount = jcount + 1 2308 END IF 2309 END_3D 2310 2311 IF ( jcount > 0 ) THEN 2312 WRITE( numout,*) 'loc_ref_tgt: stopping because jcount is non-zero; jcount = ', jcount 2313 CALL ctl_stop( 'dyn_hpg_djr : loc_ref_tgt failed to locate target points correctly' ) 2314 END IF 2315 2316 ! 4. Adjust kk_ref_for_tgt so that interpolation of simple cubic is off-centred at the bottom (does not require boundary conditions) 2317 ! This assumes that every sea point has at least 4 levels. 2318 2319 IF ( .NOT. ln_hpg_djr_ref_ccs ) THEN 2320 DO_3D(0, 0, 0, 0, 2, jpk-1) 2321 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 END_3D 2323 END IF 2324 2325 RETURN 2326 2327 END SUBROUTINE loc_ref_tgt 2328 2329 !---------------------------------------------------------------------------- 2330 2331 SUBROUTINE hpg_ffr( kt, Kmm, puu, pvv, Krhs ) 2332 !!--------------------------------------------------------------------- 2333 !! *** ROUTINE hpg_ffr *** 2334 !! 2335 !! ** Method : s-coordinate case forces on faces scheme. 2336 !! Local density subtracted using cubic or constrained cubic splines (ccs) 2337 !! Remaining densities interpolated to centre either linearly or with ccs 2338 !! Vertical representation either using quadratic density or classic 2nd order accurate 2339 !! 2340 !! ** Action : - Update puu(..,Krhs) and pvv(..,Krhs) with the now hydrastatic pressure trend 2341 !!---------------------------------------------------------------------- 2342 INTEGER , INTENT( in ) :: kt ! ocean time-step index 2343 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 2344 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 2345 !! 2346 2347 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zrhd_ref, zdrhd_k_ref ! densities (rhd) of reference profile 2348 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zz_ref ! heights of reference profile 2349 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4) :: z_rhd_pmr ! profiles minus reference 2350 2351 ! 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) 2354 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_ddepw_ij ! constrained spline horizontal differences in gdepw 2355 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: z_p_pmr, z_F_pmr ! pressures and forces on vertical faces 2356 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_F_upp ! forces on upper faces 2357 2358 INTEGER, DIMENSION(A2D(nn_hls)) :: jk_bot_ref ! bottom level in reference profiles 2359 INTEGER, DIMENSION(A2D(nn_hls),3) :: j_mbkt ! bottom levels for the 3 profiles in the cell 2360 2361 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpforces ! total of forces 2362 2363 INTEGER :: ji, jj, jk, j_uv, jr ! dummy loop indices 2364 INTEGER :: jio, jjo ! offsets for the 2 point stencil (0 or 1) 2365 INTEGER :: ji_ro, jj_ro, jr_offset ! offsets for the reference profile 2366 INTEGER :: jn_hor_pts ! number of profiles required in horizontal (4 if cubic interpolarion or 2 if not) 2367 2368 REAL(wp) :: z_r_6 ! 1/6 2369 REAL(wp) :: zr_a, zr_b, zr_c ! rhd evaluated at i, i+1/2 and i+1 2370 REAL(wp) :: za_0, za_1, za_2 ! polynomial coefficients 2371 !!---------------------------------------------------------------------- 2372 2373 IF( kt == nit000 ) THEN 2374 IF(lwp) WRITE(numout,*) 2375 IF(lwp) WRITE(numout,*) 'dyn:hpg_Lin : hydrostatic pressure gradient trend' 2376 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, forces on faces with reference removed' 2377 ENDIF 2378 ! 2379 2380 z_r_6 = 1.0_wp/6.0_wp 2381 2382 IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN 2383 jn_hor_pts = 4 2384 jr_offset = 2 2385 ELSE 2386 jn_hor_pts = 2 2387 jr_offset = 1 2388 END IF 2389 2390 ! WRITE( numout, *) ' hpg_ffr: kt, jn_hor_pts, jr_offset = ', kt, jn_hor_pts, jr_offset 2391 2392 DO j_uv = 1, 2 2393 2394 IF ( j_uv == 1 ) THEN 2395 jio = 1 ! set for the whole of the j_uv loop (loop ends near the bottom of the routine) 2396 jjo = 0 2397 ELSE 2398 jio = 0 2399 jjo = 1 2400 END IF 2401 2402 ! 1. Calculate the referene profile from all points in the stencil 2403 2404 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) 2406 IF ( ln_hpg_ffr_ref_ccs ) CALL calc_drhd_k(zrhd_ref, jk_bot_ref, zdrhd_k_ref) 2407 END IF 2408 2409 IF ( ln_dbg_hpg ) THEN 2410 CALL dbg_3dr( '2. rhd', rhd ) 2411 CALL dbg_3dr( '2. gdept', gdept(:,:,:,Kmm) ) 2412 END IF 2413 2414 ! 2. Interpolate reference profile to all points in the stencil and calculate z_rhd_pmr (profile minus reference) 2415 2416 DO jr = 1, jn_hor_pts 2417 2418 2419 IF ( j_uv == 1 ) THEN 2420 ji_ro = jr - jr_offset ! range of jio: is -1 to 2 for jn_hor_pts = 4; is 0 to 1 for jn_hor_pts = 2 2421 jj_ro = 0 2422 ELSE 2423 ji_ro = 0 2424 jj_ro = jr - jr_offset 2425 END IF 2426 2427 IF ( ln_hpg_ffr_ref ) THEN !!!! MJB should gde3w below now be gdept ?? !!!! 2428 2429 IF ( ln_hpg_ffr_ref_ccs ) THEN 2430 CALL ref_to_tgt_ccs ( ji_ro, jj_ro, gde3w, rhd, zz_ref, zrhd_ref, zdrhd_k_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2431 ELSE 2432 CALL ref_to_tgt_cub ( ji_ro, jj_ro, gde3w, rhd, zz_ref, zrhd_ref, jk_bot_ref, z_rhd_pmr(:,:,:,jr) ) 2433 END IF 2434 2435 ELSE ! no subtraction of reference profile 2436 2437 DO_3D (0,0,0,0,1,jpk-1) 2438 z_rhd_pmr(ji,jj,jk,jr) = rhd(ji+ji_ro,jj+jj_ro,jk) 2439 END_3D 2440 2441 END IF ! ln_hpg_ffr_ref ) THEN 2442 2443 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) 2445 z_depw_pmr(ji,jj,jk,jr) = gdepw(ji+ji_ro,jj+jj_ro,jk,Kmm) 2446 END_3D 2447 2448 IF ( ln_dbg_hpg ) THEN 2449 CALL dbg_3dr( '2. z_rhd_pmr', z_rhd_pmr(:,:,:,jr) ) 2450 END IF 2451 2452 END DO 2453 2454 ! 3. Do horizontal interpolation to form intermediate densities (either linear or cubic or constrained cubic) 2455 ! Transfers data points from reference stencil (2 or 4 point) to a 3 point horizontal stencil 2456 2457 IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub) THEN 2458 2459 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) 2462 CALL calc_dz_dij_ccs( j_uv, z_depw_pmr, z_ddepw_ij ) 2463 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) 2466 CALL calc_dz_dij_cub( j_uv, z_depw_pmr, z_ddepw_ij ) 2467 END IF 2468 2469 DO_3D (0,0,0,0,1,jpk-1) 2470 z_rhd_pmr(ji,jj,jk,1) = z_rhd_pmr(ji,jj,jk,2) 2471 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) 2474 END_3D 2475 2476 ELSE ! simple linear interpolation 2477 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) 2486 END_3D 2487 END IF 2488 2489 DO_2D (0,0,0,0) 2490 j_mbkt(ji,jj,1) = mbkt(ji, jj) 2491 j_mbkt(ji,jj,2) = MIN( mbkt(ji, jj), mbkt(ji+jio, jj+jjo) ) 2492 j_mbkt(ji,jj,3) = mbkt(ji+jio, jj+jjo) 2493 END_2D 2494 2495 IF ( ln_dbg_hpg ) THEN 2496 CALL dbg_3dr( '3. e3t ', e3t ) 2497 CALL dbg_3dr( '3. z_rhd_pmr: 1', z_rhd_pmr(:,:,:,1) ) 2498 CALL dbg_3dr( '3. z_rhd_pmr: 2', z_rhd_pmr(:,:,:,2) ) 2499 CALL dbg_3dr( '3. z_rhd_pmr: 3', z_rhd_pmr(:,:,:,3) ) 2500 END IF 2501 2502 ! 4. vertical interpolation of densities, calculating pressures and forces on vertical faces between w levels 2503 2504 IF ( ln_hpg_ffr_vrt_quad ) THEN 2505 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)) 2507 END DO ! jr 2508 2509 ELSE ! .NOT. ln_hpg_ffr_vrt_quad (simplest vertical integration scheme) 2510 2511 DO jr = 1, 3 2512 DO_2D(0,0,0,0) 2513 z_p_pmr(ji,jj,1,jr) = 0._wp 2514 END_2D 2515 DO jk = 1, jpk - 1 2516 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) 2518 END_2D 2519 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) ) 2521 END_2D 2522 END DO ! jk 2523 END DO ! jr 2524 2525 END IF ! ln_hpg_ffr_vrt_quad 2526 2527 IF ( ln_dbg_hpg ) THEN 2528 CALL dbg_3dr( '4. z_p_pmr: 1', z_p_pmr(:,:,:,1) ) 2529 CALL dbg_3dr( '4. z_p_pmr: 2', z_p_pmr(:,:,:,2) ) 2530 CALL dbg_3dr( '4. z_p_pmr: 3', z_p_pmr(:,:,:,3) ) 2531 CALL dbg_3dr( '4. z_F_pmr: 1', z_F_pmr(:,:,:,1) ) 2532 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 END IF 2537 2538 2539 ! 5. Calculate forces on the upper faces and hence on the total forces on the cells (zpforces) 2540 2541 DO_2D(0,0,0,0) 2542 z_F_upp(ji,jj,1) = 0._wp 2543 END_2D 2544 2545 IF ( ln_hpg_ffr_hor_ccs .OR. ln_hpg_ffr_hor_cub ) THEN ! use Simpson's rule 2546 DO_3D(0,0,0,0,2,jpk) 2547 z_F_upp(ji,jj,jk) = - z_r_6 * ( z_ddepw_ij(ji,jj,jk,1)*z_p_pmr(ji,jj,jk,1) & 2548 & + 4._wp*z_ddepw_ij(ji,jj,jk,2)*z_p_pmr(ji,jj,jk,2) & 2549 & + z_ddepw_ij(ji,jj,jk,3)*z_p_pmr(ji,jj,jk,3) ) 2550 END_3D 2551 ELSE ! use trapezoidal rule 2552 DO_3D(0,0,0,0,2,jpk) 2553 z_F_upp(ji,jj,jk) = 0.5_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji+jio,jj+jjo,jk,Kmm) ) * ( z_p_pmr(ji,jj,jk,1) + z_p_pmr(ji,jj,jk,3) ) 2554 & 2555 END_3D 2556 END IF 2557 2558 IF ( ln_dbg_hpg ) THEN 2559 CALL dbg_3dr( '4. z_F_upp: ', z_F_upp ) 2560 CALL dbg_3dr( '4. gdepw: ', gdepw ) 2561 CALL dbg_3dr( '4. z_ddepw_ij 1: ', z_ddepw_ij(:,:,:,1) ) 2562 CALL dbg_3dr( '4. z_ddepw_ij 2: ', z_ddepw_ij(:,:,:,2) ) 2563 CALL dbg_3dr( '4. z_ddepw_ij 3: ', z_ddepw_ij(:,:,:,3) ) 2564 2565 END IF 2566 2567 IF ( j_uv == 1 ) THEN 2568 DO_3D(0,0,0,0,1,jpk-1) 2569 zpforces(ji,jj,jk) = z_F_pmr(ji,jj,jk,1) - z_F_pmr(ji,jj,jk,3) + z_F_upp(ji,jj,jk) - z_F_upp(ji,jj,jk+1) 2570 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + zpforces(ji,jj,jk) / (e1u(ji,jj)*e3u(ji,jj,jk,Kmm)) 2571 END_3D 2572 IF ( ln_dbg_hpg ) CALL dbg_3dr( '5. u zpforces: ', zpforces ) 2573 2574 ELSE 2575 DO_3D(0,0,0,0,1,jpk-1) 2576 zpforces(ji,jj,jk) = z_F_pmr(ji,jj,jk,1) - z_F_pmr(ji,jj,jk,3) + z_F_upp(ji,jj,jk) - z_F_upp(ji,jj,jk+1) 2577 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + zpforces(ji,jj,jk) / (e1v(ji,jj)*e3v(ji,jj,jk,Kmm)) 2578 END_3D 2579 IF ( ln_dbg_hpg ) CALL dbg_3dr( '5. v zpforces: ', zpforces ) 2580 END IF 2581 2582 ! temporary output of fields for debugging etc. 2583 IF ( j_uv == 1) THEN 2584 CALL iom_put( "gdepw_hpg", gdepw(:,:,:,Kmm) ) 2585 CALL iom_put( "rhd_hpg", rhd ) 2586 CALL iom_put( "pressure", z_p_pmr(:,:,:,1) ) 2587 CALL iom_put( "u_force_west", z_F_pmr(:,:,:,1) ) 2588 CALL iom_put( "u_force_upper", z_F_upp ) 2589 ELSE 2590 CALL iom_put( "v_force_south", z_F_pmr(:,:,:,1) ) 2591 CALL iom_put( "v_force_upper", z_F_upp ) 2592 END IF 2593 2594 END DO ! j_uv 2595 ! 2596 END SUBROUTINE hpg_ffr 2597 2598 !------------------------------------------------------------------------------------------ 2599 2600 SUBROUTINE calc_rhd_ref (k_uv, kn_hor, prhd_ref, pz_ref, kk_bot_ref) 2601 2602 !!--------------------------------------------------------------------- 2603 !! *** ROUTINE calc_rhd_ref *** 2604 !! 2605 !! ** Method : Find the deepest cell within the stencil. 2606 !! (Later will extend to producing a reference profile that spans the highest and lowest points in the stencil) 2607 !! 2608 !! 2609 !! 2610 !! ** Action : - Set prhd_ref, pz_ref, kk_bot_ref 2611 !!---------------------------------------------------------------------- 2612 INTEGER , INTENT( in ) :: k_uv ! 1 for u-vel; 2 for v_vel 2613 INTEGER , INTENT( in ) :: kn_hor ! 4 for cubic; 2 for linear 2614 2615 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: prhd_ref ! densities of reference profile 2616 REAL(wp), DIMENSION(:,:,:), INTENT(out) :: pz_ref ! heights of " " 2617 INTEGER, DIMENSION(:,:), INTENT(out) :: kk_bot_ref ! bottom level of " " 2618 2619 INTEGER, DIMENSION(A2D(nn_hls)) :: ji_ref, jj_ref 2620 2621 INTEGER ji, jj, jk ! standard loop indices 2622 INTEGER jio, jjo ! offsets 2623 INTEGER jib, jjb ! second set of indices to deeper points 2624 INTEGER jir, jjr ! indices of reference profile 2625 REAL(wp) zhta, zhtb 2626 2627 IF (k_uv == 1) THEN 2628 jio = 1 2629 jjo = 0 2630 ELSE 2631 jio = 0 2632 jjo = 1 2633 END IF 2634 2635 2636 DO_2D( 0, 0, 0, 0 ) 2637 2638 IF ( ht_0(ji+jio,jj+jjo) >= ht_0(ji,jj) ) THEN 2639 zhta = ht_0(ji+jio,jj+jjo) 2640 ji_ref(ji,jj) = ji+jio 2641 jj_ref(ji,jj) = jj+jjo 2642 ELSE 2643 zhta = ht_0(ji,jj) 2644 ji_ref(ji,jj) = ji 2645 jj_ref(ji,jj) = jj 2646 END IF 2647 2648 IF ( kn_hor == 4 ) THEN 2649 IF ( ht_0(ji-jio,jj-jjo) >= ht_0(ji+2*jio,jj+2*jjo) ) THEN 2650 zhtb = ht_0(ji-jio,jj-jjo) 2651 jib = ji-jio 2652 jjb = jj-jjo 2653 ELSE 2654 zhtb = ht_0(ji+2*jio,jj+2*jjo) 2655 jib = ji+2*jio 2656 jjb = jj+2*jjo 2657 END IF 2658 IF ( zhta < zhtb ) THEN 2659 ji_ref(ji,jj) = jib 2660 jj_ref(ji,jj) = jjb 2661 END IF 2662 END IF 2663 2664 END_2D 2665 2666 DO_3D( 0, 0, 0, 0, 1, jpk-1) 2667 jir = ji_ref(ji,jj) 2668 jjr = jj_ref(ji,jj) 2669 prhd_ref (ji,jj,jk) = rhd(jir, jjr, jk) 2670 pz_ref (ji,jj,jk) = - gde3w(jir, jjr, jk) 2671 END_3D 2672 2673 DO_2D( 0, 0, 0, 0) 2674 jir = ji_ref(ji,jj) 2675 jjr = jj_ref(ji,jj) 2676 kk_bot_ref(ji,jj) = mbkt(jir,jjr) 2677 END_2D 2678 2679 END SUBROUTINE calc_rhd_ref 2680 2681 !------------------------------------------------------------------------------------------ 2682 2683 SUBROUTINE calc_drhd_k(p_rhd, kk_bot, p_drhd_k) 2684 2685 !!--------------------------------------------------------------------- 2686 !! *** ROUTINE calc_drhd_k *** 2687 !! 2688 !! ** Method : Calculate harmonic averages of vertical differences and apply upper and lower boundary conditions 2689 !! 2690 !! ** Action : - Set p_drhd_k 2691 !!---------------------------------------------------------------------- 2692 2693 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in) :: p_rhd ! densities of profile 2694 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: kk_bot ! bottom level of profile 2695 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_drhd_k ! harmonic mean of vertical differences of profile 2696 2697 INTEGER ji, jj, jk ! standard loop indices 2698 INTEGER jio, jjo ! offsets 2699 INTEGER iktb ! index of the bottom of ref profile 2700 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztmp ! primitive vertical differences 2701 2702 REAL(wp) cffw, z1_cff, zep 2703 2704 DO_3D( 0, 0, 0, 0, 2, jpk ) 2705 ztmp(ji,jj,jk) = p_rhd(ji,jj,jk) - p_rhd(ji,jj,jk-1) 2706 END_3D 2707 2708 zep = 1.e-15 2709 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 2710 cffw = MAX( 2._wp * ztmp(ji,jj,jk) * ztmp(ji,jj,jk+1), 0._wp ) 2711 z1_cff = ztmp(ji,jj,jk) + ztmp(ji,jj,jk+1) 2712 p_drhd_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 2713 END_3D 2714 2715 DO_2D( 0, 0, 0, 0 ) 2716 p_drhd_k(ji,jj,1) = aco_bc_rhd_srf * ( p_rhd(ji,jj,2) - p_rhd(ji,jj,1) ) - bco_bc_rhd_srf * p_drhd_k(ji,jj,2) 2717 iktb = kk_bot(ji,jj) 2718 IF ( iktb > 1 ) THEN 2719 p_drhd_k(ji,jj,iktb) = aco_bc_rhd_bot * (p_rhd(ji,jj,iktb) - p_rhd(ji,jj,iktb-1) ) - bco_bc_rhd_bot * p_drhd_k(ji,jj,iktb-1) 2720 END IF 2721 END_2D 2722 2723 END SUBROUTINE calc_drhd_k 2724 2725 !---------------------------------------------------------------------------- 2726 2727 SUBROUTINE calc_mid_ccs( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 2728 2729 !!--------------------------------------------------------------------- 2730 !! *** ROUTINE calc_mid_ccs *** 2731 !! 2732 !! ** 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 !! 2735 !! ** Action : - set p_fld_mid 2736 !!---------------------------------------------------------------------- 2737 2738 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2739 REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation) 2740 REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation) 2741 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_fld_pmr ! field in pmr form 2742 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(out) :: p_fld_mid ! field at mid-point 2743 2744 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3) :: zz_dfld_ij ! constrained spline horizontal differences (dimension 3 for consistency with calc_dfld_cen_dij) 2745 2746 INTEGER ji, jj, jk ! standard loop indices 2747 INTEGER :: j_t_levs ! indicator that field passed is on tracer levels 2748 2749 j_t_levs = 1 2750 2751 DO jk = 1, jpk 2752 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 ) 2754 2755 ! first contribution is simple centred average; p_rhd_pmr has 4 points; 2 and 3 are the central ones 2756 DO_2D( 0, 0, 0, 0) 2757 p_fld_mid(ji,jj,jk) = 0.5_wp * ( p_fld_pmr(ji,jj,jk,2) + p_fld_pmr(ji,jj,jk,3) ) 2758 END_2D 2759 2760 IF ( k_uv == 1 ) THEN 2761 DO_2D( 0, 0, 0, 0) 2762 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 2763 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 END IF 2765 END_2D 2766 ELSE ! k_uv == 2 2767 DO_2D( 0, 0, 0, 0) 2768 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 2769 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 END IF 2771 END_2D 2772 END IF ! k_uv 2773 2774 END DO ! jk 2775 2776 END SUBROUTINE calc_mid_ccs 2777 2778 !---------------------------------------------------------------------------- 2779 2780 SUBROUTINE calc_mid_cub( k_uv, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_fld_mid ) 2781 2782 !!--------------------------------------------------------------------- 2783 !! *** ROUTINE calc_mid_cub *** 2784 !! 2785 !! ** Method : Use simple cubic polynomial to interpolate 4 profiles to their central point 2786 !! The coefficients p_aco_bc_fld_hor, p_bco_bc_fld_hor are not currently used. 2787 !! The simplest form of von Neumann conditions horizontal bc is used (one could off-centred polynomials) 2788 !! ** Action : - set p_fld_mid 2789 !!---------------------------------------------------------------------- 2790 2791 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2792 REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation) 2793 REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation) 2794 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(inout) :: p_fld_pmr ! field in pmr form 2795 REAL(wp), DIMENSION(A2D(nn_hls),jpk) , INTENT(out) :: p_fld_mid ! field at mid-point 2796 2797 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences 2798 REAL(wp), DIMENSION(A2D(nn_hls),3) :: zz_dfld_ij ! constrained spline horizontal differences (dimension 3 for consistency with calc_dfld_cen_dij) 2799 2800 INTEGER ji, jj, jk ! standard loop indices 2801 REAL(wp) z_1_16, z_9_16 ! temporary sum and products 2802 REAL(wp) z_cor ! correction to the central value 2803 2804 z_1_16 = 1._wp / 16._wp ; z_9_16 = 9._wp / 16._wp 2805 2806 DO jk = 1, jpk 2807 2808 ! first contribution is simple centred average plus 1/16 of it; p_rhd_pmr has 4 points; 2 and 3 are the central ones 2809 DO_2D( 0, 0, 0, 0) 2810 p_fld_mid(ji,jj,jk) = z_9_16 * ( p_fld_pmr(ji,jj,jk,2) + p_fld_pmr(ji,jj,jk,3) ) 2811 END_2D 2812 2813 IF ( k_uv == 1 ) THEN 2814 DO_2D( 0, 0, 0, 0) 2815 IF ( umask(ji-1, jj, jk) > 0.5 ) THEN 2816 z_cor = p_fld_pmr(ji,jj,jk,1) 2817 ELSE 2818 z_cor = p_fld_pmr(ji,jj,jk,2) 2819 END IF 2820 IF ( umask(ji+1, jj, jk) > 0.5 ) THEN 2821 z_cor = z_cor + p_fld_pmr(ji,jj,jk,4) 2822 ELSE 2823 z_cor = z_cor + p_fld_pmr(ji,jj,jk,3) 2824 END IF 2825 p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - z_1_16 * z_cor 2826 END_2D 2827 ELSE ! k_uv == 2 2828 DO_2D( 0, 0, 0, 0) 2829 IF ( vmask(ji, jj-1, jk) > 0.5 ) THEN 2830 z_cor = p_fld_pmr(ji,jj,jk,1) 2831 ELSE 2832 z_cor = p_fld_pmr(ji,jj,jk,2) 2833 END IF 2834 IF ( vmask(ji, jj+1, jk) > 0.5 ) THEN 2835 z_cor = z_cor + p_fld_pmr(ji,jj,jk,4) 2836 ELSE 2837 z_cor = z_cor + p_fld_pmr(ji,jj,jk,3) 2838 END IF 2839 p_fld_mid(ji,jj,jk) = p_fld_mid(ji,jj,jk) - z_1_16 * z_cor 2840 END_2D 2841 END IF ! k_uv 2842 2843 END DO ! jk 2844 2845 END SUBROUTINE calc_mid_cub 2846 2847 !---------------------------------------------------------------------------- 2848 2849 SUBROUTINE calc_dz_dij_ccs( k_uv, p_z_pmr, p_dz_ij ) 2850 2851 !!--------------------------------------------------------------------- 2852 !! *** ROUTINE calc_dz_dij_ccs *** 2853 !! 2854 !! ** Method : Use constrained cubic spline to determine horizontal derivatives of z at 3 central points 2855 !! The routine is limited to z derivatives because the output is valid on w-levels 2856 !! ** Action : - set p_dz_ij 2857 !!---------------------------------------------------------------------- 2858 2859 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2860 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_z_pmr ! z field in pmr form 2861 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dz_ij ! constrained cubic horizontal derivatives at -1/2, 0 and 1/2 2862 2863 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences 2864 2865 INTEGER ji, jj, jk ! standard loop indices 2866 INTEGER jk_msk ! jk level to use for umask and vmask (we need z on the upper and lower faces) 2867 INTEGER j_w_levs ! indicator for data on w-levels 2868 2869 j_w_levs = 2 2870 2871 DO jk = 1, jpk 2872 2873 CALL calc_dfld_pmr_ij( k_uv, j_w_levs, jk, aco_bc_z_hor, bco_bc_z_hor, p_z_pmr, p_dz_ij ) 2874 2875 ! copy the 2nd horizontal element to the 3rd element of the output array (for an isolated velocity cell p_dz_ij is the same for all 3 elements) 2876 DO_2D( 0, 0, 0, 0) 2877 p_dz_ij(ji,jj,jk,3) = p_dz_ij(ji,jj,jk,2) 2878 END_2D 2879 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) 2881 IF ( jk == 1 ) THEN 2882 jk_msk = 1 2883 ELSE 2884 jk_msk = jk - 1 2885 END IF 2886 2887 IF ( k_uv == 1 ) THEN 2888 DO_2D( 0, 0, 0, 0) 2889 IF ( umask(ji-1, jj, jk_msk) > 0.5 .OR. umask(ji+1, jj, jk_msk) > 0.5 ) THEN 2890 p_dz_ij(ji,jj,jk,2) = 1.5_wp*( p_z_pmr(ji,jj,jk,3) - p_z_pmr(ji,jj,jk,2) ) - 0.25_wp * ( p_dz_ij(ji,jj,jk,3) + p_dz_ij(ji,jj,jk,1) ) 2891 END IF 2892 END_2D 2893 ELSE ! k_uv == 2 2894 DO_2D( 0, 0, 0, 0) 2895 IF ( vmask(ji, jj-1, jk_msk) > 0.5 .OR. vmask(ji, jj+1, jk_msk) > 0.5 ) THEN 2896 p_dz_ij(ji,jj,jk,2) = 1.5_wp*( p_z_pmr(ji,jj,jk,3) - p_z_pmr(ji,jj,jk,2) ) - 0.25_wp * ( p_dz_ij(ji,jj,jk,3) + p_dz_ij(ji,jj,jk,1) ) 2897 END IF 2898 END_2D 2899 END IF ! k_uv 2900 2901 END DO ! jk 2902 2903 END SUBROUTINE calc_dz_dij_ccs 2904 2905 !---------------------------------------------------------------------------- 2906 2907 SUBROUTINE calc_dz_dij_cub( k_uv, p_z_pmr, p_dz_ij ) 2908 2909 !!--------------------------------------------------------------------- 2910 !! *** ROUTINE calc_dz_dij_cub *** 2911 !! 2912 !! ** Method : Use simple cubic polynomial to determine horizontal derivatives at 3 central points 2913 !! based on equations (5.5) and (5.6) of SMcW 2003 2914 !! The routine is limited to z derivatives because the output is valid on w-levels 2915 !! ** Action : - set p_dz_ij 2916 !!---------------------------------------------------------------------- 2917 2918 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 2919 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_z_pmr ! z field in pmr form 2920 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dz_ij ! horizontal derivatives at -1/2, 0 and 1/2 2921 2922 REAL(wp), DIMENSION(A2D(nn_hls)) :: z_z_a, z_z_d ! z fields with boundary conditions applied 2923 REAL(wp) z_1_24 ! 1/24 2924 REAL(wp) z_a, z_b, z_c, z_d ! values of input field at the four points used (-3/2, -1/2, 1/2, 3/2) 2925 REAL(wp) z_c_m_b, z_d_m_a ! differences c - b and d - a 2926 REAL(wp) z_co_1, z_co_2, z_co_3 ! coefficients of polynomial (field = z_co_0 + z_co_1 zeta + .. ) 2927 INTEGER ji, jj, jk ! standard loop indices 2928 INTEGER jk_msk ! jk level to use for umask and vmask (we need z on the upper and lower faces) 2929 2930 z_1_24 = 1.0_wp / 24.0_wp 2931 2932 DO jk = 1, jpk 2933 2934 !------------------------------------------------------ 2935 ! 1. Use simple von Neumann conditions at the boundaries 2936 !------------------------------------------------------ 2937 IF ( jk == 1 ) THEN 2938 jk_msk = 1 2939 ELSE 2940 jk_msk = jk - 1 2941 END IF 2942 2943 IF ( k_uv == 1 ) THEN 2944 DO_2D( 0, 0, 0, 0) 2945 IF ( umask(ji-1, jj, jk_msk) > 0.5 ) THEN 2946 z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,1) 2947 ELSE 2948 z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,2) 2949 END IF 2950 IF ( umask(ji+1, jj, jk_msk) > 0.5 ) THEN 2951 z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,4) 2952 ELSE 2953 z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,3) 2954 END IF 2955 END_2D 2956 ELSE ! k_uv == 2 2957 DO_2D( 0, 0, 0, 0) 2958 IF ( vmask(ji, jj-1, jk_msk) > 0.5 ) THEN 2959 z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,1) 2960 ELSE 2961 z_z_a(ji,jj) = p_z_pmr(ji,jj,jk,2) 2962 END IF 2963 IF ( vmask(ji, jj+1, jk_msk) > 0.5 ) THEN 2964 z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,4) 2965 ELSE 2966 z_z_d(ji,jj) = p_z_pmr(ji,jj,jk,3) 2967 END IF 2968 END_2D 2969 END IF ! k_uv 2970 2971 !------------------------------------------------------ 2972 ! 2. calculate the coefficients for the polynomial fit (c_1, c_2 and c_3; we don't need c_0) 2973 !------------------------------------------------------ 2974 2975 DO_2D( 0, 0, 0, 0) 2976 z_a = z_z_a(ji,jj) 2977 z_b = p_z_pmr(ji,jj,jk,2) 2978 z_c = p_z_pmr(ji,jj,jk,3) 2979 z_d = z_z_d(ji,jj) 2980 2981 z_c_m_b = z_c - z_b ; z_d_m_a = z_d - z_a 2982 2983 ! eqn (5.6) of SMcW 2984 z_co_1 = 1.125_wp * z_c_m_b - z_1_24 * z_d_m_a 2985 z_co_2 = -0.5_wp * (z_c+z_b) + 0.5_wp * (z_d+z_a) 2986 z_co_3 = -3.0_wp * z_c_m_b + z_d_m_a 2987 2988 ! eqn (5.5) of SMcW (first derivative of it) 2989 p_dz_ij(ji,jj,jk,1) = z_co_1 - 0.5_wp*z_co_2 + 0.125_wp*z_co_3 ! zeta = -0.5 2990 p_dz_ij(ji,jj,jk,2) = z_co_1 ! zeta = 0.0 2991 p_dz_ij(ji,jj,jk,3) = z_co_1 + 0.5_wp*z_co_2 + 0.125_wp*z_co_3 ! zeta = 0.5 2992 2993 END_2D 2994 2995 END DO ! jk 2996 2997 END SUBROUTINE calc_dz_dij_cub 2998 2999 !---------------------------------------------------------------------------- 3000 3001 SUBROUTINE calc_dfld_pmr_ij( k_uv, k_lev_type, kk, p_aco_bc_fld_hor, p_bco_bc_fld_hor, p_fld_pmr, p_dfld_ij ) 3002 3003 !!--------------------------------------------------------------------- 3004 !! *** ROUTINE calc_dfld_pmr_ij *** 3005 !! 3006 !! ** Method : Calculate constrained spline horizontal derivatives 3007 !! 3008 !! ** Action : - set p_dfld_ij 3009 !!---------------------------------------------------------------------- 3010 3011 INTEGER , INTENT(in) :: k_uv ! 1 for u-vel; 2 for v_vel 3012 INTEGER , INTENT(in) :: k_lev_type ! 1 for t-level data; 2 for w-level data; used with check of land/sea mask 3013 INTEGER , INTENT(in) :: kk ! vertical level 3014 REAL(wp) , INTENT(in) :: p_aco_bc_fld_hor ! a coeff for horizontal bc (von Neumann or linear extrapolation) 3015 REAL(wp) , INTENT(in) :: p_bco_bc_fld_hor ! b coeff for horizontal bc (von Neumann or linear extrapolation) 3016 REAL(wp), DIMENSION(A2D(nn_hls),jpk,4), INTENT(in) :: p_fld_pmr ! field in pmr form 3017 REAL(wp), DIMENSION(A2D(nn_hls),jpk,3), INTENT(out) :: p_dfld_ij ! constrained spline horizontal derivatives of the field; only 1 and 2 are set! 3018 3019 REAL(wp), DIMENSION(A2D(nn_hls)) :: zdfld_21, zdfld_32, zdfld_43 ! primitive horizontal differences 3020 3021 INTEGER ji, jj ! standard loop indices 3022 INTEGER jio, jjo ! offsets 3023 INTEGER iktb ! index of the bottom of ref profile 3024 3025 REAL(wp) z1_cff, z_cff_31, z_cff_42 ! temporary sum and products 3026 REAL(wp) zep 3027 3028 INTEGER :: j_t_levs, jk_msk ! indicators that field passed in is valid on t-levels or w-levels; jk to use with masks 3029 3030 j_t_levs = 1 3031 3032 !---------------------------------------------------------------------------------------- 3033 ! 1. compute and store elementary horizontal differences zfor z_rhd_pmr arrays 3034 !---------------------------------------------------------------------------------------- 3035 3036 IF ( k_uv == 1) THEN 3037 jio = 1 3038 jjo = 0 3039 ELSE 3040 jio = 0 3041 jjo = 1 3042 END IF 3043 3044 IF ( k_lev_type == j_t_levs ) THEN 3045 jk_msk = kk 3046 ELSE 3047 IF ( kk == 1 ) THEN 3048 jk_msk = 1 3049 ELSE 3050 jk_msk = kk - 1 3051 END IF 3052 END IF 3053 3054 DO_2D( 0, 0, 0, 0 ) 3055 zdfld_21(ji,jj) = p_fld_pmr(ji,jj,kk,2) - p_fld_pmr(ji,jj,kk,1) 3056 zdfld_32(ji,jj) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) 3057 zdfld_43(ji,jj) = p_fld_pmr(ji,jj,kk,4) - p_fld_pmr(ji,jj,kk,3) 3058 END_2D 3059 3060 zep = 1.e-15 3061 DO_2D( 0, 0, 0, 0 ) 3062 z_cff_31 = MAX( 2._wp * zdfld_21(ji,jj) * zdfld_32(ji,jj), 0._wp ) 3063 z1_cff = zdfld_21(ji,jj) + zdfld_32(ji,jj) 3064 p_dfld_ij(ji,jj,kk,1) = z_cff_31 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 3065 3066 z_cff_42 = MAX( 2._wp * zdfld_32(ji,jj) * zdfld_43(ji,jj), 0._wp ) 3067 z1_cff = zdfld_32(ji,jj) + zdfld_43(ji,jj) 3068 p_dfld_ij(ji,jj,kk,2) = z_cff_42 / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 3069 END_2D 3070 3071 !---------------------------------------------------------------------------------- 3072 ! 2. apply boundary conditions at side boundaries using 5.36-5.37 3073 !---------------------------------------------------------------------------------- 3074 3075 IF ( k_uv == 1 ) THEN 3076 DO_2D( 0, 0, 0, 0) 3077 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 3078 IF ( umask(ji,jj,jk_msk) > 0.5_wp .AND. umask(ji-1,jj,jk_msk) < 0.5_wp .AND. umask(ji+1,jj,jk_msk) > 0.5_wp) THEN 3079 p_dfld_ij(ji,jj,kk,1) = p_aco_bc_fld_hor * ( p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,2) 3080 END IF 3081 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 3082 IF ( umask(ji,jj,jk_msk) < 0.5_wp .AND. umask(ji-1,jj,jk_msk) > 0.5_wp .AND. umask(ji-2,jj,jk_msk) > 0.5_wp) THEN 3083 p_dfld_ij(ji,jj,kk,2) = p_aco_bc_fld_hor * ( p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,1) 3084 END IF 3085 END_2D 3086 3087 DO_2D( 0, 0, 0, 0) 3088 ! For an isolated velocity point use the central difference (this is important) 3089 IF ( umask(ji-1, jj, jk_msk) < 0.5 .AND. umask(ji+1, jj, jk_msk) < 0.5 ) THEN 3090 p_dfld_ij(ji,jj,kk,1) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) 3091 p_dfld_ij(ji,jj,kk,2) = p_dfld_ij(ji,jj,kk,1) 3092 END IF 3093 END_2D 3094 3095 ELSE ! k_uv == 2 3096 3097 DO_2D( 0, 0, 0, 0) 3098 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 3099 IF ( vmask(ji,jj,jk_msk) > 0.5_wp .AND. vmask(ji,jj-1,jk_msk) < 0.5_wp .AND. vmask(ji,jj+1,jk_msk) > 0.5_wp) THEN 3100 p_dfld_ij(ji,jj,kk,1) = p_aco_bc_fld_hor * (p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,2) 3101 END IF 3102 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 3103 IF ( vmask(ji,jj,jk_msk) < 0.5_wp .AND. vmask(ji,jj-1,jk_msk) > 0.5_wp .AND. vmask(ji,jj-2,jk_msk) > 0.5_wp) THEN 3104 p_dfld_ij(ji,jj,kk,2) = p_aco_bc_fld_hor * (p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) ) - p_bco_bc_fld_hor * p_dfld_ij(ji,jj,kk,1) 3105 END IF 3106 END_2D 3107 3108 DO_2D( 0, 0, 0, 0) 3109 ! For an isolated velocity point use the central difference (this is important) 3110 IF ( vmask(ji, jj-1, jk_msk) < 0.5 .AND. vmask(ji, jj+1, jk_msk) < 0.5 ) THEN 3111 p_dfld_ij(ji,jj,kk,1) = p_fld_pmr(ji,jj,kk,3) - p_fld_pmr(ji,jj,kk,2) 3112 p_dfld_ij(ji,jj,kk,2) = p_dfld_ij(ji,jj,kk,1) 3113 END IF 3114 END_2D 3115 3116 END IF ! k_uv 3117 3118 END SUBROUTINE calc_dfld_pmr_ij 3119 3120 !------------------------------------------------------------------------------------------ 3121 3122 SUBROUTINE vrt_int_quad(k_mbkt, p_rhd, p_e3t, p_p, p_F) 3123 3124 !!--------------------------------------------------------------------- 3125 !! *** ROUTINE calc_rhd_k *** 3126 !! 3127 !! ** Method : Use quadratic reconstruction of p_rhd (density profile) to calculate pressure profile and 3128 !! forces on vertical faces between w levels 3129 !! 3130 !! ** Action : - set p_p and p_F 3131 !!---------------------------------------------------------------------- 3132 3133 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(in) :: k_mbkt ! bottom level 3134 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 3136 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_p ! pressures on w levels 3137 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(out) :: p_F ! force on the vertical face between w levels 3138 3139 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 3152 3153 ! Integrate densities in the vertical using quadratic fits to the densities 3154 ! zr_a is r_{-1} zr_b is r_1 ; zr_c is r_3 3155 3156 ! 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 3198 END SUBROUTINE vrt_int_quad 3199 3200 !------------------------------------------------------------------------------------------ 3201 3202 SUBROUTINE dbg_2di( cc_array, ki_2d ) 3203 CHARACTER*(*), INTENT(IN) :: cc_array 3204 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(IN) :: ki_2d 3205 3206 INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp 3207 INTEGER :: ji_prt, jj_prt 3208 3209 IF( lwp ) THEN 3210 ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min ) 3211 ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max ) 3212 jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min ) 3213 jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max ) 3214 3215 ! print out a 2D horizontal slice 3216 3217 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN 3218 IF ( jj_prt_low <= jj_prt_upp ) THEN 3219 WRITE(numout,*) 3220 WRITE(numout,*) cc_array 3221 WRITE(numout,*) ' ji_prt_low, ji_prt_upp = ', ji_prt_low, ji_prt_upp 3222 WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3223 DO jj_prt = jj_prt_low, jj_prt_upp 3224 WRITE(numout,*) jj_prt, ( ki_2d(ji_prt, jj_prt), ji_prt = ji_prt_low, ji_prt_upp ) 3225 END DO 3226 END IF 3227 END IF 3228 3229 END IF 3230 3231 RETURN 3232 END SUBROUTINE dbg_2di 3233 3234 !---------------------------------------------------------------------------- 3235 3236 SUBROUTINE dbg_2di_k( cc_array, ki_2d, kk ) 3237 CHARACTER*(*), INTENT(IN) :: cc_array 3238 INTEGER, DIMENSION(A2D(nn_hls)), INTENT(IN) :: ki_2d 3239 INTEGER, INTENT(IN) :: kk 3240 3241 INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp 3242 INTEGER :: ji_prt, jj_prt, jk_prt 3243 3244 IF( lwp ) THEN 3245 ji_prt_low = MAX( ntsi, ki_dbg_min ) 3246 ji_prt_upp = MIN( ntei, ki_dbg_max ) 3247 jj_prt_low = MAX( ntsj, kj_dbg_min ) 3248 jj_prt_upp = MIN( ntej, kj_dbg_max ) 3249 3250 ! print out a 2D (ji, jk) slice 3251 3252 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN 3253 IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN 3254 IF ( kk == kk_dbg_min ) THEN 3255 WRITE(numout,*) 3256 WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen 3257 WRITE(numout,*) 'ji_prt_low, ji_prt_upp = ', ji_prt_low, ji_prt_upp 3258 WRITE(numout,*) 'kk_dbg_min, kk_dbg_max = ', kk_dbg_min, kk_dbg_max 3259 END IF 3260 IF ( kk_dbg_min <= kk .AND. kk <= kk_dbg_max ) & 3261 & WRITE(numout,*) cc_array, kk, ( ki_2d(ji_prt, kj_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp ) 3262 END IF 3263 END IF 3264 3265 ! print out a 2D (jj, jk) slice 3266 3267 IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN 3268 IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN 3269 IF ( kk == kk_dbg_min ) THEN 3270 WRITE(numout,*) 3271 WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen 3272 WRITE(numout,*) 'jj_prt_low, jj_prt_upp = ', jj_prt_low, jj_prt_upp 3273 WRITE(numout,*) 'kk_dbg_min, kk_dbg_max = ', kk_dbg_min, kk_dbg_max 3274 END IF 3275 IF ( kk_dbg_min <= kk .AND. kk <= kk_dbg_max ) & 3276 & WRITE(numout,*) cc_array, kk, ( ki_2d(ki_dbg_cen, jj_prt), jj_prt = jj_prt_low, jj_prt_upp ) 3277 END IF 3278 END IF 3279 3280 END IF 3281 3282 RETURN 3283 END SUBROUTINE dbg_2di_k 3284 3285 !---------------------------------------------------------------------------- 3286 3287 SUBROUTINE dbg_2dr( cc_array, pr_2d ) 3288 CHARACTER*(*), INTENT(IN) :: cc_array 3289 REAL(wp), DIMENSION(A2D(nn_hls)), INTENT(IN) :: pr_2d 3290 3291 INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp 3292 INTEGER :: ji_prt, jj_prt 3293 3294 IF(lwp) THEN 3295 ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min ) 3296 ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max ) 3297 jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min ) 3298 jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max ) 3299 3300 ! print out a 2D horizontal slice 3301 3302 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN 3303 IF ( jj_prt_low <= jj_prt_upp ) THEN 3304 WRITE(numout,*) 3305 WRITE(numout,*) cc_array 3306 WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3307 DO jj_prt = jj_prt_low, jj_prt_upp 3308 WRITE(numout,*) jj_prt, ( pr_2d(ji_prt, jj_prt), ji_prt = ji_prt_low, ji_prt_upp ) 3309 END DO 3310 END IF 3311 END IF 3312 3313 END IF 3314 3315 RETURN 3316 END SUBROUTINE dbg_2dr 3317 3318 !---------------------------------------------------------------------------- 3319 3320 SUBROUTINE dbg_3di( cc_array, ki_3d ) 3321 3322 CHARACTER*(*), INTENT(IN) :: cc_array 3323 INTEGER, DIMENSION(A2D(nn_hls),jpk), INTENT(IN) :: ki_3d 3324 3325 INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp 3326 INTEGER :: ji_prt, jj_prt, jk_prt 3327 3328 ! output values 3329 3330 3331 IF(lwp) THEN 3332 ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min ) 3333 ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max ) 3334 jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min ) 3335 jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max ) 3336 3337 ! print out a 2D horizontal slice 3338 3339 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN 3340 IF ( jj_prt_low <= jj_prt_upp ) THEN 3341 WRITE(numout,*) 3342 WRITE(numout,*) cc_array, ' level = ', kk_dbg_cen 3343 WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3344 DO jj_prt = jj_prt_low, jj_prt_upp 3345 WRITE(numout,*) jj_prt, ( ki_3d(ji_prt, jj_prt, kk_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp ) 3346 END DO 3347 END IF 3348 END IF 3349 3350 ! print out a 2D (ji, jk) slice 3351 3352 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN 3353 IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN 3354 WRITE(numout,*) 3355 WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen 3356 WRITE(numout,*) ' row/lev ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3357 DO jk_prt = kk_dbg_min, kk_dbg_max 3358 WRITE(numout,*) jk_prt, ( ki_3d(ji_prt, kj_dbg_cen, jk_prt), ji_prt = ji_prt_low, ji_prt_upp ) 3359 END DO 3360 END IF 3361 END IF 3362 3363 ! print out a 2D (jj, jk) slice 3364 3365 IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN 3366 IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN 3367 WRITE(numout,*) 3368 WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen 3369 WRITE(numout,*) ' col/lev ', ( jj_prt, jj_prt = jj_prt_low, jj_prt_upp ) 3370 DO jk_prt = kk_dbg_min, kk_dbg_max 3371 WRITE(numout,*) jk_prt, ( ki_3d(ki_dbg_cen, jj_prt, jk_prt), jj_prt = jj_prt_low, jj_prt_upp ) 3372 END DO 3373 END IF 3374 END IF 3375 3376 END IF 3377 3378 RETURN 3379 END SUBROUTINE dbg_3di 3380 3381 !---------------------------------------------------------------------------- 3382 3383 SUBROUTINE dbg_3dr( cc_array, pr_3d ) 3384 3385 CHARACTER*(*), INTENT(IN) :: cc_array 3386 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(IN) :: pr_3d 3387 3388 INTEGER :: ji_prt_low, ji_prt_upp, jj_prt_low, jj_prt_upp 3389 INTEGER :: ji_prt, jj_prt, jk_prt 3390 3391 ! output values 3392 3393 3394 IF(lwp) THEN 3395 ji_prt_low = MAX( ntsi-nn_hls, ki_dbg_min ) 3396 ji_prt_upp = MIN( ntei+nn_hls, ki_dbg_max ) 3397 jj_prt_low = MAX( ntsj-nn_hls, kj_dbg_min ) 3398 jj_prt_upp = MIN( ntej+nn_hls, kj_dbg_max ) 3399 3400 ! print out a 2D horizontal slice 3401 3402 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ij ) THEN 3403 IF ( jj_prt_low <= jj_prt_upp ) THEN 3404 WRITE(numout,*) 3405 WRITE(numout,*) cc_array, ' level = ', kk_dbg_cen 3406 WRITE(numout,*) ' row/col ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3407 DO jj_prt = jj_prt_low, jj_prt_upp 3408 WRITE(numout,*) jj_prt, ( pr_3d(ji_prt, jj_prt, kk_dbg_cen), ji_prt = ji_prt_low, ji_prt_upp ) 3409 END DO 3410 END IF 3411 END IF 3412 3413 ! print out a 2D (ji, jk) slice 3414 3415 IF ( ji_prt_low <= ji_prt_upp .AND. ln_dbg_ik ) THEN 3416 IF ( ntsj <= kj_dbg_cen .AND. kj_dbg_cen <= ntej ) THEN 3417 WRITE(numout,*) 3418 WRITE(numout,*) cc_array, ' kj = ', kj_dbg_cen 3419 WRITE(numout,*) ' row/lev ', ( ji_prt, ji_prt = ji_prt_low, ji_prt_upp ) 3420 DO jk_prt = kk_dbg_min, kk_dbg_max 3421 WRITE(numout,*) jk_prt, ( pr_3d(ji_prt, kj_dbg_cen, jk_prt), ji_prt = ji_prt_low, ji_prt_upp ) 3422 END DO 3423 END IF 3424 END IF 3425 3426 ! print out a 2D (jj, jk) slice 3427 3428 IF ( jj_prt_low <= jj_prt_upp .AND. ln_dbg_jk ) THEN 3429 IF ( ntsi <= ki_dbg_cen .AND. ki_dbg_cen <= ntei ) THEN 3430 WRITE(numout,*) 3431 WRITE(numout,*) cc_array, ' ki = ', ki_dbg_cen 3432 WRITE(numout,*) ' col/lev ', ( jj_prt, jj_prt = jj_prt_low, jj_prt_upp ) 3433 DO jk_prt = kk_dbg_min, kk_dbg_max 3434 WRITE(numout,*) jk_prt, ( pr_3d(ki_dbg_cen, jj_prt, jk_prt), jj_prt = jj_prt_low, jj_prt_upp ) 3435 END DO 3436 END IF 3437 END IF 3438 3439 END IF 3440 3441 RETURN 3442 END SUBROUTINE dbg_3dr 3443 3444 1437 3445 !!====================================================================== 1438 3446 END MODULE dynhpg
Note: See TracChangeset
for help on using the changeset viewer.
