- Timestamp:
- 2015-03-31T19:58:23+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5120 r5189 44 44 USE wrk_nemo ! Memory Allocation 45 45 USE timing ! Timing 46 USE iom 46 47 47 48 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, & … … 191 198 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 192 199 ! 193 ! initialisation of ice load 194 riceload(:,:)=0.0 200 ! initialisation of ice shelf load 201 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 202 IF ( ln_isfcav ) THEN 203 CALL wrk_alloc( jpi,jpj, 2, ztstop) 204 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 205 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 206 ! 207 IF(lwp) WRITE(numout,*) 208 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 209 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 210 211 ! To use density and not density anomaly 212 znad=1._wp 213 214 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 215 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 216 217 ! compute density of the water displaced by the ice shelf 218 DO jk = 1, jpk 219 CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 220 END DO 221 222 ! compute rhd at the ice/oce interface (ice shelf side) 223 CALL eos(ztstop,risfdep,zrhdtop_isf) 224 225 ! Surface value + ice shelf gradient 226 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 227 ziceload = 0._wp 228 DO jj = 1, jpj 229 DO ji = 1, jpi ! vector opt. 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 .GE. 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 195 246 ! 196 247 END SUBROUTINE dyn_hpg_init … … 465 516 INTEGER, INTENT(in) :: kt ! ocean time-step index 466 517 !! 467 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j 468 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars469 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd518 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 519 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 520 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 470 521 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 471 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 472 !!---------------------------------------------------------------------- 473 ! 474 CALL wrk_alloc( jpi,jpj, 2, ztstop) 475 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 476 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 477 ! 478 IF( kt == nit000 ) THEN 479 IF(lwp) WRITE(numout,*) 480 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 481 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 482 ENDIF 483 522 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce, zpshpi, zpshpj 523 !!---------------------------------------------------------------------- 524 ! 525 CALL wrk_alloc( jpi,jpj, 2, ztstop) 526 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 527 CALL wrk_alloc( jpi,jpj, zrhdtop_oce, zpshpi, zpshpj) 528 ! 484 529 ! Local constant initialization 485 530 zcoef0 = - grav * 0.5_wp 531 486 532 ! To use density and not density anomaly 487 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume488 ! ELSE ; znad = 0._wp ! Fixed volume489 ! ENDIF490 533 znad=1._wp 534 491 535 ! iniitialised to 0. zhpi zhpi 492 536 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 493 537 494 !================================================================================== 495 !=====Compute iceload and contribution of the half first wet layer ================= 496 !=================================================================================== 497 498 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 499 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 500 501 ! compute density of the water displaced by the ice shelf 502 zrhd = rhd ! save rhd 503 DO jk = 1, jpk 504 zdept(:,:)=gdept_1d(jk) 505 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 506 END DO 507 WHERE ( tmask(:,:,:) == 1._wp) 508 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 509 END WHERE 510 511 ! compute rhd at the ice/oce interface (ice shelf side) 512 CALL eos(ztstop,risfdep,zrhdtop_isf) 513 514 ! compute rhd at the ice/oce interface (ocean side) 515 DO ji=1,jpi 516 DO jj=1,jpj 517 ikt=mikt(ji,jj) 518 ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 519 ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 520 END DO 521 END DO 522 CALL eos(ztstop,risfdep,zrhdtop_oce) 523 ! 524 ! Surface value + ice shelf gradient 525 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 526 ziceload = 0._wp 527 DO jj = 1, jpj 528 DO ji = 1, jpi ! vector opt. 529 ikt=mikt(ji,jj) 530 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 531 DO jk=2,ikt-1 532 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 533 & * (1._wp - tmask(ji,jj,jk)) 534 END DO 535 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 536 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 537 END DO 538 END DO 539 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 540 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 538 ! compute rhd at the ice/oce interface (ocean side) usefull to reduce unrealistic water current if homogene water 539 ! TO BE DISCUSS do we compute density at w point for the ocean/ice atmosphere 540 ! DO ji=1,jpi 541 ! DO jj=1,jpj 542 ! ikt=mikt(ji,jj) 543 ! ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 544 ! ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 545 ! END DO 546 ! END DO 547 ! CALL eos(ztstop,risfdep,zrhdtop_oce) 548 541 549 DO jj = 2, jpjm1 542 550 DO ji = fs_2, fs_jpim1 ! vector opt. … … 544 552 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 545 553 ! we assume ISF is in isostatic equilibrium 546 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) & 547 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 548 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 549 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 550 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)) ) 551 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) & 552 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 553 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 554 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 555 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 554 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,iktp1i) * ( znad + rhd(ji+1,jj ,iktp1i) ) & 555 & - fse3w(ji ,jj ,ikt ) * ( znad + rhd(ji ,jj ,ikt ) ) & 556 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 557 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,iktp1j) * ( znad + rhd(ji ,jj+1,iktp1j) ) & 558 & - fse3w(ji ,jj ,ikt ) * ( znad + rhd(ji ,jj ,ikt ) ) & 559 & + ( riceload(ji,jj+1) - riceload(ji,jj) ) ) 560 ! zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) & 561 ! & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 562 ! & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 563 ! & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 564 ! & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 565 ! zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) & 566 ! & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 567 ! & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 568 ! & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 569 ! & + ( riceload(ji,jj+1) - riceload(ji,jj) ) ) 556 570 ! s-coordinate pressure gradient correction (=0 if z coordinate) 557 571 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 564 578 END DO 565 579 END DO 566 !==================================================================================567 !===== Compute partial cell contribution for the top cell =========================568 !==================================================================================569 DO jj = 2, jpjm1570 DO ji = fs_2, fs_jpim1 ! vector opt.571 iku = miku(ji,jj) ;572 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp573 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))574 ! u direction575 IF ( iku .GT. 1 ) THEN576 ! case iku577 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu &578 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) &579 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )580 ! corrective term ( = 0 if z coordinate )581 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj)582 ! zhpi will be added in interior loop583 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap584 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure585 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku)586 587 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)588 zhpiint = zcoef0 / e1u(ji,jj) &589 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &590 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &591 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &592 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )593 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint594 END IF595 596 ! v direction597 ikv = mikv(ji,jj)598 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))599 IF ( ikv .GT. 1 ) THEN600 ! case ikv601 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv &602 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) &603 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )604 ! corrective term ( = 0 if z coordinate )605 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj)606 ! zhpi will be added in interior loop607 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap608 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure609 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv)610 611 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)612 zhpjint = zcoef0 / e2v(ji,jj) &613 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &614 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &615 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &616 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )617 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint618 END IF619 END DO620 END DO621 580 622 581 !================================================================================== 623 582 !===== Compute interior value ===================================================== 624 583 !================================================================================== 625 626 DO jj = 2, jpjm1 627 DO ji = fs_2, fs_jpim1 ! vector opt. 628 iku=miku(ji,jj); ikv=mikv(ji,jj) 629 DO jk = 2, jpkm1 584 ! interior value (2=<jk=<jpkm1) 585 DO jk = 2, jpkm1 586 DO jj = 2, jpjm1 587 DO ji = fs_2, fs_jpim1 ! vector opt. 630 588 ! hydrostatic pressure gradient along s-surfaces 631 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 632 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 633 & + zcoef0 / e1u(ji,jj) & 634 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 635 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 636 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 637 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 589 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 590 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 591 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 592 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 593 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 594 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 638 595 ! s-coordinate pressure gradient correction 639 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 640 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 641 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 642 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 643 644 ! hydrostatic pressure gradient along s-surfaces 645 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 646 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 647 & + zcoef0 / e2v(ji,jj) & 648 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 649 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 650 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 651 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 652 ! s-coordinate pressure gradient correction 653 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 654 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 655 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 596 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 597 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 598 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 599 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 656 600 ! add to the general momentum trend 657 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 601 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 602 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 658 603 END DO 659 604 END DO 660 605 END DO 661 662 !================================================================================== 663 !===== Compute bottom cell contribution (partial cell) ============================ 664 !================================================================================== 665 666 DO jj = 2, jpjm1 667 DO ji = 2, jpim1 668 iku = mbku(ji,jj) 669 ikv = mbkv(ji,jj) 670 671 IF (iku .GT. 1) THEN 672 ! remove old value (interior case) 673 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 674 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 675 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 676 ! put new value 677 ! -zpshpi to avoid double contribution of the partial step in the top layer 678 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 679 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 680 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 681 END IF 682 ! v direction 683 IF (ikv .GT. 1) THEN 684 ! remove old value (interior case) 685 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 686 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 687 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 688 ! put new value 689 ! -zpshpj to avoid double contribution of the partial step in the top layer 690 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 691 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 692 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 693 END IF 694 END DO 695 END DO 696 697 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 698 rhd = zrhd 699 ! 700 CALL wrk_dealloc( jpi,jpj,2, ztstop) 701 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 702 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 606 ! 607 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 608 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 609 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce, zpshpi, zpshpj) 703 610 ! 704 611 END SUBROUTINE hpg_isf
Note: See TracChangeset
for help on using the changeset viewer.
