- Timestamp:
- 2014-06-11T14:52:23+02:00 (10 years ago)
- Location:
- branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF
- Files:
-
- 5 edited
-
zdf_oce.F90 (modified) (2 diffs)
-
zdfbfr.F90 (modified) (18 diffs)
-
zdfddm.F90 (modified) (4 diffs)
-
zdftke.F90 (modified) (13 diffs)
-
zdftmx.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
r4147 r4666 40 40 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: avtb_2d !: horizontal shape of background Kz profile 41 41 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: bfrua, bfrva !: Bottom friction coefficients set in zdfbfr 42 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:) :: tfrua, tfrva !: top friction coefficients set in zdfbfr 42 43 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avmu , avmv !: vertical viscosity coef at uw- & vw-pts [m2/s] 43 44 REAL(wp), PUBLIC, SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avm , avt !: vertical viscosity & diffusivity coef at w-pt [m2/s] … … 57 58 ALLOCATE(avmb(jpk) , bfrua(jpi,jpj) , & 58 59 & avtb(jpk) , bfrva(jpi,jpj) , avtb_2d(jpi,jpj) , & 60 & tfrua(jpi, jpj), tfrva(jpi, jpj) , & 59 61 & avmu(jpi,jpj,jpk), avm(jpi,jpj,jpk) , & 60 62 & avmv(jpi,jpj,jpk), avt(jpi,jpj,jpk) , STAT = zdf_oce_alloc ) -
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r4624 r4666 41 41 REAL(wp), PUBLIC :: rn_bfrien ! local factor to enhance coefficient bfri (PUBLIC for TAM) 42 42 LOGICAL , PUBLIC :: ln_bfr2d ! logical switch for 2D enhancement (PUBLIC for TAM) 43 REAL(wp), PUBLIC :: rn_tfri1 ! top drag coefficient (linear case) (PUBLIC for TAM) 44 REAL(wp), PUBLIC :: rn_tfri2 ! top drag coefficient (non linear case) (PUBLIC for TAM) 45 REAL(wp), PUBLIC :: rn_tfri2_max ! Maximum top drag coefficient (non linear case and ln_loglayer=T) (PUBLIC for TAM) 46 REAL(wp), PUBLIC :: rn_tfeb2 ! background top turbulent kinetic energy [m2/s2] (PUBLIC for TAM) 47 REAL(wp), PUBLIC :: rn_tfrien ! local factor to enhance coefficient tfri (PUBLIC for TAM) 48 LOGICAL , PUBLIC :: ln_tfr2d ! logical switch for 2D enhancement (PUBLIC for TAM) 49 43 50 LOGICAL , PUBLIC :: ln_loglayer ! switch for log layer bfr coeff. (PUBLIC for TAM) 44 51 REAL(wp), PUBLIC :: rn_bfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 52 REAL(wp), PUBLIC :: rn_tfrz0 ! bottom roughness for loglayer bfr coeff (PUBLIC for TAM) 45 53 LOGICAL , PUBLIC :: ln_bfrimp ! logical switch for implicit bottom friction 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d ! 2D bottomdrag coefficient (PUBLIC for TAM)54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: bfrcoef2d, tfrcoef2d ! 2D bottom/top drag coefficient (PUBLIC for TAM) 47 55 48 56 !! * Substitutions … … 60 68 !! *** FUNCTION zdf_bfr_alloc *** 61 69 !!---------------------------------------------------------------------- 62 ALLOCATE( bfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc )70 ALLOCATE( bfrcoef2d(jpi,jpj), tfrcoef2d(jpi,jpj), STAT=zdf_bfr_alloc ) 63 71 ! 64 72 IF( lk_mpp ) CALL mpp_sum ( zdf_bfr_alloc ) … … 88 96 INTEGER :: ikbt, ikbu, ikbv ! local integers 89 97 REAL(wp) :: zvu, zuv, zecu, zecv, ztmp ! temporary scalars 90 REAL(wp), POINTER, DIMENSION(:,:) :: zbfrt 98 REAL(wp), POINTER, DIMENSION(:,:) :: zbfrt, ztfrt 91 99 !!---------------------------------------------------------------------- 92 100 ! … … 101 109 IF( nn_bfr == 2 ) THEN ! quadratic bottom friction only 102 110 ! 103 CALL wrk_alloc( jpi, jpj, zbfrt )111 CALL wrk_alloc( jpi, jpj, zbfrt, ztfrt ) 104 112 105 113 IF ( ln_loglayer.AND.lk_vvl ) THEN ! "log layer" bottom friction coefficient … … 120 128 zbfrt(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 121 129 zbfrt(ji,jj) = MIN(zbfrt(ji,jj), rn_bfri2_max) 130 ! (ISF) 131 ikbt = mikt(ji,jj) 132 ! JC: possible WAD implementation should modify line below if layers vanish 133 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_bfrz0 ))**2._wp 134 ztfrt(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 135 ztfrt(ji,jj) = MIN(ztfrt(ji,jj), rn_tfri2_max) 136 122 137 END DO 123 138 END DO … … 125 140 ELSE 126 141 zbfrt(:,:) = bfrcoef2d(:,:) 142 ztfrt(:,:) = tfrcoef2d(:,:) 127 143 ENDIF 128 144 … … 150 166 bfrua(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji+1,jj ) ) * zecu 151 167 bfrva(ji,jj) = - 0.5_wp * ( zbfrt(ji,jj) + zbfrt(ji ,jj+1) ) * zecv 168 ! 169 ! (ISF) ======================================================================== 170 ikbu = miku(ji,jj) ! ocean bottom level at u- and v-points 171 ikbv = mikv(ji,jj) ! (deepest ocean u- and v-points) 172 ! 173 zvu = 0.25 * ( vn(ji,jj ,ikbu) + vn(ji+1,jj ,ikbu) & 174 & + vn(ji,jj-1,ikbu) + vn(ji+1,jj-1,ikbu) ) 175 zuv = 0.25 * ( un(ji,jj ,ikbv) + un(ji-1,jj ,ikbv) & 176 & + un(ji,jj+1,ikbv) + un(ji-1,jj+1,ikbv) ) 177 ! 178 zecu = SQRT( un(ji,jj,ikbu) * un(ji,jj,ikbu) + zvu*zvu + rn_bfeb2 ) 179 zecv = SQRT( vn(ji,jj,ikbv) * vn(ji,jj,ikbv) + zuv*zuv + rn_bfeb2 ) 180 ! 181 tfrua(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji+1,jj ) ) * zecu 182 tfrva(ji,jj) = - 0.5_wp * ( ztfrt(ji,jj) + ztfrt(ji ,jj+1) ) * zecv 183 ! (ISF) END ==================================================================== 184 152 185 END DO 153 186 END DO … … 158 191 IF(ln_ctl) CALL prt_ctl( tab2d_1=bfrua, clinfo1=' bfr - u: ', mask1=umask, & 159 192 & tab2d_2=bfrva, clinfo2= ' v: ', mask2=vmask,ovlap=1 ) 160 CALL wrk_dealloc( jpi,jpj,zbfrt )193 CALL wrk_dealloc( jpi,jpj,zbfrt, ztfrt ) 161 194 ENDIF 162 195 ! … … 183 216 INTEGER :: ios 184 217 REAL(wp) :: zminbfr, zmaxbfr ! temporary scalars 218 REAL(wp) :: zmintfr, zmaxtfr ! temporary scalars 185 219 REAL(wp) :: ztmp, zfru, zfrv ! - - 186 220 !! 187 221 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfri2_max, rn_bfeb2, rn_bfrz0, ln_bfr2d, & 188 & rn_bfrien, ln_bfrimp, ln_loglayer 222 & rn_tfri1, rn_tfri2, rn_tfri2_max, rn_tfeb2, rn_tfrz0, ln_tfr2d, & 223 & rn_bfrien, rn_tfrien, ln_bfrimp, ln_loglayer 189 224 !!---------------------------------------------------------------------- 190 225 ! … … 215 250 bfrua(:,:) = 0.e0 216 251 bfrva(:,:) = 0.e0 252 tfrua(:,:) = 0.e0 253 tfrva(:,:) = 0.e0 217 254 ! 218 255 CASE( 1 ) 219 256 IF(lwp) WRITE(numout,*) ' linear botton friction' 220 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri1 = ', rn_bfri1257 IF(lwp) WRITE(numout,*) ' bottom friction coef. rn_bfri1 = ', rn_bfri1 221 258 IF( ln_bfr2d ) THEN 222 259 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_bfr2d = ', ln_bfr2d 223 260 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 261 ENDIF 262 IF(lwp) WRITE(numout,*) ' top friction coef. rn_bfri1 = ', rn_bfri1 263 IF( ln_tfr2d ) THEN 264 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 265 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 224 266 ENDIF 225 267 ! … … 236 278 bfrua(:,:) = - bfrcoef2d(:,:) 237 279 bfrva(:,:) = - bfrcoef2d(:,:) 280 ! 281 IF(ln_tfr2d) THEN 282 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 283 CALL iom_open('tfr_coef.nc',inum) 284 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! tfrcoef2d is used as tmp array 285 CALL iom_close(inum) 286 tfrcoef2d(:,:) = rn_tfri1 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 287 ELSE 288 tfrcoef2d(:,:) = rn_tfri1 ! initialize tfrcoef2d to the namelist variable 289 ENDIF 290 ! 291 tfrua(:,:) = - tfrcoef2d(:,:) 292 tfrva(:,:) = - tfrcoef2d(:,:) 238 293 ! 239 294 CASE( 2 ) … … 252 307 IF(lwp) WRITE(numout,*) ' coef rn_bfri2 enhancement factor rn_bfrien = ',rn_bfrien 253 308 ENDIF 309 IF(lwp) WRITE(numout,*) ' quadratic top friction' 310 IF(lwp) WRITE(numout,*) ' friction coef. rn_bfri2 = ', rn_tfri2 311 IF(lwp) WRITE(numout,*) ' Max. coef. (log case) rn_tfri2_max = ', rn_tfri2_max 312 IF(lwp) WRITE(numout,*) ' background tke rn_tfeb2 = ', rn_tfeb2 313 IF(lwp) WRITE(numout,*) ' log formulation ln_tfr2d = ', ln_loglayer 314 IF(lwp) WRITE(numout,*) ' bottom roughness rn_tfrz0 [m] = ', rn_tfrz0 315 IF( rn_tfrz0<=0.e0 ) THEN 316 WRITE(ctmp1,*) ' bottom roughness must be strictly positive' 317 CALL ctl_stop( ctmp1 ) 318 ENDIF 319 IF( ln_tfr2d ) THEN 320 IF(lwp) WRITE(numout,*) ' read coef. enhancement distribution from file ln_tfr2d = ', ln_tfr2d 321 IF(lwp) WRITE(numout,*) ' coef rn_tfri2 enhancement factor rn_tfrien = ',rn_tfrien 322 ENDIF 254 323 ! 255 324 IF(ln_bfr2d) THEN … … 263 332 bfrcoef2d(:,:) = rn_bfri2 ! initialize bfrcoef2d to the namelist variable 264 333 ENDIF 334 335 IF(ln_tfr2d) THEN 336 ! tfr_coef is a coefficient in [0,1] giving the mask where to apply the bfr enhancement 337 CALL iom_open('tfr_coef.nc',inum) 338 CALL iom_get (inum, jpdom_data, 'tfr_coef',tfrcoef2d,1) ! bfrcoef2d is used as tmp array 339 CALL iom_close(inum) 340 ! 341 tfrcoef2d(:,:) = rn_tfri2 * ( 1 + rn_tfrien * tfrcoef2d(:,:) ) 342 ELSE 343 tfrcoef2d(:,:) = rn_tfri2 ! initialize tfrcoef2d to the namelist variable 344 ENDIF 265 345 ! 266 346 IF ( ln_loglayer.AND.(.NOT.lk_vvl) ) THEN ! set "log layer" bottom friction once for all … … 279 359 bfrcoef2d(ji,jj) = MAX(bfrcoef2d(ji,jj), ztmp) 280 360 bfrcoef2d(ji,jj) = MIN(bfrcoef2d(ji,jj), rn_bfri2_max) 361 ikbt = mikt(ji,jj) 362 ztmp = tmask(ji,jj,ikbt) * ( vkarmn / LOG( 0.5_wp * fse3t_n(ji,jj,ikbt) / rn_tfrz0 ))**2._wp 363 tfrcoef2d(ji,jj) = MAX(tfrcoef2d(ji,jj), ztmp) 364 tfrcoef2d(ji,jj) = MIN(tfrcoef2d(ji,jj), rn_tfri2_max) 281 365 END DO 282 366 END DO … … 308 392 zminbfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient 309 393 zmaxbfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient 394 zmintfr = 1.e10_wp ! initialise tracker for minimum of bottom friction coefficient 395 zmaxtfr = -1.e10_wp ! initialise tracker for maximum of bottom friction coefficient 310 396 ! 311 397 # if defined key_vectopt_loop … … 339 425 zminbfr = MIN( zminbfr, MIN( zfru, ABS( bfrcoef2d(ji,jj) ) ) ) 340 426 zmaxbfr = MAX( zmaxbfr, MIN( zfrv, ABS( bfrcoef2d(ji,jj) ) ) ) 427 ! (ISF) 428 ikbu = miku(ji,jj) ! deepest ocean level at u- and v-points 429 ikbv = mikv(ji,jj) 430 zfru = 0.5 * fse3u(ji,jj,ikbu) / rdt 431 zfrv = 0.5 * fse3v(ji,jj,ikbv) / rdt 432 IF( ABS( tfrcoef2d(ji,jj) ) > zfru ) THEN 433 IF( ln_ctl ) THEN 434 WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbu 435 WRITE(numout,*) 'BFR ', ABS( tfrcoef2d(ji,jj) ), zfru 436 ENDIF 437 ictu = ictu + 1 438 ENDIF 439 IF( ABS( tfrcoef2d(ji,jj) ) > zfrv ) THEN 440 IF( ln_ctl ) THEN 441 WRITE(numout,*) 'BFR ', narea, nimpp+ji, njmpp+jj, ikbv 442 WRITE(numout,*) 'BFR ', tfrcoef2d(ji,jj), zfrv 443 ENDIF 444 ictv = ictv + 1 445 ENDIF 446 zmintfr = MIN( zmintfr, MIN( zfru, ABS( tfrcoef2d(ji,jj) ) ) ) 447 zmaxtfr = MAX( zmaxtfr, MIN( zfrv, ABS( tfrcoef2d(ji,jj) ) ) ) 448 341 449 END DO 342 450 END DO … … 346 454 CALL mpp_min( zminbfr ) 347 455 CALL mpp_max( zmaxbfr ) 456 CALL mpp_min( zmintfr ) 457 CALL mpp_max( zmaxtfr ) 348 458 ENDIF 349 459 IF( .NOT.ln_bfrimp) THEN … … 352 462 WRITE(numout,*) ' Bottom friction stability check failed at ', ictv, ' V-points ' 353 463 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zminbfr, ' to ', zmaxbfr 464 WRITE(numout,*) ' Bottom friction coefficient now ranges from: ', zmintfr, ' to ', zmaxtfr 354 465 WRITE(numout,*) ' Bottom friction coefficient will be reduced where necessary' 355 466 ENDIF -
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90
r4624 r4666 107 107 108 108 ! ! =============== 109 DO jk = 2, jpkm1! Horizontal slab109 ! ! Horizontal slab 110 110 ! ! =============== 111 DO jj = 1, jpj ! indicators: 112 DO ji = 1, jpi 113 DO jk = mikt(ji,jj)+1, jpkm1 ! Horizontal slab 111 114 ! Define the mask 112 115 ! --------------- 113 rrau(:,:,jk) = MAX( 1.e-20, rrau(:,:,jk) ) ! only retains positive value of rrau 114 115 DO jj = 1, jpj ! indicators: 116 DO ji = 1, jpi 116 rrau(ji,jj,jk) = MAX( 1.e-20, rrau(ji,jj,jk) ) ! only retains positive value of rrau 117 117 118 ! stability indicator: msks=1 if rn2>0; 0 elsewhere 118 119 IF( rn2(ji,jj,jk) + 1.e-12 <= 0. ) THEN ; zmsks(ji,jj) = 0._wp … … 136 137 ELSE ; zmskd3(ji,jj) = 1._wp 137 138 ENDIF 138 END DO139 END DO140 139 ! mask zmsk in order to have avt and avs masked 141 zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk)140 zmsks(:,:) = zmsks(:,:) * tmask(:,:,jk) 142 141 143 142 … … 145 144 ! ------------------ 146 145 ! Constant eddy coefficient: reset to the background value 147 !CDIR NOVERRCHK148 DO jj = 1, jpj149 !CDIR NOVERRCHK150 DO ji = 1, jpi151 146 zinr = 1./rrau(ji,jj,jk) 152 147 ! salt fingering … … 165 160 END DO 166 161 END DO 167 168 162 END DO 169 163 ! Increase avmu, avmv if necessary 170 164 ! -------------------------------- 171 165 !!gm to be changed following the definition of avm. 172 DO jj = 1, jpjm1 173 DO ji = 1, fs_jpim1 ! vector opt. 166 DO jj = 1, jpjm1 167 DO ji = 1, fs_jpim1 ! vector opt. 168 DO jk = miku(ji,jj)+1, jpkm1 ! Horizontal slab 174 169 avmu(ji,jj,jk) = MAX( avmu(ji,jj,jk), & 175 170 & avt(ji,jj,jk), avt(ji+1,jj,jk), & 176 171 & avs(ji,jj,jk), avs(ji+1,jj,jk) ) * umask(ji,jj,jk) 172 END DO 173 DO jk = mikv(ji,jj)+1, jpkm1 177 174 avmv(ji,jj,jk) = MAX( avmv(ji,jj,jk), & 178 175 & avt(ji,jj,jk), avt(ji,jj+1,jk), & -
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
r4624 r4666 241 241 DO jj = 2, jpjm1 ! en(1) = rn_ebb taum / rau0 (min value rn_emin0) 242 242 DO ji = fs_2, fs_jpim1 ! vector opt. 243 en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1) 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 244 248 END DO 245 249 END DO … … 342 346 END DO 343 347 ! 344 DO j k = 2, jpkm1 !* Matrix and right hand side in en345 DO j j = 2, jpjm1346 DO j i = fs_2, fs_jpim1 ! vector opt.348 DO jj = 2, jpjm1 349 DO ji = fs_2, fs_jpim1 ! vector opt. 350 DO jk = mikt(ji,jj)+1, jpkm1 !* Matrix and right hand side in en 347 351 zcof = zfact1 * tmask(ji,jj,jk) 348 352 zzd_up = zcof * ( avm (ji,jj,jk+1) + avm (ji,jj,jk ) ) & ! upper diagonal … … 360 364 ! ! right hand side in en 361 365 en(ji,jj,jk) = en(ji,jj,jk) + rdt * ( zesh2 - avt(ji,jj,jk) * rn2(ji,jj,jk) & 362 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) * tmask(ji,jj,jk) 363 END DO 364 END DO 365 END DO 366 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 367 DO jk = 3, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 368 DO jj = 2, jpjm1 369 DO ji = fs_2, fs_jpim1 ! vector opt. 366 & + zfact3 * dissl(ji,jj,jk) * en (ji,jj,jk) ) & 367 & * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) 368 END DO 369 ! !* Matrix inversion from level 2 (tke prescribed at level 1) 370 DO jk = mikt(ji,jj)+2, jpkm1 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 370 371 zdiag(ji,jj,jk) = zdiag(ji,jj,jk) - zd_lw(ji,jj,jk) * zd_up(ji,jj,jk-1) / zdiag(ji,jj,jk-1) 371 372 END DO 372 END DO 373 END DO 374 DO jj = 2, jpjm1 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 375 DO ji = fs_2, fs_jpim1 ! vector opt. 376 zd_lw(ji,jj,2) = en(ji,jj,2) - zd_lw(ji,jj,2) * en(ji,jj,1) ! Surface boudary conditions on tke 377 END DO 378 END DO 379 DO jk = 3, jpkm1 380 DO jj = 2, jpjm1 381 DO ji = fs_2, fs_jpim1 ! vector opt. 373 ! Second recurrence : Lk = RHSk - Lk / Dk-1 * Lk-1 374 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 375 ! 376 DO jk = mikt(ji,jj)+2, jpkm1 382 377 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) 383 378 END DO 384 END DO 385 END DO 386 DO jj = 2, jpjm1 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 387 DO ji = fs_2, fs_jpim1 ! vector opt. 379 ! 380 ! thrid recurrence : Ek = ( Lk - Uk * Ek+1 ) / Dk 388 381 en(ji,jj,jpkm1) = zd_lw(ji,jj,jpkm1) / zdiag(ji,jj,jpkm1) 389 END DO 390 END DO 391 DO jk = jpk-2, 2, -1 392 DO jj = 2, jpjm1 393 DO ji = fs_2, fs_jpim1 ! vector opt. 382 ! 383 DO jk = jpk-2, mikt(ji,jj)+1, -1 394 384 en(ji,jj,jk) = ( zd_lw(ji,jj,jk) - zd_up(ji,jj,jk) * en(ji,jj,jk+1) ) / zdiag(ji,jj,jk) 395 385 END DO 396 END DO 397 END DO 398 DO jk = 2, jpkm1 ! set the minimum value of tke 399 DO jj = 2, jpjm1 400 DO ji = fs_2, fs_jpim1 ! vector opt. 386 ! 387 DO jk = mikt(ji,jj), jpkm1 ! set the minimum value of tke 401 388 en(ji,jj,jk) = MAX( en(ji,jj,jk), rn_emin ) * tmask(ji,jj,jk) 402 389 END DO … … 412 399 DO ji = fs_2, fs_jpim1 ! vector opt. 413 400 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 414 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 401 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 415 402 END DO 416 403 END DO … … 421 408 jk = nmln(ji,jj) 422 409 en(ji,jj,jk) = en(ji,jj,jk) + rn_efr * en(ji,jj,1) * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 423 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 410 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,1) 424 411 END DO 425 412 END DO … … 433 420 ztx2 = utau(ji-1,jj ) + utau(ji,jj) 434 421 zty2 = vtau(ji ,jj-1) + vtau(ji,jj) 435 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) ! module of the mean stress422 ztau = 0.5_wp * SQRT( ztx2 * ztx2 + zty2 * zty2 ) * tmask(ji,jj,1) ! module of the mean stress 436 423 zdif = taum(ji,jj) - ztau ! mean of modulus - modulus of the mean 437 424 zdif = rhftau_scl * MAX( 0._wp, zdif + rhftau_add ) ! apply some modifications... 438 425 en(ji,jj,jk) = en(ji,jj,jk) + zbbrau * zdif * EXP( -fsdepw(ji,jj,jk) / htau(ji,jj) ) & 439 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) 426 & * ( 1._wp - fr_i(ji,jj) ) * tmask(ji,jj,jk) * tmask(ji,jj,jk-1) * tmask(ji,jj,1) 440 427 END DO 441 428 END DO … … 506 493 ! 507 494 IF( ln_mxl0 ) THEN ! surface mixing length = F(stress) : l=vkarmn*2.e5*taum/(rau0*g) 508 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 509 zmxlm(:,:,1) = MAX( rn_mxl0, zraug * taum(:,:) ) 510 ELSE ! surface set to the minimum value 511 zmxlm(:,:,1) = rn_mxl0 495 DO jj = 2, jpjm1 496 DO ji = fs_2, fs_jpim1 497 IF (mikt(ji,jj) .GT. 1) THEN 498 zmxlm(ji,jj,mikt(ji,jj)) = rmxl_min 499 ELSE 500 zraug = vkarmn * 2.e5_wp / ( rau0 * grav ) 501 zmxlm(ji,jj,mikt(ji,jj)) = MAX( rn_mxl0, zraug * taum(ji,jj) ) 502 END IF 503 END DO 504 END DO 505 ELSE 506 DO jj = 2, jpjm1 507 DO ji = fs_2, fs_jpim1 ! surface set to the minimum value 508 zmxlm(ji,jj,mikt(ji,jj)) = MAX( tmask(ji,jj,1) * rn_mxl0, rmxl_min) 509 END DO 510 END DO 512 511 ENDIF 513 512 zmxlm(:,:,jpk) = rmxl_min ! last level set to the interior minium value 514 513 ! 515 514 !CDIR NOVERRCHK 516 DO j k = 2, jpkm1 ! interior value : l=sqrt(2*e/n^2)517 !CDIR NOVERRCHK 518 DO j j = 2, jpjm1519 !CDIR NOVERRCHK520 DO j i = fs_2, fs_jpim1 ! vector opt.515 DO jj = 2, jpjm1 516 !CDIR NOVERRCHK 517 DO ji = fs_2, fs_jpim1 ! vector opt. 518 !CDIR NOVERRCHK 519 DO jk = mikt(ji,jj)+1, jpkm1 ! interior value : l=sqrt(2*e/n^2) 521 520 zrn2 = MAX( rn2(ji,jj,jk), rsmall ) 522 521 zmxlm(ji,jj,jk) = MAX( rmxl_min, SQRT( 2._wp * en(ji,jj,jk) / zrn2 ) ) 523 522 END DO 523 zmxld(ji,jj,mikt(ji,jj)) = zmxlm(ji,jj,mikt(ji,jj)) ! surface set to the minimum value 524 524 END DO 525 525 END DO … … 533 533 ! 534 534 CASE ( 0 ) ! bounded by the distance to surface and bottom 535 DO j k = 2, jpkm1536 DO j j = 2, jpjm1537 DO j i = fs_2, fs_jpim1 ! vector opt.538 zemxl = MIN( fsdepw(ji,jj,jk) , zmxlm(ji,jj,jk), &535 DO jj = 2, jpjm1 536 DO ji = fs_2, fs_jpim1 ! vector opt. 537 DO jk = mikt(ji,jj)+1, jpkm1 538 zemxl = MIN( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj)), zmxlm(ji,jj,jk), & 539 539 & fsdepw(ji,jj,mbkt(ji,jj)+1) - fsdepw(ji,jj,jk) ) 540 540 zmxlm(ji,jj,jk) = zemxl … … 545 545 ! 546 546 CASE ( 1 ) ! bounded by the vertical scale factor 547 DO j k = 2, jpkm1548 DO j j = 2, jpjm1549 DO j i = fs_2, fs_jpim1 ! vector opt.547 DO jj = 2, jpjm1 548 DO ji = fs_2, fs_jpim1 ! vector opt. 549 DO jk = mikt(ji,jj)+1, jpkm1 550 550 zemxl = MIN( fse3w(ji,jj,jk), zmxlm(ji,jj,jk) ) 551 551 zmxlm(ji,jj,jk) = zemxl … … 556 556 ! 557 557 CASE ( 2 ) ! |dk[xml]| bounded by e3t : 558 DO j k = 2, jpkm1 ! from the surface to the bottom :559 DO j j = 2, jpjm1560 DO j i = fs_2, fs_jpim1 ! vector opt.558 DO jj = 2, jpjm1 559 DO ji = fs_2, fs_jpim1 ! vector opt. 560 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : 561 561 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 562 562 END DO 563 END DO 564 END DO 565 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : 566 DO jj = 2, jpjm1 567 DO ji = fs_2, fs_jpim1 ! vector opt. 563 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : 568 564 zemxl = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 569 565 zmxlm(ji,jj,jk) = zemxl … … 574 570 ! 575 571 CASE ( 3 ) ! lup and ldown, |dk[xml]| bounded by e3t : 576 DO j k = 2, jpkm1 ! from the surface to the bottom : lup577 DO j j = 2, jpjm1578 DO j i = fs_2, fs_jpim1 ! vector opt.572 DO jj = 2, jpjm1 573 DO ji = fs_2, fs_jpim1 ! vector opt. 574 DO jk = mikt(ji,jj)+1, jpkm1 ! from the surface to the bottom : lup 579 575 zmxld(ji,jj,jk) = MIN( zmxld(ji,jj,jk-1) + fse3t(ji,jj,jk-1), zmxlm(ji,jj,jk) ) 580 576 END DO 581 END DO 582 END DO 583 DO jk = jpkm1, 2, -1 ! from the bottom to the surface : ldown 584 DO jj = 2, jpjm1 585 DO ji = fs_2, fs_jpim1 ! vector opt. 577 DO jk = jpkm1, mikt(ji,jj)+1, -1 ! from the bottom to the surface : ldown 586 578 zmxlm(ji,jj,jk) = MIN( zmxlm(ji,jj,jk+1) + fse3t(ji,jj,jk+1), zmxlm(ji,jj,jk) ) 587 579 END DO … … 628 620 CALL lbc_lnk( avm, 'W', 1. ) ! Lateral boundary conditions (sign unchanged) 629 621 ! 630 DO jk = 2, jpkm1 !* vertical eddy viscosity at u- and v-points 622 DO jj = 2, jpjm1 623 DO ji = fs_2, fs_jpim1 ! vector opt. 624 DO jk = miku(ji,jj)+1, jpkm1 !* vertical eddy viscosity at u- and v-points 625 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 626 END DO 627 DO jk = mikv(ji,jj)+1, jpkm1 628 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 629 END DO 630 END DO 631 END DO 632 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 633 ! 634 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 631 635 DO jj = 2, jpjm1 632 636 DO ji = fs_2, fs_jpim1 ! vector opt. 633 avmu(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji+1,jj ,jk) ) * umask(ji,jj,jk) 634 avmv(ji,jj,jk) = 0.5 * ( avm(ji,jj,jk) + avm(ji ,jj+1,jk) ) * vmask(ji,jj,jk) 635 END DO 636 END DO 637 END DO 638 CALL lbc_lnk( avmu, 'U', 1. ) ; CALL lbc_lnk( avmv, 'V', 1. ) ! Lateral boundary conditions 639 ! 640 IF( nn_pdl == 1 ) THEN !* Prandtl number case: update avt 641 DO jk = 2, jpkm1 642 DO jj = 2, jpjm1 643 DO ji = fs_2, fs_jpim1 ! vector opt. 637 DO jk = mikt(ji,jj)+1, jpkm1 644 638 zcoef = avm(ji,jj,jk) * 2._wp * fse3w(ji,jj,jk) * fse3w(ji,jj,jk) 645 639 ! ! shear … … 655 649 avt(ji,jj,jk) = MAX( zpdlr * avt(ji,jj,jk), avtb_2d(ji,jj) * avtb(jk) ) * tmask(ji,jj,jk) 656 650 # if defined key_c1d 657 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) 658 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk)! c1d config. : save Ri651 e_pdl(ji,jj,jk) = zpdlr * tmask(ji,jj,jk) ! c1d configuration : save masked Prandlt number 652 e_ric(ji,jj,jk) = zri * tmask(ji,jj,jk) ! c1d config. : save Ri 659 653 # endif 660 654 END DO -
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90
r4624 r4666 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) 128 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zav_tide(:,:,jk)* tmask(:,:,jk) * tmask(:,:,jk-1) 129 129 END DO 130 130 … … 135 135 END DO 136 136 137 DO jk = 2, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 138 zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s 137 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 138 DO ji = 1, jpi 139 DO jk = mikt(ji,jj)+1, jpkm1 !* Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zav_tide bound by 300 cm2/s 140 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s 141 END DO 142 END DO 139 143 END DO 140 144 … … 162 166 ! ! Update mixing coefs ! 163 167 ! ! ----------------------- ! 164 DO jk = 2, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 165 avt(:,:,jk) = avt(:,:,jk) + zav_tide(:,:,jk) 166 avm(:,:,jk) = avm(:,:,jk) + zav_tide(:,:,jk) 167 DO jj = 2, jpjm1 168 DO ji = fs_2, fs_jpim1 ! vector opt. 168 DO jj = 1, jpj !* Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz to recover en_tmx 169 DO ji = 1, jpi 170 DO jk = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 171 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) 173 END DO 174 END DO 175 END DO 176 177 DO jj = 2, jpjm1 178 DO ji = fs_2, fs_jpim1 ! vector opt. 179 DO jk = mikt(ji,jj)+1, jpkm1 !* update momentum & tracer diffusivity with tidal mixing 169 180 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) 170 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) … … 237 248 zsum2(:,:) = 0.e0 238 249 DO jk= 2, jpk 239 zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) 240 zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) 250 zsum1(:,:) = zsum1(:,:) + zempba_3d_1(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 251 zsum2(:,:) = zsum2(:,:) + zempba_3d_2(:,:,jk) * fse3w(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 241 252 END DO 242 253 DO jj = 1, jpj … … 274 285 zkz(:,:) = 0.e0 ! Associated potential energy consummed over the whole water column 275 286 DO jk = 2, jpkm1 276 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) 287 zkz(:,:) = zkz(:,:) + fse3w(:,:,jk) * MAX( 0.e0, rn2(:,:,jk) ) * rau0 * zavt_itf(:,:,jk) * tmask(:,:,jk) * tmask(:,:,jk-1) 277 288 END DO 278 289 … … 284 295 285 296 DO jk = 2, jpkm1 ! Mutiply by zkz to recover en_tmx, BUT bound by 30/6 ==> zavt_itf bound by 300 cm2/s 286 zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) ! kz max = 120 cm2/s297 zavt_itf(:,:,jk) = zavt_itf(:,:,jk) * MIN( zkz(:,:), 120./10. ) * tmask(:,:,jk) * tmask(:,:,jk-1) ! kz max = 120 cm2/s 287 298 END DO 288 299 … … 377 388 READ ( numnam_cfg, namzdf_tmx, IOSTAT = ios, ERR = 902 ) 378 389 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 379 IF(lwm)WRITE ( numond, namzdf_tmx )390 WRITE ( numond, namzdf_tmx ) 380 391 381 392 IF(lwp) THEN ! Control print … … 414 425 ! only the energy available for mixing is taken into account, 415 426 ! (mixing efficiency tidal dissipation efficiency) 416 en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * tmask(:,:,1)427 en_tmx(:,:) = - rn_tfe * rn_me * ( zem2(:,:) * 1.25 + zek1(:,:) ) * lmask(:,:) 417 428 418 429 ! Vertical structure (az_tmx) … … 443 454 ztpc = 0.e0 444 455 zpc(:,:,:) = MAX(rn_n2min,rn2(:,:,:)) * zav_tide(:,:,:) 445 DO j k= 2, jpkm1446 DO j j = 1, jpj447 DO j i = 1, jpi456 DO jj = 1, jpj 457 DO ji = 1, jpi 458 DO jk= mikt(ji,jj)+1, jpkm1 448 459 ztpc = ztpc + fse3w(ji,jj,jk) * e1t(ji,jj) * e2t(ji,jj) * zpc(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 449 460 END DO … … 459 470 zav_tide(:,:,:) = MIN( zav_tide(:,:,:), 60.e-4 ) 460 471 zkz(:,:) = 0.e0 461 DO j k = 2, jpkm1462 DO jj = 1, jpj463 DO ji = 1, jpi472 DO jj = 1, jpj 473 DO ji = 1, jpi 474 DO jk = mikt(ji,jj)+1, jpkm1 464 475 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) 465 END DO466 END DO476 END DO 477 END DO 467 478 END DO 468 479 ! Here zkz should be equal to en_tmx ==> multiply by en_tmx/zkz … … 484 495 WRITE(numout,*) ' Min de zkz ', ztpc, ' Max = ', maxval(zkz(:,:) ) 485 496 486 DO jk = 2, jpkm1 487 zav_tide(:,:,jk) = zav_tide(:,:,jk) * MIN( zkz(:,:), 30./6. ) !kz max = 300 cm2/s 497 DO jj = 1, jpj 498 DO ji = 1, jpi 499 DO jk = mikt(ji,jj)+1, jpkm1 500 zav_tide(ji,jj,jk) = zav_tide(ji,jj,jk) * MIN( zkz(ji,jj), 30./6. ) !kz max = 300 cm2/s 501 END DO 502 END DO 488 503 END DO 489 504 ztpc = 0.e0
Note: See TracChangeset
for help on using the changeset viewer.
