Changeset 2872
- Timestamp:
- 2011-09-28T10:18:52+02:00 (13 years ago)
- Location:
- branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM
- Files:
-
- 4 added
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynbfr.F90
r2715 r2872 13 13 USE dom_oce ! ocean space and time domain variables 14 14 USE zdf_oce ! ocean vertical physics variables 15 USE zdfbfr ! ocean bottom friction variables 15 16 USE trdmod ! ocean active dynamics and tracers trends 16 17 USE trdmod_oce ! ocean variables trends … … 51 52 !!--------------------------------------------------------------------- 52 53 ! 53 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 54 IF( .not. ln_bfrimp) THEN ! only for explicit bottom friction form 55 ! implicit bfr is implemented in dynzdf_imp 56 ! H. Liu, Sept. 2011 54 57 55 IF( l_trddyn ) THEN ! temporary save of ua and va trends 56 ztrduv(:,:,:,1) = ua(:,:,:) 57 ztrduv(:,:,:,2) = va(:,:,:) 58 ENDIF 58 zm1_2dt = - 1._wp / ( 2._wp * rdt ) 59 60 IF( l_trddyn ) THEN ! temporary save of ua and va trends 61 ztrduv(:,:,:,1) = ua(:,:,:) 62 ztrduv(:,:,:,2) = va(:,:,:) 63 ENDIF 64 59 65 60 66 # if defined key_vectopt_loop 61 DO jj = 1, 162 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)67 DO jj = 1, 1 68 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 63 69 # else 64 DO jj = 2, jpjm165 DO ji = 2, jpim170 DO jj = 2, jpjm1 71 DO ji = 2, jpim1 66 72 # endif 67 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels68 ikbv = mbkv(ji,jj)69 !70 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt)71 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu)72 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv)73 END DO74 END DO73 ikbu = mbku(ji,jj) ! deepest ocean u- & v-levels 74 ikbv = mbkv(ji,jj) 75 ! 76 ! Apply stability criteria on absolute value : abs(bfr/e3) < 1/(2dt) => bfr/e3 > -1/(2dt) 77 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) + MAX( bfrua(ji,jj) / fse3u(ji,jj,ikbu) , zm1_2dt ) * ub(ji,jj,ikbu) 78 va(ji,jj,ikbv) = va(ji,jj,ikbv) + MAX( bfrva(ji,jj) / fse3v(ji,jj,ikbv) , zm1_2dt ) * vb(ji,jj,ikbv) 79 END DO 80 END DO 75 81 76 ! 77 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 78 ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 79 ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 80 CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 81 ENDIF 82 ! ! print mean trends (used for debugging) 83 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 84 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 85 ! 82 ! 83 IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics 84 ztrduv(:,:,:,1) = ua(:,:,:) - ztrduv(:,:,:,1) 85 ztrduv(:,:,:,2) = va(:,:,:) - ztrduv(:,:,:,2) 86 CALL trd_mod( ztrduv(:,:,:,1), ztrduv(:,:,:,2), jpdyn_trd_bfr, 'DYN', kt ) 87 ENDIF 88 ! ! print mean trends (used for debugging) 89 IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' bfr - Ua: ', mask1=umask, & 90 & tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' ) 91 ! 92 ENDIF ! end explicit bottom friction 86 93 END SUBROUTINE dyn_bfr 87 94 -
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r2724 r2872 119 119 INTEGER :: ji, jj, jk, jn ! dummy loop indices 120 120 INTEGER :: icycle ! local scalar 121 INTEGER :: ikbu, ikbv ! local scalar 121 122 REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! local scalars 122 123 REAL(wp) :: z1_8, zx1, zy1 ! - - … … 124 125 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 125 126 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 127 REAL(wp) :: ua_btm, va_btm ! - - 126 128 !!---------------------------------------------------------------------- 127 129 … … 147 149 hvr_e (:,:) = hvr (:,:) 148 150 IF( ln_dynvor_een ) THEN 149 ftne(1,:) = 0. e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0151 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 150 152 DO jj = 2, jpj 151 153 DO ji = fs_2, jpi ! vector opt. 152 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 153 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 154 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 155 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 154 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 155 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 156 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 157 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 156 158 END DO 157 159 END DO … … 161 163 162 164 ! !* Local constant initialization 163 z2dt_b = 2.0 * rdt! baroclinic time step164 z1_8 = 0.5 * 0.25! coefficient for vorticity estimates165 z1_4 = 0.5 * 0.5166 zraur = 1. / rau0! 1 / volumic mass167 ! 168 zhdiv(:,:) = 0. e0! barotropic divergence169 zu_sld = 0. e0 ; zu_asp = 0.e0! tides trends (lk_tide=F)170 zv_sld = 0. e0 ; zv_asp = 0.e0165 z2dt_b = 2.0_wp * rdt ! baroclinic time step 166 z1_8 = 0.5_wp * 0.25_wp ! coefficient for vorticity estimates 167 z1_4 = 0.5_wp * 0.5_wp 168 zraur = 1._wp / rau0 ! 1 / volumic mass 169 ! 170 zhdiv(:,:) = 0._wp ! barotropic divergence 171 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 172 zv_sld = 0._wp ; zv_asp = 0._wp 171 173 172 174 ! ----------------------------------------------------------------------------- … … 176 178 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 177 179 ! ! -------------------------- 178 zua(:,:) = 0. e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0179 zva(:,:) = 0. e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0180 zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp 181 zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp 180 182 ! 181 183 DO jk = 1, jpkm1 … … 205 207 END DO 206 208 207 ! !* baroclinic momentum trend (remove the vertical mean trend)208 DO jk = 1, jpkm1 ! --------------------------209 DO jj = 2, jpjm1210 DO ji = fs_2, fs_jpim1 ! vector opt.211 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)212 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)213 END DO214 END DO215 END DO216 209 217 210 ! !* barotropic Coriolis trends * H (vorticity scheme dependent) … … 260 253 DO jj = 2, jpjm1 261 254 DO ji = fs_2, fs_jpim1 ! vector opt. 262 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 263 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 264 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 265 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 255 ! remove the rhd term according to J. Chanut's suggestion 256 zcu(ji,jj) = zcu(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 257 zcv(ji,jj) = zcv(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 266 258 END DO 267 259 END DO 268 260 ENDIF 261 262 IF(ln_bfrimp.AND.lk_vvl) THEN 263 ! Remove the bottom stress trend from 3-D sea surface level gradient and Coriolis force 264 ! H. Liu, Sept. 2011 265 DO jj = 2, jpjm1 266 DO ji = fs_2, fs_jpim1 267 ikbu = mbku(ji,jj) 268 ikbv = mbkv(ji,jj) 269 ua_btm = zcu(ji,jj) * z2dt_b * hur(ji,jj) * umask (ji,jj,ikbu) 270 va_btm = zcv(ji,jj) * z2dt_b * hvr(ji,jj) * vmask (ji,jj,ikbv) 271 272 ua(ji,jj,ikbu) = ua(ji,jj,ikbu) - bfrua(ji,jj) * ua_btm / fse3u(ji,jj,ikbu) 273 va(ji,jj,ikbv) = va(ji,jj,ikbv) - bfrva(ji,jj) * va_btm / fse3v(ji,jj,ikbv) 274 275 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 276 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 277 END DO 278 END DO 279 END IF 280 ! !* baroclinic momentum trend (remove the vertical mean trend) 281 DO jk = 1, jpkm1 ! -------------------------- 282 DO jj = 2, jpjm1 283 DO ji = fs_2, fs_jpim1 ! vector opt. 284 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 285 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 286 END DO 287 END DO 288 END DO 269 289 270 290 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend … … 279 299 ! ! from the barotropic transport trend 280 300 zcoef = -1. / z2dt_b 301 IF(.NOT.ln_bfrimp) THEN 281 302 # if defined key_vectopt_loop 282 303 DO jj = 1, 1 … … 311 332 END DO 312 333 ENDIF 334 END IF 313 335 314 336 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) … … 410 432 DO ji = fs_2, fs_jpim1 ! vector opt. 411 433 ! surface pressure gradient 412 IF( lk_vvl) THEN 413 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 414 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 415 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 416 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 417 ELSE 418 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 419 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 420 ENDIF 434 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 435 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 421 436 ! energy conserving formulation for planetary vorticity term 422 437 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) … … 427 442 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 428 443 ! after velocities with implicit bottom friction 444 445 IF(ln_bfrimp) THEN 446 ! A new method to implement the implicit bottom friction. 447 ! H. Liu 448 ! May 2010 449 ua_e(ji,jj) = zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) * umask(ji,jj,1) & 450 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 451 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 452 va_e(ji,jj) = zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) * vmask(ji,jj,1) & 453 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 454 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 455 ELSE 429 456 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 430 457 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 431 458 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 432 459 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 460 ENDIF 433 461 END DO 434 462 END DO … … 438 466 DO ji = fs_2, fs_jpim1 ! vector opt. 439 467 ! surface pressure gradient 440 IF( lk_vvl) THEN 441 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 442 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 443 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 444 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 445 ELSE 446 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 447 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 448 ENDIF 468 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 469 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 449 470 ! enstrophy conserving formulation for planetary vorticity term 450 471 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) … … 453 474 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 454 475 ! after velocities with implicit bottom friction 476 IF(ln_bfrimp) THEN 477 ! A new method to implement the implicit bottom friction. 478 ! H. Liu 479 ! May 2010 480 ua_e(ji,jj) = zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) * umask(ji,jj,1) & 481 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 482 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 483 va_e(ji,jj) = zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) * vmask(ji,jj,1) & 484 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 485 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 486 ELSE 455 487 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 456 488 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 457 489 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 458 490 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 491 ENDIF 459 492 END DO 460 493 END DO … … 464 497 DO ji = fs_2, fs_jpim1 ! vector opt. 465 498 ! surface pressure gradient 466 IF( lk_vvl) THEN 467 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 468 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 469 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 470 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 471 ELSE 472 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 473 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 474 ENDIF 499 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 500 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 475 501 ! energy/enstrophy conserving formulation for planetary vorticity term 476 502 zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & … … 479 505 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 480 506 ! after velocities with implicit bottom friction 507 IF(ln_bfrimp) THEN 508 ! A new method to implement the implicit bottom friction. 509 ! H. Liu 510 ! May 2010 511 ua_e(ji,jj) = zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) * umask(ji,jj,1) & 512 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 513 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 514 va_e(ji,jj) = zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) * vmask(ji,jj,1) & 515 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 516 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 517 ELSE 481 518 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 482 519 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 483 520 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 484 521 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 522 ENDIF 485 523 END DO 486 524 END DO -
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/DYN/dynzdf_imp.F90
r2715 r2872 20 20 USE in_out_manager ! I/O manager 21 21 USE lib_mpp ! MPP library 22 USE zdfbfr ! bottom friction setup 22 23 23 24 IMPLICIT NONE … … 61 62 REAL(wp), INTENT(in) :: p2dt ! vertical profile of tracer time-step 62 63 !! 63 INTEGER :: ji, jj, jk ! dummy loop indices 64 REAL(wp) :: z1_p2dt, zcoef, zzwi, zzws, zrhs ! local scalars 64 INTEGER :: ji, jj, jk ! dummy loop indices 65 INTEGER :: ikbum1, ikbvm1 ! local variable 66 REAL(wp) :: z1_p2dt, z2dtf, zcoef, zzwi, zzws, zrhs ! local scalars 67 68 !! * Local variables for implicit bottom friction. H. Liu 69 REAL(wp) :: zbfru, zbfrv 70 REAL(wp) :: bfr_imp = 1._wp 71 !!---------------------------------------------------------------------- 65 72 !!---------------------------------------------------------------------- 66 73 … … 73 80 IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator' 74 81 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 82 IF(.NOT.ln_bfrimp) bfr_imp = 0._wp 75 83 ENDIF 76 84 … … 80 88 81 89 ! 1. Vertical diffusion on u 90 91 ! Vertical diffusion on u&v 82 92 ! --------------------------- 83 93 ! Matrix and second member construction 84 ! bottom boundary condition: both zwi and zws must be masked as avmu can take 85 ! non zero value at the ocean bottom depending on the bottom friction 86 ! used but the bottom velocities have already been updated with the bottom 87 ! friction velocity in dyn_bfr using values from the previous timestep. There 88 ! is no need to include these in the implicit calculation. 89 ! 90 DO jk = 1, jpkm1 ! Matrix 91 DO jj = 2, jpjm1 92 DO ji = fs_2, fs_jpim1 ! vector opt. 94 !! bottom boundary condition: both zwi and zws must be masked as avmu can take 95 !! non zero value at the ocean bottom depending on the bottom friction 96 !! used but the bottom velocities have already been updated with the bottom 97 !! friction velocity in dyn_bfr using values from the previous timestep. There 98 !! is no need to include these in the implicit calculation. 99 100 ! The code has been modified here to implicitly implement bottom 101 ! friction: u(v)mask is not necessary here anymore. 102 ! H. Liu, April 2010. 103 104 ! 1. Vertical diffusion on u 105 DO jj = 2, jpjm1 106 DO ji = fs_2, fs_jpim1 ! vector opt. 107 ikbum1 = mbku(ji,jj) 108 ! Apply stability criteria on absolute value : Min abs(bfr) => Max (bfr) 109 ! zbfru = MAX( bfrua(ji,jj), fse3u(ji,jj,ikbum1)*zinv ) 110 zbfru = bfrua(ji,jj) 111 112 DO jk = 1, ikbum1 93 113 zcoef = - p2dt / fse3u(ji,jj,jk) 94 zzwi = zcoef * avmu (ji,jj,jk ) / fse3uw(ji,jj,jk ) 95 zwi(ji,jj,jk) = zzwi * umask(ji,jj,jk) 96 zzws = zcoef * avmu (ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 97 zws(ji,jj,jk) = zzws * umask(ji,jj,jk+1) 98 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 99 END DO 100 END DO 101 END DO 102 DO jj = 2, jpjm1 ! Surface boudary conditions 103 DO ji = fs_2, fs_jpim1 ! vector opt. 114 zwi(ji,jj,jk) = zcoef * avmu(ji,jj,jk ) / fse3uw(ji,jj,jk ) 115 zws(ji,jj,jk) = zcoef * avmu(ji,jj,jk+1) / fse3uw(ji,jj,jk+1) 116 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 117 END DO 118 119 ! Surface boudary conditions 104 120 zwi(ji,jj,1) = 0._wp 105 121 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 106 END DO 107 END DO 122 123 ! Bottom boudary conditions ! H. Liu, May, 2010 124 ! !commented out to be consisent with v3.2, h.liu 125 ! z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 2.0 * bfr_imp 126 z2dtf = p2dt * zbfru / fse3u(ji,jj,ikbum1) * 1.0_wp * bfr_imp 127 zws(ji,jj,ikbum1) = 0._wp 128 zwd(ji,jj,ikbum1) = 1._wp - zwi(ji,jj,ikbum1) - z2dtf 108 129 109 130 ! Matrix inversion starting from the first level … … 121 142 ! The solution (the after velocity) is in ua 122 143 !----------------------------------------------------------------------- 123 ! 124 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 125 DO jj = 2, jpjm1 126 DO ji = fs_2, fs_jpim1 ! vector opt. 144 145 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 146 DO jk = 2, ikbum1 127 147 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 128 148 END DO 129 END DO 130 END DO 131 ! 132 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 133 DO ji = fs_2, fs_jpim1 ! vector opt. 134 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ( ua(ji,jj,1) + 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) & 135 & / ( fse3u(ji,jj,1) * rau0 ) ) 136 END DO 137 END DO 138 DO jk = 2, jpkm1 139 DO jj = 2, jpjm1 140 DO ji = fs_2, fs_jpim1 ! vector opt. 149 150 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 151 z2dtf = 0.5_wp * p2dt / ( fse3u(ji,jj,1) * rau0 ) 152 ua(ji,jj,1) = ub(ji,jj,1) + p2dt * ua(ji,jj,1) + z2dtf * (utau_b(ji,jj) + utau(ji,jj)) 153 DO jk = 2, ikbum1 141 154 zrhs = ub(ji,jj,jk) + p2dt * ua(ji,jj,jk) ! zrhs=right hand side 142 155 ua(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * ua(ji,jj,jk-1) 143 156 END DO 144 END DO 145 END DO 146 ! 147 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk == 148 DO ji = fs_2, fs_jpim1 ! vector opt. 149 ua(ji,jj,jpkm1) = ua(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 150 END DO 151 END DO 152 DO jk = jpk-2, 1, -1 153 DO jj = 2, jpjm1 154 DO ji = fs_2, fs_jpim1 ! vector opt. 155 ua(ji,jj,jk) = ( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 157 158 159 ! thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk 160 ua(ji,jj,ikbum1) = ua(ji,jj,ikbum1) / zwd(ji,jj,ikbum1) 161 DO jk = ikbum1-1, 1, -1 162 ua(ji,jj,jk) =( ua(ji,jj,jk) - zws(ji,jj,jk) * ua(ji,jj,jk+1) ) / zwd(ji,jj,jk) 156 163 END DO 157 164 END DO … … 170 177 ! 2. Vertical diffusion on v 171 178 ! --------------------------- 172 ! Matrix and second member construction 173 ! bottom boundary condition: both zwi and zws must be masked as avmv can take 174 ! non zero value at the ocean bottom depending on the bottom friction 175 ! used but the bottom velocities have already been updated with the bottom 176 ! friction velocity in dyn_bfr using values from the previous timestep. There 177 ! is no need to include these in the implicit calculation. 178 ! 179 DO jk = 1, jpkm1 ! Matrix 179 180 DO ji = fs_2, fs_jpim1 ! vector opt. 180 181 DO jj = 2, jpjm1 181 DO ji = fs_2, fs_jpim1 ! vector opt. 182 ikbvm1 = mbkv(ji,jj) 183 ! Apply stability criteria on absolute value : Min abs(bfr) => Max (bfr) 184 ! zbfrv = MAX( bfrva(ji,jj), fse3v(ji,jj,ikbvm1)*zinv ) 185 zbfrv = bfrva(ji,jj) 186 187 DO jk = 1, ikbvm1 182 188 zcoef = -p2dt / fse3v(ji,jj,jk) 183 zzwi = zcoef * avmv (ji,jj,jk ) / fse3vw(ji,jj,jk ) 184 zwi(ji,jj,jk) = zzwi * vmask(ji,jj,jk) 185 zzws = zcoef * avmv (ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 186 zws(ji,jj,jk) = zzws * vmask(ji,jj,jk+1) 187 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zzws 188 END DO 189 END DO 190 END DO 191 DO jj = 2, jpjm1 ! Surface boudary conditions 192 DO ji = fs_2, fs_jpim1 ! vector opt. 189 zwi(ji,jj,jk) = zcoef * avmv(ji,jj,jk ) / fse3vw(ji,jj,jk ) 190 zws(ji,jj,jk) = zcoef * avmv(ji,jj,jk+1) / fse3vw(ji,jj,jk+1) 191 zwd(ji,jj,jk) = 1._wp - zwi(ji,jj,jk) - zws(ji,jj,jk) 192 END DO 193 194 ! Surface boudary conditions 193 195 zwi(ji,jj,1) = 0._wp 194 196 zwd(ji,jj,1) = 1._wp - zws(ji,jj,1) 195 END DO 196 END DO 197 198 ! Bottom boudary conditions 199 ! Bottom boudary conditions ! H. Liu, May, 2010 200 ! !commented out to be consisent with v3.2, h.liu 201 ! z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 2.0 * bfr_imp 202 z2dtf = p2dt * zbfrv / fse3v(ji,jj,ikbvm1) * 1.0_wp * bfr_imp 203 zws(ji,jj,ikbvm1) = 0._wp 204 zwd(ji,jj,ikbvm1) = 1._wp - zwi(ji,jj,ikbvm1) - z2dtf 197 205 198 206 ! Matrix inversion … … 206 214 ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk ) 207 215 ! 208 ! m is decomposed in the product of an upper and lower triangular matrix 216 ! m is decomposed in the product of an upper and lower triangular 217 ! matrix 209 218 ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi 210 219 ! The solution (after velocity) is in 2d array va 211 220 !----------------------------------------------------------------------- 212 ! 213 DO jk = 2, jpkm1 !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) == 214 DO jj = 2, jpjm1 215 DO ji = fs_2, fs_jpim1 ! vector opt. 221 222 ! First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) 223 DO jk = 2, ikbvm1 216 224 zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1) 217 225 END DO 218 END DO 219 END DO 220 ! 221 DO jj = 2, jpjm1 !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 == 222 DO ji = fs_2, fs_jpim1 ! vector opt. 223 va(ji,jj,1) = vb(ji,jj,1) + p2dt * ( va(ji,jj,1) + 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) & 224 & / ( fse3v(ji,jj,1) * rau0 ) ) 225 END DO 226 END DO 227 DO jk = 2, jpkm1 228 DO jj = 2, jpjm1 229 DO ji = fs_2, fs_jpim1 ! vector opt. 226 227 ! second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 228 z2dtf = 0.5_wp * p2dt / ( fse3v(ji,jj,1)*rau0 ) 229 va(ji,jj,1) = vb(ji,jj,1) + p2dt * va(ji,jj,1) + z2dtf * (vtau_b(ji,jj) + vtau(ji,jj)) 230 DO jk = 2, ikbvm1 230 231 zrhs = vb(ji,jj,jk) + p2dt * va(ji,jj,jk) ! zrhs=right hand side 231 232 va(ji,jj,jk) = zrhs - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * va(ji,jj,jk-1) 232 233 END DO 233 END DO 234 END DO 235 ! 236 DO jj = 2, jpjm1 !== thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk == 237 DO ji = fs_2, fs_jpim1 ! vector opt. 238 va(ji,jj,jpkm1) = va(ji,jj,jpkm1) / zwd(ji,jj,jpkm1) 239 END DO 240 END DO 241 DO jk = jpk-2, 1, -1 242 DO jj = 2, jpjm1 243 DO ji = fs_2, fs_jpim1 ! vector opt. 244 va(ji,jj,jk) = ( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 245 END DO 234 235 ! thrid recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk 236 va(ji,jj,ikbvm1) = va(ji,jj,ikbvm1) / zwd(ji,jj,ikbvm1) 237 238 DO jk = ikbvm1-1, 1, -1 239 va(ji,jj,jk) =( va(ji,jj,jk) - zws(ji,jj,jk) * va(ji,jj,jk+1) ) / zwd(ji,jj,jk) 240 END DO 241 246 242 END DO 247 243 END DO … … 258 254 IF( wrk_not_released(3, 3) ) CALL ctl_stop('dyn_zdf_imp: failed to release workspace array') 259 255 ! 256 260 257 END SUBROUTINE dyn_zdf_imp 261 258 -
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfbfr.F90
r2715 r2872 36 36 REAL(wp) :: rn_bfrien = 30._wp ! local factor to enhance coefficient bfri 37 37 LOGICAL :: ln_bfr2d = .false. ! logical switch for 2D enhancement 38 LOGICAL, PUBLIC :: ln_bfrimp = .false. ! logical switch for implicit bottom friction 38 39 39 40 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: bfrcoef2d ! 2D bottom drag coefficient … … 142 143 REAL(wp) :: zfru, zfrv ! - - 143 144 !! 144 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien 145 NAMELIST/nambfr/ nn_bfr, rn_bfri1, rn_bfri2, rn_bfeb2, ln_bfr2d, rn_bfrien, ln_bfrimp 145 146 !!---------------------------------------------------------------------- 146 147 -
branches/2011/dev_r2802_NOCL_bfrimp/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
r2715 r2872 81 81 !!--------------------------------------------------------------------- 82 82 # include "par_POMME_R025.h90" 83 #elif defined key_amm12 84 !!--------------------------------------------------------------------- 85 !! 'key_amm12' : Atlantic Margin Model (~12km) : AMM 86 !!--------------------------------------------------------------------- 87 # include "par_AMM12.h90" 88 #elif defined key_amm7 89 !!--------------------------------------------------------------------- 90 !! 'key_amm7' : Atlantic Margin Model (~7km) : AMM 91 !!--------------------------------------------------------------------- 92 # include "par_AMM7.h90" 83 93 #else 84 94 !!---------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.