Changeset 15603 for branches/UKMO/dev_r5518_GO6_starthour_obsoper/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
- Timestamp:
- 2021-12-16T10:39:55+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
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
Note: See TracChangeset
for help on using the changeset viewer.
