- Timestamp:
- 2019-11-22T15:29:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/LDF/ldfslp.F90
r10425 r11949 81 81 CONTAINS 82 82 83 SUBROUTINE ldf_slp( kt, prd, pn2 )83 SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm ) 84 84 !!---------------------------------------------------------------------- 85 85 !! *** ROUTINE ldf_slp *** … … 107 107 !!---------------------------------------------------------------------- 108 108 INTEGER , INTENT(in) :: kt ! ocean time-step index 109 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 109 110 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density 110 111 REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.) … … 171 172 ! 172 173 ! !== Slopes just below the mixed layer ==! 173 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml174 CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm ) ! output: uslpml, vslpml, wslpiml, wslpjml 174 175 175 176 … … 205 206 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 206 207 ! ! + 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 ) )208 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau ) ) 209 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav ) ) 209 210 ! ! uslp and vslp output in zwz and zww, resp. 210 211 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 211 212 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 212 213 ! 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)) )214 zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) ) & 215 - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) - e3u(ji,jj,miku(ji,jj),Kmm) ) 216 zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) ) & 217 - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) - e3v(ji,jj,mikv(ji,jj),Kmm) ) 217 218 ! 218 219 zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) & … … 224 225 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 225 226 ! 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. ) )227 ! zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 228 ! zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 228 229 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 229 230 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) … … 296 297 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 297 298 ! ! + 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 ) )299 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai ) ) 300 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj ) ) 300 301 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 301 302 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 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 ) 303 304 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk) 304 305 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk) … … 375 376 376 377 377 SUBROUTINE ldf_slp_triad ( kt )378 SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm ) 378 379 !!---------------------------------------------------------------------- 379 380 !! *** ROUTINE ldf_slp_triad *** … … 390 391 !!---------------------------------------------------------------------- 391 392 INTEGER, INTENT( in ) :: kt ! ocean time-step index 393 INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices 392 394 !! 393 395 INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices … … 419 421 DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set 420 422 DO ji = 1, fs_jpim1 ! vector opt. 421 zdit = ( ts b(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point422 zdis = ( ts b(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )423 zdjt = ( ts b(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point424 zdjs = ( ts b(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )423 zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! i-gradient of T & S at u-point 424 zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 425 zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! j-gradient of T & S at v-point 426 zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) ) 425 427 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 428 zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj) … … 452 454 DO ji = 1, jpi ! vector opt. 453 455 IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp 454 zdkt = ( ts b(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) )455 zdks = ( ts b(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) )456 zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) ) 457 zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) ) 456 458 ELSE 457 459 zdkt = 0._wp ! 1st level gradient set to zero … … 460 462 zdzrho_raw = ( - rab_b(ji,jj,jk+kp,jp_tem) * zdkt & 461 463 & + rab_b(ji,jj,jk+kp,jp_sal) * zdks & 462 & ) / e3w _n(ji,jj,jk+kp)464 & ) / e3w(ji,jj,jk+kp,Kmm) 463 465 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 464 466 END DO … … 470 472 DO ji = 1, jpi 471 473 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)474 z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm) 473 475 END DO 474 476 END DO … … 499 501 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 500 502 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 & - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 504 ze3_e1 = e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj) 503 505 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 506 ENDIF … … 509 511 ELSE 510 512 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 & - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 514 ze3_e2 = e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj) 513 515 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 516 ENDIF … … 542 544 ! 543 545 ! 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) ! unmasked546 zti_coord = znot_thru_surface * ( gdept(ji+1,jj ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj) 547 ztj_coord = znot_thru_surface * ( gdept(ji ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) ! unmasked 546 548 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 547 549 ztj_g_raw = ztj_raw - ztj_coord 548 550 ! 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 ze3_e1 = e3w(ji+ip,jj ,jk+kp,Kmm) * r1_e1u(ji,jj) 552 ze3_e2 = e3w(ji ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj) 551 553 ! NB: hard coded factor 5 (can be a namelist parameter...) 552 554 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) … … 561 563 zti_g_lim = ( zfacti * zti_g_lim & 562 564 & + ( 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)565 & * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 564 566 ztj_g_lim = ( zfactj * ztj_g_lim & 565 567 & + ( 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)568 & * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 567 569 ! 568 570 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim … … 596 598 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 597 599 ! 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)600 zbu = e1e2u(ji ,jj ) * e3u(ji ,jj ,jk ,Kmm) 601 zbv = e1e2v(ji ,jj ) * e3v(ji ,jj ,jk ,Kmm) 602 zbti = e1e2t(ji+ip,jj ) * e3w(ji+ip,jj ,jk+kp,Kmm) 603 zbtj = e1e2t(ji ,jj+jp) * e3w(ji ,jj+jp,jk+kp,Kmm) 602 604 ! 603 605 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked … … 618 620 619 621 620 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )622 SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm ) 621 623 !!---------------------------------------------------------------------- 622 624 !! *** ROUTINE ldf_slp_mxl *** … … 638 640 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts) 639 641 REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point) 642 INTEGER , INTENT(in) :: Kmm ! ocean time level indices 640 643 !! 641 644 INTEGER :: ji , jj , jk ! dummy loop indices … … 694 697 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 695 698 ! 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 ) )699 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau ) ) 700 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav ) ) 698 701 ! !- Slope at u- & v-points (uslpml, vslpml) 699 702 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 717 720 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 718 721 ! 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 ) )722 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai ) ) 723 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj ) ) 721 724 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 722 725 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) … … 786 789 ! DO jj = 2, jpjm1 787 790 ! 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.5791 ! uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 792 ! vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 793 ! 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 794 ! 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 795 ! END DO 793 796 ! END DO
Note: See TracChangeset
for help on using the changeset viewer.