Changeset 10222 for NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90
- Timestamp:
- 2018-10-25T11:42:23+02:00 (5 years ago)
- File:
-
- 1 edited
-
NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/TOP/PISCES/SED/sedmat.F90
r5215 r10222 1 1 MODULE sedmat 2 #if defined key_sed3 2 !!====================================================================== 4 3 !! *** MODULE sedmat *** … … 9 8 10 9 USE sed ! sediment global variable 10 USE lib_mpp ! distribued memory computing library 11 11 12 12 13 IMPLICIT NONE … … 25 26 CONTAINS 26 27 27 SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, ps ol)28 SUBROUTINE sed_mat_dsr( nvar, ndim, nlev, preac, psms, psol, dtsed_in ) 28 29 !!--------------------------------------------------------------------- 29 30 !! *** ROUTINE sed_mat_dsr *** … … 48 49 !!---------------------------------------------------------------------- 49 50 !! * Arguments 50 INTEGER , INTENT(in) :: nvar ! number of variable s51 INTEGER , INTENT(in) :: nvar ! number of variable 51 52 INTEGER , INTENT(in) :: ndim ! number of points 52 53 INTEGER , INTENT(in) :: nlev ! number of sediment levels 53 54 54 REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(in ) :: preac ! reaction rates 55 REAL(wp), DIMENSION(ndim,nlev,nvar), INTENT(inout) :: psol ! solution ( undersaturation values ) 55 REAL(wp), DIMENSION(ndim,nlev), INTENT(in ) :: preac ! reaction rates 56 REAL(wp), DIMENSION(ndim,nlev), INTENT(in ) :: psms ! reaction rates 57 REAL(wp), DIMENSION(ndim,nlev), INTENT(inout) :: psol ! solution ( undersaturation values ) 58 REAL(wp), INTENT(in) :: dtsed_in 56 59 57 60 !---Local declarations … … 67 70 !---------------------------------------------------------------------- 68 71 72 IF( ln_timing ) CALL timing_start('sed_mat_dsr') 69 73 70 74 ! Computation left hand side of linear system of … … 73 77 74 78 79 jn = nvar 75 80 ! first sediment level 76 81 DO ji = 1, ndim 77 aplus = ( ( volw3d(ji,1) / dz3d(ji,1) ) + &78 ( volw3d(ji,2) / dz3d(ji,2) ) ) / 2.82 aplus = ( ( volw3d(ji,1) / ( dz3d(ji,1) ) ) + & 83 ( volw3d(ji,2) / ( dz3d(ji,2) ) ) ) / 2. 79 84 dxplus = ( dz3d(ji,1) + dz3d(ji,2) ) / 2. 80 rplus = ( dtsed / volw3d(ji,1) ) * diff(1) * aplus / dxplus85 rplus = ( dtsed_in / ( volw3d(ji,1) ) ) * diff(ji,1,jn) * aplus / dxplus 81 86 82 87 za(ji,1) = 0. … … 87 92 DO jk = 2, nlev - 1 88 93 DO ji = 1, ndim 89 aminus = ( ( volw3d(ji,jk-1) / dz3d(ji,jk-1) ) + &90 & ( volw3d(ji,jk ) / dz3d(ji,jk) ) ) / 2.94 aminus = ( ( volw3d(ji,jk-1) / ( dz3d(ji,jk-1) ) ) + & 95 & ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) ) / 2. 91 96 dxminus = ( dz3d(ji,jk-1) + dz3d(ji,jk) ) / 2. 92 97 93 aplus = ( ( volw3d(ji,jk ) / dz3d(ji,jk) ) + &94 & ( volw3d(ji,jk+1) / dz3d(ji,jk+1) ) ) / 2.98 aplus = ( ( volw3d(ji,jk ) / ( dz3d(ji,jk ) ) ) + & 99 & ( volw3d(ji,jk+1) / ( dz3d(ji,jk+1) ) ) ) / 2. 95 100 dxplus = ( dz3d(ji,jk) + dz3d(ji,jk+1) ) / 2 96 !97 rminus = ( dtsed / volw3d(ji,jk) ) * diff(jk-1) * aminus / dxminus98 rplus = ( dtsed / volw3d(ji,jk) ) * diff(jk) * aplus / dxplus99 !101 ! 102 rminus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk-1,jn) * aminus / dxminus 103 rplus = ( dtsed_in / volw3d(ji,jk) ) * diff(ji,jk,jn) * aplus / dxplus 104 ! 100 105 za(ji,jk) = -rminus 101 106 zb(ji,jk) = 1. + rminus + rplus … … 105 110 106 111 DO ji = 1, ndim 107 aminus = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + &108 112 aminus = ( ( volw3d(ji,nlev-1) / dz3d(ji,nlev-1) ) + & 113 & ( volw3d(ji,nlev) / dz3d(ji,nlev) ) ) / 2. 109 114 dxminus = ( dz3d(ji,nlev-1) + dz3d(ji,nlev) ) / 2. 110 rminus = ( dtsed / volw3d(ji,nlev) ) * diff(nlev-1) * aminus / dxminus115 rminus = ( dtsed_in / volw3d(ji,nlev) ) * diff(ji,nlev-1,jn) * aminus / dxminus 111 116 ! 112 117 za(ji,nlev) = -rminus … … 119 124 ! ----------------------------------------------- 120 125 121 DO jn = 1, nvar122 123 zr (:,:) = psol(:,:,jn)124 zbet(: ) = zb(:,1) + preac(:,1,jn)125 psol(:,1,jn) = zr(:,1) / zbet(:) 126 zr (:,:) = psol(:,:) + psms(:,:) 127 zb (:,:) = zb(:,:) + preac(:,:) 128 zbet(: ) = zb(:,1) 129 psol(:,1) = zr(:,1) / zbet(:) 130 126 131 ! 127 128 129 zgamm(ji,jk)= zc(ji,jk-1) / zbet(ji)130 zbet(ji) = ( zb(ji,jk) + preac(ji,jk,jn)) - za(ji,jk) * zgamm(ji,jk)131 psol(ji,jk,jn) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1,jn) ) / zbet(ji)132 133 132 DO jk = 2, nlev 133 DO ji = 1, ndim 134 zgamm(ji,jk) = zc(ji,jk-1) / zbet(ji) 135 zbet(ji) = zb(ji,jk) - za(ji,jk) * zgamm(ji,jk) 136 psol(ji,jk) = ( zr(ji,jk) - za(ji,jk) * psol(ji,jk-1) ) / zbet(ji) 137 END DO 138 ENDDO 134 139 ! 135 DO jk = nlev - 1, 1, -1 136 DO ji = 1,ndim 137 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(ji,jk+1) * psol(ji,jk+1,jn) 138 END DO 139 ENDDO 140 140 DO jk = nlev - 1, 1, -1 141 DO ji = 1,ndim 142 psol(ji,jk) = psol(ji,jk) - zgamm(ji,jk+1) * psol(ji,jk+1) 143 END DO 141 144 ENDDO 142 145 146 IF( ln_timing ) CALL timing_stop('sed_mat_dsr') 147 143 148 144 149 END SUBROUTINE sed_mat_dsr 145 150 146 147 SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol ) 151 SUBROUTINE sed_mat_btb( nvar, ndim, nlev, psol, dtsed_in ) 148 152 !!--------------------------------------------------------------------- 149 153 !! *** ROUTINE sed_mat_btb *** … … 170 174 psol ! solution 171 175 176 REAL(wp), INTENT(in) :: dtsed_in 177 172 178 !---Local declarations 173 179 INTEGER :: & … … 192 198 193 199 200 IF( ln_timing ) CALL timing_start('sed_mat_btb') 201 194 202 ! first sediment level 195 aplus = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 196 dxplus = ( dz(2) + dz(3) ) / 2. 197 rplus = ( dtsed / vols(2) ) * db * aplus / dxplus 198 199 za(1) = 0. 200 zb(1) = 1. + rplus 201 zc(1) = -rplus 203 DO ji = 1, ndim 204 aplus = ( ( vols(2) / dz(2) ) + ( vols(3) / dz(3) ) ) / 2. 205 dxplus = ( dz(2) + dz(3) ) / 2. 206 rplus = ( dtsed_in / vols(2) ) * db(ji,2) * aplus / dxplus 207 208 za(1) = 0. 209 zb(1) = 1. + rplus 210 zc(1) = -rplus 202 211 203 212 204 DO jk = 2, nlev - 1 205 aminus = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 206 dxminus = ( dz(jk) + dz(jk+1) ) / 2. 207 rminus = ( dtsed / vols(jk+1) ) * db * aminus / dxminus 208 ! 209 aplus = ( ( vols(jk+1) / dz(jk+1 ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 210 dxplus = ( dz(jk+1) + dz(jk+2) ) / 2. 211 rplus = ( dtsed / vols(jk+1) ) * db * aplus / dxplus 212 ! 213 za(jk) = -rminus 214 zb(jk) = 1. + rminus + rplus 215 zc(jk) = -rplus 216 ENDDO 217 218 aminus = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 219 dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 220 rminus = ( dtsed / vols(nlev+1) ) * db * aminus / dxminus 213 DO jk = 2, nlev - 1 214 aminus = ( ( vols(jk) / dz(jk) ) + ( vols(jk+1) / dz(jk+1) ) ) / 2. 215 dxminus = ( dz(jk) + dz(jk+1) ) / 2. 216 rminus = ( dtsed_in / vols(jk+1) ) * db(ji,jk) * aminus / dxminus 217 ! 218 aplus = ( ( vols(jk+1) / dz(jk+1 ) ) + ( vols(jk+2) / dz(jk+2) ) ) / 2. 219 dxplus = ( dz(jk+1) + dz(jk+2) ) / 2. 220 rplus = ( dtsed_in / vols(jk+1) ) * db(ji,jk+1) * aplus / dxplus 221 ! 222 za(jk) = -rminus 223 zb(jk) = 1. + rminus + rplus 224 zc(jk) = -rplus 225 ENDDO 226 227 aminus = ( ( vols(nlev) / dz(nlev) ) + ( vols(nlev+1) / dz(nlev+1) ) ) / 2. 228 dxminus = ( dz(nlev) + dz(nlev+1) ) / 2. 229 rminus = ( dtsed_in / vols(nlev+1) ) * db(ji,nlev) * aminus / dxminus 230 ! 231 za(nlev) = -rminus 232 zb(nlev) = 1. + rminus 233 zc(nlev) = 0. 234 235 236 ! solves tridiagonal system of linear equations 237 ! ----------------------------------------------- 238 DO jn = 1, nvar 239 240 DO jk = 1, nlev 241 zr (ji,jk) = psol(ji,jk,jn) 242 END DO 243 zbet = zb(1) 244 psol(ji,1,jn) = zr(ji,1) / zbet 245 ! 246 DO jk = 2, nlev 247 zgamm(jk) = zc(jk-1) / zbet 248 zbet = zb(jk) - za(jk) * zgamm(jk) 249 psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 250 ENDDO 251 ! 252 DO jk = nlev - 1, 1, -1 253 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 254 ENDDO 255 256 ENDDO 257 258 END DO 221 259 ! 222 223 za(nlev) = -rminus 224 zb(nlev) = 1. + rminus 225 zc(nlev) = 0. 226 227 228 ! solves tridiagonal system of linear equations 229 ! ----------------------------------------------- 230 DO jn = 1, nvar 231 232 zr (:,:) = psol(:,:,jn) 233 zbet = zb(1) 234 psol(:,1,jn) = zr(:,1) / zbet 235 ! 236 DO jk = 2, nlev 237 zgamm(jk) = zc(jk-1) / zbet 238 zbet = zb(jk) - za(jk) * zgamm(jk) 239 DO ji = 1, ndim 240 psol(ji,jk,jn) = ( zr(ji,jk) - za(jk) * psol(ji,jk-1,jn) ) / zbet 241 END DO 242 ENDDO 243 ! 244 DO jk = nlev - 1, 1, -1 245 DO ji = 1,ndim 246 psol(ji,jk,jn) = psol(ji,jk,jn) - zgamm(jk+1) * psol(ji,jk+1,jn) 247 END DO 248 ENDDO 249 250 ENDDO 260 IF( ln_timing ) CALL timing_stop('sed_mat_btb') 251 261 252 262 253 263 END SUBROUTINE sed_mat_btb 254 264 255 256 #else257 !!======================================================================258 !! MODULE sedmat : Dummy module259 !!======================================================================260 !! $Id$261 CONTAINS262 SUBROUTINE sed_mat ! Empty routine263 END SUBROUTINE sed_mat264 !!======================================================================265 #endif266 267 265 END MODULE sedmat
Note: See TracChangeset
for help on using the changeset viewer.
