Changeset 5974 for branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2015-12-02T11:52:05+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5072_UKMO2_OBS_simplification/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5682 r5974 11 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 13 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 14 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts || defined key_esopa15 15 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts'split explicit free surface16 !! split explicit free surface 17 17 !!---------------------------------------------------------------------- 18 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 23 23 USE sbc_oce ! surface boundary condition: ocean 24 24 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables26 25 USE phycst ! physical constants 27 26 USE dynvor ! vorticity term 28 27 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 28 USE bdytides ! open boundary condition data 30 29 USE bdydyn2d ! open boundary conditions on barotropic variables 31 30 USE sbctide ! tides … … 68 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 68 70 ! Arrays below are saved to allow testing of the "no time averaging" option71 ! If this option is not retained, these could be replaced by temporary arrays72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays73 ubb_e, ub_e, &74 vbb_e, vb_e75 76 69 !! * Substitutions 77 70 # include "domzgr_substitute.h90" … … 88 81 !! *** routine dyn_spg_ts_alloc *** 89 82 !!---------------------------------------------------------------------- 90 INTEGER :: ierr( 3)83 INTEGER :: ierr(4) 91 84 !!---------------------------------------------------------------------- 92 85 ierr(:) = 0 93 86 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 87 ALLOCATE( ssha_e(jpi,jpj), sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 88 & ua_e(jpi,jpj), un_e(jpi,jpj), ub_e(jpi,jpj), ubb_e(jpi,jpj), & 89 & va_e(jpi,jpj), vn_e(jpi,jpj), vb_e(jpi,jpj), vbb_e(jpi,jpj), & 90 & hu_e(jpi,jpj), hur_e(jpi,jpj), hv_e(jpi,jpj), hvr_e(jpi,jpj), STAT= ierr(1) ) 97 91 98 92 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 99 93 100 IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 94 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 95 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 96 97 ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_adv(jpi,jpj), vn_adv(jpi,jpj), & 98 #if defined key_agrif 99 & ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj) , & 100 #endif 101 & STAT= ierr(4)) 102 102 103 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 104 105 105 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 107 ! 108 108 END FUNCTION dyn_spg_ts_alloc 109 109 110 110 111 SUBROUTINE dyn_spg_ts( kt ) … … 145 146 INTEGER :: ikbu, ikbv, noffset ! local integers 146 147 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 149 REAL(wp) :: zu_spg, zv_spg ! - -150 REAL(wp) :: zhura, zhvra 151 REAL(wp) :: za0, za1, za2, za3 152 ! 153 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e148 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 149 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 150 REAL(wp) :: zu_spg, zv_spg ! - - 151 REAL(wp) :: zhura, zhvra ! - - 152 REAL(wp) :: za0, za1, za2, za3 ! - - 153 ! 154 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 154 155 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv156 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 157 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 158 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 163 164 ! !* Allocate temporary arrays 164 165 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 165 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)166 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 167 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 167 168 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 168 169 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 187 188 ! 188 189 ! time offset in steps for bdy data update 189 IF (.NOT.ln_bt_fw) THEN ; noffset=- 2*nn_baro ; ELSE ; noffset = 0 ; ENDIF190 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 190 191 ! 191 192 IF( kt == nit000 ) THEN !* initialisation … … 219 220 ! 220 221 IF ( kt == nit000 .OR. lk_vvl ) THEN 221 IF ( ln_dynvor_een_old ) THEN 222 DO jj = 1, jpjm1 223 DO ji = 1, jpim1 224 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 225 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 226 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 227 END DO 228 END DO 222 IF ( ln_dynvor_een ) THEN !== EEN scheme ==! 223 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 224 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 225 DO jj = 1, jpjm1 226 DO ji = 1, jpim1 227 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 228 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 229 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 230 END DO 231 END DO 232 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 233 DO jj = 1, jpjm1 234 DO ji = 1, jpim1 235 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 236 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 237 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 238 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 239 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 240 END DO 241 END DO 242 END SELECT 229 243 CALL lbc_lnk( zwz, 'F', 1._wp ) 230 zwz(:,:) = ff(:,:) * zwz(:,:) 231 244 ! 232 245 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 233 246 DO jj = 2, jpj 234 DO ji = fs_2, jpi ! vector opt.247 DO ji = 2, jpi 235 248 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 236 249 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 239 252 END DO 240 253 END DO 241 ELSE IF ( ln_dynvor_een ) THEN 242 DO jj = 1, jpjm1 243 DO ji = 1, jpim1 244 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 245 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 246 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 247 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 248 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 249 END DO 250 END DO 251 CALL lbc_lnk( zwz, 'F', 1._wp ) 252 zwz(:,:) = ff(:,:) * zwz(:,:) 253 254 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 255 DO jj = 2, jpj 256 DO ji = fs_2, jpi ! vector opt. 257 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 258 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 259 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 260 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 261 END DO 262 END DO 263 ELSE 254 ! 255 ELSE !== all other schemes (ENE, ENS, MIX) 264 256 zwz(:,:) = 0._wp 265 zhf(:,:) = 0. 257 zhf(:,:) = 0._wp 266 258 IF ( .not. ln_sco ) THEN 259 260 !!gm agree the JC comment : this should be done in a much clear way 261 267 262 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 268 263 ! Set it to zero for the time being … … 276 271 277 272 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) *(1._wp- umask(:,jj,1) * umask(:,jj+1,1))273 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 279 274 END DO 280 275 … … 297 292 ! If forward start at previous time step, and centered integration, 298 293 ! then update averaging weights: 299 IF ( (.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN294 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 300 295 ll_fw_start=.FALSE. 301 296 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) … … 338 333 DO jj = 2, jpjm1 339 334 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)341 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)342 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) /e2v(ji,jj)343 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) /e2v(ji,jj)335 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 336 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 337 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 338 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 344 339 ! energy conserving formulation for planetary vorticity term 345 340 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) … … 352 347 DO ji = fs_2, fs_jpim1 ! vector opt. 353 348 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 354 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)349 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 355 350 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 356 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)351 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 357 352 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 358 353 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 360 355 END DO 361 356 ! 362 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN ! enstrophy and energy conserving scheme357 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 363 358 DO jj = 2, jpjm1 364 359 DO ji = fs_2, fs_jpim1 ! vector opt. 365 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &366 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &367 & + ftse(ji,jj ) * zwy(ji ,jj-1) &368 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )369 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &370 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &371 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &372 & + ftne(ji,jj ) * zwx(ji ,jj ) )360 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 361 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 362 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 363 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 364 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 365 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 366 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 367 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 373 368 END DO 374 369 END DO … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) /e1u(ji,jj)384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) /e2v(ji,jj)378 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 379 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 385 380 END DO 386 381 END DO … … 431 426 DO jj = 2, jpjm1 432 427 DO ji = fs_2, fs_jpim1 ! vector opt. 433 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)434 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)428 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 429 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 435 430 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 436 431 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 441 436 DO ji = fs_2, fs_jpim1 ! vector opt. 442 437 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 443 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)438 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 444 439 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 445 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)440 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 446 441 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 447 442 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 490 485 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 486 sshn_e(:,:) = sshn (:,:) 492 zun_e(:,:) = un_b (:,:)493 zvn_e(:,:) = vn_b (:,:)487 un_e (:,:) = un_b (:,:) 488 vn_e (:,:) = vn_b (:,:) 494 489 ! 495 490 hu_e (:,:) = hu (:,:) … … 499 494 ELSE ! CENTRED integration: start from BEFORE fields 500 495 sshn_e(:,:) = sshb (:,:) 501 zun_e(:,:) = ub_b (:,:)502 zvn_e(:,:) = vb_b (:,:)496 un_e (:,:) = ub_b (:,:) 497 vn_e (:,:) = vb_b (:,:) 503 498 ! 504 499 hu_e (:,:) = hu_b (:,:) … … 514 509 va_b (:,:) = 0._wp 515 510 ssha (:,:) = 0._wp ! Sum for after averaged sea level 516 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop517 zv_sum(:,:) = 0._wp511 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 512 vn_adv(:,:) = 0._wp 518 513 ! ! ==================== ! 519 514 DO jn = 1, icycle ! sub-time-step loop ! … … 523 518 ! Update only tidal forcing at open boundaries 524 519 #if defined key_tide 525 IF ( lk_bdy .AND. lk_tide ) 526 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )520 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 521 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 527 522 #endif 528 523 ! … … 539 534 540 535 ! Extrapolate barotropic velocities at step jit+0.5: 541 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)542 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)536 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 537 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 543 538 544 539 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 549 544 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 550 545 DO ji = 2, fs_jpim1 ! Vector opt. 551 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &552 & * ( e1 2t(ji ,jj) * zsshp2_e(ji ,jj) &553 & + e1 2t(ji+1,jj) * zsshp2_e(ji+1,jj) )554 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &555 & * ( e1 2t(ji,jj ) * zsshp2_e(ji,jj ) &556 & + e1 2t(ji,jj+1) * zsshp2_e(ji,jj+1) )546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 547 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 550 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 557 552 END DO 558 553 END DO … … 602 597 ! Sum over sub-time-steps to compute advective velocities 603 598 za2 = wgtbtp2(jn) 604 zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u(:,:)605 zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v(:,:)599 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 600 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 606 601 ! 607 602 ! Set next sea level: … … 609 604 DO ji = fs_2, fs_jpim1 ! vector opt. 610 605 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 611 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1 2t(ji,jj)606 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 612 607 END DO 613 608 END DO … … 627 622 DO jj = 2, jpjm1 628 623 DO ji = 2, jpim1 ! NO Vector Opt. 629 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1 2u(ji,jj) &630 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &631 & + e1 2t(ji+1,jj ) * ssha_e(ji+1,jj ) )632 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1 2v(ji,jj) &633 & * ( e1 2t(ji ,jj ) * ssha_e(ji ,jj ) &634 & + e1 2t(ji ,jj+1) * ssha_e(ji ,jj+1) )624 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 625 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 626 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 627 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 628 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 629 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 635 630 END DO 636 631 END DO … … 666 661 DO jj = 2, jpjm1 667 662 DO ji = 2, jpim1 668 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1 2u(ji ,jj) &669 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj) &670 & + e1 2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )671 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1 2v(ji ,jj ) &672 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj ) &673 & + e1 2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) & 664 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) & 667 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 674 669 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 675 670 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 … … 688 683 DO jj = 2, jpjm1 689 684 DO ji = fs_2, fs_jpim1 ! vector opt. 690 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)691 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)692 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) /e2v(ji,jj)693 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)685 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 686 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 687 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 688 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 694 689 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 695 690 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) … … 701 696 DO ji = fs_2, fs_jpim1 ! vector opt. 702 697 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)698 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 704 699 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)700 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 706 701 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 702 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 709 704 END DO 710 705 ! 711 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN !== energy and enstrophy conserving scheme ==!706 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 712 707 DO jj = 2, jpjm1 713 708 DO ji = fs_2, fs_jpim1 ! vector opt. 714 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &715 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &716 & + ftse(ji,jj ) * zwy(ji ,jj-1) &717 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )718 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &719 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &720 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &721 & + ftne(ji,jj ) * zwx(ji ,jj ) )709 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 710 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 711 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 712 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 713 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 714 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 715 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 716 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 722 717 END DO 723 718 END DO … … 729 724 DO jj = 2, jpjm1 730 725 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)732 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)726 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 727 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 733 728 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 734 729 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg … … 738 733 ! 739 734 ! Add bottom stresses: 740 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)741 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)735 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 736 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 742 737 ! 743 738 ! Surface pressure trend: … … 745 740 DO ji = fs_2, fs_jpim1 ! vector opt. 746 741 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) /e1u(ji,jj)748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) /e2v(ji,jj)742 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 743 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 749 744 zwx(ji,jj) = zu_spg 750 745 zwy(ji,jj) = zv_spg … … 756 751 DO jj = 2, jpjm1 757 752 DO ji = fs_2, fs_jpim1 ! vector opt. 758 ua_e(ji,jj) = ( zun_e(ji,jj) &753 ua_e(ji,jj) = ( un_e(ji,jj) & 759 754 & + rdtbt * ( zwx(ji,jj) & 760 755 & + zu_trd(ji,jj) & … … 762 757 & ) * umask(ji,jj,1) 763 758 764 va_e(ji,jj) = ( zvn_e(ji,jj) &759 va_e(ji,jj) = ( vn_e(ji,jj) & 765 760 & + rdtbt * ( zwy(ji,jj) & 766 761 & + zv_trd(ji,jj) & … … 777 772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 773 779 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &774 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 780 775 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 776 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 783 778 & ) * zhura 784 779 785 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &780 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 786 781 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 787 782 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 807 802 #if defined key_bdy 808 803 ! open boundaries 809 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )804 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 810 805 #endif 811 806 #if defined key_agrif … … 815 810 ! ! ---- 816 811 ubb_e (:,:) = ub_e (:,:) 817 ub_e (:,:) = zun_e(:,:)818 zun_e(:,:) = ua_e (:,:)812 ub_e (:,:) = un_e (:,:) 813 un_e (:,:) = ua_e (:,:) 819 814 ! 820 815 vbb_e (:,:) = vb_e (:,:) 821 vb_e (:,:) = zvn_e(:,:)822 zvn_e(:,:) = va_e (:,:)816 vb_e (:,:) = vn_e (:,:) 817 vn_e (:,:) = va_e (:,:) 823 818 ! 824 819 sshbb_e(:,:) = sshb_e(:,:) … … 845 840 ! ----------------------------------------------------------------------------- 846 841 ! 847 ! At this stage ssha holds a time averaged value848 ! ! Sea Surface Height at u-,v- and f-points849 IF( lk_vvl ) THEN ! (required only in key_vvl case)850 DO jj = 1, jpjm1851 DO ji = 1, jpim1 ! NO Vector Opt.852 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &853 & * ( e12t(ji ,jj) * ssha(ji ,jj) &854 & + e12t(ji+1,jj) * ssha(ji+1,jj) )855 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &856 & * ( e12t(ji,jj ) * ssha(ji,jj ) &857 & + e12t(ji,jj+1) * ssha(ji,jj+1) )858 END DO859 END DO860 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions861 ENDIF862 !863 842 ! Set advection velocity correction: 843 zwx(:,:) = un_adv(:,:) 844 zwy(:,:) = vn_adv(:,:) 864 845 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 865 un_adv(:,:) = z u_sum(:,:)*hur(:,:)866 vn_adv(:,:) = z v_sum(:,:)*hvr(:,:)846 un_adv(:,:) = zwx(:,:)*hur(:,:) 847 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 867 848 ELSE 868 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + z u_sum(:,:)) * hur(:,:)869 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + z v_sum(:,:)) * hvr(:,:)849 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 850 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 870 851 END IF 871 852 872 853 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 873 ub2_b(:,:) = z u_sum(:,:)874 vb2_b(:,:) = z v_sum(:,:)854 ub2_b(:,:) = zwx(:,:) 855 vb2_b(:,:) = zwy(:,:) 875 856 ENDIF 876 857 ! … … 882 863 END DO 883 864 ELSE 865 ! At this stage, ssha has been corrected: compute new depths at velocity points 866 DO jj = 1, jpjm1 867 DO ji = 1, jpim1 ! NO Vector Opt. 868 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 869 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 870 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 871 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 872 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 873 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 874 END DO 875 END DO 876 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 877 ! 884 878 DO jk=1,jpkm1 885 879 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 920 914 ! 921 915 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 922 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)923 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )916 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 917 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 924 918 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 925 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1069 1063 ! 1070 1064 INTEGER :: ji ,jj 1071 INTEGER :: ios ! Local integer output status for namelist read1072 1065 REAL(wp) :: zxr2, zyr2, zcmax 1073 1066 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1074 1067 !! 1075 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1076 & nn_baro, rn_bt_cmax, nn_bt_flt1077 1068 !!---------------------------------------------------------------------- 1078 1069 ! 1079 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1080 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1081 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1082 1083 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1084 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1085 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1086 IF(lwm) WRITE ( numond, namsplit ) 1087 ! 1088 ! ! Max courant number for ext. grav. waves 1070 ! Max courant number for ext. grav. waves 1089 1071 ! 1090 1072 CALL wrk_alloc( jpi, jpj, zcu ) 1091 1073 ! 1092 IF (lk_vvl) THEN 1093 DO jj = 1, jpj 1094 DO ji =1, jpi 1095 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1096 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1097 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1098 END DO 1099 END DO 1100 ELSE 1101 DO jj = 1, jpj 1102 DO ji =1, jpi 1103 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1104 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1105 zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) ) 1106 END DO 1107 END DO 1108 ENDIF 1109 1110 zcmax = MAXVAL(zcu(:,:)) 1074 DO jj = 1, jpj 1075 DO ji =1, jpi 1076 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1077 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1078 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1079 END DO 1080 END DO 1081 ! 1082 zcmax = MAXVAL( zcu(:,:) ) 1111 1083 IF( lk_mpp ) CALL mpp_max( zcmax ) 1112 1084 1113 1085 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1114 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1086 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1115 1087 1116 rdtbt = rdt / FLOAT(nn_baro)1088 rdtbt = rdt / REAL( nn_baro , wp ) 1117 1089 zcmax = zcmax * rdtbt 1118 1090 ! Print results … … 1120 1092 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1121 1093 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1122 IF( ln_bt_ nn_auto ) THEN1123 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1094 IF( ln_bt_auto ) THEN 1095 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1124 1096 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1125 1097 ELSE 1126 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1098 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1127 1099 ENDIF 1128 1100 … … 1169 1141 END SUBROUTINE dyn_spg_ts_init 1170 1142 1171 #else1172 !!---------------------------------------------------------------------------1173 !! Default case : Empty module No split explicit free surface1174 !!---------------------------------------------------------------------------1175 CONTAINS1176 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1177 dyn_spg_ts_alloc = 01178 END FUNCTION dyn_spg_ts_alloc1179 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1180 INTEGER, INTENT(in) :: kt1181 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1182 END SUBROUTINE dyn_spg_ts1183 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1184 INTEGER , INTENT(in) :: kt ! ocean time-step1185 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1186 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1187 END SUBROUTINE ts_rst1188 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1189 INTEGER , INTENT(in) :: kt ! ocean time-step1190 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1191 END SUBROUTINE dyn_spg_ts_init1192 #endif1193 1194 1143 !!====================================================================== 1195 1144 END MODULE dynspg_ts 1196 1197 1198
Note: See TracChangeset
for help on using the changeset viewer.