Changeset 592
- Timestamp:
- 2007-02-09T10:15:25+01:00 (17 years ago)
- Location:
- trunk/NEMO/OPA_SRC
- Files:
-
- 1 added
- 23 edited
-
DOM/dom_oce.F90 (modified) (1 diff)
-
DOM/domain.F90 (modified) (2 diffs)
-
DOM/domstp.F90 (modified) (1 diff)
-
DOM/domvvl.F90 (added)
-
DOM/domzgr.F90 (modified) (4 diffs)
-
DOM/domzgr_substitute.h90 (modified) (1 diff)
-
DYN/dynhpg.F90 (modified) (5 diffs)
-
DYN/dynnxt.F90 (modified) (5 diffs)
-
DYN/dynspg_exp.F90 (modified) (1 diff)
-
DYN/dynspg_flt.F90 (modified) (6 diffs)
-
DYN/dynspg_oce.F90 (modified) (1 diff)
-
DYN/dynspg_ts.F90 (modified) (12 diffs)
-
DYN/wzvmod.F90 (modified) (7 diffs)
-
LDF/ldfslp.F90 (modified) (1 diff)
-
TRA/traadv_cen2.F90 (modified) (1 diff)
-
TRA/traadv_muscl.F90 (modified) (1 diff)
-
TRA/traadv_muscl2.F90 (modified) (1 diff)
-
TRA/traadv_tvd.F90 (modified) (1 diff)
-
TRA/traadv_ubs.F90 (modified) (1 diff)
-
TRA/trabbc.F90 (modified) (8 diffs)
-
TRA/tranxt.F90 (modified) (6 diffs)
-
TRA/trasbc.F90 (modified) (1 diff)
-
TRA/trazdf.F90 (modified) (2 diffs)
-
TRA/trazdf_imp.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/dom_oce.F90
r454 r592 113 113 e3w , e3uw !: W--UW points (m) 114 114 #endif 115 #if defined key_vvl 116 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .TRUE. !: variable grid flag 117 118 !! All coordinates 119 !! --------------- 120 REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) :: & !: 121 gdep3w_1 , & !: depth of T-points (sum of e3w) (m) 122 gdept_1, gdepw_1, & !: analytical depth at T-W points (m) 123 e3v_1 , e3f_1 , & !: analytical vertical scale factors at V--F 124 e3t_1 , e3u_1 , & !: T--U points (m) 125 e3vw_1 , & !: analytical vertical scale factors at VW-- 126 e3w_1 , e3uw_1 !: W--UW points (m) 127 #else 128 LOGICAL, PUBLIC, PARAMETER :: lk_vvl = .FALSE. !: fixed grid flag 129 #endif 115 130 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & !: 116 131 hur, hvr, & !: inverse of u and v-points ocean depth (1/m) -
trunk/NEMO/OPA_SRC/DOM/domain.F90
r531 r592 27 27 USE domwri ! domain: write the meshmask file 28 28 USE closea ! closed sea or lake (dom_clo routine) 29 USE domvvl ! variable volume 29 30 30 31 IMPLICIT NONE … … 91 92 CALL dom_msk ! Masks 92 93 94 IF( lk_vvl ) CALL dom_vvl_ini ! Vertical variable mesh 93 95 94 96 ! Local depth or Inverse of the local depth of the water column at u- and v-points -
trunk/NEMO/OPA_SRC/DOM/domstp.F90
r474 r592 88 88 IF(lwp) WRITE(numout,*)' accelerating the convergence' 89 89 IF(lwp) WRITE(numout,*)' dynamics time step = ', rdt/3600., ' hours' 90 IF( ln_sco .AND. rdtmin /= rdtmax ) & 91 & CALL ctl_stop ( ' depth dependent acceleration of convergence not implemented in s-coordinates' ) 90 IF( ln_sco .AND. rdtmin /= rdtmax .AND. lk_vvl ) & 91 & CALL ctl_stop ( ' depth dependent acceleration of convergence not implemented in s-coordinates & 92 & nor in variable volume' ) 92 93 IF(lwp) WRITE(numout,*)' tracers time step : dt (hours) level' 93 94 -
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r528 r592 337 337 idta(:,:) = jpkm1 ! flat basin 338 338 zdta(:,:) = gdepw_0(jpk) 339 h_oce = gdepw_0(jpk) 339 340 340 341 ELSE ! bump … … 342 343 IF(lwp) WRITE(numout,*) ' bathymetry field: flat basin with a bump' 343 344 344 ii_bump = jpidta / 3 + 3! i-index of the bump center345 ii_bump = jpidta / 2 ! i-index of the bump center 345 346 ij_bump = jpjdta / 2 ! j-index of the bump center 346 r_bump = 6 ! bump radius (index)347 h_bump = 240.e0 ! bump height (meters)347 r_bump = 50000.e0 ! bump radius (meters) 348 h_bump = 2700.e0 ! bump height (meters) 348 349 h_oce = gdepw_0(jpk) ! background ocean depth (meters) 349 350 IF(lwp) WRITE(numout,*) ' bump characteristics: ' … … 355 356 DO jj = 1, jpjdta 356 357 DO ji = 1, jpidta 357 zi = FLOAT( ji - ii_bump ) / r_bump358 zj = FLOAT( jj - ij_bump ) / r_bump358 zi = FLOAT( ji - ii_bump ) * ppe1_m / r_bump 359 zj = FLOAT( jj - ij_bump ) * ppe2_m / r_bump 359 360 zdta(ji,jj) = h_oce - h_bump * EXP( -( zi*zi + zj*zj ) ) 360 361 END DO … … 392 393 idta( 1 , : ) = ih ; zdta( 1 , : ) = zh 393 394 idta(jpidta, : ) = ih ; zdta(jpidta, : ) = zh 394 ENDIF395 396 ! EEL R5 configuration with east and west open boundaries.397 ! Two rows of zeroes are needed at the south and north for OBCs398 399 IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN400 ih = 0 ; zh = 0.e0401 IF( ln_sco ) zh = jpkm1 ; IF( ln_sco ) zh = h_oce402 idta( : , 2 ) = jpkm1 ; zdta( : , 2 ) = h_oce403 idta( : ,jpjdta-1) = jpkm1 ; zdta( : ,jpjdta-1) = h_oce404 395 ENDIF 405 396 -
trunk/NEMO/OPA_SRC/DOM/domzgr_substitute.h90
r454 r592 40 40 # define fse3vw(i,j,k) e3vw(i,j,k) 41 41 #endif 42 43 #if defined key_vvl 44 # define fsvdept(i,j,k) gdept_1(i,j,k) 45 46 # define fsvdepw(i,j,k) gdepw_1(i,j,k) 47 # define fsvde3w(i,j,k) gdep3w_1(i,j,k) 48 49 # define fsve3t(i,j,k) e3t_1(i,j,k) 50 # define fsve3u(i,j,k) e3u_1(i,j,k) 51 # define fsve3v(i,j,k) e3v_1(i,j,k) 52 # define fsve3f(i,j,k) e3f_1(i,j,k) 53 54 # define fsve3w(i,j,k) e3w_1(i,j,k) 55 # define fsve3uw(i,j,k) e3uw_1(i,j,k) 56 # define fsve3vw(i,j,k) e3vw_1(i,j,k) 57 #else 58 # define fsvdept(i,j,k) gdept(i,j,k) 59 60 # define fsvdepw(i,j,k) gdepw(i,j,k) 61 # define fsvde3w(i,j,k) gdep3w(i,j,k) 62 63 # define fsve3t(i,j,k) e3t(i,j,k) 64 # define fsve3u(i,j,k) e3u(i,j,k) 65 # define fsve3v(i,j,k) e3v(i,j,k) 66 # define fsve3f(i,j,k) e3f(i,j,k) 67 68 # define fsve3w(i,j,k) e3w(i,j,k) 69 # define fsve3uw(i,j,k) e3uw(i,j,k) 70 # define fsve3vw(i,j,k) e3vw(i,j,k) 71 #endif 72 -
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) -
trunk/NEMO/OPA_SRC/LDF/ldfslp.F90
r461 r592 723 723 wslpjml(:,:) = 0.e0 724 724 725 IF( ln_traldf_hor .OR. ln_dynldf_hor) THEN725 IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (lk_vvl .AND. ln_rstart) ) THEN 726 726 IF(lwp) THEN 727 727 WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' -
trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90
r503 r592 304 304 305 305 ! Surface value 306 IF( lk_dynspg_rl ) THEN307 ! rigid lid : flux set to zero306 IF( lk_dynspg_rl .OR. lk_vvl ) THEN 307 ! rigid lid or variable volume: flux set to zero 308 308 zwx(:,:, 1 ) = 0.e0 ; zwy(:,:, 1 ) = 0.e0 309 309 ELSE -
trunk/NEMO/OPA_SRC/TRA/traadv_muscl.F90
r503 r592 355 355 END DO 356 356 ! surface values 357 IF( lk_dynspg_rl ) THEN ! rigid lid: flux set to zero357 IF( lk_dynspg_rl .OR. lk_vvl) THEN ! rigid lid or variable volume: flux set to zero 358 358 zt1(:,:, 1 ) = 0.e0 359 359 zs1(:,:, 1 ) = 0.e0 -
trunk/NEMO/OPA_SRC/TRA/traadv_muscl2.F90
r503 r592 428 428 429 429 ! surface values 430 IF( lk_dynspg_rl ) THEN ! rigid lid: flux set to zero430 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or variable volume: flux set to zero 431 431 zt1(:,:, 1 ) = 0.e0 432 432 zs1(:,:, 1 ) = 0.e0 -
trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90
r503 r592 130 130 ! upstream tracer flux in the k direction 131 131 ! Surface value 132 IF( lk_dynspg_rl ) THEN ! rigid lid: flux set to zero132 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid or variable volume: flux set to zero 133 133 ztw(:,:,1) = 0.e0 134 134 zsw(:,:,1) = 0.e0 -
trunk/NEMO/OPA_SRC/TRA/traadv_ubs.F90
r519 r592 373 373 ! ------------------------------------------------------------------- 374 374 ! Surface value 375 IF( lk_dynspg_rl ) THEN ! rigid lid : flux set to zero375 IF( lk_dynspg_rl .OR. lk_vvl ) THEN ! rigid lid : flux set to zero 376 376 ztw(:,:,1) = 0.e0 377 377 zsw(:,:,1) = 0.e0 -
trunk/NEMO/OPA_SRC/TRA/trabbc.F90
r541 r592 37 37 38 38 INTEGER , DIMENSION(jpi,jpj) :: nbotlevt ! ocean bottom level index at T-pt 39 REAL(wp), DIMENSION(jpi,jpj) :: qgh_trd 39 REAL(wp), DIMENSION(jpi,jpj) :: qgh_trd0 ! geothermal heating trend 40 40 41 41 !! * Substitutions … … 80 80 INTEGER :: ji, jj ! dummy loop indices 81 81 #endif 82 REAL(wp) :: zqgh_trd ! geothermal heat flux trend 82 83 !!---------------------------------------------------------------------- 83 84 … … 96 97 #if defined key_vectopt_loop && ! defined key_mpp_omp 97 98 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 98 ta(ji,1,nbotlevt(ji,1)) = ta(ji,1,nbotlevt(ji,1)) + qgh_trd(ji,1) 99 zqgh_trd = ro0cpr * qgh_trd0(ji,1) / fse3t(ji,1,nbotlevt(ji,1) ) 100 ta(ji,1,nbotlevt(ji,1)) = ta(ji,1,nbotlevt(ji,1)) + zqgh_trd 99 101 END DO 100 102 #else 101 103 DO jj = 2, jpjm1 102 104 DO ji = 2, jpim1 103 ta(ji,jj,nbotlevt(ji,jj)) = ta(ji,jj,nbotlevt(ji,jj)) + qgh_trd(ji,jj) 105 zqgh_trd = ro0cpr * qgh_trd0(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj)) 106 ta(ji,jj,nbotlevt(ji,jj)) = ta(ji,jj,nbotlevt(ji,jj)) + zqgh_trd 104 107 END DO 105 108 END DO … … 131 134 !! - NetCDF file : geothermal_heating.nc ( if necessary ) 132 135 !! 133 !! ** Action : - compute the heat geothermal trend qgh_trd136 !! ** Action : - read/fix the geothermal heat qgh_trd0 134 137 !! - compute the bottom ocean level nbotlevt 135 138 !!---------------------------------------------------------------------- … … 162 165 END DO 163 166 164 165 167 SELECT CASE ( ngeo_flux ) ! initialization of geothermal heat flux 166 168 ! … … 171 173 CASE ( 1 ) ! constant flux 172 174 IF(lwp) WRITE(numout,*) ' *** constant heat flux = ', ngeo_flux_const 173 qgh_trd (:,:) = ngeo_flux_const175 qgh_trd0(:,:) = ngeo_flux_const 174 176 ! 175 177 CASE ( 2 ) ! variable geothermal heat flux … … 178 180 IF(lwp) WRITE(numout,*) ' *** variable geothermal heat flux' 179 181 CALL iom_open ( 'geothermal_heating.nc', inum ) 180 CALL iom_get ( inum, jpdom_data, 'heatflow', qgh_trd )182 CALL iom_get ( inum, jpdom_data, 'heatflow', qgh_trd0 ) 181 183 CALL iom_close (inum) 182 184 ! 183 qgh_trd (:,:) = qgh_trd(:,:) * 1.e-3 ! conversion in W/m2185 qgh_trd0(:,:) = qgh_trd0(:,:) * 1.e-3 ! conversion in W/m2 184 186 ! 185 187 CASE DEFAULT … … 189 191 END SELECT 190 192 191 ! geothermal heat flux trend 192 193 SELECT CASE ( ngeo_flux ) 194 ! 195 CASE ( 1:2 ) ! geothermal heat flux 196 #if defined key_vectopt_loop && ! defined key_mpp_omp 197 DO ji = 1, jpij ! vector opt. (forced unrolling) 198 qgh_trd(ji,1) = ro0cpr * qgh_trd(ji,1) / fse3t(ji,1,nbotlevt(ji,1) ) 199 END DO 200 #else 201 DO jj = 1, jpj 202 DO ji = 1, jpi 203 qgh_trd(ji,jj) = ro0cpr * qgh_trd(ji,jj) / fse3t(ji,jj,nbotlevt(ji,jj)) 204 END DO 205 END DO 206 #endif 207 END SELECT 208 ! 193 209 194 END SUBROUTINE tra_bbc_init 210 195 -
trunk/NEMO/OPA_SRC/TRA/tranxt.F90
r503 r592 30 30 USE agrif_opa_interp 31 31 32 USE ocesbc ! ocean surface boundary condition 33 USE domvvl ! variable volume 34 USE dynspg_oce ! surface pressure gradient variables 35 USE phycst 36 32 37 IMPLICIT NONE 33 38 PRIVATE … … 35 40 !! * Routine accessibility 36 41 PUBLIC tra_nxt ! routine called by step.F90 42 43 REAL(wp) :: vemp ! total amount of volume added or removed by E-P forcing 44 45 !! * Substitutions 46 # include "domzgr_substitute.h90" 37 47 !!---------------------------------------------------------------------- 38 48 !! OPA 9.0 , LOCEAN-IPSL (2006) … … 79 89 REAL(wp) :: zt, zs ! temporary scalars 80 90 REAL(wp) :: zfact ! temporary scalar 91 !! Variable volume 92 REAL(wp) :: zssh ! temporary scalars 93 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3tb, zfse3tn, zfse3ta ! 3D workspace 94 81 95 !!---------------------------------------------------------------------- 96 97 !! Explicit physics with thickness weighted updates 98 IF( lk_vvl .AND. ln_zdfexp ) THEN 99 100 ! Scale factors at before and after time step 101 ! ------------------------------------------- 102 DO jk = 1, jpkm1 103 zfse3tb(:,:,jk) = fsve3t(:,:,jk) * ( 1 + sshb(:,:) * mut(:,:,jk) ) 104 zfse3ta(:,:,jk) = fsve3t(:,:,jk) * ( 1 + ssha(:,:) * mut(:,:,jk) ) 105 END DO 106 107 ! Asselin filtered scale factor at now time step 108 ! ---------------------------------------------- 109 IF( (neuler == 0 .AND. kt == nit000) .OR. lk_dynspg_ts ) THEN 110 zfse3tn(:,:,:) = fse3t(:,:,:) 111 ELSE 112 DO jk = 1, jpkm1 113 DO jj = 1, jpj 114 DO ji = 1, jpi 115 zssh = atfp * ( sshb(ji,jj) + ssha(ji,jj) ) + atfp1 * sshn(ji,jj) 116 zfse3tn(ji,jj,jk) = fsve3t(ji,jj,jk) * ( 1 + zssh * mut(ji,jj,jk) ) 117 END DO 118 END DO 119 END DO 120 ENDIF 121 122 ! Thickness weighting 123 ! ------------------- 124 ta(:,:,1:jpkm1) = ta(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1) 125 sa(:,:,1:jpkm1) = sa(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1) 126 127 tn(:,:,1:jpkm1) = tn(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1) 128 sn(:,:,1:jpkm1) = sn(:,:,1:jpkm1) * fse3t (:,:,1:jpkm1) 129 130 tb(:,:,1:jpkm1) = tb(:,:,1:jpkm1) * zfse3tb(:,:,1:jpkm1) 131 sb(:,:,1:jpkm1) = sb(:,:,1:jpkm1) * zfse3tb(:,:,1:jpkm1) 132 133 ENDIF 82 134 83 135 IF( l_trdtra ) THEN … … 85 137 ztrds(:,:,jpk) = 0.e0 86 138 ENDIF 139 87 140 ! 0. Lateral boundary conditions on ( ta, sa ) (T-point, unchanged sign) 88 141 ! ---------------------------------============ … … 165 218 ELSE ! Default case 166 219 IF( neuler == 0 .AND. kt == nit000 ) THEN 167 DO jj = 1, jpj 168 DO ji = 1, jpi 169 tb(ji,jj,jk) = tn(ji,jj,jk) 170 sb(ji,jj,jk) = sn(ji,jj,jk) 171 tn(ji,jj,jk) = ta(ji,jj,jk) 172 sn(ji,jj,jk) = sa(ji,jj,jk) 173 END DO 174 END DO 220 IF( (lk_vvl .AND. ln_zdfexp) ) THEN ! Varying levels 221 DO jj = 1, jpj 222 DO ji = 1, jpi 223 zssh = tmask(ji,jj,jk) / fse3t(ji,jj,jk) 224 tb(ji,jj,jk) = tn(ji,jj,jk) * zssh * tmask(ji,jj,jk) 225 sb(ji,jj,jk) = sn(ji,jj,jk) * zssh * tmask(ji,jj,jk) 226 zssh = tmask(ji,jj,jk) / zfse3ta(ji,jj,jk) 227 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh * tmask(ji,jj,jk) 228 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh * tmask(ji,jj,jk) 229 END DO 230 END DO 231 ELSE ! Fixed levels 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 tb(ji,jj,jk) = tn(ji,jj,jk) 235 sb(ji,jj,jk) = sn(ji,jj,jk) 236 tn(ji,jj,jk) = ta(ji,jj,jk) 237 sn(ji,jj,jk) = sa(ji,jj,jk) 238 END DO 239 END DO 240 ENDIF 175 241 IF( l_trdtra ) THEN 176 242 ztrdt(:,:,jk) = 0.e0 … … 186 252 END DO 187 253 END IF 188 DO jj = 1, jpj 189 DO ji = 1, jpi 190 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 191 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 192 tn(ji,jj,jk) = ta(ji,jj,jk) 193 sn(ji,jj,jk) = sa(ji,jj,jk) 194 END DO 195 END DO 254 IF( (lk_vvl .AND. ln_zdfexp) ) THEN ! Varying levels 255 DO jj = 1, jpj 256 DO ji = 1, jpi 257 zssh = tmask(ji,jj,jk) / zfse3tn(ji,jj,jk) 258 tb(ji,jj,jk) = ( atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) & 259 & + atfp1 * tn(ji,jj,jk) ) * zssh 260 sb(ji,jj,jk) = ( atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) & 261 & + atfp1 * sn(ji,jj,jk) ) * zssh 262 zssh = tmask(ji,jj,1) / zfse3ta(ji,jj,jk) 263 tn(ji,jj,jk) = ta(ji,jj,jk) * zssh 264 sn(ji,jj,jk) = sa(ji,jj,jk) * zssh 265 END DO 266 END DO 267 ELSE ! Fixed levels or first varying level 268 DO jj = 1, jpj 269 DO ji = 1, jpi 270 tb(ji,jj,jk) = atfp * ( tb(ji,jj,jk) + ta(ji,jj,jk) ) + atfp1 * tn(ji,jj,jk) 271 sb(ji,jj,jk) = atfp * ( sb(ji,jj,jk) + sa(ji,jj,jk) ) + atfp1 * sn(ji,jj,jk) 272 tn(ji,jj,jk) = ta(ji,jj,jk) 273 sn(ji,jj,jk) = sa(ji,jj,jk) 274 END DO 275 END DO 276 ENDIF 196 277 ENDIF 197 278 ENDIF -
trunk/NEMO/OPA_SRC/TRA/trasbc.F90
r503 r592 110 110 zse3t = 1. / fse3t(ji,jj,1) 111 111 #endif 112 zta = ro0cpr * ( qt(ji,jj) - qsr(ji,jj) ) * zse3t ! temperature : heat flux 113 zsa = emps(ji,jj) * zsrau * sn(ji,jj,1) * zse3t ! salinity : concent./dilut. effect 114 ta(ji,jj,1) = ta(ji,jj,1) + zta ! add the trend to the general tracer trend 112 IF( lk_vvl) THEN 113 zta = ro0cpr * ( qt(ji,jj) - qsr(ji,jj) ) * zse3t & ! temperature : heat flux 114 & - emp(ji,jj) * zsrau * tn(ji,jj,1) * zse3t ! & cooling/heating effet of EMP flux 115 zsa = 0.e0 ! No salinity concent./dilut. effect 116 ELSE 117 zta = ro0cpr * ( qt(ji,jj) - qsr(ji,jj) ) * zse3t ! temperature : heat flux 118 zsa = emps(ji,jj) * zsrau * sn(ji,jj,1) * zse3t ! salinity : concent./dilut. effect 119 ENDIF 120 ta(ji,jj,1) = ta(ji,jj,1) + zta ! add the trend to the general tracer trend 115 121 sa(ji,jj,1) = sa(ji,jj,1) + zsa 116 122 END DO -
trunk/NEMO/OPA_SRC/TRA/trazdf.F90
r583 r592 24 24 USE in_out_manager ! I/O manager 25 25 USE prtctl ! Print control 26 27 USE phycst 28 USE dynspg_oce 29 USE ocesbc ! ocean surface boundary condition 30 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 USE domvvl ! variable volume 26 32 27 33 IMPLICIT NONE … … 73 79 ztrds(:,:,:) = sa(:,:,:) 74 80 ENDIF 81 82 IF( lk_vvl ) CALL dom_vvl_ssh( kt ) ! compute ssha field 75 83 76 84 SELECT CASE ( nzdf ) ! compute lateral mixing trend and add it to the general trend -
trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r503 r592 31 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 32 USE prtctl ! Print control 33 USE domvvl ! variable volume 33 34 34 35 IMPLICIT NONE … … 97 98 !! * Local declarations 98 99 INTEGER :: ji, jj, jk ! dummy loop indices 99 REAL(wp) :: zavi, zrhs ! temporary scalars 100 REAL(wp) :: zavi, zrhs, znvvl, & ! temporary scalars 101 ze3tb, ze3tn, ze3ta, zvsfvvl ! variable vertical scale factors 100 102 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 101 103 zwi, zwt, zavsi ! workspace arrays … … 119 121 zavsi(:,:,jpk) = 0.e0 ; zavsi( : ,:,1) = 0.e0 120 122 123 ! I.1 Variable volume : to take into account vertical variable vertical scale factors 124 ! ------------------- 125 IF( lk_vvl ) THEN ; znvvl = 1. 126 ELSE ; znvvl = 0.e0 127 ENDIF 128 121 129 ! II. Vertical trend associated with the vertical physics 122 130 ! ======================================================= … … 145 153 END DO 146 154 155 #else 156 ! No isopycnal diffusion 157 zwt(:,:,:) = avt(:,:,:) 158 # if defined key_zdfddm 159 zavsi(:,:,:) = avs(:,:,:) 160 # endif 161 162 #endif 163 147 164 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 148 165 DO jk = 1, jpkm1 149 166 DO jj = 2, jpjm1 150 167 DO ji = fs_2, fs_jpim1 ! vector opt. 151 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 152 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 153 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 154 END DO 155 END DO 156 END DO 157 #else 158 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 159 DO jk = 1, jpkm1 160 DO jj = 2, jpjm1 161 DO ji = fs_2, fs_jpim1 ! vector opt. 162 zwi(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 163 zws(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 164 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 165 END DO 166 END DO 167 END DO 168 #endif 168 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 169 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 170 ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point 171 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 172 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 173 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 174 END DO 175 END DO 176 END DO 169 177 170 178 ! Surface boudary conditions 171 179 DO jj = 2, jpjm1 172 180 DO ji = fs_2, fs_jpim1 ! vector opt. 181 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) ) 182 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 173 183 zwi(ji,jj,1) = 0.e0 174 zwd(ji,jj,1) = 1.- zws(ji,jj,1)184 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 175 185 END DO 176 186 END DO … … 214 224 DO jj = 2, jpjm1 215 225 DO ji = fs_2, fs_jpim1 216 ta(ji,jj,1) = tb(ji,jj,1) + p2dt(1) * ta(ji,jj,1) 217 END DO 218 END DO 219 DO jk = 2, jpkm1 220 DO jj = 2, jpjm1 221 DO ji = fs_2, fs_jpim1 222 zrhs = tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ! zrhs=right hand side 226 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) ) 227 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 228 ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,1) 229 ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 230 END DO 231 END DO 232 DO jk = 2, jpkm1 233 DO jj = 2, jpjm1 234 DO ji = fs_2, fs_jpim1 235 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) 236 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl 237 ze3tn = ( 1. - znvvl ) + znvvl*fse3t (ji,jj,jk) 238 zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk) ! zrhs=right hand side 223 239 ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1) 224 240 END DO … … 249 265 250 266 ! Diagonal, inferior, superior (including the bottom boundary condition via avs masked) 251 # if defined key_ldfslp252 267 DO jk = 1, jpkm1 253 268 DO jj = 2, jpjm1 254 269 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 256 zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 257 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 258 END DO 259 END DO 260 END DO 261 # else 262 DO jk = 1, jpkm1 263 DO jj = 2, jpjm1 264 DO ji = fs_2, fs_jpim1 ! vector opt. 265 zwi(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 266 zws(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 267 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 268 END DO 269 END DO 270 END DO 271 # endif 270 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + ssha(ji,jj) * mut(ji,jj,jk) ) 271 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 272 ze3tn = ( 1. - znvvl )*fse3t(ji,jj,jk) + znvvl ! now scale factor at T-point 273 zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 274 zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 275 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 276 END DO 277 END DO 278 END DO 272 279 273 280 ! Surface boudary conditions 274 281 DO jj = 2, jpjm1 275 282 DO ji = fs_2, fs_jpim1 ! vector opt. 283 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + ssha(ji,jj) * mut(ji,jj,1) ) 284 ze3ta = ( 1. - znvvl ) + znvvl*zvsfvvl ! after scale factor at T-point 276 285 zwi(ji,jj,1) = 0.e0 277 zwd(ji,jj,1) = 1.- zws(ji,jj,1)286 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 278 287 END DO 279 288 END DO … … 316 325 DO jj = 2, jpjm1 317 326 DO ji = fs_2, fs_jpim1 318 sa(ji,jj,1) = sb(ji,jj,1) + p2dt(1) * sa(ji,jj,1) 319 END DO 320 END DO 321 DO jk = 2, jpkm1 322 DO jj = 2, jpjm1 323 DO ji = fs_2, fs_jpim1 324 zrhs = sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ! zrhs=right hand side 327 zvsfvvl = fsve3t(ji,jj,1) * ( 1 + sshb(ji,jj) * mut(ji,jj,1) ) 328 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point 329 ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,1) ! now scale factor at T-point 330 sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) 331 END DO 332 END DO 333 DO jk = 2, jpkm1 334 DO jj = 2, jpjm1 335 DO ji = fs_2, fs_jpim1 336 zvsfvvl = fsve3t(ji,jj,jk) * ( 1 + sshb(ji,jj) * mut(ji,jj,jk) ) 337 ze3tb = ( 1. - znvvl ) + znvvl*zvsfvvl ! before scale factor at T-point 338 ze3tn = ( 1. - znvvl ) + znvvl*fse3t(ji,jj,jk) ! now scale factor at T-point 339 zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk) ! zrhs=right hand side 325 340 sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) 326 341 END DO
Note: See TracChangeset
for help on using the changeset viewer.
