- Timestamp:
- 2013-11-20T17:28:04+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r4153 r4292 6 6 !! History : 2.0 ! 2006-06 (B. Levier, L. Marie) original code 7 7 !! 3.1 ! 2009-02 (G. Madec, M. Leclair, R. Benshila) pure z* coordinate 8 !! ----------------------------------------------------------------------9 #if defined key_vvl 8 !! 3.3 ! 2011-10 (M. Leclair) totally rewrote domvvl: 9 !! vvl option includes z_star and z_tilde coordinates 10 10 !!---------------------------------------------------------------------- 11 11 !! 'key_vvl' variable volume 12 12 !!---------------------------------------------------------------------- 13 !! dom_vvl : defined coefficients to distribute ssh on each layers14 13 !!---------------------------------------------------------------------- 14 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 15 !! dom_vvl_sf_nxt : Compute next vertical scale factors 16 !! dom_vvl_sf_swp : Swap vertical scale factors and update the vertical grid 17 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 18 !! dom_vvl_rst : read/write restart file 19 !! dom_vvl_ctl : Check the vvl options 20 !! dom_vvl_orca_fix : Recompute some area-weighted interpolations of vertical scale factors 21 !! : to account for manual changes to e[1,2][u,v] in some Straits 22 !!---------------------------------------------------------------------- 23 !! * Modules used 15 24 USE oce ! ocean dynamics and tracers 16 25 USE dom_oce ! ocean space and time domain 17 USE sbc_oce ! surface boundary condition: ocean 18 USE phycst ! physical constants 26 USE sbc_oce ! ocean surface boundary condition 19 27 USE in_out_manager ! I/O manager 28 USE iom ! I/O manager library 29 USE restart ! ocean restart 20 30 USE lib_mpp ! distributed memory computing library 21 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 26 36 PRIVATE 27 37 28 PUBLIC dom_vvl ! called by domain.F90 29 PUBLIC dom_vvl_2 ! called by domain.F90 30 PUBLIC dom_vvl_alloc ! called by nemogcm.F90 31 32 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mut , muu , muv , muf !: 1/H_0 at t-,u-,v-,f-points 33 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: r2dt ! vertical profile time-step, = 2 rdttra 35 ! ! except at nit000 (=rdttra) if neuler=0 38 !! * Routine accessibility 39 PUBLIC dom_vvl_init ! called by domain.F90 40 PUBLIC dom_vvl_sf_nxt ! called by step.F90 41 PUBLIC dom_vvl_sf_swp ! called by step.F90 42 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 43 PRIVATE dom_vvl_orca_fix ! called by dom_vvl_interpol 44 45 !!* Namelist nam_vvl 46 LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate 47 LOGICAL , PUBLIC :: ln_vvl_ztilde = .FALSE. ! ztilde vertical coordinate 48 LOGICAL , PUBLIC :: ln_vvl_layer = .FALSE. ! level vertical coordinate 49 LOGICAL , PUBLIC :: ln_vvl_ztilde_as_zstar = .FALSE. ! ztilde vertical coordinate 50 LOGICAL , PUBLIC :: ln_vvl_zstar_at_eqtor = .FALSE. ! ztilde vertical coordinate 51 LOGICAL , PUBLIC :: ln_vvl_kepe = .FALSE. ! kinetic/potential energy transfer 52 ! ! conservation: not used yet 53 REAL(wp) :: rn_ahe3 = 0.0_wp ! thickness diffusion coefficient 54 REAL(wp) :: rn_rst_e3t = 30._wp ! ztilde to zstar restoration timescale [days] 55 REAL(wp) :: rn_lf_cutoff = 5.0_wp ! cutoff frequency for low-pass filter [days] 56 REAL(wp) :: rn_zdef_max = 0.9_wp ! maximum fractional e3t deformation 57 LOGICAL , PUBLIC :: ln_vvl_dbg = .FALSE. ! debug control prints 58 59 !! * Module variables 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: un_td, vn_td ! thickness diffusion transport 61 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hdiv_lf ! low frequency part of hz divergence 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_b, tilde_e3t_n ! baroclinic scale factors 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: tilde_e3t_a ! baroclinic scale factors 64 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_e3t ! retoring period for scale factors 65 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: frq_rst_hdv ! retoring period for low freq. divergence 36 66 37 67 !! * Substitutions … … 39 69 # include "vectopt_loop_substitute.h90" 40 70 !!---------------------------------------------------------------------- 41 !! NEMO/OPA 4.0 , NEMO Consortium (2011)71 !! NEMO/OPA 3.3 , NEMO-Consortium (2010) 42 72 !! $Id$ 43 73 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 44 74 !!---------------------------------------------------------------------- 45 CONTAINS 75 76 CONTAINS 46 77 47 78 INTEGER FUNCTION dom_vvl_alloc() 48 79 !!---------------------------------------------------------------------- 49 !! *** ROUTINE dom_vvl_alloc *** 50 !!---------------------------------------------------------------------- 80 !! *** FUNCTION dom_vvl_alloc *** 81 !!---------------------------------------------------------------------- 82 IF( ln_vvl_zstar ) dom_vvl_alloc = 0 83 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 84 ALLOCATE( tilde_e3t_b(jpi,jpj,jpk) , tilde_e3t_n(jpi,jpj,jpk) , tilde_e3t_a(jpi,jpj,jpk) , & 85 & un_td (jpi,jpj,jpk) , vn_td (jpi,jpj,jpk) , STAT = dom_vvl_alloc ) 86 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 87 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 88 un_td = 0.0_wp 89 vn_td = 0.0_wp 90 ENDIF 91 IF( ln_vvl_ztilde ) THEN 92 ALLOCATE( frq_rst_e3t(jpi,jpj) , frq_rst_hdv(jpi,jpj) , hdiv_lf(jpi,jpj,jpk) , STAT= dom_vvl_alloc ) 93 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 94 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 95 ENDIF 96 97 END FUNCTION dom_vvl_alloc 98 99 100 SUBROUTINE dom_vvl_init 101 !!---------------------------------------------------------------------- 102 !! *** ROUTINE dom_vvl_init *** 103 !! 104 !! ** Purpose : Initialization of all scale factors, depths 105 !! and water column heights 106 !! 107 !! ** Method : - use restart file and/or initialize 108 !! - interpolate scale factors 109 !! 110 !! ** Action : - fse3t_(n/b) and tilde_e3t_(n/b) 111 !! - Regrid: fse3(u/v)_n 112 !! fse3(u/v)_b 113 !! fse3w_n 114 !! fse3(u/v)w_b 115 !! fse3(u/v)w_n 116 !! fsdept_n, fsdepw_n and fsde3w_n 117 !! - h(t/u/v)_0 118 !! - frq_rst_e3t and frq_rst_hdv 119 !! 120 !! Reference : Leclair, M., and G. Madec, 2011, Ocean Modelling. 121 !!---------------------------------------------------------------------- 122 USE phycst, ONLY : rpi, rsmall, rad 123 !! * Local declarations 124 INTEGER :: ji,jj,jk 125 INTEGER :: ii0, ii1, ij0, ij1 126 !!---------------------------------------------------------------------- 127 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_init') 128 129 IF(lwp) WRITE(numout,*) 130 IF(lwp) WRITE(numout,*) 'dom_vvl_init : Variable volume activated' 131 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~' 132 133 ! choose vertical coordinate (z_star, z_tilde or layer) 134 ! ========================== 135 CALL dom_vvl_ctl 136 137 ! Allocate module arrays 138 ! ====================== 139 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl_init : unable to allocate arrays' ) 140 141 ! Read or initialize fse3t_(b/n), tilde_e3t_(b/n) and hdiv_lf (and e3t_a(jpk)) 142 ! ============================================================================ 143 CALL dom_vvl_rst( nit000, 'READ' ) 144 fse3t_a(:,:,jpk) = e3t_0(:,:,jpk) 145 146 ! Reconstruction of all vertical scale factors at now and before time steps 147 ! ============================================================================= 148 ! Horizontal scale factor interpolations 149 ! -------------------------------------- 150 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3u_b(:,:,:), 'U' ) 151 CALL dom_vvl_interpol( fse3t_b(:,:,:), fse3v_b(:,:,:), 'V' ) 152 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3u_n(:,:,:), 'U' ) 153 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3v_n(:,:,:), 'V' ) 154 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n(:,:,:), 'F' ) 155 ! Vertical scale factor interpolations 156 ! ------------------------------------ 157 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 158 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 159 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 160 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 161 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 162 ! t- and w- points depth 163 ! ---------------------- 164 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 165 fsdepw_n(:,:,1) = 0.0_wp 166 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 167 DO jk = 2, jpk 168 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 169 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 170 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 171 END DO 172 ! Reference water column height at t-, u- and v- point 173 ! ---------------------------------------------------- 174 ht_0(:,:) = 0.0_wp 175 hu_0(:,:) = 0.0_wp 176 hv_0(:,:) = 0.0_wp 177 DO jk = 1, jpk 178 ht_0(:,:) = ht_0(:,:) + e3t_0(:,:,jk) * tmask(:,:,jk) 179 hu_0(:,:) = hu_0(:,:) + e3u_0(:,:,jk) * umask(:,:,jk) 180 hv_0(:,:) = hv_0(:,:) + e3v_0(:,:,jk) * vmask(:,:,jk) 181 END DO 182 183 ! Restoring frequencies for z_tilde coordinate 184 ! ============================================ 185 IF( ln_vvl_ztilde ) THEN 186 ! Values in days provided via the namelist; use rsmall to avoid possible division by zero errors with faulty settings 187 frq_rst_e3t(:,:) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 188 frq_rst_hdv(:,:) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 189 IF( ln_vvl_ztilde_as_zstar ) THEN 190 ! Ignore namelist settings and use these next two to emulate z-star using z-tilde 191 frq_rst_e3t(:,:) = 0.0_wp 192 frq_rst_hdv(:,:) = 1.0_wp / rdt 193 ENDIF 194 IF ( ln_vvl_zstar_at_eqtor ) THEN 195 DO jj = 1, jpj 196 DO ji = 1, jpi 197 IF( ABS(gphit(ji,jj)) >= 6.) THEN 198 ! values outside the equatorial band and transition zone (ztilde) 199 frq_rst_e3t(ji,jj) = 2.0_wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.e0_wp ) 200 frq_rst_hdv(ji,jj) = 2.0_wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.e0_wp ) 201 ELSEIF( ABS(gphit(ji,jj)) <= 2.5) THEN 202 ! values inside the equatorial band (ztilde as zstar) 203 frq_rst_e3t(ji,jj) = 0.0_wp 204 frq_rst_hdv(ji,jj) = 1.0_wp / rdt 205 ELSE 206 ! values in the transition band (linearly vary from ztilde to ztilde as zstar values) 207 frq_rst_e3t(ji,jj) = 0.0_wp + (frq_rst_e3t(ji,jj)-0.0_wp)*0.5_wp & 208 & * ( 1.0_wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 209 & * 180._wp / 3.5_wp ) ) 210 frq_rst_hdv(ji,jj) = (1.0_wp / rdt) & 211 & + ( frq_rst_hdv(ji,jj)-(1.e0_wp / rdt) )*0.5_wp & 212 & * ( 1._wp - COS( rad*(ABS(gphit(ji,jj))-2.5_wp) & 213 & * 180._wp / 3.5_wp ) ) 214 ENDIF 215 END DO 216 END DO 217 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN 218 ii0 = 103 ; ii1 = 111 ! Suppress ztilde in the Foxe Basin for ORCA2 219 ij0 = 128 ; ij1 = 135 ; 220 frq_rst_e3t( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.0_wp 221 frq_rst_hdv( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 1.e0_wp / rdt 222 ENDIF 223 ENDIF 224 ENDIF 225 226 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_init') 227 228 END SUBROUTINE dom_vvl_init 229 230 231 SUBROUTINE dom_vvl_sf_nxt( kt ) 232 !!---------------------------------------------------------------------- 233 !! *** ROUTINE dom_vvl_sf_nxt *** 234 !! 235 !! ** Purpose : - compute the after scale factors used in tra_zdf, dynnxt, 236 !! tranxt and dynspg routines 237 !! 238 !! ** Method : - z_star case: Repartition of ssh INCREMENT proportionnaly to the level thickness. 239 !! - z_tilde_case: after scale factor increment = 240 !! high frequency part of horizontal divergence 241 !! + retsoring towards the background grid 242 !! + thickness difusion 243 !! Then repartition of ssh INCREMENT proportionnaly 244 !! to the "baroclinic" level thickness. 245 !! 246 !! ** Action : - hdiv_lf : restoring towards full baroclinic divergence in z_tilde case 247 !! - tilde_e3t_a: after increment of vertical scale factor 248 !! in z_tilde case 249 !! - fse3(t/u/v)_a 250 !! 251 !! Reference : Leclair, M., and Madec, G. 2011, Ocean Modelling. 252 !!---------------------------------------------------------------------- 253 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze3t 254 REAL(wp), POINTER, DIMENSION(:,: ) :: zht, z_scale, zwu, zwv, zhdiv 255 !! * Arguments 256 INTEGER, INTENT( in ) :: kt ! time step 257 !! * Local declarations 258 INTEGER :: ji, jj, jk ! dummy loop indices 259 INTEGER , DIMENSION(3) :: ijk_max, ijk_min ! temporary integers 260 REAL(wp) :: z2dt ! temporary scalars 261 REAL(wp) :: z_tmin, z_tmax ! temporary scalars 262 !!---------------------------------------------------------------------- 263 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 264 CALL wrk_alloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 265 CALL wrk_alloc( jpi, jpj, jpk, ze3t ) 266 267 IF(kt == nit000) THEN 268 IF(lwp) WRITE(numout,*) 269 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_nxt : compute after scale factors' 270 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~' 271 ENDIF 272 273 ! ******************************* ! 274 ! After acale factors at t-points ! 275 ! ******************************* ! 276 277 ! ! ----------------- ! 278 IF( ln_vvl_zstar ) THEN ! z_star coordinate ! 279 ! ! ----------------- ! 280 281 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * tmask(:,:,1) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 282 DO jk = 1, jpkm1 283 fse3t_a(:,:,jk) = fse3t_b(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 284 END DO 285 286 ! ! --------------------------- ! 287 ELSEIF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 288 ! ! --------------------------- ! 289 290 ! I - initialization 291 ! ================== 292 293 ! 1 - barotropic divergence 294 ! ------------------------- 295 zhdiv(:,:) = 0. 296 zht(:,:) = 0. 297 DO jk = 1, jpkm1 298 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk) 299 zht (:,:) = zht (:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 300 END DO 301 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask(:,:,1) ) 302 303 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) 304 ! -------------------------------------------------- 305 IF( ln_vvl_ztilde ) THEN 306 IF( kt .GT. nit000 ) THEN 307 DO jk = 1, jpkm1 308 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 309 & * ( hdiv_lf(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 310 END DO 311 ENDIF 312 END IF 313 314 ! II - after z_tilde increments of vertical scale factors 315 ! ======================================================= 316 tilde_e3t_a(:,:,:) = 0.0_wp ! tilde_e3t_a used to store tendency terms 317 318 ! 1 - High frequency divergence term 319 ! ---------------------------------- 320 IF( ln_vvl_ztilde ) THEN ! z_tilde case 321 DO jk = 1, jpkm1 322 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 323 END DO 324 ELSE ! layer case 325 DO jk = 1, jpkm1 326 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - fse3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) 327 END DO 328 END IF 329 330 ! 2 - Restoring term (z-tilde case only) 331 ! ------------------ 332 IF( ln_vvl_ztilde ) THEN 333 DO jk = 1, jpk 334 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 335 END DO 336 END IF 337 338 ! 3 - Thickness diffusion term 339 ! ---------------------------- 340 zwu(:,:) = 0.0_wp 341 zwv(:,:) = 0.0_wp 342 ! a - first derivative: diffusive fluxes 343 DO jk = 1, jpkm1 344 DO jj = 1, jpjm1 345 DO ji = 1, fs_jpim1 ! vector opt. 346 un_td(ji,jj,jk) = rn_ahe3 * umask(ji,jj,jk) * re2u_e1u(ji,jj) & 347 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji+1,jj ,jk) ) 348 vn_td(ji,jj,jk) = rn_ahe3 * vmask(ji,jj,jk) * re1v_e2v(ji,jj) & 349 & * ( tilde_e3t_b(ji,jj,jk) - tilde_e3t_b(ji ,jj+1,jk) ) 350 zwu(ji,jj) = zwu(ji,jj) + un_td(ji,jj,jk) 351 zwv(ji,jj) = zwv(ji,jj) + vn_td(ji,jj,jk) 352 END DO 353 END DO 354 END DO 355 ! b - correction for last oceanic u-v points 356 DO jj = 1, jpj 357 DO ji = 1, jpi 358 un_td(ji,jj,mbku(ji,jj)) = un_td(ji,jj,mbku(ji,jj)) - zwu(ji,jj) 359 vn_td(ji,jj,mbkv(ji,jj)) = vn_td(ji,jj,mbkv(ji,jj)) - zwv(ji,jj) 360 END DO 361 END DO 362 ! c - second derivative: divergence of diffusive fluxes 363 DO jk = 1, jpkm1 364 DO jj = 2, jpjm1 365 DO ji = fs_2, fs_jpim1 ! vector opt. 366 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) + ( un_td(ji-1,jj ,jk) - un_td(ji,jj,jk) & 367 & + vn_td(ji ,jj-1,jk) - vn_td(ji,jj,jk) & 368 & ) * r1_e12t(ji,jj) 369 END DO 370 END DO 371 END DO 372 ! d - thickness diffusion transport: boundary conditions 373 ! (stored for tracer advction and continuity equation) 374 CALL lbc_lnk( un_td , 'U' , -1.) 375 CALL lbc_lnk( vn_td , 'V' , -1.) 376 377 ! 4 - Time stepping of baroclinic scale factors 378 ! --------------------------------------------- 379 ! Leapfrog time stepping 380 ! ~~~~~~~~~~~~~~~~~~~~~~ 381 IF( neuler == 0 .AND. kt == nit000 ) THEN 382 z2dt = rdt 383 ELSE 384 z2dt = 2.0_wp * rdt 385 ENDIF 386 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1. ) 387 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 388 389 ! Maximum deformation control 390 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 391 ze3t(:,:,jpk) = 0.0_wp 392 DO jk = 1, jpkm1 393 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 394 END DO 395 z_tmax = MAXVAL( ze3t(:,:,:) ) 396 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 397 z_tmin = MINVAL( ze3t(:,:,:) ) 398 IF( lk_mpp ) CALL mpp_min( z_tmin ) ! min over the global domain 399 ! - ML - test: for the moment, stop simulation for too large e3_t variations 400 IF( ( z_tmax .GT. rn_zdef_max ) .OR. ( z_tmin .LT. - rn_zdef_max ) ) THEN 401 IF( lk_mpp ) THEN 402 CALL mpp_maxloc( ze3t, tmask, z_tmax, ijk_max(1), ijk_max(2), ijk_max(3) ) 403 CALL mpp_minloc( ze3t, tmask, z_tmin, ijk_min(1), ijk_min(2), ijk_min(3) ) 404 ELSE 405 ijk_max = MAXLOC( ze3t(:,:,:) ) 406 ijk_max(1) = ijk_max(1) + nimpp - 1 407 ijk_max(2) = ijk_max(2) + njmpp - 1 408 ijk_min = MINLOC( ze3t(:,:,:) ) 409 ijk_min(1) = ijk_min(1) + nimpp - 1 410 ijk_min(2) = ijk_min(2) + njmpp - 1 411 ENDIF 412 IF (lwp) THEN 413 WRITE(numout, *) 'MAX( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmax 414 WRITE(numout, *) 'at i, j, k=', ijk_max 415 WRITE(numout, *) 'MIN( tilde_e3t_a(:,:,:) / e3t_0(:,:,:) ) =', z_tmin 416 WRITE(numout, *) 'at i, j, k=', ijk_min 417 CALL ctl_warn('MAX( ABS( tilde_e3t_a(:,:,:) ) / e3t_0(:,:,:) ) too high') 418 ENDIF 419 ENDIF 420 ! - ML - end test 421 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 422 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 423 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 424 425 ! Add "tilda" part to the after scale factor 426 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 427 fse3t_a(:,:,:) = e3t_0(:,:,:) + tilde_e3t_a(:,:,:) 428 429 ! III - Barotropic repartition of the sea surface height over the baroclinic profile 430 ! ================================================================================== 431 ! add e3t(n-1) "star" Asselin-filtered 432 DO jk = 1, jpkm1 433 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_b(:,:,jk) - e3t_0(:,:,jk) - tilde_e3t_b(:,:,jk) 434 END DO 435 ! add ( ssh increment + "baroclinicity error" ) proportionnaly to e3t(n) 436 ! - ML - baroclinicity error should be better treated in the future 437 ! i.e. locally and not spread over the water column. 438 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 439 zht(:,:) = 0. 440 DO jk = 1, jpkm1 441 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 442 END DO 443 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) - zht(:,:) ) / ( ht_0(:,:) + sshn(:,:) + 1. - tmask(:,:,1) ) 444 DO jk = 1, jpkm1 445 fse3t_a(:,:,jk) = fse3t_a(:,:,jk) + fse3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 446 END DO 447 448 ENDIF 449 450 IF( ln_vvl_dbg ) THEN ! - ML - test: control prints for debuging 451 ! 452 IF( lwp ) WRITE(numout, *) 'kt =', kt 453 IF ( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 454 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( zht(:,:) ) ) 455 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 456 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(SUM(tilde_e3t_a))) =', z_tmax 457 END IF 458 ! 459 zht(:,:) = 0.0_wp 460 DO jk = 1, jpkm1 461 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 462 END DO 463 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 464 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 465 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(fse3t_n))) =', z_tmax 466 ! 467 zht(:,:) = 0.0_wp 468 DO jk = 1, jpkm1 469 zht(:,:) = zht(:,:) + fse3t_a(:,:,jk) * tmask(:,:,jk) 470 END DO 471 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 472 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 473 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(fse3t_a))) =', z_tmax 474 ! 475 zht(:,:) = 0.0_wp 476 DO jk = 1, jpkm1 477 zht(:,:) = zht(:,:) + fse3t_b(:,:,jk) * tmask(:,:,jk) 478 END DO 479 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 480 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 481 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshb-SUM(fse3t_b))) =', z_tmax 482 ! 483 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshb(:,:) ) ) 484 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 485 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshb))) =', z_tmax 486 ! 487 z_tmax = MAXVAL( tmask(:,:,1) * ABS( sshn(:,:) ) ) 488 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 489 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(sshn))) =', z_tmax 490 ! 491 z_tmax = MAXVAL( tmask(:,:,1) * ABS( ssha(:,:) ) ) 492 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 493 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ssha))) =', z_tmax 494 END IF 495 496 ! *********************************** ! 497 ! After scale factors at u- v- points ! 498 ! *********************************** ! 499 500 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3u_a(:,:,:), 'U' ) 501 CALL dom_vvl_interpol( fse3t_a(:,:,:), fse3v_a(:,:,:), 'V' ) 502 503 CALL wrk_dealloc( jpi, jpj, zht, z_scale, zwu, zwv, zhdiv ) 504 CALL wrk_dealloc( jpi, jpj, jpk, ze3t ) 505 506 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_nxt') 507 508 END SUBROUTINE dom_vvl_sf_nxt 509 510 511 SUBROUTINE dom_vvl_sf_swp( kt ) 512 !!---------------------------------------------------------------------- 513 !! *** ROUTINE dom_vvl_sf_swp *** 514 !! 515 !! ** Purpose : compute time filter and swap of scale factors 516 !! compute all depths and related variables for next time step 517 !! write outputs and restart file 518 !! 519 !! ** Method : - swap of e3t with trick for volume/tracer conservation 520 !! - reconstruct scale factor at other grid points (interpolate) 521 !! - recompute depths and water height fields 522 !! 523 !! ** Action : - fse3t_(b/n), tilde_e3t_(b/n) and fse3(u/v)_n ready for next time step 524 !! - Recompute: 525 !! fse3(u/v)_b 526 !! fse3w_n 527 !! fse3(u/v)w_b 528 !! fse3(u/v)w_n 529 !! fsdept_n, fsdepw_n and fsde3w_n 530 !! h(u/v) and h(u/v)r 531 !! 532 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 533 !! Leclair, M., and G. Madec, 2011, Ocean Modelling. 534 !!---------------------------------------------------------------------- 535 !! * Arguments 536 INTEGER, INTENT( in ) :: kt ! time step 537 !! * Local declarations 538 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_e3t_def 539 INTEGER :: jk ! dummy loop indices 540 !!---------------------------------------------------------------------- 541 542 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_sf_swp') 51 543 ! 52 ALLOCATE( mut (jpi,jpj,jpk) , muu (jpi,jpj,jpk) , muv (jpi,jpj,jpk) , muf (jpi,jpj,jpk) , & 53 & r2dt (jpk) , STAT=dom_vvl_alloc ) 54 ! 55 IF( lk_mpp ) CALL mpp_sum ( dom_vvl_alloc ) 56 IF( dom_vvl_alloc /= 0 ) CALL ctl_warn('dom_vvl_alloc: failed to allocate arrays') 544 CALL wrk_alloc( jpi, jpj, jpk, z_e3t_def ) 57 545 ! 58 END FUNCTION dom_vvl_alloc 59 60 61 SUBROUTINE dom_vvl 62 !!---------------------------------------------------------------------- 63 !! *** ROUTINE dom_vvl *** 64 !! 65 !! ** Purpose : compute mu coefficients at t-, u-, v- and f-points to 66 !! spread ssh over the whole water column (scale factors) 67 !! set the before and now ssh at u- and v-points 68 !! (also f-point in now case) 69 !!---------------------------------------------------------------------- 546 IF( kt == nit000 ) THEN 547 IF(lwp) WRITE(numout,*) 548 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_swp : - time filter and swap of scale factors' 549 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step' 550 ENDIF 70 551 ! 71 INTEGER :: ji, jj, jk ! dummy loop indices 72 REAL(wp) :: zcoefu, zcoefv , zcoeff ! local scalars 73 REAL(wp) :: zvt , zvt_ip1, zvt_jp1, zvt_ip1jp1 ! - - 74 REAL(wp), POINTER, DIMENSION(:,:) :: zee_t, zee_u, zee_v, zee_f ! 2D workspace 75 !!---------------------------------------------------------------------- 552 ! Time filter and swap of scale factors 553 ! ===================================== 554 ! - ML - fse3(t/u/v)_b are allready computed in dynnxt. 555 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 556 IF( neuler == 0 .AND. kt == nit000 ) THEN 557 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 558 ELSE 559 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 560 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 561 ENDIF 562 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 563 ENDIF 564 fse3t_n(:,:,:) = fse3t_a(:,:,:) 565 fse3u_n(:,:,:) = fse3u_a(:,:,:) 566 fse3v_n(:,:,:) = fse3v_a(:,:,:) 567 568 ! Compute all missing vertical scale factor and depths 569 ! ==================================================== 570 ! Horizontal scale factor interpolations 571 ! -------------------------------------- 572 ! - ML - fse3u_b and fse3v_b are allready computed in dynnxt 573 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3f_n (:,:,:), 'F' ) 574 ! Vertical scale factor interpolations 575 ! ------------------------------------ 576 CALL dom_vvl_interpol( fse3t_n(:,:,:), fse3w_n (:,:,:), 'W' ) 577 CALL dom_vvl_interpol( fse3u_n(:,:,:), fse3uw_n(:,:,:), 'UW' ) 578 CALL dom_vvl_interpol( fse3v_n(:,:,:), fse3vw_n(:,:,:), 'VW' ) 579 CALL dom_vvl_interpol( fse3u_b(:,:,:), fse3uw_b(:,:,:), 'UW' ) 580 CALL dom_vvl_interpol( fse3v_b(:,:,:), fse3vw_b(:,:,:), 'VW' ) 581 ! t- and w- points depth 582 ! ---------------------- 583 fsdept_n(:,:,1) = 0.5_wp * fse3w_n(:,:,1) 584 fsdepw_n(:,:,1) = 0.0_wp 585 fsde3w_n(:,:,1) = fsdept_n(:,:,1) - sshn(:,:) 586 DO jk = 2, jpk 587 fsdept_n(:,:,jk) = fsdept_n(:,:,jk-1) + fse3w_n(:,:,jk) 588 fsdepw_n(:,:,jk) = fsdepw_n(:,:,jk-1) + fse3t_n(:,:,jk-1) 589 fsde3w_n(:,:,jk) = fsdept_n(:,:,jk ) - sshn (:,:) 590 END DO 591 ! Local depth and Inverse of the local depth of the water column at u- and v- points 592 ! ---------------------------------------------------------------------------------- 593 hu(:,:) = 0. 594 hv(:,:) = 0. 595 DO jk = 1, jpk 596 hu(:,:) = hu(:,:) + fse3u_n(:,:,jk) * umask(:,:,jk) 597 hv(:,:) = hv(:,:) + fse3v_n(:,:,jk) * vmask(:,:,jk) 598 END DO 599 ! Inverse of the local depth 600 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1. - umask(:,:,1) ) 601 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1. - vmask(:,:,1) ) 602 603 ! Write outputs 604 ! ============= 605 z_e3t_def(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 606 CALL iom_put( "e3t_n" , fse3t_n (:,:,:) ) 607 CALL iom_put( "dept_n" , fsde3w_n (:,:,:) ) 608 CALL iom_put( "e3tdef" , z_e3t_def(:,:,:) ) 609 610 ! write restart file 611 ! ================== 612 IF( lrst_oce ) CALL dom_vvl_rst( kt, 'WRITE' ) 76 613 ! 77 IF( nn_timing == 1 ) CALL timing_start('dom_vvl')614 CALL wrk_dealloc( jpi, jpj, jpk, z_e3t_def ) 78 615 ! 79 CALL wrk_alloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 616 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_sf_swp') 617 618 END SUBROUTINE dom_vvl_sf_swp 619 620 621 SUBROUTINE dom_vvl_interpol( pe3_in, pe3_out, pout ) 622 !!--------------------------------------------------------------------- 623 !! *** ROUTINE dom_vvl__interpol *** 624 !! 625 !! ** Purpose : interpolate scale factors from one grid point to another 626 !! 627 !! ** Method : e3_out = e3_0 + interpolation(e3_in - e3_0) 628 !! - horizontal interpolation: grid cell surface averaging 629 !! - vertical interpolation: simple averaging 630 !!---------------------------------------------------------------------- 631 !! * Arguments 632 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 633 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 634 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors 635 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 636 !! * Local declarations 637 INTEGER :: ji, jj, jk ! dummy loop indices 638 LOGICAL :: l_is_orca ! local logical 639 !!---------------------------------------------------------------------- 640 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_interpol') 641 ! 642 l_is_orca = .FALSE. 643 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) l_is_orca = .TRUE. ! ORCA R2 configuration - will need to correct some locations 644 645 SELECT CASE ( pout ) 646 ! ! ------------------------------------- ! 647 CASE( 'U' ) ! interpolation from T-point to U-point ! 648 ! ! ------------------------------------- ! 649 ! horizontal surface weighted interpolation 650 DO jk = 1, jpk 651 DO jj = 1, jpjm1 652 DO ji = 1, fs_jpim1 ! vector opt. 653 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * r1_e12u(ji,jj) & 654 & * ( e12t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 655 & + e12t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) ) 656 END DO 657 END DO 658 END DO 659 ! 660 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 661 ! boundary conditions 662 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1. ) 663 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 664 ! ! ------------------------------------- ! 665 CASE( 'V' ) ! interpolation from T-point to V-point ! 666 ! ! ------------------------------------- ! 667 ! horizontal surface weighted interpolation 668 DO jk = 1, jpk 669 DO jj = 1, jpjm1 670 DO ji = 1, fs_jpim1 ! vector opt. 671 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) * r1_e12v(ji,jj) & 672 & * ( e12t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 673 & + e12t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) ) 674 END DO 675 END DO 676 END DO 677 ! 678 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 679 ! boundary conditions 680 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1. ) 681 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 682 ! ! ------------------------------------- ! 683 CASE( 'F' ) ! interpolation from U-point to F-point ! 684 ! ! ------------------------------------- ! 685 ! horizontal surface weighted interpolation 686 DO jk = 1, jpk 687 DO jj = 1, jpjm1 688 DO ji = 1, fs_jpim1 ! vector opt. 689 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) * r1_e12f(ji,jj) & 690 & * ( e12u(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3u_0(ji,jj ,jk) ) & 691 & + e12u(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3u_0(ji,jj+1,jk) ) ) 692 END DO 693 END DO 694 END DO 695 ! 696 IF( l_is_orca ) CALL dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 697 ! boundary conditions 698 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1. ) 699 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 700 ! ! ------------------------------------- ! 701 CASE( 'W' ) ! interpolation from T-point to W-point ! 702 ! ! ------------------------------------- ! 703 ! vertical simple interpolation 704 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 705 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 706 DO jk = 2, jpk 707 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * tmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 708 & + 0.5_wp * tmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 709 END DO 710 ! ! -------------------------------------- ! 711 CASE( 'UW' ) ! interpolation from U-point to UW-point ! 712 ! ! -------------------------------------- ! 713 ! vertical simple interpolation 714 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 715 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 716 DO jk = 2, jpk 717 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * umask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 718 & + 0.5_wp * umask(:,:,jk) * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 719 END DO 720 ! ! -------------------------------------- ! 721 CASE( 'VW' ) ! interpolation from V-point to VW-point ! 722 ! ! -------------------------------------- ! 723 ! vertical simple interpolation 724 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 725 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 726 DO jk = 2, jpk 727 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * vmask(:,:,jk) ) * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 728 & + 0.5_wp * vmask(:,:,jk) * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 729 END DO 730 END SELECT 80 731 ! 81 IF(lwp) THEN 732 733 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_interpol') 734 735 END SUBROUTINE dom_vvl_interpol 736 737 SUBROUTINE dom_vvl_rst( kt, cdrw ) 738 !!--------------------------------------------------------------------- 739 !! *** ROUTINE dom_vvl_rst *** 740 !! 741 !! ** Purpose : Read or write VVL file in restart file 742 !! 743 !! ** Method : use of IOM library 744 !! if the restart does not contain vertical scale factors, 745 !! they are set to the _0 values 746 !! if the restart does not contain vertical scale factors increments (z_tilde), 747 !! they are set to 0. 748 !!---------------------------------------------------------------------- 749 !! * Arguments 750 INTEGER , INTENT(in) :: kt ! ocean time-step 751 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 752 !! * Local declarations 753 INTEGER :: id1, id2, id3, id4, id5 ! local integers 754 !!---------------------------------------------------------------------- 755 ! 756 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_rst') 757 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 758 ! ! =============== 759 IF( ln_rstart ) THEN !* Read the restart file 760 CALL rst_read_open ! open the restart file if necessary 761 id1 = iom_varid( numror, 'fse3t_b', ldstop = .FALSE. ) 762 id2 = iom_varid( numror, 'fse3t_n', ldstop = .FALSE. ) 763 id3 = iom_varid( numror, 'tilde_e3t_b', ldstop = .FALSE. ) 764 id4 = iom_varid( numror, 'tilde_e3t_n', ldstop = .FALSE. ) 765 id5 = iom_varid( numror, 'hdif_lf', ldstop = .FALSE. ) 766 ! ! --------- ! 767 ! ! all cases ! 768 ! ! --------- ! 769 IF( MIN( id1, id2 ) > 0 ) THEN ! all required arrays exist 770 CALL iom_get( numror, jpdom_autoglo, 'fse3t_b', fse3t_b(:,:,:) ) 771 CALL iom_get( numror, jpdom_autoglo, 'fse3t_n', fse3t_n(:,:,:) ) 772 IF( neuler == 0 ) THEN 773 fse3t_b(:,:,:) = fse3t_n(:,:,:) 774 ENDIF 775 ELSE IF( id1 > 0 ) THEN 776 IF(lwp) write(numout,*) 'dom_vvl_rst WARNING : fse3t_n not found in restart files' 777 IF(lwp) write(numout,*) 'fse3t_n set equal to fse3t_b.' 778 fse3t_b(:,:,:) = fse3t_n(:,:,:) 779 ELSE ! one at least array is missing 780 CALL ctl_stop( 'dom_vvl_rst: vvl cannot restart from a non vvl run' ) 781 ENDIF 782 ! ! ----------- ! 783 IF( ln_vvl_zstar ) THEN ! z_star case ! 784 ! ! ----------- ! 785 IF( MIN( id3, id4 ) > 0 ) THEN 786 CALL ctl_stop( 'dom_vvl_rst: z_star cannot restart from a z_tilde or layer run' ) 787 ENDIF 788 ! ! ----------------------- ! 789 ELSE ! z_tilde and layer cases ! 790 ! ! ----------------------- ! 791 IF( MIN( id3, id4 ) > 0 ) THEN ! all required arrays exist 792 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 793 CALL iom_get( numror, jpdom_autoglo, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 794 ELSE ! one at least array is missing 795 tilde_e3t_b(:,:,:) = 0.0_wp 796 tilde_e3t_n(:,:,:) = 0.0_wp 797 ENDIF 798 ! ! ------------ ! 799 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 800 ! ! ------------ ! 801 IF( id5 > 0 ) THEN ! required array exists 802 CALL iom_get( numror, jpdom_autoglo, 'hdiv_lf', hdiv_lf(:,:,:) ) 803 ELSE ! array is missing 804 hdiv_lf(:,:,:) = 0.0_wp 805 ENDIF 806 ENDIF 807 ENDIF 808 ! 809 ELSE !* Initialize at "rest" 810 fse3t_b(:,:,:) = e3t_0(:,:,:) 811 fse3t_n(:,:,:) = e3t_0(:,:,:) 812 IF( ln_vvl_ztilde .OR. ln_vvl_layer) THEN 813 tilde_e3t_b(:,:,:) = 0.0_wp 814 tilde_e3t_n(:,:,:) = 0.0_wp 815 IF( ln_vvl_ztilde ) hdiv_lf(:,:,:) = 0.0_wp 816 END IF 817 ENDIF 818 819 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 820 ! ! =================== 821 IF(lwp) WRITE(numout,*) '---- dom_vvl_rst ----' 822 ! ! --------- ! 823 ! ! all cases ! 824 ! ! --------- ! 825 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_b', fse3t_b(:,:,:) ) 826 CALL iom_rstput( kt, nitrst, numrow, 'fse3t_n', fse3t_n(:,:,:) ) 827 ! ! ----------------------- ! 828 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde and layer cases ! 829 ! ! ----------------------- ! 830 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_b', tilde_e3t_b(:,:,:) ) 831 CALL iom_rstput( kt, nitrst, numrow, 'tilde_e3t_n', tilde_e3t_n(:,:,:) ) 832 END IF 833 ! ! -------------! 834 IF( ln_vvl_ztilde ) THEN ! z_tilde case ! 835 ! ! ------------ ! 836 CALL iom_rstput( kt, nitrst, numrow, 'hdiv_lf', hdiv_lf(:,:,:) ) 837 ENDIF 838 839 ENDIF 840 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_rst') 841 842 END SUBROUTINE dom_vvl_rst 843 844 845 SUBROUTINE dom_vvl_ctl 846 !!--------------------------------------------------------------------- 847 !! *** ROUTINE dom_vvl_ctl *** 848 !! 849 !! ** Purpose : Control the consistency between namelist options 850 !! for vertical coordinate 851 !!---------------------------------------------------------------------- 852 INTEGER :: ioptio 853 854 NAMELIST/nam_vvl/ ln_vvl_zstar, ln_vvl_ztilde, ln_vvl_layer, ln_vvl_ztilde_as_zstar, & 855 & ln_vvl_zstar_at_eqtor , rn_ahe3 , rn_rst_e3t , & 856 & rn_lf_cutoff , rn_zdef_max , ln_vvl_dbg ! not yet implemented: ln_vvl_kepe 857 !!---------------------------------------------------------------------- 858 859 REWIND ( numnam ) ! Read Namelist nam_vvl : vertical coordinate 860 READ ( numnam, nam_vvl ) 861 862 IF(lwp) THEN ! Namelist print 82 863 WRITE(numout,*) 83 WRITE(numout,*) 'dom_vvl : Variable volume initialization' 84 WRITE(numout,*) '~~~~~~~~ compute coef. used to spread ssh over each layers' 864 WRITE(numout,*) 'dom_vvl_ctl : choice/control of the variable vertical coordinate' 865 WRITE(numout,*) '~~~~~~~~~~~' 866 WRITE(numout,*) ' Namelist nam_vvl : chose a vertical coordinate' 867 WRITE(numout,*) ' zstar ln_vvl_zstar = ', ln_vvl_zstar 868 WRITE(numout,*) ' ztilde ln_vvl_ztilde = ', ln_vvl_ztilde 869 WRITE(numout,*) ' layer ln_vvl_layer = ', ln_vvl_layer 870 WRITE(numout,*) ' ztilde as zstar ln_vvl_ztilde_as_zstar = ', ln_vvl_ztilde_as_zstar 871 WRITE(numout,*) ' ztilde near the equator ln_vvl_zstar_at_eqtor = ', ln_vvl_zstar_at_eqtor 872 ! WRITE(numout,*) ' Namelist nam_vvl : chose kinetic-to-potential energy conservation' 873 ! WRITE(numout,*) ' ln_vvl_kepe = ', ln_vvl_kepe 874 WRITE(numout,*) ' Namelist nam_vvl : thickness diffusion coefficient' 875 WRITE(numout,*) ' rn_ahe3 = ', rn_ahe3 876 WRITE(numout,*) ' Namelist nam_vvl : maximum e3t deformation fractional change' 877 WRITE(numout,*) ' rn_zdef_max = ', rn_zdef_max 878 IF( ln_vvl_ztilde_as_zstar ) THEN 879 WRITE(numout,*) ' ztilde running in zstar emulation mode; ' 880 WRITE(numout,*) ' ignoring namelist timescale parameters and using:' 881 WRITE(numout,*) ' hard-wired : z-tilde to zstar restoration timescale (days)' 882 WRITE(numout,*) ' rn_rst_e3t = 0.0' 883 WRITE(numout,*) ' hard-wired : z-tilde cutoff frequency of low-pass filter (days)' 884 WRITE(numout,*) ' rn_lf_cutoff = 1.0/rdt' 885 ELSE 886 WRITE(numout,*) ' Namelist nam_vvl : z-tilde to zstar restoration timescale (days)' 887 WRITE(numout,*) ' rn_rst_e3t = ', rn_rst_e3t 888 WRITE(numout,*) ' Namelist nam_vvl : z-tilde cutoff frequency of low-pass filter (days)' 889 WRITE(numout,*) ' rn_lf_cutoff = ', rn_lf_cutoff 890 ENDIF 891 WRITE(numout,*) ' Namelist nam_vvl : debug prints' 892 WRITE(numout,*) ' ln_vvl_dbg = ', ln_vvl_dbg 85 893 ENDIF 86 87 IF( dom_vvl_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'dom_vvl : unable to allocate arrays' ) 88 89 fsdept(:,:,:) = gdept (:,:,:) 90 fsdepw(:,:,:) = gdepw (:,:,:) 91 fsde3w(:,:,:) = gdep3w(:,:,:) 92 fse3t (:,:,:) = e3t (:,:,:) 93 fse3u (:,:,:) = e3u (:,:,:) 94 fse3v (:,:,:) = e3v (:,:,:) 95 fse3f (:,:,:) = e3f (:,:,:) 96 fse3w (:,:,:) = e3w (:,:,:) 97 fse3uw(:,:,:) = e3uw (:,:,:) 98 fse3vw(:,:,:) = e3vw (:,:,:) 99 100 ! !== mu computation ==! 101 zee_t(:,:) = fse3t_0(:,:,1) ! Lower bound : thickness of the first model level 102 zee_u(:,:) = fse3u_0(:,:,1) 103 zee_v(:,:) = fse3v_0(:,:,1) 104 zee_f(:,:) = fse3f_0(:,:,1) 105 DO jk = 2, jpkm1 ! Sum of the masked vertical scale factors 106 zee_t(:,:) = zee_t(:,:) + fse3t_0(:,:,jk) * tmask(:,:,jk) 107 zee_u(:,:) = zee_u(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 108 zee_v(:,:) = zee_v(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 109 DO jj = 1, jpjm1 ! f-point : fmask=shlat at coasts, use the product of umask 110 zee_f(:,jj) = zee_f(:,jj) + fse3f_0(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 111 END DO 112 END DO 113 ! ! Compute and mask the inverse of the local depth at T, U, V and F points 114 zee_t(:,:) = 1._wp / zee_t(:,:) * tmask(:,:,1) 115 zee_u(:,:) = 1._wp / zee_u(:,:) * umask(:,:,1) 116 zee_v(:,:) = 1._wp / zee_v(:,:) * vmask(:,:,1) 117 DO jj = 1, jpjm1 ! f-point case fmask cannot be used 118 zee_f(:,jj) = 1._wp / zee_f(:,jj) * umask(:,jj,1) * umask(:,jj+1,1) 119 END DO 120 CALL lbc_lnk( zee_f, 'F', 1. ) ! lateral boundary condition on ee_f 121 ! 122 DO jk = 1, jpk ! mu coefficients 123 mut(:,:,jk) = zee_t(:,:) * tmask(:,:,jk) ! T-point at T levels 124 muu(:,:,jk) = zee_u(:,:) * umask(:,:,jk) ! U-point at T levels 125 muv(:,:,jk) = zee_v(:,:) * vmask(:,:,jk) ! V-point at T levels 126 END DO 127 DO jk = 1, jpk ! F-point : fmask=shlat at coasts, use the product of umask 128 DO jj = 1, jpjm1 129 muf(:,jj,jk) = zee_f(:,jj) * umask(:,jj,jk) * umask(:,jj+1,jk) ! at T levels 130 END DO 131 muf(:,jpj,jk) = 0._wp 132 END DO 133 CALL lbc_lnk( muf, 'F', 1. ) ! lateral boundary condition 134 135 136 hu_0(:,:) = 0.e0 ! Reference ocean depth at U- and V-points 137 hv_0(:,:) = 0.e0 138 DO jk = 1, jpk 139 hu_0(:,:) = hu_0(:,:) + fse3u_0(:,:,jk) * umask(:,:,jk) 140 hv_0(:,:) = hv_0(:,:) + fse3v_0(:,:,jk) * vmask(:,:,jk) 141 END DO 142 143 DO jj = 1, jpjm1 ! initialise before and now Sea Surface Height at u-, v-, f-points 144 DO ji = 1, jpim1 ! NO vector opt. 145 zcoefu = 0.50_wp / ( e1u(ji,jj) * e2u(ji,jj) ) * umask(ji,jj,1) 146 zcoefv = 0.50_wp / ( e1v(ji,jj) * e2v(ji,jj) ) * vmask(ji,jj,1) 147 zcoeff = 0.25_wp / ( e1f(ji,jj) * e2f(ji,jj) ) * umask(ji,jj,1) * umask(ji,jj+1,1) 148 ! 149 zvt = e1e2t(ji ,jj ) * sshb(ji ,jj ) ! before fields 150 zvt_ip1 = e1e2t(ji+1,jj ) * sshb(ji+1,jj ) 151 zvt_jp1 = e1e2t(ji ,jj+1) * sshb(ji ,jj+1) 152 sshu_b(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 153 sshv_b(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 154 ! 155 zvt = e1e2t(ji ,jj ) * sshn(ji ,jj ) ! now fields 156 zvt_ip1 = e1e2t(ji+1,jj ) * sshn(ji+1,jj ) 157 zvt_jp1 = e1e2t(ji ,jj+1) * sshn(ji ,jj+1) 158 zvt_ip1jp1 = e1e2t(ji+1,jj+1) * sshn(ji+1,jj+1) 159 sshu_n(ji,jj) = zcoefu * ( zvt + zvt_ip1 ) 160 sshv_n(ji,jj) = zcoefv * ( zvt + zvt_jp1 ) 161 sshf_n(ji,jj) = zcoeff * ( zvt + zvt_ip1 + zvt_jp1 + zvt_ip1jp1 ) 162 END DO 163 END DO 164 CALL lbc_lnk( sshu_n, 'U', 1. ) ; CALL lbc_lnk( sshu_b, 'U', 1. ) ! lateral boundary conditions 165 CALL lbc_lnk( sshv_n, 'V', 1. ) ; CALL lbc_lnk( sshv_b, 'V', 1. ) 166 CALL lbc_lnk( sshf_n, 'F', 1. ) 167 ! 168 CALL wrk_dealloc( jpi, jpj, zee_t, zee_u, zee_v, zee_f ) 169 ! 170 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl') 171 ! 172 END SUBROUTINE dom_vvl 173 174 175 SUBROUTINE dom_vvl_2( kt, pe3u_b, pe3v_b ) 176 !!---------------------------------------------------------------------- 177 !! *** ROUTINE dom_vvl_2 *** 178 !! 179 !! ** Purpose : compute the vertical scale factors at u- and v-points 180 !! in variable volume case. 181 !! 182 !! ** Method : In variable volume case (non linear sea surface) the 183 !! the vertical scale factor at velocity points is computed 184 !! as the average of the cell surface weighted e3t. 185 !! It uses the sea surface heigth so it have to be initialized 186 !! after ssh is read/set 187 !!---------------------------------------------------------------------- 188 INTEGER , INTENT(in ) :: kt ! ocean time-step index 189 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe3u_b, pe3v_b ! before vertical scale factor at u- & v-pts 190 ! 191 INTEGER :: ji, jj, jk ! dummy loop indices 192 INTEGER :: iku, ikv ! local integers 193 INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers 194 REAL(wp) :: zvt, zvtip1, zvtjp1 ! local scalars 195 !!---------------------------------------------------------------------- 196 ! 197 IF( nn_timing == 1 ) CALL timing_start('dom_vvl_2') 198 ! 199 IF( lwp .AND. kt == nit000 ) THEN 894 895 ioptio = 0 ! Parameter control 896 IF( ln_vvl_ztilde_as_zstar ) ln_vvl_ztilde = .true. 897 IF( ln_vvl_zstar ) ioptio = ioptio + 1 898 IF( ln_vvl_ztilde ) ioptio = ioptio + 1 899 IF( ln_vvl_layer ) ioptio = ioptio + 1 900 901 IF( ioptio /= 1 ) CALL ctl_stop( 'Choose ONE vertical coordinate in namelist nam_vvl' ) 902 903 IF(lwp) THEN ! Print the choice 200 904 WRITE(numout,*) 201 WRITE(numout,*) 'dom_vvl_2 : Variable volume, fse3t_b initialization' 202 WRITE(numout,*) '~~~~~~~~~ ' 203 pe3u_b(:,:,jpk) = fse3u_0(:,:,jpk) 204 pe3v_b(:,:,jpk) = fse3v_0(:,:,jpk) 905 IF( ln_vvl_zstar ) WRITE(numout,*) ' zstar vertical coordinate is used' 906 IF( ln_vvl_ztilde ) WRITE(numout,*) ' ztilde vertical coordinate is used' 907 IF( ln_vvl_layer ) WRITE(numout,*) ' layer vertical coordinate is used' 908 IF( ln_vvl_ztilde_as_zstar ) WRITE(numout,*) ' to emulate a zstar coordinate' 909 ! - ML - Option not developed yet 910 ! IF( ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option used' 911 ! IF( .NOT. ln_vvl_kepe ) WRITE(numout,*) ' kinetic to potential energy transfer : option not used' 205 912 ENDIF 206 207 DO jk = 1, jpkm1 ! set the before scale factors at u- & v-points 208 DO jj = 2, jpjm1 209 DO ji = fs_2, fs_jpim1 210 zvt = ( fse3t_b(ji ,jj ,jk) - fse3t_0(ji ,jj ,jk) ) * e1e2t(ji ,jj ) 211 zvtip1 = ( fse3t_b(ji+1,jj ,jk) - fse3t_0(ji+1,jj ,jk) ) * e1e2t(ji+1,jj ) 212 zvtjp1 = ( fse3t_b(ji ,jj+1,jk) - fse3t_0(ji ,jj+1,jk) ) * e1e2t(ji ,jj+1) 213 pe3u_b(ji,jj,jk) = fse3u_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtip1 ) / ( e1u(ji,jj) * e2u(ji,jj) ) 214 pe3v_b(ji,jj,jk) = fse3v_0(ji,jj,jk) + 0.5_wp * ( zvt + zvtjp1 ) / ( e1v(ji,jj) * e2v(ji,jj) ) 215 END DO 216 END DO 217 END DO 218 219 ! Correct scale factors at locations that have been individually modified in domhgr 220 ! Such modifications break the relationship between e1e2t and e1u*e2u etc. Recompute 221 ! scale factors ignoring the modified metric. 913 914 END SUBROUTINE dom_vvl_ctl 915 916 SUBROUTINE dom_vvl_orca_fix( pe3_in, pe3_out, pout ) 917 !!--------------------------------------------------------------------- 918 !! *** ROUTINE dom_vvl_orca_fix *** 919 !! 920 !! ** Purpose : Correct surface weighted, horizontally interpolated, 921 !! scale factors at locations that have been individually 922 !! modified in domhgr. Such modifications break the 923 !! relationship between e12t and e1u*e2u etc. 924 !! Recompute some scale factors ignoring the modified metric. 925 !!---------------------------------------------------------------------- 926 !! * Arguments 927 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( in ) :: pe3_in ! input e3 to be interpolated 928 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( inout ) :: pe3_out ! output interpolated e3 929 CHARACTER(LEN=*), INTENT( in ) :: pout ! grid point of out scale factors 930 ! ! = 'U', 'V', 'W, 'F', 'UW' or 'VW' 931 !! * Local declarations 932 INTEGER :: ji, jj, jk ! dummy loop indices 933 INTEGER :: ij0, ij1, ii0, ii1 ! dummy loop indices 934 !! acc 935 !! Hmm with the time splitting these "fixes" seem to do more harm than good. Temporarily disabled for 936 !! the ORCA2 tests (by changing jp_cfg test from 2 to 3) pending further investigations 937 !! 222 938 ! ! ===================== 223 939 IF( cp_cfg == "orca" .AND. jp_cfg == 2 ) THEN ! ORCA R2 configuration 224 940 ! ! ===================== 941 !! acc 225 942 IF( nn_cla == 0 ) THEN 226 943 ! 227 944 ii0 = 139 ; ii1 = 140 ! Gibraltar Strait (e2u was modified) 228 ij0 = 102 ; ij1 = 102 229 DO jk = 1, jpkm1 ! set the before scale factors at u-points945 ij0 = 102 ; ij1 = 102 946 DO jk = 1, jpkm1 230 947 DO jj = mj0(ij0), mj1(ij1) 231 948 DO ji = mi0(ii0), mi1(ii1) 232 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 233 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 949 SELECT CASE ( pout ) 950 CASE( 'U' ) 951 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 952 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 953 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 954 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 955 CASE( 'F' ) 956 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 957 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 958 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 959 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 960 END SELECT 234 961 END DO 235 962 END DO … … 237 964 ! 238 965 ii0 = 160 ; ii1 = 160 ! Bab el Mandeb (e2u and e1v were modified) 239 ij0 = 88 ; ij1 = 88 240 DO jk = 1, jpkm1 ! set the before scale factors at u-points966 ij0 = 88 ; ij1 = 88 967 DO jk = 1, jpkm1 241 968 DO jj = mj0(ij0), mj1(ij1) 242 969 DO ji = mi0(ii0), mi1(ii1) 243 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 244 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 970 SELECT CASE ( pout ) 971 CASE( 'U' ) 972 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 973 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 974 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 975 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 976 CASE( 'V' ) 977 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 978 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 979 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 980 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 981 CASE( 'F' ) 982 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 983 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 984 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 985 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 986 END SELECT 245 987 END DO 246 988 END DO 247 989 END DO 248 DO jk = 1, jpkm1 ! set the before scale factors at v-points249 DO jj = mj0(ij0), mj1(ij1)250 DO ji = mi0(ii0), mi1(ii1)251 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj)252 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) )253 END DO254 END DO255 END DO256 990 ENDIF 257 991 258 992 ii0 = 145 ; ii1 = 146 ! Danish Straits (e2u was modified) 259 ij0 = 116 ; ij1 = 116 260 DO jk = 1, jpkm1 ! set the before scale factors at u-points 261 DO jj = mj0(ij0), mj1(ij1) 262 DO ji = mi0(ii0), mi1(ii1) 263 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 264 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 265 END DO 266 END DO 267 END DO 268 ! 993 ij0 = 116 ; ij1 = 116 994 DO jk = 1, jpkm1 995 DO jj = mj0(ij0), mj1(ij1) 996 DO ji = mi0(ii0), mi1(ii1) 997 SELECT CASE ( pout ) 998 CASE( 'U' ) 999 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1000 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1001 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1002 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1003 CASE( 'F' ) 1004 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1005 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1006 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1007 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1008 END SELECT 1009 END DO 1010 END DO 1011 END DO 269 1012 ENDIF 1013 ! 270 1014 ! ! ===================== 271 1015 IF( cp_cfg == "orca" .AND. jp_cfg == 1 ) THEN ! ORCA R1 configuration 272 1016 ! ! ===================== 273 1017 ! 274 1018 ii0 = 281 ; ii1 = 282 ! Gibraltar Strait (e2u was modified) 275 ij0 = 200 ; ij1 = 200 276 DO jk = 1, jpkm1 ! set the before scale factors at u-points 277 DO jj = mj0(ij0), mj1(ij1) 278 DO ji = mi0(ii0), mi1(ii1) 279 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 280 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 281 END DO 282 END DO 283 END DO 284 1019 ij0 = 200 ; ij1 = 200 1020 DO jk = 1, jpkm1 1021 DO jj = mj0(ij0), mj1(ij1) 1022 DO ji = mi0(ii0), mi1(ii1) 1023 SELECT CASE ( pout ) 1024 CASE( 'U' ) 1025 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1026 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1027 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1028 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1029 CASE( 'F' ) 1030 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1031 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1032 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1033 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1034 END SELECT 1035 END DO 1036 END DO 1037 END DO 1038 ! 285 1039 ii0 = 314 ; ii1 = 315 ! Bhosporus Strait (e2u was modified) 286 ij0 = 208 ; ij1 = 208 287 DO jk = 1, jpkm1 ! set the before scale factors at u-points 288 DO jj = mj0(ij0), mj1(ij1) 289 DO ji = mi0(ii0), mi1(ii1) 290 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 291 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 292 END DO 293 END DO 294 END DO 295 1040 ij0 = 208 ; ij1 = 208 1041 DO jk = 1, jpkm1 1042 DO jj = mj0(ij0), mj1(ij1) 1043 DO ji = mi0(ii0), mi1(ii1) 1044 SELECT CASE ( pout ) 1045 CASE( 'U' ) 1046 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1047 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1048 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1049 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1050 CASE( 'F' ) 1051 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1052 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1053 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1054 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1055 END SELECT 1056 END DO 1057 END DO 1058 END DO 1059 ! 296 1060 ii0 = 44 ; ii1 = 44 ! Lombok Strait (e1v was modified) 297 ij0 = 124 ; ij1 = 125 298 DO jk = 1, jpkm1 ! set the before scale factors at v-points 299 DO jj = mj0(ij0), mj1(ij1) 300 DO ji = mi0(ii0), mi1(ii1) 301 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 302 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 303 END DO 304 END DO 305 END DO 306 1061 ij0 = 124 ; ij1 = 125 1062 DO jk = 1, jpkm1 1063 DO jj = mj0(ij0), mj1(ij1) 1064 DO ji = mi0(ii0), mi1(ii1) 1065 SELECT CASE ( pout ) 1066 CASE( 'V' ) 1067 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1068 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1069 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1070 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1071 END SELECT 1072 END DO 1073 END DO 1074 END DO 1075 ! 307 1076 ii0 = 48 ; ii1 = 48 ! Sumba Strait (e1v was modified) [closed from bathy_11 on] 308 ij0 = 124 ; ij1 = 125 309 DO jk = 1, jpkm1 ! set the before scale factors at v-points 310 DO jj = mj0(ij0), mj1(ij1) 311 DO ji = mi0(ii0), mi1(ii1) 312 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 313 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 314 END DO 315 END DO 316 END DO 317 1077 ij0 = 124 ; ij1 = 125 1078 DO jk = 1, jpkm1 1079 DO jj = mj0(ij0), mj1(ij1) 1080 DO ji = mi0(ii0), mi1(ii1) 1081 SELECT CASE ( pout ) 1082 CASE( 'V' ) 1083 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1084 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1085 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1086 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1087 END SELECT 1088 END DO 1089 END DO 1090 END DO 1091 ! 318 1092 ii0 = 53 ; ii1 = 53 ! Ombai Strait (e1v was modified) 319 ij0 = 124 ; ij1 = 125 320 DO jk = 1, jpkm1 ! set the before scale factors at v-points 321 DO jj = mj0(ij0), mj1(ij1) 322 DO ji = mi0(ii0), mi1(ii1) 323 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 324 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 325 END DO 326 END DO 327 END DO 328 1093 ij0 = 124 ; ij1 = 125 1094 DO jk = 1, jpkm1 1095 DO jj = mj0(ij0), mj1(ij1) 1096 DO ji = mi0(ii0), mi1(ii1) 1097 SELECT CASE ( pout ) 1098 CASE( 'V' ) 1099 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1100 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1101 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1102 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1103 END SELECT 1104 END DO 1105 END DO 1106 END DO 1107 ! 329 1108 ii0 = 56 ; ii1 = 56 ! Timor Passage (e1v was modified) 330 ij0 = 124 ; ij1 = 125 331 DO jk = 1, jpkm1 ! set the before scale factors at v-points 332 DO jj = mj0(ij0), mj1(ij1) 333 DO ji = mi0(ii0), mi1(ii1) 334 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 335 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 336 END DO 337 END DO 338 END DO 339 1109 ij0 = 124 ; ij1 = 125 1110 DO jk = 1, jpkm1 1111 DO jj = mj0(ij0), mj1(ij1) 1112 DO ji = mi0(ii0), mi1(ii1) 1113 SELECT CASE ( pout ) 1114 CASE( 'V' ) 1115 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1116 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1117 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1118 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1119 END SELECT 1120 END DO 1121 END DO 1122 END DO 1123 ! 340 1124 ii0 = 55 ; ii1 = 55 ! West Halmahera Strait (e1v was modified) 341 ij0 = 141 ; ij1 = 142 342 DO jk = 1, jpkm1 ! set the before scale factors at v-points 343 DO jj = mj0(ij0), mj1(ij1) 344 DO ji = mi0(ii0), mi1(ii1) 345 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 346 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 347 END DO 348 END DO 349 END DO 350 1125 ij0 = 141 ; ij1 = 142 1126 DO jk = 1, jpkm1 1127 DO jj = mj0(ij0), mj1(ij1) 1128 DO ji = mi0(ii0), mi1(ii1) 1129 SELECT CASE ( pout ) 1130 CASE( 'V' ) 1131 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1132 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1133 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1134 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1135 END SELECT 1136 END DO 1137 END DO 1138 END DO 1139 ! 351 1140 ii0 = 58 ; ii1 = 58 ! East Halmahera Strait (e1v was modified) 352 ij0 = 141 ; ij1 = 142 353 DO jk = 1, jpkm1 ! set the before scale factors at v-points 354 DO jj = mj0(ij0), mj1(ij1) 355 DO ji = mi0(ii0), mi1(ii1) 356 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 357 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 358 END DO 359 END DO 360 END DO 361 362 ! 1141 ij0 = 141 ; ij1 = 142 1142 DO jk = 1, jpkm1 1143 DO jj = mj0(ij0), mj1(ij1) 1144 DO ji = mi0(ii0), mi1(ii1) 1145 SELECT CASE ( pout ) 1146 CASE( 'V' ) 1147 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1148 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1149 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1150 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1151 END SELECT 1152 END DO 1153 END DO 1154 END DO 363 1155 ENDIF 364 ! ! ======================1156 ! ! ===================== 365 1157 IF( cp_cfg == "orca" .AND. jp_cfg == 05 ) THEN ! ORCA R05 configuration 366 ! ! ====================== 1158 ! ! ===================== 1159 ! 367 1160 ii0 = 563 ; ii1 = 564 ! Gibraltar Strait (e2u was modified) 368 ij0 = 327 ; ij1 = 327 369 DO jk = 1, jpkm1 ! set the before scale factors at u-points 370 DO jj = mj0(ij0), mj1(ij1) 371 DO ji = mi0(ii0), mi1(ii1) 372 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 373 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 374 END DO 375 END DO 376 END DO 377 ! 378 ii0 = 627 ; ii1 = 628 ! Bosphore Strait (e2u was modified) 379 ij0 = 343 ; ij1 = 343 380 DO jk = 1, jpkm1 ! set the before scale factors at u-points 381 DO jj = mj0(ij0), mj1(ij1) 382 DO ji = mi0(ii0), mi1(ii1) 383 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 384 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 1161 ij0 = 327 ; ij1 = 327 1162 DO jk = 1, jpkm1 1163 DO jj = mj0(ij0), mj1(ij1) 1164 DO ji = mi0(ii0), mi1(ii1) 1165 SELECT CASE ( pout ) 1166 CASE( 'U' ) 1167 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1168 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1169 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1170 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1171 CASE( 'F' ) 1172 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1173 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1174 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1175 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1176 END SELECT 1177 END DO 1178 END DO 1179 END DO 1180 ! 1181 ii0 = 627 ; ii1 = 628 ! Bosphorus Strait (e2u was modified) 1182 ij0 = 343 ; ij1 = 343 1183 DO jk = 1, jpkm1 1184 DO jj = mj0(ij0), mj1(ij1) 1185 DO ji = mi0(ii0), mi1(ii1) 1186 SELECT CASE ( pout ) 1187 CASE( 'U' ) 1188 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1189 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1190 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1191 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1192 CASE( 'F' ) 1193 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1194 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1195 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1196 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1197 END SELECT 385 1198 END DO 386 1199 END DO … … 388 1201 ! 389 1202 ii0 = 93 ; ii1 = 94 ! Sumba Strait (e2u was modified) 390 ij0 = 232 ; ij1 = 232 391 DO jk = 1, jpkm1 ! set the before scale factors at u-points 392 DO jj = mj0(ij0), mj1(ij1) 393 DO ji = mi0(ii0), mi1(ii1) 394 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 395 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 1203 ij0 = 232 ; ij1 = 232 1204 DO jk = 1, jpkm1 1205 DO jj = mj0(ij0), mj1(ij1) 1206 DO ji = mi0(ii0), mi1(ii1) 1207 SELECT CASE ( pout ) 1208 CASE( 'U' ) 1209 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1210 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1211 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1212 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1213 CASE( 'F' ) 1214 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1215 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1216 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1217 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1218 END SELECT 396 1219 END DO 397 1220 END DO … … 399 1222 ! 400 1223 ii0 = 103 ; ii1 = 103 ! Ombai Strait (e2u was modified) 401 ij0 = 232 ; ij1 = 232 402 DO jk = 1, jpkm1 ! set the before scale factors at u-points 403 DO jj = mj0(ij0), mj1(ij1) 404 DO ji = mi0(ii0), mi1(ii1) 405 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 406 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 1224 ij0 = 232 ; ij1 = 232 1225 DO jk = 1, jpkm1 1226 DO jj = mj0(ij0), mj1(ij1) 1227 DO ji = mi0(ii0), mi1(ii1) 1228 SELECT CASE ( pout ) 1229 CASE( 'U' ) 1230 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1231 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1232 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1233 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1234 CASE( 'F' ) 1235 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1236 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1237 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1238 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1239 END SELECT 407 1240 END DO 408 1241 END DO … … 410 1243 ! 411 1244 ii0 = 15 ; ii1 = 15 ! Palk Strait (e2u was modified) 412 ij0 = 270 ; ij1 = 270 413 DO jk = 1, jpkm1 ! set the before scale factors at u-points 414 DO jj = mj0(ij0), mj1(ij1) 415 DO ji = mi0(ii0), mi1(ii1) 416 zvt = fse3t_b(ji,jj,jk) * e1t(ji,jj) 417 pe3u_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji+1,jj,jk) * e1t(ji+1,jj) ) / ( e1u(ji,jj) ) 1245 ij0 = 270 ; ij1 = 270 1246 DO jk = 1, jpkm1 1247 DO jj = mj0(ij0), mj1(ij1) 1248 DO ji = mi0(ii0), mi1(ii1) 1249 SELECT CASE ( pout ) 1250 CASE( 'U' ) 1251 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) & 1252 & * ( e1t(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3t_0(ji ,jj,jk) ) & 1253 & + e1t(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3t_0(ji+1,jj,jk) ) & 1254 & ) / e1u(ji,jj) + e3u_0(ji,jj,jk) 1255 CASE( 'F' ) 1256 pe3_out(ji,jj,jk) = 0.5_wp * umask(ji,jj,jk) * umask(ji,jj+1,jk) & 1257 & * ( e1u(ji ,jj) * ( pe3_in(ji ,jj,jk) - e3u_0(ji ,jj,jk) ) & 1258 & + e1u(ji+1,jj) * ( pe3_in(ji+1,jj,jk) - e3u_0(ji+1,jj,jk) ) & 1259 & ) / e1f(ji,jj) + e3f_0(ji,jj,jk) 1260 END SELECT 418 1261 END DO 419 1262 END DO … … 421 1264 ! 422 1265 ii0 = 87 ; ii1 = 87 ! Lombok Strait (e1v was modified) 423 ij0 = 232 ; ij1 = 233 424 DO jk = 1, jpkm1 ! set the before scale factors at v-points 425 DO jj = mj0(ij0), mj1(ij1) 426 DO ji = mi0(ii0), mi1(ii1) 427 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 428 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 1266 ij0 = 232 ; ij1 = 233 1267 DO jk = 1, jpkm1 1268 DO jj = mj0(ij0), mj1(ij1) 1269 DO ji = mi0(ii0), mi1(ii1) 1270 SELECT CASE ( pout ) 1271 CASE( 'V' ) 1272 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1273 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1274 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1275 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1276 END SELECT 429 1277 END DO 430 1278 END DO … … 432 1280 ! 433 1281 ii0 = 662 ; ii1 = 662 ! Bab el Mandeb (e1v was modified) 434 ij0 = 276 ; ij1 = 276 435 DO jk = 1, jpkm1 ! set the before scale factors at v-points 436 DO jj = mj0(ij0), mj1(ij1) 437 DO ji = mi0(ii0), mi1(ii1) 438 zvt = fse3t_b(ji,jj,jk) * e2t(ji,jj) 439 pe3v_b(ji,jj,jk) = 0.5_wp * ( zvt + fse3t_b(ji,jj+1,jk) * e2t(ji,jj+1) ) / ( e2v(ji,jj) ) 440 END DO 441 END DO 442 END DO 443 ! 1282 ij0 = 276 ; ij1 = 276 1283 DO jk = 1, jpkm1 1284 DO jj = mj0(ij0), mj1(ij1) 1285 DO ji = mi0(ii0), mi1(ii1) 1286 SELECT CASE ( pout ) 1287 CASE( 'V' ) 1288 pe3_out(ji,jj,jk) = 0.5_wp * vmask(ji,jj,jk) & 1289 & * ( e2t(ji,jj ) * ( pe3_in(ji,jj ,jk) - e3t_0(ji,jj ,jk) ) & 1290 & + e2t(ji,jj+1) * ( pe3_in(ji,jj+1,jk) - e3t_0(ji,jj+1,jk) ) & 1291 & ) / e2v(ji,jj) + e3v_0(ji,jj,jk) 1292 END SELECT 1293 END DO 1294 END DO 1295 END DO 444 1296 ENDIF 445 ! End of individual corrections to scale factors 446 447 IF( ln_zps ) THEN ! minimum of the e3t at partial cell level 448 DO jj = 2, jpjm1 449 DO ji = fs_2, fs_jpim1 450 iku = mbku(ji,jj) 451 ikv = mbkv(ji,jj) 452 pe3u_b(ji,jj,iku) = MIN( fse3t_b(ji,jj,iku), fse3t_b(ji+1,jj ,iku) ) 453 pe3v_b(ji,jj,ikv) = MIN( fse3t_b(ji,jj,ikv), fse3t_b(ji ,jj+1,ikv) ) 454 END DO 455 END DO 456 ENDIF 457 458 pe3u_b(:,:,:) = pe3u_b(:,:,:) - fse3u_0(:,:,:) ! anomaly to avoid zero along closed boundary/extra halos 459 pe3v_b(:,:,:) = pe3v_b(:,:,:) - fse3v_0(:,:,:) 460 CALL lbc_lnk( pe3u_b(:,:,:), 'U', 1. ) ! lateral boundary conditions 461 CALL lbc_lnk( pe3v_b(:,:,:), 'V', 1. ) 462 pe3u_b(:,:,:) = pe3u_b(:,:,:) + fse3u_0(:,:,:) ! recover the full scale factor 463 pe3v_b(:,:,:) = pe3v_b(:,:,:) + fse3v_0(:,:,:) 464 ! 465 IF( nn_timing == 1 ) CALL timing_stop('dom_vvl_2') 466 ! 467 END SUBROUTINE dom_vvl_2 468 469 #else 470 !!---------------------------------------------------------------------- 471 !! Default option : Empty routine 472 !!---------------------------------------------------------------------- 473 CONTAINS 474 SUBROUTINE dom_vvl 475 END SUBROUTINE dom_vvl 476 SUBROUTINE dom_vvl_2(kdum, pudum, pvdum ) 477 USE par_kind 478 INTEGER , INTENT(in ) :: kdum 479 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pudum, pvdum 480 END SUBROUTINE dom_vvl_2 481 #endif 1297 END SUBROUTINE dom_vvl_orca_fix 482 1298 483 1299 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.
