- Timestamp:
- 2019-11-22T15:29:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/dynzdf.F90
r11281 r11949 45 45 CONTAINS 46 46 47 SUBROUTINE dyn_zdf( kt )47 SUBROUTINE dyn_zdf( kt, Kbb, Kmm, Krhs, puu, pvv, Kaa ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE dyn_zdf *** … … 54 54 !! 55 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! u a = ub + 2*dt * uavector form or linear free surf.57 !! u a = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_aotherwise56 !! u(after) = u(before) + 2*dt * u(rhs) vector form or linear free surf. 57 !! u(after) = ( e3u_b*u(before) + 2*dt * e3u_n*u(rhs) ) / e3u(after) otherwise 58 58 !! - update the after velocity with the implicit vertical mixing. 59 59 !! This requires to solver the following system: 60 !! u a = ua + 1/e3u_a dk+1[ mi(avm) / e3uw_adk[ua] ]60 !! u(after) = u(after) + 1/e3u(after) dk+1[ mi(avm) / e3uw(after) dk[ua] ] 61 61 !! with the following surface/top/bottom boundary condition: 62 62 !! surface: wind stress input (averaged over kt-1/2 & kt+1/2) 63 63 !! top & bottom : top stress (iceshelf-ocean) & bottom stress (cf zdfdrg.F90) 64 64 !! 65 !! ** Action : ( ua,va) after velocity65 !! ** Action : (puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) after velocity 66 66 !!--------------------------------------------------------------------- 67 INTEGER, INTENT(in) :: kt ! ocean time-step index 67 INTEGER , INTENT( in ) :: kt ! ocean time-step index 68 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 69 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 68 70 ! 69 71 INTEGER :: ji, jj, jk ! dummy loop indices … … 96 98 ! 97 99 ! !* explicit top/bottom drag case 98 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, ub, vb, ua, va ) ! add top/bottom friction trend to (ua,va)100 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa)) 99 101 ! 100 102 ! 101 103 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 102 104 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 103 ztrdu(:,:,:) = ua(:,:,:)104 ztrdv(:,:,:) = va(:,:,:)105 ENDIF 106 ! 107 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in ua,va)105 ztrdu(:,:,:) = puu(:,:,:,Krhs) 106 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 107 ENDIF 108 ! 109 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa)) 108 110 ! 109 111 ! ! time stepping except vertical diffusion 110 112 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 111 113 DO jk = 1, jpkm1 112 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)113 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)114 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kbb) + r2dt * puu(:,:,jk,Krhs) ) * umask(:,:,jk) 115 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kbb) + r2dt * pvv(:,:,jk,Krhs) ) * vmask(:,:,jk) 114 116 END DO 115 117 ELSE ! applied on thickness weighted velocity 116 118 DO jk = 1, jpkm1 117 ua(:,:,jk) = ( e3u_b(:,:,jk) * ub(:,:,jk) &118 & + r2dt * e3u _n(:,:,jk) * ua(:,:,jk) ) / e3u_a(:,:,jk) * umask(:,:,jk)119 va(:,:,jk) = ( e3v_b(:,:,jk) * vb(:,:,jk) &120 & + r2dt * e3v _n(:,:,jk) * va(:,:,jk) ) / e3v_a(:,:,jk) * vmask(:,:,jk)119 puu(:,:,jk,Kaa) = ( e3u(:,:,jk,Kbb) * puu(:,:,jk,Kbb) & 120 & + r2dt * e3u(:,:,jk,Kmm) * puu(:,:,jk,Krhs) ) / e3u(:,:,jk,Kaa) * umask(:,:,jk) 121 pvv(:,:,jk,Kaa) = ( e3v(:,:,jk,Kbb) * pvv(:,:,jk,Kbb) & 122 & + r2dt * e3v(:,:,jk,Kmm) * pvv(:,:,jk,Krhs) ) / e3v(:,:,jk,Kaa) * vmask(:,:,jk) 121 123 END DO 122 124 ENDIF … … 125 127 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 126 128 ! not lead to the effective stress seen over the whole barotropic loop. 127 ! G. Madec : in linear free surface, e3u _a = e3u_n = e3u_0, so systematic use of e3u_a129 ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa) 128 130 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 129 131 DO jk = 1, jpkm1 ! remove barotropic velocities 130 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)131 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)132 puu(:,:,jk,Kaa) = ( puu(:,:,jk,Kaa) - uu_b(:,:,Kaa) ) * umask(:,:,jk) 133 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - vv_b(:,:,Kaa) ) * vmask(:,:,jk) 132 134 END DO 133 135 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only … … 135 137 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 136 138 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 137 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)138 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)139 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua140 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va139 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 140 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 141 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 142 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 141 143 END DO 142 144 END DO … … 146 148 iku = miku(ji,jj) ! top ocean level at u- and v-points 147 149 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 148 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku)149 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv)150 ua(ji,jj,iku) = ua(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua151 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va150 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) 151 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) 152 puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua 153 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va 152 154 END DO 153 155 END DO … … 165 167 DO jj = 2, jpjm1 166 168 DO ji = fs_2, fs_jpim1 ! vector opt. 167 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point169 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 168 170 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 169 & / ( ze3ua * e3uw _n(ji,jj,jk) ) * wumask(ji,jj,jk )171 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 170 172 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 171 & / ( ze3ua * e3uw _n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)173 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 172 174 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 173 175 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua … … 182 184 DO jj = 2, jpjm1 183 185 DO ji = fs_2, fs_jpim1 ! vector opt. 184 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point185 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw _n(ji,jj,jk) ) * wumask(ji,jj,jk )186 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw _n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)186 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 187 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 188 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 187 189 zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua 188 190 zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua … … 197 199 DO ji = fs_2, fs_jpim1 ! vector opt. 198 200 zwi(ji,jj,1) = 0._wp 199 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)200 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw _n(ji,jj,2) ) * wumask(ji,jj,2)201 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 202 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2) 201 203 zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua 202 204 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 210 212 DO jj = 2, jpjm1 211 213 DO ji = fs_2, fs_jpim1 ! vector opt. 212 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point214 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 213 215 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 214 & / ( ze3ua * e3uw _n(ji,jj,jk) ) * wumask(ji,jj,jk )216 & / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 215 217 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 216 & / ( ze3ua * e3uw _n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)218 & / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 217 219 zwi(ji,jj,jk) = zzwi 218 220 zws(ji,jj,jk) = zzws … … 225 227 DO jj = 2, jpjm1 226 228 DO ji = fs_2, fs_jpim1 ! vector opt. 227 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point228 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw _n(ji,jj,jk) ) * wumask(ji,jj,jk )229 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw _n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1)229 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point 230 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk ) 231 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1) 230 232 zwi(ji,jj,jk) = zzwi 231 233 zws(ji,jj,jk) = zzws … … 254 256 DO ji = 2, jpim1 255 257 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 256 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point258 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 257 259 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 258 260 END DO … … 263 265 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 264 266 iku = miku(ji,jj) ! ocean top level at u- and v-points 265 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point267 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point 266 268 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 267 269 END DO … … 282 284 ! m is decomposed in the product of an upper and a lower triangular matrix 283 285 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 284 ! The solution (the after velocity) is in ua286 ! The solution (the after velocity) is in puu(:,:,:,Kaa) 285 287 !----------------------------------------------------------------------- 286 288 ! … … 295 297 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 296 298 DO ji = fs_2, fs_jpim1 ! vector opt. 297 ze3ua = ( 1._wp - r_vvl ) * e3u _n(ji,jj,1) + r_vvl * e3u_a(ji,jj,1)298 ua(ji,jj,1) = ua(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &299 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) + r_vvl * e3u(ji,jj,1,Kaa) 300 puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 299 301 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 300 302 END DO … … 303 305 DO jj = 2, jpjm1 304 306 DO ji = fs_2, fs_jpim1 305 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1)307 puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa) 306 308 END DO 307 309 END DO … … 310 312 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 311 313 DO ji = fs_2, fs_jpim1 ! vector opt. 312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)314 puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 313 315 END DO 314 316 END DO … … 316 318 DO jj = 2, jpjm1 317 319 DO ji = fs_2, fs_jpim1 318 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk)320 puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 319 321 END DO 320 322 END DO … … 331 333 DO jj = 2, jpjm1 332 334 DO ji = fs_2, fs_jpim1 ! vector opt. 333 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point335 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 334 336 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 335 & / ( ze3va * e3vw _n(ji,jj,jk) ) * wvmask(ji,jj,jk )337 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 336 338 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 337 & / ( ze3va * e3vw _n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)339 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 338 340 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 339 341 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va … … 348 350 DO jj = 2, jpjm1 349 351 DO ji = fs_2, fs_jpim1 ! vector opt. 350 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point351 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw _n(ji,jj,jk) ) * wvmask(ji,jj,jk )352 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw _n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)352 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 353 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 354 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 353 355 zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va 354 356 zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va … … 363 365 DO ji = fs_2, fs_jpim1 ! vector opt. 364 366 zwi(ji,jj,1) = 0._wp 365 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)366 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw _n(ji,jj,2) ) * wvmask(ji,jj,2)367 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 368 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2) 367 369 zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va 368 370 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 376 378 DO jj = 2, jpjm1 377 379 DO ji = fs_2, fs_jpim1 ! vector opt. 378 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point380 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 379 381 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 380 & / ( ze3va * e3vw _n(ji,jj,jk) ) * wvmask(ji,jj,jk )382 & / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 381 383 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 382 & / ( ze3va * e3vw _n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)384 & / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 383 385 zwi(ji,jj,jk) = zzwi 384 386 zws(ji,jj,jk) = zzws … … 391 393 DO jj = 2, jpjm1 392 394 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point394 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw _n(ji,jj,jk) ) * wvmask(ji,jj,jk )395 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw _n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1)395 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point 396 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk ) 397 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1) 396 398 zwi(ji,jj,jk) = zzwi 397 399 zws(ji,jj,jk) = zzws … … 419 421 DO ji = 2, jpim1 420 422 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 421 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point423 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 422 424 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 423 425 END DO … … 427 429 DO ji = 2, jpim1 428 430 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 429 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point431 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point 430 432 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 431 433 END DO … … 446 448 ! m is decomposed in the product of an upper and lower triangular matrix 447 449 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 448 ! The solution (after velocity) is in 2d array va450 ! The solution (after velocity) is in 2d array pvv(:,:,:,Kaa) 449 451 !----------------------------------------------------------------------- 450 452 ! … … 459 461 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 460 462 DO ji = fs_2, fs_jpim1 ! vector opt. 461 ze3va = ( 1._wp - r_vvl ) * e3v _n(ji,jj,1) + r_vvl * e3v_a(ji,jj,1)462 va(ji,jj,1) = va(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &463 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) + r_vvl * e3v(ji,jj,1,Kaa) 464 pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 463 465 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 464 466 END DO … … 467 469 DO jj = 2, jpjm1 468 470 DO ji = fs_2, fs_jpim1 ! vector opt. 469 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1)471 pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa) 470 472 END DO 471 473 END DO … … 474 476 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 475 477 DO ji = fs_2, fs_jpim1 ! vector opt. 476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)478 pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1) 477 479 END DO 478 480 END DO … … 480 482 DO jj = 2, jpjm1 481 483 DO ji = fs_2, fs_jpim1 482 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk)484 pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk) 483 485 END DO 484 486 END DO … … 486 488 ! 487 489 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 488 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)489 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)490 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt )490 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:) 491 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:) 492 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm ) 491 493 DEALLOCATE( ztrdu, ztrdv ) 492 494 ENDIF 493 495 ! ! print mean trends (used for debugging) 494 IF(ln_ctl) CALL prt_ctl( tab3d_1= ua, clinfo1=' zdf - Ua: ', mask1=umask, &495 & tab3d_2= va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )496 IF(ln_ctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Kaa), clinfo1=' zdf - Ua: ', mask1=umask, & 497 & tab3d_2=pvv(:,:,:,Kaa), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 496 498 ! 497 499 IF( ln_timing ) CALL timing_stop('dyn_zdf')
Note: See TracChangeset
for help on using the changeset viewer.