Changeset 10874 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
- Timestamp:
- 2019-04-15T15:57:37+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
r10825 r10874 45 45 CONTAINS 46 46 47 SUBROUTINE dyn_zdf( kt , ktlev1, ktlev2, ktlev3, kt2lev, pu_rhs, pv_rhs)47 SUBROUTINE dyn_zdf( kt ) 48 48 !!---------------------------------------------------------------------- 49 49 !! *** ROUTINE dyn_zdf *** … … 54 54 !! 55 55 !! ** Method : - Leap-Frog time stepping on all trends but the vertical mixing 56 !! pu_rhs = uu(:,:,:,ktlev1) + 2*dt * pu_rhsvector form or linear free surf.57 !! pu_rhs = ( e3u_b*uu(:,:,:,ktlev1) + 2*dt * e3u_n*pu_rhs ) / e3u(:,:,:,ktlev3)otherwise56 !! ua = ub + 2*dt * ua vector form or linear free surf. 57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_a otherwise 58 58 !! - update the after velocity with the implicit vertical mixing. 59 59 !! This requires to solver the following system: 60 !! pu_rhs = pu_rhs + 1/e3u(:,:,:,ktlev3)dk+1[ mi(avm) / e3uw_a dk[ua] ]60 !! ua = ua + 1/e3u_a 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 : ( pu_rhs,pv_rhs) after velocity65 !! ** Action : (ua,va) after velocity 66 66 !!--------------------------------------------------------------------- 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 67 INTEGER, INTENT(in) :: kt ! ocean time-step index 71 68 ! 72 69 INTEGER :: ji, jj, jk ! dummy loop indices … … 99 96 ! 100 97 ! !* explicit top/bottom drag case 101 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, u u(:,:,:,ktlev1), vv(:,:,:,ktlev1), pu_rhs, pv_rhs ) ! add top/bottom friction trend to (pu_rhs,pv_rhs)98 IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, ub, vb, ua, va ) ! add top/bottom friction trend to (ua,va) 102 99 ! 103 100 ! 104 101 IF( l_trddyn ) THEN !* temporary save of ta and sa trends 105 102 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 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)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) 111 108 ! 112 109 ! ! time stepping except vertical diffusion 113 110 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity 114 111 DO jk = 1, jpkm1 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)112 ua(:,:,jk) = ( ub(:,:,jk) + r2dt * ua(:,:,jk) ) * umask(:,:,jk) 113 va(:,:,jk) = ( vb(:,:,jk) + r2dt * va(:,:,jk) ) * vmask(:,:,jk) 117 114 END DO 118 115 ELSE ! applied on thickness weighted velocity 119 116 DO jk = 1, jpkm1 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)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) 124 121 END DO 125 122 ENDIF … … 128 125 ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does 129 126 ! not lead to the effective stress seen over the whole barotropic loop. 130 ! G. Madec : in linear free surface, e3u (:,:,:,ktlev3) = e3u(:,:,:,ktlev2) = e3u_0, so systematic use of e3u(:,:,:,ktlev3)127 ! G. Madec : in linear free surface, e3u_a = e3u_n = e3u_0, so systematic use of e3u_a 131 128 IF( ln_drgimp .AND. ln_dynspg_ts ) THEN 132 129 DO jk = 1, jpkm1 ! remove barotropic velocities 133 pu_rhs(:,:,jk) = ( pu_rhs(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk)134 pv_rhs(:,:,jk) = ( pv_rhs(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk)130 ua(:,:,jk) = ( ua(:,:,jk) - ua_b(:,:) ) * umask(:,:,jk) 131 va(:,:,jk) = ( va(:,:,jk) - va_b(:,:) ) * vmask(:,:,jk) 135 132 END DO 136 133 DO jj = 2, jpjm1 ! Add bottom/top stress due to barotropic component only … … 138 135 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 139 136 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 140 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) / ze3ua143 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) / ze3va137 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) / ze3ua 140 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / ze3va 144 141 END DO 145 142 END DO … … 149 146 iku = miku(ji,jj) ! top ocean level at u- and v-points 150 147 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 151 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) / ze3ua154 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) / ze3va148 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) / ze3ua 151 va(ji,jj,ikv) = va(ji,jj,ikv) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / ze3va 155 152 END DO 156 153 END DO … … 168 165 DO jj = 2, jpjm1 169 166 DO ji = fs_2, fs_jpim1 ! vector opt. 170 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point167 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 171 168 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 172 & / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )169 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 173 170 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 174 & / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)171 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 175 172 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 176 173 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 185 182 DO jj = 2, jpjm1 186 183 DO ji = fs_2, fs_jpim1 ! vector opt. 187 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point188 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)184 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 185 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) 190 187 zWui = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 191 188 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 200 197 DO ji = fs_2, fs_jpim1 ! vector opt. 201 198 zwi(ji,jj,1) = 0._wp 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)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) 204 201 zWus = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) 205 202 zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp ) … … 213 210 DO jj = 2, jpjm1 214 211 DO ji = fs_2, fs_jpim1 ! vector opt. 215 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point212 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 216 213 zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) & 217 & / ( ze3ua * e3uw (ji,jj,jk ,kt2lev) ) * wumask(ji,jj,jk )214 & / ( ze3ua * e3uw_n(ji,jj,jk ) ) * wumask(ji,jj,jk ) 218 215 zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) & 219 & / ( ze3ua * e3uw (ji,jj,jk+1,kt2lev) ) * wumask(ji,jj,jk+1)216 & / ( ze3ua * e3uw_n(ji,jj,jk+1) ) * wumask(ji,jj,jk+1) 220 217 zwi(ji,jj,jk) = zzwi 221 218 zws(ji,jj,jk) = zzws … … 228 225 DO jj = 2, jpjm1 229 226 DO ji = fs_2, fs_jpim1 ! vector opt. 230 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,jk,ktlev2) + r_vvl * e3u(ji,jj,jk,ktlev3) ! after scale factor at U-point231 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)227 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,jk) + r_vvl * e3u_a(ji,jj,jk) ! after scale factor at U-point 228 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) 233 230 zwi(ji,jj,jk) = zzwi 234 231 zws(ji,jj,jk) = zzws … … 257 254 DO ji = 2, jpim1 258 255 iku = mbku(ji,jj) ! ocean bottom level at u- and v-points 259 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point256 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 260 257 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua 261 258 END DO … … 266 263 !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed 267 264 iku = miku(ji,jj) ! ocean top level at u- and v-points 268 ze3ua = ( 1._wp - r_vvl ) * e3u (ji,jj,iku,ktlev2) + r_vvl * e3u(ji,jj,iku,ktlev3) ! after scale factor at T-point265 ze3ua = ( 1._wp - r_vvl ) * e3u_n(ji,jj,iku) + r_vvl * e3u_a(ji,jj,iku) ! after scale factor at T-point 269 266 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua 270 267 END DO … … 285 282 ! m is decomposed in the product of an upper and a lower triangular matrix 286 283 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 287 ! The solution (the after velocity) is in pu_rhs284 ! The solution (the after velocity) is in ua 288 285 !----------------------------------------------------------------------- 289 286 ! … … 298 295 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 299 296 DO ji = fs_2, fs_jpim1 ! vector opt. 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) ) &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) ) & 302 299 & / ( ze3ua * rau0 ) * umask(ji,jj,1) 303 300 END DO … … 306 303 DO jj = 2, jpjm1 307 304 DO ji = fs_2, fs_jpim1 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)305 ua(ji,jj,jk) = ua(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 309 306 END DO 310 307 END DO … … 313 310 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==! 314 311 DO ji = fs_2, fs_jpim1 ! vector opt. 315 pu_rhs(ji,jj,jpkm1) = pu_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)312 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 316 313 END DO 317 314 END DO … … 319 316 DO jj = 2, jpjm1 320 317 DO ji = fs_2, fs_jpim1 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)318 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 322 319 END DO 323 320 END DO … … 334 331 DO jj = 2, jpjm1 335 332 DO ji = fs_2, fs_jpim1 ! vector opt. 336 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point333 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 337 334 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 338 & / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )335 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 339 336 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 340 & / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)337 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 341 338 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 342 339 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 351 348 DO jj = 2, jpjm1 352 349 DO ji = fs_2, fs_jpim1 ! vector opt. 353 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point354 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)350 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 351 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) 356 353 zWvi = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 357 354 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 366 363 DO ji = fs_2, fs_jpim1 ! vector opt. 367 364 zwi(ji,jj,1) = 0._wp 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)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) 370 367 zWvs = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) 371 368 zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp ) … … 379 376 DO jj = 2, jpjm1 380 377 DO ji = fs_2, fs_jpim1 ! vector opt. 381 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point378 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 382 379 zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) & 383 & / ( ze3va * e3vw (ji,jj,jk ,kt2lev) ) * wvmask(ji,jj,jk )380 & / ( ze3va * e3vw_n(ji,jj,jk ) ) * wvmask(ji,jj,jk ) 384 381 zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) & 385 & / ( ze3va * e3vw (ji,jj,jk+1,kt2lev) ) * wvmask(ji,jj,jk+1)382 & / ( ze3va * e3vw_n(ji,jj,jk+1) ) * wvmask(ji,jj,jk+1) 386 383 zwi(ji,jj,jk) = zzwi 387 384 zws(ji,jj,jk) = zzws … … 394 391 DO jj = 2, jpjm1 395 392 DO ji = fs_2, fs_jpim1 ! vector opt. 396 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,jk,ktlev2) + r_vvl * e3v(ji,jj,jk,ktlev3) ! after scale factor at V-point397 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)393 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,jk) + r_vvl * e3v_a(ji,jj,jk) ! after scale factor at V-point 394 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) 399 396 zwi(ji,jj,jk) = zzwi 400 397 zws(ji,jj,jk) = zzws … … 422 419 DO ji = 2, jpim1 423 420 ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 424 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point421 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 425 422 zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va 426 423 END DO … … 430 427 DO ji = 2, jpim1 431 428 ikv = mikv(ji,jj) ! (first wet ocean u- and v-points) 432 ze3va = ( 1._wp - r_vvl ) * e3v (ji,jj,ikv,ktlev2) + r_vvl * e3v(ji,jj,ikv,ktlev3) ! after scale factor at T-point429 ze3va = ( 1._wp - r_vvl ) * e3v_n(ji,jj,ikv) + r_vvl * e3v_a(ji,jj,ikv) ! after scale factor at T-point 433 430 zwd(ji,jj,iku) = zwd(ji,jj,iku) - r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3va 434 431 END DO … … 449 446 ! m is decomposed in the product of an upper and lower triangular matrix 450 447 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 451 ! The solution (after velocity) is in 2d array pv_rhs448 ! The solution (after velocity) is in 2d array va 452 449 !----------------------------------------------------------------------- 453 450 ! … … 462 459 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==! 463 460 DO ji = fs_2, fs_jpim1 ! vector opt. 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) ) &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) ) & 466 463 & / ( ze3va * rau0 ) * vmask(ji,jj,1) 467 464 END DO … … 470 467 DO jj = 2, jpjm1 471 468 DO ji = fs_2, fs_jpim1 ! vector opt. 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)469 va(ji,jj,jk) = va(ji,jj,jk) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 473 470 END DO 474 471 END DO … … 477 474 DO jj = 2, jpjm1 !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==! 478 475 DO ji = fs_2, fs_jpim1 ! vector opt. 479 pv_rhs(ji,jj,jpkm1) = pv_rhs(ji,jj,jpkm1) / zwd(ji,jj,jpkm1)476 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 480 477 END DO 481 478 END DO … … 483 480 DO jj = 2, jpjm1 484 481 DO ji = fs_2, fs_jpim1 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)482 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 486 483 END DO 487 484 END DO … … 489 486 ! 490 487 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 491 ztrdu(:,:,:) = ( pu_rhs(:,:,:) - uu(:,:,:,ktlev1) ) / r2dt - ztrdu(:,:,:)492 ztrdv(:,:,:) = ( pv_rhs(:,:,:) - vv(:,:,:,ktlev1) ) / r2dt - ztrdv(:,:,:)488 ztrdu(:,:,:) = ( ua(:,:,:) - ub(:,:,:) ) / r2dt - ztrdu(:,:,:) 489 ztrdv(:,:,:) = ( va(:,:,:) - vb(:,:,:) ) / r2dt - ztrdv(:,:,:) 493 490 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 494 491 DEALLOCATE( ztrdu, ztrdv )
Note: See TracChangeset
for help on using the changeset viewer.
