- Timestamp:
- 2014-08-05T12:24:00+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r3294 r4736 27 27 PRIVATE 28 28 29 PUBLIC dyn_zad ! routine called by step.F90 29 PUBLIC dyn_zad ! routine called by dynadv.F90 30 PUBLIC dyn_zad_zts ! routine called by dynadv.F90 30 31 31 32 !! * Substitutions … … 83 84 DO jj = 2, jpj ! vertical fluxes 84 85 DO ji = fs_2, jpi ! vector opt. 85 zww(ji,jj) = 0.25 * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk)86 zww(ji,jj) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 86 87 END DO 87 88 END DO … … 95 96 DO jj = 2, jpjm1 ! Surface and bottom values set to zero 96 97 DO ji = fs_2, fs_jpim1 ! vector opt. 97 zwuw(ji,jj, 1 ) = 0. e098 zwvw(ji,jj, 1 ) = 0. e099 zwuw(ji,jj,jpk) = 0. e0100 zwvw(ji,jj,jpk) = 0. e098 zwuw(ji,jj, 1 ) = 0._wp 99 zwvw(ji,jj, 1 ) = 0._wp 100 zwuw(ji,jj,jpk) = 0._wp 101 zwvw(ji,jj,jpk) = 0._wp 101 102 END DO 102 103 END DO … … 132 133 END SUBROUTINE dyn_zad 133 134 135 SUBROUTINE dyn_zad_zts ( kt ) 136 !!---------------------------------------------------------------------- 137 !! *** ROUTINE dynzad_zts *** 138 !! 139 !! ** Purpose : Compute the now vertical momentum advection trend and 140 !! add it to the general trend of momentum equation. This version 141 !! uses sub-timesteps for improved numerical stability with small 142 !! vertical grid sizes. This is especially relevant when using 143 !! embedded ice with thin surface boxes. 144 !! 145 !! ** Method : The now vertical advection of momentum is given by: 146 !! w dz(u) = ua + 1/(e1u*e2u*e3u) mk+1[ mi(e1t*e2t*wn) dk(un) ] 147 !! w dz(v) = va + 1/(e1v*e2v*e3v) mk+1[ mj(e1t*e2t*wn) dk(vn) ] 148 !! Add this trend to the general trend (ua,va): 149 !! (ua,va) = (ua,va) + w dz(u,v) 150 !! 151 !! ** Action : - Update (ua,va) with the vert. momentum adv. trends 152 !! - Save the trends in (ztrdu,ztrdv) ('key_trddyn') 153 !!---------------------------------------------------------------------- 154 INTEGER, INTENT(in) :: kt ! ocean time-step inedx 155 ! 156 INTEGER :: ji, jj, jk, jl ! dummy loop indices 157 INTEGER :: jnzts = 5 ! number of sub-timesteps for vertical advection 158 INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps 159 REAL(wp) :: zua, zva ! temporary scalars 160 REAL(wp) :: zr_rdt ! temporary scalar 161 REAL(wp) :: z2dtzts ! length of Euler forward sub-timestep for vertical advection 162 REAL(wp) :: zts ! length of sub-timestep for vertical advection 163 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwuw , zwvw, zww 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv 165 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zus , zvs 166 !!---------------------------------------------------------------------- 167 ! 168 IF( nn_timing == 1 ) CALL timing_start('dyn_zad_zts') 169 ! 170 CALL wrk_alloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 171 CALL wrk_alloc( jpi,jpj,jpk,3, zus, zvs ) 172 ! 173 IF( kt == nit000 ) THEN 174 IF(lwp)WRITE(numout,*) 175 IF(lwp)WRITE(numout,*) 'dyn_zad_zts : arakawa advection scheme with sub-timesteps' 176 ENDIF 177 178 IF( l_trddyn ) THEN ! Save ua and va trends 179 CALL wrk_alloc( jpi, jpj, jpk, ztrdu, ztrdv ) 180 ztrdu(:,:,:) = ua(:,:,:) 181 ztrdv(:,:,:) = va(:,:,:) 182 ENDIF 183 184 IF( neuler == 0 .AND. kt == nit000 ) THEN 185 z2dtzts = rdt / REAL( jnzts, wp ) ! = rdt (restart with Euler time stepping) 186 ELSE 187 z2dtzts = 2._wp * rdt / REAL( jnzts, wp ) ! = 2 rdt (leapfrog) 188 ENDIF 189 190 DO jk = 2, jpkm1 ! Calculate and store vertical fluxes 191 DO jj = 2, jpj 192 DO ji = fs_2, jpi ! vector opt. 193 zww(ji,jj,jk) = 0.25_wp * e1t(ji,jj) * e2t(ji,jj) * wn(ji,jj,jk) 194 END DO 195 END DO 196 END DO 197 198 DO jj = 2, jpjm1 ! Surface and bottom advective fluxes set to zero 199 DO ji = fs_2, fs_jpim1 ! vector opt. 200 zwuw(ji,jj, 1 ) = 0._wp 201 zwvw(ji,jj, 1 ) = 0._wp 202 zwuw(ji,jj,jpk) = 0._wp 203 zwvw(ji,jj,jpk) = 0._wp 204 END DO 205 END DO 206 207 ! Start with before values and use sub timestepping to reach after values 208 209 zus(:,:,:,1) = ub(:,:,:) 210 zvs(:,:,:,1) = vb(:,:,:) 211 212 DO jl = 1, jnzts ! Start of sub timestepping loop 213 214 IF( jl == 1 ) THEN ! Euler forward to kick things off 215 jtb = 1 ; jtn = 1 ; jta = 2 216 zts = z2dtzts 217 ELSEIF( jl == 2 ) THEN ! First leapfrog step 218 jtb = 1 ; jtn = 2 ; jta = 3 219 zts = 2._wp * z2dtzts 220 ELSE ! Shuffle pointers for subsequent leapfrog steps 221 jtb = MOD(jtb,3) + 1 222 jtn = MOD(jtn,3) + 1 223 jta = MOD(jta,3) + 1 224 ENDIF 225 226 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 227 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 228 DO ji = fs_2, fs_jpim1 ! vector opt. 229 zwuw(ji,jj,jk) = ( zww(ji+1,jj ,jk) + zww(ji,jj,jk) ) * ( zus(ji,jj,jk-1,jtn)-zus(ji,jj,jk,jtn) ) 230 zwvw(ji,jj,jk) = ( zww(ji ,jj+1,jk) + zww(ji,jj,jk) ) * ( zvs(ji,jj,jk-1,jtn)-zvs(ji,jj,jk,jtn) ) 231 END DO 232 END DO 233 END DO 234 DO jk = 1, jpkm1 ! Vertical momentum advection at u- and v-points 235 DO jj = 2, jpjm1 236 DO ji = fs_2, fs_jpim1 ! vector opt. 237 ! ! vertical momentum advective trends 238 zua = - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) ) 239 zva = - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) ) 240 zus(ji,jj,jk,jta) = zus(ji,jj,jk,jtb) + zua * zts 241 zvs(ji,jj,jk,jta) = zvs(ji,jj,jk,jtb) + zva * zts 242 END DO 243 END DO 244 END DO 245 246 END DO ! End of sub timestepping loop 247 248 zr_rdt = 1._wp / ( REAL( jnzts, wp ) * z2dtzts ) 249 DO jk = 1, jpkm1 ! Recover trends over the outer timestep 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 ! vector opt. 252 ! ! vertical momentum advective trends 253 ! ! add the trends to the general momentum trends 254 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zus(ji,jj,jk,jta) - ub(ji,jj,jk)) * zr_rdt 255 va(ji,jj,jk) = va(ji,jj,jk) + ( zvs(ji,jj,jk,jta) - vb(ji,jj,jk)) * zr_rdt 256 END DO 257 END DO 258 END DO 259 260 IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic 261 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 262 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 263 CALL trd_mod(ztrdu, ztrdv, jpdyn_trd_zad, 'DYN', kt) 264 CALL wrk_dealloc( jpi, jpj, jpk, ztrdu, ztrdv ) 265 ENDIF 266 ! ! Control print 267 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' zad - Ua: ', mask1=umask, & 268 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 269 ! 270 CALL wrk_dealloc( jpi,jpj,jpk, zwuw , zwvw, zww ) 271 CALL wrk_dealloc( jpi,jpj,jpk,3, zus, zvs ) 272 ! 273 IF( nn_timing == 1 ) CALL timing_stop('dyn_zad_zts') 274 ! 275 END SUBROUTINE dyn_zad_zts 276 134 277 !!====================================================================== 135 278 END MODULE dynzad
Note: See TracChangeset
for help on using the changeset viewer.
