Changeset 10825 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
- Timestamp:
- 2019-04-02T16:40:54+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
r10364 r10825 45 45 CONTAINS 46 46 47 SUBROUTINE dyn_zdf( kt )47 SUBROUTINE dyn_zdf( kt, ktlev1, ktlev2, ktlev3, kt2lev, pu_rhs, pv_rhs ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE dyn_zdf *** … … 54 54 !! 55 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! ua = ub + 2*dt * uavector form or linear free surf.57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_aotherwise56 !! pu_rhs = uu(:,:,:,ktlev1) + 2*dt * pu_rhs vector form or linear free surf. 57 !! pu_rhs = ( e3u_b*uu(:,:,:,ktlev1) + 2*dt * e3u_n*pu_rhs ) / e3u(:,:,:,ktlev3) otherwise 58 58 !! - update the after velocity with the implicit vertical mixing. 59 59 !! This requires to solver the following system: 60 !! ua = ua + 1/e3u_adk+1[ mi(avm) / e3uw_a dk[ua] ]60 !! pu_rhs = pu_rhs + 1/e3u(:,:,:,ktlev3) dk+1[ mi(avm) / e3uw_a 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 : (pu_rhs,pv_rhs) 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) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms 69 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms 70 REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pu_rhs, pv_rhs ! momentum trends -> momentum after fields 68 71 ! 69 72 INTEGER :: ji, jj, jk ! dummy loop indices … … 96 99 ! 97 100 ! !* explicit top/bottom drag case 98 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, u b, vb, ua, va ) ! add top/bottom friction trend to (ua,va)101 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, uu(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs ) ! add top/bottom friction trend to (pu_rhs,pv_rhs) 99 102 ! 100 103 ! 101 104 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 102 105 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)106 ztrdu(:,:,:) = pu_rhs(:,:,:) 107 ztrdv(:,:,:) = pv_rhs(:,:,:) 108 ENDIF 109 ! 110 ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in pu_rhs,pv_rhs) 108 111 ! 109 112 ! ! time stepping except vertical diffusion 110 113 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 111 114 DO jk = 1, jpkm1 112 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk)113 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk)115 pu_rhs(:,:,jk) = ( uu(:,:,jk,ktlev1) + r2dt * pu_rhs(:,:,jk) ) * umask(:,:,jk) 116 pv_rhs(:,:,jk) = ( vv(:,:,jk,ktlev1) + r2dt * pv_rhs(:,:,jk) ) * vmask(:,:,jk) 114 117 END DO 115 118 ELSE ! applied on thickness weighted velocity 116 119 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)120 pu_rhs(:,:,jk) = ( e3u(:,:,jk,ktlev1) * uu(:,:,jk,ktlev1) & 121 & + r2dt * e3u(:,:,jk,ktlev2) * pu_rhs(:,:,jk) ) / e3u(:,:,jk,ktlev3) * umask(:,:,jk) 122 pv_rhs(:,:,jk) = ( e3v(:,:,jk,ktlev1) * vv(:,:,jk,ktlev1) & 123 & + r2dt * e3v(:,:,jk,ktlev2) * pv_rhs(:,:,jk) ) / e3v(:,:,jk,ktlev3) * vmask(:,:,jk) 121 124 END DO 122 125 ENDIF … … 125 128 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 126 129 ! 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_a130 ! G. Madec : in linear free surface, e3u(:,:,:,ktlev3) = e3u(:,:,:,ktlev2) = e3u_0, so systematic use of e3u(:,:,:,ktlev3) 128 131 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 129 132 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)133 pu_rhs(:,:,jk) = ( pu_rhs(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 134 pv_rhs(:,:,jk) = ( pv_rhs(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 132 135 END DO 133 136 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only … … 135 138 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 136 139 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) / ze3va140 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) 141 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) 142 pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * ua_b(ji,jj) / ze3ua 143 pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 141 144 END DO 142 145 END DO … … 146 149 iku = miku(ji,jj) ! top ocean level at u- and v-points 147 150 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) / ze3va151 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) 152 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) 153 pu_rhs(ji,jj,iku) = pu_rhs(ji,jj,iku) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * ua_b(ji,jj) / ze3ua 154 pv_rhs(ji,jj,ikv) = pv_rhs(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 152 155 END DO 153 156 END DO … … 165 168 DO jj = 2, jpjm1 166 169 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-point170 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point 168 171 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 )172 & / ( ze3ua * e3uw(ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk ) 170 173 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)174 & / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 172 175 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 173 176 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 182 185 DO jj = 2, jpjm1 183 186 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)187 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point 188 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk ) 189 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 187 190 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 188 191 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 197 200 DO ji = fs_2, fs_jpim1 ! vector opt. 198 201 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)202 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3) 203 zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) / ( ze3ua * e3uw(ji,jj,2,kt2lev) ) * wumask(ji,jj,2) 201 204 zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) 202 205 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 210 213 DO jj = 2, jpjm1 211 214 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-point215 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point 213 216 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 )217 & / ( ze3ua * e3uw(ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk ) 215 218 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)219 & / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 217 220 zwi(ji,jj,jk) = zzwi 218 221 zws(ji,jj,jk) = zzws … … 225 228 DO jj = 2, jpjm1 226 229 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)230 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point 231 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) / ( ze3ua * e3uw(ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk ) 232 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) / ( ze3ua * e3uw(ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1) 230 233 zwi(ji,jj,jk) = zzwi 231 234 zws(ji,jj,jk) = zzws … … 254 257 DO ji = 2, jpim1 255 258 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-point259 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point 257 260 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 258 261 END DO … … 263 266 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 264 267 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-point268 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point 266 269 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 267 270 END DO … … 282 285 ! m is decomposed in the product of an upper and a lower triangular matrix 283 286 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 284 ! The solution (the after velocity) is in ua287 ! The solution (the after velocity) is in pu_rhs 285 288 !----------------------------------------------------------------------- 286 289 ! … … 295 298 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 296 299 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) ) &300 ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,ktlev2) + r_vvl * e3u(ji,jj,1,ktlev3) 301 pu_rhs(ji,jj,1) = pu_rhs(ji,jj,1) + r2dt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 299 302 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 300 303 END DO … … 303 306 DO jj = 2, jpjm1 304 307 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)308 pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pu_rhs(ji,jj,jk-1) 306 309 END DO 307 310 END DO … … 310 313 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 311 314 DO ji = fs_2, fs_jpim1 ! vector opt. 312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)315 pu_rhs(ji,jj,jpkm1) = pu_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 313 316 END DO 314 317 END DO … … 316 319 DO jj = 2, jpjm1 317 320 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)321 pu_rhs(ji,jj,jk) = ( pu_rhs(ji,jj,jk) - zws(ji,jj,jk) * pu_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk) 319 322 END DO 320 323 END DO … … 331 334 DO jj = 2, jpjm1 332 335 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-point336 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point 334 337 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 )338 & / ( ze3va * e3vw(ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk ) 336 339 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)340 & / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 338 341 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 339 342 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 348 351 DO jj = 2, jpjm1 349 352 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)353 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point 354 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk ) 355 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 353 356 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 354 357 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 363 366 DO ji = fs_2, fs_jpim1 ! vector opt. 364 367 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)368 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3) 369 zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) / ( ze3va * e3vw(ji,jj,2,kt2lev) ) * wvmask(ji,jj,2) 367 370 zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) 368 371 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 376 379 DO jj = 2, jpjm1 377 380 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-point381 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point 379 382 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 )383 & / ( ze3va * e3vw(ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk ) 381 384 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)385 & / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 383 386 zwi(ji,jj,jk) = zzwi 384 387 zws(ji,jj,jk) = zzws … … 391 394 DO jj = 2, jpjm1 392 395 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)396 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point 397 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) / ( ze3va * e3vw(ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk ) 398 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) / ( ze3va * e3vw(ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1) 396 399 zwi(ji,jj,jk) = zzwi 397 400 zws(ji,jj,jk) = zzws … … 419 422 DO ji = 2, jpim1 420 423 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-point424 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point 422 425 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 423 426 END DO … … 427 430 DO ji = 2, jpim1 428 431 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-point432 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point 430 433 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 431 434 END DO … … 446 449 ! m is decomposed in the product of an upper and lower triangular matrix 447 450 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 448 ! The solution (after velocity) is in 2d array va451 ! The solution (after velocity) is in 2d array pv_rhs 449 452 !----------------------------------------------------------------------- 450 453 ! … … 459 462 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 460 463 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) ) &464 ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,ktlev2) + r_vvl * e3v(ji,jj,1,ktlev3) 465 pv_rhs(ji,jj,1) = pv_rhs(ji,jj,1) + r2dt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 463 466 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 464 467 END DO … … 467 470 DO jj = 2, jpjm1 468 471 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)472 pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pv_rhs(ji,jj,jk-1) 470 473 END DO 471 474 END DO … … 474 477 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 475 478 DO ji = fs_2, fs_jpim1 ! vector opt. 476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)479 pv_rhs(ji,jj,jpkm1) = pv_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 477 480 END DO 478 481 END DO … … 480 483 DO jj = 2, jpjm1 481 484 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)485 pv_rhs(ji,jj,jk) = ( pv_rhs(ji,jj,jk) - zws(ji,jj,jk) * pv_rhs(ji,jj,jk+1) ) / zwd(ji,jj,jk) 483 486 END DO 484 487 END DO … … 486 489 ! 487 490 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 488 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:)489 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:)491 ztrdu(:,:,:) = ( pu_rhs(:,:,:) - uu(:,:,:,ktlev1) ) / r2dt - ztrdu(:,:,:) 492 ztrdv(:,:,:) = ( pv_rhs(:,:,:) - vv(:,:,:,ktlev1) ) / r2dt - ztrdv(:,:,:) 490 493 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 491 494 DEALLOCATE( ztrdu, ztrdv )
Note: See TracChangeset
for help on using the changeset viewer.