Changeset 2602
- Timestamp:
- 2011-02-21T16:23:38+01:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r2528 r2602 15 15 !! 3.2 ! 2009-03 (G. Madec) heat and salt content trends 16 16 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 17 !! - ! 2011-02 (A. Coward, C. Ethe, G. Madec) improvment of surface boundary condition 17 18 !!---------------------------------------------------------------------- 18 19 … … 39 40 PUBLIC tra_zdf_imp ! routine called by step.F90 40 41 42 REAL(wp) :: r_vvl ! variable volume indicator, =1 if lk_vvl=T, =0 otherwise 43 41 44 !! * Substitutions 42 45 # include "domzgr_substitute.h90" … … 55 58 !! *** ROUTINE tra_zdf_imp *** 56 59 !! 57 !! ** Purpose : Compute the trend due to the vertical tracer diffusion 58 !! including the vertical component of lateral mixing (only for 2nd 59 !! order operator, for fourth order it is already computed and add 60 !! to the general trend in traldf.F) and add it to the general trend 61 !! of the tracer equations. 62 !! 63 !! ** Method : The vertical component of the lateral diffusive trends 64 !! is provided by a 2nd order operator rotated along neutral or geo- 65 !! potential surfaces to which an eddy induced advection can be 66 !! added. It is computed using before fields (forward in time) and 67 !! isopycnal or geopotential slopes computed in routine ldfslp. 68 !! 69 !! Second part: vertical trend associated with the vertical physics 70 !! =========== (including the vertical flux proportional to dk[t] 71 !! associated with the lateral mixing, through the 72 !! update of avt) 73 !! The vertical diffusion of the tracer t is given by: 74 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 60 !! ** Purpose : Compute the after tracer through a implicit computation 61 !! of the vertical tracer diffusion (including the vertical component 62 !! of lateral mixing (only for 2nd order operator, for fourth order 63 !! it is already computed and add to the general trend in traldf) 64 !! 65 !! ** Method : The vertical diffusion of the tracer t is given by: 66 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 75 67 !! It is computed using a backward time scheme (t=ta). 68 !! If lk_zdfddm=T, use avs for salinity or for passive tracers 76 69 !! Surface and bottom boundary conditions: no diffusive flux on 77 70 !! both tracers (bottom, applied through the masked field avt). 78 !! Add this trend to the general trend ta,sa : 79 !! ta = ta + dz( avt dz(t) ) 80 !! if lk_zdfddm=T, use avs for salinity or for passive tracers 81 !! (sa = sa + dz( avs dz(t) ) 82 !! 83 !! Third part: recover avt resulting from the vertical physics 84 !! ========== alone, for further diagnostics (for example to 85 !! compute the turbocline depth in zdfmxl.F90). 86 !! if lk_zdfddm=T, use avt = zavt 87 !! (avs = zavs if lk_zdfddm=T ) 88 !! 89 !! ** Action : - Update (ta) with before vertical diffusion trend 71 !! If iso-neutral mixing, add to avt the contribution due to lateral mixing. 72 !! 73 !! ** Action : - pta becomes the after tracer 90 74 !!--------------------------------------------------------------------- 91 75 USE oce , ONLY : zwd => ua ! ua used as workspace … … 100 84 !! 101 85 INTEGER :: ji, jj, jk, jn ! dummy loop indices 102 REAL(wp) :: z avi, zrhs, znvvl! local scalars86 REAL(wp) :: zrhs ! local scalars 103 87 REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors 104 88 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt ! workspace arrays … … 109 93 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 110 94 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 111 zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F 95 ! 96 IF( lk_vvl ) THEN ; r_vvl = 1._wp ! Variable volume indicator 97 ELSE ; r_vvl = 0._wp 98 ENDIF 112 99 ENDIF 113 100 ! 114 ! I. Local initialization 115 ! ----------------------- 116 zwd(1,:, : ) = 0._wp ; zwd(jpi,:,:) = 0._wp 117 zws(1,:, : ) = 0._wp ; zws(jpi,:,:) = 0._wp 118 zwi(1,:, : ) = 0._wp ; zwi(jpi,:,:) = 0._wp 119 zwt(1,:, : ) = 0._wp ; zwt(jpi,:,:) = 0._wp 120 zwt(:,:,jpk) = 0._wp ; zwt( : ,:,1) = 0._wp 121 122 ! I.1 Variable volume : to take into account vertical variable vertical scale factors 123 ! ------------------- 124 IF( lk_vvl ) THEN ; znvvl = 1._wp 125 ELSE ; znvvl = 0._wp 126 ENDIF 127 128 ! II. Vertical trend associated with the vertical physics 129 ! ======================================================= 130 ! (including the vertical flux proportional to dk[t] associated 131 ! with the lateral mixing, through the avt update) 132 ! dk[ avt dk[ (t,s) ] ] diffusive trends 133 134 ! 135 ! II.0 Matrix construction 136 ! ------------------------ 137 DO jn = 1, kjpt 101 ! ! ============= ! 102 DO jn = 1, kjpt ! tracer loop ! 103 ! ! ============= ! 138 104 ! 139 105 ! Matrix construction 140 ! ------------------------ 141 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 106 ! -------------------- 107 ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer 108 ! 109 IF( ( cdtype == 'TRA' .AND. ( ( jn == jp_tem ) .OR. ( jn == jp_sal .AND. lk_zdfddm ) ) ) .OR. & 110 & ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 111 ! 112 ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers 113 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ; zwt(:,:,2:jpk) = avt (:,:,2:jpk) 114 ELSE ; zwt(:,:,2:jpk) = fsavs(:,:,2:jpk) 115 ENDIF 116 zwt(:,:,1) = 0._wp 117 ! 142 118 #if defined key_ldfslp 143 IF( ln_traldf_grif ) THEN 119 ! isoneutral diffusion: add the contribution 120 IF( ln_traldf_grif ) THEN ! Griffies isoneutral diff 144 121 DO jk = 2, jpkm1 145 122 DO jj = 2, jpjm1 146 123 DO ji = fs_2, fs_jpim1 ! vector opt. 147 ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 148 zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 149 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 124 zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk) 150 125 END DO 151 126 END DO 152 127 END DO 153 ! update and save of avt (and avs if double diffusive mixing) 154 ELSE IF( l_traldf_rot ) THEN 128 ELSE IF( l_traldf_rot ) THEN ! standard isoneutral diff 155 129 DO jk = 2, jpkm1 156 130 DO jj = 2, jpjm1 157 131 DO ji = fs_2, fs_jpim1 ! vector opt. 158 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing 159 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 160 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 161 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 132 zwt(ji,jj,jk) = zwt(ji,jj,jk) + fsahtw(ji,jj,jk) & 133 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 134 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 162 135 END DO 163 136 END DO 164 137 END DO 165 ELSE ! no rotation but key_ldfslp defined166 zwt(:,:,:) = avt(:,:,:)167 138 ENDIF 168 #else169 ! No isopycnal diffusion170 zwt(:,:,:) = avt(:,:,:)171 139 #endif 172 ! Diagonal, inferior, superior (including the bottom boundary condition via avtmasked)140 ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked) 173 141 DO jk = 1, jpkm1 174 142 DO jj = 2, jpjm1 175 143 DO ji = fs_2, fs_jpim1 ! vector opt. 176 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point177 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point144 ze3ta = ( 1. - r_vvl ) + r_vvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 145 ze3tn = r_vvl + ( 1. - r_vvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 178 146 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 179 147 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) … … 182 150 END DO 183 151 END DO 152 ! 153 !! Matrix inversion from the first level 154 !!---------------------------------------------------------------------- 155 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 156 ! 157 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 158 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 159 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 160 ! ( ... )( ... ) ( ... ) 161 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 162 ! 163 ! m is decomposed in the product of an upper and lower triangular matrix. 164 ! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi. 165 ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal 166 ! and "superior" (above diagonal) components of the tridiagonal system. 167 ! The solution will be in the 4d array pta. 168 ! The 3d array zwt is used as a work space array. 169 ! En route to the solution pta is used a to evaluate the rhs and then 170 ! used as a work space array: its value is modified. 171 ! 184 172 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 173 ! done once for all passive tracers (so included in the IF instruction) 185 174 DO jj = 2, jpjm1 186 175 DO ji = fs_2, fs_jpim1 … … 196 185 END DO 197 186 ! 198 ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN199 #if defined key_ldfslp200 IF( ln_traldf_grif ) THEN201 DO jk = 2, jpkm1202 DO jj = 2, jpjm1203 DO ji = fs_2, fs_jpim1 ! vector opt.204 zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing205 ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing206 zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature)207 END DO208 END DO209 END DO210 ELSE IF( l_traldf_rot ) THEN211 DO jk = 2, jpkm1212 DO jj = 2, jpjm1213 DO ji = fs_2, fs_jpim1 ! vector opt.214 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing215 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) &216 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) )217 zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on salinity)218 END DO219 END DO220 END DO221 ELSE ! no rotation but key_ldfslp defined222 zwt(:,:,:) = fsavs(:,:,:)223 ENDIF224 #else225 ! No isopycnal diffusion226 zwt(:,:,:) = fsavs(:,:,:)227 #endif228 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked)229 DO jk = 1, jpkm1230 DO jj = 2, jpjm1231 DO ji = fs_2, fs_jpim1 ! vector opt.232 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point233 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point234 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) )235 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) )236 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk)237 END DO238 END DO239 END DO240 ! Surface boudary conditions241 DO jj = 2, jpjm1242 DO ji = fs_2, fs_jpim1 ! vector opt.243 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) ! after scale factor at T-point244 zwi(ji,jj,1) = 0._wp245 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1)246 END DO247 END DO248 !249 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k)250 DO jj = 2, jpjm1251 DO ji = fs_2, fs_jpim1252 zwt(ji,jj,1) = zwd(ji,jj,1)253 END DO254 END DO255 DO jk = 2, jpkm1256 DO jj = 2, jpjm1257 DO ji = fs_2, fs_jpim1258 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)259 END DO260 END DO261 END DO262 !263 187 END IF 264 265 ! II.1. Vertical diffusion on tracer 266 ! --------------------------- 267 188 ! 268 189 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 269 190 DO jj = 2, jpjm1 270 191 DO ji = fs_2, fs_jpim1 271 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1)272 ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1)192 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,1) 193 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t(ji,jj,1) 273 194 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 274 195 END DO … … 277 198 DO jj = 2, jpjm1 278 199 DO ji = fs_2, fs_jpim1 279 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk)280 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk)200 ze3tb = ( 1. - r_vvl ) + r_vvl * fse3t_b(ji,jj,jk) 201 ze3tn = ( 1. - r_vvl ) + r_vvl * fse3t (ji,jj,jk) 281 202 zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side 282 203 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) … … 285 206 END DO 286 207 287 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 208 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer) 288 209 DO jj = 2, jpjm1 289 210 DO ji = fs_2, fs_jpim1 … … 299 220 END DO 300 221 END DO 301 ! 302 END DO 222 ! ! ================= ! 223 END DO ! end tracer loop ! 224 ! ! ================= ! 303 225 ! 304 226 END SUBROUTINE tra_zdf_imp
Note: See TracChangeset
for help on using the changeset viewer.
