Changeset 592 for trunk/NEMO/OPA_SRC/DYN
- Timestamp:
- 2007-02-09T10:15:25+01:00 (17 years ago)
- Location:
- trunk/NEMO/OPA_SRC/DYN
- Files:
-
- 7 edited
-
dynhpg.F90 (modified) (5 diffs)
-
dynnxt.F90 (modified) (5 diffs)
-
dynspg_exp.F90 (modified) (1 diff)
-
dynspg_flt.F90 (modified) (6 diffs)
-
dynspg_oce.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (12 diffs)
-
wzvmod.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynhpg.F90
r541 r592 160 160 WRITE(numout,*) ' s-coord. (ROTated axes scheme) ln_hpg_rot = ', ln_hpg_rot 161 161 WRITE(numout,*) ' weighting coeff. (wdj scheme) gamm = ', gamm 162 ENDIF 163 164 IF( lk_vvl .AND. .NOT. ln_hpg_sco ) THEN 165 CALL ctl_stop( 'hpg_ctl : variable volume key_vvl compatible only with the standard jacobian formulation hpg_sco') 162 166 ENDIF 163 167 … … 384 388 INTEGER, INTENT(in) :: kt ! ocean time-step index 385 389 !! 386 INTEGER :: ji, jj, jk ! dummy loop indices387 REAL(wp) :: zcoef0, zuap, zvap ! temporary scalars390 INTEGER :: ji, jj, jk ! dummy loop indices 391 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 388 392 !!---------------------------------------------------------------------- 389 393 … … 396 400 ! Local constant initialization 397 401 zcoef0 = - grav * 0.5 402 ! To use density and not density anomaly 403 IF ( lk_vvl ) THEN ; znad = 1. ! Variable volume 404 ELSE ; znad = 0.e0 ! Fixed volume 405 ENDIF 398 406 399 407 ! Surface value … … 401 409 DO ji = fs_2, fs_jpim1 ! vector opt. 402 410 ! hydrostatic pressure gradient along s-surfaces 403 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * rhd(ji+1,jj ,1) &404 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )405 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * rhd(ji ,jj+1,1) &406 & - fse3w(ji ,jj ,1) * rhd(ji ,jj ,1) )411 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 412 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 413 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 414 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 407 415 ! s-coordinate pressure gradient correction 408 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) &416 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2*znad ) & 409 417 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 410 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) ) &418 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2*znad ) & 411 419 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 412 420 ! add to the general momentum trend … … 422 430 ! hydrostatic pressure gradient along s-surfaces 423 431 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 424 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) &425 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) )432 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 433 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 426 434 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 427 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) &428 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) )435 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 436 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 429 437 ! s-coordinate pressure gradient correction 430 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) ) &438 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2*znad ) & 431 439 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 432 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) ) &440 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2*znad ) & 433 441 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 434 442 ! add to the general momentum trend -
trunk/NEMO/OPA_SRC/DYN/dynnxt.F90
r392 r592 21 21 USE agrif_opa_update 22 22 USE agrif_opa_interp 23 USE domvvl ! variable volume 23 24 24 25 IMPLICIT NONE … … 27 28 !! * Accessibility 28 29 PUBLIC dyn_nxt ! routine called by step.F90 30 !! * Substitutions 31 # include "domzgr_substitute.h90" 29 32 !!---------------------------------------------------------------------- 30 33 … … 70 73 INTEGER :: ji, jj, jk ! dummy loop indices 71 74 REAL(wp) :: z2dt ! temporary scalar 75 !! Variable volume 76 REAL(wp) :: zsshun, zsshvn ! temporary scalars 77 REAL(wp), DIMENSION(jpi,jpj) :: & 78 zsshub, zsshua, zsshvb, zsshva ! 2D workspace 79 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 80 zfse3ub, zfse3un, zfse3ua, & ! 3D workspace 81 zfse3vb, zfse3vn, zfse3va 72 82 !!---------------------------------------------------------------------- 73 83 !! OPA 9.0 , LOCEAN-IPSL (2005) … … 85 95 z2dt = 2. * rdt 86 96 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt 97 98 !! Explicit physics with thickness weighted updates 99 IF( lk_vvl .AND. .NOT. lk_dynspg_flt ) THEN 100 101 ! Sea surface elevation time stepping 102 ! ----------------------------------- 103 ! 104 DO jj = 1, jpjm1 105 DO ji = 1,jpim1 106 107 ! Sea Surface Height at u-point before 108 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 109 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 110 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * sshbb(ji+1,jj ) ) 111 112 ! Sea Surface Height at v-point before 113 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 114 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * sshbb(ji ,jj ) & 115 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * sshbb(ji ,jj+1) ) 116 117 ! Sea Surface Height at u-point after 118 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 119 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 120 & + e1t(ji+1,jj ) * e2t(ji+1,jj ) * ssha(ji+1,jj ) ) 121 122 ! Sea Surface Height at v-point after 123 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 124 & * ( e1t(ji ,jj ) * e2t(ji ,jj ) * ssha(ji ,jj ) & 125 & + e1t(ji ,jj+1) * e2t(ji ,jj+1) * ssha(ji ,jj+1) ) 126 127 END DO 128 END DO 129 130 ! Boundaries conditions 131 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. ) 132 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. ) 133 134 ! Scale factors at before and after time step 135 ! ------------------------------------------- 136 DO jk = 1, jpkm1 137 zfse3ub(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshub(:,:) * muu(:,:,jk) ) 138 zfse3ua(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshua(:,:) * muu(:,:,jk) ) 139 zfse3vb(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshvb(:,:) * muv(:,:,jk) ) 140 zfse3va(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshva(:,:) * muv(:,:,jk) ) 141 END DO 142 143 ! Asselin filtered scale factor at now time step 144 ! ---------------------------------------------- 145 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 146 zfse3un(:,:,:) = fse3u(:,:,:) 147 zfse3vn(:,:,:) = fse3v(:,:,:) 148 ELSE 149 DO jk = 1, jpkm1 150 DO jj = 1, jpj 151 DO ji = 1, jpi 152 zsshun = atfp * ( zsshub(ji,jj) + zsshua(ji,jj) ) + atfp1 * sshu(ji,jj) 153 zsshvn = atfp * ( zsshvb(ji,jj) + zsshva(ji,jj) ) + atfp1 * sshv(ji,jj) 154 zfse3un(ji,jj,jk) = fsve3u(ji,jj,jk) * ( 1 + zsshun * muu(ji,jj,jk) ) 155 zfse3vn(ji,jj,jk) = fsve3v(ji,jj,jk) * ( 1 + zsshvn * muv(ji,jj,jk) ) 156 END DO 157 END DO 158 END DO 159 ENDIF 160 161 ! Thickness weighting 162 ! ------------------- 163 ua(:,:,1:jpkm1) = ua(:,:,1:jpkm1) * fse3u (:,:,1:jpkm1) 164 va(:,:,1:jpkm1) = va(:,:,1:jpkm1) * fse3v (:,:,1:jpkm1) 165 166 un(:,:,1:jpkm1) = un(:,:,1:jpkm1) * fse3u (:,:,1:jpkm1) 167 vn(:,:,1:jpkm1) = vn(:,:,1:jpkm1) * fse3v (:,:,1:jpkm1) 168 169 ub(:,:,1:jpkm1) = ub(:,:,1:jpkm1) * zfse3ub(:,:,1:jpkm1) 170 vb(:,:,1:jpkm1) = vb(:,:,1:jpkm1) * zfse3vb(:,:,1:jpkm1) 171 172 ENDIF 87 173 88 174 ! Lateral boundary conditions on ( ua, va ) … … 146 232 # endif 147 233 #endif 234 148 235 ! Time filter and swap of dynamics arrays 149 236 ! ------------------------------------------ 150 237 IF( neuler == 0 .AND. kt == nit000 ) THEN 151 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 152 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 153 ! Euler (forward) time stepping 154 ub(ji,jj,jk) = un(ji,jj,jk) 155 vb(ji,jj,jk) = vn(ji,jj,jk) 156 un(ji,jj,jk) = ua(ji,jj,jk) 157 vn(ji,jj,jk) = va(ji,jj,jk) 158 END DO 159 END DO 238 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 239 ! caution: don't use (:,:) for this loop 240 ! it causes optimization problems on NEC in auto-tasking 241 DO jj = 1, jpj 242 DO ji = 1, jpi 243 zsshun = umask(ji,jj,jk) / fse3u(ji,jj,jk) 244 zsshvn = vmask(ji,jj,jk) / fse3v(ji,jj,jk) 245 ub(ji,jj,jk) = un(ji,jj,jk) * zsshun * umask(ji,jj,jk) 246 vb(ji,jj,jk) = vn(ji,jj,jk) * zsshvn * vmask(ji,jj,jk) 247 zsshun = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 248 zsshvn = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 249 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun * umask(ji,jj,jk) 250 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn * vmask(ji,jj,jk) 251 END DO 252 END DO 253 ELSE ! Fixed levels 254 DO jj = 1, jpj 255 DO ji = 1, jpi 256 ! Euler (forward) time stepping 257 ub(ji,jj,jk) = un(ji,jj,jk) 258 vb(ji,jj,jk) = vn(ji,jj,jk) 259 un(ji,jj,jk) = ua(ji,jj,jk) 260 vn(ji,jj,jk) = va(ji,jj,jk) 261 END DO 262 END DO 263 ENDIF 160 264 ELSE 161 DO jj = 1, jpj ! caution: don't use (:,:) for this loop 162 DO ji = 1, jpi ! it causes optimization problems on NEC in auto-tasking 163 ! Leap-frog time stepping 164 ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 165 vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 166 un(ji,jj,jk) = ua(ji,jj,jk) 167 vn(ji,jj,jk) = va(ji,jj,jk) 168 END DO 169 END DO 265 IF( (lk_vvl .AND. .NOT. lk_dynspg_flt) ) THEN ! Varying levels 266 ! caution: don't use (:,:) for this loop 267 ! it causes optimization problems on NEC in auto-tasking 268 DO jj = 1, jpj 269 DO ji = 1, jpi 270 zsshun = umask(ji,jj,jk) / zfse3un(ji,jj,jk) 271 zsshvn = vmask(ji,jj,jk) / zfse3vn(ji,jj,jk) 272 ub(ji,jj,jk) = ( atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) & 273 & + atfp1 * un(ji,jj,jk) ) * zsshun 274 vb(ji,jj,jk) = ( atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) & 275 & + atfp1 * vn(ji,jj,jk) ) * zsshvn 276 zsshun = umask(ji,jj,jk) / zfse3ua(ji,jj,jk) 277 zsshvn = vmask(ji,jj,jk) / zfse3va(ji,jj,jk) 278 un(ji,jj,jk) = ua(ji,jj,jk) * zsshun 279 vn(ji,jj,jk) = va(ji,jj,jk) * zsshvn 280 END DO 281 END DO 282 ELSE ! Fixed levels 283 DO jj = 1, jpj 284 DO ji = 1, jpi 285 ! Leap-frog time stepping 286 ub(ji,jj,jk) = atfp * ( ub(ji,jj,jk) + ua(ji,jj,jk) ) + atfp1 * un(ji,jj,jk) 287 vb(ji,jj,jk) = atfp * ( vb(ji,jj,jk) + va(ji,jj,jk) ) + atfp1 * vn(ji,jj,jk) 288 un(ji,jj,jk) = ua(ji,jj,jk) 289 vn(ji,jj,jk) = va(ji,jj,jk) 290 END DO 291 END DO 292 ENDIF 170 293 ENDIF 171 294 ! ! =============== -
trunk/NEMO/OPA_SRC/DYN/dynspg_exp.F90
r575 r592 110 110 ENDIF 111 111 112 ! 0. Initialization 113 ! ----------------- 114 ! read or estimate sea surface height and vertically integrated velocities 115 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 116 z2dt = 2. * rdt ! time step: leap-frog 117 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 118 zraur = 1.e0 / rauw 119 120 121 ! 1. Surface pressure gradient (now) 122 ! ---------------------------- 123 DO jj = 2, jpjm1 124 DO ji = fs_2, fs_jpim1 ! vector opt. 125 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 126 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 127 END DO 128 END DO 129 130 ! 2. Add the surface pressure trend to the general trend 131 ! ------------------------------------------------------ 132 DO jk = 1, jpkm1 112 IF( .NOT. lk_vvl ) THEN 113 114 ! 0. Initialization 115 ! ----------------- 116 ! read or estimate sea surface height and vertically integrated velocities 117 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 118 z2dt = 2. * rdt ! time step: leap-frog 119 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 120 zraur = 1.e0 / rauw 121 122 ! 1. Surface pressure gradient (now) 123 ! ---------------------------- 133 124 DO jj = 2, jpjm1 134 125 DO ji = fs_2, fs_jpim1 ! vector opt. 135 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)136 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)126 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 127 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 137 128 END DO 138 129 END DO 139 END DO 130 131 ! 2. Add the surface pressure trend to the general trend 132 ! ------------------------------------------------------ 133 DO jk = 1, jpkm1 134 DO jj = 2, jpjm1 135 DO ji = fs_2, fs_jpim1 ! vector opt. 136 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 137 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 138 END DO 139 END DO 140 END DO 140 141 141 ! 3. Vertical integration of horizontal divergence of velocities 142 ! -------------------------------- 143 zhdiv(:,:) = 0.e0 144 DO jk = jpkm1, 1, -1 145 DO jj = 2, jpjm1 146 DO ji = fs_2, fs_jpim1 ! vector opt. 147 zhdiv(ji,jj) = zhdiv(ji,jj) + ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * un(ji ,jj ,jk) & 148 & - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 149 & + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * vn(ji ,jj ,jk) & 150 & - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 151 & / ( e1t(ji,jj) * e2t(ji,jj) ) 142 ! 3. Vertical integration of horizontal divergence of velocities 143 ! -------------------------------- 144 zhdiv(:,:) = 0.e0 145 DO jk = jpkm1, 1, -1 146 DO jj = 2, jpjm1 147 DO ji = fs_2, fs_jpim1 ! vector opt. 148 zhdiv(ji,jj) = zhdiv(ji,jj) + ( e2u(ji ,jj ) * fse3u(ji ,jj ,jk) * un(ji ,jj ,jk) & 149 & - e2u(ji-1,jj ) * fse3u(ji-1,jj ,jk) * un(ji-1,jj ,jk) & 150 & + e1v(ji ,jj ) * fse3v(ji ,jj ,jk) * vn(ji ,jj ,jk) & 151 & - e1v(ji ,jj-1) * fse3v(ji ,jj-1,jk) * vn(ji ,jj-1,jk) ) & 152 & / ( e1t(ji,jj) * e2t(ji,jj) ) 153 END DO 152 154 END DO 153 155 END DO 154 END DO155 156 156 157 #if defined key_obc 157 ! open boundaries (div must be zero behind the open boundary)158 ! mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column159 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east160 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west161 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north162 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south158 ! open boundaries (div must be zero behind the open boundary) 159 ! mpp remark: The zeroing of zhdiv can probably be extended to 1->jpi/jpj for the correct row/column 160 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east 161 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west 162 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 163 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south 163 164 #endif 164 165 165 166 ! 4. Sea surface elevation time stepping 167 ! -------------------------------------- 168 ! Euler (forward) time stepping, no time filter 169 IF( neuler == 0 .AND. kt == nit000 ) THEN 170 DO jj = 1, jpj 171 DO ji = 1, jpi 172 ! after free surface elevation 173 zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 174 ! swap of arrays 175 sshb(ji,jj) = sshn(ji,jj) 176 sshn(ji,jj) = zssha 177 END DO 178 END DO 179 ELSE 180 ! Leap-frog time stepping and time filter 181 DO jj = 1, jpj 182 DO ji = 1, jpi 183 ! after free surface elevation 184 zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 185 ! time filter and array swap 186 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj) 187 sshn(ji,jj) = zssha 188 END DO 189 END DO 166 ! 4. Sea surface elevation time stepping 167 ! -------------------------------------- 168 ! Euler (forward) time stepping, no time filter 169 IF( neuler == 0 .AND. kt == nit000 ) THEN 170 DO jj = 1, jpj 171 DO ji = 1, jpi 172 ! after free surface elevation 173 zssha = sshb(ji,jj) - rdt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 174 ! swap of arrays 175 sshb(ji,jj) = sshn(ji,jj) 176 sshn(ji,jj) = zssha 177 END DO 178 END DO 179 ELSE 180 ! Leap-frog time stepping and time filter 181 DO jj = 1, jpj 182 DO ji = 1, jpi 183 ! after free surface elevation 184 zssha = sshb(ji,jj) - z2dt * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) * tmask(ji,jj,1) 185 ! time filter and array swap 186 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj) 187 sshn(ji,jj) = zssha 188 END DO 189 END DO 190 ENDIF 191 192 ELSE !! Variable volume, ssh time-stepping already done 193 194 ! read or estimate sea surface height and vertically integrated velocities 195 IF( lk_obc ) CALL obc_dta_bt( kt, 0 ) 196 197 ! Sea surface elevation swap 198 ! ----------------------------- 199 ! 200 sshbb(:,:) = sshb(:,:) ! Needed for the dynamics time-stepping 201 ! 202 IF( neuler == 0 .AND. kt == nit000 ) THEN 203 ! No time filter 204 sshb(:,:) = sshn(:,:) 205 sshn(:,:) = ssha(:,:) 206 ELSE 207 ! Time filter 208 sshb(:,:) = atfp * ( sshb(:,:) + ssha(:,:) ) + atfp1 * sshn(:,:) 209 sshn(:,:) = ssha(:,:) 210 ENDIF 211 190 212 ENDIF 191 213 -
trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90
r575 r592 46 46 USE iom 47 47 USE restart ! only for lrst_oce 48 USE domvvl ! variable volume 48 49 49 50 IMPLICIT NONE … … 109 110 !! References : Roullet and Madec 1999, JGR. 110 111 !!--------------------------------------------------------------------- 112 !! * Modules used 113 USE oce , ONLY : zub => ta, & ! ta used as workspace 114 zvb => sa ! sa " " 115 111 116 INTEGER, INTENT( in ) :: kt ! ocean time-step index 112 117 INTEGER, INTENT( out ) :: kindic ! solver convergence flag (<0 if not converge) … … 114 119 INTEGER :: ji, jj, jk ! dummy loop indices 115 120 REAL(wp) :: z2dt, z2dtg, zraur, znugdt, & ! temporary scalars 116 & znurau, z ssha, zgcb, zbtd,& ! " "121 & znurau, zgcb, zbtd, & ! " " 117 122 & ztdgu, ztdgv ! " " 123 !! Variable volume 124 REAL(wp), DIMENSION(jpi,jpj) :: & 125 zsshub, zsshua, zsshvb, zsshva, zssha ! 2D workspace 126 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! 3D workspace 127 zfse3ub, zfse3ua, zfse3vb, zfse3va ! 3D workspace 118 128 !!---------------------------------------------------------------------- 119 129 ! … … 142 152 znurau = znugdt * zraur 143 153 144 ! Surface pressure gradient (now) 145 DO jj = 2, jpjm1 146 DO ji = fs_2, fs_jpim1 ! vector opt. 147 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 148 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 149 END DO 150 END DO 151 152 ! Add the surface pressure trend to the general trend 153 DO jk = 1, jpkm1 154 !! Explicit physics with thickness weighted updates 155 IF( lk_vvl ) THEN ! variable volume 156 157 DO jj = 1, jpjm1 158 DO ji = 1,jpim1 159 160 ! Sea Surface Height at u-point before 161 zsshub(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 162 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & 163 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) 164 165 ! Sea Surface Height at v-point before 166 zsshvb(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 167 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & 168 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) 169 170 ! Sea Surface Height at u-point after 171 zsshua(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 172 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 173 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 174 175 ! Sea Surface Height at v-point after 176 zsshva(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 177 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 178 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 179 180 END DO 181 END DO 182 183 ! Boundaries conditions 184 CALL lbc_lnk( zsshub, 'U', 1. ) ; CALL lbc_lnk( zsshvb, 'V', 1. ) 185 CALL lbc_lnk( zsshua, 'U', 1. ) ; CALL lbc_lnk( zsshva, 'V', 1. ) 186 187 ! Scale factors at before and after time step 188 ! ------------------------------------------- 189 DO jk = 1, jpkm1 190 zfse3ua(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshua(:,:) * muu(:,:,jk) ) 191 zfse3va(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshva(:,:) * muv(:,:,jk) ) 192 zfse3ub(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshub(:,:) * muu(:,:,jk) ) 193 zfse3vb(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshvb(:,:) * muv(:,:,jk) ) 194 END DO 195 196 ! Thickness weighting 197 ! ------------------- 198 ua(:,:,1:jpkm1) = ua(:,:,1:jpkm1) * fse3u(:,:,1:jpkm1) 199 va(:,:,1:jpkm1) = va(:,:,1:jpkm1) * fse3v(:,:,1:jpkm1) 200 201 zub(:,:,1:jpkm1) = ub(:,:,1:jpkm1) * zfse3ub(:,:,1:jpkm1) 202 zvb(:,:,1:jpkm1) = vb(:,:,1:jpkm1) * zfse3vb(:,:,1:jpkm1) 203 204 ! Evaluate the masked next velocity (effect of the additional force not included) 205 DO jk = 1, jpkm1 206 DO jj = 2, jpjm1 207 DO ji = fs_2, fs_jpim1 ! vector opt. 208 ua(ji,jj,jk) = ( zub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) /zfse3ua(ji,jj,jk) * umask(ji,jj,jk) 209 va(ji,jj,jk) = ( zvb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) /zfse3va(ji,jj,jk) * vmask(ji,jj,jk) 210 END DO 211 END DO 212 END DO 213 214 ELSE ! fixed volume 215 216 ! Surface pressure gradient (now) 154 217 DO jj = 2, jpjm1 155 218 DO ji = fs_2, fs_jpim1 ! vector opt. 156 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 157 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 158 END DO 159 END DO 160 END DO 161 162 ! Evaluate the masked next velocity (effect of the additional force not included) 163 DO jk = 1, jpkm1 164 DO jj = 2, jpjm1 165 DO ji = fs_2, fs_jpim1 ! vector opt. 166 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 167 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 168 END DO 169 END DO 170 END DO 219 spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj) 220 spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj) 221 END DO 222 END DO 223 224 ! Add the surface pressure trend to the general trend 225 DO jk = 1, jpkm1 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 ! vector opt. 228 ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj) 229 va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj) 230 END DO 231 END DO 232 END DO 233 234 ! Evaluate the masked next velocity (effect of the additional force not included) 235 DO jk = 1, jpkm1 236 DO jj = 2, jpjm1 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 ua(ji,jj,jk) = ( ub(ji,jj,jk) + z2dt * ua(ji,jj,jk) ) * umask(ji,jj,jk) 239 va(ji,jj,jk) = ( vb(ji,jj,jk) + z2dt * va(ji,jj,jk) ) * vmask(ji,jj,jk) 240 END DO 241 END DO 242 END DO 243 244 ENDIF 171 245 172 246 #if defined key_obc … … 221 295 CALL lbc_lnk( spgu, 'U', -1. ) 222 296 CALL lbc_lnk( spgv, 'V', -1. ) 297 298 IF( lk_vvl ) CALL sol_mat( kt ) 223 299 224 300 ! Right hand side of the elliptic equation and first guess … … 334 410 ! Sea surface elevation time stepping 335 411 ! ----------------------------------- 412 IF( lk_vvl ) THEN ! after free surface elevation 413 zssha(:,:) = ssha(:,:) 414 ELSE 415 zssha(:,:) = sshb(:,:) + z2dt * ( wn(:,:,1) - emp(:,:) * zraur ) * tmask(:,:,1) 416 ENDIF 417 336 418 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping, no time filter 337 DO jj = 1, jpj 338 DO ji = 1, jpi 339 ! after free surface elevation 340 zssha = sshb(ji,jj) + rdt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1) 341 ! swap of arrays 342 sshb(ji,jj) = sshn(ji,jj) 343 sshn(ji,jj) = zssha 344 END DO 345 END DO 419 ! swap of arrays 420 sshb(:,:) = sshn (:,:) 421 sshn(:,:) = zssha(:,:) 346 422 ELSE ! Leap-frog time stepping and time filter 347 DO jj = 1, jpj 348 DO ji = 1, jpi 349 ! after free surface elevation 350 zssha = sshb(ji,jj) + z2dt * ( wn(ji,jj,1) - emp(ji,jj) * zraur ) * tmask(ji,jj,1) 351 ! time filter and array swap 352 sshb(ji,jj) = atfp * ( sshb(ji,jj) + zssha ) + atfp1 * sshn(ji,jj) 353 sshn(ji,jj) = zssha 354 END DO 355 END DO 423 ! time filter and array swap 424 sshb(:,:) = atfp * ( sshb(:,:) + zssha(:,:) ) + atfp1 * sshn(:,:) 425 sshn(:,:) = zssha(:,:) 356 426 ENDIF 357 427 -
trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90
r367 r592 39 39 #endif 40 40 41 #if defined key_dynspg_ts || defined key_ esopa41 #if defined key_dynspg_ts || defined key_vvl || defined key_esopa 42 42 !! Time splitting variables 43 43 !! ------------------------ -
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r575 r592 34 34 USE iom 35 35 USE restart ! only for lrst_oce 36 USE domvvl ! variable volume 36 37 37 38 IMPLICIT NONE … … 102 103 zsshb_e, zub_e, zvb_e, & ! " " 103 104 zun_e, zvn_e ! " " 105 !! Variable volume 106 REAL(wp), DIMENSION(jpi,jpj) :: & 107 zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e ! 2D workspace 108 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace 104 109 !!---------------------------------------------------------------------- 105 110 … … 137 142 ENDIF 138 143 ! 144 ENDIF 145 146 ! Substract the surface pressure gradient from the momentum 147 ! --------------------------------------------------------- 148 IF( lk_vvl ) THEN ! Variable volume 149 150 ! 0. Local initialization 151 ! ----------------------- 152 zspgu_1(:,:) = 0.e0 153 zspgv_1(:,:) = 0.e0 154 155 ! 1. Surface pressure gradient (now) 156 ! ---------------------------- 157 DO jj = 2, jpjm1 158 DO ji = fs_2, fs_jpim1 ! vector opt. 159 zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 160 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj) 161 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 162 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 163 END DO 164 END DO 165 166 ! 2. Substract the surface pressure trend from the general trend 167 ! ------------------------------------------------------ 168 DO jk = 1, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 ! vector opt. 171 ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj) 172 va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) 173 END DO 174 END DO 175 END DO 176 139 177 ENDIF 140 178 … … 269 307 zua_b (:,:) = un_b (:,:) 270 308 zva_b (:,:) = vn_b (:,:) 309 IF( lk_vvl ) THEN 310 zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point 311 zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point 312 hu_e(:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 313 hv_e(:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 314 zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point 315 zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point 316 ENDIF 271 317 272 318 ! set ssh corrections to 0 … … 330 376 DO ji = fs_2, fs_jpim1 ! vector opt. 331 377 ! surface pressure gradient 332 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 333 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 378 IF( lk_vvl) THEN 379 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 380 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 381 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 382 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 383 ELSE 384 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 385 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 386 ENDIF 334 387 ! energy conserving formulation for planetary vorticity term 335 388 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) … … 349 402 DO ji = fs_2, fs_jpim1 ! vector opt. 350 403 ! surface pressure gradient 351 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 352 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 404 IF( lk_vvl) THEN 405 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 406 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 407 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 408 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 409 ELSE 410 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 411 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 412 ENDIF 353 413 ! enstrophy conserving formulation for planetary vorticity term 354 414 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & … … 369 429 DO ji = fs_2, fs_jpim1 ! vector opt. 370 430 ! surface pressure gradient 371 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 372 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 431 IF( lk_vvl) THEN 432 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & 433 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) 434 zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & 435 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) 436 ELSE 437 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 438 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 439 ENDIF 373 440 ! energy/enstrophy conserving formulation for planetary vorticity term 374 441 zcubt = + zfac25 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & … … 403 470 ! Time filter and swap of dynamics arrays 404 471 ! --------------------------------------- 405 IF( neuler == 0 .AND. kt == nit000) THEN ! Euler (forward) time stepping472 IF( neuler == 0 .AND. jit == 1 ) THEN ! Euler (forward) time stepping 406 473 zsshb_e(:,:) = sshn_e(:,:) 407 474 zub_e (:,:) = zun_e (:,:) … … 419 486 ENDIF 420 487 488 IF( lk_vvl ) THEN ! Variable volume 489 490 ! Sea surface elevation 491 ! --------------------- 492 DO jj = 1, jpjm1 493 DO ji = 1,jpim1 494 495 ! Sea Surface Height at u-point before 496 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 497 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 498 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 499 500 ! Sea Surface Height at v-point before 501 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 502 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 503 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 504 505 END DO 506 END DO 507 508 ! Boundaries conditions 509 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 510 511 ! Scale factors at before and after time step 512 ! ------------------------------------------- 513 DO jk = 1, jpkm1 514 zfse3un_e(:,:,jk) = fsve3u(:,:,jk) * ( 1 + zsshun_e(:,:) * muu(:,:,jk) ) 515 zfse3vn_e(:,:,jk) = fsve3v(:,:,jk) * ( 1 + zsshvn_e(:,:) * muv(:,:,jk) ) 516 END DO 517 518 ! Ocean depth at U- and V-points 519 hu_e(:,:) = 0.e0 520 hv_e(:,:) = 0.e0 521 522 DO jk = 1, jpk 523 hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 524 hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 525 END DO 526 527 ENDIF ! End variable volume case 528 421 529 ! ! ==================== ! 422 530 END DO ! end loop ! … … 444 552 ! Horizontal divergence of time averaged barotropic transports 445 553 !------------------------------------------------------------- 446 DO jj = 2, jpjm1 447 DO ji = fs_2, fs_jpim1 ! vector opt. 448 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) & 449 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) & 450 & / ( e1t(ji,jj) * e2t(ji,jj) ) 451 END DO 452 END DO 453 454 #if defined key_obc 554 IF( .NOT. lk_vvl ) THEN 555 DO jj = 2, jpjm1 556 DO ji = fs_2, fs_jpim1 ! vector opt. 557 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) & 558 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) & 559 & / ( e1t(ji,jj) * e2t(ji,jj) ) 560 END DO 561 END DO 562 ENDIF 563 564 #if defined key_obc && ! defined key_vvl 455 565 ! open boundaries (div must be zero behind the open boundary) 456 566 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column … … 463 573 ! sea surface height 464 574 !------------------- 465 sshb(:,:) = sshn(:,:) 466 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 575 IF( lk_vvl ) THEN 576 sshbb(:,:) = sshb(:,:) 577 sshb (:,:) = sshn(:,:) 578 sshn (:,:) = ssha(:,:) 579 ELSE 580 sshb (:,:) = sshn(:,:) 581 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 582 ENDIF 467 583 468 584 ! ... Boundary conditions on sshn … … 477 593 ! ------------------------------------------ 478 594 sshb_b(:,:) = sshn_b (:,:) 595 IF ( kt == nit000 ) zssha_b(:,:) = sshn(:,:) 479 596 sshn_b(:,:) = zssha_b(:,:) 480 597 un_b (:,:) = zua_b (:,:) -
trunk/NEMO/OPA_SRC/DYN/wzvmod.F90
r455 r592 4 4 !! Ocean diagnostic variable : vertical velocity 5 5 !!============================================================================== 6 6 !! History : 5.0 ! 90-10 (C. Levy, G. Madec) Original code 7 !! 7.0 ! 96-01 (G. Madec) Statement function for e3 8 !! 8.5 ! 02-07 (G. Madec) Free form, F90 7 9 !!---------------------------------------------------------------------- 8 10 !! wzv : Compute the vertical velocity … … 14 16 USE prtctl ! Print control 15 17 18 USE domvvl ! Variable volume 19 USE phycst 20 USE ocesbc ! ocean surface boundary condition 21 USE lbclnk ! ocean lateral boundary condition (or mpp link) 22 16 23 IMPLICIT NONE 17 24 PRIVATE 18 25 19 26 !! * Routine accessibility 20 PUBLIC wzv ! routine called by step.F90 and inidtr.F9027 PUBLIC wzv ! routine called by step.F90 and inidtr.F90 21 28 22 29 !! * Substitutions 23 30 # include "domzgr_substitute.h90" 31 !!---------------------------------------------------------------------- 32 !! OPA 9.0 , LOCEAN-IPSL (2005) 33 !! $Header$ 34 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 24 35 !!---------------------------------------------------------------------- 25 36 … … 40 51 !! velocity is computed by integrating the horizontal divergence 41 52 !! from the bottom to the surface. 42 !! 43 !! in r egid-lid case, w=0 at the sea surface.53 !! The boundary conditions are w=0 at the bottom (no flux) and, 54 !! in rigid-lid case, w=0 at the sea surface. 44 55 !! 45 56 !! ** action : wn array : the now vertical velocity 46 !!47 !! History :48 !! 5.0 ! 90-10 (C. Levy, G. Madec) Original code49 !! 7.0 ! 96-01 (G. Madec) Statement function for e350 !! 8.5 ! 02-07 (G. Madec) Free form, F9051 57 !!---------------------------------------------------------------------- 52 58 !! * Arguments … … 55 61 !! * Local declarations 56 62 INTEGER :: jj, jk ! dummy loop indices 57 !!----------------------------------------------------------------------58 !! OPA 9.0 , LOCEAN-IPSL (2005)59 !! $Header$60 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt61 63 !!---------------------------------------------------------------------- 62 64 … … 103 105 !! 104 106 !! ** action : wn array : the now vertical velocity 105 !!106 !! History :107 !! 9.0 ! 02-07 (G. Madec) Vector optimization108 107 !!---------------------------------------------------------------------- 109 108 !! * Arguments … … 111 110 112 111 !! * Local declarations 113 INTEGER :: jk ! dummy loop indices 114 !!---------------------------------------------------------------------- 115 !! OPA 8.5, LODYC-IPSL (2002) 112 INTEGER :: jk ! dummy loop indices 113 !! Variable volume 114 INTEGER :: ji, jj ! dummy loop indices 115 REAL(wp) :: z2dt, zraur ! temporary scalar 116 REAL(wp), DIMENSION (jpi,jpj) :: zssha, zun, zvn, zhdiv 116 117 !!---------------------------------------------------------------------- 117 118 … … 125 126 ENDIF 126 127 127 ! Computation from the bottom 128 DO jk = jpkm1, 1, -1 129 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 130 END DO 128 IF( lk_vvl ) THEN ! Variable volume 129 ! 130 z2dt = 2. * rdt ! time step: leap-frog 131 IF( neuler == 0 .AND. kt == nit000 ) z2dt = rdt ! time step: Euler if restart from rest 132 zraur = 1. / rauw 133 134 ! Vertically integrated quantities 135 ! -------------------------------- 136 zun(:,:) = 0.e0 137 zvn(:,:) = 0.e0 138 ! 139 DO jk = 1, jpkm1 ! Vertically integrated transports (now) 140 zun(:,:) = zun(:,:) + fse3u(:,:,jk) * un(:,:,jk) 141 zvn(:,:) = zvn(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 142 END DO 143 144 ! Horizontal divergence of barotropic transports 145 !-------------------------------------------------- 146 DO jj = 2, jpjm1 147 DO ji = 2, jpim1 ! vector opt. 148 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun(ji ,jj ) & 149 & - e2u(ji-1,jj ) * zun(ji-1,jj ) & 150 & + e1v(ji ,jj ) * zvn(ji ,jj ) & 151 & - e1v(ji ,jj-1) * zvn(ji ,jj-1) ) & 152 & / ( e1t(ji,jj) * e2t(ji,jj) ) 153 END DO 154 END DO 155 156 #if defined key_obc && ( key_dynspg_exp || key_dynspg_ts ) 157 ! open boundaries (div must be zero behind the open boundary) 158 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 159 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 160 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 161 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 162 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 163 #endif 164 165 CALL lbc_lnk( zhdiv, 'T', 1. ) 166 167 ! Sea surface elevation time stepping 168 ! ----------------------------------- 169 zssha(:,:) = sshb(:,:) - z2dt * ( zraur * emp(:,:) + zhdiv(:,:) ) * tmask(:,:,1) 170 171 ! Vertical velocity computed from bottom 172 ! -------------------------------------- 173 DO jk = jpkm1, 1, -1 174 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) & 175 & - ( zssha(:,:) - sshb(:,:) ) * fsve3t(:,:,jk) * mut(:,:,jk) / z2dt 176 END DO 177 178 ELSE ! Fixed volume 179 180 ! Vertical velocity computed from bottom 181 ! -------------------------------------- 182 DO jk = jpkm1, 1, -1 183 wn(:,:,jk) = wn(:,:,jk+1) - fse3t(:,:,jk) * hdivn(:,:,jk) 184 END DO 185 186 ENDIF 131 187 132 188 IF(ln_ctl) CALL prt_ctl(tab3d_1=wn, clinfo1=' w**2 - : ', mask1=wn)
Note: See TracChangeset
for help on using the changeset viewer.
