Changeset 2148 for branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN
- Timestamp:
- 2010-10-04T15:53:42+02:00 (14 years ago)
- Location:
- branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN
- Files:
-
- 5 edited
-
divcur.F90 (modified) (3 diffs)
-
dynnxt.F90 (modified) (6 diffs)
-
dynzdf_exp.F90 (modified) (8 diffs)
-
dynzdf_imp.F90 (modified) (17 diffs)
-
sshwzv.F90 (modified) (15 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/divcur.F90
r1953 r2148 14 14 USE in_out_manager ! I/O manager 15 15 USE obc_oce ! ocean lateral open boundary condition 16 USE bdy_oce 16 USE bdy_oce ! Unstructured open boundaries variables 17 17 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 18 18 … … 87 87 INTEGER :: ii, ij, jl ! temporary integer 88 88 INTEGER :: ijt, iju ! temporary integer 89 REAL(wp) :: zraur, zdep 89 90 REAL(wp), DIMENSION( jpi ,1:jpj+2) :: zwu ! workspace 90 91 REAL(wp), DIMENSION(-1:jpi+2, jpj ) :: zwv ! workspace … … 301 302 !! * Local declarations 302 303 INTEGER :: ji, jj, jk ! dummy loop indices 304 REAL(wp) :: zraur, zdep 303 305 !!---------------------------------------------------------------------- 304 306 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynnxt.F90
r1970 r2148 22 22 USE oce ! ocean dynamics and tracers 23 23 USE dom_oce ! ocean space and time domain 24 USE sbc_oce ! Surface boundary condition: ocean fields 25 USE phycst ! physical constants 24 26 USE dynspg_oce ! type of surface pressure gradient 25 27 USE dynadv ! dynamics: vector invariant versus flux form … … 87 89 !! un,vn now horizontal velocity of next time-step 88 90 !!---------------------------------------------------------------------- 91 USE oce, ONLY : ze3u_f => ta ! use ta as 3D workspace 92 USE oce, ONLY : ze3v_f => sa ! use sa as 3D workspace 89 93 INTEGER, INTENT( in ) :: kt ! ocean time-step index 90 94 !! … … 95 99 REAL(wp) :: zue3a , zue3n , zue3b ! temporary scalar 96 100 REAL(wp) :: zve3a , zve3n , zve3b ! - - 97 REAL(wp) :: ze3u_b, ze3u_n, ze3u_a ! - -98 REAL(wp) :: ze3v_b, ze3v_n, ze3v_a ! - -99 101 REAL(wp) :: zuf , zvf ! - - 102 REAL(wp) :: zec ! - - 103 REAL(wp) :: zv_t_ij , zv_t_ip1j ! - - 104 REAL(wp) :: zv_t_ijp1 ! - - 105 REAL(wp), DIMENSION(jpi,jpj) :: zs_t, zs_u_1, zs_v_1 ! temporary 2D workspace 100 106 !!---------------------------------------------------------------------- 101 107 … … 146 152 # if defined key_obc 147 153 ! !* OBC open boundaries 148 IF( lk_obc )CALL obc_dyn( kt )154 CALL obc_dyn( kt ) 149 155 ! 150 156 IF ( lk_dynspg_exp .OR. lk_dynspg_ts ) THEN … … 212 218 END DO 213 219 ELSE !* Leap-Frog : Asselin filter and swap 214 IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN ! applied on velocity 220 ! ! =============! 221 IF( .NOT. lk_vvl ) THEN ! Fixed volume ! 222 ! ! =============! 215 223 DO jk = 1, jpkm1 216 224 DO jj = 1, jpj 217 225 DO ji = 1, jpi 218 zuf = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk)219 zvf = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk)226 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 227 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 220 228 ! 221 229 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity … … 226 234 END DO 227 235 END DO 228 ELSE ! applied on thickness weighted velocity 236 ! ! ================! 237 ELSE ! Variable volume ! 238 ! ! ================! 239 ! Before scale factor at t-points 240 ! ------------------------------- 229 241 DO jk = 1, jpkm1 230 DO jj = 1, jpj 231 DO ji = 1, jpi 232 ze3u_a = fse3u_a(ji,jj,jk) 233 ze3v_a = fse3v_a(ji,jj,jk) 234 ze3u_n = fse3u_n(ji,jj,jk) 235 ze3v_n = fse3v_n(ji,jj,jk) 236 ze3u_b = fse3u_b(ji,jj,jk) 237 ze3v_b = fse3v_b(ji,jj,jk) 238 ! 239 zue3a = ua(ji,jj,jk) * ze3u_a 240 zve3a = va(ji,jj,jk) * ze3v_a 241 zue3n = un(ji,jj,jk) * ze3u_n 242 zve3n = vn(ji,jj,jk) * ze3v_n 243 zue3b = ub(ji,jj,jk) * ze3u_b 244 zve3b = vb(ji,jj,jk) * ze3v_b 245 ! 246 zuf = ( atfp * ( zue3b + zue3a ) + atfp1 * zue3n ) & 247 & / ( atfp * ( ze3u_b + ze3u_a ) + atfp1 * ze3u_n ) * umask(ji,jj,jk) 248 zvf = ( atfp * ( zve3b + zve3a ) + atfp1 * zve3n ) & 249 & / ( atfp * ( ze3v_b + ze3v_a ) + atfp1 * ze3v_n ) * vmask(ji,jj,jk) 250 ! 251 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 252 vb(ji,jj,jk) = zvf 253 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 254 vn(ji,jj,jk) = va(ji,jj,jk) 255 END DO 256 END DO 257 END DO 242 fse3t_b(:,:,jk) = fse3t_n(:,:,jk) & 243 & + atfp * ( fse3t_b(:,:,jk) + fse3t_a(:,:,jk) & 244 & - 2.e0 * fse3t_n(:,:,jk) ) 245 ENDDO 246 ! Add volume filter correction only at the first level of t-point scale factors 247 zec = atfp * rdt / rau0 248 fse3t_b(:,:,1) = fse3t_b(:,:,1) - zec * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 249 ! surface at t-points and inverse surface at (u/v)-points used in surface averaging computations 250 zs_t (:,:) = e1t(:,:) * e2t(:,:) 251 zs_u_1(:,:) = 0.5 / e1u(:,:) * e2u(:,:) 252 zs_v_1(:,:) = 0.5 / e1v(:,:) * e2v(:,:) 253 ! 254 IF( ln_dynadv_vec ) THEN 255 ! Before scale factor at (u/v)-points 256 ! ----------------------------------- 257 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 258 DO jk = 1, jpkm1 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 261 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 262 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 263 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 264 fse3u_b(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 265 fse3v_b(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 266 END DO 267 END DO 268 END DO 269 ! lateral boundary conditions 270 CALL lbc_lnk( fse3u_b(:,:,:), 'U', 1. ) 271 CALL lbc_lnk( fse3v_b(:,:,:), 'V', 1. ) 272 ! Add initial scale factor to scale factor anomaly 273 fse3u_b(:,:,:) = fse3u_b(:,:,:) + fse3u_0(:,:,:) 274 fse3v_b(:,:,:) = fse3v_b(:,:,:) + fse3v_0(:,:,:) 275 ! Leap-Frog - Asselin filter and swap: applied on velocity 276 ! ----------------------------------- 277 DO jk = 1, jpkm1 278 DO jj = 1, jpj 279 DO ji = 1, jpi 280 zuf = un(ji,jj,jk) + atfp * ( ub(ji,jj,jk) - 2.e0 * un(ji,jj,jk) + ua(ji,jj,jk) ) 281 zvf = vn(ji,jj,jk) + atfp * ( vb(ji,jj,jk) - 2.e0 * vn(ji,jj,jk) + va(ji,jj,jk) ) 282 ! 283 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 284 vb(ji,jj,jk) = zvf 285 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 286 vn(ji,jj,jk) = va(ji,jj,jk) 287 END DO 288 END DO 289 END DO 290 ! 291 ELSE 292 ! Temporary filered scale factor at (u/v)-points (will become before scale factor) 293 !----------------------------------------------- 294 ! Scale factor anomaly at (u/v)-points: surface averaging of scale factor at t-points 295 DO jk = 1, jpkm1 296 DO jj = 1, jpjm1 297 DO ji = 1, jpim1 298 zv_t_ij = zs_t(ji ,jj ) * fse3t_b(ji ,jj ,jk) 299 zv_t_ip1j = zs_t(ji+1,jj ) * fse3t_b(ji+1,jj ,jk) 300 zv_t_ijp1 = zs_t(ji ,jj+1) * fse3t_b(ji ,jj+1,jk) 301 ze3u_f(ji,jj,jk) = umask(ji,jj,jk) * ( zs_u_1(ji,jj) * ( zv_t_ij + zv_t_ip1j ) - fse3u_0(ji,jj,jk) ) 302 ze3v_f(ji,jj,jk) = vmask(ji,jj,jk) * ( zs_v_1(ji,jj) * ( zv_t_ij + zv_t_ijp1 ) - fse3v_0(ji,jj,jk) ) 303 END DO 304 END DO 305 END DO 306 ! lateral boundary conditions 307 CALL lbc_lnk( ze3u_f, 'U', 1. ) 308 CALL lbc_lnk( ze3v_f, 'V', 1. ) 309 ! Add initial scale factor to scale factor anomaly 310 ze3u_f(:,:,:) = ze3u_f(:,:,:) + fse3u_0(:,:,:) 311 ze3v_f(:,:,:) = ze3v_f(:,:,:) + fse3v_0(:,:,:) 312 ! Leap-Frog - Asselin filter and swap: applied on thickness weighted velocity 313 ! ----------------------------------- =========================== 314 DO jk = 1, jpkm1 315 DO jj = 1, jpj 316 DO ji = 1, jpim1 317 zue3a = ua(ji,jj,jk) * fse3u_a(ji,jj,jk) 318 zve3a = va(ji,jj,jk) * fse3v_a(ji,jj,jk) 319 zue3n = un(ji,jj,jk) * fse3u_n(ji,jj,jk) 320 zve3n = vn(ji,jj,jk) * fse3v_n(ji,jj,jk) 321 zue3b = ub(ji,jj,jk) * fse3u_b(ji,jj,jk) 322 zve3b = vb(ji,jj,jk) * fse3v_b(ji,jj,jk) 323 ! 324 zuf = ( zue3n + atfp * ( zue3b - 2.e0 * zue3n + zue3a ) ) / ze3u_f(ji,jj,jk) 325 zvf = ( zve3n + atfp * ( zve3b - 2.e0 * zve3n + zve3a ) ) / ze3v_f(ji,jj,jk) 326 ! 327 ub(ji,jj,jk) = zuf ! ub <-- filtered velocity 328 vb(ji,jj,jk) = zvf 329 un(ji,jj,jk) = ua(ji,jj,jk) ! un <-- ua 330 vn(ji,jj,jk) = va(ji,jj,jk) 331 END DO 332 END DO 333 END DO 334 fse3u_b(:,:,:) = ze3u_f(:,:,:) ! e3u_b <-- filtered scale factor 335 fse3v_b(:,:,:) = ze3v_f(:,:,:) 336 CALL lbc_lnk( ub, 'U', -1. ) ! lateral boundary conditions 337 CALL lbc_lnk( vb, 'V', -1. ) 338 ENDIF 339 ! 258 340 ENDIF 341 ! 259 342 ENDIF 260 343 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_exp.F90
r1537 r2148 4 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 5 5 !!============================================================================== 6 !! History : ! 90-10 (B. Blanke) Original code 7 !! ! 97-05 (G. Madec) vertical component of isopycnal 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal 8 !! NEMO 1.0 ! 1002-08 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! sion using an explicit time-stepping scheme. 14 15 !!---------------------------------------------------------------------- 15 !! * Modules used16 16 USE oce ! ocean dynamics and tracers 17 17 USE dom_oce ! ocean space and time domain … … 24 24 PRIVATE 25 25 26 !! * Routine accessibility 27 PUBLIC dyn_zdf_exp ! called by step.F90 26 PUBLIC dyn_zdf_exp ! called by step.F90 28 27 29 28 !! * Substitutions … … 31 30 # include "vectopt_loop_substitute.h90" 32 31 !!---------------------------------------------------------------------- 33 !! OPA 9.0 , LOCEAN-IPSL (2005)32 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 34 33 !! $Id$ 35 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 47 46 !! technique. The vertical diffusion of momentum is given by: 48 47 !! diffu = dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ub) ) 49 !! Surface boundary conditions: wind stress input 48 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 50 49 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F90) 51 50 !! Add this trend to the general trend ua : … … 54 53 !! ** Action : - Update (ua,va) with the vertical diffusive trend 55 54 !!--------------------------------------------------------------------- 56 !! * Arguments 57 INTEGER , INTENT( in ) :: kt ! ocean time-step index 58 REAL(wp), INTENT( in ) :: p2dt ! time-step 59 60 !! * Local declarations 61 INTEGER :: ji, jj, jk, jl ! dummy loop indices 62 REAL(wp) :: zrau0r, zlavmr, zua, zva ! temporary scalars 63 REAL(wp), DIMENSION(jpi,jpk) :: zwx, zwy, zwz, zww ! temporary workspace arrays 55 INTEGER , INTENT(in) :: kt ! ocean time-step index 56 REAL(wp), INTENT(in) :: p2dt ! time-step 57 !! 58 INTEGER :: ji, jj, jk, jl ! dummy loop indices 59 REAL(wp) :: zrau0r, zlavmr, zua, zva ! temporary scalars 60 REAL(wp), DIMENSION(jpi,jpk) :: zwx, zwy, zwz, zww ! 2D workspace 64 61 !!---------------------------------------------------------------------- 65 62 66 63 IF( kt == nit000 ) THEN 67 64 IF(lwp) WRITE(numout,*) 68 IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion explicit operator'65 IF(lwp) WRITE(numout,*) 'dyn_zdf_exp : vertical momentum diffusion - explicit operator' 69 66 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 70 67 ENDIF 71 68 72 ! Local constant initialization 73 ! ----------------------------- 74 zrau0r = 1. / rau0 ! inverse of the reference density 75 zlavmr = 1. / float( nn_zdfexp ) ! inverse of the number of sub time step 69 zrau0r = 1. / rau0 ! Local constant initialization 70 zlavmr = 1. / REAL( nn_zdfexp ) 76 71 77 ! ! =============== 78 DO jj = 2, jpjm1 ! Vertical slab 79 ! ! =============== 80 81 ! Surface boundary condition 82 DO ji = 2, jpim1 83 zwy(ji,1) = utau(ji,jj) * zrau0r 84 zww(ji,1) = vtau(ji,jj) * zrau0r 72 ! ! =============== 73 DO jj = 2, jpjm1 ! Vertical slab 74 ! ! =============== 75 DO ji = 2, jpim1 ! Surface boundary condition 76 zwy(ji,1) = ( utau_b(ji,jj) + utau(ji,jj) ) * zrau0r 77 zww(ji,1) = ( vtau_b(ji,jj) + vtau(ji,jj) ) * zrau0r 85 78 END DO 86 87 ! Initialization of x, z and contingently trends array 88 DO jk = 1, jpk 79 DO jk = 1, jpk ! Initialization of x, z and contingently trends array 89 80 DO ji = 2, jpim1 90 81 zwx(ji,jk) = ub(ji,jj,jk) … … 92 83 END DO 93 84 END DO 94 95 ! Time splitting loop 96 DO jl = 1, nn_zdfexp 97 98 ! First vertical derivative 99 DO jk = 2, jpk 85 ! 86 DO jl = 1, nn_zdfexp ! Time splitting loop 87 ! 88 DO jk = 2, jpk ! First vertical derivative 100 89 DO ji = 2, jpim1 101 90 zwy(ji,jk) = avmu(ji,jj,jk) * ( zwx(ji,jk-1) - zwx(ji,jk) ) / fse3uw(ji,jj,jk) … … 103 92 END DO 104 93 END DO 105 106 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 107 DO jk = 1, jpkm1 94 DO jk = 1, jpkm1 ! Second vertical derivative and trend estimation at kt+l*rdt/nn_zdfexp 108 95 DO ji = 2, jpim1 109 zua = zlavmr *( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk)110 zva = zlavmr *( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk)96 zua = zlavmr * ( zwy(ji,jk) - zwy(ji,jk+1) ) / fse3u(ji,jj,jk) 97 zva = zlavmr * ( zww(ji,jk) - zww(ji,jk+1) ) / fse3v(ji,jj,jk) 111 98 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 112 99 va(ji,jj,jk) = va(ji,jj,jk) + zva 113 114 zwx(ji,jk) = zwx(ji,jk) + p2dt *zua*umask(ji,jj,jk)115 zwz(ji,jk) = zwz(ji,jk) + p2dt *zva*vmask(ji,jj,jk)100 ! 101 zwx(ji,jk) = zwx(ji,jk) + p2dt * zua * umask(ji,jj,jk) 102 zwz(ji,jk) = zwz(ji,jk) + p2dt * zva * vmask(ji,jj,jk) 116 103 END DO 117 104 END DO 118 105 ! 119 106 END DO 120 121 ! ! =============== 122 END DO ! End of slab 123 ! ! =============== 124 107 ! ! =============== 108 END DO ! End of slab 109 ! ! =============== 125 110 END SUBROUTINE dyn_zdf_exp 126 111 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r1662 r2148 4 4 !! Ocean dynamics: vertical component(s) of the momentum mixing trend 5 5 !!============================================================================== 6 !! History : ! 90-10 (B. Blanke) Original code 7 !! ! 97-05 (G. Madec) vertical component of isopycnal 8 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 6 !! History : OPA ! 1990-10 (B. Blanke) Original code 7 !! 8.0 ! 1997-05 (G. Madec) vertical component of isopycnal 8 !! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 9 10 !!---------------------------------------------------------------------- 10 11 … … 13 14 !! sion using a implicit time-stepping. 14 15 !!---------------------------------------------------------------------- 15 !! OPA 9.0 , LOCEAN-IPSL (2005)16 !! $Id$17 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)18 !!----------------------------------------------------------------------19 !! * Modules used20 16 USE oce ! ocean dynamics and tracers 21 17 USE dom_oce ! ocean space and time domain … … 28 24 PRIVATE 29 25 30 !! * Routine accessibility 31 PUBLIC dyn_zdf_imp ! called by step.F90 26 PUBLIC dyn_zdf_imp ! called by step.F90 32 27 33 28 !! * Substitutions … … 35 30 # include "vectopt_loop_substitute.h90" 36 31 !!---------------------------------------------------------------------- 37 !! OPA 9.0 , LOCEAN-IPSL (2005)32 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 38 33 !! $Id$ 39 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 40 35 !!---------------------------------------------------------------------- 41 36 42 37 CONTAINS 43 44 38 45 39 SUBROUTINE dyn_zdf_imp( kt, p2dt ) … … 54 48 !! dz( avmu dz(u) ) = 1/e3u dk+1( avmu/e3uw dk(ua) ) 55 49 !! backward time stepping 56 !! Surface boundary conditions: wind stress input 50 !! Surface boundary conditions: wind stress input (averaged over kt-1/2 & kt+1/2) 57 51 !! Bottom boundary conditions : bottom stress (cf zdfbfr.F) 58 52 !! Add this trend to the general trend ua : 59 53 !! ua = ua + dz( avmu dz(u) ) 60 54 !! 61 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive 62 !! mixing trend. 55 !! ** Action : - Update (ua,va) arrays with the after vertical diffusive mixing trend. 63 56 !!--------------------------------------------------------------------- 64 !! * Modules used 65 USE oce, ONLY : zwd => ta, & ! use ta as workspace 66 zws => sa ! use sa as workspace 67 68 !! * Arguments 69 INTEGER , INTENT( in ) :: kt ! ocean time-step index 70 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 71 72 !! * Local declarations 73 INTEGER :: ji, jj, jk ! dummy loop indices 74 REAL(wp) :: zrau0r, z2dtf, zcoef, zzws, zrhs ! temporary scalars 75 REAL(wp) :: zzwi ! temporary scalars 76 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! temporary workspace arrays 57 USE oce, ONLY : zwd => ta ! use ta as workspace 58 USE oce, ONLY : zws => sa ! use sa as workspace 59 !! 60 INTEGER , INTENT( in ) :: kt ! ocean time-step index 61 REAL(wp), INTENT( in ) :: p2dt ! vertical profile of tracer time-step 62 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: zrau0r, zcoef ! temporary scalars 65 REAL(wp) :: zzwi, zzws, zrhs ! temporary scalars 66 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi ! 3D workspace 77 67 !!---------------------------------------------------------------------- 78 68 … … 95 85 ! friction velocity in dyn_bfr using values from the previous timestep. There 96 86 ! is no need to include these in the implicit calculation. 97 DO jk = 1, jpkm1 87 ! 88 DO jk = 1, jpkm1 ! Matrix 98 89 DO jj = 2, jpjm1 99 90 DO ji = fs_2, fs_jpim1 ! vector opt. 100 91 zcoef = - p2dt / fse3u(ji,jj,jk) 101 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk )102 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk)103 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1)104 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1)92 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 93 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 94 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 95 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 105 96 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws 106 97 END DO 107 98 END DO 108 99 END DO 109 110 ! Surface boudary conditions 111 DO jj = 2, jpjm1 100 DO jj = 2, jpjm1 ! Surface boudary conditions 112 101 DO ji = fs_2, fs_jpim1 ! vector opt. 113 102 zwi(ji,jj,1) = 0. … … 130 119 ! The solution (the after velocity) is in ua 131 120 !----------------------------------------------------------------------- 132 133 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 134 DO jk = 2, jpkm1 121 ! 122 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 135 123 DO jj = 2, jpjm1 136 124 DO ji = fs_2, fs_jpim1 ! vector opt. … … 139 127 END DO 140 128 END DO 141 142 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 143 DO jj = 2, jpjm1 144 DO ji = fs_2, fs_jpim1 ! vector opt. 145 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 146 !!! ua(ji,jj,1) = ub(ji,jj,1) & 147 !!! + p2dt * ( ua(ji,jj,1) + utau(ji,jj) / ( fse3u(ji,jj,1)*rau0 ) ) 148 z2dtf = p2dt / ( fse3u(ji,jj,1)*rau0 ) 149 ua(ji,jj,1) = ub(ji,jj,1) & 150 + p2dt * ua(ji,jj,1) + z2dtf * utau(ji,jj) 129 ! 130 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 131 DO ji = fs_2, fs_jpim1 ! vector opt. 132 ua(ji,jj,1) = ub(ji,jj,1) & 133 & + p2dt * ( ua(ji,jj,1) + 0.5 * ( utau_b(ji,jj) + utau(ji,jj) ) / ( fse3u(ji,jj,1) * rau0 ) ) 151 134 END DO 152 135 END DO … … 159 142 END DO 160 143 END DO 161 162 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 163 DO jj = 2, jpjm1 144 ! 145 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 164 146 DO ji = fs_2, fs_jpim1 ! vector opt. 165 147 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 173 155 END DO 174 156 END DO 175 176 ! Normalization to obtain the general momentum trend ua 177 DO jk = 1, jpkm1 157 ! 158 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend ua == 178 159 DO jj = 2, jpjm1 179 160 DO ji = fs_2, fs_jpim1 ! vector opt. … … 192 173 ! friction velocity in dyn_bfr using values from the previous timestep. There 193 174 ! is no need to include these in the implicit calculation. 194 DO jk = 1, jpkm1 175 ! 176 DO jk = 1, jpkm1 ! Matrix 195 177 DO jj = 2, jpjm1 196 178 DO ji = fs_2, fs_jpim1 ! vector opt. 197 179 zcoef = -p2dt / fse3v(ji,jj,jk) 198 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk )180 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 199 181 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 200 zzws = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1)182 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 201 183 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 202 184 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zzws … … 204 186 END DO 205 187 END DO 206 207 ! Surface boudary conditions 208 DO jj = 2, jpjm1 188 DO jj = 2, jpjm1 ! Surface boudary conditions 209 189 DO ji = fs_2, fs_jpim1 ! vector opt. 210 190 zwi(ji,jj,1) = 0.e0 … … 223 203 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 224 204 ! 225 ! m is decomposed in the product of an upper and lower triangular 226 ! matrix 205 ! m is decomposed in the product of an upper and lower triangular matrix 227 206 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 228 207 ! The solution (after velocity) is in 2d array va 229 208 !----------------------------------------------------------------------- 230 231 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 232 DO jk = 2, jpkm1 209 ! 210 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 233 211 DO jj = 2, jpjm1 234 212 DO ji = fs_2, fs_jpim1 ! vector opt. … … 237 215 END DO 238 216 END DO 239 240 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 241 DO jj = 2, jpjm1 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 !!! change les resultats (derniers digit, pas significativement + rapide 1* de moins) 244 !!! va(ji,jj,1) = vb(ji,jj,1) & 245 !!! + p2dt * ( va(ji,jj,1) + vtau(ji,jj) / ( fse3v(ji,jj,1)*rau0 ) ) 246 z2dtf = p2dt / ( fse3v(ji,jj,1)*rau0 ) 247 va(ji,jj,1) = vb(ji,jj,1) & 248 + p2dt * va(ji,jj,1) + z2dtf * vtau(ji,jj) 217 ! 218 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 va(ji,jj,1) = vb(ji,jj,1) & 221 & + p2dt * ( va(ji,jj,1) + 0.5 * ( vtau_b(ji,jj) + vtau(ji,jj) ) / ( fse3v(ji,jj,1) * rau0 ) ) 249 222 END DO 250 223 END DO … … 257 230 END DO 258 231 END DO 259 260 ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 261 DO jj = 2, jpjm1 232 ! 233 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 262 234 DO ji = fs_2, fs_jpim1 ! vector opt. 263 235 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) … … 271 243 END DO 272 244 END DO 273 274 ! Normalization to obtain the general momentum trend va 275 DO jk = 1, jpkm1 245 ! 246 DO jk = 1, jpkm1 !== Normalization to obtain the general momentum trend va == 276 247 DO jj = 2, jpjm1 277 248 DO ji = fs_2, fs_jpim1 ! vector opt. … … 280 251 END DO 281 252 END DO 282 253 ! 283 254 END SUBROUTINE dyn_zdf_imp 284 255 -
branches/DEV_r2106_LOCEAN2010/NEMO/OPA_SRC/DYN/sshwzv.F90
r2113 r2148 5 5 !!============================================================================== 6 6 !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code 7 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) modified LF-RA 7 8 !!---------------------------------------------------------------------- 8 9 … … 27 28 USE diaar5, ONLY : lk_diaar5 28 29 USE iom 29 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 USE sbcrnf, ONLY : rnf_dep, rnf_mod_dep ! River runoff 30 31 31 32 IMPLICIT NONE … … 38 39 # include "domzgr_substitute.h90" 39 40 # include "vectopt_loop_substitute.h90" 40 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 41 !!---------------------------------------------------------------------- 42 !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010) 43 43 !! $Id$ 44 44 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 54 54 !! and update the now vertical coordinate (lk_vvl=T). 55 55 !! 56 !! ** Method : - 57 !! 58 !! - Using the incompressibility hypothesis, the vertical 56 !! ** Method : - Using the incompressibility hypothesis, the vertical 59 57 !! velocity is computed by integrating the horizontal divergence 60 58 !! from the bottom to the surface minus the scale factor evolution. … … 63 61 !! ** action : ssha : after sea surface height 64 62 !! wn : now vertical velocity 65 !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height 66 !! at u-, v-, f-point s 67 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 63 !! sshu_a, sshv_a, sshf_a : after sea surface height (lk_vvl=T) 64 !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points 65 !! 66 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 68 67 !!---------------------------------------------------------------------- 69 68 USE oce, ONLY : z3d => ta ! use ta as 3D workspace … … 73 72 INTEGER :: ji, jj, jk ! dummy loop indices 74 73 REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars 75 REAL(wp) :: z2dt, z raur! temporary scalars74 REAL(wp) :: z2dt, z1_2dt, z1_rau0 ! temporary scalars 76 75 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace 77 76 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace … … 79 78 80 79 IF( kt == nit000 ) THEN 80 ! 81 81 IF(lwp) WRITE(numout,*) 82 82 IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' … … 95 95 sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 96 96 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 97 sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) &98 & + sshb(ji+1,jj) + sshb(ji+1,jj+1) )99 97 sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & 100 98 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) 101 99 sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & 102 100 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) 103 sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) &104 & + sshn(ji+1,jj) + sshn(ji+1,jj+1) )105 101 END DO 106 102 END DO 107 103 CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) 108 104 CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) 109 CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) 105 DO jj = 1, jpjm1 106 DO ji = 1, jpim1 ! NO Vector Opt. 107 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 108 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 109 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 110 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 111 END DO 112 END DO 113 CALL lbc_lnk( sshf_n, 'F', 1. ) 110 114 ENDIF 111 115 ! 112 116 ENDIF 113 117 114 ! !------------------------------ !115 IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case)116 ! !------------------------------!118 ! !------------------------------------------! 119 IF( lk_vvl ) THEN ! Regridding: Update Now Vertical coord. ! (only in vvl case) 120 ! !------------------------------------------! 117 121 DO jk = 1, jpkm1 118 fsdept(:,:,jk) = fsdept_n(:,:,jk) 122 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays 119 123 fsdepw(:,:,jk) = fsdepw_n(:,:,jk) 120 124 fsde3w(:,:,jk) = fsde3w_n(:,:,jk) 121 125 ! 122 fse3t (:,:,jk) = fse3t_n (:,:,jk) 126 fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays 123 127 fse3u (:,:,jk) = fse3u_n (:,:,jk) 124 128 fse3v (:,:,jk) = fse3v_n (:,:,jk) … … 128 132 fse3vw(:,:,jk) = fse3vw_n(:,:,jk) 129 133 END DO 130 ! ! now ocean depth (at u- and v-points)131 hu(:,:) = hu_0(:,:) + sshu_n(:,:) 134 ! 135 hu(:,:) = hu_0(:,:) + sshu_n(:,:) ! now ocean depth (at u- and v-points) 132 136 hv(:,:) = hv_0(:,:) + sshv_n(:,:) 133 ! 137 ! ! now masked inverse of the ocean depth (at u- and v-points) 134 138 hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) 135 139 hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) 136 140 ! 137 ! 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 rnf_dep(ji,jj) = 0. 141 DO jk = 1, rnf_mod_dep(ji,jj) ! recalculates rnf_dep to be the depth 142 rnf_dep(ji,jj) = rnf_dep(ji,jj) + fse3t(ji,jj,jk) ! in metres to the bottom of the relevant grid box 143 ENDDO 144 ENDDO 145 ENDDO 146 ! 147 ENDIF 148 149 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 150 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 151 152 ! set time step size (Euler/Leapfrog) 153 z2dt = 2. * rdt 141 ENDIF 142 ! 143 144 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 145 IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) 146 147 z2dt = 2. * rdt ! set time step size (Euler/Leapfrog) 154 148 IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt 155 156 zraur = 1. / rau0157 149 158 150 ! !------------------------------! … … 163 155 zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) 164 156 END DO 165 166 157 ! ! Sea surface elevation time stepping 167 ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 158 ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used 159 ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b -2*rnf ) = emp 160 z1_rau0 = 0.5 / rau0 161 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) - 2 * rnf(:,:) ) + zhdiv(:,:) ) ) & 162 & * tmask(:,:,1) 168 163 169 164 #if defined key_obc 170 IF 165 IF( Agrif_Root() ) THEN 171 166 ssha(:,:) = ssha(:,:) * obctmsk(:,:) 172 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm)167 CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) 173 168 ENDIF 174 169 #endif 175 176 170 ! ! Sea Surface Height at u-,v- and f-points (vvl case only) 177 171 IF( lk_vvl ) THEN ! (required only in key_vvl case) … … 184 178 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 185 179 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 186 sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) &187 & * ( ssha(ji ,jj) + ssha(ji ,jj+1) &188 & + ssha(ji+1,jj) + ssha(ji+1,jj+1) )189 180 END DO 190 181 END DO 191 CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions 182 ! Boundaries conditions 183 CALL lbc_lnk( sshu_a, 'U', 1. ) 192 184 CALL lbc_lnk( sshv_a, 'V', 1. ) 193 CALL lbc_lnk( sshf_a, 'F', 1. ) 194 ENDIF 195 185 ENDIF 196 186 ! !------------------------------! 197 187 ! ! Now Vertical Velocity ! 198 188 ! !------------------------------! 199 ! ! integrate from the bottom the hor. divergence 200 DO jk = jpkm1, 1, -1 201 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 202 & - ( fse3t_a(:,:,jk) & 203 & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt 189 z1_2dt = 1.e0 / z2dt 190 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise 192 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & 193 & - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) ) & 194 & * tmask(:,:,jk) * z1_2dt 204 195 END DO 205 ! 196 197 ! !------------------------------! 198 ! ! outputs ! 199 ! !------------------------------! 206 200 CALL iom_put( "woce", wn ) ! vertical velocity 207 201 CALL iom_put( "ssh" , sshn ) ! sea surface height 208 202 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height 209 IF( lk_diaar5 ) THEN 203 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value 204 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport. 210 205 z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) 211 206 DO jk = 1, jpk 212 207 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) 213 208 END DO 214 CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport 215 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport 216 ENDIF 209 CALL iom_put( "w_masstr" , z3d ) 210 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) 211 ENDIF 212 ! 213 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) 217 214 ! 218 215 END SUBROUTINE ssh_wzv … … 227 224 !! ssha already computed in ssh_wzv 228 225 !! 229 !! ** Method : - apply Asselin time fiter to now ssh and swap : 230 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 231 !! sshn = ssha 226 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 227 !! from the filter, see Leclair and Madec 2010) and swap : 228 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) 229 !! - atfp * rdt * ( emp_b - emp ) / rau0 230 !! sshn = ssha 232 231 !! 233 232 !! ** action : - sshb, sshn : before & now sea surface height 234 233 !! ready for the next time step 235 !!---------------------------------------------------------------------- 236 INTEGER, INTENT( in ) :: kt ! ocean time-step index 237 !! 238 INTEGER :: ji, jj ! dummy loop indices 234 !! 235 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 236 !!---------------------------------------------------------------------- 237 INTEGER, INTENT(in) :: kt ! ocean time-step index 238 !! 239 INTEGER :: ji, jj ! dummy loop indices 240 REAL(wp) :: zec ! temporary scalar 239 241 !!---------------------------------------------------------------------- 240 242 … … 245 247 ENDIF 246 248 247 ! Time filter and swap of the ssh 248 ! ------------------------------- 249 250 IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points 251 ! ! ---------------------- ! 249 ! !--------------------------! 250 IF( lk_vvl ) THEN ! Variable volume levels ! 251 ! !--------------------------! 252 ! 253 ! ssh at t-, u-, v, f-points 254 !=========================== 252 255 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 253 256 sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) 254 257 sshu_n(:,:) = sshu_a(:,:) 255 258 sshv_n(:,:) = sshv_a(:,:) 256 sshf_n(:,:) = sshf_a(:,:) 259 DO jj = 1, jpjm1 260 DO ji = 1, jpim1 ! NO Vector Opt. 261 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 262 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 263 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 264 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 265 END DO 266 END DO 267 ! Boundaries conditions 268 CALL lbc_lnk( sshf_n, 'F', 1. ) 257 269 ELSE ! Leap-Frog time-stepping: Asselin filter + swap 270 zec = atfp * rdt / rau0 258 271 DO jj = 1, jpj 259 272 DO ji = 1, jpi ! before <-- now filtered 260 sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) 261 sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) 262 sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) 263 sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) 273 sshb (ji,jj) = sshn (ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) & 274 & - zec * ( emp_b(ji,jj) - emp(ji,jj) ) * tmask(ji,jj,1) 264 275 sshn (ji,jj) = ssha (ji,jj) ! now <-- after 265 276 sshu_n(ji,jj) = sshu_a(ji,jj) 266 277 sshv_n(ji,jj) = sshv_a(ji,jj) 267 sshf_n(ji,jj) = sshf_a(ji,jj) 268 END DO 269 END DO 278 END DO 279 END DO 280 DO jj = 1, jpjm1 281 DO ji = 1, jpim1 ! NO Vector Opt. 282 sshf_n(ji,jj) = 0.5 * umask(ji,jj,1) * umask(ji,jj+1,1) & 283 & / ( e1f(ji,jj ) * e2f(ji,jj ) ) & 284 & * ( e1u(ji,jj ) * e2u(ji,jj ) * sshu_n(ji,jj ) & 285 & + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) ) 286 END DO 287 END DO 288 ! Boundaries conditions 289 CALL lbc_lnk( sshf_n, 'F', 1. ) 290 DO jj = 1, jpjm1 291 DO ji = 1, jpim1 ! NO Vector Opt. 292 sshu_b(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 293 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 294 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 295 sshv_b(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 296 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 297 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 298 END DO 299 END DO 300 ! Boundaries conditions 301 CALL lbc_lnk( sshu_b, 'U', 1. ) 302 CALL lbc_lnk( sshv_b, 'V', 1. ) 270 303 ENDIF 271 ! 272 ELSE ! fixed levels : ssh at t-point only 273 ! ! ------------ ! 304 ! !--------------------------! 305 ELSE ! fixed levels ! 306 ! !--------------------------! 307 ! 308 ! ssh at t-point only 309 !==================== 274 310 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter 275 311 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) … … 278 314 DO jj = 1, jpj 279 315 DO ji = 1, jpi ! before <-- now filtered 280 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 316 sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) 281 317 sshn(ji,jj) = ssha(ji,jj) ! now <-- after 282 318 END DO … … 286 322 ENDIF 287 323 ! 288 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 )324 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 289 325 ! 290 326 END SUBROUTINE ssh_nxt
Note: See TracChangeset
for help on using the changeset viewer.
