Changeset 10884
- Timestamp:
- 2019-04-18T17:00:30+02:00 (6 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src
- Files:
-
- 6 edited
-
OCE/DYN/dynzdf.F90 (modified) (31 diffs)
-
OCE/TRA/traadv_ubs.F90 (modified) (1 diff)
-
OCE/TRA/trazdf.F90 (modified) (13 diffs)
-
OCE/step.F90 (modified) (2 diffs)
-
TOP/TRP/trctrp.F90 (modified) (2 diffs)
-
TOP/TRP/trczdf.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynzdf.F90
r10874 r10884 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 !! ua = ub + 2*dt * uavector form or linear free surf.57 !! ua = ( e3u_b*ub + 2*dt * e3u_n*ua ) / e3u_aotherwise56 !! puu(:,:,:,Kaa) = puu(:,:,:,Kbb) + 2*dt * puu(:,:,:,Krhs) vector form or linear free surf. 57 !! puu(:,:,:,Kaa) = ( e3u_b*puu(:,:,:,Kbb) + 2*dt * e3u_n*puu(:,:,:,Krhs) ) / e3u(:,:,:,Kaa) 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_a dk+1[ mi(avm) / e3uw_adk[ua] ]60 !! puu(:,:,:,Kaa) = puu(:,:,:,Kaa) + 1/e3u(:,:,:,Kaa) dk+1[ mi(avm) / e3uw(:,:,:,Kaa) 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, 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) - ua_b(:,:) ) * umask(:,:,jk) 133 pvv(:,:,jk,Kaa) = ( pvv(:,:,jk,Kaa) - va_b(:,:) ) * 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) ) * ua_b(ji,jj) / ze3ua 142 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * va_b(ji,jj) / 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) ) * ua_b(ji,jj) / ze3ua 153 pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + r2dt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * va_b(ji,jj) / 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 = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 173 175 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 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 = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) 188 190 zWus = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) … … 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 = 0.5_wp * ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) 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 = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 339 341 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 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 = 0.5_wp * ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * wvmask(ji,jj,jk ) 354 356 zWvs = 0.5_wp * ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * wvmask(ji,jj,jk+1) … … 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 = 0.5_wp * ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) 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 ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / r2dt - ztrdu(:,:,:) 491 ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / r2dt - ztrdv(:,:,:) 490 492 CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt ) 491 493 DEALLOCATE( ztrdu, ztrdv ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90
r10880 r10884 92 92 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 93 93 REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu_mm, pv_mm, pww ! 3 ocean transport components 94 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! activetracers and RHS of tracer equation94 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 95 95 ! 96 96 INTEGER :: ji, jj, jk, jn ! dummy loop indices -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trazdf.F90
r10874 r10884 44 44 CONTAINS 45 45 46 SUBROUTINE tra_zdf( kt )46 SUBROUTINE tra_zdf( kt, Kbb, Kmm, Krhs, pts, Kaa ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_zdf *** … … 50 50 !! ** Purpose : compute the vertical ocean tracer physics. 51 51 !!--------------------------------------------------------------------- 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 52 INTEGER , INTENT(in) :: kt ! ocean time-step index 53 INTEGER , INTENT(in) :: Kbb, Kmm, Krhs, Kaa ! time level indices 54 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 53 55 ! 54 56 INTEGER :: jk ! Dummy loop indices … … 70 72 IF( l_trdtra ) THEN !* Save ta and sa trends 71 73 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)74 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Kaa) 75 ztrds(:,:,:) = pts(:,:,:,jp_sal,Kaa) 74 76 ENDIF 75 77 ! 76 78 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )79 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, Kbb, Kmm, pts, Kaa, jpts ) 78 80 79 81 !!gm WHY here ! and I don't like that ! … … 81 83 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 84 ! JMM : restore negative salinities to small salinities: 83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp85 WHERE( pts(:,:,:,jp_sal,Kaa) < 0._wp ) pts(:,:,:,jp_sal,Kaa) = 0.1_wp 84 86 !!gm 85 87 86 88 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 89 DO jk = 1, jpkm1 88 ztrdt(:,:,jk) = ( ( tsa(:,:,jk,jp_tem)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_tem)*e3t_b(:,:,jk) ) &89 & / (e3t _n(:,:,jk)*r2dt) ) - ztrdt(:,:,jk)90 ztrds(:,:,jk) = ( ( tsa(:,:,jk,jp_sal)*e3t_a(:,:,jk) - tsb(:,:,jk,jp_sal)*e3t_b(:,:,jk) ) &91 & / (e3t _n(:,:,jk)*r2dt) ) - ztrds(:,:,jk)90 ztrdt(:,:,jk) = ( ( pts(:,:,jk,jp_tem,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_tem,Kbb)*e3t(:,:,jk,Kbb) ) & 91 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrdt(:,:,jk) 92 ztrds(:,:,jk) = ( ( pts(:,:,jk,jp_sal,Kaa)*e3t(:,:,jk,Kaa) - pts(:,:,jk,jp_sal,Kbb)*e3t(:,:,jk,Kbb) ) & 93 & / (e3t(:,:,jk,Kmm)*r2dt) ) - ztrds(:,:,jk) 92 94 END DO 93 95 !!gm this should be moved in trdtra.F90 and done on all trends … … 107 109 108 110 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )111 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, pt, Kaa, kjpt ) 110 112 !!---------------------------------------------------------------------- 111 113 !! *** ROUTINE tra_zdf_imp *** … … 125 127 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 128 !! 127 !! ** Action : - pt abecomes the after tracer128 !!--------------------------------------------------------------------- 129 INTEGER , INTENT(in ) :: kt ! ocean time-step index130 INTEGER , INTENT(in ) :: kit000 ! first time step index131 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)132 INTEGER , INTENT(in ) :: kjpt ! number of tracers133 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step134 REAL(wp) , DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields135 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt ), INTENT(inout) :: pta ! in: tracer trend ; out: after tracer field129 !! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer 130 !!--------------------------------------------------------------------- 131 INTEGER , INTENT(in ) :: kt ! ocean time-step index 132 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 133 INTEGER , INTENT(in ) :: kit000 ! first time step index 134 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 135 INTEGER , INTENT(in ) :: kjpt ! number of tracers 136 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 137 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation 136 138 ! 137 139 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 181 183 DO jj = 2, jpjm1 182 184 DO ji = fs_2, fs_jpim1 ! vector opt. (ensure same order of calculation as below if wi=0.) 183 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w _n(ji,jj,jk)184 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w _n(ji,jj,jk+1)185 zwd(ji,jj,jk) = e3t _a(ji,jj,jk) - zzwi - zzws &185 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm) 186 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 187 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws & 186 188 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 187 189 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) … … 194 196 DO jj = 2, jpjm1 195 197 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w _n(ji,jj,jk)197 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w _n(ji,jj,jk+1)198 zwd(ji,jj,jk) = e3t _a(ji,jj,jk) - zwi(ji,jj,jk) - zws(ji,jj,jk)198 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm) 199 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm) 200 zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk) 199 201 END DO 200 202 END DO … … 218 220 ! The solution will be in the 4d array pta. 219 221 ! The 3d array zwt is used as a work space array. 220 ! En route to the solution pt ais used a to evaluate the rhs and then222 ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then 221 223 ! used as a work space array: its value is modified. 222 224 ! … … 238 240 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 239 241 DO ji = fs_2, fs_jpim1 240 pt a(ji,jj,1,jn) = e3t_b(ji,jj,1) * ptb(ji,jj,1,jn) + p2dt * e3t_n(ji,jj,1) * pta(ji,jj,1,jn)242 pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Kaa) 241 243 END DO 242 244 END DO … … 244 246 DO jj = 2, jpjm1 245 247 DO ji = fs_2, fs_jpim1 246 zrhs = e3t _b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * e3t_n(ji,jj,jk) * pta(ji,jj,jk,jn) ! zrhs=right hand side247 pt a(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn)248 zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Kaa) ! zrhs=right hand side 249 pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa) 248 250 END DO 249 251 END DO … … 252 254 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 253 255 DO ji = fs_2, fs_jpim1 254 pt a(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)256 pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 255 257 END DO 256 258 END DO … … 258 260 DO jj = 2, jpjm1 259 261 DO ji = fs_2, fs_jpim1 260 pt a(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) &262 pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) & 261 263 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 262 264 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
r10883 r10884 208 208 ENDIF 209 209 210 CALL dyn_zdf ( kstp ) ! vertical diffusion210 CALL dyn_zdf( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion ==> after 211 211 212 212 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> … … 261 261 IF( ln_diaptr ) CALL dia_ptr ! Poleward adv/ldf TRansports diagnostics 262 262 !!gm 263 CALL tra_zdf ( kstp ) ! vertical mixing and after tracer fields263 CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vert. mixing & after tracer ==> after 264 264 IF( ln_zdfnpc ) CALL tra_npc ( kstp ) ! update after fields by non-penetrative convection 265 265 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
r10880 r10884 77 77 IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_trc ! tracers sponge 78 78 #endif 79 CALL trc_zdf ( kt ) ! vertical mixing and after tracer fields79 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa ) ! vert. mixing & after tracer ==> after 80 80 CALL trc_nxt ( kt ) ! tracer fields at next time step 81 81 IF( ln_trcrad ) CALL trc_rad ( kt ) ! Correct artificial negative concentrations … … 86 86 CALL trc_sbc( kt ) ! surface boundary condition 87 87 IF( ln_trcdmp ) CALL trc_dmp( kt ) ! internal damping trends 88 CALL trc_zdf( kt ) ! vertical mixing and after tracer fields88 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa ) ! vert. mixing & after tracer ==> after 89 89 CALL trc_nxt( kt ) ! tracer fields at next time step 90 90 IF( ln_trcrad ) CALL trc_rad( kt ) ! Correct artificial negative concentrations -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trczdf.F90
r10068 r10884 36 36 CONTAINS 37 37 38 SUBROUTINE trc_zdf( kt )38 SUBROUTINE trc_zdf( kt, Kbb, Kmm, Krhs, ptr, Kaa ) 39 39 !!---------------------------------------------------------------------- 40 40 !! *** ROUTINE trc_zdf *** … … 43 43 !! an implicit time-stepping scheme. 44 44 !!--------------------------------------------------------------------- 45 INTEGER, INTENT( in ) :: kt ! ocean time-step index 45 INTEGER , INTENT(in ) :: kt ! ocean time-step index 46 INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 47 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation 46 48 ! 47 49 INTEGER :: jk, jn … … 52 54 IF( ln_timing ) CALL timing_start('trc_zdf') 53 55 ! 54 IF( l_trdtrc ) ztrtrd(:,:,:,:) = tra(:,:,:,:)56 IF( l_trdtrc ) ztrtrd(:,:,:,:) = ptr(:,:,:,:,Krhs) 55 57 ! 56 CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, trb, tra, jptra ) ! implicit scheme58 CALL tra_zdf_imp( kt, nittrc000, 'TRC', r2dttrc, Kbb, Kmm, ptr, Kaa, jptra ) ! implicit scheme 57 59 ! 58 60 IF( l_trdtrc ) THEN ! save the vertical diffusive trends for further diagnostics 59 61 DO jn = 1, jptra 60 62 DO jk = 1, jpkm1 61 ztrtrd(:,:,jk,jn) = ( ( tra(:,:,jk,jn) - trb(:,:,jk,jn) ) / r2dttrc ) - ztrtrd(:,:,jk,jn)63 ztrtrd(:,:,jk,jn) = ( ( ptr(:,:,jk,jn,Kaa) - ptr(:,:,jk,jn,Kbb) ) / r2dttrc ) - ztrtrd(:,:,jk,jn) 62 64 END DO 63 65 CALL trd_tra( kt, 'TRC', jn, jptra_zdf, ztrtrd(:,:,:,jn) )
Note: See TracChangeset
for help on using the changeset viewer.
