Changeset 2528 for trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
- Timestamp:
- 2010-12-27T18:33:53+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/TRA/trazdf_imp.F90
- Property svn:eol-style deleted
r1517 r2528 2 2 !!====================================================================== 3 3 !! *** MODULE trazdf_imp *** 4 !! Ocean activetracers: vertical component of the tracer mixing trend4 !! Ocean tracers: vertical component of the tracer mixing trend 5 5 !!====================================================================== 6 6 !! History : OPA ! 1990-10 (B. Blanke) Original code … … 14 14 !! 2.0 ! 2006-11 (G. Madec) New step reorganisation 15 15 !! 3.2 ! 2009-03 (G. Madec) heat and salt content trends 16 !! 3.3 ! 2010-06 (C. Ethe, G. Madec) Merge TRA-TRC 16 17 !!---------------------------------------------------------------------- 17 18 … … 23 24 USE dom_oce ! ocean space and time domain variables 24 25 USE zdf_oce ! ocean vertical physics variables 26 USE trc_oce ! share passive tracers/ocean variables 27 USE domvvl ! variable volume 25 28 USE ldftra_oce ! ocean active tracers: lateral physics 29 USE ldftra ! lateral mixing type 26 30 USE ldfslp ! lateral physics: slope of diffusion 27 USE trdmod ! ocean active tracers trends28 USE trdmod_oce ! ocean variables trends29 31 USE zdfddm ! ocean vertical physics: double diffusion 32 USE traldf_iso_grif ! active tracers: Griffies operator 30 33 USE in_out_manager ! I/O manager 31 34 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 32 USE prtctl ! Print control33 USE domvvl ! variable volume34 USE ldftra ! lateral mixing type35 35 36 36 IMPLICIT NONE … … 45 45 # include "vectopt_loop_substitute.h90" 46 46 !!---------------------------------------------------------------------- 47 !! NEMO/OPA 3. 2 , LOCEAN-IPSL (2009)47 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 48 48 !! $Id$ 49 !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)49 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 50 50 !!---------------------------------------------------------------------- 51 51 CONTAINS 52 53 SUBROUTINE tra_zdf_imp( kt, p2dt )52 53 SUBROUTINE tra_zdf_imp( kt, cdtype, p2dt, ptb, pta, kjpt ) 54 54 !!---------------------------------------------------------------------- 55 55 !! *** ROUTINE tra_zdf_imp *** … … 71 71 !! associated with the lateral mixing, through the 72 72 !! update of avt) 73 !! The vertical diffusion of t racers (t & s)is given by:73 !! The vertical diffusion of the tracer t is given by: 74 74 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 75 75 !! It is computed using a backward time scheme (t=ta). … … 78 78 !! Add this trend to the general trend ta,sa : 79 79 !! ta = ta + dz( avt dz(t) ) 80 !! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T ) 80 !! if lk_zdfddm=T, use avs for salinity or for passive tracers 81 !! (sa = sa + dz( avs dz(t) ) 81 82 !! 82 83 !! Third part: recover avt resulting from the vertical physics 83 84 !! ========== alone, for further diagnostics (for example to 84 85 !! compute the turbocline depth in zdfmxl.F90). 85 !! avt = zavt86 !! if lk_zdfddm=T, use avt = zavt 86 87 !! (avs = zavs if lk_zdfddm=T ) 87 88 !! 88 !! ** Action : - Update (ta,sa) with before vertical diffusion trend 89 !! 89 !! ** Action : - Update (ta) with before vertical diffusion trend 90 90 !!--------------------------------------------------------------------- 91 91 USE oce , ONLY : zwd => ua ! ua used as workspace 92 92 USE oce , ONLY : zws => va ! va - - 93 !! 94 INTEGER , INTENT(in) :: kt ! ocean time-step index 95 REAL(wp), DIMENSION(jpk), INTENT(in) :: p2dt ! vertical profile of tracer time-step 96 !! 97 INTEGER :: ji, jj, jk ! dummy loop indices 98 REAL(wp) :: zavi, zrhs, znvvl ! temporary scalars 99 REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors 100 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt, zavsi ! workspace arrays 93 !! 94 INTEGER , INTENT(in ) :: kt ! ocean time-step index 95 CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator) 96 INTEGER , INTENT(in ) :: kjpt ! number of tracers 97 REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step 98 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before and now tracer fields 99 REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend 100 !! 101 INTEGER :: ji, jj, jk, jn ! dummy loop indices 102 REAL(wp) :: zavi, zrhs, znvvl ! local scalars 103 REAL(wp) :: ze3tb, ze3tn, ze3ta ! variable vertical scale factors 104 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwi, zwt ! workspace arrays 101 105 !!--------------------------------------------------------------------- 102 106 103 IF( kt == nit000 ) THEN107 IF( kt == nit000 ) THEN 104 108 IF(lwp)WRITE(numout,*) 105 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing '109 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing on ', cdtype 106 110 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 107 zavi = 0. e0! avoid warning at compilation phase when lk_ldfslp=F111 zavi = 0._wp ! avoid warning at compilation phase when lk_ldfslp=F 108 112 ENDIF 109 113 ! 110 114 ! I. Local initialization 111 115 ! ----------------------- 112 zwd (1,:, : ) = 0.e0 ; zwd (jpi,:,:) = 0.e0 113 zws (1,:, : ) = 0.e0 ; zws (jpi,:,:) = 0.e0 114 zwi (1,:, : ) = 0.e0 ; zwi (jpi,:,:) = 0.e0 115 zwt (1,:, : ) = 0.e0 ; zwt (jpi,:,:) = 0.e0 116 zavsi(1,:, : ) = 0.e0 ; zavsi(jpi,:,:) = 0.e0 117 zwt (:,:,jpk) = 0.e0 ; zwt ( : ,:,1) = 0.e0 118 zavsi(:,:,jpk) = 0.e0 ; zavsi( : ,:,1) = 0.e0 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 119 121 120 122 ! I.1 Variable volume : to take into account vertical variable vertical scale factors 121 123 ! ------------------- 122 IF( lk_vvl ) THEN ; znvvl = 1. 123 ELSE ; znvvl = 0. e0124 IF( lk_vvl ) THEN ; znvvl = 1._wp 125 ELSE ; znvvl = 0._wp 124 126 ENDIF 125 127 … … 130 132 ! dk[ avt dk[ (t,s) ] ] diffusive trends 131 133 132 134 ! 133 135 ! II.0 Matrix construction 134 136 ! ------------------------ 135 137 DO jn = 1, kjpt 138 ! 139 ! Matrix construction 140 ! ------------------------ 141 IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN 136 142 #if defined key_ldfslp 137 ! update and save of avt (and avs if double diffusive mixing) 138 IF( l_traldf_rot ) THEN 139 DO jk = 2, jpkm1 143 IF( ln_traldf_grif ) THEN 144 DO jk = 2, jpkm1 145 DO jj = 2, jpjm1 146 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) 150 END DO 151 END DO 152 END DO 153 ! update and save of avt (and avs if double diffusive mixing) 154 ELSE IF( l_traldf_rot ) THEN 155 DO jk = 2, jpkm1 156 DO jj = 2, jpjm1 157 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) 162 END DO 163 END DO 164 END DO 165 ELSE ! no rotation but key_ldfslp defined 166 zwt(:,:,:) = avt(:,:,:) 167 ENDIF 168 #else 169 ! No isopycnal diffusion 170 zwt(:,:,:) = avt(:,:,:) 171 #endif 172 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 173 DO jk = 1, jpkm1 174 DO jj = 2, jpjm1 175 DO ji = fs_2, fs_jpim1 ! vector opt. 176 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 177 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 178 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 179 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 180 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 181 END DO 182 END DO 183 END DO 184 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 185 DO jj = 2, jpjm1 186 DO ji = fs_2, fs_jpim1 187 zwt(ji,jj,1) = zwd(ji,jj,1) 188 END DO 189 END DO 190 DO jk = 2, jpkm1 191 DO jj = 2, jpjm1 192 DO ji = fs_2, fs_jpim1 193 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 194 END DO 195 END DO 196 END DO 197 ! 198 ELSE IF( ( cdtype == 'TRA' .AND. jn == jp_sal ) .OR. ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN 199 #if defined key_ldfslp 200 IF( ln_traldf_grif ) THEN 201 DO jk = 2, jpkm1 202 DO jj = 2, jpjm1 203 DO ji = fs_2, fs_jpim1 ! vector opt. 204 zavi = ah_wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 205 ! zavi = fsahtw(ji,jj,jk) * wslp2(ji,jj,jk) ! vertical mixing coef. due to lateral mixing 206 zwt(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 207 END DO 208 END DO 209 END DO 210 ELSE IF( l_traldf_rot ) THEN 211 DO jk = 2, jpkm1 212 DO jj = 2, jpjm1 213 DO ji = fs_2, fs_jpim1 ! vector opt. 214 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing 215 & * ( 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 DO 219 END DO 220 END DO 221 ELSE ! no rotation but key_ldfslp defined 222 zwt(:,:,:) = fsavs(:,:,:) 223 ENDIF 224 #else 225 ! No isopycnal diffusion 226 zwt(:,:,:) = fsavs(:,:,:) 227 #endif 228 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 229 DO jk = 1, jpkm1 230 DO jj = 2, jpjm1 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,jk) ! after scale factor at T-point 233 ze3tn = znvvl + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) ! now scale factor at T-point 234 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 DO 238 END DO 239 END DO 240 ! Surface boudary conditions 140 241 DO jj = 2, jpjm1 141 242 DO ji = fs_2, fs_jpim1 ! vector opt. 142 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing 143 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 144 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 145 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 146 # if defined key_zdfddm 147 zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity 148 # endif 149 END DO 150 END DO 151 END DO 152 ELSE ! no rotation but key_ldfslp defined 153 zwt (:,:,:) = avt(:,:,:) 154 # if defined key_zdfddm 155 zavsi(:,:,:) = avs(:,:,:) ! avs /= avt (double diffusion mixing) 156 # endif 157 ENDIF 158 #else 159 ! No isopycnal diffusion 160 zwt(:,:,:) = avt(:,:,:) 161 # if defined key_zdfddm 162 zavsi(:,:,:) = avs(:,:,:) 163 # endif 164 165 #endif 166 167 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 168 DO jk = 1, jpkm1 169 DO jj = 2, jpjm1 170 DO ji = fs_2, fs_jpim1 ! vector opt. 171 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 172 & + znvvl * fse3t_a(ji,jj,jk) 173 ze3tn = znvvl & ! now scale factor at T-point 174 & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 175 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 176 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 177 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 178 END DO 179 END DO 180 END DO 181 182 ! Surface boudary conditions 183 DO jj = 2, jpjm1 184 DO ji = fs_2, fs_jpim1 ! vector opt. 185 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 186 & + znvvl * fse3t_a(ji,jj,1) 187 zwi(ji,jj,1) = 0.e0 188 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 189 END DO 190 END DO 191 192 193 ! II.1. Vertical diffusion on t 194 ! --------------------------- 195 196 !! Matrix inversion from the first level 197 !!---------------------------------------------------------------------- 198 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 199 ! 200 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 201 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 202 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 203 ! ( ... )( ... ) ( ... ) 204 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 205 ! 206 ! m is decomposed in the product of an upper and lower triangular matrix 207 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 208 ! The second member is in 2d array zwy 209 ! The solution is in 2d array zwx 210 ! The 3d arry zwt is a work space array 211 ! zwy is used and then used as a work space array : its value is modified! 212 213 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 214 DO jj = 2, jpjm1 215 DO ji = fs_2, fs_jpim1 216 zwt(ji,jj,1) = zwd(ji,jj,1) 217 END DO 218 END DO 219 DO jk = 2, jpkm1 243 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) ! after scale factor at T-point 244 zwi(ji,jj,1) = 0._wp 245 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 246 END DO 247 END DO 248 ! 249 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 250 DO jj = 2, jpjm1 251 DO ji = fs_2, fs_jpim1 252 zwt(ji,jj,1) = zwd(ji,jj,1) 253 END DO 254 END DO 255 DO jk = 2, jpkm1 256 DO jj = 2, jpjm1 257 DO ji = fs_2, fs_jpim1 258 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1) 259 END DO 260 END DO 261 END DO 262 ! 263 END IF 264 265 ! II.1. Vertical diffusion on tracer 266 ! --------------------------- 267 268 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 220 269 DO jj = 2, jpjm1 221 270 DO ji = fs_2, fs_jpim1 222 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 271 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 272 ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 273 pta(ji,jj,1,jn) = ze3tb * ptb(ji,jj,1,jn) + p2dt(1) * ze3tn * pta(ji,jj,1,jn) 223 274 END DO 224 275 END DO 225 END DO 226 227 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 230 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,1) 231 ze3tn = ( 1. - znvvl ) + znvvl * fse3t(ji,jj,1) 232 ta(ji,jj,1) = ze3tb * tb(ji,jj,1) + p2dt(1) * ze3tn * ta(ji,jj,1) 276 DO jk = 2, jpkm1 277 DO jj = 2, jpjm1 278 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) 281 zrhs = ze3tb * ptb(ji,jj,jk,jn) + p2dt(jk) * ze3tn * pta(ji,jj,jk,jn) ! zrhs=right hand side 282 pta(ji,jj,jk,jn) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pta(ji,jj,jk-1,jn) 283 END DO 284 END DO 233 285 END DO 234 END DO 235 DO jk = 2, jpkm1286 287 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 236 288 DO jj = 2, jpjm1 237 289 DO ji = fs_2, fs_jpim1 238 ze3tb = ( 1. - znvvl ) + znvvl * fse3t_b(ji,jj,jk) 239 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) 240 zrhs = ze3tb * tb(ji,jj,jk) + p2dt(jk) * ze3tn * ta(ji,jj,jk) ! zrhs=right hand side 241 ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * ta(ji,jj,jk-1) 290 pta(ji,jj,jpkm1,jn) = pta(ji,jj,jpkm1,jn) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 242 291 END DO 243 292 END DO 244 END DO 245 246 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 247 ! Save the masked temperature after in ta 248 ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 249 DO jj = 2, jpjm1 250 DO ji = fs_2, fs_jpim1 251 ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 293 DO jk = jpk-2, 1, -1 294 DO jj = 2, jpjm1 295 DO ji = fs_2, fs_jpim1 296 pta(ji,jj,jk,jn) = ( pta(ji,jj,jk,jn) - zws(ji,jj,jk) * pta(ji,jj,jk+1,jn) ) & 297 & / zwt(ji,jj,jk) * tmask(ji,jj,jk) 298 END DO 299 END DO 252 300 END DO 253 END DO 254 DO jk = jpk-2, 1, -1 255 DO jj = 2, jpjm1 256 DO ji = fs_2, fs_jpim1 257 ta(ji,jj,jk) = ( ta(ji,jj,jk) - zws(ji,jj,jk) * ta(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) 258 END DO 259 END DO 260 END DO 261 262 ! II.2 Vertical diffusion on salinity 263 ! ----------------------------------- 264 265 #if defined key_zdfddm 266 ! Rebuild the Matrix as avt /= avs 267 268 ! Diagonal, inferior, superior (including the bottom boundary condition via avs masked) 269 DO jk = 1, jpkm1 270 DO jj = 2, jpjm1 271 DO ji = fs_2, fs_jpim1 ! vector opt. 272 ze3ta = ( 1. - znvvl ) & ! after scale factor at T-point 273 & + znvvl * fse3t_a(ji,jj,jk) 274 ze3tn = znvvl & ! now scale factor at T-point 275 & + ( 1. - znvvl ) * fse3t_n(ji,jj,jk) 276 zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( ze3tn * fse3w(ji,jj,jk ) ) 277 zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( ze3tn * fse3w(ji,jj,jk+1) ) 278 zwd(ji,jj,jk) = ze3ta - zwi(ji,jj,jk) - zws(ji,jj,jk) 279 END DO 280 END DO 281 END DO 282 283 ! Surface boudary conditions 284 DO jj = 2, jpjm1 285 DO ji = fs_2, fs_jpim1 ! vector opt. 286 ze3ta = ( 1. - znvvl ) + znvvl * fse3t_a(ji,jj,1) 287 zwi(ji,jj,1) = 0.e0 288 zwd(ji,jj,1) = ze3ta - zws(ji,jj,1) 289 END DO 290 END DO 291 #endif 292 293 294 !! Matrix inversion from the first level 295 !!---------------------------------------------------------------------- 296 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 297 ! 298 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 299 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 300 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 301 ! ( ... )( ... ) ( ... ) 302 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 303 ! 304 ! m is decomposed in the product of an upper and lower triangular 305 ! matrix 306 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 307 ! The second member is in 2d array zwy 308 ! The solution is in 2d array zwx 309 ! The 3d arry zwt is a work space array 310 ! zwy is used and then used as a work space array : its value is modified! 311 312 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 313 DO jj = 2, jpjm1 314 DO ji = fs_2, fs_jpim1 315 zwt(ji,jj,1) = zwd(ji,jj,1) 316 END DO 317 END DO 318 DO jk = 2, jpkm1 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 321 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 322 END DO 323 END DO 324 END DO 325 326 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 329 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point 330 & + znvvl * fse3t_b(ji,jj,1) 331 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,1) ! now scale factor at T-point 332 sa(ji,jj,1) = ze3tb * sb(ji,jj,1) + p2dt(1) * ze3tn * sa(ji,jj,1) 333 END DO 334 END DO 335 DO jk = 2, jpkm1 336 DO jj = 2, jpjm1 337 DO ji = fs_2, fs_jpim1 338 ze3tb = ( 1. - znvvl ) & ! before scale factor at T-point 339 & + znvvl * fse3t_b(ji,jj,jk) 340 ze3tn = ( 1. - znvvl ) + znvvl * fse3t (ji,jj,jk) ! now scale factor at T-point 341 zrhs = ze3tb * sb(ji,jj,jk) + p2dt(jk) * ze3tn * sa(ji,jj,jk) ! zrhs=right hand side 342 sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) 343 END DO 344 END DO 345 END DO 346 347 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 348 ! Save the masked temperature after in ta 349 ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 350 DO jj = 2, jpjm1 351 DO ji = fs_2, fs_jpim1 352 sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 353 END DO 354 END DO 355 DO jk = jpk-2, 1, -1 356 DO jj = 2, jpjm1 357 DO ji = fs_2, fs_jpim1 358 sa(ji,jj,jk) = ( sa(ji,jj,jk) - zws(ji,jj,jk) * sa(ji,jj,jk+1) ) / zwt(ji,jj,jk) * tmask(ji,jj,jk) 359 END DO 360 END DO 301 ! 361 302 END DO 362 303 !
Note: See TracChangeset
for help on using the changeset viewer.
