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