- Timestamp:
- 2020-06-10T13:13:39+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_1d_bugfixes/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r11442 r13088 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 USE stopack 35 35 … … 62 62 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 63 63 REAL(wp) :: rn_hsro ! Minimum surface roughness 64 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) 65 65 66 66 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 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 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 69 69 REAL(wp) :: rghmin = -0.28_wp 70 70 REAL(wp) :: rgh0 = 0.0329_wp … … 73 73 REAL(wp) :: ra2 = 0.74_wp 74 74 REAL(wp) :: rb1 = 16.60_wp 75 REAL(wp) :: rb2 = 10.10_wp 76 REAL(wp) :: re2 = 1.33_wp 75 REAL(wp) :: rb2 = 10.10_wp 76 REAL(wp) :: re2 = 1.33_wp 77 77 REAL(wp) :: rl1 = 0.107_wp 78 78 REAL(wp) :: rl2 = 0.0032_wp … … 134 134 INTEGER :: ji, jj, jk, ibot, ibotm1, dir ! dummy loop arguments 135 135 REAL(wp) :: zesh2, zsigpsi, zcoef, zex1, zex2 ! local scalars 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 136 REAL(wp) :: ztx2, zty2, zup, zdown, zcof ! - - 137 137 REAL(wp) :: zratio, zrn2, zflxb, sh ! - - 138 138 REAL(wp) :: prod, buoy, diss, zdiss, sm ! - - … … 140 140 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 141 141 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 142 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 143 143 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) 144 144 REAL(wp), POINTER, DIMENSION(:,:,:) :: eb ! tke at time before … … 157 157 CALL wrk_alloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 158 158 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 159 159 160 160 ! Preliminary computing 161 161 … … 166 166 avm (:,:,:) = avm_k (:,:,:) 167 167 avmu(:,:,:) = avmu_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 168 avmv(:,:,:) = avmv_k(:,:,:) 169 169 ENDIF 170 170 171 171 ! Compute surface and bottom friction at T-points 172 !CDIR NOVERRCHK 173 DO jj = 2, jpjm1 174 !CDIR NOVERRCHK 175 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. 176 176 ! 177 177 ! surface friction 178 178 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 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 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 189 189 190 190 ! Set surface roughness length 191 191 SELECT CASE ( nn_z0_met ) 192 192 ! 193 CASE ( 0 ) ! Constant roughness 193 CASE ( 0 ) ! Constant roughness 194 194 zhsro(:,:) = rn_hsro 195 195 CASE ( 1 ) ! Standard Charnock formula … … 227 227 IF( nn_clos == 0 ) THEN ! Mellor-Yamada 228 228 DO jk = 2, jpkm1 229 DO jj = 2, jpjm1 229 DO jj = 2, jpjm1 230 230 DO ji = fs_2, fs_jpim1 ! vector opt. 231 231 zup = mxln(ji,jj,jk) * fsdepw(ji,jj,mbkt(ji,jj)+1) … … 276 276 IF( ln_sigpsi ) THEN 277 277 zsigpsi = MIN( 1._wp, zesh2 / eps(ji,jj,jk) ) ! 0. <= zsigpsi <= 1. 278 zwall_psi(ji,jj,jk) = rsc_psi / & 278 zwall_psi(ji,jj,jk) = rsc_psi / & 279 279 & ( zsigpsi * rsc_psi + (1._wp-zsigpsi) * rsc_psi0 / MAX( zwall(ji,jj,jk), 1._wp ) ) 280 280 ELSE … … 295 295 ! diagonal 296 296 z_elem_b(ji,jj,jk) = 1._wp - z_elem_a(ji,jj,jk) - z_elem_c(ji,jj,jk) & 297 & + rdt * zdiss * tmask(ji,jj,jk) 297 & + rdt * zdiss * tmask(ji,jj,jk) 298 298 ! 299 299 ! right hand side in en … … 317 317 ! First level 318 318 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 319 en(:,:,1) = MAX(en(:,:,1), rn_emin) 319 en(:,:,1) = MAX(en(:,:,1), rn_emin) 320 320 z_elem_a(:,:,1) = en(:,:,1) 321 321 z_elem_c(:,:,1) = 0._wp 322 322 z_elem_b(:,:,1) = 1._wp 323 ! 323 ! 324 324 ! One level below 325 325 en(:,:,2) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1 * ((zhsro(:,:)+fsdepw(:,:,2)) & 326 326 & / zhsro(:,:) )**(1.5_wp*ra_sf))**(2._wp/3._wp) 327 327 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 328 z_elem_a(:,:,2) = 0._wp 328 z_elem_a(:,:,2) = 0._wp 329 329 z_elem_c(:,:,2) = 0._wp 330 330 z_elem_b(:,:,2) = 1._wp … … 335 335 ! Dirichlet conditions at k=1 336 336 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 337 en(:,:,1) = MAX(en(:,:,1), rn_emin) 337 en(:,:,1) = MAX(en(:,:,1), rn_emin) 338 338 z_elem_a(:,:,1) = en(:,:,1) 339 339 z_elem_c(:,:,1) = 0._wp … … 358 358 SELECT CASE ( nn_bc_bot ) 359 359 ! 360 CASE ( 0 ) ! Dirichlet 360 CASE ( 0 ) ! Dirichlet 361 361 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = rn_lmin 362 362 ! ! Balance between the production and the dissipation terms … … 378 378 z_elem_c(ji,jj,ibotm1) = 0._wp 379 379 z_elem_b(ji,jj,ibotm1) = 1._wp 380 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 380 en(ji,jj,ibotm1) = MAX( rc02r * ustarb2(ji,jj), rn_emin ) 381 381 END DO 382 382 END DO 383 383 ! 384 384 CASE ( 1 ) ! Neumman boundary condition 385 ! 385 ! 386 386 !CDIR NOVERRCHK 387 387 DO jj = 2, jpjm1 … … 429 429 END DO 430 430 END DO 431 ! ! set the minimum value of tke 431 ! ! set the minimum value of tke 432 432 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 433 433 … … 471 471 DO jj = 2, jpjm1 472 472 DO ji = fs_2, fs_jpim1 ! vector opt. 473 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 474 474 END DO 475 475 END DO … … 490 490 ! 491 491 ! psi / k 492 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 492 zratio = psi(ji,jj,jk) / eb(ji,jj,jk) 493 493 ! 494 494 ! psi3+ : stable : B=-KhN²<0 => N²>0 if rn2>0 dir = 1 (stable) otherwise dir = 0 (unstable) … … 510 510 zesh2 = dir * ( prod + buoy ) + (1._wp - dir ) * prod ! production term 511 511 zdiss = dir * ( diss / psi(ji,jj,jk) ) + (1._wp - dir ) * (diss-buoy) / psi(ji,jj,jk) ! dissipation term 512 ! 512 ! 513 513 ! building the matrix 514 514 zcof = rfact_psi * zwall_psi(ji,jj,jk) * tmask(ji,jj,jk) … … 552 552 z_elem_c(:,:,2) = 0._wp 553 553 z_elem_b(:,:,2) = 1._wp 554 ! 554 ! 555 555 ! 556 556 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz … … 576 576 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 577 577 578 ! 578 ! 579 579 ! 580 580 END SELECT … … 586 586 ! 587 587 ! 588 CASE ( 0 ) ! Dirichlet 588 CASE ( 0 ) ! Dirichlet 589 589 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 590 590 ! ! Balance between the production and the dissipation terms … … 611 611 ! 612 612 CASE ( 1 ) ! Neumman boundary condition 613 ! 613 ! 614 614 !CDIR NOVERRCHK 615 615 DO jj = 2, jpjm1 … … 693 693 DO jj = 2, jpjm1 694 694 DO ji = fs_2, fs_jpim1 ! vector opt. 695 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) 696 696 END DO 697 697 END DO … … 720 720 eps(ji,jj,jk) = MAX( eps(ji,jj,jk), rn_epsmin ) 721 721 mxln(ji,jj,jk) = rc03 * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / eps(ji,jj,jk) 722 ! 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) 723 723 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 724 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) ) 725 725 END DO 726 726 END DO 727 END DO 727 END DO 728 728 729 729 ! … … 807 807 END DO 808 808 END DO 809 #if defined key_traldf_c2d || key_traldf_c3d 809 810 IF( ln_stopack) THEN 810 811 IF( nn_spp_avt > 0 ) CALL spp_gen(kt, avt(:,:,jk), nn_spp_avt, rn_avt_sd,jk) 811 812 IF( nn_spp_avm > 0 ) CALL spp_gen(kt, avm(:,:,jk), nn_spp_avm, rn_avm_sd,jk) 812 ENDIF 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 813 819 END DO 814 820 ! … … 851 857 !!---------------------------------------------------------------------- 852 858 !! *** ROUTINE zdf_gls_init *** 853 !! 854 !! ** Purpose : Initialization of the vertical eddy diffivity and 859 !! 860 !! ** Purpose : Initialization of the vertical eddy diffivity and 855 861 !! viscosity when using a gls turbulent closure scheme 856 862 !! … … 1071 1077 ! 1072 1078 END SELECT 1073 1079 1074 1080 ! !* Set Schmidt number for psi diffusion in the wave breaking case 1075 1081 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1076 1082 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1077 1083 IF( ln_sigpsi ) THEN 1078 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1084 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1079 1085 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1080 1086 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work … … 1084 1090 rsc_psi0 = rsc_psi 1085 1091 ENDIF 1086 1092 1087 1093 ! !* Shear free turbulence parameters 1088 1094 ! … … 1130 1136 rc04 = rc03 * rc0 1131 1137 rsbc_tke1 = -3._wp/2._wp*rn_crban*ra_sf*rl_sf ! Dirichlet + Wave breaking 1132 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1138 rsbc_tke2 = rdt * rn_crban / rl_sf ! Neumann + Wave breaking 1133 1139 zcr = MAX(rsmall, rsbc_tke1**(1./(-ra_sf*3._wp/2._wp))-1._wp ) 1134 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1140 rtrans = 0.2_wp / zcr ! Ad. inverse transition length between log and wave layer 1135 1141 rsbc_zs1 = rn_charn/grav ! Charnock formula for surface roughness 1136 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1142 rsbc_zs2 = rn_frac_hs / 0.85_wp / grav * 665._wp ! Rascle formula for surface roughness 1137 1143 rsbc_psi1 = -0.5_wp * rdt * rc0**(rpp-2._wp*rmm) / rsc_psi 1138 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1144 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1139 1145 1140 1146 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke … … 1151 1157 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk) 1152 1158 END DO 1153 ! 1159 ! 1154 1160 CALL gls_rst( nit000, 'READ' ) !* read or initialize all required files 1155 1161 ! … … 1162 1168 !!--------------------------------------------------------------------- 1163 1169 !! *** ROUTINE ts_rst *** 1164 !! 1170 !! 1165 1171 !! ** Purpose : Read or write TKE file (en) in restart file 1166 1172 !! 1167 1173 !! ** Method : use of IOM library 1168 !! if the restart does not contain TKE, en is either 1174 !! if the restart does not contain TKE, en is either 1169 1175 !! set to rn_emin or recomputed (nn_igls/=0) 1170 1176 !!---------------------------------------------------------------------- … … 1178 1184 !!---------------------------------------------------------------------- 1179 1185 ! 1180 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1186 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1181 1187 ! ! --------------- 1182 1188 IF( ln_rstart ) THEN !* Read the restart file … … 1197 1203 CALL iom_get( numror, jpdom_autoglo, 'mxln' , mxln ) 1198 1204 IF(nn_timing == 2) CALL timing_stop('iom_rstget') 1199 ELSE 1205 ELSE 1200 1206 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1201 1207 en (:,:,:) = rn_emin 1202 mxln(:,:,:) = 0.05 1208 mxln(:,:,:) = 0.05 1203 1209 avt_k (:,:,:) = avt (:,:,:) 1204 1210 avm_k (:,:,:) = avm (:,:,:) … … 1210 1216 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1211 1217 en (:,:,:) = rn_emin 1212 mxln(:,:,:) = 0.05 1218 mxln(:,:,:) = 0.05 1213 1219 ENDIF 1214 1220 ! … … 1217 1223 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1218 1224 IF(nn_timing == 2) CALL timing_start('iom_rstput') 1219 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1225 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1220 1226 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1221 1227 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1222 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1228 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1223 1229 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1224 1230 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) … … 1239 1245 WRITE(*,*) 'zdf_gls_init: You should not have seen this print! error?' 1240 1246 END SUBROUTINE zdf_gls_init 1241 1247 1242 1248 SUBROUTINE zdf_gls( kt ) ! Empty routine 1243 1249 IMPLICIT NONE 1244 INTEGER, INTENT(in) :: kt ! ocean time step 1250 INTEGER, INTENT(in) :: kt ! ocean time step 1245 1251 WRITE(*,*) 'zdf_gls: You should not have seen this print! error?', kt 1246 1252 END SUBROUTINE zdf_gls 1247 1253 1248 1254 SUBROUTINE gls_rst( kt, cdrw ) ! Empty routine 1249 1255 IMPLICIT NONE … … 1256 1262 !!====================================================================== 1257 1263 END MODULE zdfgls 1258
Note: See TracChangeset
for help on using the changeset viewer.
