- Timestamp:
- 2011-09-28T10:18:52+02:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r2715 r2872 20 20 USE in_out_manager ! I/O manager 21 21 USE lib_mpp ! MPP library 22 USE zdfbfr ! bottom friction setup 22 23 23 24 IMPLICIT NONE … … 61 62 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 62 63 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 INTEGER :: ikbum1, ikbvm1 ! local variable 66 REAL(wp) :: z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 67 68 !! * Local variables for implicit bottom friction. H. Liu 69 REAL(wp) :: zbfru, zbfrv 70 REAL(wp) :: bfr_imp = 1._wp 71 !!---------------------------------------------------------------------- 65 72 !!---------------------------------------------------------------------- 66 73 … … 73 80 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 74 81 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 82 IF(.NOT.ln_bfrimp) bfr_imp = 0._wp 75 83 ENDIF 76 84 … … 80 88 81 89 ! 1. Vertical diffusion on u 90 91 ! Vertical diffusion on u&v 82 92 ! --------------------------- 83 93 ! Matrix and second member construction 84 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 85 ! non zero value at the ocean bottom depending on the bottom friction 86 ! used but the bottom velocities have already been updated with the bottom 87 ! friction velocity in dyn_bfr using values from the previous timestep. There 88 ! is no need to include these in the implicit calculation. 89 ! 90 DO jk = 1, jpkm1 ! Matrix 91 DO jj = 2, jpjm1 92 DO ji = fs_2, fs_jpim1 ! vector opt. 94 !! bottom boundary condition: both zwi and zws must be masked as avmu can take 95 !! non zero value at the ocean bottom depending on the bottom friction 96 !! used but the bottom velocities have already been updated with the bottom 97 !! friction velocity in dyn_bfr using values from the previous timestep. There 98 !! is no need to include these in the implicit calculation. 99 100 ! The code has been modified here to implicitly implement bottom 101 ! friction: u(v)mask is not necessary here anymore. 102 ! H. Liu, April 2010. 103 104 ! 1. Vertical diffusion on u 105 DO jj = 2, jpjm1 106 DO ji = fs_2, fs_jpim1 ! vector opt. 107 ikbum1 = mbku(ji,jj) 108 ! Apply stability criteria on absolute value : Min abs(bfr) => Max (bfr) 109 ! zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv ) 110 zbfru = bfrua(ji,jj) 111 112 DO jk = 1, ikbum1 93 113 zcoef = - p2dt / fse3u(ji,jj,jk) 94 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 95 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 96 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 97 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 98 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 99 END DO 100 END DO 101 END DO 102 DO jj = 2, jpjm1 ! Surface boudary conditions 103 DO ji = fs_2, fs_jpim1 ! vector opt. 114 zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk ) / fse3uw(ji,jj,jk ) 115 zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 116 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 117 END DO 118 119 ! Surface boudary conditions 104 120 zwi(ji,jj,1) = 0._wp 105 121 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 106 END DO 107 END DO 122 123 ! Bottom boudary conditions ! H. Liu, May, 2010 124 ! !commented out to be consisent with v3.2, h.liu 125 ! z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp 126 z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp 127 zws(ji,jj,ikbum1) = 0._wp 128 zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 108 129 109 130 ! Matrix inversion starting from the first level … … 121 142 ! The solution (the after velocity) is in ua 122 143 !----------------------------------------------------------------------- 123 ! 124 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 144 145 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 146 DO jk = 2, ikbum1 127 147 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 128 148 END DO 129 END DO 130 END DO 131 ! 132 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 133 DO ji = fs_2, fs_jpim1 ! vector opt. 134 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 135 & / ( fse3u(ji,jj,1) * rau0 ) ) 136 END DO 137 END DO 138 DO jk = 2, jpkm1 139 DO jj = 2, jpjm1 140 DO ji = fs_2, fs_jpim1 ! vector opt. 149 150 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 151 z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 152 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 153 DO jk = 2, ikbum1 141 154 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side 142 155 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 143 156 END DO 144 END DO 145 END DO 146 ! 147 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 148 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 150 END DO 151 END DO 152 DO jk = jpk-2, 1, -1 153 DO jj = 2, jpjm1 154 DO ji = fs_2, fs_jpim1 ! vector opt. 155 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 157 158 159 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 160 ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 161 DO jk = ikbum1-1, 1, -1 162 ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 156 163 END DO 157 164 END DO … … 170 177 ! 2. Vertical diffusion on v 171 178 ! --------------------------- 172 ! Matrix and second member construction 173 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 174 ! non zero value at the ocean bottom depending on the bottom friction 175 ! used but the bottom velocities have already been updated with the bottom 176 ! friction velocity in dyn_bfr using values from the previous timestep. There 177 ! is no need to include these in the implicit calculation. 178 ! 179 DO jk = 1, jpkm1 ! Matrix 179 180 DO ji = fs_2, fs_jpim1 ! vector opt. 180 181 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 ikbvm1 = mbkv(ji,jj) 183 ! Apply stability criteria on absolute value : Min abs(bfr) => Max (bfr) 184 ! zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv ) 185 zbfrv = bfrva(ji,jj) 186 187 DO jk = 1, ikbvm1 182 188 zcoef = -p2dt / fse3v(ji,jj,jk) 183 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 184 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 185 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 186 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 187 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 188 END DO 189 END DO 190 END DO 191 DO jj = 2, jpjm1 ! Surface boudary conditions 192 DO ji = fs_2, fs_jpim1 ! vector opt. 189 zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk ) / fse3vw(ji,jj,jk ) 190 zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 191 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 192 END DO 193 194 ! Surface boudary conditions 193 195 zwi(ji,jj,1) = 0._wp 194 196 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 195 END DO 196 END DO 197 198 ! Bottom boudary conditions 199 ! Bottom boudary conditions ! H. Liu, May, 2010 200 ! !commented out to be consisent with v3.2, h.liu 201 ! z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp 202 z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp 203 zws(ji,jj,ikbvm1) = 0._wp 204 zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 197 205 198 206 ! Matrix inversion … … 206 214 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 207 215 ! 208 ! m is decomposed in the product of an upper and lower triangular matrix 216 ! m is decomposed in the product of an upper and lower triangular 217 ! matrix 209 218 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 210 219 ! The solution (after velocity) is in 2d array va 211 220 !----------------------------------------------------------------------- 212 ! 213 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 214 DO jj = 2, jpjm1 215 DO ji = fs_2, fs_jpim1 ! vector opt. 221 222 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 223 DO jk = 2, ikbvm1 216 224 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 217 225 END DO 218 END DO 219 END DO 220 ! 221 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 224 & / ( fse3v(ji,jj,1) * rau0 ) ) 225 END DO 226 END DO 227 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 226 227 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 228 z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 229 va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 230 DO jk = 2, ikbvm1 230 231 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side 231 232 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 232 233 END DO 233 END DO 234 END DO 235 ! 236 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 239 END DO 240 END DO 241 DO jk = jpk-2, 1, -1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 245 END DO 234 235 ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 236 va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 237 238 DO jk = ikbvm1-1, 1, -1 239 va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 240 END DO 241 246 242 END DO 247 243 END DO … … 258 254 IF( wrk_not_released(3, 3) ) CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 259 255 ! 256 260 257 END SUBROUTINE dyn_zdf_imp 261 258
Note: See TracChangeset
for help on using the changeset viewer.
