- Timestamp:
- 2015-07-09T12:14:37+02:00 (9 years ago)
- Location:
- branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF
- Files:
-
- 6 edited
-
zdfbfr.F90 (modified) (9 diffs)
-
zdfddm.F90 (modified) (3 diffs)
-
zdfgls.F90 (modified) (31 diffs)
-
zdfini.F90 (modified) (3 diffs)
-
zdftke.F90 (modified) (23 diffs)
-
zdftmx.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r5477 r5572 120 120 zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 121 121 zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 122 ! (ISF)123 ikbt = mikt(ji,jj)124 ! JC: possible WAD implementation should modify line below if layers vanish125 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp126 ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp)127 ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max)128 129 122 END DO 130 123 END DO 124 ! (ISF) 125 IF ( ln_isfcav ) THEN 126 DO jj = 1, jpj 127 DO ji = 1, jpi 128 ikbt = mikt(ji,jj) 129 ! JC: possible WAD implementation should modify line below if layers vanish 130 ztmp = (1-tmask(ji,jj,1)) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 131 ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 132 ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 133 END DO 134 END DO 135 END IF 131 136 ! 132 137 ELSE … … 152 157 ! 153 158 ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 154 IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 155 bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) & 156 & + ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) ) & 157 & * zecu * (1._wp - umask(ji,jj,1)) 158 END IF 159 IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 160 bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) & 161 & + ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) ) & 162 & * zecv * (1._wp - vmask(ji,jj,1)) 163 END IF 164 ! (ISF) ======================================================================== 165 ikbu = miku(ji,jj) ! ocean bottom level at u- and v-points 166 ikbv = mikv(ji,jj) ! (deepest ocean u- and v-points) 167 ! 168 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & 169 & + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) ) 170 zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) & 171 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 172 ! 173 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 174 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 175 ! 176 tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) * zecu * (1._wp - umask(ji,jj,1)) 177 tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 178 ! (ISF) END ==================================================================== 179 ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 180 IF ( miku(ji,jj) + 2 .GE. mbku(ji,jj) ) THEN 181 tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) & 182 & + ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) ) & 183 & * zecu * (1._wp - umask(ji,jj,1)) 184 END IF 185 IF ( mikv(ji,jj) + 2 .GE. mbkv(ji,jj) ) THEN 186 tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) & 187 & + ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) ) & 188 & * zecv * (1._wp - vmask(ji,jj,1)) 159 IF ( ln_isfcav ) THEN 160 IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 161 bfrua(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) & 162 & + ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) ) & 163 & * zecu * (1._wp - umask(ji,jj,1)) 164 END IF 165 IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 166 bfrva(ji,jj) = - 0.5_wp * ( ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) & 167 & + ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) ) & 168 & * zecv * (1._wp - vmask(ji,jj,1)) 169 END IF 189 170 END IF 190 171 END DO 191 172 END DO 192 !193 173 CALL lbc_lnk( bfrua, 'U', 1. ) ; CALL lbc_lnk( bfrva, 'V', 1. ) ! Lateral boundary condition 174 175 IF ( ln_isfcav ) THEN 176 DO jj = 2, jpjm1 177 DO ji = 2, jpim1 178 ! (ISF) ======================================================================== 179 ikbu = miku(ji,jj) ! ocean top level at u- and v-points 180 ikbv = mikv(ji,jj) ! (1st wet ocean u- and v-points) 181 ! 182 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & 183 & + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) ) 184 zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) & 185 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 186 ! 187 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_tfeb2 ) 188 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_tfeb2 ) 189 ! 190 tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) * zecu * (1._wp - umask(ji,jj,1)) 191 tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) * zecv * (1._wp - vmask(ji,jj,1)) 192 ! (ISF) END ==================================================================== 193 ! in case of 2 cell water column, we assume each cell feels the top and bottom friction 194 IF ( miku(ji,jj) + 1 .GE. mbku(ji,jj) ) THEN 195 tfrua(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) & 196 & + ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) ) & 197 & * zecu * (1._wp - umask(ji,jj,1)) 198 END IF 199 IF ( mikv(ji,jj) + 1 .GE. mbkv(ji,jj) ) THEN 200 tfrva(ji,jj) = - 0.5_wp * ( ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) & 201 & + ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) ) & 202 & * zecv * (1._wp - vmask(ji,jj,1)) 203 END IF 204 END DO 205 END DO 206 CALL lbc_lnk( tfrua, 'U', 1. ) ; CALL lbc_lnk( tfrva, 'V', 1. ) ! Lateral boundary condition 207 END IF 208 ! 194 209 ! 195 210 IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, & … … 264 279 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 265 280 ENDIF 266 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_bfri1 267 IF( ln_tfr2d ) THEN 268 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 269 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 270 ENDIF 281 IF ( ln_isfcav ) THEN 282 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_tfri1 283 IF( ln_tfr2d ) THEN 284 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 285 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 286 ENDIF 287 END IF 271 288 ! 272 289 IF(ln_bfr2d) THEN … … 282 299 bfrua(:,:) = - bfrcoef2d(:,:) 283 300 bfrva(:,:) = - bfrcoef2d(:,:) 301 ! 302 IF ( ln_isfcav ) THEN 303 IF(ln_tfr2d) THEN 304 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 305 CALL iom_open('tfr_coef.nc',inum) 306 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 307 CALL iom_close(inum) 308 tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 309 ELSE 310 tfrcoef2d(:,:) = rn_tfri1 ! initialize tfrcoef2d to the namelist variable 311 ENDIF 312 ! 313 tfrua(:,:) = - tfrcoef2d(:,:) 314 tfrva(:,:) = - tfrcoef2d(:,:) 315 END IF 284 316 ! 285 317 CASE( 2 ) … … 298 330 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 299 331 ENDIF 300 IF(lwp) WRITE(numout,*) ' quadratic top friction' 301 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_tfri2 302 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 303 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 304 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 305 IF(lwp) WRITE(numout,*) ' bottom roughness rn_tfrz0 [m] = ', rn_tfrz0 306 IF( rn_tfrz0<=0.e0 ) THEN 307 WRITE(ctmp1,*) ' bottom roughness must be strictly positive' 308 CALL ctl_stop( ctmp1 ) 309 ENDIF 310 IF( ln_tfr2d ) THEN 311 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 312 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 313 ENDIF 332 IF ( ln_isfcav ) THEN 333 IF(lwp) WRITE(numout,*) ' quadratic top friction' 334 IF(lwp) WRITE(numout,*) ' friction coef. rn_tfri2 = ', rn_tfri2 335 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 336 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 337 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 338 IF(lwp) WRITE(numout,*) ' top roughness rn_tfrz0 [m] = ', rn_tfrz0 339 IF( rn_tfrz0<=0.e0 ) THEN 340 WRITE(ctmp1,*) ' top roughness must be strictly positive' 341 CALL ctl_stop( ctmp1 ) 342 ENDIF 343 IF( ln_tfr2d ) THEN 344 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 345 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 346 ENDIF 347 END IF 314 348 ! 315 349 IF(ln_bfr2d) THEN … … 323 357 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 324 358 ENDIF 359 360 IF ( ln_isfcav ) THEN 361 IF(ln_tfr2d) THEN 362 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 363 CALL iom_open('tfr_coef.nc',inum) 364 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 365 CALL iom_close(inum) 366 ! 367 tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 368 ELSE 369 tfrcoef2d(:,:) = rn_tfri2 ! initialize tfrcoef2d to the namelist variable 370 ENDIF 371 END IF 325 372 ! 326 373 IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all … … 333 380 END DO 334 381 END DO 382 IF ( ln_isfcav ) THEN 383 DO jj = 1, jpj 384 DO ji = 1, jpi 385 ikbt = mikt(ji,jj) 386 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 387 tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 388 tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 389 END DO 390 END DO 391 END IF 335 392 ENDIF 336 393 ! … … 385 442 zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) ) 386 443 zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) ) 444 ! (ISF) 445 IF ( ln_isfcav ) THEN 446 ikbu = miku(ji,jj) ! 1st wet ocean level at u- and v-points 447 ikbv = mikv(ji,jj) 448 zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 449 zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 450 IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 451 IF( ln_ctl ) THEN 452 WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbu 453 WRITE(numout,*) 'TFR ', ABS( tfrcoef2d(ji,jj) ), zfru 454 ENDIF 455 ictu = ictu + 1 456 ENDIF 457 IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 458 IF( ln_ctl ) THEN 459 WRITE(numout,*) 'TFR ', narea, nimpp+ji, njmpp+jj, ikbv 460 WRITE(numout,*) 'TFR ', tfrcoef2d(ji,jj), zfrv 461 ENDIF 462 ictv = ictv + 1 463 ENDIF 464 zmintfr = MIN( zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) ) ) 465 zmaxtfr = MAX( zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) ) ) 466 END IF 467 ! END ISF 387 468 END DO 388 469 END DO … … 392 473 CALL mpp_min( zminbfr ) 393 474 CALL mpp_max( zmaxbfr ) 475 IF ( ln_isfcav) CALL mpp_min( zmintfr ) 476 IF ( ln_isfcav) CALL mpp_max( zmaxtfr ) 394 477 ENDIF 395 478 IF( .NOT.ln_bfrimp) THEN 396 479 IF( lwp .AND. ictu + ictv > 0 ) THEN 397 WRITE(numout,*) ' Bottom friction stability check failed at ', ictu, ' U-points '398 WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points '480 WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictu, ' U-points ' 481 WRITE(numout,*) ' Bottom/Top friction stability check failed at ', ictv, ' V-points ' 399 482 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 400 WRITE(numout,*) ' Bottomfriction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr401 WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary'483 IF ( ln_isfcav ) WRITE(numout,*) ' Top friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 484 WRITE(numout,*) ' Bottom/Top friction coefficient will be reduced where necessary' 402 485 ENDIF 403 486 ENDIF -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r5477 r5572 156 156 END DO 157 157 ! mask zmsk in order to have avt and avs masked 158 zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)158 zmsks(:,:) = zmsks(:,:) * wmask(:,:,jk) 159 159 160 160 … … 191 191 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), & 192 192 & avt(ji,jj,jk), avt(ji+1,jj,jk), & 193 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * umask(ji,jj,jk)193 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * wumask(ji,jj,jk) 194 194 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), & 195 195 & avt(ji,jj,jk), avt(ji,jj+1,jk), & 196 & avs(ji,jj,jk), avs(ji,jj+1,jk) ) * vmask(ji,jj,jk)196 & avs(ji,jj,jk), avs(ji,jj+1,jk) ) * wvmask(ji,jj,jk) 197 197 END DO 198 198 END DO … … 255 255 IF( zdf_ddm_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'zdf_ddm_init : unable to allocate arrays' ) 256 256 ! ! initialization to masked Kz 257 avs(:,:,:) = rn_avt0 * tmask(:,:,:)257 avs(:,:,:) = rn_avt0 * wmask(:,:,:) 258 258 ! 259 259 END SUBROUTINE zdf_ddm_init -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfgls.F90
r5477 r5572 20 20 USE domvvl ! ocean space and time domain : variable volume layer 21 21 USE zdf_oce ! ocean vertical physics 22 USE zdfbfr ! bottom friction (only for rn_bfrz0) 22 23 USE sbc_oce ! surface boundary condition: ocean 23 24 USE phycst ! physical constants … … 52 53 53 54 ! !! ** Namelist namzdf_gls ** 54 LOGICAL :: ln_crban ! =T use Craig and Banner scheme55 55 LOGICAL :: ln_length_lim ! use limit on the dissipation rate under stable stratification (Galperin et al. 1988) 56 56 LOGICAL :: ln_sigpsi ! Activate Burchard (2003) modification for k-eps closure & wave breaking mixing 57 INTEGER :: nn_tkebc_surf ! TKE surface boundary condition (=0/1) 58 INTEGER :: nn_tkebc_bot ! TKE bottom boundary condition (=0/1) 59 INTEGER :: nn_psibc_surf ! PSI surface boundary condition (=0/1) 60 INTEGER :: nn_psibc_bot ! PSI bottom boundary condition (=0/1) 57 INTEGER :: nn_bc_surf ! surface boundary condition (=0/1) 58 INTEGER :: nn_bc_bot ! bottom boundary condition (=0/1) 59 INTEGER :: nn_z0_met ! Method for surface roughness computation 61 60 INTEGER :: nn_stab_func ! stability functions G88, KC or Canuto (=0/1/2) 62 61 INTEGER :: nn_clos ! closure 0/1/2/3 MY82/k-eps/k-w/gen … … 66 65 REAL(wp) :: rn_charn ! Charnock constant for surface breaking waves mixing : 1400. (standard) or 2.e5 (Stacey value) 67 66 REAL(wp) :: rn_crban ! Craig and Banner constant for surface breaking waves mixing 68 69 REAL(wp) :: hsro = 0.003_wp ! Minimum surface roughness70 REAL(wp) :: hbro = 0.003_wp ! Bottom roughness (m) 67 REAL(wp) :: rn_hsro ! Minimum surface roughness 68 REAL(wp) :: rn_frac_hs ! Fraction of wave height as surface roughness (if nn_z0_met > 1) 69 71 70 REAL(wp) :: rcm_sf = 0.73_wp ! Shear free turbulence parameters 72 71 REAL(wp) :: ra_sf = -2.0_wp ! Must be negative -2 < ra_sf < -1 … … 96 95 REAL(wp) :: rm7 = 0.0_wp 97 96 REAL(wp) :: rm8 = 0.318_wp 98 97 REAL(wp) :: rtrans = 0.1_wp 99 98 REAL(wp) :: rc02, rc02r, rc03, rc04 ! coefficients deduced from above parameters 100 REAL(wp) :: rc03_sqrt2_galp ! - - - - 101 REAL(wp) :: rsbc_tke1, rsbc_tke2, rsbc_tke3, rfact_tke ! - - - - 102 REAL(wp) :: rsbc_psi1, rsbc_psi2, rsbc_psi3, rfact_psi ! - - - - 103 REAL(wp) :: rsbc_mb , rsbc_std , rsbc_zs ! - - - - 99 REAL(wp) :: rsbc_tke1, rsbc_tke2, rfact_tke ! - - - - 100 REAL(wp) :: rsbc_psi1, rsbc_psi2, rfact_psi ! - - - - 101 REAL(wp) :: rsbc_zs1, rsbc_zs2 ! - - - - 104 102 REAL(wp) :: rc0, rc2, rc3, rf6, rcff, rc_diff ! - - - - 105 103 REAL(wp) :: rs0, rs1, rs2, rs4, rs5, rs6 ! - - - - … … 147 145 REAL(wp) :: gh, gm, shr, dif, zsqen, zav ! - - 148 146 REAL(wp), POINTER, DIMENSION(:,: ) :: zdep 147 REAL(wp), POINTER, DIMENSION(:,: ) :: zkar 149 148 REAL(wp), POINTER, DIMENSION(:,: ) :: zflxs ! Turbulence fluxed induced by internal waves 150 149 REAL(wp), POINTER, DIMENSION(:,: ) :: zhsro ! Surface roughness (surface waves) … … 153 152 REAL(wp), POINTER, DIMENSION(:,:,:) :: shear ! vertical shear 154 153 REAL(wp), POINTER, DIMENSION(:,:,:) :: eps ! dissipation rate 155 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi.AND.ln_crban=T) 156 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a, z_elem_b, z_elem_c, psi 154 REAL(wp), POINTER, DIMENSION(:,:,:) :: zwall_psi ! Wall function use in the wb case (ln_sigpsi) 155 REAL(wp), POINTER, DIMENSION(:,:,:) :: psi ! psi at time now 156 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_a ! element of the first matrix diagonal 157 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_b ! element of the second matrix diagonal 158 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_elem_c ! element of the third matrix diagonal 157 159 !!-------------------------------------------------------------------- 158 160 ! 159 161 IF( nn_timing == 1 ) CALL timing_start('zdf_gls') 160 162 ! 161 CALL wrk_alloc( jpi,jpj, zdep, z flxs, zhsro )162 CALL wrk_alloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi )163 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 ) 165 164 166 ! Preliminary computing 165 167 … … 174 176 175 177 ! Compute surface and bottom friction at T-points 176 !CDIR NOVERRCHK 177 DO jj = 2, jpjm1 178 !CDIR NOVERRCHK 179 DO ji = fs_2, fs_jpim1 ! vector opt. 180 ! 181 ! surface friction 178 !CDIR NOVERRCHK 179 DO jj = 2, jpjm1 180 !CDIR NOVERRCHK 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 ! 183 ! surface friction 182 184 ustars2(ji,jj) = r1_rau0 * taum(ji,jj) * tmask(ji,jj,1) 183 ! 184 ! bottom friction (explicit before friction) 185 ! Note that we chose here not to bound the friction as in dynbfr) 186 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 187 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 188 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 189 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 190 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 191 END DO 192 END DO 193 194 ! In case of breaking surface waves mixing, 195 ! Compute surface roughness length according to Charnock formula: 196 IF( ln_crban ) THEN ; zhsro(:,:) = MAX(rsbc_zs * ustars2(:,:), hsro) 197 ELSE ; zhsro(:,:) = hsro 198 ENDIF 185 ! 186 ! bottom friction (explicit before friction) 187 ! Note that we chose here not to bound the friction as in dynbfr) 188 ztx2 = ( bfrua(ji,jj) * ub(ji,jj,mbku(ji,jj)) + bfrua(ji-1,jj) * ub(ji-1,jj,mbku(ji-1,jj)) ) & 189 & * ( 1._wp - 0.5_wp * umask(ji,jj,1) * umask(ji-1,jj,1) ) 190 zty2 = ( bfrva(ji,jj) * vb(ji,jj,mbkv(ji,jj)) + bfrva(ji,jj-1) * vb(ji,jj-1,mbkv(ji,jj-1)) ) & 191 & * ( 1._wp - 0.5_wp * vmask(ji,jj,1) * vmask(ji,jj-1,1) ) 192 ustarb2(ji,jj) = SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) 193 END DO 194 END DO 195 196 ! Set surface roughness length 197 SELECT CASE ( nn_z0_met ) 198 ! 199 CASE ( 0 ) ! Constant roughness 200 zhsro(:,:) = rn_hsro 201 CASE ( 1 ) ! Standard Charnock formula 202 zhsro(:,:) = MAX(rsbc_zs1 * ustars2(:,:), rn_hsro) 203 CASE ( 2 ) ! Roughness formulae according to Rascle et al., Ocean Modelling (2008) 204 zdep(:,:) = 30.*TANH(2.*0.3/(28.*SQRT(MAX(ustars2(:,:),rsmall)))) ! Wave age (eq. 10) 205 zhsro(:,:) = MAX(rsbc_zs2 * ustars2(:,:) * zdep(:,:)**1.5, rn_hsro) ! zhsro = rn_frac_hs * Hsw (eq. 11) 206 ! 207 END SELECT 199 208 200 209 ! Compute shear and dissipation rate … … 303 312 ! 304 313 ! Set surface condition on zwall_psi (1 at the bottom) 305 IF( ln_sigpsi ) THEN 306 zcoef = rsc_psi / rsc_psi0 307 DO jj = 2, jpjm1 308 DO ji = fs_2, fs_jpim1 ! vector opt. 309 zwall_psi(ji,jj,1) = zcoef 310 END DO 311 END DO 312 ENDIF 313 314 zwall_psi(:,:,1) = zwall_psi(:,:,2) 315 zwall_psi(:,:,jpk) = 1. 316 ! 314 317 ! Surface boundary condition on tke 315 318 ! --------------------------------- 316 319 ! 317 SELECT CASE ( nn_ tkebc_surf )320 SELECT CASE ( nn_bc_surf ) 318 321 ! 319 322 CASE ( 0 ) ! Dirichlet case 320 ! 321 IF (ln_crban) THEN ! Wave induced mixing case 322 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 323 ! ! balance between the production and the dissipation terms including the wave effect 324 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 325 z_elem_a(:,:,1) = en(:,:,1) 326 z_elem_c(:,:,1) = 0._wp 327 z_elem_b(:,:,1) = 1._wp 328 ! 329 ! one level below 330 en(:,:,2) = MAX( rsbc_tke1 * ustars2(:,:) * ( (zhsro(:,:)+fsdepw(:,:,2))/zhsro(:,:) )**ra_sf, rn_emin ) 331 z_elem_a(:,:,2) = 0._wp 332 z_elem_c(:,:,2) = 0._wp 333 z_elem_b(:,:,2) = 1._wp 334 ! 335 ELSE ! No wave induced mixing case 336 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 337 ! ! balance between the production and the dissipation terms 338 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 339 z_elem_a(:,:,1) = en(:,:,1) 340 z_elem_c(:,:,1) = 0._wp 341 z_elem_b(:,:,1) = 1._wp 342 ! 343 ! one level below 344 en(:,:,2) = MAX( rc02r * ustars2(:,:), rn_emin ) 345 z_elem_a(:,:,2) = 0._wp 346 z_elem_c(:,:,2) = 0._wp 347 z_elem_b(:,:,2) = 1._wp 348 ! 349 ENDIF 350 ! 323 ! First level 324 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 325 en(:,:,1) = MAX(en(:,:,1), rn_emin) 326 z_elem_a(:,:,1) = en(:,:,1) 327 z_elem_c(:,:,1) = 0._wp 328 z_elem_b(:,:,1) = 1._wp 329 ! 330 ! 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) 332 en(:,:,2) = MAX(en(:,:,2), rn_emin ) 333 z_elem_a(:,:,2) = 0._wp 334 z_elem_c(:,:,2) = 0._wp 335 z_elem_b(:,:,2) = 1._wp 336 ! 337 ! 351 338 CASE ( 1 ) ! Neumann boundary condition on d(e)/dz 352 ! 353 IF (ln_crban) THEN ! Shear free case: d(e)/dz= Fw 354 ! 355 ! Dirichlet conditions at k=1 (Cosmetic) 356 en(:,:,1) = MAX( rsbc_tke1 * ustars2(:,:), rn_emin ) 357 z_elem_a(:,:,1) = en(:,:,1) 358 z_elem_c(:,:,1) = 0._wp 359 z_elem_b(:,:,1) = 1._wp 360 ! at k=2, set de/dz=Fw 361 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 362 z_elem_a(:,:,2) = 0._wp 363 zflxs(:,:) = rsbc_tke3 * ustars2(:,:)**1.5_wp * ( (zhsro(:,:)+fsdept(:,:,1) ) / zhsro(:,:) )**(1.5*ra_sf) 364 en(:,:,2) = en(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 365 ! 366 ELSE ! No wave induced mixing case: d(e)/dz=0. 367 ! 368 ! Dirichlet conditions at k=1 (Cosmetic) 369 en(:,:,1) = MAX( rc02r * ustars2(:,:), rn_emin ) 370 z_elem_a(:,:,1) = en(:,:,1) 371 z_elem_c(:,:,1) = 0._wp 372 z_elem_b(:,:,1) = 1._wp 373 ! at k=2 set de/dz=0.: 374 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 375 z_elem_a(:,:,2) = 0._wp 376 ! 377 ENDIF 378 ! 339 ! 340 ! Dirichlet conditions at k=1 341 en(:,:,1) = rc02r * ustars2(:,:) * (1._wp + rsbc_tke1)**(2._wp/3._wp) 342 en(:,:,1) = MAX(en(:,:,1), rn_emin) 343 z_elem_a(:,:,1) = en(:,:,1) 344 z_elem_c(:,:,1) = 0._wp 345 z_elem_b(:,:,1) = 1._wp 346 ! 347 ! at k=2, set de/dz=Fw 348 !cbr 349 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 350 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) 355 ! 356 ! 379 357 END SELECT 380 358 … … 382 360 ! -------------------------------- 383 361 ! 384 SELECT CASE ( nn_ tkebc_bot )362 SELECT CASE ( nn_bc_bot ) 385 363 ! 386 364 CASE ( 0 ) ! Dirichlet … … 457 435 ! ! set the minimum value of tke 458 436 en(:,:,:) = MAX( en(:,:,:), rn_emin ) 459 437 460 438 !!----------------------------------------!! 461 439 !! Solve prognostic equation for psi !! … … 560 538 ! --------------------------------- 561 539 ! 562 SELECT CASE ( nn_ psibc_surf )540 SELECT CASE ( nn_bc_surf ) 563 541 ! 564 542 CASE ( 0 ) ! Dirichlet boundary conditions 565 ! 566 IF( ln_crban ) THEN ! Wave induced mixing case 567 ! ! en(1) = q2(1) = 0.5 * (15.8 * Ccb)^(2/3) * u*^2 568 ! ! balance between the production and the dissipation terms including the wave effect 569 zdep(:,:) = rl_sf * zhsro(:,:) 570 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 571 z_elem_a(:,:,1) = psi(:,:,1) 572 z_elem_c(:,:,1) = 0._wp 573 z_elem_b(:,:,1) = 1._wp 574 ! 575 ! one level below 576 zex1 = (rmm*ra_sf+rnn) 577 zex2 = (rmm*ra_sf) 578 zdep(:,:) = ( (zhsro(:,:) + fsdepw(:,:,2))**zex1 ) / zhsro(:,:)**zex2 579 psi (:,:,2) = rsbc_psi1 * ustars2(:,:)**rmm * zdep(:,:) * tmask(:,:,1) 580 z_elem_a(:,:,2) = 0._wp 581 z_elem_c(:,:,2) = 0._wp 582 z_elem_b(:,:,2) = 1._wp 583 ! 584 ELSE ! No wave induced mixing case 585 ! ! en(1) = u*^2/C0^2 & l(1) = K*zs 586 ! ! balance between the production and the dissipation terms 587 ! 588 zdep(:,:) = vkarmn * zhsro(:,:) 589 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 590 z_elem_a(:,:,1) = psi(:,:,1) 591 z_elem_c(:,:,1) = 0._wp 592 z_elem_b(:,:,1) = 1._wp 593 ! 594 ! one level below 595 zdep(:,:) = vkarmn * ( zhsro(:,:) + fsdepw(:,:,2) ) 596 psi (:,:,2) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 597 z_elem_a(:,:,2) = 0._wp 598 z_elem_c(:,:,2) = 0._wp 599 z_elem_b(:,:,2) = 1. 600 ! 601 ENDIF 602 ! 543 ! 544 ! Surface value 545 zdep(:,:) = zhsro(:,:) * rl_sf ! Cosmetic 546 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 547 z_elem_a(:,:,1) = psi(:,:,1) 548 z_elem_c(:,:,1) = 0._wp 549 z_elem_b(:,:,1) = 1._wp 550 ! 551 ! One level below 552 zkar(:,:) = (rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdepw(:,:,2)/zhsro(:,:) ))) 553 zdep(:,:) = (zhsro(:,:) + fsdepw(:,:,2)) * zkar(:,:) 554 psi (:,:,2) = rc0**rpp * en(:,:,2)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 555 z_elem_a(:,:,2) = 0._wp 556 z_elem_c(:,:,2) = 0._wp 557 z_elem_b(:,:,2) = 1._wp 558 ! 559 ! 603 560 CASE ( 1 ) ! Neumann boundary condition on d(psi)/dz 604 ! 605 IF( ln_crban ) THEN ! Wave induced mixing case 606 ! 607 zdep(:,:) = rl_sf * zhsro(:,:) 608 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 609 z_elem_a(:,:,1) = psi(:,:,1) 610 z_elem_c(:,:,1) = 0._wp 611 z_elem_b(:,:,1) = 1._wp 612 ! 613 ! Neumann condition at k=2 614 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 615 z_elem_a(:,:,2) = 0._wp 616 ! 617 ! Set psi vertical flux at the surface: 618 zdep(:,:) = (zhsro(:,:) + fsdept(:,:,1))**(rmm*ra_sf+rnn-1._wp) / zhsro(:,:)**(rmm*ra_sf) 619 zflxs(:,:) = rsbc_psi3 * ( zwall_psi(:,:,1)*avm(:,:,1) + zwall_psi(:,:,2)*avm(:,:,2) ) & 620 & * en(:,:,1)**rmm * zdep 621 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 622 ! 623 ELSE ! No wave induced mixing 624 ! 625 zdep(:,:) = vkarmn * zhsro(:,:) 626 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 627 z_elem_a(:,:,1) = psi(:,:,1) 628 z_elem_c(:,:,1) = 0._wp 629 z_elem_b(:,:,1) = 1._wp 630 ! 631 ! Neumann condition at k=2 632 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 633 z_elem_a(ji,jj,2) = 0._wp 634 ! 635 ! Set psi vertical flux at the surface: 636 zdep(:,:) = zhsro(:,:) + fsdept(:,:,1) 637 zflxs(:,:) = rsbc_psi2 * ( avm(:,:,1) + avm(:,:,2) ) * en(:,:,1)**rmm * zdep**(rnn-1._wp) 638 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 639 ! 640 ENDIF 641 ! 561 ! 562 ! Surface value: Dirichlet 563 zdep(:,:) = zhsro(:,:) * rl_sf 564 psi (:,:,1) = rc0**rpp * en(:,:,1)**rmm * zdep(:,:)**rnn * tmask(:,:,1) 565 z_elem_a(:,:,1) = psi(:,:,1) 566 z_elem_c(:,:,1) = 0._wp 567 z_elem_b(:,:,1) = 1._wp 568 ! 569 ! Neumann condition at k=2 570 z_elem_b(:,:,2) = z_elem_b(:,:,2) + z_elem_a(:,:,2) ! Remove z_elem_a from z_elem_b 571 z_elem_a(:,:,2) = 0._wp 572 ! 573 ! Set psi vertical flux at the surface: 574 zkar(:,:) = rl_sf + (vkarmn-rl_sf)*(1._wp-exp(-rtrans*fsdept(:,:,1)/zhsro(:,:) )) ! Lengh scale slope 575 zdep(:,:) = ((zhsro(:,:) + fsdept(:,:,1)) / zhsro(:,:))**(rmm*ra_sf) 576 zflxs(:,:) = (rnn + rsbc_tke1 * (rnn + rmm*ra_sf) * zdep(:,:))*(1._wp + rsbc_tke1*zdep(:,:))**(2._wp*rmm/3._wp-1_wp) 577 zdep(:,:) = rsbc_psi1 * (zwall_psi(:,:,1)*avm(:,:,1)+zwall_psi(:,:,2)*avm(:,:,2)) * & 578 & ustars2(:,:)**rmm * zkar(:,:)**rnn * (zhsro(:,:) + fsdept(:,:,1))**(rnn-1.) 579 zflxs(:,:) = zdep(:,:) * zflxs(:,:) 580 psi(:,:,2) = psi(:,:,2) + zflxs(:,:) / fse3w(:,:,2) 581 582 ! 583 ! 642 584 END SELECT 643 585 … … 645 587 ! -------------------------------- 646 588 ! 647 SELECT CASE ( nn_psibc_bot ) 589 SELECT CASE ( nn_bc_bot ) 590 ! 648 591 ! 649 592 CASE ( 0 ) ! Dirichlet 650 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * hbro593 ! ! en(ibot) = u*^2 / Co2 and mxln(ibot) = vkarmn * rn_bfrz0 651 594 ! ! Balance between the production and the dissipation terms 652 595 !CDIR NOVERRCHK … … 656 599 ibot = mbkt(ji,jj) + 1 ! k bottom level of w-point 657 600 ibotm1 = mbkt(ji,jj) ! k-1 bottom level of w-point but >=1 658 zdep(ji,jj) = vkarmn * hbro601 zdep(ji,jj) = vkarmn * rn_bfrz0 659 602 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 660 603 z_elem_a(ji,jj,ibot) = 0._wp … … 663 606 ! 664 607 ! Just above last level, Dirichlet condition again (GOTM like) 665 zdep(ji,jj) = vkarmn * ( hbro+ fse3t(ji,jj,ibotm1) )608 zdep(ji,jj) = vkarmn * ( rn_bfrz0 + fse3t(ji,jj,ibotm1) ) 666 609 psi (ji,jj,ibotm1) = rc0**rpp * en(ji,jj,ibot )**rmm * zdep(ji,jj)**rnn 667 610 z_elem_a(ji,jj,ibotm1) = 0._wp … … 681 624 ! 682 625 ! Bottom level Dirichlet condition: 683 zdep(ji,jj) = vkarmn * hbro626 zdep(ji,jj) = vkarmn * rn_bfrz0 684 627 psi (ji,jj,ibot) = rc0**rpp * en(ji,jj,ibot)**rmm * zdep(ji,jj)**rnn 685 628 ! … … 693 636 ! 694 637 ! Set psi vertical flux at the bottom: 695 zdep(ji,jj) = hbro+ 0.5_wp*fse3t(ji,jj,ibotm1)638 zdep(ji,jj) = rn_bfrz0 + 0.5_wp*fse3t(ji,jj,ibotm1) 696 639 zflxb = rsbc_psi2 * ( avm(ji,jj,ibot) + avm(ji,jj,ibotm1) ) & 697 640 & * (0.5_wp*(en(ji,jj,ibot)+en(ji,jj,ibotm1)))**rmm * zdep(ji,jj)**(rnn-1._wp) … … 736 679 DO jj = 2, jpjm1 737 680 DO ji = fs_2, fs_jpim1 ! vector opt. 738 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / psi(ji,jj,jk)681 eps(ji,jj,jk) = rc03 * en(ji,jj,jk) * en(ji,jj,jk) * SQRT( en(ji,jj,jk) ) / MAX( psi(ji,jj,jk), rn_epsmin) 739 682 END DO 740 683 END DO … … 783 726 ! Galperin criterium (NOTE : Not required if the proper value of C3 in stable cases is calculated) 784 727 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 785 mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk))728 IF (ln_length_lim) mxln(ji,jj,jk) = MIN( rn_clim_galp * SQRT( 2._wp * en(ji,jj,jk) / zrn2 ), mxln(ji,jj,jk) ) 786 729 END DO 787 730 END DO … … 847 790 ! Boundary conditions on stability functions for momentum (Neumann): 848 791 ! Lines below are useless if GOTM style Dirichlet conditions are used 849 zcoef = rcm_sf / SQRT( 2._wp ) 792 793 avmv(:,:,1) = avmv(:,:,2) 794 850 795 DO jj = 2, jpjm1 851 796 DO ji = fs_2, fs_jpim1 ! vector opt. 852 avmv(ji,jj,1) = zcoef 853 END DO 854 END DO 855 zcoef = rc0 / SQRT( 2._wp ) 856 DO jj = 2, jpjm1 857 DO ji = fs_2, fs_jpim1 ! vector opt. 858 avmv(ji,jj,mbkt(ji,jj)+1) = zcoef 797 avmv(ji,jj,mbkt(ji,jj)+1) = avmv(ji,jj,mbkt(ji,jj)) 859 798 END DO 860 799 END DO … … 900 839 avmv_k(:,:,:) = avmv(:,:,:) 901 840 ! 902 CALL wrk_dealloc( jpi,jpj, zdep, z flxs, zhsro )841 CALL wrk_dealloc( jpi,jpj, zdep, zkar, zflxs, zhsro ) 903 842 CALL wrk_dealloc( jpi,jpj,jpk, eb, mxlb, shear, eps, zwall_psi, z_elem_a, z_elem_b, z_elem_c, psi ) 904 843 ! … … 932 871 !! 933 872 NAMELIST/namzdf_gls/rn_emin, rn_epsmin, ln_length_lim, & 934 & rn_clim_galp, ln_crban, ln_sigpsi, & 935 & rn_crban, rn_charn, & 936 & nn_tkebc_surf, nn_tkebc_bot, & 937 & nn_psibc_surf, nn_psibc_bot, & 873 & rn_clim_galp, ln_sigpsi, rn_hsro, & 874 & rn_crban, rn_charn, rn_frac_hs, & 875 & nn_bc_surf, nn_bc_bot, nn_z0_met, & 938 876 & nn_stab_func, nn_clos 939 877 !!---------------------------------------------------------- … … 955 893 WRITE(numout,*) '~~~~~~~~~~~~' 956 894 WRITE(numout,*) ' Namelist namzdf_gls : set gls mixing parameters' 957 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 958 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 959 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 960 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 961 WRITE(numout,*) ' TKE Surface boundary condition nn_tkebc_surf = ', nn_tkebc_surf 962 WRITE(numout,*) ' TKE Bottom boundary condition nn_tkebc_bot = ', nn_tkebc_bot 963 WRITE(numout,*) ' PSI Surface boundary condition nn_psibc_surf = ', nn_psibc_surf 964 WRITE(numout,*) ' PSI Bottom boundary condition nn_psibc_bot = ', nn_psibc_bot 965 WRITE(numout,*) ' Craig and Banner scheme ln_crban = ', ln_crban 966 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 895 WRITE(numout,*) ' minimum value of en rn_emin = ', rn_emin 896 WRITE(numout,*) ' minimum value of eps rn_epsmin = ', rn_epsmin 897 WRITE(numout,*) ' Limit dissipation rate under stable stratif. ln_length_lim = ', ln_length_lim 898 WRITE(numout,*) ' Galperin limit (Standard: 0.53, Holt: 0.26) rn_clim_galp = ', rn_clim_galp 899 WRITE(numout,*) ' TKE Surface boundary condition nn_bc_surf = ', nn_bc_surf 900 WRITE(numout,*) ' TKE Bottom boundary condition nn_bc_bot = ', nn_bc_bot 901 WRITE(numout,*) ' Modify psi Schmidt number (wb case) ln_sigpsi = ', ln_sigpsi 967 902 WRITE(numout,*) ' Craig and Banner coefficient rn_crban = ', rn_crban 968 903 WRITE(numout,*) ' Charnock coefficient rn_charn = ', rn_charn 904 WRITE(numout,*) ' Surface roughness formula nn_z0_met = ', nn_z0_met 905 WRITE(numout,*) ' Wave height frac. (used if nn_z0_met=2) rn_frac_hs = ', rn_frac_hs 969 906 WRITE(numout,*) ' Stability functions nn_stab_func = ', nn_stab_func 970 907 WRITE(numout,*) ' Type of closure nn_clos = ', nn_clos 971 WRITE(numout,*) ' Hard coded parameters' 972 WRITE(numout,*) ' Surface roughness (m) hsro = ', hsro 973 WRITE(numout,*) ' Bottom roughness (m) hbro = ', hbro 908 WRITE(numout,*) ' Surface roughness (m) rn_hsro = ', rn_hsro 909 WRITE(numout,*) ' Bottom roughness (m) (nambfr namelist) rn_bfrz0 = ', rn_bfrz0 974 910 ENDIF 975 911 … … 978 914 979 915 ! !* Check of some namelist values 980 IF( nn_tkebc_surf < 0 .OR. nn_tkebc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_surf is 0 or 1' ) 981 IF( nn_psibc_surf < 0 .OR. nn_psibc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_surf is 0 or 1' ) 982 IF( nn_tkebc_bot < 0 .OR. nn_tkebc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_tkebc_bot is 0 or 1' ) 983 IF( nn_psibc_bot < 0 .OR. nn_psibc_bot > 1 ) CALL ctl_stop( 'bad flag: nn_psibc_bot is 0 or 1' ) 916 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 917 IF( nn_bc_surf < 0 .OR. nn_bc_surf > 1 ) CALL ctl_stop( 'bad flag: nn_bc_surf is 0 or 1' ) 918 IF( nn_z0_met < 0 .OR. nn_z0_met > 2 ) CALL ctl_stop( 'bad flag: nn_z0_met is 0, 1 or 2' ) 984 919 IF( nn_stab_func < 0 .OR. nn_stab_func > 3 ) CALL ctl_stop( 'bad flag: nn_stab_func is 0, 1, 2 and 3' ) 985 920 IF( nn_clos < 0 .OR. nn_clos > 3 ) CALL ctl_stop( 'bad flag: nn_clos is 0, 1, 2 or 3' ) … … 1001 936 SELECT CASE ( nn_stab_func ) 1002 937 CASE( 0, 1 ) ; rpsi3m = 2.53_wp ! G88 or KC stability functions 1003 CASE( 2 ) ; rpsi3m = 2. 38_wp ! Canuto A stability functions938 CASE( 2 ) ; rpsi3m = 2.62_wp ! Canuto A stability functions 1004 939 CASE( 3 ) ; rpsi3m = 2.38 ! Canuto B stability functions (caution : constant not identified) 1005 940 END SELECT … … 1012 947 rnn = -1._wp 1013 948 rsc_tke = 1._wp 1014 rsc_psi = 1. 3_wp ! Schmidt number for psi949 rsc_psi = 1.2_wp ! Schmidt number for psi 1015 950 rpsi1 = 1.44_wp 1016 951 rpsi3p = 1._wp … … 1140 1075 ! ! See Eq. (13) of Carniel et al, OM, 30, 225-239, 2009 1141 1076 ! ! or Eq. (17) of Burchard, JPO, 31, 3133-3145, 2001 1142 IF( ln_sigpsi .AND. ln_crban ) THEN 1143 zcr = SQRT( 1.5_wp*rsc_tke ) * rcm_sf / vkarmn 1144 rsc_psi0 = vkarmn*vkarmn / ( rpsi2 * rcm_sf*rcm_sf ) & 1145 & * ( rnn*rnn - 4._wp/3._wp * zcr*rnn*rmm - 1._wp/3._wp * zcr*rnn & 1146 & + 2._wp/9._wp * rmm * zcr*zcr + 4._wp/9._wp * zcr*zcr * rmm*rmm ) 1077 IF( ln_sigpsi ) THEN 1078 ra_sf = -1.5 ! Set kinetic energy slope, then deduce rsc_psi and rl_sf 1079 ! Verification: retrieve Burchard (2001) results by uncomenting the line below: 1080 ! Note that the results depend on the value of rn_cm_sf which is constant (=rc0) in his work 1081 ! ra_sf = -SQRT(2./3.*rc0**3./rn_cm_sf*rn_sc_tke)/vkarmn 1082 rsc_psi0 = rsc_tke/(24.*rpsi2)*(-1.+(4.*rnn + ra_sf*(1.+4.*rmm))**2./(ra_sf**2.)) 1147 1083 ELSE 1148 1084 rsc_psi0 = rsc_psi … … 1151 1087 ! !* Shear free turbulence parameters 1152 1088 ! 1153 ra_sf = -4._wp * rnn * SQRT( rsc_tke ) / ( (1._wp+4._wp*rmm) * SQRT( rsc_tke ) & 1154 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1155 rl_sf = rc0 * SQRT( rc0 / rcm_sf ) & 1156 & * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm*rmm) * rsc_tke & 1157 & + 12._wp * rsc_psi0 * rpsi2 & 1158 & - (1._wp + 4._wp*rmm) * SQRT( rsc_tke*(rsc_tke+ 24._wp*rsc_psi0*rpsi2) ) ) & 1159 & / ( 12._wp*rnn*rnn ) ) 1089 ra_sf = -4._wp*rnn*SQRT(rsc_tke) / ( (1._wp+4._wp*rmm)*SQRT(rsc_tke) & 1090 & - SQRT(rsc_tke + 24._wp*rsc_psi0*rpsi2 ) ) 1091 1092 IF ( rn_crban==0._wp ) THEN 1093 rl_sf = vkarmn 1094 ELSE 1095 rl_sf = rc0 * SQRT(rc0/rcm_sf) * SQRT( ( (1._wp + 4._wp*rmm + 8._wp*rmm**2_wp)*rsc_tke & 1096 & + 12._wp * rsc_psi0*rpsi2 - (1._wp + 4._wp*rmm) & 1097 & *SQRT(rsc_tke*(rsc_tke & 1098 & + 24._wp*rsc_psi0*rpsi2)) ) & 1099 & /(12._wp*rnn**2.) & 1100 & ) 1101 ENDIF 1160 1102 1161 1103 ! … … 1187 1129 rc03 = rc02 * rc0 1188 1130 rc04 = rc03 * rc0 1189 rc03_sqrt2_galp = rc03 / SQRT(2._wp) / rn_clim_galp 1190 rsbc_mb = 0.5_wp * (15.8_wp*rn_crban)**(2._wp/3._wp) ! Surf. bound. cond. from Mellor and Blumberg 1191 rsbc_std = 3.75_wp ! Surf. bound. cond. standard (prod=diss) 1192 rsbc_tke1 = (-rsc_tke*rn_crban/(rcm_sf*ra_sf*rl_sf))**(2._wp/3._wp) ! k_eps = 53. Dirichlet + Wave breaking 1193 rsbc_tke2 = 0.5_wp / rau0 1194 rsbc_tke3 = rdt * rn_crban ! Neumann + Wave breaking 1195 rsbc_zs = rn_charn / grav ! Charnock formula 1196 rsbc_psi1 = rc0**rpp * rsbc_tke1**rmm * rl_sf**rnn ! Dirichlet + Wave breaking 1197 rsbc_psi2 = -0.5_wp * rdt * rc0**rpp * rnn * vkarmn**rnn / rsc_psi ! Neumann + NO Wave breaking 1198 rsbc_psi3 = -0.5_wp * rdt * rc0**rpp * rl_sf**rnn / rsc_psi * (rnn + rmm*ra_sf) ! Neumann + Wave breaking 1199 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1200 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1131 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 1133 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 1135 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 1137 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 1139 1140 rfact_tke = -0.5_wp / rsc_tke * rdt ! Cst used for the Diffusion term of tke 1141 rfact_psi = -0.5_wp / rsc_psi * rdt ! Cst used for the Diffusion term of tke 1201 1142 1202 1143 ! !* Wall proximity function … … 1257 1198 IF(lwp) WRITE(numout,*) ' ===>>>> : previous run without gls scheme, en and mxln computed by iterative loop' 1258 1199 en (:,:,:) = rn_emin 1259 mxln(:,:,:) = 0.0 011200 mxln(:,:,:) = 0.05 1260 1201 avt_k (:,:,:) = avt (:,:,:) 1261 1202 avm_k (:,:,:) = avm (:,:,:) … … 1267 1208 IF(lwp) WRITE(numout,*) ' ===>>>> : Initialisation of en and mxln by background values' 1268 1209 en (:,:,:) = rn_emin 1269 mxln(:,:,:) = 0.0 011210 mxln(:,:,:) = 0.05 1270 1211 ENDIF 1271 1212 ! … … 1273 1214 ! ! ------------------- 1274 1215 IF(lwp) WRITE(numout,*) '---- gls-rst ----' 1275 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1216 CALL iom_rstput( kt, nitrst, numrow, 'en' , en ) 1276 1217 CALL iom_rstput( kt, nitrst, numrow, 'avt' , avt_k ) 1277 1218 CALL iom_rstput( kt, nitrst, numrow, 'avm' , avm_k ) 1278 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1219 CALL iom_rstput( kt, nitrst, numrow, 'avmu' , avmu_k ) 1279 1220 CALL iom_rstput( kt, nitrst, numrow, 'avmv' , avmv_k ) 1280 1221 CALL iom_rstput( kt, nitrst, numrow, 'mxln' , mxln ) -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfini.F90
r5477 r5572 14 14 !!---------------------------------------------------------------------- 15 15 USE par_oce ! mesh and scale factors 16 USE sbc_oce ! surface module (only for nn_isf in the option compatibility test)17 16 USE ldftra_oce ! ocean active tracers: lateral physics 18 17 USE ldfdyn_oce ! ocean dynamics lateral physics … … 118 117 IF( ioptio == 0 .OR. ioptio > 1 .AND. .NOT. lk_esopa ) & 119 118 & CALL ctl_stop( ' one and only one vertical diffusion option has to be defined ' ) 120 IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. nn_isf .NE. 0) &119 IF( ( lk_zdfric .OR. lk_zdfgls .OR. lk_zdfkpp ) .AND. ln_isfcav ) & 121 120 & CALL ctl_stop( ' only zdfcst and zdftke were tested with ice shelves cavities ' ) 122 121 ! … … 125 124 IF(lwp) WRITE(numout,*) ' convection :' 126 125 ! 127 IF( ln_zdfnpc ) CALL ctl_stop( ' zdf_init: non penetrative convective scheme is not working', & 128 & ' set ln_zdfnpc to FALSE' ) 126 #if defined key_top 127 IF( ln_zdfnpc ) CALL ctl_stop( ' zdf_init: npc scheme is not working with key_top' ) 128 #endif 129 129 ! 130 130 ioptio = 0 -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r5477 r5572 26 26 !! ! + cleaning of the parameters + bugs correction 27 27 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 28 !! 3.6 ! 2014-11 (P. Mathiot) add ice shelf capability 28 29 !!---------------------------------------------------------------------- 29 30 #if defined key_zdftke || defined key_esopa … … 236 237 zfact3 = 0.5_wp * rn_ediss 237 238 ! 239 ! 238 240 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 239 241 ! ! Surface boundary condition on tke 240 242 ! !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 243 IF ( ln_isfcav ) THEN 244 DO jj = 2, jpjm1 ! en(mikt(ji,jj)) = rn_emin 245 DO ji = fs_2, fs_jpim1 ! vector opt. 246 en(ji,jj,mikt(ji,jj))=rn_emin * tmask(ji,jj,1) 247 END DO 248 END DO 249 END IF 241 250 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 251 DO ji = fs_2, fs_jpim1 ! vector opt. 243 IF (mikt(ji,jj) .GT. 1) THEN 244 en(ji,jj,mikt(ji,jj))=rn_emin 245 ELSE 246 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 247 END IF 252 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 248 253 END DO 249 254 END DO … … 301 306 END DO 302 307 zcof = 0.016 / SQRT( zrhoa * zcdrag ) 308 !CDIR NOVERRCHK 303 309 DO jk = 2, jpkm1 !* TKE Langmuir circulation source term added to en 304 DO jj = 2, jpjm1 310 !CDIR NOVERRCHK 311 DO jj = 2, jpjm1 312 !CDIR NOVERRCHK 305 313 DO ji = fs_2, fs_jpim1 ! vector opt. 306 314 zus = zcof * SQRT( taum(ji,jj) ) ! Stokes drift … … 309 317 zwlc = zind * rn_lc * zus * SIN( rpi * fsdepw(ji,jj,jk) / zhlc(ji,jj) ) 310 318 ! ! TKE Langmuir circulation source term 311 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * tmask(ji,jj,jk)319 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 312 320 END DO 313 321 END DO … … 328 336 avmu(ji,jj,jk) = avmu(ji,jj,jk) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) & 329 337 & * ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) & 330 & / ( fse3uw_n(ji,jj,jk)&331 & * fse3uw_b(ji,jj,jk))338 & / ( fse3uw_n(ji,jj,jk) & 339 & * fse3uw_b(ji,jj,jk) ) 332 340 avmv(ji,jj,jk) = avmv(ji,jj,jk) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) & 333 341 & * ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) & … … 338 346 END DO 339 347 ! 340 DO j j = 2, jpjm1341 DO j i = fs_2, fs_jpim1 ! vector opt.342 DO j k = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en348 DO jk = 2, jpkm1 !* Matrix and right hand side in en 349 DO jj = 2, jpjm1 350 DO ji = fs_2, fs_jpim1 ! vector opt. 343 351 zcof = zfact1 * tmask(ji,jj,jk) 344 352 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 357 365 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 358 366 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 359 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 360 END DO 361 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 362 DO jk = mikt(ji,jj)+2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 367 & * wmask(ji,jj,jk) 368 END DO 369 END DO 370 END DO 371 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 372 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 373 DO jj = 2, jpjm1 374 DO ji = fs_2, fs_jpim1 ! vector opt. 363 375 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 364 376 END DO 365 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 366 zd_lw(ji,jj,mikt(ji,jj)+1) = en(ji,jj,mikt(ji,jj)+1) - zd_lw(ji,jj,mikt(ji,jj)+1) * en(ji,jj,mikt(ji,jj)) ! Surface boudary conditions on tke 367 ! 368 DO jk = mikt(ji,jj)+2, jpkm1 377 END DO 378 END DO 379 ! 380 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 384 END DO 385 END DO 386 DO jk = 3, jpkm1 387 DO jj = 2, jpjm1 388 DO ji = fs_2, fs_jpim1 ! vector opt. 369 389 zd_lw(ji,jj,jk) = en(ji,jj,jk) - zd_lw(ji,jj,jk) / zdiag(ji,jj,jk-1) *zd_lw(ji,jj,jk-1) 370 390 END DO 371 ! 372 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 391 END DO 392 END DO 393 ! 394 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 395 DO jj = 2, jpjm1 396 DO ji = fs_2, fs_jpim1 ! vector opt. 373 397 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 374 ! 375 DO jk = jpk-2, mikt(ji,jj)+1, -1 398 END DO 399 END DO 400 DO jk = jpk-2, 2, -1 401 DO jj = 2, jpjm1 402 DO ji = fs_2, fs_jpim1 ! vector opt. 376 403 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 377 404 END DO 378 ! 379 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 380 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 405 END DO 406 END DO 407 DO jk = 2, jpkm1 ! set the minimum value of tke 408 DO jj = 2, jpjm1 409 DO ji = fs_2, fs_jpim1 ! vector opt. 410 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * wmask(ji,jj,jk) 381 411 END DO 382 412 END DO … … 391 421 DO ji = fs_2, fs_jpim1 ! vector opt. 392 422 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 393 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)423 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 394 424 END DO 395 425 END DO … … 400 430 jk = nmln(ji,jj) 401 431 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 402 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1)432 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 403 433 END DO 404 434 END DO … … 416 446 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 417 447 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 418 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1)448 & * ( 1._wp - fr_i(ji,jj) ) * wmask(ji,jj,jk) * tmask(ji,jj,1) 419 449 END DO 420 450 END DO … … 484 514 ! !* Buoyancy length scale: l=sqrt(2*e/n**2) 485 515 ! 516 ! initialisation of interior minimum value (avoid a 2d loop with mikt) 517 zmxlm(:,:,:) = rmxl_min 518 zmxld(:,:,:) = rmxl_min 519 ! 486 520 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 487 521 DO jj = 2, jpjm1 488 522 DO ji = fs_2, fs_jpim1 489 IF (mikt(ji,jj) .GT. 1) THEN 490 zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 491 ELSE 492 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 493 zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 494 END IF 523 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 524 zmxlm(ji,jj,1) = MAX( rn_mxl0, zraug * taum(ji,jj) * tmask(ji,jj,1) ) 495 525 END DO 496 526 END DO 497 527 ELSE 498 DO jj = 2, jpjm1 499 DO ji = fs_2, fs_jpim1 ! surface set to the minimum value 500 zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 501 END DO 502 END DO 528 zmxlm(:,:,1) = rn_mxl0 503 529 ENDIF 504 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 505 ! 506 !CDIR NOVERRCHK 507 DO jj = 2, jpjm1 508 !CDIR NOVERRCHK 509 DO ji = fs_2, fs_jpim1 ! vector opt. 510 !CDIR NOVERRCHK 511 DO jk = mikt(ji,jj)+1, jpkm1 ! interior value : l=sqrt(2*e/n^2) 530 ! 531 !CDIR NOVERRCHK 532 DO jk = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2) 533 !CDIR NOVERRCHK 534 DO jj = 2, jpjm1 535 !CDIR NOVERRCHK 536 DO ji = fs_2, fs_jpim1 ! vector opt. 512 537 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 513 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 514 END DO 515 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 538 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 539 END DO 516 540 END DO 517 541 END DO … … 519 543 ! !* Physical limits for the mixing length 520 544 ! 521 zmxld(:,:, 1 ) = zmxlm(:,:,1) ! surface set to the zmxlm value545 zmxld(:,:,1 ) = zmxlm(:,:,1) ! surface set to the minimum value 522 546 zmxld(:,:,jpk) = rmxl_min ! last level set to the minimum value 523 547 ! 524 548 SELECT CASE ( nn_mxl ) 525 549 ! 550 ! where wmask = 0 set zmxlm == fse3w 526 551 CASE ( 0 ) ! bounded by the distance to surface and bottom 527 DO j j = 2, jpjm1528 DO j i = fs_2, fs_jpim1 ! vector opt.529 DO j k = mikt(ji,jj)+1, jpkm1552 DO jk = 2, jpkm1 553 DO jj = 2, jpjm1 554 DO ji = fs_2, fs_jpim1 ! vector opt. 530 555 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 531 556 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 532 zmxlm(ji,jj,jk) = zemxl 533 zmxld(ji,jj,jk) = zemxl 557 ! wmask prevent zmxlm = 0 if jk = mikt(ji,jj) 558 zmxlm(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 559 zmxld(ji,jj,jk) = zemxl * wmask(ji,jj,jk) + MIN(zmxlm(ji,jj,jk),fse3w(ji,jj,jk)) * (1 - wmask(ji,jj,jk)) 534 560 END DO 535 561 END DO … … 537 563 ! 538 564 CASE ( 1 ) ! bounded by the vertical scale factor 539 DO j j = 2, jpjm1540 DO j i = fs_2, fs_jpim1 ! vector opt.541 DO j k = mikt(ji,jj)+1, jpkm1565 DO jk = 2, jpkm1 566 DO jj = 2, jpjm1 567 DO ji = fs_2, fs_jpim1 ! vector opt. 542 568 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 543 569 zmxlm(ji,jj,jk) = zemxl … … 548 574 ! 549 575 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 550 DO j j = 2, jpjm1551 DO j i = fs_2, fs_jpim1 ! vector opt.552 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom :576 DO jk = 2, jpkm1 ! from the surface to the bottom : 577 DO jj = 2, jpjm1 578 DO ji = fs_2, fs_jpim1 ! vector opt. 553 579 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 554 580 END DO 555 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 581 END DO 582 END DO 583 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 556 586 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 557 587 zmxlm(ji,jj,jk) = zemxl … … 562 592 ! 563 593 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 564 DO j j = 2, jpjm1565 DO j i = fs_2, fs_jpim1 ! vector opt.566 DO j k = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup594 DO jk = 2, jpkm1 ! from the surface to the bottom : lup 595 DO jj = 2, jpjm1 596 DO ji = fs_2, fs_jpim1 ! vector opt. 567 597 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 568 598 END DO 569 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 599 END DO 600 END DO 601 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 602 DO jj = 2, jpjm1 603 DO ji = fs_2, fs_jpim1 ! vector opt. 570 604 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 571 605 END DO … … 604 638 zsqen = SQRT( en(ji,jj,jk) ) 605 639 zav = rn_ediff * zmxlm(ji,jj,jk) * zsqen 606 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * tmask(ji,jj,jk)607 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)640 avm (ji,jj,jk) = MAX( zav, avmb(jk) ) * wmask(ji,jj,jk) 641 avt (ji,jj,jk) = MAX( zav, avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 608 642 dissl(ji,jj,jk) = zsqen / zmxld(ji,jj,jk) 609 643 END DO … … 612 646 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 613 647 ! 614 DO jj = 2, jpjm1 615 DO ji = fs_2, fs_jpim1 ! vector opt. 616 DO jk = miku(ji,jj)+1, jpkm1 !* vertical eddy viscosity at u- and v-points 617 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 618 END DO 619 DO jk = mikv(ji,jj)+1, jpkm1 620 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 648 DO jk = 2, jpkm1 !* vertical eddy viscosity at wu- and wv-points 649 DO jj = 2, jpjm1 650 DO ji = fs_2, fs_jpim1 ! vector opt. 651 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 652 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 621 653 END DO 622 654 END DO … … 625 657 ! 626 658 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 627 DO j j = 2, jpjm1628 DO j i = fs_2, fs_jpim1 ! vector opt.629 DO j k = mikt(ji,jj)+1, jpkm1659 DO jk = 2, jpkm1 660 DO jj = 2, jpjm1 661 DO ji = fs_2, fs_jpim1 ! vector opt. 630 662 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 631 663 ! ! shear … … 639 671 !!gm and even better with the use of the "true" ri_crit=0.22222... (this change the results!) 640 672 !!gm zpdlr = MAX( 0.1_wp, ri_crit / MAX( ri_crit , zri ) ) 641 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk)673 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * wmask(ji,jj,jk) 642 674 # if defined key_c1d 643 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number644 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri675 e_pdl(ji,jj,jk) = zpdlr * wmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 676 e_ric(ji,jj,jk) = zri * wmask(ji,jj,jk) ! c1d config. : save Ri 645 677 # endif 646 678 END DO … … 729 761 IF( nn_pdl < 0 .OR. nn_pdl > 1 ) CALL ctl_stop( 'bad flag: nn_pdl is 0 or 1 ' ) 730 762 IF( nn_htau < 0 .OR. nn_htau > 1 ) CALL ctl_stop( 'bad flag: nn_htau is 0, 1 or 2 ' ) 731 IF( nn_etau == 3 .AND. .NOT. l k_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' )763 IF( nn_etau == 3 .AND. .NOT. ln_cpl ) CALL ctl_stop( 'nn_etau == 3 : HF taum only known in coupled mode' ) 732 764 733 765 IF( ln_mxl0 ) THEN … … 749 781 ! !* set vertical eddy coef. to the background value 750 782 DO jk = 1, jpk 751 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)752 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)753 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)754 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)783 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 784 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 785 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 786 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 755 787 END DO 756 788 dissl(:,:,:) = 1.e-12_wp … … 803 835 en (:,:,:) = rn_emin * tmask(:,:,:) 804 836 CALL tke_avn ! recompute avt, avm, avmu, avmv and dissl (approximation) 837 ! 838 avt_k (:,:,:) = avt (:,:,:) 839 avm_k (:,:,:) = avm (:,:,:) 840 avmu_k(:,:,:) = avmu(:,:,:) 841 avmv_k(:,:,:) = avmv(:,:,:) 842 ! 805 843 DO jit = nit000 + 1, nit000 + 10 ; CALL zdf_tke( jit ) ; END DO 806 844 ENDIF … … 808 846 en(:,:,:) = rn_emin * tmask(:,:,:) 809 847 DO jk = 1, jpk ! set the Kz to the background value 810 avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)811 avm (:,:,jk) = avmb(jk) * tmask(:,:,jk)812 avmu(:,:,jk) = avmb(jk) * umask(:,:,jk)813 avmv(:,:,jk) = avmb(jk) * vmask(:,:,jk)848 avt (:,:,jk) = avtb(jk) * wmask (:,:,jk) 849 avm (:,:,jk) = avmb(jk) * wmask (:,:,jk) 850 avmu(:,:,jk) = avmb(jk) * wumask(:,:,jk) 851 avmv(:,:,jk) = avmb(jk) * wvmask(:,:,jk) 814 852 END DO 815 853 ENDIF -
branches/UKMO/dev_r5107_hadgem3_cplseq/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r5477 r5572 126 126 zkz(:,:) = 0.e0 !* Associated potential energy consummed over the whole water column 127 127 DO jk = 2, jpkm1 128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1)128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk) * wmask(:,:,jk) 129 129 END DO 130 130 … … 135 135 END DO 136 136 137 DO j j = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx138 DO j i = 1, jpi139 DO j k = mikt(ji,jj)+1, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s137 DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 138 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 139 DO ji = 1, jpi 140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s 141 141 END DO 142 142 END DO … … 166 166 ! ! Update mixing coefs ! 167 167 ! ! ----------------------- ! 168 DO j j = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx169 DO j i = 1, jpi170 DO j k = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing171 avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) 172 avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) 168 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 169 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 170 DO ji = 1, jpi 171 avt(ji,jj,jk) = avt(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 172 avm(ji,jj,jk) = avm(ji,jj,jk) + zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 173 173 END DO 174 174 END DO 175 175 END DO 176 176 177 DO j j = 2, jpjm1178 DO j i = fs_2, fs_jpim1 ! vector opt.179 DO j k = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing180 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * umask(ji,jj,jk)181 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * vmask(ji,jj,jk)177 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 178 DO jj = 2, jpjm1 179 DO ji = fs_2, fs_jpim1 ! vector opt. 180 avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji+1,jj ,jk) ) * wumask(ji,jj,jk) 181 avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5 * ( zav_tide(ji,jj,jk) + zav_tide(ji ,jj+1,jk) ) * wvmask(ji,jj,jk) 182 182 END DO 183 183 END DO … … 457 457 ztpc = 0.e0 458 458 zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 459 DO j j = 1, jpj460 DO j i = 1, jpi461 DO j k= mikt(ji,jj)+1, jpkm1462 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)459 DO jk= 2, jpkm1 460 DO jj = 1, jpj 461 DO ji = 1, jpi 462 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 463 463 END DO 464 464 END DO … … 473 473 zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 ) 474 474 zkz(:,:) = 0.e0 475 DO j j = 1, jpj476 DO j i = 1, jpi477 DO j k = mikt(ji,jj)+1, jpkm1478 zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX( 0.e0, rn2(ji,jj,jk) ) * rau0 * zav_tide(ji,jj,jk)* tmask(ji,jj,jk)475 DO jk = 2, jpkm1 476 DO jj = 1, jpj 477 DO ji = 1, jpi 478 zkz(ji,jj) = zkz(ji,jj) + fse3w(ji,jj,jk) * MAX(0.e0, rn2(ji,jj,jk)) * rau0 * zav_tide(ji,jj,jk) * wmask(ji,jj,jk) 479 479 END DO 480 480 END DO … … 498 498 WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 499 499 500 DO j j = 1, jpj501 DO j i = 1, jpi502 DO j k = mikt(ji,jj)+1, jpkm1503 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s500 DO jk = 2, jpkm1 501 DO jj = 1, jpj 502 DO ji = 1, jpi 503 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) * wmask(ji,jj,jk) !kz max = 300 cm2/s 504 504 END DO 505 505 END DO … … 510 510 DO jj = 1, jpj 511 511 DO ji = 1, jpi 512 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj)512 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 513 513 END DO 514 514 END DO … … 519 519 DO jk = 1, jpk 520 520 ze_z = SUM( e1t(:,:) * e2t(:,:) * zav_tide(:,:,jk) * tmask_i(:,:) ) & 521 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )521 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 522 522 ztpc = 1.E50 523 523 DO jj = 1, jpj … … 540 540 END DO 541 541 ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) & 542 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )542 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 543 543 WRITE(numout,*) ' jk= ', jk,' ', ze_z * 1.e4,' cm2/s' 544 544 END DO … … 546 546 zkz(:,:) = az_tmx(:,:,jk) /rn_n2min 547 547 ze_z = SUM( e1t(:,:) * e2t(:,:) * zkz(:,:) * tmask_i(:,:) ) & 548 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * tmask (:,:,jk) * tmask_i(:,:) ) )548 & / MAX( 1.e-20, SUM( e1t(:,:) * e2t(:,:) * wmask (:,:,jk) * tmask_i(:,:) ) ) 549 549 WRITE(numout,*) 550 550 WRITE(numout,*) ' N2 min - jk= ', jk,' ', ze_z * 1.e4,' cm2/s min= ',MINVAL(zkz)*1.e4, &
Note: See TracChangeset
for help on using the changeset viewer.
