- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r5109 r6808 8 8 !! 3.3 ! 2010-10 (C. Bricaud) Add in the reference 9 9 !!---------------------------------------------------------------------- 10 #if defined key_zdfgls || defined key_esopa10 #if defined key_zdfgls 11 11 !!---------------------------------------------------------------------- 12 12 !! 'key_zdfgls' Generic Length Scale vertical physics … … 42 42 LOGICAL , PUBLIC, PARAMETER :: lk_zdfgls = .TRUE. !: TKE vertical mixing flag 43 43 ! 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: en !: now turbulent kinetic energy45 44 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: mxln !: now mixing length 46 45 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: zwall !: wall function 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avt_k ! not enhanced Kz48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avm_k ! not enhanced Kz49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmu_k ! not enhanced Kz50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: avmv_k ! not enhanced Kz51 46 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustars2 !: Squared surface velocity scale at T-points 52 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ustarb2 !: Squared bottom velocity scale at T-points … … 107 102 108 103 !! * Substitutions 109 # include "domzgr_substitute.h90"110 104 # include "vectopt_loop_substitute.h90" 111 105 !!---------------------------------------------------------------------- … … 120 114 !! *** FUNCTION zdf_gls_alloc *** 121 115 !!---------------------------------------------------------------------- 122 ALLOCATE( en(jpi,jpj,jpk), mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 123 & avt_k (jpi,jpj,jpk) , avm_k (jpi,jpj,jpk), & 124 & avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), & 125 & ustars2(jpi,jpj), ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 116 ALLOCATE( mxln(jpi,jpj,jpk), zwall(jpi,jpj,jpk) , & 117 & ustars2(jpi,jpj) , ustarb2(jpi,jpj) , STAT= zdf_gls_alloc ) 126 118 ! 127 119 IF( lk_mpp ) CALL mpp_sum ( zdf_gls_alloc ) … … 161 153 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 162 154 ! 163 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro )164 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )155 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 156 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 165 157 166 158 ! Preliminary computing … … 176 168 177 169 ! Compute surface and bottom friction at T-points 178 !CDIR NOVERRCHK179 170 DO jj = 2, jpjm1 180 !CDIR NOVERRCHK181 171 DO ji = fs_2, fs_jpim1 ! vector opt. 182 172 ! … … 213 203 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 214 204 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 215 & / ( fse3uw_n(ji,jj,jk) & 216 & * fse3uw_b(ji,jj,jk) ) 205 & / ( e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) ) 217 206 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 218 207 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & 219 & / ( fse3vw_n(ji,jj,jk) & 220 & * fse3vw_b(ji,jj,jk) ) 208 & / ( e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) ) 221 209 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT(en(ji,jj,jk)) / mxln(ji,jj,jk) 222 210 END DO … … 235 223 DO jj = 2, jpjm1 236 224 DO ji = fs_2, fs_jpim1 ! vector opt. 237 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1)238 zdown = vkarmn * fsdepw(ji,jj,jk) * ( -fsdepw(ji,jj,jk) + fsdepw(ji,jj,mbkt(ji,jj)+1) )225 zup = mxln(ji,jj,jk) * gdepw_n(ji,jj,mbkt(ji,jj)+1) 226 zdown = vkarmn * gdepw_n(ji,jj,jk) * ( -gdepw_n(ji,jj,jk) + gdepw_n(ji,jj,mbkt(ji,jj)+1) ) 239 227 zcoef = ( zup / MAX( zdown, rsmall ) ) 240 228 zwall (ji,jj,jk) = ( 1._wp + re2 * zcoef*zcoef ) * tmask(ji,jj,jk) … … 293 281 ! lower diagonal 294 282 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 295 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )283 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 296 284 ! 297 285 ! upper diagonal 298 286 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 299 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )287 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 300 288 ! 301 289 ! diagonal … … 329 317 ! 330 318 ! One level below 331 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 319 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+gdepw_n(:,:,2)) & 320 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 332 321 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 333 322 z_elem_a(:,:,2) = 0._wp … … 349 338 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 350 339 z_elem_a(:,:,2) = 0._wp 351 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:)) )) 352 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) * ((zhsro(:,:)+fsdept(:,:,1))/zhsro(:,:) )**(1.5_wp*ra_sf) 353 354 en(:,:,2) = en(:,:,2) + zflxs(:,:)/fse3w(:,:,2) 340 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1.-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:)) )) 341 zflxs(:,:) = rsbc_tke2 * ustars2(:,:)**1.5_wp * zkar(:,:) & 342 & * ((zhsro(:,:)+gdept_n(:,:,1)) / zhsro(:,:) )**(1.5_wp*ra_sf) 343 344 en(:,:,2) = en(:,:,2) + zflxs(:,:)/e3w_n(:,:,2) 355 345 ! 356 346 ! … … 365 355 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 366 356 ! ! Balance between the production and the dissipation terms 367 !CDIR NOVERRCHK 368 DO jj = 2, jpjm1 369 !CDIR NOVERRCHK 357 DO jj = 2, jpjm1 370 358 DO ji = fs_2, fs_jpim1 ! vector opt. 371 359 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 388 376 CASE ( 1 ) ! Neumman boundary condition 389 377 ! 390 !CDIR NOVERRCHK 391 DO jj = 2, jpjm1 392 !CDIR NOVERRCHK 378 DO jj = 2, jpjm1 393 379 DO ji = fs_2, fs_jpim1 ! vector opt. 394 380 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 519 505 ! lower diagonal 520 506 z_elem_a(ji,jj,jk) = zcof * ( avm (ji,jj,jk ) + avm (ji,jj,jk-1) ) & 521 & / ( fse3t(ji,jj,jk-1) * fse3w(ji,jj,jk ) )507 & / ( e3t_n(ji,jj,jk-1) * e3w_n(ji,jj,jk ) ) 522 508 ! upper diagonal 523 509 z_elem_c(ji,jj,jk) = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & 524 & / ( fse3t(ji,jj,jk ) * fse3w(ji,jj,jk) )510 & / ( e3t_n(ji,jj,jk ) * e3w_n(ji,jj,jk) ) 525 511 ! diagonal 526 512 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & … … 550 536 ! 551 537 ! One level below 552 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans* fsdepw(:,:,2)/zhsro(:,:) )))553 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:)538 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdepw_n(:,:,2)/zhsro(:,:) ))) 539 zdep(:,:) = (zhsro(:,:) + gdepw_n(:,:,2)) * zkar(:,:) 554 540 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 555 541 z_elem_a(:,:,2) = 0._wp … … 572 558 ! 573 559 ! Set psi vertical flux at the surface: 574 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans* fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope575 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf)560 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*gdept_n(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 561 zdep(:,:) = ((zhsro(:,:) + gdept_n(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 576 562 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 577 563 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 578 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.)564 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + gdept_n(:,:,1))**(rnn-1.) 579 565 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 580 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2)566 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / e3w_n(:,:,2) 581 567 582 568 ! … … 593 579 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 594 580 ! ! Balance between the production and the dissipation terms 595 !CDIR NOVERRCHK 596 DO jj = 2, jpjm1 597 !CDIR NOVERRCHK 581 DO jj = 2, jpjm1 598 582 DO ji = fs_2, fs_jpim1 ! vector opt. 599 583 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 606 590 ! 607 591 ! Just above last level, Dirichlet condition again (GOTM like) 608 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) )592 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + e3t_n(ji,jj,ibotm1) ) 609 593 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 610 594 z_elem_a(ji,jj,ibotm1) = 0._wp … … 616 600 CASE ( 1 ) ! Neumman boundary condition 617 601 ! 618 !CDIR NOVERRCHK 619 DO jj = 2, jpjm1 620 !CDIR NOVERRCHK 602 DO jj = 2, jpjm1 621 603 DO ji = fs_2, fs_jpim1 ! vector opt. 622 604 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point … … 636 618 ! 637 619 ! Set psi vertical flux at the bottom: 638 zdep(ji,jj) = rn_bfrz0 + 0.5_wp* fse3t(ji,jj,ibotm1)620 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*e3t_n(ji,jj,ibotm1) 639 621 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 640 622 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) 641 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / fse3w(ji,jj,ibotm1)623 psi(ji,jj,ibotm1) = psi(ji,jj,ibotm1) + zflxb / e3w_n(ji,jj,ibotm1) 642 624 END DO 643 625 END DO … … 839 821 avmv_k(:,:,:) = avmv(:,:,:) 840 822 ! 841 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro )842 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )823 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 824 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 843 825 ! 844 826 IF( nn_timing == 1 ) CALL timing_stop('zdf_gls')
Note: See TracChangeset
for help on using the changeset viewer.
