- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90
r5836 r6060 73 73 74 74 !! * Substitutions 75 # include "domzgr_substitute.h90"76 75 # include "vectopt_loop_substitute.h90" 77 76 !!---------------------------------------------------------------------- … … 181 180 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0 182 181 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 183 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/ fse3u(ji,jj,jk)* ABS( zau ) )184 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/ fse3v(ji,jj,jk)* ABS( zav ) )182 zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,jk)* ABS( zau ) ) 183 zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,jk)* ABS( zav ) ) 185 184 ! ! uslp and vslp output in zwz and zww, resp. 186 185 zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) ) 187 186 zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) ) 188 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) &189 & + zfi * uslpml(ji,jj) &190 & * 0.5_wp * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk)-fse3u(ji,jj,1) ) &187 zwz(ji,jj,jk) = ( ( 1. - zfi) * zau / ( zbu - zeps ) & 188 & + zfi * uslpml(ji,jj) & 189 & * 0.5_wp * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk)-e3u_n(ji,jj,1) ) & 191 190 & / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 5._wp ) ) * umask(ji,jj,jk) 192 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) &193 & + zfj * vslpml(ji,jj) &194 & * 0.5_wp * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk)-fse3v(ji,jj,1) ) &191 zww(ji,jj,jk) = ( ( 1. - zfj) * zav / ( zbv - zeps ) & 192 & + zfj * vslpml(ji,jj) & 193 & * 0.5_wp * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk)-e3v_n(ji,jj,1) ) & 195 194 & / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 5. ) ) * vmask(ji,jj,jk) 196 195 !!gm modif to suppress omlmask.... (as in Griffies case) … … 198 197 ! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp ) 199 198 ! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp ) 200 ! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )201 ! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )199 ! zci = 0.5 * ( gdept_n(ji+1,jj,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) ) 200 ! zcj = 0.5 * ( gdept_n(ji,jj+1,jk)+gdept_n(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) ) 202 201 ! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk) 203 202 ! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk) … … 270 269 ! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0. 271 270 ! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 272 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/ fse3w(ji,jj,jk)* ABS( zai ) )273 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/ fse3w(ji,jj,jk)* ABS( zaj ) )271 zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zai ) ) 272 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,jk)* ABS( zaj ) ) 274 273 ! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.) 275 274 zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0 276 zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp )275 zck = gdepw_n(ji,jj,jk) / MAX( hmlp(ji,jj), 10._wp ) 277 276 zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * tmask(ji,jj,jk) 278 277 zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * tmask(ji,jj,jk) … … 281 280 ! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0. 282 281 ! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp ) 283 ! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )282 ! zck = gdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. ) 284 283 ! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk) 285 284 ! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk) … … 441 440 zdks = 0._wp 442 441 ENDIF 443 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / fse3w(ji,jj,jk+kp)442 zdzrho_raw = ( - zalbet(ji,jj,jk) * zdkt + zbeta0*zdks ) / e3w_n(ji,jj,jk+kp) 444 443 zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln 445 444 END DO … … 451 450 DO ji = 1, jpi 452 451 jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth 453 z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)452 z1_mlbw(ji,jj) = 1._wp / gdepw_n(ji,jj,jk) 454 453 END DO 455 454 END DO … … 480 479 ! Add s-coordinate slope at t-points (do this by *subtracting* gradient of depth) 481 480 zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) & 482 & - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk)483 ze3_e1 = fse3w(ji+ip,jj,jk-kp) * r1_e1u(ji,jj)481 & - ( gdept_n(ji+1,jj,jk-kp) - gdept_n(ji,jj,jk-kp) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk) 482 ze3_e1 = e3w_n(ji+ip,jj,jk-kp) * r1_e1u(ji,jj) 484 483 zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw ) 485 484 ENDIF … … 490 489 ELSE 491 490 ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) & 492 & - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk)493 ze3_e2 = fse3w(ji,jj+jp,jk-kp) / e2v(ji,jj)491 & - ( gdept_n(ji,jj+1,jk-kp) - gdept_n(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk) 492 ze3_e2 = e3w_n(ji,jj+jp,jk-kp) / e2v(ji,jj) 494 493 ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw ) 495 494 ENDIF … … 523 522 ! 524 523 ! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface 525 zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) * r1_e1u(ji,jj)526 ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked524 zti_coord = znot_thru_surface * ( gdept_n(ji+1,jj ,jk) - gdept_n(ji,jj,jk) ) * r1_e1u(ji,jj) 525 ztj_coord = znot_thru_surface * ( gdept_n(ji ,jj+1,jk) - gdept_n(ji,jj,jk) ) * r1_e2v(ji,jj) ! unmasked 527 526 zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces 528 527 ztj_g_raw = ztj_raw - ztj_coord 529 528 ! additional limit required in bilaplacian case 530 ze3_e1 = fse3w(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj)531 ze3_e2 = fse3w(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj)529 ze3_e1 = e3w_n(ji+ip,jj ,jk+kp) * r1_e1u(ji,jj) 530 ze3_e2 = e3w_n(ji ,jj+jp,jk+kp) * r1_e2v(ji,jj) 532 531 ! NB: hard coded factor 5 (can be a namelist parameter...) 533 532 zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw ) … … 542 541 zti_g_lim = ( zfacti * zti_g_lim & 543 542 & + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) & 544 & * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)543 & * gdepw_n(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp) 545 544 ztj_g_lim = ( zfactj * ztj_g_lim & 546 545 & + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) & 547 & * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)546 & * gdepw_n(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp) 548 547 ! 549 548 triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim … … 577 576 triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw 578 577 ! 579 zbu = e1e2u(ji ,jj ) * fse3u(ji ,jj ,jk )580 zbv = e1e2v(ji ,jj ) * fse3v(ji ,jj ,jk )581 zbti = e1e2t(ji+ip,jj ) * fse3w(ji+ip,jj ,jk+kp)582 zbtj = e1e2t(ji ,jj+jp) * fse3w(ji ,jj+jp,jk+kp)578 zbu = e1e2u(ji ,jj ) * e3u_n(ji ,jj ,jk ) 579 zbv = e1e2v(ji ,jj ) * e3v_n(ji ,jj ,jk ) 580 zbti = e1e2t(ji+ip,jj ) * e3w_n(ji+ip,jj ,jk+kp) 581 zbtj = e1e2t(ji ,jj+jp) * e3w_n(ji ,jj+jp,jk+kp) 583 582 ! 584 583 wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked … … 682 681 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0 683 682 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 684 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/ fse3u(ji,jj,iku)* ABS( zau ) )685 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/ fse3v(ji,jj,ikv)* ABS( zav ) )683 zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u_n(ji,jj,iku)* ABS( zau ) ) 684 zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v_n(ji,jj,ikv)* ABS( zav ) ) 686 685 ! !- Slope at u- & v-points (uslpml, vslpml) 687 686 uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku) … … 705 704 ! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0. 706 705 ! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt) 707 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/ fse3w(ji,jj,ik)* ABS( zai ) )708 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/ fse3w(ji,jj,ik)* ABS( zaj ) )706 zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zai ) ) 707 zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w_n(ji,jj,ik)* ABS( zaj ) ) 709 708 ! !- i- & j-slope at w-points (wslpiml, wslpjml) 710 709 wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik) … … 767 766 768 767 !!gm I no longer understand this..... 769 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. ( lk_vvl.AND. ln_rstart) ) THEN768 !!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (.NOT.ln_linssh .AND. ln_rstart) ) THEN 770 769 ! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces' 771 770 ! … … 775 774 ! 776 775 ! ! set the slope of diffusion to the slope of s-surfaces 777 ! ! ( c a u t i o n : minus sign as fsdep has positive value )776 ! ! ( c a u t i o n : minus sign as dep has positive value ) 778 777 ! DO jk = 1, jpk 779 778 ! DO jj = 2, jpjm1 780 779 ! DO ji = fs_2, fs_jpim1 ! vector opt. 781 ! uslp (ji,jj,jk) = - ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)782 ! vslp (ji,jj,jk) = - ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)783 ! wslpi(ji,jj,jk) = - ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5784 ! wslpj(ji,jj,jk) = - ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5780 ! uslp (ji,jj,jk) = - ( gdept_n(ji+1,jj,jk) - gdept_n(ji ,jj ,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk) 781 ! vslp (ji,jj,jk) = - ( gdept_n(ji,jj+1,jk) - gdept_n(ji ,jj ,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk) 782 ! 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.5 783 ! 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.5 785 784 ! END DO 786 785 ! END DO
Note: See TracChangeset
for help on using the changeset viewer.
