- Timestamp:
- 2015-01-07T19:03:53+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r4826_NOC_WAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r4770 r5014 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.? ! 2014-09 (H. Liu) Add Wetting/Drying components 13 14 !!--------------------------------------------------------------------- 14 15 #if defined key_dynspg_ts || defined key_esopa … … 17 18 !!---------------------------------------------------------------------- 18 19 !! dyn_spg_ts : compute surface pressure gradient trend using a time- 19 !! splitting scheme and add to the general trend 20 !! splitting scheme and add to the general trend 20 21 !!---------------------------------------------------------------------- 21 22 USE oce ! ocean dynamics and tracers … … 26 27 USE dynvor ! vorticity term 27 28 USE bdy_par ! for lk_bdy 28 USE bdytides ! open boundary condition data 29 USE bdytides ! open boundary condition data 29 30 USE bdydyn2d ! open boundary conditions on barotropic variables 30 31 USE sbctide ! tides … … 38 39 USE zdf_oce ! Bottom friction coefts 39 40 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 41 USE timing ! Timing 41 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE wadlmt ! wetting/drying flux limter 42 44 USE dynadv, ONLY: ln_dynadv_vec 43 45 #if defined key_agrif 44 46 USE agrif_opa_interp ! agrif 45 47 #endif 46 #if defined key_asminc 48 #if defined key_asminc 47 49 USE asminc ! Assimilation increment 48 50 #endif … … 51 53 PRIVATE 52 54 53 PUBLIC dyn_spg_ts ! routine called in dynspg.F90 55 PUBLIC dyn_spg_ts ! routine called in dynspg.F90 54 56 PUBLIC dyn_spg_ts_alloc ! " " " " 55 57 PUBLIC dyn_spg_ts_init ! " " " " … … 97 99 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 98 100 99 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 100 102 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 101 103 … … 110 112 !!---------------------------------------------------------------------- 111 113 !! 112 !! ** Purpose : 114 !! ** Purpose : 113 115 !! -Compute the now trend due to the explicit time stepping 114 !! of the quasi-linear barotropic system. 116 !! of the quasi-linear barotropic system. 115 117 !! 116 !! ** Method : 118 !! ** Method : 117 119 !! Barotropic variables are advanced from internal time steps 118 120 !! "n" to "n+1" if ln_bt_fw=T 119 !! or from 121 !! or from 120 122 !! "n-1" to "n+1" if ln_bt_fw=F 121 123 !! thanks to a generalized forward-backward time stepping (see ref. below). … … 126 128 !! -Compute barotropic advective velocities at step "n" : un_adv, vn_adv 127 129 !! These are used to advect tracers and are compliant with discrete 128 !! continuity equation taken at the baroclinic time steps. This 130 !! continuity equation taken at the baroclinic time steps. This 129 131 !! ensures tracers conservation. 130 132 !! -Update 3d trend (ua, va) with barotropic component. 131 133 !! 132 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 133 !! The regional oceanic modeling system (ROMS): 134 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 135 !! The regional oceanic modeling system (ROMS): 134 136 !! a split-explicit, free-surface, 135 !! topography-following-coordinate oceanic model. 136 !! Ocean Modelling, 9, 347-404. 137 !! topography-following-coordinate oceanic model. 138 !! Ocean Modelling, 9, 347-404. 137 139 !!--------------------------------------------------------------------- 138 140 ! 139 141 INTEGER, INTENT(in) :: kt ! ocean time-step index 140 142 ! 141 LOGICAL :: ll_fw_start ! if true, forward integration 142 LOGICAL :: ll_init ! if true, special startup of 2d equations 143 INTEGER :: ji, jj, jk, jn ! dummy loop indices 144 INTEGER :: ikbu, ikbv, noffset ! local integers 145 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 143 LOGICAL :: ll_fw_start ! if true, forward integration 144 LOGICAL :: ll_init ! if true, special startup of 2d equations 145 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 146 INTEGER :: ji, jj, jk, jn ! dummy loop indices 147 INTEGER :: ikbu, ikbv, noffset ! local integers 148 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 146 149 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 147 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 150 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 148 151 REAL(wp) :: zu_spg, zv_spg ! - - 149 REAL(wp) :: zhura, zhvra ! - - 150 REAL(wp) :: za0, za1, za2, za3 ! - - 152 REAL(wp) :: zhura, zhvra ! - - 153 REAL(wp) :: za0, za1, za2, za3 ! - - 154 155 REAL(wp) :: ztmp ! temporary vaialbe 151 156 ! 152 157 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e … … 156 161 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 157 162 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 163 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 164 !!---------------------------------------------------------------------- 159 165 ! … … 167 173 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 168 174 CALL wrk_alloc( jpi, jpj, zhf ) 175 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 169 176 ! 170 177 ! !* Local constant initialization 171 z1_12 = 1._wp / 12._wp 172 z1_8 = 0.125_wp 178 z1_12 = 1._wp / 12._wp 179 z1_8 = 0.125_wp 173 180 z1_4 = 0.25_wp 174 z1_2 = 0.5_wp 181 z1_2 = 0.5_wp 175 182 zraur = 1._wp / rau0 176 183 ! 177 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step184 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 178 185 z2dt_bf = rdt 179 186 ELSE 180 187 z2dt_bf = 2.0_wp * rdt 181 188 ENDIF 182 z1_2dt_b = 1.0_wp / z2dt_bf 183 ! 184 ll_init = ln_bt_av ! if no time averaging, then no specific restart189 z1_2dt_b = 1.0_wp / z2dt_bf 190 ! 191 ll_init = ln_bt_av ! if no time averaging, then no specific restart 185 192 ll_fw_start = .FALSE. 186 193 ! 187 194 ! time offset in steps for bdy data update 188 195 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 189 196 ! 190 IF( kt == nit000 ) THEN 197 IF( kt == nit000 ) THEN !* initialisation 191 198 ! 192 199 IF(lwp) WRITE(numout,*) … … 245 252 IF ( .not. ln_sco ) THEN 246 253 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 247 ! Set it to zero for the time being 254 ! Set it to zero for the time being 248 255 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 249 256 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth … … 264 271 END DO 265 272 CALL lbc_lnk( zhf, 'F', 1._wp ) 266 ! JC: TBC. hf should be greater than 0 273 ! JC: TBC. hf should be greater than 0 267 274 DO jj = 1, jpj 268 275 DO ji = 1, jpi … … 274 281 ENDIF 275 282 ! 276 ! If forward start at previous time step, and centered integration, 283 ! If forward start at previous time step, and centered integration, 277 284 ! then update averaging weights: 278 285 IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN … … 280 287 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 281 288 ENDIF 282 289 283 290 ! ----------------------------------------------------------------------------- 284 291 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 285 292 ! ----------------------------------------------------------------------------- 286 ! 293 ! 287 294 ! 288 295 ! !* e3*d/dt(Ua) (Vertically integrated) … … 293 300 DO jk = 1, jpkm1 294 301 zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 295 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 302 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 296 303 END DO 297 304 ! … … 311 318 ! !* barotropic Coriolis trends (vorticity scheme dependent) 312 319 ! ! -------------------------------------------------------- 313 zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes 320 zwx(:,:) = un_b(:,:) * hu(:,:) * e2u(:,:) ! now fluxes 314 321 zwy(:,:) = vn_b(:,:) * hv(:,:) * e1v(:,:) 315 322 ! … … 353 360 END DO 354 361 ! 355 ENDIF 362 ENDIF 356 363 ! 357 364 ! !* Right-Hand-Side of the barotropic momentum equation 358 365 ! ! ---------------------------------------------------- 359 366 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 360 DO jj = 2, jpjm1 361 DO ji = fs_2, fs_jpim1 ! vector opt. 362 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 363 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 364 END DO 365 END DO 367 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 368 DO jj = 2, jpjm1 369 DO ji = 2, jpim1 370 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 371 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 372 & rn_wdmin1 + rn_wdmin2 373 IF(ll_tmp1) THEN 374 zcpx(ji,jj) = 1.0_wp 375 ELSE IF(ll_tmp2) THEN 376 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen ! here 377 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) /& 378 & (sshn(ji+1,jj) - sshn(ji,jj))) 379 ELSE 380 zcpx(ji,jj) = 0._wp 381 END IF 382 383 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 384 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 385 & rn_wdmin1 + rn_wdmin2 386 IF(ll_tmp1) THEN 387 zcpy(ji,jj) = 1.0_wp 388 ELSE IF(ll_tmp2) THEN 389 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen ! here 390 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) + bathy(ji,jj)) /& 391 & (sshn(ji,jj+1) - sshn(ji,jj))) 392 ELSE 393 zcpy(ji,jj) = 0._wp 394 END IF 395 END DO 396 END DO 397 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 398 ENDIF 399 400 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 401 DO jj = 2, jpjm1 402 DO ji = 2, jpim1 403 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / & 404 & e1u(ji,jj) * zcpx(ji,jj) 405 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / & 406 & e2v(ji,jj) * zcpy(ji,jj) 407 END DO 408 END DO 409 ELSE 410 DO jj = 2, jpjm1 411 DO ji = fs_2, fs_jpim1 ! vector opt. 412 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 413 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 414 END DO 415 END DO 416 END IF 417 366 418 ENDIF 367 419 … … 371 423 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 372 424 END DO 373 END DO 374 ! 375 ! ! Add bottom stress contribution from baroclinic velocities:425 END DO 426 ! 427 ! ! Add bottom stress contribution from baroclinic velocities: 376 428 IF (ln_bt_fw) THEN 377 DO jj = 2, jpjm1 429 DO jj = 2, jpjm1 378 430 DO ji = fs_2, fs_jpim1 ! vector opt. 379 ikbu = mbku(ji,jj) 380 ikbv = mbkv(ji,jj) 431 ikbu = mbku(ji,jj) 432 ikbv = mbkv(ji,jj) 381 433 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 382 434 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) … … 386 438 DO jj = 2, jpjm1 387 439 DO ji = fs_2, fs_jpim1 ! vector opt. 388 ikbu = mbku(ji,jj) 389 ikbv = mbkv(ji,jj) 440 ikbu = mbku(ji,jj) 441 ikbv = mbkv(ji,jj) 390 442 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 391 443 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) … … 397 449 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 398 450 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 399 ! 451 ! 400 452 IF (ln_bt_fw) THEN ! Add wind forcing 401 453 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) … … 404 456 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 405 457 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 406 ENDIF 458 ENDIF 407 459 ! 408 460 IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing 409 461 IF (ln_bt_fw) THEN 410 DO jj = 2, jpjm1 462 DO jj = 2, jpjm1 411 463 DO ji = fs_2, fs_jpim1 ! vector opt. 412 464 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj) … … 417 469 END DO 418 470 ELSE 419 DO jj = 2, jpjm1 471 DO jj = 2, jpjm1 420 472 DO ji = fs_2, fs_jpim1 ! vector opt. 421 473 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & … … 427 479 END DO 428 480 END DO 429 ENDIF 481 ENDIF 430 482 ENDIF 431 483 ! !* Right-Hand-Side of the barotropic ssh equation … … 450 502 ! 451 503 ! ----------------------------------------------------------------------- 452 ! Phase 2 : Integration of the barotropic equations 504 ! Phase 2 : Integration of the barotropic equations 453 505 ! ----------------------------------------------------------------------- 454 506 ! 455 507 ! ! ==================== ! 456 508 ! ! Initialisations ! 457 ! ! ==================== ! 458 ! Initialize barotropic variables: 509 ! ! ==================== ! 510 ! Initialize barotropic variables: 459 511 IF( ll_init )THEN 460 512 sshbb_e(:,:) = 0._wp … … 465 517 vb_e (:,:) = 0._wp 466 518 ENDIF 467 ! 468 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 469 sshn_e(:,:) = sshn (:,:) 470 zun_e (:,:) = un_b (:,:) 519 520 IF(ln_wd) THEN !preserve the positivity of water depth 521 !ssh[b,n,a] should have already been processed for this 522 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 523 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 524 ENDIF 525 ! 526 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 527 sshn_e(:,:) = sshn (:,:) 528 zun_e (:,:) = un_b (:,:) 471 529 zvn_e (:,:) = vn_b (:,:) 472 530 ! 473 hu_e (:,:) = hu (:,:) 474 hv_e (:,:) = hv (:,:) 475 hur_e (:,:) = hur (:,:) 531 hu_e (:,:) = hu (:,:) 532 hv_e (:,:) = hv (:,:) 533 hur_e (:,:) = hur (:,:) 476 534 hvr_e (:,:) = hvr (:,:) 477 535 ELSE ! CENTRED integration: start from BEFORE fields 478 536 sshn_e(:,:) = sshb (:,:) 479 zun_e (:,:) = ub_b (:,:) 537 zun_e (:,:) = ub_b (:,:) 480 538 zvn_e (:,:) = vb_b (:,:) 481 539 ! 482 hu_e (:,:) = hu_b (:,:) 483 hv_e (:,:) = hv_b (:,:) 484 hur_e (:,:) = hur_b(:,:) 540 hu_e (:,:) = hu_b (:,:) 541 hv_e (:,:) = hv_b (:,:) 542 hur_e (:,:) = hur_b(:,:) 485 543 hvr_e (:,:) = hvr_b(:,:) 486 544 ENDIF … … 489 547 ! 490 548 ! Initialize sums: 491 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 549 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 492 550 va_b (:,:) = 0._wp 493 551 ssha (:,:) = 0._wp ! Sum for after averaged sea level … … 506 564 ! 507 565 ! Set extrapolation coefficients for predictor step: 508 IF ((jn<3).AND.ll_init) THEN ! Forward 509 za1 = 1._wp 510 za2 = 0._wp 511 za3 = 0._wp 512 ELSE ! AB3-AM4 Coefficients: bet=0.281105 566 IF ((jn<3).AND.ll_init) THEN ! Forward 567 za1 = 1._wp 568 za2 = 0._wp 569 za3 = 0._wp 570 ELSE ! AB3-AM4 Coefficients: bet=0.281105 513 571 za1 = 1.781105_wp ! za1 = 3/2 + bet 514 572 za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet) … … 524 582 ! Extrapolate Sea Level at step jit+0.5: 525 583 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 584 585 526 586 ! 527 587 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points … … 535 595 END DO 536 596 END DO 597 537 598 CALL lbc_lnk( zwx, 'U', 1._wp ) ; CALL lbc_lnk( zwy, 'V', 1._wp ) 538 599 ! 539 600 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 540 601 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 602 IF(ln_wd) THEN 603 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 604 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 605 END IF 541 606 ELSE 542 607 zhup2_e (:,:) = hu(:,:) … … 552 617 ! 553 618 #if defined key_agrif 554 ! Set fluxes during predictor step to ensure 619 ! Set fluxes during predictor step to ensure 555 620 ! volume conservation 556 621 IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN … … 577 642 ENDIF 578 643 #endif 644 IF(ln_wd) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 579 645 ! 580 646 ! Sum over sub-time-steps to compute advective velocities … … 584 650 ! 585 651 ! Set next sea level: 586 DO jj = 2, jpjm1 652 DO jj = 2, jpjm1 587 653 DO ji = fs_2, fs_jpim1 ! vector opt. 588 654 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & … … 590 656 END DO 591 657 END DO 658 !IF(ln_wd) THEN 659 592 660 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 661 IF(ln_wd) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 593 662 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 594 663 … … 600 669 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 601 670 #endif 602 ! 671 ! 603 672 ! Sea Surface Height at u-,v-points (vvl case only) 604 IF ( lk_vvl ) THEN 673 IF ( lk_vvl ) THEN 605 674 DO jj = 2, jpjm1 606 675 DO ji = 2, jpim1 ! NO Vector Opt. … … 613 682 END DO 614 683 END DO 684 615 685 CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 616 ENDIF 617 ! 686 ENDIF 687 ! 618 688 ! Half-step back interpolation of SSH for surface pressure computation: 619 689 !---------------------------------------------------------------------- 620 690 IF ((jn==1).AND.ll_init) THEN 621 691 za0=1._wp ! Forward-backward 622 za1=0._wp 692 za1=0._wp 623 693 za2=0._wp 624 694 za3=0._wp … … 627 697 za1=-0.1666666666666_wp ! za1 = gam 628 698 za2= 0.0833333333333_wp ! za2 = eps 629 za3= 0._wp 630 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 631 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 632 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 699 za3= 0._wp 700 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 701 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 702 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 633 703 za2=0.088_wp ! za2 = gam 634 704 za3=0.013_wp ! za3 = eps … … 640 710 ! 641 711 ! Compute associated depths at U and V points: 642 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN 643 ! 644 DO jj = 2, jpjm1 712 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN 713 ! 714 DO jj = 2, jpjm1 645 715 DO ji = 2, jpim1 646 716 zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) & … … 650 720 & * ( e12t(ji ,jj ) * zsshp2_e(ji ,jj ) & 651 721 & + e12t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 652 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 722 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 653 723 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 654 724 END DO 655 725 END DO 726 727 IF(ln_wd) THEN 728 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 729 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 730 END IF 731 !shall we call lbc_lnk for zhu(v)st_e() here? 732 656 733 ENDIF 657 734 ! … … 692 769 zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 693 770 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 694 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 771 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 695 772 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 696 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 773 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 697 774 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 698 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 775 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 699 776 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 700 777 END DO 701 778 END DO 702 ! 779 ! 703 780 ENDIF 704 781 ! … … 720 797 ! 721 798 ! Surface pressure trend: 722 DO jj = 2, jpjm1 723 DO ji = fs_2, fs_jpim1 ! vector opt. 724 ! Add surface pressure gradient 725 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 726 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 727 zwx(ji,jj) = zu_spg 728 zwy(ji,jj) = zv_spg 729 END DO 730 END DO 799 IF(ln_wd) THEN ! Calculating and applying W/D gravity filters 800 DO jj = 2, jpjm1 801 DO ji = 2, jpim1 802 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) 803 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) + & 804 & rn_wdmin1 + rn_wdmin2 805 IF(ll_tmp1) THEN 806 zcpx(ji,jj) = 1.0_wp 807 ELSE IF(ll_tmp2) THEN 808 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen ! here 809 zcpx(ji,jj) = ABS((zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 810 & (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj))) 811 ELSE 812 zcpx(ji,jj) = 0._wp 813 END IF 814 815 ll_tmp1 = MIN(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) 816 ll_tmp2 = MAX(zsshp2_e(ji,jj), zsshp2_e(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) + & 817 & rn_wdmin1 + rn_wdmin2 818 IF(ll_tmp1) THEN 819 zcpy(ji,jj) = 1.0_wp 820 ELSE IF(ll_tmp2) THEN 821 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen ! here 822 zcpy(ji,jj) = ABS((zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) /& 823 & (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj))) 824 ELSE 825 zcpy(ji,jj) = 0._wp 826 END IF 827 END DO 828 END DO 829 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 830 ENDIF 831 832 IF(ln_wd) THEN 833 DO jj = 2, jpjm1 834 DO ji = 2, jpim1 835 ! Add surface pressure gradient 836 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 837 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 838 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 839 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 840 END DO 841 END DO 842 ELSE 843 DO jj = 2, jpjm1 844 DO ji = fs_2, fs_jpim1 ! vector opt. 845 ! Add surface pressure gradient 846 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 847 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 848 zwx(ji,jj) = zu_spg 849 zwy(ji,jj) = zv_spg 850 END DO 851 END DO 852 END IF 853 731 854 ! 732 855 ! Set next velocities: 733 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN 856 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form 734 857 DO jj = 2, jpjm1 735 858 DO ji = fs_2, fs_jpim1 ! vector opt. 736 ua_e(ji,jj) = ( zun_e(ji,jj) & 859 ua_e(ji,jj) = ( zun_e(ji,jj) & 737 860 & + rdtbt * ( zwx(ji,jj) & 738 861 & + zu_trd(ji,jj) & 739 & + zu_frc(ji,jj) ) & 862 & + zu_frc(ji,jj) ) & 740 863 & ) * umask(ji,jj,1) 741 864 … … 748 871 END DO 749 872 750 ELSE 873 ELSE ! Flux form 751 874 DO jj = 2, jpjm1 752 875 DO ji = fs_2, fs_jpim1 ! vector opt. 753 876 754 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 755 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 756 757 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 758 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 877 IF(ln_wd) THEN 878 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 879 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 880 ELSE 881 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 882 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 883 END IF 884 885 zhura = umask(ji,jj,1)/(zhura + 1._wp - umask(ji,jj,1)) 886 zhvra = vmask(ji,jj,1)/(zhvra + 1._wp - vmask(ji,jj,1)) 887 888 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 889 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 759 890 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 760 891 & + hu(ji,jj) * zu_frc(ji,jj) ) & … … 771 902 ! 772 903 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 773 ! ! ---------------------------------------------- 774 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 775 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 904 ! ! ---------------------------------------------- 905 IF(ln_wd) THEN 906 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 907 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 908 ELSE 909 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 910 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 911 END IF 912 776 913 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 777 914 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) … … 781 918 ! ! ----------------------- 782 919 ! 783 CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries 920 CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries 784 921 CALL lbc_lnk( va_e , 'V', -1._wp ) 785 922 786 #if defined key_bdy 923 #if defined key_bdy 787 924 ! open boundaries 788 925 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e ) 789 926 #endif 790 #if defined key_agrif 927 #if defined key_agrif 791 928 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( jn ) ! Agrif 792 929 #endif … … 807 944 ! !* Sum over whole bt loop 808 945 ! ! ---------------------- 809 za1 = wgtbtp1(jn) 946 za1 = wgtbtp1(jn) 810 947 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities 811 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 812 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 813 ELSE 948 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 949 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 950 ELSE ! Sum transports 814 951 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 815 952 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 816 953 ENDIF 817 ! 954 ! ! Sum sea level 818 955 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 819 956 ! ! ==================== ! … … 841 978 ! 842 979 ! Set advection velocity correction: 843 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 980 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 844 981 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 845 982 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) … … 894 1031 ! 895 1032 ! 896 #endif 1033 #endif 897 1034 ! 898 1035 ! !* write time-spliting arrays in the restart … … 905 1042 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 906 1043 CALL wrk_dealloc( jpi, jpj, zhf ) 1044 IF(ln_wd) CALL wrk_alloc( jpi, jpj, zcpx, zcpy ) 907 1045 ! 908 1046 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 918 1056 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 919 1057 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 920 INTEGER, INTENT(inout) :: jpit ! cycle length 1058 INTEGER, INTENT(inout) :: jpit ! cycle length 921 1059 REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights 922 1060 zwgt2 ! Secondary weights 923 1061 924 1062 INTEGER :: jic, jn, ji ! temporary integers 925 1063 REAL(wp) :: za1, za2 … … 930 1068 931 1069 ! Set time index when averaged value is requested 932 IF (ll_fw) THEN 1070 IF (ll_fw) THEN 933 1071 jic = nn_baro 934 1072 ELSE … … 938 1076 ! Set primary weights: 939 1077 IF (ll_av) THEN 940 ! Define simple boxcar window for primary weights 941 ! (width = nn_baro, centered around jic) 1078 ! Define simple boxcar window for primary weights 1079 ! (width = nn_baro, centered around jic) 942 1080 SELECT CASE ( nn_bt_flt ) 943 1081 CASE( 0 ) ! No averaging … … 947 1085 CASE( 1 ) ! Boxcar, width = nn_baro 948 1086 DO jn = 1, 3*nn_baro 949 za1 = ABS(float(jn-jic))/float(nn_baro) 1087 za1 = ABS(float(jn-jic))/float(nn_baro) 950 1088 IF (za1 < 0.5_wp) THEN 951 1089 zwgt1(jn) = 1._wp … … 956 1094 CASE( 2 ) ! Boxcar, width = 2 * nn_baro 957 1095 DO jn = 1, 3*nn_baro 958 za1 = ABS(float(jn-jic))/float(nn_baro) 1096 za1 = ABS(float(jn-jic))/float(nn_baro) 959 1097 IF (za1 < 1._wp) THEN 960 1098 zwgt1(jn) = 1._wp … … 969 1107 jpit = jic 970 1108 ENDIF 971 1109 972 1110 ! Set secondary weights 973 1111 DO jn = 1, jpit … … 999 1137 ! 1000 1138 IF( TRIM(cdrw) == 'READ' ) THEN 1001 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) ) 1002 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) ) 1139 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) ) 1140 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) ) 1003 1141 IF( .NOT.ln_bt_av ) THEN 1004 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) 1005 CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) ) 1142 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) 1143 CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) ) 1006 1144 CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) ) 1007 CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) ) 1008 CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) ) 1145 CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) ) 1146 CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) ) 1009 1147 CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) ) 1010 1148 ENDIF … … 1012 1150 ! Read time integrated fluxes 1013 1151 IF ( .NOT.Agrif_Root() ) THEN 1014 CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) ) 1152 CALL iom_get( numror, jpdom_autoglo, 'ub2_i_b' , ub2_i_b(:,:) ) 1015 1153 CALL iom_get( numror, jpdom_autoglo, 'vb2_i_b' , vb2_i_b(:,:) ) 1016 1154 ENDIF … … 1022 1160 ! 1023 1161 IF (.NOT.ln_bt_av) THEN 1024 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 1162 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 1025 1163 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 1026 1164 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) … … 1070 1208 CALL wrk_alloc( jpi, jpj, zcu ) 1071 1209 ! 1072 IF (lk_vvl) THEN 1210 IF (lk_vvl) THEN 1073 1211 DO jj = 1, jpj 1074 1212 DO ji =1, jpi … … 1093 1231 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1094 1232 IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1095 1233 1096 1234 rdtbt = rdt / FLOAT(nn_baro) 1097 1235 zcmax = zcmax * rdtbt 1098 1236 ! Print results 1099 1237 IF(lwp) WRITE(numout,*) 1100 1238 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' … … 1129 1267 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1130 1268 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1131 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1269 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1132 1270 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1133 1271 END SELECT … … 1142 1280 ENDIF 1143 1281 IF ( zcmax>0.9_wp ) THEN 1144 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1282 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1145 1283 ENDIF 1146 1284 ! … … 1165 1303 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 1166 1304 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw 1167 END SUBROUTINE ts_rst 1305 END SUBROUTINE ts_rst 1168 1306 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine 1169 1307 INTEGER , INTENT(in) :: kt ! ocean time-step … … 1171 1309 END SUBROUTINE dyn_spg_ts_init 1172 1310 #endif 1173 1311 1174 1312 !!====================================================================== 1175 1313 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
