- Timestamp:
- 2021-12-16T10:39:55+01:00 (3 years ago)
- Location:
- branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF
- Files:
-
- 4 edited
-
zdfbfr.F90 (modified) (8 diffs)
-
zdfevd.F90 (modified) (7 diffs)
-
zdfgls.F90 (modified) (40 diffs)
-
zdftke.F90 (modified) (43 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r12555 r15603 26 26 USE wrk_nemo ! Memory Allocation 27 27 USE phycst, ONLY: vkarmn 28 USE stopack 28 29 29 30 IMPLICIT NONE … … 52 53 REAL(wp), PUBLIC :: rn_tfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 53 54 LOGICAL , PUBLIC :: ln_bfrimp ! logical switch for implicit bottom friction 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d, tfrcoef2d ! 2D bottom/top drag coefficient (PUBLIC for TAM)55 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d, tfrcoef2d,bfrcoef2d0 ! 2D bottom/top drag coefficient (PUBLIC for TAM) 55 56 56 57 !! * Substitutions … … 68 69 !! *** FUNCTION zdf_bfr_alloc *** 69 70 !!---------------------------------------------------------------------- 70 ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )71 ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), bfrcoef2d0(jpi,jpj),STAT=zdf_bfr_alloc ) 71 72 ! 72 73 IF( lk_mpp ) CALL mpp_sum ( zdf_bfr_alloc ) … … 107 108 IF(lflush) CALL flush(numout) 108 109 ENDIF 110 ! 111 #if defined key_traldf_c2d || key_traldf_c3d 112 IF( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) THEN 113 bfrcoef2d = bfrcoef2d0 114 CALL spp_gen(kt, bfrcoef2d, nn_spp_bfr, rn_bfr_sd, jk_spp_bfr) 115 ENDIF 116 #else 117 IF ( ln_stopack .AND. ( nn_spp_bfr > 0 ) ) & 118 & CALL ctl_stop( 'zdf_bfr: parameter perturbation will only work with '// & 119 'key_traldf_c2d or key_traldf_c3d') 120 #endif 109 121 ! 110 122 IF( nn_bfr == 2 ) THEN ! quadratic bottom friction only … … 135 147 END DO 136 148 END IF 137 ! 149 ! 138 150 ELSE 139 151 zbfrt(:,:) = bfrcoef2d(:,:) … … 178 190 DO ji = 2, jpim1 179 191 ! (ISF) ======================================================================== 180 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 192 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 181 193 ikbv = mikv(ji,jj) ! (1st wet ocean u- and v-points) 182 194 ! … … 359 371 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 360 372 ENDIF 361 373 362 374 IF ( ln_isfcav ) THEN 363 375 IF(ln_tfr2d) THEN … … 494 506 ENDIF 495 507 ! 508 bfrcoef2d0(:,:) = bfrcoef2d(:,:) 496 509 IF( nn_timing == 1 ) CALL timing_stop('zdf_bfr_init') 497 510 ! -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfevd.F90
r12555 r15603 20 20 USE zdfkpp ! KPP vertical mixing 21 21 USE trd_oce ! trends: ocean variables 22 USE trdtra ! trends manager: tracers 22 USE trdtra ! trends manager: tracers 23 23 USE in_out_manager ! I/O manager 24 24 USE iom ! for iom_put 25 25 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 26 26 USE timing ! Timing 27 USE stopack 27 28 28 29 IMPLICIT NONE … … 30 31 31 32 PUBLIC zdf_evd ! called by step.F90 33 REAL(wp), SAVE, ALLOCATABLE, DIMENSION(:,:) :: rn_avevd0 32 34 33 35 !! * Substitutions … … 43 45 !!---------------------------------------------------------------------- 44 46 !! *** ROUTINE zdf_evd *** 45 !! 47 !! 46 48 !! ** Purpose : Local increased the vertical eddy viscosity and diffu- 47 49 !! sivity coefficients when a static instability is encountered. 48 50 !! 49 51 !! ** Method : avt, avm, and the 4 neighbouring avmu, avmv coefficients 50 !! are set to avevd (namelist parameter) if the water column is 52 !! are set to avevd (namelist parameter) if the water column is 51 53 !! statically unstable (i.e. if rn2 < -1.e-12 ) 52 54 !! … … 70 72 IF(lwp) WRITE(numout,*) 71 73 IF(lwp .AND. lflush) CALL flush(numout) 74 ALLOCATE ( rn_avevd0(jpi,jpj) ) 75 rn_avevd0(:,:) = rn_avevd 72 76 ENDIF 73 77 74 78 zavt_evd(:,:,:) = avt(:,:,:) ! set avt prior to evd application 79 80 #if defined key_traldf_c2d || key_traldf_c3d 81 IF( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) THEN 82 rn_avevd0(:,:) = rn_avevd 83 CALL spp_gen(kt, rn_avevd0, nn_spp_aevd, rn_aevd_sd, jk_spp_aevd) 84 ENDIF 85 #else 86 IF ( ln_stopack .AND. ( nn_spp_aevd > 0 ) ) & 87 & CALL ctl_stop( 'zdf_evd: parameter perturbation will only work with '// & 88 'key_traldf_c2d or key_traldf_c3d') 89 #endif 75 90 76 91 SELECT CASE ( nn_evdm ) … … 80 95 zavm_evd(:,:,:) = avm(:,:,:) ! set avm prior to evd application 81 96 ! 82 DO jk = 1, jpkm1 97 DO jk = 1, jpkm1 83 98 DO jj = 2, jpj ! no vector opt. 84 99 DO ji = 2, jpi … … 89 104 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) THEN 90 105 #endif 91 avt (ji ,jj ,jk) = rn_avevd * tmask(ji ,jj ,jk)92 avm (ji ,jj ,jk) = rn_avevd * tmask(ji ,jj ,jk)93 avmu(ji ,jj ,jk) = rn_avevd * umask(ji ,jj ,jk)94 avmu(ji-1,jj ,jk) = rn_avevd * umask(ji-1,jj ,jk)95 avmv(ji ,jj ,jk) = rn_avevd * vmask(ji ,jj ,jk)96 avmv(ji ,jj-1,jk) = rn_avevd * vmask(ji ,jj-1,jk)106 avt (ji ,jj ,jk) = rn_avevd0(ji,jj) * tmask(ji ,jj ,jk) 107 avm (ji ,jj ,jk) = rn_avevd0(ji,jj) * tmask(ji ,jj ,jk) 108 avmu(ji ,jj ,jk) = rn_avevd0(ji,jj) * umask(ji ,jj ,jk) 109 avmu(ji-1,jj ,jk) = rn_avevd0(ji,jj) * umask(ji-1,jj ,jk) 110 avmv(ji ,jj ,jk) = rn_avevd0(ji,jj) * vmask(ji ,jj ,jk) 111 avmv(ji ,jj-1,jk) = rn_avevd0(ji,jj) * vmask(ji ,jj-1,jk) 97 112 ENDIF 98 113 END DO 99 114 END DO 100 END DO 115 END DO 101 116 CALL lbc_lnk( avt , 'W', 1. ) ; CALL lbc_lnk( avm , 'W', 1. ) ! Lateral boundary conditions 102 117 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) … … 105 120 CALL iom_put( "avm_evd", zavm_evd ) ! output this change 106 121 ! 107 CASE DEFAULT ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 122 CASE DEFAULT ! enhance vertical eddy diffusivity only (if rn2<-1.e-12) 108 123 DO jk = 1, jpkm1 109 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 124 !!! WHERE( rn2(:,:,jk) <= -1.e-12 ) avt(:,:,jk) = tmask(:,:,jk) * avevd ! agissant sur T SEUL! 110 125 DO jj = 1, jpj ! loop over the whole domain (no lbc_lnk call) 111 126 DO ji = 1, jpi 112 127 #if defined key_zdfkpp 113 128 ! no evd mixing in the boundary layer with KPP 114 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 .AND. fsdepw(ji,jj,jk) > hkpp(ji,jj) ) & 129 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 .AND. fsdepw(ji,jj,jk) > hkpp(ji,jj) ) & 115 130 #else 116 131 IF( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) & 117 132 #endif 118 avt(ji,jj,jk) = rn_avevd * tmask(ji,jj,jk)133 avt(ji,jj,jk) = rn_avevd0(ji,jj) * tmask(ji,jj,jk) 119 134 END DO 120 135 END DO 121 136 END DO 122 137 ! 123 END SELECT 138 END SELECT 124 139 125 140 zavt_evd(:,:,:) = avt(:,:,:) - zavt_evd(:,:,:) ! change in avt due to evd -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r12555 r15603 2 2 !!====================================================================== 3 3 !! *** MODULE zdfgls *** 4 !! Ocean physics: vertical mixing coefficient computed from the gls 4 !! Ocean physics: vertical mixing coefficient computed from the gls 5 5 !! turbulent closure parameterization 6 6 !!====================================================================== … … 16 16 !! gls_rst : read/write gls restart in ocean restart file 17 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and active tracers 18 USE oce ! ocean dynamics and active tracers 19 19 USE dom_oce ! ocean space and time domain 20 20 USE domvvl ! ocean space and time domain : variable volume layer … … 31 31 USE iom ! I/O manager library 32 32 USE timing ! Timing 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 34 USE stopack 34 35 35 36 IMPLICIT NONE … … 61 62 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 62 63 REAL(wp) :: rn_hsro ! Minimum surface roughness 63 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 64 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 64 65 65 66 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 66 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 67 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 67 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 68 REAL(wp) :: rl_sf = 0.2_wp ! 0 <rl_sf<vkarmn 68 69 REAL(wp) :: rghmin = -0.28_wp 69 70 REAL(wp) :: rgh0 = 0.0329_wp … … 72 73 REAL(wp) :: ra2 = 0.74_wp 73 74 REAL(wp) :: rb1 = 16.60_wp 74 REAL(wp) :: rb2 = 10.10_wp 75 REAL(wp) :: re2 = 1.33_wp 75 REAL(wp) :: rb2 = 10.10_wp 76 REAL(wp) :: re2 = 1.33_wp 76 77 REAL(wp) :: rl1 = 0.107_wp 77 78 REAL(wp) :: rl2 = 0.0032_wp … … 133 134 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 134 135 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 135 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 137 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 137 138 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - … … 139 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 140 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 143 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 143 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before … … 156 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 157 158 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 158 159 159 160 ! Preliminary computing 160 161 … … 165 166 avm (:,:,:) = avm_k (:,:,:) 166 167 avmu(:,:,:) = avmu_k(:,:,:) 167 avmv(:,:,:) = avmv_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 168 169 ENDIF 169 170 170 171 ! Compute surface and bottom friction at T-points 171 !CDIR NOVERRCHK 172 DO jj = 2, jpjm1 173 !CDIR NOVERRCHK 174 DO ji = fs_2, fs_jpim1 ! vector opt. 172 !CDIR NOVERRCHK 173 DO jj = 2, jpjm1 174 !CDIR NOVERRCHK 175 DO ji = fs_2, fs_jpim1 ! vector opt. 175 176 ! 176 177 ! surface friction 177 178 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 178 ! 179 ! bottom friction (explicit before friction) 180 ! Note that we chose here not to bound the friction as in dynbfr) 181 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 182 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 183 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 184 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 185 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 186 END DO 187 END DO 179 ! 180 ! bottom friction (explicit before friction) 181 ! Note that we chose here not to bound the friction as in dynbfr) 182 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 183 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 184 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 185 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 186 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 187 END DO 188 END DO 188 189 189 190 ! Set surface roughness length 190 191 SELECT CASE ( nn_z0_met ) 191 192 ! 192 CASE ( 0 ) ! Constant roughness 193 CASE ( 0 ) ! Constant roughness 193 194 zhsro(:,:) = rn_hsro 194 195 CASE ( 1 ) ! Standard Charnock formula … … 226 227 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 227 228 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 229 DO jj = 2, jpjm1 229 230 DO ji = fs_2, fs_jpim1 ! vector opt. 230 231 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) … … 275 276 IF( ln_sigpsi ) THEN 276 277 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 277 zwall_psi(ji,jj,jk) = rsc_psi / & 278 zwall_psi(ji,jj,jk) = rsc_psi / & 278 279 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 279 280 ELSE … … 294 295 ! diagonal 295 296 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 296 & + rdt * zdiss * tmask(ji,jj,jk) 297 & + rdt * zdiss * tmask(ji,jj,jk) 297 298 ! 298 299 ! right hand side in en … … 316 317 ! First level 317 318 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 318 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 320 z_elem_a(:,:,1) = en(:,:,1) 320 321 z_elem_c(:,:,1) = 0._wp 321 322 z_elem_b(:,:,1) = 1._wp 322 ! 323 ! 323 324 ! One level below 324 325 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 325 326 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 326 327 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 327 z_elem_a(:,:,2) = 0._wp 328 z_elem_a(:,:,2) = 0._wp 328 329 z_elem_c(:,:,2) = 0._wp 329 330 z_elem_b(:,:,2) = 1._wp … … 334 335 ! Dirichlet conditions at k=1 335 336 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 336 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 338 z_elem_a(:,:,1) = en(:,:,1) 338 339 z_elem_c(:,:,1) = 0._wp … … 357 358 SELECT CASE ( nn_bc_bot ) 358 359 ! 359 CASE ( 0 ) ! Dirichlet 360 CASE ( 0 ) ! Dirichlet 360 361 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 361 362 ! ! Balance between the production and the dissipation terms … … 377 378 z_elem_c(ji,jj,ibotm1) = 0._wp 378 379 z_elem_b(ji,jj,ibotm1) = 1._wp 379 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 380 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 380 381 END DO 381 382 END DO 382 383 ! 383 384 CASE ( 1 ) ! Neumman boundary condition 384 ! 385 ! 385 386 !CDIR NOVERRCHK 386 387 DO jj = 2, jpjm1 … … 428 429 END DO 429 430 END DO 430 ! ! set the minimum value of tke 431 ! ! set the minimum value of tke 431 432 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 432 433 … … 470 471 DO jj = 2, jpjm1 471 472 DO ji = fs_2, fs_jpim1 ! vector opt. 472 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 473 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 473 474 END DO 474 475 END DO … … 489 490 ! 490 491 ! psi / k 491 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 492 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 492 493 ! 493 494 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) … … 509 510 zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term 510 511 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 511 ! 512 ! 512 513 ! building the matrix 513 514 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) … … 551 552 z_elem_c(:,:,2) = 0._wp 552 553 z_elem_b(:,:,2) = 1._wp 553 ! 554 ! 554 555 ! 555 556 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz … … 575 576 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 576 577 577 ! 578 ! 578 579 ! 579 580 END SELECT … … 585 586 ! 586 587 ! 587 CASE ( 0 ) ! Dirichlet 588 CASE ( 0 ) ! Dirichlet 588 589 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 589 590 ! ! Balance between the production and the dissipation terms … … 610 611 ! 611 612 CASE ( 1 ) ! Neumman boundary condition 612 ! 613 ! 613 614 !CDIR NOVERRCHK 614 615 DO jj = 2, jpjm1 … … 692 693 DO jj = 2, jpjm1 693 694 DO ji = fs_2, fs_jpim1 ! vector opt. 694 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 695 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 695 696 END DO 696 697 END DO … … 719 720 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 720 721 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 721 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 722 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 722 723 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 723 724 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 724 725 END DO 725 726 END DO 726 END DO 727 END DO 727 728 728 729 ! … … 806 807 END DO 807 808 END DO 809 #if defined key_traldf_c2d || key_traldf_c3d 810 IF( ln_stopack) THEN 811 IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk) 812 IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk) 813 ENDIF 814 #else 815 IF ( ln_stopack ) & 816 & CALL ctl_stop( 'zdf_gls: parameter perturbation will only work with '// & 817 'key_traldf_c2d or key_traldf_c3d') 818 #endif 808 819 END DO 809 820 ! … … 846 857 !!---------------------------------------------------------------------- 847 858 !! *** ROUTINE zdf_gls_init *** 848 !! 849 !! ** Purpose : Initialization of the vertical eddy diffivity and 859 !! 860 !! ** Purpose : Initialization of the vertical eddy diffivity and 850 861 !! viscosity when using a gls turbulent closure scheme 851 862 !! … … 881 892 READ ( numnam_cfg, namzdf_gls, IOSTAT = ios, ERR = 902 ) 882 893 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_gls in configuration namelist', lwp ) 883 IF(lwm .AND. nprint > 2) WRITE ( numond, namzdf_gls )894 IF(lwm) WRITE ( numond, namzdf_gls ) 884 895 885 896 IF(lwp) THEN !* Control print … … 1069 1080 END SELECT 1070 1081 ! 1071 IF(lwp .AND. lflush) CALL flush(numout) 1082 IF(lwp .AND. lflush) CALL flush(numout) 1072 1083 ! !* Set Schmidt number for psi diffusion in the wave breaking case 1073 1084 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1074 1085 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1075 1086 IF( ln_sigpsi ) THEN 1076 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1087 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1077 1088 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1078 1089 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work … … 1082 1093 rsc_psi0 = rsc_psi 1083 1094 ENDIF 1084 1095 1085 1096 ! !* Shear free turbulence parameters 1086 1097 ! … … 1129 1140 rc04 = rc03 * rc0 1130 1141 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1131 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1142 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1132 1143 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1133 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1144 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1134 1145 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1135 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1146 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1136 1147 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1137 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1148 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1138 1149 1139 1150 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke … … 1150 1161 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1151 1162 END DO 1152 ! 1163 ! 1153 1164 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1154 1165 ! … … 1161 1172 !!--------------------------------------------------------------------- 1162 1173 !! *** ROUTINE ts_rst *** 1163 !! 1174 !! 1164 1175 !! ** Purpose : Read or write TKE file (en) in restart file 1165 1176 !! 1166 1177 !! ** Method : use of IOM library 1167 !! if the restart does not contain TKE, en is either 1178 !! if the restart does not contain TKE, en is either 1168 1179 !! set to rn_emin or recomputed (nn_igls/=0) 1169 1180 !!---------------------------------------------------------------------- … … 1177 1188 !!---------------------------------------------------------------------- 1178 1189 ! 1179 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1190 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1180 1191 ! ! --------------- 1181 1192 IF( ln_rstart ) THEN !* Read the restart file … … 1196 1207 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1197 1208 IF(nn_timing == 2) CALL timing_stop('iom_rstget') 1198 ELSE 1209 ELSE 1199 1210 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1200 1211 IF(lwp .AND. lflush) CALL flush(numout) 1201 1212 en (:,:,:) = rn_emin 1202 mxln(:,:,:) = 0.05 1213 mxln(:,:,:) = 0.05 1203 1214 avt_k (:,:,:) = avt (:,:,:) 1204 1215 avm_k (:,:,:) = avm (:,:,:) … … 1211 1222 IF(lwp .AND. lflush) CALL flush(numout) 1212 1223 en (:,:,:) = rn_emin 1213 mxln(:,:,:) = 0.05 1224 mxln(:,:,:) = 0.05 1214 1225 ENDIF 1215 1226 ! … … 1219 1230 IF(lwp .AND. lflush) CALL flush(numout) 1220 1231 IF(nn_timing == 2) CALL timing_start('iom_rstput') 1221 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1232 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1222 1233 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1223 1234 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1224 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1235 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1225 1236 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1226 1237 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) … … 1241 1252 WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?' 1242 1253 END SUBROUTINE zdf_gls_init 1243 1254 1244 1255 SUBROUTINE zdf_gls( kt ) ! Empty routine 1245 1256 IMPLICIT NONE 1246 INTEGER, INTENT(in) :: kt ! ocean time step 1257 INTEGER, INTENT(in) :: kt ! ocean time step 1247 1258 WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 1248 1259 END SUBROUTINE zdf_gls 1249 1260 1250 1261 SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine 1251 1262 IMPLICIT NONE … … 1258 1269 !!====================================================================== 1259 1270 END MODULE zdfgls 1260 -
branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r12555 r15603 2 2 !!====================================================================== 3 3 !! *** MODULE zdftke *** 4 !! Ocean physics: vertical mixing coefficient computed from the tke 4 !! Ocean physics: vertical mixing coefficient computed from the tke 5 5 !! turbulent closure parameterization 6 6 !!===================================================================== … … 22 22 !! - ! 2008-05 (J.-M. Molines, G. Madec) 2D form of avtb 23 23 !! - ! 2008-06 (G. Madec) style + DOCTOR name for namelist parameters 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 24 !! - ! 2008-12 (G. Reffray) stable discretization of the production term 25 !! 3.2 ! 2009-06 (G. Madec, S. Masson) TKE restart compatible with key_cpl 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase … … 52 52 USE wrk_nemo ! work arrays 53 53 USE timing ! Timing 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 54 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 55 55 #if defined key_agrif 56 56 USE agrif_opa_interp 57 57 USE agrif_opa_update 58 58 #endif 59 60 59 USE stopack 61 60 62 61 IMPLICIT NONE … … 75 74 INTEGER :: nn_pdl ! Prandtl number or not (ratio avt/avm) (=0/1) 76 75 REAL(wp) :: rn_ediff ! coefficient for avt: avt=rn_ediff*mxl*sqrt(e) 77 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 76 REAL(wp) :: rn_ediss ! coefficient of the Kolmogoroff dissipation 78 77 REAL(wp) :: rn_ebb ! coefficient of the surface input of tke 79 78 REAL(wp) :: rn_emin ! minimum value of tke [m2/s2] … … 94 93 95 94 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: htau ! depth of tke penetration (nn_htau) 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term97 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 498 95 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dissl ! now mixing lenght of dissipation 99 96 #if defined key_c1d … … 102 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_pdl, e_ric !: prandl and local Richardson numbers 103 100 #endif 101 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: e_niw !: TKE budget- near-inertial waves term 103 REAL(wp) , ALLOCATABLE, SAVE, DIMENSION(:,:) :: efr ! surface boundary condition for nn_etau = 4 104 104 105 105 !! * Substitutions … … 118 118 !!---------------------------------------------------------------------- 119 119 ALLOCATE( & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 120 & efr (jpi,jpj) , e_niw(jpi,jpj,jpk) , & 121 121 #if defined key_c1d 122 122 & e_dis(jpi,jpj,jpk) , e_mix(jpi,jpj,jpk) , & … … 147 147 !! surface: en = max( rn_emin0, rn_ebb * taum ) 148 148 !! bottom : en = rn_emin 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 149 !! The associated critical Richardson number is: ri_cri = 2/(2+rn_ediss/rn_ediff) 150 !! 151 !! The now Turbulent kinetic energy is computed using the following 152 152 !! time stepping: implicit for vertical diffusion term, linearized semi 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 153 !! implicit for kolmogoroff dissipation term, and explicit forward for 154 !! both buoyancy and shear production terms. Therefore a tridiagonal 155 155 !! linear system is solved. Note that buoyancy and shear terms are 156 156 !! discretized in a energy conserving form (Bruchard 2002). … … 160 160 !! 161 161 !! The now vertical eddy vicosity and diffusivity coefficients are 162 !! given by: 162 !! given by: 163 163 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 164 !! avt = max( avmb, pdl * avm ) 164 !! avt = max( avmb, pdl * avm ) 165 165 !! eav = max( avmb, avm ) 166 166 !! where pdl, the inverse of the Prandtl number is 1 if nn_pdl=0 and 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 167 !! given by an empirical funtion of the localRichardson number if nn_pdl=1 168 168 !! 169 169 !! ** Action : compute en (now turbulent kinetic energy) … … 180 180 ! 181 181 IF( kt /= nit000 ) THEN ! restore before value to compute tke 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 182 avt (:,:,:) = avt_k (:,:,:) 183 avm (:,:,:) = avm_k (:,:,:) 184 avmu(:,:,:) = avmu_k(:,:,:) 185 avmv(:,:,:) = avmv_k(:,:,:) 186 ENDIF 187 ! 188 #if defined key_traldf_c2d || key_traldf_c3d 189 IF( ln_stopack) THEN 190 IF( nn_spp_tkelc > 0 ) THEN 191 rn_lc0 = rn_lc 192 CALL spp_gen( kt, rn_lc0, nn_spp_tkelc, rn_tkelc_sd, jk_spp_tkelc ) 193 ENDIF 194 IF( nn_spp_tkedf > 0 ) THEN 195 rn_ediff0 = rn_ediff 196 CALL spp_gen( kt, rn_ediff0, nn_spp_tkedf, rn_tkedf_sd, jk_spp_tkedf ) 197 ENDIF 198 IF( nn_spp_tkeds > 0 ) THEN 199 rn_ediss0 = rn_ediss 200 CALL spp_gen( kt, rn_ediss0, nn_spp_tkeds, rn_tkeds_sd, jk_spp_tkeds ) 201 ENDIF 202 IF( nn_spp_tkebb > 0 ) THEN 203 rn_ebb0 = rn_ebb 204 CALL spp_gen( kt, rn_ebb0, nn_spp_tkebb, rn_tkebb_sd, jk_spp_tkebb ) 205 ENDIF 206 IF( nn_spp_tkefr > 0 ) THEN 207 rn_efr0 = rn_efr 208 CALL spp_gen( kt, rn_efr0, nn_spp_tkefr, rn_tkefr_sd, jk_spp_tkefr ) 209 ENDIF 210 ENDIF 211 #else 212 IF ( ln_stopack ) & 213 & CALL ctl_stop( 'zdf_tke: parameter perturbation will only work with '// & 214 'key_traldf_c2d or key_traldf_c3d') 215 #endif 187 216 ! 188 217 CALL tke_tke ! now tke (en) … … 190 219 CALL tke_avn ! now avt, avm, avmu, avmv 191 220 ! 192 avt_k (:,:,:) = avt (:,:,:) 193 avm_k (:,:,:) = avm (:,:,:) 194 avmu_k(:,:,:) = avmu(:,:,:) 195 avmv_k(:,:,:) = avmv(:,:,:) 221 avt_k (:,:,:) = avt (:,:,:) 222 avm_k (:,:,:) = avm (:,:,:) 223 avmu_k(:,:,:) = avmu(:,:,:) 224 avmv_k(:,:,:) = avmv(:,:,:) 196 225 ! 197 226 #if defined key_agrif 198 ! Update child grid f => parent grid 227 ! Update child grid f => parent grid 199 228 IF( .NOT.Agrif_Root() ) CALL Agrif_Update_Tke( kt ) ! children only 200 #endif 201 ! 229 #endif 230 IF ( kt == nitend ) THEN 231 DEALLOCATE ( rn_lc0, rn_ediff0, rn_ediss0, rn_ebb0, rn_efr0 ) 232 ENDIF 233 ! 234 202 235 END SUBROUTINE zdf_tke 203 236 … … 225 258 REAL(wp) :: zrhoa = 1.22 ! Air density kg/m3 226 259 REAL(wp) :: zcdrag = 1.5e-3 ! drag coefficient 227 REAL(wp) :: z bbrau, zesh2! temporary scalars228 REAL(wp) :: zfact1 , zfact2, zfact3! - -260 REAL(wp) :: zesh2 ! temporary scalars 261 REAL(wp) :: zfact1 ! - - 229 262 REAL(wp) :: ztx2 , zty2 , zcof ! - - 230 263 REAL(wp) :: ztau , zdif ! - - … … 233 266 !!bfr REAL(wp) :: zebot ! - - 234 267 INTEGER , POINTER, DIMENSION(:,: ) :: imlc 235 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc 268 REAL(wp), POINTER, DIMENSION(:,: ) :: zhlc, zbbrau,zfact2,zfact3 236 269 REAL(wp), POINTER, DIMENSION(:,:,:) :: zpelc, zdiag, zd_up, zd_lw 237 270 !!-------------------------------------------------------------------- … … 240 273 ! 241 274 CALL wrk_alloc( jpi,jpj, imlc ) ! integer 242 CALL wrk_alloc( jpi,jpj, zhlc ) 243 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 244 ! 245 zbbrau = rn_ebb / rau0 ! Local constant initialisation 246 zfact1 = -.5_wp * rdt 247 zfact2 = 1.5_wp * rdt * rn_ediss 248 zfact3 = 0.5_wp * rn_ediss 275 CALL wrk_alloc( jpi,jpj, zhlc ) 276 CALL wrk_alloc( jpi,jpj, zbbrau ) 277 CALL wrk_alloc( jpi,jpj, zfact2 ) 278 CALL wrk_alloc( jpi,jpj, zfact3 ) 279 CALL wrk_alloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 280 ! 281 zbbrau = rn_ebb0 / rau0 ! Local constant initialisation 282 zfact1 = -.5_wp * rdt 283 zfact2 = 1.5_wp * rdt * rn_ediss0 284 zfact3 = 0.5_wp * rn_ediss0 249 285 ! 250 286 ! … … 261 297 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 262 298 DO ji = fs_2, fs_jpim1 ! vector opt. 263 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)264 END DO 265 END DO 266 299 en(ji,jj,1) = MAX( rn_emin0, zbbrau(ji,jj) * taum(ji,jj) ) * tmask(ji,jj,1) 300 END DO 301 END DO 302 267 303 !!bfr - start commented area 268 304 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 303 339 imlc(:,:) = mbkt(:,:) + 1 ! Initialization to the number of w ocean point (=2 over land) 304 340 DO jk = jpkm1, 2, -1 305 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 341 DO jj = 1, jpj ! Last w-level at which zpelc>=0.5*us*us 306 342 DO ji = 1, jpi ! with us=0.016*wind(starting from jpk-1) 307 343 zus = zcof * taum(ji,jj) … … 311 347 END DO 312 348 ! ! finite LC depth 313 DO jj = 1, jpj 349 DO jj = 1, jpj 314 350 DO ji = 1, jpi 315 351 zhlc(ji,jj) = fsdepw(ji,jj,imlc(ji,jj)) … … 326 362 ! ! vertical velocity due to LC 327 363 zind = 0.5 - SIGN( 0.5, fsdepw(ji,jj,jk) - zhlc(ji,jj) ) 328 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) )364 zwlc = zind * rn_lc0(ji,jj) * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 329 365 ! ! TKE Langmuir circulation source term 330 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) ) 366 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( 1._wp - fr_i(ji,jj) )* ( zwlc * zwlc * zwlc ) / & 331 367 & zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 368 332 369 END DO 333 370 END DO … … 347 384 DO ji = 1, jpi 348 385 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 349 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 386 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 350 387 & / ( fse3uw_n(ji,jj,jk) & 351 388 & * fse3uw_b(ji,jj,jk) ) … … 376 413 ! ! shear prod. at w-point weightened by mask 377 414 zesh2 = ( avmu(ji-1,jj,jk) + avmu(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) ) & 378 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 415 & + ( avmv(ji,jj-1,jk) + avmv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) ) 379 416 ! 380 417 zd_up(ji,jj,jk) = zzd_up ! Matrix (zdiag, zd_up, zd_lw) 381 418 zd_lw(ji,jj,jk) = zzd_lw 382 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2 * dissl(ji,jj,jk) * tmask(ji,jj,jk)419 zdiag(ji,jj,jk) = 1._wp - zzd_lw - zzd_up + zfact2(ji,jj) * dissl(ji,jj,jk) * tmask(ji,jj,jk) 383 420 ! 384 421 ! ! right hand side in en 385 422 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 386 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) &423 & + zfact3(ji,jj) * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 387 424 & * wmask(ji,jj,jk) 388 425 END DO … … 433 470 END DO 434 471 435 ! ! Save TKE prior to nn_etau addition 436 e_niw(:,:,:) = en(:,:,:) 437 ! 472 ! ! Save TKE prior to nn_etau addition 473 e_niw(:,:,:) = en(:,:,:) 474 ! 438 475 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 439 476 ! ! TKE due to surface and internal wave breaking … … 455 492 DO jj = 2, jpjm1 456 493 DO ji = fs_2, fs_jpim1 ! vector opt. 457 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &494 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 458 495 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 459 496 END DO … … 464 501 DO ji = fs_2, fs_jpim1 ! vector opt. 465 502 jk = nmln(ji,jj) 466 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &503 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr0(ji,jj) * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 467 504 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 468 505 END DO … … 477 514 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 478 515 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 479 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 480 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 516 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 517 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 481 518 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 482 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) &519 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau(ji,jj) * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 483 520 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 484 521 END DO … … 486 523 END DO 487 524 ELSEIF( nn_etau == 4 ) THEN !* column integral independant of htau (rn_efr must be scaled up) 488 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 525 IF( nn_htau == 2 ) THEN ! efr dependant on time-varying htau 489 526 DO jj = 2, jpjm1 490 527 DO ji = fs_2, fs_jpim1 ! vector opt. … … 504 541 CALL lbc_lnk( en, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 505 542 ! 506 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 507 DO jj = 2, jpjm1 508 DO ji = fs_2, fs_jpim1 ! vector opt. 509 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 510 END DO 511 END DO 512 END DO 513 ! 514 CALL lbc_lnk( e_niw, 'W', 1. ) 543 DO jk = 2, jpkm1 ! TKE budget: near-inertial waves term 544 DO jj = 2, jpjm1 545 DO ji = fs_2, fs_jpim1 ! vector opt. 546 e_niw(ji,jj,jk) = en(ji,jj,jk) - e_niw(ji,jj,jk) 547 END DO 548 END DO 549 END DO 550 ! 551 CALL lbc_lnk( e_niw, 'W', 1. ) 515 552 ! 516 553 CALL wrk_dealloc( jpi,jpj, imlc ) ! integer 517 CALL wrk_dealloc( jpi,jpj, zhlc ) 518 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 554 CALL wrk_dealloc( jpi,jpj, zhlc ) 555 CALL wrk_dealloc( jpi,jpj, zbbrau ) 556 CALL wrk_dealloc( jpi,jpj, zfact2 ) 557 CALL wrk_dealloc( jpi,jpj, zfact3 ) 558 CALL wrk_dealloc( jpi,jpj,jpk, zpelc, zdiag, zd_up, zd_lw ) 519 559 ! 520 560 IF( nn_timing == 1 ) CALL timing_stop('tke_tke') … … 529 569 !! ** Purpose : Compute the vertical eddy viscosity and diffusivity 530 570 !! 531 !! ** Method : At this stage, en, the now TKE, is known (computed in 532 !! the tke_tke routine). First, the now mixing lenth is 571 !! ** Method : At this stage, en, the now TKE, is known (computed in 572 !! the tke_tke routine). First, the now mixing lenth is 533 573 !! computed from en and the strafification (N^2), then the mixings 534 574 !! coefficients are computed. 535 575 !! - Mixing length : a first evaluation of the mixing lengh 536 576 !! scales is: 537 !! mxl = sqrt(2*en) / N 577 !! mxl = sqrt(2*en) / N 538 578 !! where N is the brunt-vaisala frequency, with a minimum value set 539 579 !! to rmxl_min (rn_mxl0) in the interior (surface) ocean. 540 !! The mixing and dissipative length scale are bound as follow : 580 !! The mixing and dissipative length scale are bound as follow : 541 581 !! nn_mxl=0 : mxl bounded by the distance to surface and bottom. 542 582 !! zmxld = zmxlm = mxl 543 583 !! nn_mxl=1 : mxl bounded by the e3w and zmxld = zmxlm = mxl 544 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 584 !! nn_mxl=2 : mxl bounded such that the vertical derivative of mxl is 545 585 !! less than 1 (|d/dz(mxl)|<1) and zmxld = zmxlm = mxl 546 586 !! nn_mxl=3 : mxl is bounded from the surface to the bottom usings 547 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 548 !! the surface to obtain ldown. the resulting length 587 !! |d/dz(xml)|<1 to obtain lup, and from the bottom to 588 !! the surface to obtain ldown. the resulting length 549 589 !! scales are: 550 !! zmxld = sqrt( lup * ldown ) 590 !! zmxld = sqrt( lup * ldown ) 551 591 !! zmxlm = min ( lup , ldown ) 552 592 !! - Vertical eddy viscosity and diffusivity: 553 593 !! avm = max( avtb, rn_ediff * zmxlm * en^1/2 ) 554 !! avt = max( avmb, pdlr * avm ) 594 !! avt = max( avmb, pdlr * avm ) 555 595 !! with pdlr=1 if nn_pdl=0, pdlr=1/pdl=F(Ri) otherwise. 556 596 !! … … 567 607 IF( nn_timing == 1 ) CALL timing_start('tke_avn') 568 608 569 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 609 CALL wrk_alloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 570 610 571 611 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< … … 576 616 ! 577 617 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 578 zmxlm(:,:,:) = rmxl_min 618 zmxlm(:,:,:) = rmxl_min 579 619 zmxld(:,:,:) = rmxl_min 580 620 ! … … 586 626 END DO 587 627 END DO 588 ELSE 628 ELSE 589 629 zmxlm(:,:,1) = rn_mxl0 590 630 ENDIF … … 604 644 ! !* Physical limits for the mixing length 605 645 ! 606 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 646 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 607 647 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 608 648 ! … … 698 738 DO ji = fs_2, fs_jpim1 ! vector opt. 699 739 zsqen = SQRT( en(ji,jj,jk) ) 700 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen740 zav = rn_ediff0(ji,jj) * zmxlm(ji,jj,jk) * zsqen 701 741 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 702 742 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) … … 704 744 END DO 705 745 END DO 746 #if defined key_traldf_c2d || key_traldf_c3d 747 IF( ln_stopack) THEN 748 IF(nn_spp_avt > 0 ) CALL spp_gen( 1, avt(:,:,jk), nn_spp_avt, rn_avt_sd, jk_spp_avt, jk) 749 IF(nn_spp_avm > 0 ) CALL spp_gen( 1, avm(:,:,jk), nn_spp_avm, rn_avm_sd, jk_spp_avm, jk) 750 ENDIF 751 #else 752 IF ( ln_stopack ) & 753 & CALL ctl_stop( 'tke_avn: parameter perturbation will only work with '// & 754 'key_traldf_c2d or key_traldf_c3d') 755 #endif 706 756 END DO 707 757 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) … … 749 799 ENDIF 750 800 ! 751 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 801 CALL wrk_dealloc( jpi,jpj,jpk, zmpdl, zmxlm, zmxld ) 752 802 ! 753 803 IF( nn_timing == 1 ) CALL timing_stop('tke_avn') … … 759 809 !!---------------------------------------------------------------------- 760 810 !! *** ROUTINE zdf_tke_init *** 761 !! 762 !! ** Purpose : Initialization of the vertical eddy diffivity and 811 !! 812 !! ** Purpose : Initialization of the vertical eddy diffivity and 763 813 !! viscosity when using a tke turbulent closure scheme 764 814 !! … … 776 826 & rn_emin0, rn_bshear, nn_mxl , ln_mxl0 , & 777 827 & rn_mxl0 , nn_pdl , ln_lc , rn_lc , & 778 & nn_etau , nn_htau , rn_efr , rn_c 828 & nn_etau , nn_htau , rn_efr , rn_c 779 829 !!---------------------------------------------------------------------- 780 830 … … 843 893 rn_mxl0 = rmxl_min 844 894 ENDIF 845 895 896 ALLOCATE( rn_lc0 (jpi,jpj) ) ; rn_lc0 = rn_lc 897 ALLOCATE( rn_ediff0(jpi,jpj) ) ; rn_ediff0 = rn_ediff 898 ALLOCATE( rn_ediss0(jpi,jpj) ) ; rn_ediss0 = rn_ediss 899 ALLOCATE( rn_ebb0 (jpi,jpj) ) ; rn_ebb0 = rn_ebb 900 ALLOCATE( rn_efr0 (jpi,jpj) ) ; rn_efr0 = rn_efr 901 846 902 IF( nn_etau == 2 ) THEN 847 903 ierr = zdf_mxl_alloc() … … 855 911 856 912 ! !* depth of penetration of surface tke 857 IF( nn_etau /= 0 ) THEN 913 IF( nn_etau /= 0 ) THEN 858 914 htau(:,:) = 0._wp 859 915 SELECT CASE( nn_htau ) ! Choice of the depth of penetration … … 861 917 htau(:,:) = 10._wp 862 918 CASE( 1 ) ! F(latitude) : 0.5m to 30m poleward of 40 degrees 863 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 919 htau(:,:) = MAX( 0.5_wp, MIN( 30._wp, 45._wp* ABS( SIN( rpi/180._wp * gphit(:,:) ) ) ) ) 864 920 CASE( 2 ) ! fraction of depth-integrated TKE within mixed-layer 865 921 rhtau = -1._wp / LOG( 1._wp - rn_c ) … … 904 960 END DO 905 961 dissl(:,:,:) = 1.e-12_wp 906 ! 962 ! 907 963 CALL tke_rst( nit000, 'READ' ) !* read or initialize all required files 908 964 ! … … 913 969 !!--------------------------------------------------------------------- 914 970 !! *** ROUTINE tke_rst *** 915 !! 971 !! 916 972 !! ** Purpose : Read or write TKE file (en) in restart file 917 973 !! 918 974 !! ** Method : use of IOM library 919 !! if the restart does not contain TKE, en is either 920 !! set to rn_emin or recomputed 975 !! if the restart does not contain TKE, en is either 976 !! set to rn_emin or recomputed 921 977 !!---------------------------------------------------------------------- 922 978 INTEGER , INTENT(in) :: kt ! ocean time-step … … 927 983 !!---------------------------------------------------------------------- 928 984 ! 929 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 985 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 930 986 ! ! --------------- 931 987 IF( ln_rstart ) THEN !* Read the restart file … … 956 1012 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 957 1013 ! 958 avt_k (:,:,:) = avt (:,:,:)959 avm_k (:,:,:) = avm (:,:,:)960 avmu_k(:,:,:) = avmu(:,:,:)961 avmv_k(:,:,:) = avmv(:,:,:)962 !963 1014 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 964 1015 ENDIF 965 1016 ELSE !* Start from rest 966 1017 en(:,:,:) = rn_emin * tmask(:,:,:) 967 DO jk = 1, jpk ! set the Kz to the background value968 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk)969 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk)970 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk)971 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk)972 END DO973 1018 ENDIF 974 1019 ! 1020 avt_k (:,:,:) = avt (:,:,:) 1021 avm_k (:,:,:) = avm (:,:,:) 1022 avmu_k(:,:,:) = avmu(:,:,:) 1023 avmv_k(:,:,:) = avmv(:,:,:) 1024 975 1025 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file 976 1026 ! ! -------------------
Note: See TracChangeset
for help on using the changeset viewer.
