Changeset 10883 for NEMO/branches/2019
- Timestamp:
- 2019-04-18T14:29:58+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE
- Files:
-
- 8 edited
-
ZDF/zdfgls.F90 (modified) (13 diffs)
-
ZDF/zdfosm.F90 (modified) (61 diffs)
-
ZDF/zdfphy.F90 (modified) (7 diffs)
-
ZDF/zdfric.F90 (modified) (4 diffs)
-
ZDF/zdfsh2.F90 (modified) (4 diffs)
-
ZDF/zdftke.F90 (modified) (21 diffs)
-
nemogcm.F90 (modified) (1 diff)
-
step.F90 (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfgls.F90
r10425 r10883 124 124 125 125 126 SUBROUTINE zdf_gls( kt, p_sh2, p_avm, p_avt )126 SUBROUTINE zdf_gls( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 127 127 !!---------------------------------------------------------------------- 128 128 !! *** ROUTINE zdf_gls *** … … 134 134 !! 135 135 INTEGER , INTENT(in ) :: kt ! ocean time step 136 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 136 137 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 137 138 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) … … 176 177 zmsku = ( 2._wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) ) 177 178 zmskv = ( 2._wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) ! (CAUTION: CdU<0) 178 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( u b(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 &179 & + ( zmskv*( v b(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 )179 ustar2_bot(ji,jj) = - rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 180 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 180 181 END DO 181 182 END DO … … 185 186 zmsku = ( 2. - umask(ji-1,jj,mikt(ji,jj)) * umask(ji,jj,mikt(ji,jj)) ) 186 187 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) ! (CAUTION: CdU<0) 187 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( u b(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 &188 & + ( zmskv*( v b(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 )188 ustar2_top(ji,jj) = - rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 189 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 189 190 END DO 190 191 END DO … … 222 223 DO jj = 2, jpjm1 223 224 DO ji = fs_2, fs_jpim1 ! vector opt. 224 zup = hmxl_n(ji,jj,jk) * gdepw _n(ji,jj,mbkt(ji,jj)+1)225 zdown = vkarmn * gdepw _n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) )225 zup = hmxl_n(ji,jj,jk) * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) 226 zdown = vkarmn * gdepw(ji,jj,jk,Kmm) * ( -gdepw(ji,jj,jk,Kmm) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) ) 226 227 zcoef = ( zup / MAX( zdown, rsmall ) ) 227 228 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) … … 276 277 zcof = rfact_tke * tmask(ji,jj,jk) 277 278 ! ! lower diagonal, in fact not used for jk = 2 (see surface conditions) 278 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t _n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )279 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 279 280 ! ! upper diagonal, in fact not used for jk = ibotm1 (see bottom conditions) 280 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t _n(ji,jj,jk ) * e3w_n(ji,jj,jk) )281 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 281 282 ! ! diagonal 282 283 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 306 307 ! 307 308 ! One level below 308 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw _n(:,:,2)) &309 en (:,:,2) = MAX( rc02r * ustar2_surf(:,:) * ( 1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw(:,:,2,Kmm)) & 309 310 & / zhsro(:,:) )**(1.5_wp*ra_sf) )**(2._wp/3._wp) , rn_emin ) 310 311 zd_lw(:,:,2) = 0._wp … … 325 326 zdiag(:,:,2) = zdiag(:,:,2) + zd_lw(:,:,2) ! Remove zd_lw from zdiag 326 327 zd_lw(:,:,2) = 0._wp 327 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept _n(:,:,1)/zhsro(:,:)) ))328 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:)) )) 328 329 zflxs(:,:) = rsbc_tke2 * ustar2_surf(:,:)**1.5_wp * zkar(:,:) & 329 & * ( ( zhsro(:,:)+gdept _n(:,:,1) ) / zhsro(:,:) )**(1.5_wp*ra_sf)330 !!gm why not : * ( 1._wp + gdept _n(:,:,1) / zhsro(:,:) )**(1.5_wp*ra_sf)331 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w _n(:,:,2)330 & * ( ( zhsro(:,:)+gdept(:,:,1,Kmm) ) / zhsro(:,:) )**(1.5_wp*ra_sf) 331 !!gm why not : * ( 1._wp + gdept(:,:,1,Kmm) / zhsro(:,:) )**(1.5_wp*ra_sf) 332 en(:,:,2) = en(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 332 333 ! 333 334 ! … … 526 527 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) 527 528 ! ! lower diagonal 528 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t _n(ji,jj,jk-1) * e3w_n(ji,jj,jk) )529 zd_lw(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) ) / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk,Kmm) ) 529 530 ! ! upper diagonal 530 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t _n(ji,jj,jk ) * e3w_n(ji,jj,jk) )531 zd_up(ji,jj,jk) = zcof * ( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) ) / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk,Kmm) ) 531 532 ! ! diagonal 532 533 zdiag(ji,jj,jk) = 1._wp - zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) + rdt * zdiss * wmask(ji,jj,jk) … … 554 555 ! 555 556 ! One level below 556 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw _n(:,:,2)/zhsro(:,:) )))557 zdep (:,:) = (zhsro(:,:) + gdepw _n(:,:,2)) * zkar(:,:)557 zkar (:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdepw(:,:,2,Kmm)/zhsro(:,:) ))) 558 zdep (:,:) = (zhsro(:,:) + gdepw(:,:,2,Kmm)) * zkar(:,:) 558 559 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 559 560 zd_lw(:,:,2) = 0._wp … … 575 576 ! 576 577 ! Set psi vertical flux at the surface: 577 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept _n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope578 zdep (:,:) = ((zhsro(:,:) + gdept _n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)578 zkar (:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-EXP(-rtrans*gdept(:,:,1,Kmm)/zhsro(:,:) )) ! Lengh scale slope 579 zdep (:,:) = ((zhsro(:,:) + gdept(:,:,1,Kmm)) / zhsro(:,:))**(rmm*ra_sf) 579 580 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 580 581 zdep (:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*p_avm(:,:,1)+zwall_psi(:,:,2)*p_avm(:,:,2)) * & 581 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept _n(:,:,1))**(rnn-1.)582 & ustar2_surf(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept(:,:,1,Kmm))**(rnn-1.) 582 583 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 583 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w _n(:,:,2)584 psi (:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w(:,:,2,Kmm) 584 585 ! 585 586 END SELECT … … 607 608 ! 608 609 ! Just above last level, Dirichlet condition again (GOTM like) 609 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t _n(ji,jj,ibotm1) )610 zdep(ji,jj) = vkarmn * ( r_z0_bot + e3t(ji,jj,ibotm1,Kmm) ) 610 611 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 611 612 zd_lw(ji,jj,ibotm1) = 0._wp … … 635 636 ! 636 637 ! Set psi vertical flux at the bottom: 637 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t _n(ji,jj,ibotm1)638 zdep(ji,jj) = r_z0_bot + 0.5_wp*e3t(ji,jj,ibotm1,Kmm) 638 639 zflxb = rsbc_psi2 * ( p_avm(ji,jj,ibot) + p_avm(ji,jj,ibotm1) ) & 639 640 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 640 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w _n(ji,jj,ibotm1)641 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w(ji,jj,ibotm1,Kmm) 641 642 END DO 642 643 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90
r10425 r10883 42 42 !!---------------------------------------------------------------------- 43 43 USE oce ! ocean dynamics and active tracers 44 ! uses w nfrom previous time step (which is now wb) to calculate hbl44 ! uses ww from previous time step (which is now wb) to calculate hbl 45 45 USE dom_oce ! ocean space and time domain 46 46 USE zdf_oce ! ocean vertical physics … … 122 122 123 123 124 SUBROUTINE zdf_osm( kt, p_avm, p_avt )124 SUBROUTINE zdf_osm( kt, Kbb, Kmm, p_avm, p_avt ) 125 125 !!---------------------------------------------------------------------- 126 126 !! *** ROUTINE zdf_osm *** … … 157 157 !! the equation number. (LMD94, here after) 158 158 !!---------------------------------------------------------------------- 159 INTEGER , INTENT(in ) :: kt ! ocean time step 159 INTEGER , INTENT(in ) :: kt ! ocean time step 160 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 160 161 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 161 162 !! … … 314 315 zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1) 315 316 ! Upwards surface salinity flux for non-local term 316 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts n(ji,jj,1,jp_sal) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)317 zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1) 317 318 ! Non radiative upwards surface buoyancy flux 318 319 zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) - grav * zbeta * zws0(ji,jj) … … 407 408 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 408 409 ! BL must be always 2 levels deep. 409 hbl(:,:) = MAX(hbl(:,:), gdepw _n(:,:,3) )410 hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) ) 410 411 ibld(:,:) = 3 411 412 DO jk = 4, jpkm1 412 413 DO jj = 2, jpjm1 413 414 DO ji = 2, jpim1 414 IF ( hbl(ji,jj) >= gdepw _n(ji,jj,jk) ) THEN415 IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 415 416 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 416 417 ENDIF … … 430 431 zthick=0._wp 431 432 DO jm = 2, ibld(ji,jj) 432 zthick=zthick+e3t _n(ji,jj,jm)433 zt = zt + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)434 zs = zs + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)435 zu = zu + e3t _n(ji,jj,jm) &436 & * ( u b(ji,jj,jm) + ub(ji - 1,jj,jm) ) &433 zthick=zthick+e3t(ji,jj,jm,Kmm) 434 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 435 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 436 zu = zu + e3t(ji,jj,jm,Kmm) & 437 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 437 438 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 438 zv = zv + e3t _n(ji,jj,jm) &439 & * ( v b(ji,jj,jm) + vb(ji,jj - 1,jm) ) &439 zv = zv + e3t(ji,jj,jm,Kmm) & 440 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 440 441 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 441 442 END DO … … 444 445 zu_bl(ji,jj) = zu / zthick 445 446 zv_bl(ji,jj) = zv / zthick 446 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_tem)447 zds_bl(ji,jj) = zs_bl(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_sal)448 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( u b(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj)) ) &447 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 448 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 449 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 449 450 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 450 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( v b(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj)) ) &451 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 451 452 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 452 453 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) … … 487 488 ibld(:,:) = 3 488 489 489 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - w n(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it490 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 490 491 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:)) 491 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + w n(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom492 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 492 493 493 494 DO jk = 4, jpkm1 494 495 DO jj = 2, jpjm1 495 496 DO ji = 2, jpim1 496 IF ( zhbl_t(ji,jj) >= gdepw _n(ji,jj,jk) ) THEN497 IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN 497 498 ibld(ji,jj) = MIN(mbkt(ji,jj), jk) 498 499 ENDIF … … 520 521 521 522 DO jk = imld(ji,jj), ibld(ji,jj) 522 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts n(ji,jj,jm,jp_tem) ) &523 & - zbeta * ( zs_bl(ji,jj) - ts n(ji,jj,jm,jp_sal) ) ), 0.0 ) + zvel_max524 525 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w _n(ji,jj,jk) )523 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 524 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max 525 526 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 526 527 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 527 528 528 IF ( zhbl_s >= gdepw _n(ji,jj,jm+1) ) jm = jm + 1529 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 529 530 END DO 530 531 hbl(ji,jj) = zhbl_s … … 534 535 ! stable 535 536 DO jk = imld(ji,jj), ibld(ji,jj) 536 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts n(ji,jj,jm,jp_tem) ) &537 & - zbeta * ( zs_bl(ji,jj) - ts n(ji,jj,jm,jp_sal) ) ), 0.0 ) &537 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) & 538 & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) & 538 539 & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s 539 540 … … 543 544 & + ( ( 0.32 / 3.0 ) * EXP( - 2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) ) & 544 545 & - ( 0.32 / 3.0 - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s ) ) ) & 545 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w _n(ji,jj,jk) / zdhdt(ji,jj) ! ALMG to investigate whether need to include wnhere546 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj) ! ALMG to investigate whether need to include ww here 546 547 547 548 zhbl_s = MIN(zhbl_s, ht_n(ji,jj)) 548 IF ( zhbl_s >= gdepw _n(ji,jj,jm) ) jm = jm + 1549 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 549 550 END DO 550 hbl(ji,jj) = MAX(zhbl_s, gdepw _n(ji,jj,3) )551 hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) ) 551 552 ibld(ji,jj) = MAX(jm, 3 ) 552 553 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) … … 558 559 hbli(ji,jj) = hbl(ji,jj) 559 560 ELSE 560 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw _n(ji,jj,3) )561 hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) ) 561 562 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj) 562 563 ENDIF 563 564 ENDIF 564 zhbl(ji,jj) = gdepw _n(ji,jj,ibld(ji,jj))565 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 565 566 END DO 566 567 END DO … … 581 582 zthick=0._wp 582 583 DO jm = 2, ibld(ji,jj) 583 zthick=zthick+e3t _n(ji,jj,jm)584 zt = zt + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)585 zs = zs + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)586 zu = zu + e3t _n(ji,jj,jm) &587 & * ( u b(ji,jj,jm) + ub(ji - 1,jj,jm) ) &584 zthick=zthick+e3t(ji,jj,jm,Kmm) 585 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 586 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 587 zu = zu + e3t(ji,jj,jm,Kmm) & 588 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 588 589 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 589 zv = zv + e3t _n(ji,jj,jm) &590 & * ( v b(ji,jj,jm) + vb(ji,jj - 1,jm) ) &590 zv = zv + e3t(ji,jj,jm,Kmm) & 591 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 591 592 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 592 593 END DO … … 595 596 zu_bl(ji,jj) = zu / zthick 596 597 zv_bl(ji,jj) = zv / zthick 597 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_tem)598 zds_bl(ji,jj) = zs_bl(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_sal)599 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( u b(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj)) ) &598 zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 599 zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 600 zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 600 601 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 601 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( v b(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj)) ) &602 zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 602 603 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 603 604 zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj) 604 zhbl(ji,jj) = gdepw _n(ji,jj,ibld(ji,jj))605 zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm) 605 606 IF ( lconv(ji,jj) ) THEN 606 607 IF ( zdb_bl(ji,jj) > 0._wp )THEN … … 616 617 zwb_ent(ji,jj) = 0._wp 617 618 ENDIF 618 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t _n(ji,jj,ibld(ji,jj)) ) , 1 )619 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 619 620 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 620 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))621 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 621 622 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 622 623 ELSE ! IF (zdb_bl) 623 624 imld(ji,jj) = ibld(ji,jj) - 1 624 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))625 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 625 626 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 626 627 ENDIF … … 632 633 zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) & 633 634 & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01 , 0.2 ) 634 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t _n(ji,jj,ibld(ji,jj)) ) , 1 )635 inhml = MAX( INT( zari * zhbl(ji,jj) / e3t(ji,jj,ibld(ji,jj),Kmm) ) , 1 ) 635 636 imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1) 636 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))637 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 637 638 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 638 639 ELSE 639 640 imld(ji,jj) = ibld(ji,jj) - 1 640 zhml(ji,jj) = gdepw _n(ji,jj,imld(ji,jj))641 zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm) 641 642 zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj) 642 643 ENDIF ! IF (zdb_bl > 0.0) … … 665 666 zthick=0._wp 666 667 DO jm = 2, imld(ji,jj) 667 zthick=zthick+e3t _n(ji,jj,jm)668 zt = zt + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)669 zs = zs + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)670 zu = zu + e3t _n(ji,jj,jm) &671 & * ( u b(ji,jj,jm) + ub(ji - 1,jj,jm) ) &668 zthick=zthick+e3t(ji,jj,jm,Kmm) 669 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 670 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 671 zu = zu + e3t(ji,jj,jm,Kmm) & 672 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 672 673 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 673 zv = zv + e3t _n(ji,jj,jm) &674 & * ( v b(ji,jj,jm) + vb(ji,jj - 1,jm) ) &674 zv = zv + e3t(ji,jj,jm,Kmm) & 675 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 675 676 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 676 677 END DO … … 679 680 zu_ml(ji,jj) = zu / zthick 680 681 zv_ml(ji,jj) = zv / zthick 681 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_tem)682 zds_ml(ji,jj) = zs_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_sal)683 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( u b(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj)) ) &682 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 683 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 684 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 684 685 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 685 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( v b(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj)) ) &686 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 686 687 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 687 688 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) … … 696 697 zthick=0._wp 697 698 DO jm = 2, imld(ji,jj) 698 zthick=zthick+e3t _n(ji,jj,jm)699 zt = zt + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)700 zs = zs + e3t _n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)701 zu = zu + e3t _n(ji,jj,jm) &702 & * ( u b(ji,jj,jm) + ub(ji - 1,jj,jm) ) &699 zthick=zthick+e3t(ji,jj,jm,Kmm) 700 zt = zt + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm) 701 zs = zs + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm) 702 zu = zu + e3t(ji,jj,jm,Kmm) & 703 & * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) & 703 704 & / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) ) 704 zv = zv + e3t _n(ji,jj,jm) &705 & * ( v b(ji,jj,jm) + vb(ji,jj - 1,jm) ) &705 zv = zv + e3t(ji,jj,jm,Kmm) & 706 & * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) & 706 707 & / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) ) 707 708 END DO … … 710 711 zu_ml(ji,jj) = zu / zthick 711 712 zv_ml(ji,jj) = zv / zthick 712 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_tem)713 zds_ml(ji,jj) = zs_ml(ji,jj) - ts n(ji,jj,ibld(ji,jj),jp_sal)714 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( u b(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj)) ) &713 zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm) 714 zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm) 715 zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) & 715 716 & / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) ) 716 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( v b(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj)) ) &717 zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) & 717 718 & / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) ) 718 719 zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj) … … 775 776 zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) ) 776 777 DO jk = 2 , ibld(ji,jj) 777 znd = -( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)778 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 778 779 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 779 780 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) … … 791 792 zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj) 792 793 DO jk = 2, ibld(ji,jj) 793 znd = gdepw _n(ji,jj,jk) / zhbl(ji,jj)794 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 794 795 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) 795 796 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 ) … … 801 802 zbgrad = zdb_bl(ji,jj) / zdh(ji,jj) 802 803 DO jk = 2, ibld(ji,jj) 803 znd = -( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)804 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 804 805 zdtdz_pyc(ji,jj,jk) = ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 805 806 zdbdz_pyc(ji,jj,jk) = zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) … … 824 825 & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird 825 826 DO jk = 2 , ibld(ji,jj)-1 826 znd = -( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)827 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 827 828 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -1.75 * ( znd + 0.75 )**2 ) 828 829 zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 ) … … 833 834 zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj) 834 835 DO jk = 2, ibld(ji,jj) 835 znd = gdepw _n(ji,jj,jk) / zhbl(ji,jj)836 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 836 837 IF ( znd < 1.0 ) THEN 837 838 zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 ) … … 880 881 IF ( lconv(ji,jj) ) THEN 881 882 DO jk = 2, imld(ji,jj) ! mixed layer diffusivity 882 zznd_ml = gdepw _n(ji,jj,jk) / zhml(ji,jj)883 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 883 884 ! 884 885 zdiffut(ji,jj,jk) = 0.8 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5 … … 890 891 IF ( zdh(ji,jj) > 0._wp ) THEN 891 892 DO jk = imld(ji,jj)+1 , ibld(ji,jj) 892 zznd_pyc = -( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)893 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 893 894 ! 894 895 zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc ) … … 897 898 END DO 898 899 ENDIF 899 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with w ntaken out900 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t _n(ji,jj,ibld(ji,jj))900 ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out 901 zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm) 901 902 ! could be taken out, take account of entrainment represents as a diffusivity 902 903 ! should remove w from here, represents entrainment … … 904 905 ! stable conditions 905 906 DO jk = 2, ibld(ji,jj) 906 zznd_ml = gdepw _n(ji,jj,jk) / zhbl(ji,jj)907 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 907 908 zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5 908 909 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 ) … … 932 933 IF ( lconv(ji,jj) ) THEN 933 934 DO jk = 2, imld(ji,jj) 934 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)935 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 935 936 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj) 936 937 ! … … 967 968 IF ( lconv(ji,jj) ) THEN 968 969 DO jk = 2, imld(ji,jj) 969 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)970 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 970 971 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) & 971 972 & + 0.00125 * EXP ( - zznd_d ) * zsc_uw_2(ji,jj) ) & … … 978 979 ! Stable conditions 979 980 DO jk = 2, ibld(ji,jj) ! corrected to ibld 980 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)981 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 981 982 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) & 982 983 & * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj) … … 1001 1002 IF (lconv(ji,jj) ) THEN 1002 1003 DO jk = 2, imld(ji,jj) 1003 zznd_ml = gdepw _n(ji,jj,jk) / zhml(ji,jj)1004 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1004 1005 ! calculate turbulent length scale 1005 1006 zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) & … … 1035 1036 IF ( lconv(ji,jj) ) THEN 1036 1037 DO jk = 2 , imld(ji,jj) 1037 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)1038 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1038 1039 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) & 1039 1040 & * ( 1.0 - EXP( -0.5 * zznd_d ) ) & … … 1078 1079 ELSE 1079 1080 DO jk = 2, ibld(ji,jj) 1080 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)1081 znd = gdepw _n(ji,jj,jk) / zhbl(ji,jj)1081 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1082 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1082 1083 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + & 1083 1084 & 7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj) … … 1104 1105 IF ( lconv(ji,jj) ) THEN 1105 1106 DO jk = 2, imld(ji,jj) 1106 zznd_ml = gdepw _n(ji,jj,jk) / zhml(ji,jj)1107 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)1107 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj) 1108 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1108 1109 ghamu(ji,jj,jk) = ghamu(ji,jj,jk)& 1109 1110 & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj) … … 1114 1115 ELSE 1115 1116 DO jk = 2, ibld(ji,jj) 1116 znd = gdepw _n(ji,jj,jk) / zhbl(ji,jj)1117 zznd_d = gdepw _n(ji,jj,jk) / dstokes(ji,jj)1117 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1118 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj) 1118 1119 IF ( zznd_d <= 2.0 ) THEN 1119 1120 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 & … … 1141 1142 IF ( lconv(ji,jj) ) THEN 1142 1143 DO jk = 2, ibld(ji,jj) 1143 znd = ( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about1144 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1144 1145 IF ( znd >= 0.0 ) THEN 1145 1146 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) ) … … 1152 1153 ELSE 1153 1154 DO jk = 2, ibld(ji,jj) 1154 znd = ( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about1155 znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about 1155 1156 IF ( znd >= 0.0 ) THEN 1156 1157 ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) ) … … 1171 1172 DO ji = 2, jpim1 1172 1173 DO jk= 2, ibld(ji,jj) 1173 znd = gdepw _n(ji,jj,jk) / zhbl(ji,jj)1174 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj) 1174 1175 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk) 1175 1176 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk) … … 1194 1195 END DO 1195 1196 DO jk = imld(ji,jj), ibld(ji,jj) 1196 znd = -( gdepw _n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)1197 znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj) 1197 1198 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd ) 1198 1199 ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd ) … … 1239 1240 DO jj = 1, jpjm1 1240 1241 DO ji = 1, jpim1 ! vector opt. 1241 z3du(ji,jj,jk) = 0.5 * ( u n(ji,jj,jk-1) - un(ji ,jj,jk) ) &1242 & * ( u b(ji,jj,jk-1) - ub(ji ,jj,jk) ) * wumask(ji,jj,jk) &1243 & / ( e3uw _n(ji,jj,jk) * e3uw_b(ji,jj,jk) )1244 z3dv(ji,jj,jk) = 0.5 * ( v n(ji,jj,jk-1) - vn(ji,jj ,jk) ) &1245 & * ( v b(ji,jj,jk-1) - vb(ji,jj ,jk) ) * wvmask(ji,jj,jk) &1246 & / ( e3vw _n(ji,jj,jk) * e3vw_b(ji,jj,jk) )1242 z3du(ji,jj,jk) = 0.5 * ( uu(ji,jj,jk-1,Kmm) - uu(ji ,jj,jk,Kmm) ) & 1243 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji ,jj,jk,Kbb) ) * wumask(ji,jj,jk) & 1244 & / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) 1245 z3dv(ji,jj,jk) = 0.5 * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj ,jk,Kmm) ) & 1246 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj ,jk,Kbb) ) * wvmask(ji,jj,jk) & 1247 & / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) 1247 1248 END DO 1248 1249 END DO … … 1364 1365 1365 1366 1366 SUBROUTINE zdf_osm_init 1367 SUBROUTINE zdf_osm_init( Kmm ) 1367 1368 !!---------------------------------------------------------------------- 1368 1369 !! *** ROUTINE zdf_osm_init *** … … 1376 1377 !! ** input : Namlist namosm 1377 1378 !!---------------------------------------------------------------------- 1379 INTEGER, INTENT(in) :: Kmm ! time level index (middle) 1380 ! 1378 1381 INTEGER :: ios ! local integer 1379 1382 INTEGER :: ji, jj, jk ! dummy loop indices … … 1423 1426 IF( zdf_osm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' ) 1424 1427 1425 call osm_rst( nit000, 'READ' ) !* read or initialize hbl1428 call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl 1426 1429 1427 1430 IF( ln_zdfddm) THEN … … 1517 1520 1518 1521 1519 SUBROUTINE osm_rst( kt, cdrw )1522 SUBROUTINE osm_rst( kt, Kmm, cdrw ) 1520 1523 !!--------------------------------------------------------------------- 1521 1524 !! *** ROUTINE osm_rst *** … … 1527 1530 !!---------------------------------------------------------------------- 1528 1531 1529 INTEGER, INTENT(in) :: kt 1532 INTEGER , INTENT(in) :: kt ! ocean time step index 1533 INTEGER , INTENT(in) :: Kmm ! ocean time level index (middle) 1530 1534 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1531 1535 … … 1545 1549 id1 = iom_varid( numror, 'wn' , ldstop = .FALSE. ) 1546 1550 IF( id1 > 0 ) THEN ! 'wn' exists; read 1547 CALL iom_get( numror, jpdom_autoglo, 'wn', w n, ldxios = lrxios )1548 WRITE(numout,*) ' ===>>>> : w nread from restart file'1551 CALL iom_get( numror, jpdom_autoglo, 'wn', ww, ldxios = lrxios ) 1552 WRITE(numout,*) ' ===>>>> : ww read from restart file' 1549 1553 ELSE 1550 w n(:,:,:) = 0._wp1551 WRITE(numout,*) ' ===>>>> : w nnot in restart file, set to zero initially'1554 ww(:,:,:) = 0._wp 1555 WRITE(numout,*) ' ===>>>> : ww not in restart file, set to zero initially' 1552 1556 END IF 1553 1557 id1 = iom_varid( numror, 'hbl' , ldstop = .FALSE. ) … … 1568 1572 IF( TRIM(cdrw) == 'WRITE') THEN !* Write hbli into the restart file, then return 1569 1573 IF(lwp) WRITE(numout,*) '---- osm-rst ----' 1570 CALL iom_rstput( kt, nitrst, numrow, 'wn' , w n, ldxios = lwxios )1574 CALL iom_rstput( kt, nitrst, numrow, 'wn' , ww , ldxios = lwxios ) 1571 1575 CALL iom_rstput( kt, nitrst, numrow, 'hbl' , hbl , ldxios = lwxios ) 1572 1576 CALL iom_rstput( kt, nitrst, numrow, 'hbli' , hbli, ldxios = lwxios ) … … 1580 1584 ALLOCATE( imld_rst(jpi,jpj) ) 1581 1585 ! w-level of the mixing and mixed layers 1582 CALL eos_rab( ts n, rab_n )1583 CALL bn2(ts n, rab_n, rn2)1586 CALL eos_rab( ts(:,:,:,:,Kmm), rab_n ) 1587 CALL bn2(ts(:,:,:,:,Kmm), rab_n, rn2) 1584 1588 imld_rst(:,:) = nlb10 ! Initialization to the number of w ocean point 1585 1589 hbl(:,:) = 0._wp ! here hbl used as a dummy variable, integrating vertically N^2 … … 1591 1595 DO ji = 1, jpi 1592 1596 ikt = mbkt(ji,jj) 1593 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w _n(ji,jj,jk)1597 hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm) 1594 1598 IF( hbl(ji,jj) < zN2_c ) imld_rst(ji,jj) = MIN( jk , ikt ) + 1 ! Mixed layer level 1595 1599 END DO … … 1600 1604 DO ji = 1, jpi 1601 1605 iiki = imld_rst(ji,jj) 1602 hbl (ji,jj) = gdepw _n(ji,jj,iiki) * ssmask(ji,jj) ! Turbocline depth1606 hbl (ji,jj) = gdepw(ji,jj,iiki ,Kmm) * ssmask(ji,jj) ! Turbocline depth 1603 1607 END DO 1604 1608 END DO … … 1610 1614 1611 1615 1612 SUBROUTINE tra_osm( kt )1616 SUBROUTINE tra_osm( kt, Kmm, pts, Krhs ) 1613 1617 !!---------------------------------------------------------------------- 1614 1618 !! *** ROUTINE tra_osm *** … … 1620 1624 REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdt, ztrds ! 3D workspace 1621 1625 !!---------------------------------------------------------------------- 1622 INTEGER, INTENT(in) :: kt 1626 INTEGER , INTENT(in) :: kt ! time step index 1627 INTEGER , INTENT(in) :: Kmm, Krhs ! time level indices 1628 REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts ! active tracers and RHS of tracer equation 1629 ! 1623 1630 INTEGER :: ji, jj, jk 1624 1631 ! … … 1630 1637 1631 1638 IF( l_trdtra ) THEN !* Save ta and sa trends 1632 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = tsa(:,:,:,jp_tem)1633 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = tsa(:,:,:,jp_sal)1639 ALLOCATE( ztrdt(jpi,jpj,jpk) ) ; ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) 1640 ALLOCATE( ztrds(jpi,jpj,jpk) ) ; ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) 1634 1641 ENDIF 1635 1642 … … 1638 1645 DO jj = 2, jpjm1 1639 1646 DO ji = 2, jpim1 1640 tsa(ji,jj,jk,jp_tem) = tsa(ji,jj,jk,jp_tem) &1647 pts(ji,jj,jk,jp_tem,Krhs) = pts(ji,jj,jk,jp_tem,Krhs) & 1641 1648 & - ( ghamt(ji,jj,jk ) & 1642 & - ghamt(ji,jj,jk+1) ) /e3t _n(ji,jj,jk)1643 tsa(ji,jj,jk,jp_sal) = tsa(ji,jj,jk,jp_sal) &1649 & - ghamt(ji,jj,jk+1) ) /e3t(ji,jj,jk,Kmm) 1650 pts(ji,jj,jk,jp_sal,Krhs) = pts(ji,jj,jk,jp_sal,Krhs) & 1644 1651 & - ( ghams(ji,jj,jk ) & 1645 & - ghams(ji,jj,jk+1) ) / e3t _n(ji,jj,jk)1652 & - ghams(ji,jj,jk+1) ) / e3t(ji,jj,jk,Kmm) 1646 1653 END DO 1647 1654 END DO … … 1651 1658 ! save the non-local tracer flux trends for diagnostic 1652 1659 IF( l_trdtra ) THEN 1653 ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)1654 ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)1660 ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:) 1661 ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:) 1655 1662 !!bug gm jpttdzdf ==> jpttosm 1656 1663 CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt ) … … 1660 1667 1661 1668 IF(ln_ctl) THEN 1662 CALL prt_ctl( tab3d_1= tsa(:,:,:,jp_tem), clinfo1=' osm - Ta: ', mask1=tmask, &1663 & tab3d_2= tsa(:,:,:,jp_sal), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )1669 CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' osm - Ta: ', mask1=tmask, & 1670 & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' ) 1664 1671 ENDIF 1665 1672 ! … … 1684 1691 1685 1692 1686 SUBROUTINE dyn_osm( kt )1693 SUBROUTINE dyn_osm( kt, Kmm, puu, pvv, Krhs ) 1687 1694 !!---------------------------------------------------------------------- 1688 1695 !! *** ROUTINE dyn_osm *** … … 1693 1700 !! ** Method : ??? 1694 1701 !!---------------------------------------------------------------------- 1695 INTEGER, INTENT(in) :: kt ! 1702 INTEGER , INTENT( in ) :: kt ! ocean time step index 1703 INTEGER , INTENT( in ) :: Kmm, Krhs ! ocean time level indices 1704 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 1696 1705 ! 1697 1706 INTEGER :: ji, jj, jk ! dummy loop indices … … 1708 1717 DO jj = 2, jpjm1 1709 1718 DO ji = 2, jpim1 1710 ua(ji,jj,jk) = ua(ji,jj,jk) &1719 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) & 1711 1720 & - ( ghamu(ji,jj,jk ) & 1712 & - ghamu(ji,jj,jk+1) ) / e3u _n(ji,jj,jk)1713 va(ji,jj,jk) = va(ji,jj,jk) &1721 & - ghamu(ji,jj,jk+1) ) / e3u(ji,jj,jk,Kmm) 1722 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) & 1714 1723 & - ( ghamv(ji,jj,jk ) & 1715 & - ghamv(ji,jj,jk+1) ) / e3v _n(ji,jj,jk)1724 & - ghamv(ji,jj,jk+1) ) / e3v(ji,jj,jk,Kmm) 1716 1725 END DO 1717 1726 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfphy.F90
r10874 r10883 61 61 CONTAINS 62 62 63 SUBROUTINE zdf_phy_init 63 SUBROUTINE zdf_phy_init( Kmm ) 64 64 !!---------------------------------------------------------------------- 65 65 !! *** ROUTINE zdf_phy_init *** … … 70 70 !! set horizontal shape and vertical profile of background mixing coef. 71 71 !!---------------------------------------------------------------------- 72 INTEGER, INTENT(in) :: Kmm ! time level index (middle) 73 ! 72 74 INTEGER :: jk ! dummy loop indices 73 75 INTEGER :: ioptio, ios ! local integers … … 193 195 IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init ; ENDIF 194 196 IF( ln_zdfgls ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_GLS ; CALL zdf_gls_init ; ENDIF 195 IF( ln_zdfosm ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_OSM ; CALL zdf_osm_init ; ENDIF197 IF( ln_zdfosm ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_OSM ; CALL zdf_osm_init( Kmm ) ; ENDIF 196 198 ! 197 199 IF( ioptio /= 1 ) CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' ) … … 218 220 219 221 220 SUBROUTINE zdf_phy( kt )222 SUBROUTINE zdf_phy( kt, Kbb, Kmm ) 221 223 !!---------------------------------------------------------------------- 222 224 !! *** ROUTINE zdf_phy *** … … 230 232 !! bottom stress..... <<<<====verifier ! 231 233 !!---------------------------------------------------------------------- 232 INTEGER, INTENT(in) :: kt ! ocean time-step index 234 INTEGER, INTENT(in) :: kt ! ocean time-step index 235 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 233 236 ! 234 237 INTEGER :: ji, jj, jk ! dummy loop indice … … 254 257 ! 255 258 IF( l_zdfsh2 ) & !* shear production at w-points (energy conserving form) 256 CALL zdf_sh2( ub, vb, un, vn, avm_k, & ! <<== in257 & 259 CALL zdf_sh2( Kbb, Kmm, avm_k, & ! <<== in 260 & zsh2 ) ! ==>> out : shear production 258 261 ! 259 262 SELECT CASE ( nzdf_phy ) !* Vertical eddy viscosity and diffusivity coefficients at w-points 260 CASE( np_RIC ) ; CALL zdf_ric( kt, gdept_n, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz261 CASE( np_TKE ) ; CALL zdf_tke( kt 262 CASE( np_GLS ) ; CALL zdf_gls( kt 263 CASE( np_OSM ) ; CALL zdf_osm( kt 263 CASE( np_RIC ) ; CALL zdf_ric( kt, Kmm, zsh2, avm_k, avt_k ) ! Richardson number dependent Kz 264 CASE( np_TKE ) ; CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! TKE closure scheme for Kz 265 CASE( np_GLS ) ; CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k ) ! GLS closure scheme for Kz 266 CASE( np_OSM ) ; CALL zdf_osm( kt, Kbb, Kmm , avm_k, avt_k ) ! OSMOSIS closure scheme for Kz 264 267 ! CASE( np_CST ) ! Constant Kz (reset avt, avm to the background value) 265 268 ! ! avt_k and avm_k set one for all at initialisation phase … … 318 321 IF( ln_zdfgls ) CALL gls_rst( kt, 'WRITE' ) 319 322 IF( ln_zdfric ) CALL ric_rst( kt, 'WRITE' ) 320 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after w nhas been updated323 ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after ww has been updated 321 324 ENDIF 322 325 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfric.F90
r10068 r10883 112 112 113 113 114 SUBROUTINE zdf_ric( kt, pdept, p_sh2, p_avm, p_avt )114 SUBROUTINE zdf_ric( kt, Kmm, p_sh2, p_avm, p_avt ) 115 115 !!---------------------------------------------------------------------- 116 116 !! *** ROUTINE zdfric *** … … 125 125 !! avt = avm0 / (1 + rn_alp*ri) 126 126 !! with ri = N^2 / dz(u)**2 127 !! = e3w**2 * rn2/[ mi( dk(u b) )+mj( dk(vb) ) ]127 !! = e3w**2 * rn2/[ mi( dk(uu(:,:,:,Kbb)) )+mj( dk(vv(:,:,:,Kbb)) ) ] 128 128 !! avm0= rn_avmri / (1 + rn_alp*Ri)**nn_ric 129 129 !! where ri is the before local Richardson number, … … 152 152 !!---------------------------------------------------------------------- 153 153 INTEGER , INTENT(in ) :: kt ! ocean time-step 154 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdept ! depth of t-point [m]154 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 155 155 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 156 156 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) … … 189 189 DO jj = 2, jpjm1 190 190 DO ji = 2, jpim1 191 IF( pdept(ji,jj,jk) < zh_ekm(ji,jj) ) THEN191 IF( gdept(ji,jj,jk,Kmm) < zh_ekm(ji,jj) ) THEN 192 192 p_avm(ji,jj,jk) = MAX( p_avm(ji,jj,jk), rn_wvmix ) * wmask(ji,jj,jk) 193 193 p_avt(ji,jj,jk) = MAX( p_avt(ji,jj,jk), rn_wtmix ) * wmask(ji,jj,jk) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfsh2.F90
r10069 r10883 11 11 !! zdf_sh2 : compute mixing the shear production term of TKE 12 12 !!---------------------------------------------------------------------- 13 USE oce 13 14 USE dom_oce ! domain: ocean 14 15 ! … … 28 29 CONTAINS 29 30 30 SUBROUTINE zdf_sh2( pub, pvb, pun, pvn, p_avm, p_sh2 )31 SUBROUTINE zdf_sh2( Kbb, Kmm, p_avm, p_sh2 ) 31 32 !!---------------------------------------------------------------------- 32 33 !! *** ROUTINE zdf_sh2 *** … … 47 48 !! References : Bruchard, OM 2002 48 49 !! --------------------------------------------------------------------- 49 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pub, pvb, pun, pvn ! before, now horizontal velocities50 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 50 51 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm ! vertical eddy viscosity (w-points) 51 52 REAL(wp), DIMENSION(:,:,:) , INTENT( out) :: p_sh2 ! shear production of TKE (w-points) … … 59 60 DO ji = 1, jpim1 60 61 zsh2u(ji,jj) = ( p_avm(ji+1,jj,jk) + p_avm(ji,jj,jk) ) & 61 & * ( pun(ji,jj,jk-1) - pun(ji,jj,jk) ) &62 & * ( pub(ji,jj,jk-1) - pub(ji,jj,jk) ) / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) * wumask(ji,jj,jk)62 & * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) & 63 & * ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) / ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) 63 64 zsh2v(ji,jj) = ( p_avm(ji,jj+1,jk) + p_avm(ji,jj,jk) ) & 64 & * ( pvn(ji,jj,jk-1) - pvn(ji,jj,jk) ) &65 & * ( pvb(ji,jj,jk-1) - pvb(ji,jj,jk) ) / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) * wvmask(ji,jj,jk)65 & * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) & 66 & * ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) 66 67 END DO 67 68 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdftke.F90
r10874 r10883 109 109 110 110 111 SUBROUTINE zdf_tke( kt, p_sh2, p_avm, p_avt )111 SUBROUTINE zdf_tke( kt, Kbb, Kmm, p_sh2, p_avm, p_avt ) 112 112 !!---------------------------------------------------------------------- 113 113 !! *** ROUTINE zdf_tke *** … … 155 155 !!---------------------------------------------------------------------- 156 156 INTEGER , INTENT(in ) :: kt ! ocean time step 157 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 157 158 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_sh2 ! shear production term 158 159 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: p_avm, p_avt ! momentum and tracer Kz (w-points) 159 160 !!---------------------------------------------------------------------- 160 161 ! 161 CALL tke_tke( gdepw_n, e3t_n, e3w_n, p_sh2, p_avm, p_avt ) ! now tke (en)162 ! 163 CALL tke_avn( gdepw_n, e3t_n, e3w_n, p_avm, p_avt ) ! now avt, avm, dissl162 CALL tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) ! now tke (en) 163 ! 164 CALL tke_avn( Kbb, Kmm, p_avm, p_avt ) ! now avt, avm, dissl 164 165 ! 165 166 END SUBROUTINE zdf_tke 166 167 167 168 168 SUBROUTINE tke_tke( pdepw, p_e3t, p_e3w, p_sh2, p_avm, p_avt )169 SUBROUTINE tke_tke( Kbb, Kmm, p_sh2, p_avm, p_avt ) 169 170 !!---------------------------------------------------------------------- 170 171 !! *** ROUTINE tke_tke *** … … 186 187 USE zdf_oce , ONLY : en ! ocean vertical physics 187 188 !! 188 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pdepw ! depth of w-points 189 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 189 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 190 190 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_sh2 ! shear production term 191 191 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) … … 243 243 zmskv = ( 2. - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) ) 244 244 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 245 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( u b(ji,jj,mbkt(ji,jj))+ub(ji-1,jj,mbkt(ji,jj)) ) )**2 &246 & + ( zmskv*( v b(ji,jj,mbkt(ji,jj))+vb(ji,jj-1,mbkt(ji,jj)) ) )**2 )245 zebot = - 0.001875_wp * rCdU_bot(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2 & 246 & + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2 ) 247 247 en(ji,jj,mbkt(ji,jj)+1) = MAX( zebot, rn_emin ) * ssmask(ji,jj) 248 248 END DO … … 254 254 zmskv = ( 2. - vmask(ji,jj-1,mikt(ji,jj)) * vmask(ji,jj,mikt(ji,jj)) ) 255 255 ! ! where 0.001875 = (rn_ebb0/rau0) * 0.5 = 3.75*0.5/1000. (CAUTION CdU<0) 256 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( u b(ji,jj,mikt(ji,jj))+ub(ji-1,jj,mikt(ji,jj)) ) )**2 &257 & + ( zmskv*( v b(ji,jj,mikt(ji,jj))+vb(ji,jj-1,mikt(ji,jj)) ) )**2 )256 zetop = - 0.001875_wp * rCdU_top(ji,jj) * SQRT( ( zmsku*( uu(ji,jj,mikt(ji,jj),Kbb)+uu(ji-1,jj,mikt(ji,jj),Kbb) ) )**2 & 257 & + ( zmskv*( vv(ji,jj,mikt(ji,jj),Kbb)+vv(ji,jj-1,mikt(ji,jj),Kbb) ) )**2 ) 258 258 en(ji,jj,mikt(ji,jj)) = MAX( zetop, rn_emin ) * (1._wp - tmask(ji,jj,1)) ! masked at ocean surface 259 259 END DO … … 268 268 ! 269 269 ! !* total energy produce by LC : cumulative sum over jk 270 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * pdepw(:,:,1) * p_e3w(:,:,1)270 zpelc(:,:,1) = MAX( rn2b(:,:,1), 0._wp ) * gdepw(:,:,1,Kmm) * e3w(:,:,1,Kmm) 271 271 DO jk = 2, jpk 272 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * pdepw(:,:,jk) * p_e3w(:,:,jk)272 zpelc(:,:,jk) = zpelc(:,:,jk-1) + MAX( rn2b(:,:,jk), 0._wp ) * gdepw(:,:,jk,Kmm) * e3w(:,:,jk,Kmm) 273 273 END DO 274 274 ! !* finite Langmuir Circulation depth … … 286 286 DO jj = 1, jpj 287 287 DO ji = 1, jpi 288 zhlc(ji,jj) = pdepw(ji,jj,imlc(ji,jj))288 zhlc(ji,jj) = gdepw(ji,jj,imlc(ji,jj),Kmm) 289 289 END DO 290 290 END DO … … 302 302 IF ( zfr_i(ji,jj) /= 0. ) THEN 303 303 ! vertical velocity due to LC 304 IF ( pdepw(ji,jj,jk) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN304 IF ( gdepw(ji,jj,jk,Kmm) - zhlc(ji,jj) < 0 .AND. wmask(ji,jj,jk) /= 0. ) THEN 305 305 ! ! vertical velocity due to LC 306 zwlc = rn_lc * SIN( rpi * pdepw(ji,jj,jk) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i306 zwlc = rn_lc * SIN( rpi * gdepw(ji,jj,jk,Kmm) / zhlc(ji,jj) ) ! warning: optimization: zus^3 is in zfr_i 307 307 ! ! TKE Langmuir circulation source term 308 308 en(ji,jj,jk) = en(ji,jj,jk) + rdt * zfr_i(ji,jj) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) … … 342 342 ! ! eddy coefficient (ensure numerical stability) 343 343 zzd_up = zcof * MAX( p_avm(ji,jj,jk+1) + p_avm(ji,jj,jk ) , 2.e-5_wp ) & ! upper diagonal 344 & / ( p_e3t(ji,jj,jk ) * p_e3w(ji,jj,jk) )344 & / ( e3t(ji,jj,jk ,Kmm) * e3w(ji,jj,jk ,Kmm) ) 345 345 zzd_lw = zcof * MAX( p_avm(ji,jj,jk ) + p_avm(ji,jj,jk-1) , 2.e-5_wp ) & ! lower diagonal 346 & / ( p_e3t(ji,jj,jk-1) * p_e3w(ji,jj,jk) )346 & / ( e3t(ji,jj,jk-1,Kmm) * e3w(ji,jj,jk ,Kmm) ) 347 347 ! 348 348 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) … … 402 402 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 403 403 !!gm BUG : in the exp remove the depth of ssh !!! 404 !!gm i.e. use gde3w in argument ( pdepw)404 !!gm i.e. use gde3w in argument (gdepw(:,:,:,Kmm)) 405 405 406 406 … … 409 409 DO jj = 2, jpjm1 410 410 DO ji = fs_2, fs_jpim1 ! vector opt. 411 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - pdepw(ji,jj,jk) / htau(ji,jj) ) &411 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 412 412 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 413 413 END DO … … 418 418 DO ji = fs_2, fs_jpim1 ! vector opt. 419 419 jk = nmln(ji,jj) 420 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( - pdepw(ji,jj,jk) / htau(ji,jj) ) &420 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 421 421 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 422 422 END DO … … 431 431 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 432 432 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 433 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( - pdepw(ji,jj,jk) / htau(ji,jj) ) &433 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -gdepw(ji,jj,jk,Kmm) / htau(ji,jj) ) & 434 434 & * MAX(0.,1._wp - rn_eice *fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 435 435 END DO … … 441 441 442 442 443 SUBROUTINE tke_avn( pdepw, p_e3t, p_e3w, p_avm, p_avt )443 SUBROUTINE tke_avn( Kbb, Kmm, p_avm, p_avt ) 444 444 !!---------------------------------------------------------------------- 445 445 !! *** ROUTINE tke_avn *** … … 477 477 USE zdf_oce , ONLY : en, avtb, avmb, avtb_2d ! ocean vertical physics 478 478 !! 479 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pdepw ! depth (w-points) 480 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: p_e3t, p_e3w ! level thickness (t- & w-points) 479 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 481 480 REAL(wp), DIMENSION(:,:,:), INTENT( out) :: p_avm, p_avt ! vertical eddy viscosity & diffusivity (w-points) 482 481 ! … … 526 525 ! 527 526 !!gm Not sure of that coding for ISF.... 528 ! where wmask = 0 set zmxlm == p_e3w527 ! where wmask = 0 set zmxlm == e3w(:,:,:,Kmm) 529 528 CASE ( 0 ) ! bounded by the distance to surface and bottom 530 529 DO jk = 2, jpkm1 531 530 DO jj = 2, jpjm1 532 531 DO ji = fs_2, fs_jpim1 ! vector opt. 533 zemxl = MIN( pdepw(ji,jj,jk) - pdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), &534 & pdepw(ji,jj,mbkt(ji,jj)+1) - pdepw(ji,jj,jk) )532 zemxl = MIN( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm), zmxlm(ji,jj,jk), & 533 & gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) - gdepw(ji,jj,jk,Kmm) ) 535 534 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 536 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))537 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , p_e3w(ji,jj,jk) ) * (1 - wmask(ji,jj,jk))535 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 536 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN( zmxlm(ji,jj,jk) , e3w(ji,jj,jk,Kmm) ) * (1 - wmask(ji,jj,jk)) 538 537 END DO 539 538 END DO … … 544 543 DO jj = 2, jpjm1 545 544 DO ji = fs_2, fs_jpim1 ! vector opt. 546 zemxl = MIN( p_e3w(ji,jj,jk), zmxlm(ji,jj,jk) )545 zemxl = MIN( e3w(ji,jj,jk,Kmm), zmxlm(ji,jj,jk) ) 547 546 zmxlm(ji,jj,jk) = zemxl 548 547 zmxld(ji,jj,jk) = zemxl … … 555 554 DO jj = 2, jpjm1 556 555 DO ji = fs_2, fs_jpim1 ! vector opt. 557 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )556 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 558 557 END DO 559 558 END DO … … 562 561 DO jj = 2, jpjm1 563 562 DO ji = fs_2, fs_jpim1 ! vector opt. 564 zemxl = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )563 zemxl = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 565 564 zmxlm(ji,jj,jk) = zemxl 566 565 zmxld(ji,jj,jk) = zemxl … … 573 572 DO jj = 2, jpjm1 574 573 DO ji = fs_2, fs_jpim1 ! vector opt. 575 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + p_e3t(ji,jj,jk-1), zmxlm(ji,jj,jk) )574 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + e3t(ji,jj,jk-1,Kmm), zmxlm(ji,jj,jk) ) 576 575 END DO 577 576 END DO … … 580 579 DO jj = 2, jpjm1 581 580 DO ji = fs_2, fs_jpim1 ! vector opt. 582 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + p_e3t(ji,jj,jk+1), zmxlm(ji,jj,jk) )581 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + e3t(ji,jj,jk+1,Kmm), zmxlm(ji,jj,jk) ) 583 582 END DO 584 583 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
r10876 r10883 448 448 449 449 ! ! Ocean physics 450 CALL zdf_phy_init ! Vertical physics450 CALL zdf_phy_init( Nnn ) ! Vertical physics 451 451 452 452 ! ! Lateral physics -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
r10880 r10883 137 137 138 138 ! VERTICAL PHYSICS 139 CALL zdf_phy( kstp )! vertical physics update (top/bot drag, avt, avs, avm + MLD)139 CALL zdf_phy( kstp, Nbb, Nnn ) ! vertical physics update (top/bot drag, avt, avs, avm + MLD) 140 140 141 141 ! LATERAL PHYSICS … … 196 196 CALL dyn_vor( kstp, Nnn , uu, vv, Nrhs ) ! vorticity ==> RHS 197 197 CALL dyn_ldf ( kstp ) ! lateral mixing 198 IF( ln_zdfosm ) CALL dyn_osm ( kstp ) ! OSMOSIS non-local velocity fluxes198 IF( ln_zdfosm ) CALL dyn_osm( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes ==> RHS 199 199 CALL dyn_hpg ( kstp ) ! horizontal gradient of Hydrostatic pressure 200 200 CALL dyn_spg ( kstp ) ! surface pressure gradient … … 253 253 #endif 254 254 CALL tra_adv( kstp, Nbb, Nnn , ts, Nrhs ) ! hor. + vert. advection ==> RHS 255 IF( ln_zdfosm ) CALL tra_osm ( kstp ) ! OSMOSIS non-local tracer fluxes255 IF( ln_zdfosm ) CALL tra_osm( kstp, Nnn , ts, Nrhs ) ! OSMOSIS non-local tracer fluxes ==> RHS 256 256 IF( lrst_oce .AND. ln_zdfosm ) & 257 & CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts257 & CALL osm_rst( kstp, Nnn, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts 258 258 CALL tra_ldf ( kstp ) ! lateral mixing 259 259
Note: See TracChangeset
for help on using the changeset viewer.
