- Timestamp:
- 2016-11-28T17:04:10+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_INGV_UKMO_2016/NEMOGCM/NEMO/OPA_SRC/DYN/dynadv_ubs.F90
r5836 r7351 35 35 36 36 !! * Substitutions 37 # include "domzgr_substitute.h90"38 37 # include "vectopt_loop_substitute.h90" 39 38 !!---------------------------------------------------------------------- … … 71 70 !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling. 72 71 !!---------------------------------------------------------------------- 73 INTEGER, INTENT(in) :: kt ! ocean time-step index 74 ! 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 72 INTEGER, INTENT(in) :: kt ! ocean time-step index 73 ! 74 INTEGER :: ji, jj, jk ! dummy loop indices 75 REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars 78 76 REAL(wp), POINTER, DIMENSION(:,:,: ) :: zfu, zfv 79 77 REAL(wp), POINTER, DIMENSION(:,:,: ) :: zfu_t, zfv_t, zfu_f, zfv_f, zfu_uw, zfv_vw, zfw … … 83 81 IF( nn_timing == 1 ) CALL timing_start('dyn_adv_ubs') 84 82 ! 85 CALL wrk_alloc( jpi, jpj, jpk,zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )86 CALL wrk_alloc( jpi, jpj, jpk, jpts,zlu_uu, zlv_vv, zlu_uv, zlv_vu )83 CALL wrk_alloc( jpi,jpj,jpk, zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 84 CALL wrk_alloc( jpi,jpj,jpk,jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu ) 87 85 ! 88 86 IF( kt == nit000 ) THEN … … 101 99 zlu_uv(:,:,:,:) = 0._wp 102 100 zlv_vu(:,:,:,:) = 0._wp 103 104 IF( l_trddyn ) THEN ! Save ua and vatrends101 ! 102 IF( l_trddyn ) THEN ! trends: store the input trends 105 103 zfu_uw(:,:,:) = ua(:,:,:) 106 104 zfv_vw(:,:,:) = va(:,:,:) 107 105 ENDIF 108 109 106 ! ! =========================== ! 110 107 DO jk = 1, jpkm1 ! Laplacian of the velocity ! 111 108 ! ! =========================== ! 112 109 ! ! horizontal volume fluxes 113 zfu(:,:,jk) = e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)114 zfv(:,:,jk) = e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)110 zfu(:,:,jk) = e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 111 zfv(:,:,jk) = e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 115 112 ! 116 113 DO jj = 2, jpjm1 ! laplacian 117 114 DO ji = fs_2, fs_jpim1 ! vector opt. 118 !119 115 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) 120 116 zlv_vv(ji,jj,jk,1) = ( vb (ji ,jj+1,jk) - 2.*vb (ji,jj,jk) + vb (ji ,jj-1,jk) ) * vmask(ji,jj,jk) … … 137 133 CALL lbc_lnk( zlv_vv(:,:,:,1), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,1), 'V', 1. ) 138 134 CALL lbc_lnk( zlv_vv(:,:,:,2), 'V', 1. ) ; CALL lbc_lnk( zlv_vu(:,:,:,2), 'V', 1. ) 139 135 ! 140 136 ! ! ====================== ! 141 137 ! ! Horizontal advection ! 142 138 DO jk = 1, jpkm1 ! ====================== ! 143 139 ! ! horizontal volume fluxes 144 zfu(:,:,jk) = 0.25 * e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk)145 zfv(:,:,jk) = 0.25 * e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk)140 zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u_n(:,:,jk) * un(:,:,jk) 141 zfv(:,:,jk) = 0.25_wp * e1v(:,:) * e3v_n(:,:,jk) * vn(:,:,jk) 146 142 ! 147 143 DO jj = 1, jpjm1 ! horizontal momentum fluxes at T- and F-point … … 150 146 zvj = ( vn(ji,jj,jk) + vn(ji ,jj+1,jk) ) 151 147 ! 152 IF (zui > 0) THEN ; zl_u = zlu_uu(ji ,jj,jk,1)153 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1)154 ENDIF 155 IF (zvj > 0) THEN ; zl_v = zlv_vv(ji,jj ,jk,1)156 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1)148 IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk,1) 149 ELSE ; zl_u = zlu_uu(ji+1,jj,jk,1) 150 ENDIF 151 IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,jk,1) 152 ELSE ; zl_v = zlv_vv(ji,jj+1,jk,1) 157 153 ENDIF 158 154 ! … … 166 162 zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) ) 167 163 zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) ) 168 IF (zfuj > 0) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1)169 ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1)170 ENDIF 171 IF (zfvi > 0) THEN ; zl_u = zlu_uv( ji,jj ,jk,1)172 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1)164 IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk,1) 165 ELSE ; zl_v = zlv_vu( ji+1,jj,jk,1) 166 ENDIF 167 IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk,1) 168 ELSE ; zl_u = zlu_uv( ji,jj+1,jk,1) 173 169 ENDIF 174 170 ! … … 181 177 DO jj = 2, jpjm1 ! divergence of horizontal momentum fluxes 182 178 DO ji = fs_2, fs_jpim1 ! vector opt. 183 zbu = e1e2u(ji,jj) * fse3u(ji,jj,jk) 184 zbv = e1e2v(ji,jj) * fse3v(ji,jj,jk) 185 ! 186 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj ,jk) - zfu_t(ji ,jj ,jk) & 187 & + zfv_f(ji ,jj ,jk) - zfv_f(ji ,jj-1,jk) ) / zbu 188 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji ,jj ,jk) - zfu_f(ji-1,jj ,jk) & 189 & + zfv_t(ji ,jj+1,jk) - zfv_t(ji ,jj ,jk) ) / zbv 190 END DO 191 END DO 192 END DO 193 IF( l_trddyn ) THEN ! save the horizontal advection trend for diagnostic 179 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) & 180 & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 181 va(ji,jj,jk) = va(ji,jj,jk) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) & 182 & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 183 END DO 184 END DO 185 END DO 186 IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic 194 187 zfu_uw(:,:,:) = ua(:,:,:) - zfu_uw(:,:,:) 195 188 zfv_vw(:,:,:) = va(:,:,:) - zfv_vw(:,:,:) … … 198 191 zfv_t(:,:,:) = va(:,:,:) 199 192 ENDIF 200 201 193 ! ! ==================== ! 202 194 ! ! Vertical advection ! 203 DO jk = 1, jpkm1! ==================== !204 ! ! Vertical volume fluxesÊ205 zfw(:,:,jk) = 0.25 * e1e2t(:,:) * wn(:,:,jk)206 !207 IF( jk == 1 ) THEN ! surface/bottom advective fluxes208 zfu_uw( :,:,jpk) = 0.e0 ! Bottom value : flux set to zero209 zfv_vw( :,:,jpk) = 0.e0210 ! ! Surface value :211 IF( lk_vvl ) THEN ! variable volume : flux set to zero212 zfu_uw(:,:, 1 ) = 0.e0213 zfv_vw(:,:, 1 ) = 0.e0214 ELSE ! constant volume : advection through the surface215 DO jj = 2, jpjm1216 DO ji = fs_2, fs_jpim1217 zfu_uw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji+1,jj ,1) ) * un(ji,jj,1)218 zfv_vw(ji,jj, 1 ) = 2.e0 * ( zfw(ji,jj,1) + zfw(ji ,jj+1,1) ) * vn(ji,jj,1)219 END DO220 END DO221 ENDIF222 ELSE ! interior fluxes223 DO jj = 2, jpjm1224 DO ji = fs_2, fs_jpim1 ! vector opt.225 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj ,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) )226 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji ,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) )227 END DO228 END DO229 ENDIF230 END DO231 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence232 DO jj = 2, jpjm1233 DO ji = fs_2, fs_jpim1 ! vector opt.234 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) &235 & / ( e1e2u(ji,jj) * fse3u(ji,jj,jk) )236 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) &237 & / ( e1e2v(ji,jj) * fse3v(ji,jj,jk))238 END DO 239 END DO 240 END DO 241 ! 242 IF( l_trddyn ) THEN 195 ! ! ==================== ! 196 DO jj = 2, jpjm1 ! surface/bottom advective fluxes set to zero 197 DO ji = fs_2, fs_jpim1 198 zfu_uw(ji,jj,jpk) = 0._wp 199 zfv_vw(ji,jj,jpk) = 0._wp 200 zfu_uw(ji,jj, 1 ) = 0._wp 201 zfv_vw(ji,jj, 1 ) = 0._wp 202 END DO 203 END DO 204 IF( ln_linssh ) THEN ! constant volume : advection through the surface 205 DO jj = 2, jpjm1 206 DO ji = fs_2, fs_jpim1 207 zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji+1,jj) * wn(ji+1,jj,1) ) * un(ji,jj,1) 208 zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * wn(ji,jj,1) + e1e2t(ji,jj+1) * wn(ji,jj+1,1) ) * vn(ji,jj,1) 209 END DO 210 END DO 211 ENDIF 212 DO jk = 2, jpkm1 ! interior fluxes 213 DO jj = 2, jpj 214 DO ji = 2, jpi 215 zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * wn(ji,jj,jk) 216 END DO 217 END DO 218 DO jj = 2, jpjm1 219 DO ji = fs_2, fs_jpim1 ! vector opt. 220 zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( un(ji,jj,jk) + un(ji,jj,jk-1) ) 221 zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( vn(ji,jj,jk) + vn(ji,jj,jk-1) ) 222 END DO 223 END DO 224 END DO 225 DO jk = 1, jpkm1 ! divergence of vertical momentum flux divergence 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 ! vector opt. 228 ua(ji,jj,jk) = ua(ji,jj,jk) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk) 229 va(ji,jj,jk) = va(ji,jj,jk) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk) 230 END DO 231 END DO 232 END DO 233 ! 234 IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic 243 235 zfu_t(:,:,:) = ua(:,:,:) - zfu_t(:,:,:) 244 236 zfv_t(:,:,:) = va(:,:,:) - zfv_t(:,:,:) 245 237 CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt ) 246 238 ENDIF 247 ! 239 ! ! Control print 248 240 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ubs2 adv - Ua: ', mask1=umask, & 249 241 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 250 242 ! 251 CALL wrk_dealloc( jpi, jpj, jpk,zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw )252 CALL wrk_dealloc( jpi, jpj, jpk, jpts,zlu_uu, zlv_vv, zlu_uv, zlv_vu )243 CALL wrk_dealloc( jpi,jpj,jpk, zfu_t , zfv_t , zfu_f , zfv_f, zfu_uw, zfv_vw, zfu, zfv, zfw ) 244 CALL wrk_dealloc( jpi,jpj,jpk,jpts, zlu_uu, zlv_vv, zlu_uv, zlv_vu ) 253 245 ! 254 246 IF( nn_timing == 1 ) CALL timing_stop('dyn_adv_ubs')
Note: See TracChangeset
for help on using the changeset viewer.
