Changeset 457 for trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90
- Timestamp:
- 2006-05-10T19:01:19+02:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/TRA/trazdf_imp.F90
r258 r457 1 1 MODULE trazdf_imp 2 2 !!============================================================================== 3 !! 3 !! *** MODULE trazdf_imp *** 4 4 !! Ocean active tracers: vertical component of the tracer mixing trend 5 5 !!============================================================================== 6 7 !!---------------------------------------------------------------------- 8 !! tra_zdf_imp : update the tracer trend with the vertical diffusion 9 !! using an implicit time-stepping. 6 !! History : 7 !! 6.0 ! 90-10 (B. Blanke) Original code 8 !! 7.0 ! 91-11 (G. Madec) 9 !! ! 92-06 (M. Imbard) correction on tracer trend loops 10 !! ! 96-01 (G. Madec) statement function for e3 11 !! ! 97-05 (G. Madec) vertical component of isopycnal 12 !! ! 97-07 (G. Madec) geopotential diffusion in s-coord 13 !! ! 00-08 (G. Madec) double diffusive mixing 14 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 15 !! 9.0 ! 04-08 (C. Talandier) New trends organization 16 !! " " ! 06-11 (G. Madec) New step reorganisation 17 !!---------------------------------------------------------------------- 18 !! tra_zdf_imp : Update the tracer trend with the diagonal vertical 19 !! part of the mixing tensor. 20 !! Vector optimization, use k-j-i loops. 10 21 !!---------------------------------------------------------------------- 11 22 !! * Modules used 12 USE oce ! ocean dynamics and active tracers13 USE dom_oce ! ocean space and time domain 14 USE zdf_oce ! ocean vertical physics 23 USE oce ! ocean dynamics and tracers variables 24 USE dom_oce ! ocean space and time domain variables 25 USE zdf_oce ! ocean vertical physics variables 15 26 USE ldftra_oce ! ocean active tracers: lateral physics 16 USE zdfddm ! ocean vertical physics: double diffusion 17 USE zdfkpp ! KPP parameterisation 27 USE ldfslp ! lateral physics: slope of diffusion 18 28 USE trdmod ! ocean active tracers trends 19 29 USE trdmod_oce ! ocean variables trends 30 USE zdfddm ! ocean vertical physics: double diffusion 20 31 USE in_out_manager ! I/O manager 32 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 21 33 USE prtctl ! Print control 22 23 34 24 35 IMPLICIT NONE … … 26 37 27 38 !! * Routine accessibility 28 PUBLIC tra_zdf_imp ! routine called by step.F90 29 30 !! * Module variable 31 REAL(wp), DIMENSION(jpk) :: & 32 r2dt ! vertical profile of 2 x tracer time-step 39 PUBLIC tra_zdf_imp ! routine called by step.F90 33 40 34 41 !! * Substitutions 35 42 # include "domzgr_substitute.h90" 43 # include "ldftra_substitute.h90" 36 44 # include "zdfddm_substitute.h90" 37 !!---------------------------------------------------------------------- 38 !! OPA 9.0 , LOCEAN-IPSL (2005) 39 !! $Header$ 40 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 41 !!---------------------------------------------------------------------- 42 45 # include "vectopt_loop_substitute.h90" 46 !!---------------------------------------------------------------------- 47 !!---------------------------------------------------------------------- 48 !! OPA 9.0 , LOCEAN-IPSL (2005) 49 !!---------------------------------------------------------------------- 43 50 CONTAINS 44 45 SUBROUTINE tra_zdf_imp( kt )51 52 SUBROUTINE tra_zdf_imp( kt, p2dt ) 46 53 !!---------------------------------------------------------------------- 47 54 !! *** ROUTINE tra_zdf_imp *** 48 55 !! 49 !! ** Purpose : Compute the trend due to the vertical tracer mixing 50 !! using an implicit time stepping and add it to the general trend 51 !! of the tracer equations. 52 !! 53 !! ** Method : The vertical diffusion of tracers (t & s) is given by: 54 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(ta) ) 55 !! It is thus evaluated using a backward time scheme 56 !! ** Purpose : Compute the trend due to the vertical tracer diffusion 57 !! including the vertical component of lateral mixing (only for 2nd 58 !! order operator, for fourth order it is already computed and add 59 !! to the general trend in traldf.F) and add it to the general trend 60 !! of the tracer equations. 61 !! 62 !! ** Method : The vertical component of the lateral diffusive trends 63 !! is provided by a 2nd order operator rotated along neural or geo- 64 !! potential surfaces to which an eddy induced advection can be 65 !! added. It is computed using before fields (forward in time) and 66 !! isopycnal or geopotential slopes computed in routine ldfslp. 67 !! 68 !! Second part: vertical trend associated with the vertical physics 69 !! =========== (including the vertical flux proportional to dk[t] 70 !! associated with the lateral mixing, through the 71 !! update of avt) 72 !! The vertical diffusion of tracers (t & s) is given by: 73 !! difft = dz( avt dz(t) ) = 1/e3t dk+1( avt/e3w dk(t) ) 74 !! It is computed using a backward time scheme (t=ta). 56 75 !! Surface and bottom boundary conditions: no diffusive flux on 57 76 !! both tracers (bottom, applied through the masked field avt). 58 77 !! Add this trend to the general trend ta,sa : 59 !! ta = ta + dz( avt dz(t) ) 60 !! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T) 61 !! 62 !! ** Action : - Update (ta,sa) with the before vertical diffusion trend 63 !! - save the trends in (ttrd,strd) ('key_trdtra') 64 !! 65 !! History : 66 !! 6.0 ! 90-10 (B. Blanke) Original code 67 !! 7.0 ! 91-11 (G. Madec) 68 !! ! 92-06 (M. Imbard) correction on tracer trend loops 69 !! ! 96-01 (G. Madec) statement function for e3 70 !! ! 97-05 (G. Madec) vertical component of isopycnal 71 !! ! 97-07 (G. Madec) geopotential diffusion in s-coord 72 !! ! 00-08 (G. Madec) double diffusive mixing 73 !! 8.5 ! 02-08 (G. Madec) F90: Free form and module 74 !! 9.0 ! 04-08 (C. Talandier) New trends organization 75 !! 9.0 ! 05-01 (C. Ethe ) non-local flux in KPP vertical mixing scheme 78 !! ta = ta + dz( avt dz(t) ) 79 !! (sa = sa + dz( avs dz(t) ) if lk_zdfddm=T ) 80 !! 81 !! Third part: recover avt resulting from the vertical physics 82 !! ========== alone, for further diagnostics (for example to 83 !! compute the turbocline depth in zdfmxl.F90). 84 !! avt = zavt 85 !! (avs = zavs if lk_zdfddm=T ) 86 !! 87 !! ** Action : - Update (ta,sa) with before vertical diffusion trend 88 !! 76 89 !!--------------------------------------------------------------------- 77 !! * Modules used 78 USE oce, ONLY : ztdta => ua, & ! use ua as 3D workspace 79 ztdsa => va ! use va as 3D workspace 80 90 !! * Modules used 91 USE oce , ONLY : zwd => ua, & ! ua used as workspace 92 zws => va ! va " " 81 93 !! * Arguments 82 INTEGER, INTENT( in ) :: kt ! ocean time-step index 94 INTEGER, INTENT( in ) :: kt ! ocean time-step index 95 REAL(wp), DIMENSION(jpk), INTENT( in ) :: & 96 p2dt ! vertical profile of tracer time-step 83 97 84 98 !! * Local declarations 85 INTEGER :: ji, jj, jk ! dummy loop indices 86 INTEGER :: ikst, ikenm2, ikstp1 87 REAL(wp), DIMENSION(jpi,jpk) :: & 88 zwd, zws, zwi, & ! ??? 89 zwx, zwy, zwz, zwt ! ??? 99 INTEGER :: ji, jj, jk ! dummy loop indices 100 REAL(wp) :: zavi, zrhs ! temporary scalars 101 REAL(wp), DIMENSION(jpi,jpj,jpk) :: & 102 zwi, zwt, zavsi ! workspace arrays 90 103 !!--------------------------------------------------------------------- 91 104 92 93 ! 0. Local constant initialization 94 ! -------------------------------- 95 96 ! time step = 2 rdttra ex 97 IF( neuler == 0 .AND. kt == nit000 ) THEN 98 r2dt(:) = rdttra(:) ! restarting with Euler time stepping 99 ELSEIF( kt <= nit000 + 1) THEN 100 r2dt(:) = 2. * rdttra(:) ! leapfrog 105 IF( kt == nit000 ) THEN 106 IF(lwp)WRITE(numout,*) 107 IF(lwp)WRITE(numout,*) 'tra_zdf_imp : implicit vertical mixing (k-j-i loops)' 108 IF(lwp)WRITE(numout,*) '~~~~~~~~~~~ ' 109 zavi = 0.e0 ! avoid warning at compilation phase when lk_ldfslp=F 101 110 ENDIF 102 111 103 ! Save ta and sa trends 104 IF( l_trdtra ) THEN 105 ztdta(:,:,:) = ta(:,:,:) 106 ztdsa(:,:,:) = sa(:,:,:) 107 ENDIF 108 109 ! ! =============== 110 DO jj = 2, jpjm1 ! Vertical slab 111 ! ! =============== 112 ! 0. Matrix construction 113 ! ---------------------- 114 115 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 116 DO jk = 1, jpkm1 117 DO ji = 2, jpim1 118 zwi(ji,jk) = - r2dt(jk) * avt(ji,jj,jk ) & 119 / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 120 zws(ji,jk) = - r2dt(jk) * avt(ji,jj,jk+1) & 121 / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 122 zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk) 123 END DO 124 END DO 125 126 ! Surface boudary conditions 127 DO ji = 2, jpim1 128 zwi(ji,1) = 0.e0 129 zwd(ji,1) = 1. - zws(ji,1) 130 END DO 131 132 ! 1. Vertical diffusion on temperature 133 ! -------------------------=========== 134 135 ! Second member construction 136 #if defined key_zdfkpp 137 ! add non-local temperature flux ( in convective case only) 138 DO jk = 1, jpkm1 139 DO ji = 2, jpim1 140 zwy(ji,jk) = tb(ji,jj,jk) + r2dt(jk) * ta(ji,jj,jk) & 141 & - r2dt(jk) * ( ghats(ji,jj,jk) * avt(ji,jj,jk) - ghats(ji,jj,jk+1) * avt(ji,jj,jk+1) ) & 142 & * wt0(ji,jj) / fse3t(ji,jj,jk) 143 END DO 144 END DO 112 ! I. Local initialization 113 ! ----------------------- 114 zwd (1,:, : ) = 0.e0 ; zwd (jpi,:,:) = 0.e0 115 zws (1,:, : ) = 0.e0 ; zws (jpi,:,:) = 0.e0 116 zwi (1,:, : ) = 0.e0 ; zwi (jpi,:,:) = 0.e0 117 zwt (1,:, : ) = 0.e0 ; zwt (jpi,:,:) = 0.e0 118 zavsi(1,:, : ) = 0.e0 ; zavsi(jpi,:,:) = 0.e0 119 zwt (:,:,jpk) = 0.e0 ; zwt ( : ,:,1) = 0.e0 120 zavsi(:,:,jpk) = 0.e0 ; zavsi( : ,:,1) = 0.e0 121 122 ! II. Vertical trend associated with the vertical physics 123 ! ======================================================= 124 ! (including the vertical flux proportional to dk[t] associated 125 ! with the lateral mixing, through the avt update) 126 ! dk[ avt dk[ (t,s) ] ] diffusive trends 127 128 129 ! II.0 Matrix construction 130 ! ------------------------ 131 132 #if defined key_ldfslp 133 ! update and save of avt (and avs if double diffusive mixing) 134 DO jk = 2, jpkm1 135 DO jj = 2, jpjm1 136 DO ji = fs_2, fs_jpim1 ! vector opt. 137 zavi = fsahtw(ji,jj,jk) & ! vertical mixing coef. due to lateral mixing 138 & * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk) & 139 & + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) 140 zwt(ji,jj,jk) = avt(ji,jj,jk) + zavi ! zwt=avt+zavi (total vertical mixing coef. on temperature) 141 # if defined key_zdfddm 142 zavsi(ji,jj,jk) = fsavs(ji,jj,jk) + zavi ! dd mixing: zavsi = total vertical mixing coef. on salinity 143 # endif 144 END DO 145 END DO 146 END DO 147 148 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 149 DO jk = 1, jpkm1 150 DO jj = 2, jpjm1 151 DO ji = fs_2, fs_jpim1 ! vector opt. 152 zwi(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 153 zws(ji,jj,jk) = - p2dt(jk) * zwt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 154 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 155 END DO 156 END DO 157 END DO 145 158 #else 146 DO jk = 1, jpkm1 147 DO ji = 2, jpim1 148 zwy(ji,jk) = tb(ji,jj,jk) + r2dt(jk) * ta(ji,jj,jk) 149 END DO 150 END DO 159 ! Diagonal, inferior, superior (including the bottom boundary condition via avt masked) 160 DO jk = 1, jpkm1 161 DO jj = 2, jpjm1 162 DO ji = fs_2, fs_jpim1 ! vector opt. 163 zwi(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 164 zws(ji,jj,jk) = - p2dt(jk) * avt(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 165 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 166 END DO 167 END DO 168 END DO 151 169 #endif 152 170 153 ! Matrix inversion from the first level 154 ikst = 1 155 156 # include "zdf.matrixsolver.h90" 157 158 ! Save the masked temperature after in ta 159 ! (c a u t i o n: temperature not its trend, Leap-frog scheme done 160 ! it will not be done in tranxt) 161 DO jk = 1, jpkm1 162 DO ji = 2, jpim1 163 ta(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk) 164 END DO 165 END DO 166 167 168 ! 2. Vertical diffusion on salinity 169 ! -------------------------======== 171 ! Surface boudary conditions 172 DO jj = 2, jpjm1 173 DO ji = fs_2, fs_jpim1 ! vector opt. 174 zwi(ji,jj,1) = 0.e0 175 zwd(ji,jj,1) = 1. - zws(ji,jj,1) 176 END DO 177 END DO 178 179 180 ! II.1. Vertical diffusion on t 181 ! --------------------------- 182 183 !! Matrix inversion from the first level 184 !!---------------------------------------------------------------------- 185 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 186 ! 187 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 188 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 189 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 190 ! ( ... )( ... ) ( ... ) 191 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 192 ! 193 ! m is decomposed in the product of an upper and lower triangular matrix 194 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 195 ! The second member is in 2d array zwy 196 ! The solution is in 2d array zwx 197 ! The 3d arry zwt is a work space array 198 ! zwy is used and then used as a work space array : its value is modified! 199 200 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 201 DO jj = 2, jpjm1 202 DO ji = fs_2, fs_jpim1 203 zwt(ji,jj,1) = zwd(ji,jj,1) 204 END DO 205 END DO 206 DO jk = 2, jpkm1 207 DO jj = 2, jpjm1 208 DO ji = fs_2, fs_jpim1 209 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 210 END DO 211 END DO 212 END DO 213 214 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 215 DO jj = 2, jpjm1 216 DO ji = fs_2, fs_jpim1 217 ta(ji,jj,1) = tb(ji,jj,1) + p2dt(1) * ta(ji,jj,1) 218 END DO 219 END DO 220 DO jk = 2, jpkm1 221 DO jj = 2, jpjm1 222 DO ji = fs_2, fs_jpim1 223 zrhs = tb(ji,jj,jk) + p2dt(jk) * ta(ji,jj,jk) ! zrhs=right hand side 224 ta(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *ta(ji,jj,jk-1) 225 END DO 226 END DO 227 END DO 228 229 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 230 ! Save the masked temperature after in ta 231 ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 232 DO jj = 2, jpjm1 233 DO ji = fs_2, fs_jpim1 234 ta(ji,jj,jpkm1) = ta(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 235 END DO 236 END DO 237 DO jk = jpk-2, 1, -1 238 DO jj = 2, jpjm1 239 DO ji = fs_2, fs_jpim1 240 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) 241 END DO 242 END DO 243 END DO 244 245 ! II.2 Vertical diffusion on salinity 246 ! ----------------------------------- 170 247 171 248 #if defined key_zdfddm 172 ! Rebuild the Matrix as avt /= avs 173 174 ! Diagonal, inferior, superior 175 ! (including the bottom boundary condition via avs masked) 176 DO jk = 1, jpkm1 177 DO ji = 2, jpim1 178 zwi(ji,jk) = - r2dt(jk) * fsavs(ji,jj,jk ) & 179 /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 180 zws(ji,jk) = - r2dt(jk) * fsavs(ji,jj,jk+1) & 181 /( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 182 zwd(ji,jk) = 1. - zwi(ji,jk) - zws(ji,jk) 183 END DO 184 END DO 185 186 ! Surface boudary conditions 187 DO ji = 2, jpim1 188 zwi(ji,1) = 0.e0 189 zwd(ji,1) = 1. - zws(ji,1) 190 END DO 249 ! Rebuild the Matrix as avt /= avs 250 251 ! Diagonal, inferior, superior (including the bottom boundary condition via avs masked) 252 # if defined key_ldfslp 253 DO jk = 1, jpkm1 254 DO jj = 2, jpjm1 255 DO ji = fs_2, fs_jpim1 ! vector opt. 256 zwi(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 257 zws(ji,jj,jk) = - p2dt(jk) * zavsi(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 258 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 259 END DO 260 END DO 261 END DO 262 # else 263 DO jk = 1, jpkm1 264 DO jj = 2, jpjm1 265 DO ji = fs_2, fs_jpim1 ! vector opt. 266 zwi(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk ) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk ) ) 267 zws(ji,jj,jk) = - p2dt(jk) * avs(ji,jj,jk+1) / ( fse3t(ji,jj,jk) * fse3w(ji,jj,jk+1) ) 268 zwd(ji,jj,jk) = 1. - zwi(ji,jj,jk) - zws(ji,jj,jk) 269 END DO 270 END DO 271 END DO 272 # endif 273 274 ! Surface boudary conditions 275 DO jj = 2, jpjm1 276 DO ji = fs_2, fs_jpim1 ! vector opt. 277 zwi(ji,jj,1) = 0.e0 278 zwd(ji,jj,1) = 1. - zws(ji,jj,1) 279 END DO 280 END DO 191 281 #endif 192 ! Second member construction 193 #if defined key_zdfkpp 194 ! add non-local salinity flux ( in convective case only) 195 DO jk = 1, jpkm1 196 DO ji = 2, jpim1 197 zwy(ji,jk) = sb(ji,jj,jk) + r2dt(jk) * sa(ji,jj,jk) & 198 & - r2dt(jk) * ( ghats(ji,jj,jk) * fsavs(ji,jj,jk) - ghats(ji,jj,jk+1) * fsavs(ji,jj,jk+1) ) & 199 & * ws0(ji,jj) / fse3t(ji,jj,jk) 200 END DO 201 END DO 202 #else 203 DO jk = 1, jpkm1 204 DO ji = 2, jpim1 205 zwy(ji,jk) = sb(ji,jj,jk) + r2dt(jk) * sa(ji,jj,jk) 206 END DO 207 END DO 208 #endif 209 210 ! Matrix inversion from the first level 211 ikst = 1 212 213 # include "zdf.matrixsolver.h90" 214 215 216 ! Save the masked salinity after in sa 217 ! (c a u t i o n: salinity not its trend, Leap-frog scheme done 218 ! it will not be done in tranxt) 219 DO jk = 1, jpkm1 220 DO ji = 2, jpim1 221 sa(ji,jj,jk) = zwx(ji,jk) * tmask(ji,jj,jk) 222 END DO 223 END DO 224 225 ! ! =============== 226 END DO ! End of slab 227 ! ! =============== 228 229 ! save the trends for diagnostic 230 ! Compute and save the vertical diffusive temperature & salinity trends 231 IF( l_trdtra ) THEN 232 ! compute the vertical diffusive trends in substracting the previous 233 ! trends ztdta()/ztdsa() to the new one computed (dT/dt or dS/dt) 234 ! with the new temperature/salinity ta/sa 235 DO jk = 1, jpkm1 236 ztdta(:,:,jk) = ( ( ta(:,:,jk) - tb(:,:,jk) ) / r2dt(jk) ) & ! new trend 237 & - ztdta(:,:,jk) ! old trend 238 ztdsa(:,:,jk) = ( ( sa(:,:,jk) - sb(:,:,jk) ) / r2dt(jk) ) & ! new trend 239 & - ztdsa(:,:,jk) ! old trend 240 END DO 241 242 CALL trd_mod(ztdta, ztdsa, jpttdzdf, 'TRA', kt) 243 ENDIF 244 245 IF(ln_ctl) THEN ! print mean trends (used for debugging) 246 CALL prt_ctl(tab3d_1=ta, clinfo1=' zdf - Ta: ', mask1=tmask, & 247 & tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra') 248 ENDIF 282 283 284 !! Matrix inversion from the first level 285 !!---------------------------------------------------------------------- 286 ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk ) 287 ! 288 ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 ) 289 ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 ) 290 ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 ) 291 ! ( ... )( ... ) ( ... ) 292 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 293 ! 294 ! m is decomposed in the product of an upper and lower triangular 295 ! matrix 296 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 297 ! The second member is in 2d array zwy 298 ! The solution is in 2d array zwx 299 ! The 3d arry zwt is a work space array 300 ! zwy is used and then used as a work space array : its value is modified! 301 302 ! first recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) 303 DO jj = 2, jpjm1 304 DO ji = fs_2, fs_jpim1 305 zwt(ji,jj,1) = zwd(ji,jj,1) 306 END DO 307 END DO 308 DO jk = 2, jpkm1 309 DO jj = 2, jpjm1 310 DO ji = fs_2, fs_jpim1 311 zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1) 312 END DO 313 END DO 314 END DO 315 316 ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1 317 DO jj = 2, jpjm1 318 DO ji = fs_2, fs_jpim1 319 sa(ji,jj,1) = sb(ji,jj,1) + p2dt(1) * sa(ji,jj,1) 320 END DO 321 END DO 322 DO jk = 2, jpkm1 323 DO jj = 2, jpjm1 324 DO ji = fs_2, fs_jpim1 325 zrhs = sb(ji,jj,jk) + p2dt(jk) * sa(ji,jj,jk) ! zrhs=right hand side 326 sa(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *sa(ji,jj,jk-1) 327 END DO 328 END DO 329 END DO 330 331 ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk 332 ! Save the masked temperature after in ta 333 ! (c a u t i o n: temperature not its trend, Leap-frog scheme done it will not be done in tranxt) 334 DO jj = 2, jpjm1 335 DO ji = fs_2, fs_jpim1 336 sa(ji,jj,jpkm1) = sa(ji,jj,jpkm1) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1) 337 END DO 338 END DO 339 DO jk = jpk-2, 1, -1 340 DO jj = 2, jpjm1 341 DO ji = fs_2, fs_jpim1 342 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) 343 END DO 344 END DO 345 END DO 249 346 250 347 END SUBROUTINE tra_zdf_imp
Note: See TracChangeset
for help on using the changeset viewer.
