Changeset 1566 for trunk/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
- Timestamp:
- 2009-07-31T16:34:08+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynadv_ubs.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r1528 r1566 5 5 !! trend using a 3rd order upstream biased scheme 6 6 !!====================================================================== 7 !! History : 9.0 ! 06-08 (R. Benshila, L. Debreu) Original code 7 !! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code 8 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 9 !!---------------------------------------------------------------------- 9 10 … … 27 28 REAL(wp), PARAMETER :: gamma2 = 1._wp/8._wp ! =0 2nd order ; =1/8 4th order centred 28 29 29 !! * Routine accessibility 30 PUBLIC dyn_adv_ubs ! routine called by step.F90 30 PUBLIC dyn_adv_ubs ! routine called by step.F90 31 31 32 32 !! * Substitutions … … 34 34 # include "vectopt_loop_substitute.h90" 35 35 !!---------------------------------------------------------------------- 36 !! OPA 9.0 , LODYC-IPSL (2006)36 !! NEMO/OPA 3.2 , LODYC-IPSL (2009) 37 37 !! $Id$ 38 38 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) … … 46 46 !! 47 47 !! ** Purpose : Compute the now momentum advection trend in flux form 48 !! and the general trend of the momentum equation.48 !! and the general trend of the momentum equation. 49 49 !! 50 50 !! ** Method : The scheme is the one implemeted in ROMS. It depends … … 64 64 !! gamma1=1/4 and gamma2=1/8. 65 65 !! 66 !! ** Action : - update (ua,va)with the 3D advective momentum trends66 !! ** Action : - (ua,va) updated with the 3D advective momentum trends 67 67 !! 68 68 !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 69 69 !!---------------------------------------------------------------------- 70 USE oce, ONLY: zfu => ta , &! use ta as 3D workspace71 zfv => sa! use sa as 3D workspace72 70 USE oce, ONLY: zfu => ta ! use ta as 3D workspace 71 USE oce, ONLY: zfv => sa ! use sa as 3D workspace 72 !! 73 73 INTEGER, INTENT(in) :: kt ! ocean time-step index 74 75 INTEGER :: ji, jj, jk 76 REAL(wp) :: z ua, zva, zbu, zbv! temporary scalars77 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! temporary scalars74 !! 75 INTEGER :: ji, jj, jk ! dummy loop indices 76 REAL(wp) :: zbu, zbv ! temporary scalars 77 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! temporary scalars 78 78 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu_t, zfu_f ! temporary workspace 79 79 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfv_t, zfv_f ! " " … … 92 92 zfu_f(:,:,:) = 0.e0 93 93 zfv_f(:,:,:) = 0.e0 94 94 ! 95 95 zlu_uu(:,:,:,:) = 0.e0 96 96 zlv_vv(:,:,:,:) = 0.e0 … … 103 103 ENDIF 104 104 105 ! ! =============== 106 DO jk = 1, jpkm1 ! Horizontal slab 107 ! ! =============== 108 ! Laplacian 109 ! --------- 105 ! ! =========================== ! 106 DO jk = 1, jpkm1 ! Laplacian of the velocity ! 107 ! ! =========================== ! 108 ! ! horizontal volume fluxes 110 109 zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 111 110 zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 112 113 DO jj = 2, jpjm1 111 ! 112 DO jj = 2, jpjm1 ! laplacian 114 113 DO ji = fs_2, fs_jpim1 ! vector opt. 115 114 zlu_uu(ji,jj,jk,1) = ( ub (ji+1,jj,jk)-2.*ub (ji,jj,jk)+ub (ji-1,jj,jk) ) * umask(ji,jj,jk) … … 124 123 END DO 125 124 END DO 126 ! ! =============== 127 END DO ! End of slab 128 ! ! =============== 129 130 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) 131 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) 132 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) 133 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) 134 135 CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 136 CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 137 CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 138 CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 139 140 ! I. Horizontal advection 141 ! ----------------------- 142 ! ! =============== 143 DO jk = 1, jpkm1 ! Horizontal slab 144 ! ! =============== 145 ! horizontal volume fluxes 125 END DO 126 !!gm BUG !!! just below this should be +1 in all the communications 127 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', -1.) 128 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', -1.) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', -1.) 129 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', -1.) 130 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', -1.) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', -1.) 131 132 !!gm corrected: 133 CALL lbc_lnk( zlu_uu(:,:,:,1), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,1), 'U', 1. ) 134 CALL lbc_lnk( zlu_uu(:,:,:,2), 'U', 1. ) ; CALL lbc_lnk( zlu_uv(:,:,:,2), 'U', 1. ) 135 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 136 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 137 !!gm end 138 139 ! ! ====================== ! 140 ! ! Horizontal advection ! 141 DO jk = 1, jpkm1 ! ====================== ! 142 ! ! horizontal volume fluxes 146 143 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 147 144 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 148 149 ! horizontal momentum fluxes at T- and F-point 150 DO jj = 1, jpjm1 145 ! 146 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point 151 147 DO ji = 1, fs_jpim1 ! vector opt. 152 148 zui = ( un(ji,jj,jk) + un(ji+1,jj ,jk) ) 153 149 zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 154 150 ! 155 151 IF (zui > 0) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) 156 152 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) … … 159 155 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1) 160 156 ENDIF 161 157 ! 162 158 zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) & 163 159 & - gamma2 * ( zlu_uu(ji,jj,jk,2) + zlu_uu(ji+1,jj ,jk,2) ) ) & … … 166 162 & - gamma2 * ( zlv_vv(ji,jj,jk,2) + zlv_vv(ji ,jj+1,jk,2) ) ) & 167 163 & * ( zvj - gamma1 * zl_v) 168 164 ! 169 165 zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) ) 170 166 zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) ) … … 175 171 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1) 176 172 ENDIF 177 173 ! 178 174 zfv_f(ji ,jj ,jk) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,jk,2) + zlv_vu(ji+1,jj ,jk,2) ) ) & 179 175 & * ( un(ji,jj,jk) + un(ji ,jj+1,jk) - gamma1 * zl_u ) … … 182 178 END DO 183 179 END DO 184 185 ! divergence of horizontal momentum fluxes 186 DO jj = 2, jpjm1 180 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 187 181 DO ji = fs_2, fs_jpim1 ! vector opt. 188 182 zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) 189 183 zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) 190 ! horizontal advective trends 191 zua = - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 192 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 193 zva = - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 194 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 195 ! add it to the general tracer trends 196 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 197 va(ji,jj,jk) = va(ji,jj,jk) + zva 198 END DO 199 END DO 200 ! ! =============== 201 END DO ! End of slab 202 ! ! =============== 203 204 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 184 ! 185 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 186 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 187 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 188 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 189 END DO 190 END DO 191 END DO 192 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 205 193 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 206 194 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) 207 195 CALL trd_mod( zfu_uw, zfv_vw, jpdyn_trd_had, 'DYN', kt ) 208 ENDIF209 210 ! II. Vertical advection211 ! ----------------------212 213 IF( l_trddyn ) THEN ! Save ua and va trends214 196 zfu_t(:,:,:) = ua(:,:,:) 215 197 zfv_t(:,:,:) = va(:,:,:) 216 198 ENDIF 217 199 218 ! Second order centered tracer flux at w-point 219 DO jk = 1, jpkm1 220 ! Vertical volume fluxes 200 ! ! ==================== ! 201 ! ! Vertical advection ! 202 DO jk = 1, jpkm1 ! ==================== ! 203 ! ! Vertical volume fluxesÊ 221 204 zfw(:,:,jk) = 0.25 * e1t(:,:) * e2t(:,:) * wn(:,:,jk) 222 ! Vertical advective fluxes223 IF( jk == 1 ) THEN 224 zfu_uw(:,:,jpk) = 0.e0 ! Bottomvalue : flux set to zero205 ! 206 IF( jk == 1 ) THEN ! surface/bottom advective fluxes 207 zfu_uw(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero 225 208 zfv_vw(:,:,jpk) = 0.e0 226 ! ! Surface value 227 IF( lk_vvl ) THEN 228 ! variable volume : flux set to zero 209 ! ! Surface value : 210 IF( lk_vvl ) THEN ! variable volume : flux set to zero 229 211 zfu_uw(:,:, 1 ) = 0.e0 230 212 zfv_vw(:,:, 1 ) = 0.e0 231 ELSE 232 ! free surface-constant volume 213 ELSE ! constant volume : advection through the surface 233 214 DO jj = 2, jpjm1 234 215 DO ji = fs_2, fs_jpim1 … … 238 219 END DO 239 220 ENDIF 240 ELSE 241 ! ! interior fluxes 221 ELSE ! interior fluxes 242 222 DO jj = 2, jpjm1 243 223 DO ji = fs_2, fs_jpim1 ! vector opt. … … 248 228 ENDIF 249 229 END DO 250 251 ! momentum flux divergence at u-, v-points added to the general trend 252 DO jk = 1, jpkm1 230 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 253 231 DO jj = 2, jpjm1 254 232 DO ji = fs_2, fs_jpim1 ! vector opt. 255 zua =- ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) &233 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) & 256 234 & / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 257 zva =- ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) &235 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) & 258 236 & / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 259 ! add it to the general tracer trends 260 ua(ji,jj,jk) = ua(ji,jj,jk) + zua 261 va(ji,jj,jk) = va(ji,jj,jk) + zva 262 END DO 263 END DO 264 END DO 265 266 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 237 END DO 238 END DO 239 END DO 240 ! 241 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 267 242 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 268 243 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 269 244 CALL trd_mod( zfu_t, zfv_t, jpdyn_trd_zad, 'DYN', kt ) 270 245 ENDIF 271 272 ! ! Control print 246 ! ! Control print 273 247 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, & 274 248 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 275 249 ! 276 250 END SUBROUTINE dyn_adv_ubs 277 251
Note: See TracChangeset
for help on using the changeset viewer.
