Changeset 10825 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA
- 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/TRA/trazdf.F90
r10425 r10825 44 44 CONTAINS 45 45 46 SUBROUTINE tra_zdf( kt )46 SUBROUTINE tra_zdf( kt, ktlev1, ktlev2, ktlev3, kt2lev, pts_rhs ) 47 47 !!---------------------------------------------------------------------- 48 48 !! *** ROUTINE tra_zdf *** … … 51 51 !!--------------------------------------------------------------------- 52 52 INTEGER, INTENT(in) :: kt ! ocean time-step index 53 INTEGER, INTENT(in) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms 54 INTEGER, INTENT(in) :: kt2lev ! time level index for 2-time-level source terms 55 REAL(wp), INTENT( inout), DIMENSION(jpi,jpj,jpk,jpts) :: pts_rhs ! temperature and salinity trends 53 56 ! 54 57 INTEGER :: jk ! Dummy loop indices … … 70 73 IF( l_trdtra ) THEN !* Save ta and sa trends 71 74 ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) ) 72 ztrdt(:,:,:) = tsa(:,:,:,jp_tem)73 ztrds(:,:,:) = tsa(:,:,:,jp_sal)75 ztrdt(:,:,:) = pts_rhs(:,:,:,jp_tem) 76 ztrds(:,:,:) = pts_rhs(:,:,:,jp_sal) 74 77 ENDIF 75 78 ! 76 79 ! !* compute lateral mixing trend and add it to the general trend 77 CALL tra_zdf_imp( kt, nit000, 'TRA', r2dt, tsb, tsa, jpts )80 CALL tra_zdf_imp( kt, nit000, ktlev1, ktlev2, ktlev3, kt2lev, 'TRA', r2dt, ts(:,:,:,:,ktlev1), pts_rhs, jpts ) 78 81 79 82 !!gm WHY here ! and I don't like that ! … … 81 84 ! JMM avoid negative salinities near river outlet ! Ugly fix 82 85 ! JMM : restore negative salinities to small salinities: 83 WHERE( tsa(:,:,:,jp_sal) < 0._wp ) tsa(:,:,:,jp_sal) = 0.1_wp86 WHERE( pts_rhs(:,:,:,jp_sal) < 0._wp ) pts_rhs(:,:,:,jp_sal) = 0.1_wp 84 87 !!gm 85 88 86 89 IF( l_trdtra ) THEN ! save the vertical diffusive trends for further diagnostics 87 90 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)91 ztrdt(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_tem)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_tem,ktlev1)*e3t(:,:,jk,ktlev1) ) & 92 & / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrdt(:,:,jk) 93 ztrds(:,:,jk) = ( ( pts_rhs(:,:,jk,jp_sal)*e3t(:,:,jk,ktlev3) - ts(:,:,jk,jp_sal,ktlev1)*e3t(:,:,jk,ktlev1) ) & 94 & / (e3t(:,:,jk,ktlev2)*r2dt) ) - ztrds(:,:,jk) 92 95 END DO 93 96 !!gm this should be moved in trdtra.F90 and done on all trends … … 107 110 108 111 109 SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, ptb, pta, kjpt )112 SUBROUTINE tra_zdf_imp( kt, kit000, ktlev1, ktlev2, ktlev3, kt2lev, cdtype, p2dt, pt, pt_rhs, kjpt ) 110 113 !!---------------------------------------------------------------------- 111 114 !! *** ROUTINE tra_zdf_imp *** … … 125 128 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 126 129 !! 127 !! ** Action : - pt abecomes the after tracer130 !! ** Action : - pt_rhs becomes the after tracer 128 131 !!--------------------------------------------------------------------- 129 132 INTEGER , INTENT(in ) :: kt ! ocean time-step index 130 133 INTEGER , INTENT(in ) :: kit000 ! first time step index 134 INTEGER , INTENT(in ) :: ktlev1, ktlev2, ktlev3 ! time level indices for 3-time-level source terms 135 INTEGER , INTENT(in ) :: kt2lev ! time level index for 2-time-level source terms 131 136 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 132 137 INTEGER , INTENT(in ) :: kjpt ! number of tracers 133 138 REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step 134 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt b! before and now tracer fields135 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt a! in: tracer trend ; out: after tracer field139 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt ! before and now tracer fields 140 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! in: tracer trend ; out: after tracer field 136 141 ! 137 142 INTEGER :: ji, jj, jk, jn ! dummy loop indices … … 181 186 DO jj = 2, jpjm1 182 187 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 &188 zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,kt2lev) 189 zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev) 190 zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zzwi - zzws & 186 191 & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) ) 187 192 zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp ) … … 194 199 DO jj = 2, jpjm1 195 200 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)201 zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,kt2lev) 202 zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,kt2lev) 203 zwd(ji,jj,jk) = e3t(ji,jj,jk,ktlev3) - zwi(ji,jj,jk) - zws(ji,jj,jk) 199 204 END DO 200 205 END DO … … 216 221 ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal 217 222 ! and "superior" (above diagonal) components of the tridiagonal system. 218 ! The solution will be in the 4d array pt a.223 ! The solution will be in the 4d array pt_rhs. 219 224 ! 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 then225 ! En route to the solution pt_rhs is used a to evaluate the rhs and then 221 226 ! used as a work space array: its value is modified. 222 227 ! … … 238 243 DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1 239 244 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)245 pt_rhs(ji,jj,1,jn) = e3t(ji,jj,1,ktlev1) * pt(ji,jj,1,jn) + p2dt * e3t(ji,jj,1,ktlev2) * pt_rhs(ji,jj,1,jn) 241 246 END DO 242 247 END DO … … 244 249 DO jj = 2, jpjm1 245 250 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)251 zrhs = e3t(ji,jj,jk,ktlev1) * pt(ji,jj,jk,jn) + p2dt * e3t(ji,jj,jk,ktlev2) * pt_rhs(ji,jj,jk,jn) ! zrhs=right hand side 252 pt_rhs(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt_rhs(ji,jj,jk-1,jn) 248 253 END DO 249 254 END DO … … 252 257 DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 253 258 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)259 pt_rhs(ji,jj,jpkm1,jn) = pt_rhs(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 255 260 END DO 256 261 END DO … … 258 263 DO jj = 2, jpjm1 259 264 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) ) &265 pt_rhs(ji,jj,jk,jn) = ( pt_rhs(ji,jj,jk,jn) - zws(ji,jj,jk) * pt_rhs(ji,jj,jk+1,jn) ) & 261 266 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 262 267 END DO
Note: See TracChangeset
for help on using the changeset viewer.
