- Timestamp:
- 2015-11-30T20:55:41+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5944 r5956 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 ) … … 146 147 INTEGER :: iktu, iktv ! local integers 147 148 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 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(:,:) :: z un_e, zvn_e, zsshp2_e149 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 150 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 151 REAL(wp) :: zu_spg, zv_spg ! - - 152 REAL(wp) :: zhura, zhvra ! - - 153 REAL(wp) :: za0, za1, za2, za3 ! - - 154 ! 155 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 155 156 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 156 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv157 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 157 158 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 158 159 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a … … 164 165 ! !* Allocate temporary arrays 165 166 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 166 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)167 CALL wrk_alloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc)167 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd) 168 CALL wrk_alloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 168 169 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 169 170 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) … … 188 189 ! 189 190 ! time offset in steps for bdy data update 190 IF (.NOT.ln_bt_fw) THEN ; noffset=- 2*nn_baro ; ELSE ; noffset = 0 ; ENDIF191 IF (.NOT.ln_bt_fw) THEN ; noffset=-nn_baro ; ELSE ; noffset = 0 ; ENDIF 191 192 ! 192 193 IF( kt == nit000 ) THEN !* initialisation … … 220 221 ! 221 222 IF ( kt == nit000 .OR. lk_vvl ) THEN 222 IF ( ln_dynvor_een_old ) THEN 223 DO jj = 1, jpjm1 224 DO ji = 1, jpim1 225 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 226 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 227 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 228 END DO 229 END DO 223 IF ( ln_dynvor_een ) THEN !== EEN scheme ==! 224 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 225 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 226 DO jj = 1, jpjm1 227 DO ji = 1, jpim1 228 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 229 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 230 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 231 END DO 232 END DO 233 END DO 234 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 235 DO jj = 1, jpjm1 236 DO ji = 1, jpim1 237 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 238 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 239 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 240 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 241 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 242 END DO 243 END DO 244 END SELECT 230 245 CALL lbc_lnk( zwz, 'F', 1._wp ) 231 246 zwz(:,:) = ff(:,:) * zwz(:,:) … … 233 248 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 234 249 DO jj = 2, jpj 235 DO ji = fs_2, jpi ! vector opt.250 DO ji = 2, jpi 236 251 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 237 252 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 240 255 END DO 241 256 END DO 242 ELSE IF ( ln_dynvor_een ) THEN 243 DO jj = 1, jpjm1 244 DO ji = 1, jpim1 245 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 246 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 247 & / ( MAX( 1.0_wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) + & 248 & ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 249 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 250 END DO 251 END DO 252 CALL lbc_lnk( zwz, 'F', 1._wp ) 253 zwz(:,:) = ff(:,:) * zwz(:,:) 254 255 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 256 DO jj = 2, jpj 257 DO ji = fs_2, jpi ! vector opt. 258 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 259 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 260 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 261 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 262 END DO 263 END DO 264 ELSE 257 ! 258 ELSE !== all other schemes (ENE, ENS, MIX) 265 259 zwz(:,:) = 0._wp 266 zhf(:,:) = 0. 260 zhf(:,:) = 0._wp 267 261 IF ( .not. ln_sco ) THEN 262 263 !!gm agree the JC comment : this should be done in a much clear way 264 268 265 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 269 266 ! Set it to zero for the time being … … 277 274 278 275 DO jj = 1, jpjm1 279 zhf(:,jj) = zhf(:,jj) *(1._wp- umask(:,jj,1) * umask(:,jj+1,1))276 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 280 277 END DO 281 278 … … 298 295 ! If forward start at previous time step, and centered integration, 299 296 ! then update averaging weights: 300 IF ( (.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN297 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 301 298 ll_fw_start=.FALSE. 302 299 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) … … 339 336 DO jj = 2, jpjm1 340 337 DO ji = fs_2, fs_jpim1 ! vector opt. 341 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)342 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)343 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) /e2v(ji,jj)344 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) /e2v(ji,jj)338 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 339 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 340 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 341 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 345 342 ! energy conserving formulation for planetary vorticity term 346 343 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) … … 353 350 DO ji = fs_2, fs_jpim1 ! vector opt. 354 351 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 355 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)352 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 356 353 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 357 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)354 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 358 355 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 359 356 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 361 358 END DO 362 359 ! 363 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN ! enstrophy and energy conserving scheme360 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 364 361 DO jj = 2, jpjm1 365 362 DO ji = fs_2, fs_jpim1 ! vector opt. 366 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &367 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &368 & + ftse(ji,jj ) * zwy(ji ,jj-1) &369 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )370 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &371 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &372 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &373 & + ftne(ji,jj ) * zwx(ji ,jj ) )363 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 364 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 365 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 366 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 367 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 368 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 369 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 370 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 374 371 END DO 375 372 END DO … … 382 379 DO jj = 2, jpjm1 383 380 DO ji = fs_2, fs_jpim1 ! vector opt. 384 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) /e1u(ji,jj)385 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) /e2v(ji,jj)381 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 382 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 386 383 END DO 387 384 END DO … … 457 454 DO jj = 2, jpjm1 458 455 DO ji = fs_2, fs_jpim1 ! vector opt. 459 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)460 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)456 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 457 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 461 458 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 462 459 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 467 464 DO ji = fs_2, fs_jpim1 ! vector opt. 468 465 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 469 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)466 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 470 467 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 471 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)468 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 472 469 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 473 470 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 491 488 ENDIF 492 489 #endif 493 ! !* Fill boundary data arrays withAGRIF494 ! ! ------------------------------------ -490 ! !* Fill boundary data arrays for AGRIF 491 ! ! ------------------------------------ 495 492 #if defined key_agrif 496 493 IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) … … 516 513 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 517 514 sshn_e(:,:) = sshn (:,:) 518 zun_e(:,:) = un_b (:,:)519 zvn_e(:,:) = vn_b (:,:)515 un_e (:,:) = un_b (:,:) 516 vn_e (:,:) = vn_b (:,:) 520 517 ! 521 518 hu_e (:,:) = hu (:,:) … … 525 522 ELSE ! CENTRED integration: start from BEFORE fields 526 523 sshn_e(:,:) = sshb (:,:) 527 zun_e(:,:) = ub_b (:,:)528 zvn_e(:,:) = vb_b (:,:)524 un_e (:,:) = ub_b (:,:) 525 vn_e (:,:) = vb_b (:,:) 529 526 ! 530 527 hu_e (:,:) = hu_b (:,:) … … 540 537 va_b (:,:) = 0._wp 541 538 ssha (:,:) = 0._wp ! Sum for after averaged sea level 542 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop543 zv_sum(:,:) = 0._wp539 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 540 vn_adv(:,:) = 0._wp 544 541 ! ! ==================== ! 545 542 DO jn = 1, icycle ! sub-time-step loop ! … … 549 546 ! Update only tidal forcing at open boundaries 550 547 #if defined key_tide 551 IF ( lk_bdy .AND. lk_tide ) 552 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset )548 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 549 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, time_offset=noffset ) 553 550 #endif 554 551 ! … … 565 562 566 563 ! Extrapolate barotropic velocities at step jit+0.5: 567 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)568 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)564 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 565 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 569 566 570 567 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) … … 628 625 ! Sum over sub-time-steps to compute advective velocities 629 626 za2 = wgtbtp2(jn) 630 zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u(:,:)631 zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v(:,:)627 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 628 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 632 629 ! 633 630 ! Set next sea level: … … 635 632 DO ji = fs_2, fs_jpim1 ! vector opt. 636 633 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 637 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1 2t(ji,jj)634 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 638 635 END DO 639 636 END DO … … 693 690 DO ji = 2, jpim1 694 691 zx1 = z1_2 * ssumask(ji ,jj) * r1_e12u(ji ,jj) & 695 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj) &696 & + e1 2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )692 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 693 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 697 694 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e12v(ji ,jj ) & 698 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj ) &699 & + e1 2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )695 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 696 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 700 697 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 701 698 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 … … 714 711 DO jj = 2, jpjm1 715 712 DO ji = fs_2, fs_jpim1 ! vector opt. 716 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)717 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)718 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) /e2v(ji,jj)719 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)713 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 714 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 715 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 716 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 720 717 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 721 718 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) … … 727 724 DO ji = fs_2, fs_jpim1 ! vector opt. 728 725 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 729 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)726 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 730 727 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 731 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)728 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 732 729 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 733 730 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 735 732 END DO 736 733 ! 737 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN !== energy and enstrophy conserving scheme ==!734 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 738 735 DO jj = 2, jpjm1 739 736 DO ji = fs_2, fs_jpim1 ! vector opt. 740 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &741 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &742 & + ftse(ji,jj ) * zwy(ji ,jj-1) &743 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )744 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &745 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &746 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &747 & + ftne(ji,jj ) * zwx(ji ,jj ) )737 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 738 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 739 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 740 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 741 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 742 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 743 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 744 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 748 745 END DO 749 746 END DO … … 755 752 DO jj = 2, jpjm1 756 753 DO ji = fs_2, fs_jpim1 ! vector opt. 757 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)758 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)754 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 755 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 759 756 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 760 757 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg … … 764 761 ! 765 762 ! Add bottom stresses: 766 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:)767 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:)763 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 764 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 768 765 ! 769 766 ! Add top stresses: … … 775 772 DO ji = fs_2, fs_jpim1 ! vector opt. 776 773 ! Add surface pressure gradient 777 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) /e1u(ji,jj)778 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) /e2v(ji,jj)774 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 775 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 779 776 zwx(ji,jj) = zu_spg 780 777 zwy(ji,jj) = zv_spg … … 786 783 DO jj = 2, jpjm1 787 784 DO ji = fs_2, fs_jpim1 ! vector opt. 788 ua_e(ji,jj) = ( zun_e(ji,jj) &785 ua_e(ji,jj) = ( un_e(ji,jj) & 789 786 & + rdtbt * ( zwx(ji,jj) & 790 787 & + zu_trd(ji,jj) & … … 792 789 & ) * ssumask(ji,jj) 793 790 794 va_e(ji,jj) = ( zvn_e(ji,jj) &791 va_e(ji,jj) = ( vn_e(ji,jj) & 795 792 & + rdtbt * ( zwy(ji,jj) & 796 793 & + zv_trd(ji,jj) & … … 807 804 zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 808 805 809 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) &806 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 810 807 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 811 808 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & … … 813 810 & ) * zhura 814 811 815 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &812 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 816 813 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 817 814 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & … … 837 834 #if defined key_bdy 838 835 ! open boundaries 839 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )836 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 840 837 #endif 841 838 #if defined key_agrif … … 845 842 ! ! ---- 846 843 ubb_e (:,:) = ub_e (:,:) 847 ub_e (:,:) = zun_e(:,:)848 zun_e(:,:) = ua_e (:,:)844 ub_e (:,:) = un_e (:,:) 845 un_e (:,:) = ua_e (:,:) 849 846 ! 850 847 vbb_e (:,:) = vb_e (:,:) 851 vb_e (:,:) = zvn_e(:,:)852 zvn_e(:,:) = va_e (:,:)848 vb_e (:,:) = vn_e (:,:) 849 vn_e (:,:) = va_e (:,:) 853 850 ! 854 851 sshbb_e(:,:) = sshb_e(:,:) … … 875 872 ! ----------------------------------------------------------------------------- 876 873 ! 877 ! At this stage ssha holds a time averaged value878 ! ! Sea Surface Height at u-,v- and f-points879 IF( lk_vvl ) THEN ! (required only in key_vvl case)880 DO jj = 1, jpjm1881 DO ji = 1, jpim1 ! NO Vector Opt.882 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e12u(ji,jj) &883 & * ( e12t(ji ,jj) * ssha(ji ,jj) &884 & + e12t(ji+1,jj) * ssha(ji+1,jj) )885 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e12v(ji,jj) &886 & * ( e12t(ji,jj ) * ssha(ji,jj ) &887 & + e12t(ji,jj+1) * ssha(ji,jj+1) )888 END DO889 END DO890 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions891 ENDIF892 !893 874 ! Set advection velocity correction: 875 zwx(:,:) = un_adv(:,:) 876 zwy(:,:) = vn_adv(:,:) 894 877 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 895 un_adv(:,:) = z u_sum(:,:)*hur(:,:)896 vn_adv(:,:) = z v_sum(:,:)*hvr(:,:)878 un_adv(:,:) = zwx(:,:)*hur(:,:) 879 vn_adv(:,:) = zwy(:,:)*hvr(:,:) 897 880 ELSE 898 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + z u_sum(:,:)) * hur(:,:)899 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + z v_sum(:,:)) * hvr(:,:)881 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:)) * hur(:,:) 882 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:)) * hvr(:,:) 900 883 END IF 901 884 902 885 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 903 ub2_b(:,:) = z u_sum(:,:)904 vb2_b(:,:) = z v_sum(:,:)886 ub2_b(:,:) = zwx(:,:) 887 vb2_b(:,:) = zwy(:,:) 905 888 ENDIF 906 889 ! … … 912 895 END DO 913 896 ELSE 897 ! At this stage, ssha has been corrected: compute new depths at velocity points 898 DO jj = 1, jpjm1 899 DO ji = 1, jpim1 ! NO Vector Opt. 900 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 901 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 902 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 903 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 904 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 905 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 906 END DO 907 END DO 908 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 909 ! 914 910 DO jk=1,jpkm1 915 911 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b … … 930 926 #if defined key_agrif 931 927 ! Save time integrated fluxes during child grid integration 932 ! (used to update coarse grid transports) 933 ! Useless with 2nd order momentum schemes 928 ! (used to update coarse grid transports at next time step) 934 929 ! 935 930 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN … … 951 946 ! 952 947 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 953 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd , zun_e, zvn_e)954 CALL wrk_dealloc( jpi, jpj, zwx, zwy, z u_sum, zv_sum, zssh_frc, zu_frc, zv_frc )948 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd ) 949 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 955 950 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 956 951 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) … … 1100 1095 ! 1101 1096 INTEGER :: ji ,jj 1102 INTEGER :: ios ! Local integer output status for namelist read1103 1097 REAL(wp) :: zxr2, zyr2, zcmax 1104 1098 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1105 1099 !! 1106 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, &1107 & nn_baro, rn_bt_cmax, nn_bt_flt1108 1100 !!---------------------------------------------------------------------- 1109 1101 ! 1110 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1111 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1112 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1113 1114 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1115 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1116 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1117 IF(lwm) WRITE ( numond, namsplit ) 1118 ! 1119 ! ! Max courant number for ext. grav. waves 1102 ! Max courant number for ext. grav. waves 1120 1103 ! 1121 1104 CALL wrk_alloc( jpi, jpj, zcu ) 1122 1105 ! 1123 IF (lk_vvl) THEN 1124 DO jj = 1, jpj 1125 DO ji =1, jpi 1126 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1127 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1128 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1129 END DO 1130 END DO 1131 ELSE 1132 DO jj = 1, jpj 1133 DO ji =1, jpi 1134 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1135 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1136 zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) ) 1137 END DO 1138 END DO 1139 ENDIF 1140 1141 zcmax = MAXVAL(zcu(:,:)) 1106 DO jj = 1, jpj 1107 DO ji =1, jpi 1108 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1109 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1110 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1111 END DO 1112 END DO 1113 ! 1114 zcmax = MAXVAL( zcu(:,:) ) 1142 1115 IF( lk_mpp ) CALL mpp_max( zcmax ) 1143 1116 1144 1117 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1145 IF (ln_bt_ nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1118 IF (ln_bt_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1146 1119 1147 rdtbt = rdt / FLOAT(nn_baro)1120 rdtbt = rdt / REAL( nn_baro , wp ) 1148 1121 zcmax = zcmax * rdtbt 1149 1122 ! Print results … … 1151 1124 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1152 1125 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1153 IF( ln_bt_ nn_auto ) THEN1154 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1126 IF( ln_bt_auto ) THEN 1127 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1155 1128 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1156 1129 ELSE 1157 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1130 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1158 1131 ENDIF 1159 1132 … … 1200 1173 END SUBROUTINE dyn_spg_ts_init 1201 1174 1202 #else1203 !!---------------------------------------------------------------------------1204 !! Default case : Empty module No split explicit free surface1205 !!---------------------------------------------------------------------------1206 CONTAINS1207 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1208 dyn_spg_ts_alloc = 01209 END FUNCTION dyn_spg_ts_alloc1210 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1211 INTEGER, INTENT(in) :: kt1212 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1213 END SUBROUTINE dyn_spg_ts1214 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1215 INTEGER , INTENT(in) :: kt ! ocean time-step1216 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1217 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1218 END SUBROUTINE ts_rst1219 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1220 INTEGER , INTENT(in) :: kt ! ocean time-step1221 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1222 END SUBROUTINE dyn_spg_ts_init1223 #endif1224 1225 1175 !!====================================================================== 1226 1176 END MODULE dynspg_ts 1227 1228 1229
Note: See TracChangeset
for help on using the changeset viewer.
