Changeset 15797
- Timestamp:
- 2022-04-25T11:39:18+02:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/CO6_shelfclimate/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r7566 r15797 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 34 35 35 IMPLICIT NONE … … 61 61 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 62 62 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) 63 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 64 64 65 65 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 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 68 68 REAL(wp) :: rghmin = -0.28_wp 69 69 REAL(wp) :: rgh0 = 0.0329_wp … … 72 72 REAL(wp) :: ra2 = 0.74_wp 73 73 REAL(wp) :: rb1 = 16.60_wp 74 REAL(wp) :: rb2 = 10.10_wp 75 REAL(wp) :: re2 = 1.33_wp 74 REAL(wp) :: rb2 = 10.10_wp 75 REAL(wp) :: re2 = 1.33_wp 76 76 REAL(wp) :: rl1 = 0.107_wp 77 77 REAL(wp) :: rl2 = 0.0032_wp … … 133 133 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 134 134 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 135 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 135 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 136 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 137 137 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - … … 139 139 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 140 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 143 143 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before … … 156 156 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 157 157 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 158 158 159 159 ! Preliminary computing 160 160 … … 165 165 avm (:,:,:) = avm_k (:,:,:) 166 166 avmu(:,:,:) = avmu_k(:,:,:) 167 avmv(:,:,:) = avmv_k(:,:,:) 167 avmv(:,:,:) = avmv_k(:,:,:) 168 168 ENDIF 169 169 170 170 ! Compute surface and bottom friction at T-points 171 !CDIR NOVERRCHK 172 DO jj = 2, jpjm1 173 !CDIR NOVERRCHK 171 !CDIR NOVERRCHK 172 DO jj = 2, jpjm1 173 !CDIR NOVERRCHK 174 174 DO ji = fs_2, fs_jpim1 ! vector opt. 175 175 ! 176 176 ! surface friction 177 177 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 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 188 188 189 189 ! Set surface roughness length 190 190 SELECT CASE ( nn_z0_met ) 191 191 ! 192 CASE ( 0 ) ! Constant roughness 192 CASE ( 0 ) ! Constant roughness 193 193 zhsro(:,:) = rn_hsro 194 194 CASE ( 1 ) ! Standard Charnock formula … … 226 226 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 227 227 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 228 DO jj = 2, jpjm1 229 229 DO ji = fs_2, fs_jpim1 ! vector opt. 230 230 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) … … 275 275 IF( ln_sigpsi ) THEN 276 276 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 277 zwall_psi(ji,jj,jk) = rsc_psi / & 277 zwall_psi(ji,jj,jk) = rsc_psi / & 278 278 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 279 279 ELSE … … 294 294 ! diagonal 295 295 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) 296 & + rdt * zdiss * tmask(ji,jj,jk) 297 297 ! 298 298 ! right hand side in en … … 316 316 ! First level 317 317 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 318 en(:,:,1) = MAX(en(:,:,1), rn_emin) 318 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 319 z_elem_a(:,:,1) = en(:,:,1) 320 320 z_elem_c(:,:,1) = 0._wp 321 321 z_elem_b(:,:,1) = 1._wp 322 ! 322 ! 323 323 ! One level below 324 324 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 325 325 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 326 326 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 327 z_elem_a(:,:,2) = 0._wp 327 z_elem_a(:,:,2) = 0._wp 328 328 z_elem_c(:,:,2) = 0._wp 329 329 z_elem_b(:,:,2) = 1._wp … … 334 334 ! Dirichlet conditions at k=1 335 335 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 336 en(:,:,1) = MAX(en(:,:,1), rn_emin) 336 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 337 z_elem_a(:,:,1) = en(:,:,1) 338 338 z_elem_c(:,:,1) = 0._wp … … 357 357 SELECT CASE ( nn_bc_bot ) 358 358 ! 359 CASE ( 0 ) ! Dirichlet 359 CASE ( 0 ) ! Dirichlet 360 360 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 361 361 ! ! Balance between the production and the dissipation terms … … 377 377 z_elem_c(ji,jj,ibotm1) = 0._wp 378 378 z_elem_b(ji,jj,ibotm1) = 1._wp 379 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 379 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 380 380 END DO 381 381 END DO 382 382 ! 383 383 CASE ( 1 ) ! Neumman boundary condition 384 ! 384 ! 385 385 !CDIR NOVERRCHK 386 386 DO jj = 2, jpjm1 … … 428 428 END DO 429 429 END DO 430 ! ! set the minimum value of tke 430 ! ! set the minimum value of tke 431 431 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 432 432 … … 470 470 DO jj = 2, jpjm1 471 471 DO ji = fs_2, fs_jpim1 ! vector opt. 472 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 472 psi(ji,jj,jk) = rc02 * eb(ji,jj,jk) * mxlb(ji,jj,jk)**rnn 473 473 END DO 474 474 END DO … … 489 489 ! 490 490 ! psi / k 491 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 491 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 492 492 ! 493 493 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) … … 509 509 zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term 510 510 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 511 ! 511 ! 512 512 ! building the matrix 513 513 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) … … 551 551 z_elem_c(:,:,2) = 0._wp 552 552 z_elem_b(:,:,2) = 1._wp 553 ! 553 ! 554 554 ! 555 555 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz … … 575 575 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 576 576 577 ! 577 ! 578 578 ! 579 579 END SELECT … … 585 585 ! 586 586 ! 587 CASE ( 0 ) ! Dirichlet 587 CASE ( 0 ) ! Dirichlet 588 588 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 589 589 ! ! Balance between the production and the dissipation terms … … 610 610 ! 611 611 CASE ( 1 ) ! Neumman boundary condition 612 ! 612 ! 613 613 !CDIR NOVERRCHK 614 614 DO jj = 2, jpjm1 … … 692 692 DO jj = 2, jpjm1 693 693 DO ji = fs_2, fs_jpim1 ! vector opt. 694 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 694 eps(ji,jj,jk) = rc04 * en(ji,jj,jk) * psi(ji,jj,jk) 695 695 END DO 696 696 END DO … … 719 719 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 720 720 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) 721 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 722 722 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 723 723 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 724 END DO 725 725 END DO 726 END DO 726 END DO 727 727 728 728 ! … … 846 846 !!---------------------------------------------------------------------- 847 847 !! *** ROUTINE zdf_gls_init *** 848 !! 849 !! ** Purpose : Initialization of the vertical eddy diffivity and 848 !! 849 !! ** Purpose : Initialization of the vertical eddy diffivity and 850 850 !! viscosity when using a gls turbulent closure scheme 851 851 !! … … 1066 1066 ! 1067 1067 END SELECT 1068 1068 1069 1069 ! !* Set Schmidt number for psi diffusion in the wave breaking case 1070 1070 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1071 1071 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1072 1072 IF( ln_sigpsi ) THEN 1073 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1073 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1074 1074 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1075 1075 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work … … 1079 1079 rsc_psi0 = rsc_psi 1080 1080 ENDIF 1081 1081 1082 1082 ! !* Shear free turbulence parameters 1083 1083 ! … … 1125 1125 rc04 = rc03 * rc0 1126 1126 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1127 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1127 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1128 1128 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1129 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1129 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1130 1130 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1131 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1131 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1132 1132 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1133 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1133 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1134 1134 1135 1135 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke … … 1146 1146 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1147 1147 END DO 1148 ! 1148 ! 1149 1149 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1150 1150 ! … … 1157 1157 !!--------------------------------------------------------------------- 1158 1158 !! *** ROUTINE ts_rst *** 1159 !! 1159 !! 1160 1160 !! ** Purpose : Read or write TKE file (en) in restart file 1161 1161 !! 1162 1162 !! ** Method : use of IOM library 1163 !! if the restart does not contain TKE, en is either 1163 !! if the restart does not contain TKE, en is either 1164 1164 !! set to rn_emin or recomputed (nn_igls/=0) 1165 1165 !!---------------------------------------------------------------------- … … 1173 1173 !!---------------------------------------------------------------------- 1174 1174 ! 1175 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1175 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1176 1176 ! ! --------------- 1177 1177 IF( ln_rstart ) THEN !* Read the restart file … … 1190 1190 CALL iom_get( numror, jpdom_autoglo, 'avmv' , avmv ) 1191 1191 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1192 ELSE 1192 ELSE 1193 1193 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1194 1194 en (:,:,:) = rn_emin 1195 mxln(:,:,:) = 0.05 1195 mxln(:,:,:) = 0.05 1196 1196 avt_k (:,:,:) = avt (:,:,:) 1197 1197 avm_k (:,:,:) = avm (:,:,:) … … 1203 1203 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1204 1204 en (:,:,:) = rn_emin 1205 mxln(:,:,:) = 0.05 1205 mxln(:,:,:) = 0.05 1206 1206 ENDIF 1207 1207 ! … … 1209 1209 ! ! ------------------- 1210 1210 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1211 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1211 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1212 1212 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1213 1213 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1214 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1214 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1215 1215 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1216 1216 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) … … 1241 1241 !!====================================================================== 1242 1242 END MODULE zdfgls 1243
Note: See TracChangeset
for help on using the changeset viewer.
