- Timestamp:
- 2015-12-07T16:11:45+01:00 (9 years ago)
- Location:
- branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN
- Files:
-
- 6 edited
-
divhor.F90 (modified) (2 diffs)
-
dynhpg.F90 (modified) (8 diffs)
-
dynnxt.F90 (modified) (2 diffs)
-
dynspg.F90 (modified) (1 diff)
-
dynspg_ts.F90 (modified) (14 diffs)
-
dynzad.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/divhor.F90
r6006 r6012 20 20 USE oce ! ocean dynamics and tracers 21 21 USE dom_oce ! ocean space and time domain 22 USE sbc_oce, ONLY : ln_rnf, nn_isf ! surface boundary condition: ocean22 USE sbc_oce, ONLY : ln_rnf, ln_isf ! surface boundary condition: ocean 23 23 USE sbcrnf ! river runoff 24 24 USE sbcisf ! ice shelf … … 91 91 END DO 92 92 ! 93 IF( ln_rnf 93 IF( ln_rnf ) CALL sbc_rnf_div( hdivn ) !== runoffs ==! (update hdivn field) 94 94 ! 95 IF( ln_ divisf .AND. nn_isf > 0) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field)95 IF( ln_isf ) CALL sbc_isf_div( hdivn ) !== ice shelf ==! (update hdivn field) 96 96 ! 97 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn )!== ice sheet ==! (update hdivn field)97 IF( ln_iscpl .AND. ln_hsb ) CALL iscpl_div( hdivn ) !== ice sheet ==! (update hdivn field) 98 98 ! 99 CALL lbc_lnk( hdivn, 'T', 1. ) 99 CALL lbc_lnk( hdivn, 'T', 1. ) !== lateral boundary cond. ==! (no sign change) 100 100 ! 101 101 IF( nn_timing == 1 ) CALL timing_stop('div_hor') -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5930 r6012 45 45 USE wrk_nemo ! Memory Allocation 46 46 USE timing ! Timing 47 USE iom 47 48 48 49 IMPLICIT NONE … … 130 131 INTEGER :: ioptio = 0 ! temporary integer 131 132 INTEGER :: ios ! Local integer output status for namelist read 133 !! 134 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 135 REAL(wp) :: znad 136 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 137 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 138 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 132 139 !! 133 140 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & … … 190 197 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 191 198 ! 192 ! initialisation of ice load 193 riceload(:,:)=0.0 199 ! initialisation of ice shelf load 200 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 201 IF ( ln_isfcav ) THEN 202 CALL wrk_alloc( jpi,jpj, 2, ztstop) 203 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 204 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 205 ! 206 IF(lwp) WRITE(numout,*) 207 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 208 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 209 210 ! To use density and not density anomaly 211 znad=1._wp 212 213 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 214 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 215 216 ! compute density of the water displaced by the ice shelf 217 DO jk = 1, jpk 218 CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 219 END DO 220 221 ! compute rhd at the ice/oce interface (ice shelf side) 222 CALL eos(ztstop,risfdep,zrhdtop_isf) 223 224 ! Surface value + ice shelf gradient 225 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 226 ! divided by 2 later 227 ziceload = 0._wp 228 DO jj = 1, jpj 229 DO ji = 1, jpi 230 ikt=mikt(ji,jj) 231 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 232 DO jk=2,ikt-1 233 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 234 & * (1._wp - tmask(ji,jj,jk)) 235 END DO 236 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 237 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 238 END DO 239 END DO 240 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 241 242 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 243 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 244 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 245 END IF 194 246 ! 195 247 END SUBROUTINE dyn_hpg_init … … 447 499 SUBROUTINE hpg_isf( kt ) 448 500 !!--------------------------------------------------------------------- 449 !! *** ROUTINE hpg_ sco***501 !! *** ROUTINE hpg_isf *** 450 502 !! 451 503 !! ** Method : s-coordinate case. Jacobian scheme. … … 466 518 INTEGER, INTENT(in) :: kt ! ocean time-step index 467 519 !! 468 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices469 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars470 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd520 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 521 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 522 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 471 523 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 472 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 473 !!---------------------------------------------------------------------- 474 ! 475 CALL wrk_alloc( jpi,jpj, 2, ztstop) 476 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 477 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 478 ! 479 IF( kt == nit000 ) THEN 480 IF(lwp) WRITE(numout,*) 481 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 482 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 483 ENDIF 484 524 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 525 !!---------------------------------------------------------------------- 526 ! 527 CALL wrk_alloc( jpi,jpj, 2, ztstop) 528 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 529 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 530 ! 485 531 ! Local constant initialization 486 532 zcoef0 = - grav * 0.5_wp 533 487 534 ! To use density and not density anomaly 488 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume489 ! ELSE ; znad = 0._wp ! Fixed volume490 ! ENDIF491 535 znad=1._wp 536 492 537 ! iniitialised to 0. zhpi zhpi 493 538 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 494 539 495 ! Partial steps: top & bottom before horizontal gradient496 !jc CALL zps_hde_isf( kt, jpts, tsn, gtsu, gtsv, gtui, gtvi, &497 !jc & rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv , &498 !jc & grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi )499 500 !==================================================================================501 !=====Compute iceload and contribution of the half first wet layer =================502 !===================================================================================503 504 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude)505 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp506 507 ! compute density of the water displaced by the ice shelf508 zrhd = rhd ! save rhd509 DO jk = 1, jpk510 zdept(:,:)=gdept_1d(jk)511 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk))512 END DO513 WHERE ( tmask(:,:,:) == 1._wp)514 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd515 END WHERE516 517 ! compute rhd at the ice/oce interface (ice shelf side)518 CALL eos(ztstop,risfdep,zrhdtop_isf)519 520 540 ! compute rhd at the ice/oce interface (ocean side) 541 ! usefull to reduce residual current in the test case ISOMIP with no melting 521 542 DO ji=1,jpi 522 543 DO jj=1,jpj … … 527 548 END DO 528 549 CALL eos(ztstop,risfdep,zrhdtop_oce) 529 ! 530 ! Surface value + ice shelf gradient 531 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 532 ziceload = 0._wp 533 DO jj = 1, jpj 534 DO ji = 1, jpi ! vector opt. 535 ikt=mikt(ji,jj) 536 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 537 DO jk=2,ikt-1 538 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 539 & * (1._wp - tmask(ji,jj,jk)) 540 END DO 541 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 542 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 543 END DO 544 END DO 545 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 546 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 550 551 !================================================================================== 552 !===== Compute surface value ===================================================== 553 !================================================================================== 547 554 DO jj = 2, jpjm1 548 555 DO ji = fs_2, fs_jpim1 ! vector opt. … … 550 557 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 551 558 ! we assume ISF is in isostatic equilibrium 552 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj 553 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj) ) &554 & - 0.5_wp * fse3w(ji ,jj ,ikt )&555 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&556 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)))557 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji 558 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) &559 & - 0.5_wp * fse3w(ji ,jj ,ikt )&560 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) )&561 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ))559 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj,iktp1i) & 560 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 561 & - 0.5_wp * fse3w(ji,jj,ikt) & 562 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 563 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 564 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji,jj+1,iktp1j) & 565 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 566 & - 0.5_wp * fse3w(ji,jj,ikt) & 567 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 568 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 562 569 ! s-coordinate pressure gradient correction (=0 if z coordinate) 563 570 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 570 577 END DO 571 578 END DO 572 !==================================================================================573 !===== Compute partial cell contribution for the top cell =========================574 !==================================================================================575 DO jj = 2, jpjm1576 DO ji = fs_2, fs_jpim1 ! vector opt.577 iku = miku(ji,jj) ;578 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp579 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku))580 ! u direction581 IF ( iku .GT. 1 ) THEN582 ! case iku583 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &584 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &585 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )586 ! corrective term ( = 0 if z coordinate )587 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)588 ! zhpi will be added in interior loop589 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap590 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure591 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)592 593 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)594 zhpiint = zcoef0 / e1u(ji,jj) &595 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &596 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &597 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &598 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )599 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint600 END IF601 602 ! v direction603 ikv = mikv(ji,jj)604 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv))605 IF ( ikv .GT. 1 ) THEN606 ! case ikv607 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &608 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &609 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )610 ! corrective term ( = 0 if z coordinate )611 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)612 ! zhpi will be added in interior loop613 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap614 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure615 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)616 617 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)618 zhpjint = zcoef0 / e2v(ji,jj) &619 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &620 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &621 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &622 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )623 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint624 END IF625 END DO626 END DO627 579 628 580 !================================================================================== 629 581 !===== Compute interior value ===================================================== 630 582 !================================================================================== 631 632 DO jj = 2, jpjm1 633 DO ji = fs_2, fs_jpim1 ! vector opt. 634 iku=miku(ji,jj); ikv=mikv(ji,jj) 635 DO jk = 2, jpkm1 583 ! interior value (2=<jk=<jpkm1) 584 DO jk = 2, jpkm1 585 DO jj = 2, jpjm1 586 DO ji = fs_2, fs_jpim1 ! vector opt. 636 587 ! hydrostatic pressure gradient along s-surfaces 637 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 638 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 639 & + zcoef0 / e1u(ji,jj) & 640 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 641 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 642 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 643 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 588 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 589 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 590 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 591 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 592 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 593 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 644 594 ! s-coordinate pressure gradient correction 645 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 646 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 647 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 648 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 649 650 ! hydrostatic pressure gradient along s-surfaces 651 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 652 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 653 & + zcoef0 / e2v(ji,jj) & 654 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 655 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 656 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 657 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 658 ! s-coordinate pressure gradient correction 659 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 660 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 661 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 595 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 596 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 597 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 598 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 662 599 ! add to the general momentum trend 663 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 600 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 601 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 664 602 END DO 665 603 END DO 666 604 END DO 667 668 !================================================================================== 669 !===== Compute bottom cell contribution (partial cell) ============================ 670 !================================================================================== 671 672 DO jj = 2, jpjm1 673 DO ji = 2, jpim1 674 iku = mbku(ji,jj) 675 ikv = mbkv(ji,jj) 676 677 IF (iku .GT. 1) THEN 678 ! remove old value (interior case) 679 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 680 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 681 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 682 ! put new value 683 ! -zpshpi to avoid double contribution of the partial step in the top layer 684 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 685 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 686 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 687 END IF 688 ! v direction 689 IF (ikv .GT. 1) THEN 690 ! remove old value (interior case) 691 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 692 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 693 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 694 ! put new value 695 ! -zpshpj to avoid double contribution of the partial step in the top layer 696 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 697 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 698 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 699 END IF 700 END DO 701 END DO 702 703 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 704 rhd = zrhd 705 ! 706 CALL wrk_dealloc( jpi,jpj,2, ztstop) 707 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 708 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 605 ! 606 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 607 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 608 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 709 609 ! 710 610 END SUBROUTINE hpg_isf -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynnxt.F90
r6006 r6012 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices 97 INTEGER :: ikt ! local integers 97 98 REAL(wp) :: zue3a, zue3n, zue3b, zuf ! local scalars 98 99 REAL(wp) :: zve3a, zve3n, zve3b, zvf, z1_2dt ! - - … … 220 221 ! Add volume filter correction: compatibility with tracer advection scheme 221 222 ! => time filter + conservation correction (only at the first level) 222 IF ( nn_isf == 0) THEN ! if no ice shelf melting223 IF ( .NOT. ln_isf ) THEN ! if no ice shelf melting 223 224 fse3t_b(:,:,1) = fse3t_b(:,:,1) - atfp * rdt * r1_rau0 * ( emp_b(:,:) - emp(:,:) & 224 & -rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1)225 & - rnf_b(:,:) + rnf(:,:) ) * tmask(:,:,1) 225 226 ELSE ! if ice shelf melting 226 227 DO jj = 1,jpj 227 228 DO ji = 1,jpi 228 jk= mikt(ji,jj)229 fse3t_b(ji,jj, jk) = fse3t_b(ji,jj,jk) - atfp * rdt * r1_rau0&229 ikt = mikt(ji,jj) 230 fse3t_b(ji,jj,ikt) = fse3t_b(ji,jj,ikt) - atfp * rdt * r1_rau0 & 230 231 & * ( (emp_b(ji,jj) - emp(ji,jj) ) & 231 232 & - (rnf_b(ji,jj) - rnf(ji,jj) ) & 232 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj, jk)233 & + (fwfisf_b(ji,jj) - fwfisf(ji,jj)) ) * tmask(ji,jj,ikt) 233 234 END DO 234 235 END DO -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
r5930 r6012 231 231 IF( ioptio > 1 .OR. ( ioptio == 0 .AND. .NOT. lk_c1d ) ) & 232 232 & CALL ctl_stop( ' Choose only one surface pressure gradient scheme' ) 233 IF( ln_dynspg_ ts.AND. ln_isfcav ) &234 & CALL ctl_stop( ' dynspg_ tsnot tested with ice shelf cavity ' )233 IF( ln_dynspg_exp .AND. ln_isfcav ) & 234 & CALL ctl_stop( ' dynspg_exp not tested with ice shelf cavity ' ) 235 235 ! 236 236 IF( ln_dynspg_exp) nspg = 0 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5930 r6012 145 145 INTEGER :: ji, jj, jk, jn ! dummy loop indices 146 146 INTEGER :: ikbu, ikbv, noffset ! local integers 147 INTEGER :: iktu, iktv ! local integers 147 148 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 148 149 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - … … 384 385 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 385 386 DO ji = fs_2, fs_jpim1 386 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)387 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)387 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 388 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 388 389 END DO 389 390 END DO … … 413 414 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 414 415 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 416 ! 417 ! ! Add top stress contribution from baroclinic velocities: 418 IF (ln_bt_fw) THEN 419 DO jj = 2, jpjm1 420 DO ji = fs_2, fs_jpim1 ! vector opt. 421 iktu = miku(ji,jj) 422 iktv = mikv(ji,jj) 423 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 424 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 425 END DO 426 END DO 427 ELSE 428 DO jj = 2, jpjm1 429 DO ji = fs_2, fs_jpim1 ! vector opt. 430 iktu = miku(ji,jj) 431 iktv = mikv(ji,jj) 432 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 433 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 434 END DO 435 END DO 436 ENDIF 437 ! 438 ! Note that the "unclipped" top friction parameter is used even with explicit drag 439 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * tfrua(:,:) * zwx(:,:) 440 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * tfrva(:,:) * zwy(:,:) 415 441 ! 416 442 IF (ln_bt_fw) THEN ! Add wind forcing … … 544 570 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 545 571 DO ji = 2, fs_jpim1 ! Vector opt. 546 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) &572 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 547 573 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 548 574 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 549 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) &575 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 550 576 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 551 577 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) … … 607 633 END DO 608 634 END DO 609 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)635 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 610 636 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 611 637 … … 622 648 DO jj = 2, jpjm1 623 649 DO ji = 2, jpim1 ! NO Vector Opt. 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) )650 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 651 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 652 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 653 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 654 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 655 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 630 656 END DO 631 657 END DO … … 661 687 DO jj = 2, jpjm1 662 688 DO ji = 2, jpim1 663 zx1 = z1_2 * umask(ji ,jj,1) * r1_e1e2u(ji ,jj) &689 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 664 690 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 665 691 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 666 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e1e2v(ji ,jj ) &692 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 667 693 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 668 694 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) … … 736 762 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 737 763 ! 764 ! Add top stresses: 765 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 766 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 767 ! 738 768 ! Surface pressure trend: 739 769 DO jj = 2, jpjm1 … … 755 785 & + zu_trd(ji,jj) & 756 786 & + zu_frc(ji,jj) ) & 757 & ) * umask(ji,jj,1)787 & ) * ssumask(ji,jj) 758 788 759 789 va_e(ji,jj) = ( vn_e(ji,jj) & … … 761 791 & + zv_trd(ji,jj) & 762 792 & + zv_frc(ji,jj) ) & 763 & ) * vmask(ji,jj,1)764 END DO 765 END DO 766 767 ELSE 793 & ) * ssvmask(ji,jj) 794 END DO 795 END DO 796 797 ELSE ! Flux form 768 798 DO jj = 2, jpjm1 769 799 DO ji = fs_2, fs_jpim1 ! vector opt. 770 800 771 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1))772 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1))801 zhura = ssumask(ji,jj)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj)) 802 zhvra = ssvmask(ji,jj)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj)) 773 803 774 804 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & … … 791 821 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 792 822 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 793 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) )794 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) )823 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 824 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 795 825 ! 796 826 ENDIF … … 827 857 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 828 858 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 829 ELSE 859 ELSE ! Sum transports 830 860 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 831 861 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) … … 881 911 END DO 882 912 ! Save barotropic velocities not transport: 883 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )884 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )913 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 914 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 885 915 ENDIF 886 916 ! … … 897 927 ! 898 928 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 899 IF ( Agrif_NbStepint() .EQ.0 ) THEN929 IF ( Agrif_NbStepint() == 0 ) THEN 900 930 ub2_i_b(:,:) = 0.e0 901 931 vb2_i_b(:,:) = 0.e0 -
branches/2015/dev_MetOffice_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynzad.F90
r5836 r6012 90 90 DO jj = 2, jpjm1 ! vertical momentum advection at w-point 91 91 DO ji = fs_2, fs_jpim1 ! vector opt. 92 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) -un(ji,jj,jk) )93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) -vn(ji,jj,jk) )92 zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( un(ji,jj,jk-1) - un(ji,jj,jk) ) 93 zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( vn(ji,jj,jk-1) - vn(ji,jj,jk) ) 94 94 END DO 95 95 END DO
Note: See TracChangeset
for help on using the changeset viewer.
