- Timestamp:
- 2020-05-14T21:46:00+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/LDF/ldfslp.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser
- Property svn:externals
-
old new 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@HEAD sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/LDF/ldfslp.F90
r10425 r12928 21 21 !!---------------------------------------------------------------------- 22 22 USE oce ! ocean dynamics and tracers 23 USE isf_oce ! ice shelf 23 24 USE dom_oce ! ocean space and time domain 24 25 ! USE ldfdyn ! lateral diffusion: eddy viscosity coef. … … 73 74 74 75 !! * Substitutions 75 # include " vectopt_loop_substitute.h90"76 # include "do_loop_substitute.h90" 76 77 !!---------------------------------------------------------------------- 77 78 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 81 82 CONTAINS 82 83 83 SUBROUTINE ldf_slp( kt, prd, pn2 )84 SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm ) 84 85 !!---------------------------------------------------------------------- 85 86 !! *** ROUTINE ldf_slp *** … … 107 108 !!---------------------------------------------------------------------- 108 109 INTEGER , INTENT(in) :: kt ! ocean time-step index 110 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 109 111 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density 110 112 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.) … … 134 136 zwz(:,:,:) = 0._wp 135 137 ! 136 DO jk = 1, jpk !== i- & j-gradient of density ==! 137 DO jj = 1, jpjm1 138 DO ji = 1, fs_jpim1 ! vector opt. 139 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 140 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 141 END DO 142 END DO 143 END DO 138 DO_3D_10_10( 1, jpk ) 139 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 140 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 141 END_3D 144 142 IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level 145 DO jj = 1, jpjm1 146 DO ji = 1, jpim1 147 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 148 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 149 END DO 150 END DO 143 DO_2D_10_10 144 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 145 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 146 END_2D 151 147 ENDIF 152 148 IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level 153 DO jj = 1, jpjm1 154 DO ji = 1, jpim1 155 IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 156 IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 157 END DO 158 END DO 149 DO_2D_10_10 150 IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 151 IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 152 END_2D 159 153 ENDIF 160 154 ! … … 171 165 ! 172 166 ! !== Slopes just below the mixed layer ==! 173 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml167 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm ) ! output: uslpml, vslpml, wslpiml, wslpjml 174 168 175 169 … … 178 172 ! 179 173 IF ( ln_isfcav ) THEN 180 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) & 183 & - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) ) 184 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) & 185 & - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) ) 186 END DO 187 END DO 174 DO_2D_00_00 175 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) & 176 & - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) ) 177 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) & 178 & - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) ) 179 END_2D 188 180 ELSE 189 DO jj = 2, jpjm1 190 DO ji = fs_2, fs_jpim1 ! vector opt. 191 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 192 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 193 END DO 194 END DO 181 DO_2D_00_00 182 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 183 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 184 END_2D 195 185 END IF 196 186 197 DO jk = 2, jpkm1 !* Slopes at u and v points 198 DO jj = 2, jpjm1 199 DO ji = fs_2, fs_jpim1 ! vector opt. 200 ! ! horizontal and vertical density gradient at u- and v-points 201 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 202 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 203 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 204 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 205 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 206 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 207 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau ) ) 208 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav ) ) 209 ! ! uslp and vslp output in zwz and zww, resp. 210 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 211 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 212 ! thickness of water column between surface and level k at u/v point 213 zdepu = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji+1,jj,jk) ) & 214 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u_n(ji,jj,miku(ji,jj)) ) 215 zdepv = 0.5_wp * ( ( gdept_n (ji,jj,jk) + gdept_n (ji,jj+1,jk) ) & 216 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v_n(ji,jj,mikv(ji,jj)) ) 217 ! 218 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 219 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 220 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 221 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 187 DO_3D_00_00( 2, jpkm1 ) 188 ! ! horizontal and vertical density gradient at u- and v-points 189 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 190 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 191 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 192 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 193 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 194 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 195 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 196 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 197 ! ! Fred Dupont: add a correction for bottom partial steps: 198 ! ! max slope = 1/2 * e3 / e1 199 IF (ln_zps .AND. jk==mbku(ji,jj)) & 200 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , - 2._wp * e1u(ji,jj) / e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 201 IF (ln_zps .AND. jk==mbkv(ji,jj)) & 202 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , - 2._wp * e2v(ji,jj) / e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 203 ! ! uslp and vslp output in zwz and zww, resp. 204 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 205 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 206 ! thickness of water column between surface and level k at u/v point 207 zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) ) & 208 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u(ji,jj,miku(ji,jj),Kmm) ) 209 zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) ) & 210 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v(ji,jj,mikv(ji,jj),Kmm) ) 211 ! 212 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 213 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 214 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 215 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 222 216 !!gm modif to suppress omlmask.... (as in Griffies case) 223 217 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 224 218 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 225 219 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 226 ! zci = 0.5 * ( gdept _n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )227 ! zcj = 0.5 * ( gdept _n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )220 ! zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 221 ! zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 228 222 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 229 223 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 230 224 !!gm end modif 231 END DO 232 END DO 233 END DO 225 END_3D 234 226 CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1., zww, 'V', -1. ) ! lateral boundary conditions 235 227 ! 236 228 ! !* horizontal Shapiro filter 237 229 DO jk = 2, jpkm1 238 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 239 DO ji = 2, jpim1 240 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 241 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 242 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 243 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 244 & + 4.* zwz(ji ,jj ,jk) ) 245 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 246 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 247 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 248 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 249 & + 4.* zww(ji,jj ,jk) ) 250 END DO 251 END DO 230 DO_2D_00_00 231 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 232 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 233 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 234 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 235 & + 4.* zwz(ji ,jj ,jk) ) 236 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 237 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 238 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 239 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 240 & + 4.* zww(ji,jj ,jk) ) 241 END_2D 252 242 DO jj = 3, jpj-2 ! other rows 253 DO ji = fs_2, fs_jpim1 ! vector opt.243 DO ji = 2, jpim1 ! vector opt. 254 244 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 255 245 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 265 255 END DO 266 256 ! !* decrease along coastal boundaries 267 DO jj = 2, jpjm1 268 DO ji = fs_2, fs_jpim1 ! vector opt. 269 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 270 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 271 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 272 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 273 END DO 274 END DO 257 DO_2D_00_00 258 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 259 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 260 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 261 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 262 END_2D 275 263 END DO 276 264 … … 279 267 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 280 268 ! 281 DO jk = 2, jpkm1 282 DO jj = 2, jpjm1 283 DO ji = fs_2, fs_jpim1 ! vector opt. 284 ! !* Local vertical density gradient evaluated from N^2 285 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 286 ! !* Slopes at w point 287 ! ! i- & j-gradient of density at w-points 288 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 289 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 290 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 291 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 292 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 293 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 294 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 295 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 296 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 297 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 298 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai ) ) 299 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj ) ) 300 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 301 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 302 zck = ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj) - gdepw_n(ji,jj,mikt(ji,jj)), 10._wp ) 303 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 304 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 269 DO_3D_00_00( 2, jpkm1 ) 270 ! !* Local vertical density gradient evaluated from N^2 271 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 272 ! !* Slopes at w point 273 ! ! i- & j-gradient of density at w-points 274 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 275 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 276 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 277 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 278 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 279 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 280 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 281 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 282 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 283 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 284 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai ) ) 285 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj ) ) 286 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 287 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 288 zck = ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) ) / MAX( hmlp(ji,jj) - gdepw(ji,jj,mikt(ji,jj),Kmm), 10._wp ) 289 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 290 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 305 291 306 292 !!gm modif to suppress omlmask.... (as in Griffies operator) … … 311 297 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 312 298 !!gm end modif 313 END DO 314 END DO 315 END DO 299 END_3D 316 300 CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1., zww, 'T', -1. ) ! lateral boundary conditions 317 301 ! 318 302 ! !* horizontal Shapiro filter 319 303 DO jk = 2, jpkm1 320 DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only 321 DO ji = 2, jpim1 322 zcofw = wmask(ji,jj,jk) * z1_16 323 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 324 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 325 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 326 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 327 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 328 329 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 330 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 331 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 332 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 333 & + 4.* zww(ji ,jj ,jk) ) * zcofw 334 END DO 335 END DO 304 DO_2D_00_00 305 zcofw = wmask(ji,jj,jk) * z1_16 306 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 307 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 308 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 309 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 310 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 311 312 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 313 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 314 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 315 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 316 & + 4.* zww(ji ,jj ,jk) ) * zcofw 317 END_2D 336 318 DO jj = 3, jpj-2 ! other rows 337 DO ji = fs_2, fs_jpim1 ! vector opt.319 DO ji = 2, jpim1 ! vector opt. 338 320 zcofw = wmask(ji,jj,jk) * z1_16 339 321 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & … … 351 333 END DO 352 334 ! !* decrease in vicinity of topography 353 DO jj = 2, jpjm1 354 DO ji = fs_2, fs_jpim1 ! vector opt. 355 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 356 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 357 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 358 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 359 END DO 360 END DO 335 DO_2D_00_00 336 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 337 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 338 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 339 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 340 END_2D 361 341 END DO 362 342 … … 365 345 CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1. , vslp , 'V', -1. , wslpi, 'W', -1., wslpj, 'W', -1. ) 366 346 367 IF( ln_ctl) THEN347 IF(sn_cfctl%l_prtctl) THEN 368 348 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 369 349 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) … … 375 355 376 356 377 SUBROUTINE ldf_slp_triad ( kt )357 SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm ) 378 358 !!---------------------------------------------------------------------- 379 359 !! *** ROUTINE ldf_slp_triad *** … … 390 370 !!---------------------------------------------------------------------- 391 371 INTEGER, INTENT( in ) :: kt ! ocean time-step index 372 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 392 373 !! 393 374 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 402 383 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 403 384 REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw 404 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalbet405 385 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 406 386 REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb ! for Griffies operator only … … 416 396 ! 417 397 ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln) 418 DO jk = 1, jpkm1 ! done each pair of triad 419 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 420 DO ji = 1, fs_jpim1 ! vector opt. 421 zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point 422 zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 423 zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point 424 zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) ) 425 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 426 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 427 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 428 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 429 END DO 430 END DO 431 END DO 398 DO_3D_10_10( 1, jpkm1 ) 399 zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! i-gradient of T & S at u-point 400 zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 401 zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! j-gradient of T & S at v-point 402 zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 403 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 404 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 405 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 406 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 407 END_3D 432 408 ! 433 409 IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom 434 DO jj = 1, jpjm1 435 DO ji = 1, jpim1 436 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 437 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 438 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 439 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 440 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 441 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 442 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 443 END DO 444 END DO 410 DO_2D_10_10 411 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 412 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 413 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 414 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 415 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 416 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 417 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 418 END_2D 445 419 ENDIF 446 420 ! … … 448 422 449 423 DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==! 450 DO jk = 1, jpkm1 ! done each pair of triad 451 DO jj = 1, jpj ! NB: not masked ==> a minimum value is set 452 DO ji = 1, jpi ! vector opt. 453 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 454 zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) ) 455 zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) ) 456 ELSE 457 zdkt = 0._wp ! 1st level gradient set to zero 458 zdks = 0._wp 459 ENDIF 460 zdzrho_raw = ( - rab_b(ji,jj,jk+kp,jp_tem) * zdkt & 461 & + rab_b(ji,jj,jk+kp,jp_sal) * zdks & 462 & ) / e3w_n(ji,jj,jk+kp) 463 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 464 END DO 465 END DO 466 END DO 424 DO_3D_11_11( 1, jpkm1 ) 425 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 426 zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) ) 427 zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) ) 428 ELSE 429 zdkt = 0._wp ! 1st level gradient set to zero 430 zdks = 0._wp 431 ENDIF 432 zdzrho_raw = ( - rab_b(ji,jj,jk ,jp_tem) * zdkt & 433 & + rab_b(ji,jj,jk ,jp_sal) * zdks & 434 & ) / e3w(ji,jj,jk+kp,Kmm) 435 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 436 END_3D 467 437 END DO 468 438 ! 469 DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==! 470 DO ji = 1, jpi 471 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 472 z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk) 473 END DO 474 END DO 439 DO_2D_11_11 440 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 441 z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm) 442 END_2D 475 443 ! 476 444 ! !== intialisations to zero ==! … … 489 457 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 490 458 DO kp = 0, 1 ! with only the slope-max limit and MASKED 491 DO jj = 1, jpjm1 492 DO ji = 1, fs_jpim1 493 ip = jl ; jp = jl 494 ! 495 jk = nmln(ji+ip,jj) + 1 496 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 497 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 498 ELSE 499 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 500 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 501 & - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 502 ze3_e1 = e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 503 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 504 ENDIF 505 ! 506 jk = nmln(ji,jj+jp) + 1 507 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 508 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 509 ELSE 510 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 511 & - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 512 ze3_e2 = e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj) 513 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 514 ENDIF 515 END DO 516 END DO 459 DO_2D_10_10 460 ip = jl ; jp = jl 461 ! 462 jk = nmln(ji+ip,jj) + 1 463 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 464 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 465 ELSE 466 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 467 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 468 & - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 469 ze3_e1 = e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj) 470 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 471 ENDIF 472 ! 473 jk = nmln(ji,jj+jp) + 1 474 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 475 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 476 ELSE 477 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 478 & - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 479 ze3_e2 = e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj) 480 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 481 ENDIF 482 END_2D 517 483 END DO 518 484 END DO … … 528 494 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 529 495 znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0 530 DO jj = 1, jpjm1 531 DO ji = 1, fs_jpim1 ! vector opt. 532 ! 533 ! Calculate slope relative to geopotentials used for GM skew fluxes 534 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 535 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 536 ! masked by umask taken at the level of dz(rho) 537 ! 538 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 539 ! 540 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 541 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 542 ! 543 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 544 zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) 545 ztj_coord = znot_thru_surface * ( gdept_n(ji ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 546 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 547 ztj_g_raw = ztj_raw - ztj_coord 548 ! additional limit required in bilaplacian case 549 ze3_e1 = e3w_n(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 550 ze3_e2 = e3w_n(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 551 ! NB: hard coded factor 5 (can be a namelist parameter...) 552 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 553 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 554 ! 555 ! Below ML use limited zti_g as is & mask 556 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 557 ! 558 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 559 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 560 ! ! otherwise zfact=0 561 zti_g_lim = ( zfacti * zti_g_lim & 562 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 563 & * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 564 ztj_g_lim = ( zfactj * ztj_g_lim & 565 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 566 & * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 567 ! 568 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 569 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 570 ! 571 ! Get coefficients of isoneutral diffusion tensor 572 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 573 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 574 ! i.e. 33 term = (real slope* 31, 13 terms) 575 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 576 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 577 ! 578 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 579 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 580 ! 581 IF( ln_triad_iso ) THEN 582 zti_raw = zti_lim*zti_lim / zti_raw 583 ztj_raw = ztj_lim*ztj_lim / ztj_raw 584 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 585 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 586 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 587 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 588 ENDIF 589 ! ! switching triad scheme 590 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 591 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 592 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 593 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 594 ! 595 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 596 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 597 ! 598 zbu = e1e2u(ji ,jj ) * e3u_n(ji ,jj ,jk ) 599 zbv = e1e2v(ji ,jj ) * e3v_n(ji ,jj ,jk ) 600 zbti = e1e2t(ji+ip,jj ) * e3w_n(ji+ip,jj ,jk+kp) 601 zbtj = e1e2t(ji ,jj+jp) * e3w_n(ji ,jj+jp,jk+kp) 602 ! 603 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 604 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 605 END DO 606 END DO 496 DO_2D_10_10 497 ! 498 ! Calculate slope relative to geopotentials used for GM skew fluxes 499 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 500 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 501 ! masked by umask taken at the level of dz(rho) 502 ! 503 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 504 ! 505 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 506 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 507 ! 508 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 509 zti_coord = znot_thru_surface * ( gdept(ji+1,jj ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 510 ztj_coord = znot_thru_surface * ( gdept(ji ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) ! unmasked 511 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 512 ztj_g_raw = ztj_raw - ztj_coord 513 ! additional limit required in bilaplacian case 514 ze3_e1 = e3w(ji+ip,jj ,jk+kp,Kmm) * r1_e1u(ji,jj) 515 ze3_e2 = e3w(ji ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj) 516 ! NB: hard coded factor 5 (can be a namelist parameter...) 517 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 518 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 519 ! 520 ! Below ML use limited zti_g as is & mask 521 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 522 ! 523 zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1 524 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 525 ! ! otherwise zfact=0 526 zti_g_lim = ( zfacti * zti_g_lim & 527 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 528 & * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 529 ztj_g_lim = ( zfactj * ztj_g_lim & 530 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 531 & * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 532 ! 533 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 534 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 535 ! 536 ! Get coefficients of isoneutral diffusion tensor 537 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 538 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 539 ! i.e. 33 term = (real slope* 31, 13 terms) 540 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 541 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 542 ! 543 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 544 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 545 ! 546 IF( ln_triad_iso ) THEN 547 zti_raw = zti_lim*zti_lim / zti_raw 548 ztj_raw = ztj_lim*ztj_lim / ztj_raw 549 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 550 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 551 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 552 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 553 ENDIF 554 ! ! switching triad scheme 555 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 556 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 557 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 558 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 559 ! 560 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 561 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 562 ! 563 zbu = e1e2u(ji ,jj ) * e3u(ji ,jj ,jk ,Kmm) 564 zbv = e1e2v(ji ,jj ) * e3v(ji ,jj ,jk ,Kmm) 565 zbti = e1e2t(ji+ip,jj ) * e3w(ji+ip,jj ,jk+kp,Kmm) 566 zbtj = e1e2t(ji ,jj+jp) * e3w(ji ,jj+jp,jk+kp,Kmm) 567 ! 568 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 569 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 570 END_2D 607 571 END DO 608 572 END DO … … 618 582 619 583 620 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )584 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm ) 621 585 !!---------------------------------------------------------------------- 622 586 !! *** ROUTINE ldf_slp_mxl *** … … 638 602 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 639 603 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 604 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 640 605 !! 641 606 INTEGER :: ji , jj , jk ! dummy loop indices … … 658 623 ! 659 624 ! !== surface mixed layer mask ! 660 DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise 661 DO jj = 1, jpj 662 DO ji = 1, jpi 663 ik = nmln(ji,jj) - 1 664 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 665 ELSE ; omlmask(ji,jj,jk) = 0._wp 666 ENDIF 667 END DO 668 END DO 669 END DO 625 DO_3D_11_11( 1, jpk ) 626 ik = nmln(ji,jj) - 1 627 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 628 ELSE ; omlmask(ji,jj,jk) = 0._wp 629 ENDIF 630 END_3D 670 631 671 632 … … 680 641 !----------------------------------------------------------------------- 681 642 ! 682 DO jj = 2, jpjm1 683 DO ji = 2, jpim1 684 ! !== Slope at u- & v-points just below the Mixed Layer ==! 685 ! 686 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 687 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 688 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 689 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 690 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 691 ! !- horizontal density gradient at u- & v-points 692 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 693 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 694 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 695 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 696 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau ) ) 697 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav ) ) 698 ! !- Slope at u- & v-points (uslpml, vslpml) 699 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 700 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 701 ! 702 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 703 ! 704 ik = MIN( nmln(ji,jj) + 1, jpk ) 705 ikm1 = MAX( 1, ik-1 ) 706 ! !- vertical density gradient for w-slope (from N^2) 707 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 708 ! !- horizontal density i- & j-gradient at w-points 709 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 710 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 711 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 712 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 713 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 714 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 715 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 716 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 717 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 718 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 719 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai ) ) 720 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj ) ) 721 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 722 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 723 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 724 END DO 725 END DO 643 DO_2D_00_00 644 ! !== Slope at u- & v-points just below the Mixed Layer ==! 645 ! 646 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 647 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 648 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 649 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 650 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 651 ! !- horizontal density gradient at u- & v-points 652 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 653 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 654 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 655 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 656 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau ) ) 657 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav ) ) 658 ! !- Slope at u- & v-points (uslpml, vslpml) 659 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 660 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 661 ! 662 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 663 ! 664 ik = MIN( nmln(ji,jj) + 1, jpk ) 665 ikm1 = MAX( 1, ik-1 ) 666 ! !- vertical density gradient for w-slope (from N^2) 667 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 668 ! !- horizontal density i- & j-gradient at w-points 669 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 670 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 671 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 672 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 673 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 674 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 675 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 676 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 677 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 678 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 679 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai ) ) 680 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj ) ) 681 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 682 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 683 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 684 END_2D 726 685 !!gm this lbc_lnk should be useless.... 727 686 CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1. , vslpml , 'V', -1. , wslpiml, 'W', -1. , wslpjml, 'W', -1. ) … … 785 744 ! DO jk = 1, jpk 786 745 ! DO jj = 2, jpjm1 787 ! DO ji = fs_2, fs_jpim1 ! vector opt.788 ! uslp (ji,jj,jk) = - ( gdept _n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)789 ! vslp (ji,jj,jk) = - ( gdept _n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)790 ! wslpi(ji,jj,jk) = - ( gdepw _n(ji+1,jj,jk) - gdepw_n(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5791 ! wslpj(ji,jj,jk) = - ( gdepw _n(ji,jj+1,jk) - gdepw_n(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5746 ! DO ji = 2, jpim1 ! vector opt. 747 ! uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 748 ! vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 749 ! wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kmm) - gdepw(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5 750 ! wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kmm) - gdepw(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5 792 751 ! END DO 793 752 ! END DO
Note: See TracChangeset
for help on using the changeset viewer.
