Changeset 3259 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM
- Timestamp:
- 2012-01-13T16:24:15+01:00 (12 years ago)
- Location:
- branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC
- Files:
-
- 3 edited
-
DYN/dynbfr.F90 (modified) (2 diffs)
-
DYN/dynzdf_imp.F90 (modified) (8 diffs)
-
ZDF/zdfbfr.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r3231 r3259 5 5 !!============================================================================== 6 6 !! History : 3.2 ! 2008-11 (A. C. Coward) Original code 7 !! 3.4 ! 2011-09 (H. Liu) Make it consistent with semi-implicit 8 !! Bottom friction (ln_bfrimp = .true.) 7 9 !!---------------------------------------------------------------------- 8 10 … … 57 59 IF( .NOT.ln_bfrimp) THEN ! only for explicit bottom friction form 58 60 ! implicit bfr is implemented in dynzdf_imp 59 ! H. Liu, Sept. 201160 61 61 62 zm1_2dt = - 1._wp / ( 2._wp * rdt ) -
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r3231 r3259 8 8 !! NEMO 0.5 ! 2002-08 (G. Madec) F90: Free form and module 9 9 !! 3.3 ! 2010-04 (M. Leclair, G. Madec) Forcing averaged over 2 time steps 10 !! 3.4 ! 2012-01 (H. Liu) Semi-implicit bottom friction 10 11 !!---------------------------------------------------------------------- 11 12 … … 20 21 USE in_out_manager ! I/O manager 21 22 USE lib_mpp ! MPP library 22 USE zdfbfr ! bottom friction setup23 USE zdfbfr ! Bottom friction setup 23 24 USE wrk_nemo ! Memory Allocation 24 25 USE timing ! Timing … … 60 61 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 61 62 !! 62 INTEGER :: ji, jj, jk ! dummy loop indices 63 INTEGER :: ikbum1, ikbvm1 ! local variable 64 REAL(wp) :: z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 65 66 !! * Local variables for implicit bottom friction. H. Liu 67 REAL(wp) :: zbfru, zbfrv 68 REAL(wp) :: zbfr_imp = 0._wp ! toggle (SAVE'd by assignment) 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 INTEGER :: ikbu, ikbv ! local integers 65 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 66 !!---------------------------------------------------------------------- 67 69 68 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwd, zws 70 69 !!---------------------------------------------------------------------- … … 78 77 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 79 78 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 80 IF(ln_bfrimp) zbfr_imp = 1._wp81 79 ENDIF 82 80 … … 85 83 z1_p2dt = 1._wp / p2dt ! inverse of the timestep 86 84 87 ! 1. Vertical diffusion on u 88 89 ! Vertical diffusion on u&v 85 ! 1. Apply semi-implicit bottom friction 86 ! -------------------------------------- 87 ! Only needed for semi-implicit bottom friction setup. The explicit 88 ! bottom friction has been included in "u(v)a" which act as the R.H.S 89 ! column vector of the tri-diagonal matrix equation 90 ! 91 IF( ln_bfrimp ) THEN 92 # if defined key_vectopt_loop 93 DO jj = 1, 1 94 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 95 # else 96 DO jj = 2, jpjm1 97 DO ji = 2, jpim1 98 # endif 99 ikbu = mbku(ji,jj) ! ocean bottom level at u- and v-points 100 ikbv = mbkv(ji,jj) ! (deepest ocean u- and v-points) 101 avmu(ji,jj,ikbu+1) = -bfrua(ji,jj) * fse3uw(ji,jj,ikbu+1) 102 avmv(ji,jj,ikbv+1) = -bfrva(ji,jj) * fse3vw(ji,jj,ikbv+1) 103 END DO 104 END DO 105 ENDIF 106 107 ! 2. Vertical diffusion on u 90 108 ! --------------------------- 91 109 ! Matrix and second member construction 92 !! bottom boundary condition: both zwi and zws must be masked as avmu can take 93 !! non zero value at the ocean bottom depending on the bottom friction 94 !! used but the bottom velocities have already been updated with the bottom 95 !! friction velocity in dyn_bfr using values from the previous timestep. There 96 !! is no need to include these in the implicit calculation. 97 98 ! The code has been modified here to implicitly implement bottom 99 ! friction: u(v)mask is not necessary here anymore. 100 ! H. Liu, April 2010. 101 102 ! 1. Vertical diffusion on u 103 DO jj = 2, jpjm1 104 DO ji = fs_2, fs_jpim1 ! vector opt. 105 ikbum1 = mbku(ji,jj) 106 zbfru = bfrua(ji,jj) 107 108 DO jk = 1, ikbum1 110 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 111 ! non zero value at the ocean bottom depending on the bottom friction used. 112 ! 113 DO jk = 1, jpkm1 ! Matrix 114 DO jj = 2, jpjm1 115 DO ji = fs_2, fs_jpim1 ! vector opt. 109 116 zcoef = - p2dt / fse3u(ji,jj,jk) 110 zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk ) / fse3uw(ji,jj,jk ) 111 zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 112 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 113 END DO 114 115 ! Surface boundary conditions 117 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 118 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 119 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 120 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 121 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 122 END DO 123 END DO 124 END DO 125 DO jj = 2, jpjm1 ! Surface boudary conditions 126 DO ji = fs_2, fs_jpim1 ! vector opt. 116 127 zwi(ji,jj,1) = 0._wp 117 128 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 118 119 ! Bottom boundary conditions ! H. Liu, May, 2010 120 ! !commented out to be consistent with v3.2, h.liu 121 ! z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0_wp * zbfr_imp 122 z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * zbfr_imp 123 zws(ji,jj,ikbum1) = 0._wp 124 zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 129 END DO 130 END DO 125 131 126 132 ! Matrix inversion starting from the first level … … 138 144 ! The solution (the after velocity) is in ua 139 145 !----------------------------------------------------------------------- 140 141 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 142 DO jk = 2, ikbum1 146 ! 147 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 148 DO jj = 2, jpjm1 149 DO ji = fs_2, fs_jpim1 ! vector opt. 143 150 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 144 151 END DO 145 146 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 147 z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 148 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 149 DO jk = 2, ikbum1 152 END DO 153 END DO 154 ! 155 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 156 DO ji = fs_2, fs_jpim1 ! vector opt. 157 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 158 & / ( fse3u(ji,jj,1) * rau0 ) ) 159 END DO 160 END DO 161 DO jk = 2, jpkm1 162 DO jj = 2, jpjm1 163 DO ji = fs_2, fs_jpim1 ! vector opt. 150 164 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side 151 165 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 152 166 END DO 153 154 155 ! third recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 156 ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 157 DO jk = ikbum1-1, 1, -1 158 ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 167 END DO 168 END DO 169 ! 170 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 171 DO ji = fs_2, fs_jpim1 ! vector opt. 172 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 173 END DO 174 END DO 175 DO jk = jpk-2, 1, -1 176 DO jj = 2, jpjm1 177 DO ji = fs_2, fs_jpim1 ! vector opt. 178 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 159 179 END DO 160 180 END DO … … 171 191 172 192 173 ! 2. Vertical diffusion on v193 ! 3. Vertical diffusion on v 174 194 ! --------------------------- 175 176 DO ji = fs_2, fs_jpim1 ! vector opt.177 DO jj = 2, jpjm1178 ikbvm1 = mbkv(ji,jj)179 zbfrv = bfrva(ji,jj)180 181 DO j k = 1, ikbvm1195 ! Matrix and second member construction 196 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 197 ! non zero value at the ocean bottom depending on the bottom friction used 198 ! 199 DO jk = 1, jpkm1 ! Matrix 200 DO jj = 2, jpjm1 201 DO ji = fs_2, fs_jpim1 ! vector opt. 182 202 zcoef = -p2dt / fse3v(ji,jj,jk) 183 zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk ) / fse3vw(ji,jj,jk ) 184 zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 185 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 186 END DO 187 188 ! Surface boundary conditions 203 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 204 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 205 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 206 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 207 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 208 END DO 209 END DO 210 END DO 211 DO jj = 2, jpjm1 ! Surface boudary conditions 212 DO ji = fs_2, fs_jpim1 ! vector opt. 189 213 zwi(ji,jj,1) = 0._wp 190 214 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 191 192 ! Bottom boundary conditions ! H. Liu, May, 2010 193 ! !commented out to be consistent with v3.2, h.liu 194 ! z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0_wp * zbfr_imp 195 z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * zbfr_imp 196 zws(ji,jj,ikbvm1) = 0._wp 197 zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 215 END DO 216 END DO 198 217 199 218 ! Matrix inversion … … 207 226 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 208 227 ! 209 ! m is decomposed in the product of an upper and lower triangular 210 ! matrix 228 ! m is decomposed in the product of an upper and lower triangular matrix 211 229 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 212 230 ! The solution (after velocity) is in 2d array va 213 231 !----------------------------------------------------------------------- 214 215 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 216 DO jk = 2, ikbvm1 232 ! 233 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 234 DO jj = 2, jpjm1 235 DO ji = fs_2, fs_jpim1 ! vector opt. 217 236 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 218 237 END DO 219 220 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 221 z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 222 va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 223 DO jk = 2, ikbvm1 238 END DO 239 END DO 240 ! 241 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 244 & / ( fse3v(ji,jj,1) * rau0 ) ) 245 END DO 246 END DO 247 DO jk = 2, jpkm1 248 DO jj = 2, jpjm1 249 DO ji = fs_2, fs_jpim1 ! vector opt. 224 250 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side 225 251 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 226 252 END DO 227 228 ! third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 229 va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 230 231 DO jk = ikbvm1-1, 1, -1 232 va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 233 END DO 234 253 END DO 254 END DO 255 ! 256 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 257 DO ji = fs_2, fs_jpim1 ! vector opt. 258 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 259 END DO 260 END DO 261 DO jk = jpk-2, 1, -1 262 DO jj = 2, jpjm1 263 DO ji = fs_2, fs_jpim1 ! vector opt. 264 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 265 END DO 235 266 END DO 236 267 END DO -
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r3229 r3259 123 123 & tab2d_2=bfrva, clinfo2= ' v: ', mask2=vmask,ovlap=1 ) 124 124 ENDIF 125 125 126 ! 126 127 IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr') … … 170 171 WRITE(numout,*) 'Implicit bottom friction can only be used when ln_zdfexp=.false.' 171 172 WRITE(numout,*) ' but you set: ln_bfrimp=.true. and ln_zdfexp=.true.!!!!' 172 WRITE(ctmp1,*) ' bad ln_bfrimp value = .true.'173 WRITE(ctmp1,*) ' set either "ln_zdfexp = .false" or "ln_bfrimp = .false." 173 174 CALL ctl_stop( ctmp1 ) 174 175 END IF … … 272 273 CALL mpp_max( zmaxbfr ) 273 274 ENDIF 275 IF( .NOT.ln_bfrimp) THEN 274 276 IF( lwp .AND. ictu + ictv > 0 ) THEN 275 277 WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points ' … … 278 280 WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 279 281 ENDIF 282 ENDIF 280 283 ! 281 284 IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr_init')
Note: See TracChangeset
for help on using the changeset viewer.
