- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/LDF/ldfslp.F90 (modified) (23 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/LDF/ldfslp.F90
r10425 r13463 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" 77 # include "domzgr_substitute.h90" 76 78 !!---------------------------------------------------------------------- 77 79 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 81 83 CONTAINS 82 84 83 SUBROUTINE ldf_slp( kt, prd, pn2 )85 SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm ) 84 86 !!---------------------------------------------------------------------- 85 87 !! *** ROUTINE ldf_slp *** … … 107 109 !!---------------------------------------------------------------------- 108 110 INTEGER , INTENT(in) :: kt ! ocean time-step index 111 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 109 112 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density 110 113 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.) … … 134 137 zwz(:,:,:) = 0._wp 135 138 ! 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 139 DO_3D( 1, 0, 1, 0, 1, jpk ) 140 zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) ) 141 zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) ) 142 END_3D 144 143 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 144 DO_2D( 1, 0, 1, 0 ) 145 zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj) 146 zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj) 147 END_2D 151 148 ENDIF 152 149 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 150 DO_2D( 1, 0, 1, 0 ) 151 IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj) 152 IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj) 153 END_2D 159 154 ENDIF 160 155 ! … … 171 166 ! 172 167 ! !== Slopes just below the mixed layer ==! 173 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml168 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm ) ! output: uslpml, vslpml, wslpiml, wslpjml 174 169 175 170 … … 178 173 ! 179 174 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 175 DO_2D( 0, 0, 0, 0 ) 176 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) & 177 & - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) ) 178 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) & 179 & - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) ) 180 END_2D 188 181 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 182 DO_2D( 0, 0, 0, 0 ) 183 zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp) 184 zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp) 185 END_2D 195 186 END IF 196 187 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) 188 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 189 ! ! horizontal and vertical density gradient at u- and v-points 190 zau = zgru(ji,jj,jk) * r1_e1u(ji,jj) 191 zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj) 192 zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) ) 193 zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) ) 194 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 195 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 196 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 197 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 198 ! ! Fred Dupont: add a correction for bottom partial steps: 199 ! ! max slope = 1/2 * e3 / e1 200 IF (ln_zps .AND. jk==mbku(ji,jj)) & 201 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , & 202 & - 2._wp * e1u(ji,jj) / e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 203 IF (ln_zps .AND. jk==mbkv(ji,jj)) & 204 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , & 205 & - 2._wp * e2v(ji,jj) / e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 206 ! ! uslp and vslp output in zwz and zww, resp. 207 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 208 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 209 ! thickness of water column between surface and level k at u/v point 210 zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) ) & 211 & - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) & 212 & - e3u(ji,jj,miku(ji,jj),Kmm) ) 213 zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) ) & 214 & - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) & 215 & - e3v(ji,jj,mikv(ji,jj),Kmm) ) 216 ! 217 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & 218 & + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk) 219 zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) & 220 & + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk) 222 221 !!gm modif to suppress omlmask.... (as in Griffies case) 223 222 ! ! ! jk must be >= ML level for zf=1. otherwise zf=0. 224 223 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 225 224 ! 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. ) )225 ! zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 226 ! zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 228 227 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 229 228 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) 230 229 !!gm end modif 231 END DO 232 END DO 233 END DO 234 CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1., zww, 'V', -1. ) ! lateral boundary conditions 230 END_3D 231 CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1.0_wp, zww, 'V', -1.0_wp ) ! lateral boundary conditions 235 232 ! 236 233 ! !* horizontal Shapiro filter 237 234 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 235 DO_2D( 0, 0, 0, 0 ) 236 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 237 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 238 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 239 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 240 & + 4.* zwz(ji ,jj ,jk) ) 241 vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 242 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 243 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 244 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 245 & + 4.* zww(ji,jj ,jk) ) 246 END_2D 252 247 DO jj = 3, jpj-2 ! other rows 253 DO ji = fs_2, fs_jpim1 ! vector opt.248 DO ji = 2, jpim1 ! vector opt. 254 249 uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 255 250 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & … … 265 260 END DO 266 261 ! !* 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 262 DO_2D( 0, 0, 0, 0 ) 263 uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp & 264 & * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp 265 vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp & 266 & * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp 267 END_2D 275 268 END DO 276 269 … … 279 272 ! =========================== | wslpj = mij( d/dj( prd ) / d/dz( prd ) 280 273 ! 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) 274 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 275 ! !* Local vertical density gradient evaluated from N^2 276 zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) 277 ! !* Slopes at w point 278 ! ! i- & j-gradient of density at w-points 279 zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) & 280 & + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj) 281 zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) & 282 & + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj) 283 zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) & 284 & + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk) 285 zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) & 286 & + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk) 287 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 288 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 289 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai ) ) 290 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj ) ) 291 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 292 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 293 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 ) 294 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 295 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) 305 296 306 297 !!gm modif to suppress omlmask.... (as in Griffies operator) 307 298 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 308 299 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 309 ! zck = gdepw(ji,jj,jk ) / MAX( hmlp(ji,jj), 10. )300 ! zck = gdepw(ji,jj,jk,Kmm) / MAX( hmlp(ji,jj), 10. ) 310 301 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 311 302 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) 312 303 !!gm end modif 313 END DO 314 END DO 315 END DO 316 CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1., zww, 'T', -1. ) ! lateral boundary conditions 304 END_3D 305 CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1.0_wp, zww, 'T', -1.0_wp ) ! lateral boundary conditions 317 306 ! 318 307 ! !* horizontal Shapiro filter 319 308 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 309 DO_2D( 0, 0, 0, 0 ) 310 zcofw = wmask(ji,jj,jk) * z1_16 311 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & 312 & + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) & 313 & + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) & 314 & + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) & 315 & + 4.* zwz(ji ,jj ,jk) ) * zcofw 316 317 wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) & 318 & + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) & 319 & + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) & 320 & + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) & 321 & + 4.* zww(ji ,jj ,jk) ) * zcofw 322 END_2D 336 323 DO jj = 3, jpj-2 ! other rows 337 DO ji = fs_2, fs_jpim1 ! vector opt.324 DO ji = 2, jpim1 ! vector opt. 338 325 zcofw = wmask(ji,jj,jk) * z1_16 339 326 wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) & … … 351 338 END DO 352 339 ! !* 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 340 DO_2D( 0, 0, 0, 0 ) 341 zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) & 342 & * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25 343 wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck 344 wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck 345 END_2D 361 346 END DO 362 347 363 348 ! IV. Lateral boundary conditions 364 349 ! =============================== 365 CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1. , vslp , 'V', -1. , wslpi, 'W', -1., wslpj, 'W', -1.)366 367 IF( ln_ctl) THEN350 CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1.0_wp , vslp , 'V', -1.0_wp , wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp ) 351 352 IF(sn_cfctl%l_prtctl) THEN 368 353 CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk) 369 354 CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk) … … 375 360 376 361 377 SUBROUTINE ldf_slp_triad ( kt )362 SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm ) 378 363 !!---------------------------------------------------------------------- 379 364 !! *** ROUTINE ldf_slp_triad *** … … 390 375 !!---------------------------------------------------------------------- 391 376 INTEGER, INTENT( in ) :: kt ! ocean time-step index 377 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 392 378 !! 393 379 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 402 388 REAL(wp) :: zbeta0, ze3_e1, ze3_e2 403 389 REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw 404 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zalbet405 390 REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients 406 391 REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb ! for Griffies operator only … … 416 401 ! 417 402 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 403 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 404 zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! i-gradient of T & S at u-point 405 zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 406 zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! j-gradient of T & S at v-point 407 zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 408 zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj) 409 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) 410 zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 411 zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 412 END_3D 432 413 ! 433 414 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 415 DO_2D( 1, 0, 1, 0 ) 416 iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points) 417 zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature 418 zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity 419 zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj) 420 zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj) 421 zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign 422 zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw ) 423 END_2D 445 424 ENDIF 446 425 ! … … 448 427 449 428 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 429 DO_3D( 1, 1, 1, 1, 1, jpkm1 ) 430 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 431 zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) ) 432 zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) ) 433 ELSE 434 zdkt = 0._wp ! 1st level gradient set to zero 435 zdks = 0._wp 436 ENDIF 437 zdzrho_raw = ( - rab_b(ji,jj,jk ,jp_tem) * zdkt & 438 & + rab_b(ji,jj,jk ,jp_sal) * zdks & 439 & ) / e3w(ji,jj,jk+kp,Kmm) 440 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 441 END_3D 467 442 END DO 468 443 ! 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 444 DO_2D( 1, 1, 1, 1 ) 445 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 446 z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm) 447 END_2D 475 448 ! 476 449 ! !== intialisations to zero ==! … … 489 462 DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base 490 463 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 464 DO_2D( 1, 0, 1, 0 ) 465 ip = jl ; jp = jl 466 ! 467 jk = nmln(ji+ip,jj) + 1 468 IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom 469 zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp 470 ELSE 471 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 472 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 473 & - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 474 ze3_e1 = e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj) 475 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 476 ENDIF 477 ! 478 jk = nmln(ji,jj+jp) + 1 479 IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom 480 ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp 481 ELSE 482 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 483 & - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 484 ze3_e2 = e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj) 485 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 486 ENDIF 487 END_2D 517 488 END DO 518 489 END DO … … 528 499 ! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface 529 500 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 501 DO_2D( 1, 0, 1, 0 ) 502 ! 503 ! Calculate slope relative to geopotentials used for GM skew fluxes 504 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 505 ! Limit by slope *relative to geopotentials* by rn_slpmax, and mask by psi-point 506 ! masked by umask taken at the level of dz(rho) 507 ! 508 ! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti) 509 ! 510 zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked 511 ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp) 512 ! 513 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 514 zti_coord = znot_thru_surface * ( gdept(ji+1,jj ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 515 ztj_coord = znot_thru_surface * ( gdept(ji ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) ! unmasked 516 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 517 ztj_g_raw = ztj_raw - ztj_coord 518 ! additional limit required in bilaplacian case 519 ze3_e1 = e3w(ji+ip,jj ,jk+kp,Kmm) * r1_e1u(ji,jj) 520 ze3_e2 = e3w(ji ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj) 521 ! NB: hard coded factor 5 (can be a namelist parameter...) 522 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) 523 ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw ) 524 ! 525 ! Below ML use limited zti_g as is & mask 526 ! Inside ML replace by linearly reducing sx_mlb towards surface & mask 527 ! 528 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 529 zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1 530 ! ! otherwise zfact=0 531 zti_g_lim = ( zfacti * zti_g_lim & 532 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 533 & * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 534 ztj_g_lim = ( zfactj * ztj_g_lim & 535 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 536 & * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 537 ! 538 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim 539 triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim 540 ! 541 ! Get coefficients of isoneutral diffusion tensor 542 ! 1. Utilise gradients *relative* to s-coordinate, so add t-point slopes (*subtract* depth gradients) 543 ! 2. We require that isoneutral diffusion gives no vertical buoyancy flux 544 ! i.e. 33 term = (real slope* 31, 13 terms) 545 ! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms 546 ! Equivalent to tapering A_iso = sx_limited**2/(real slope)**2 547 ! 548 zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces 549 ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp) 550 ! 551 IF( ln_triad_iso ) THEN 552 zti_raw = zti_lim*zti_lim / zti_raw 553 ztj_raw = ztj_lim*ztj_lim / ztj_raw 554 zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw ) 555 ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw ) 556 zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw 557 ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw 558 ENDIF 559 ! ! switching triad scheme 560 zisw = (1._wp - rn_sw_triad ) + rn_sw_triad & 561 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) ) 562 zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad & 563 & * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) ) 564 ! 565 triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw 566 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 567 ! 568 zbu = e1e2u(ji ,jj ) * e3u(ji ,jj ,jk ,Kmm) 569 zbv = e1e2v(ji ,jj ) * e3v(ji ,jj ,jk ,Kmm) 570 zbti = e1e2t(ji+ip,jj ) * e3w(ji+ip,jj ,jk+kp,Kmm) 571 zbtj = e1e2t(ji ,jj+jp) * e3w(ji ,jj+jp,jk+kp,Kmm) 572 ! 573 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked 574 wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim 575 END_2D 607 576 END DO 608 577 END DO … … 611 580 wslp2(:,:,1) = 0._wp ! force the surface wslp to zero 612 581 613 CALL lbc_lnk( 'ldfslp', wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked582 CALL lbc_lnk( 'ldfslp', wslp2, 'W', 1.0_wp ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked 614 583 ! 615 584 IF( ln_timing ) CALL timing_stop('ldf_slp_triad') … … 618 587 619 588 620 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )589 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm ) 621 590 !!---------------------------------------------------------------------- 622 591 !! *** ROUTINE ldf_slp_mxl *** … … 638 607 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 639 608 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 609 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 640 610 !! 641 611 INTEGER :: ji , jj , jk ! dummy loop indices … … 658 628 ! 659 629 ! !== 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 630 DO_3D( 1, 1, 1, 1, 1, jpk ) 631 ik = nmln(ji,jj) - 1 632 IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp 633 ELSE ; omlmask(ji,jj,jk) = 0._wp 634 ENDIF 635 END_3D 670 636 671 637 … … 680 646 !----------------------------------------------------------------------- 681 647 ! 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 648 DO_2D( 0, 0, 0, 0 ) 649 ! !== Slope at u- & v-points just below the Mixed Layer ==! 650 ! 651 ! !- vertical density gradient for u- and v-slopes (from dzr at T-point) 652 iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1) 653 ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) ! 654 zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) ) 655 zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) ) 656 ! !- horizontal density gradient at u- & v-points 657 zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj) 658 zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj) 659 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 660 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 661 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau ) ) 662 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav ) ) 663 ! !- Slope at u- & v-points (uslpml, vslpml) 664 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) 665 vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv) 666 ! 667 ! !== i- & j-slopes at w-points just below the Mixed Layer ==! 668 ! 669 ik = MIN( nmln(ji,jj) + 1, jpk ) 670 ikm1 = MAX( 1, ik-1 ) 671 ! !- vertical density gradient for w-slope (from N^2) 672 zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. ) 673 ! !- horizontal density i- & j-gradient at w-points 674 zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) & 675 & + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj) 676 zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) & 677 & + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj) 678 zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) & 679 & + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik) 680 zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) & 681 & + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik) 682 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 683 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 684 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai ) ) 685 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj ) ) 686 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 687 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) 688 wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik) 689 END_2D 726 690 !!gm this lbc_lnk should be useless.... 727 CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1. , vslpml , 'V', -1. , wslpiml, 'W', -1. , wslpjml, 'W', -1.)691 CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1.0_wp , vslpml , 'V', -1.0_wp , wslpiml, 'W', -1.0_wp , wslpjml, 'W', -1.0_wp ) 728 692 ! 729 693 END SUBROUTINE ldf_slp_mxl … … 785 749 ! DO jk = 1, jpk 786 750 ! 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.5751 ! DO ji = 2, jpim1 ! vector opt. 752 ! uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 753 ! vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 754 ! 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 755 ! 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 756 ! END DO 793 757 ! END DO
Note: See TracChangeset
for help on using the changeset viewer.
