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