- Timestamp:
- 2015-12-16T16:44:35+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r6060 r6069 45 45 USE wrk_nemo ! Memory Allocation 46 46 USE timing ! Timing 47 USE iom 47 48 48 49 IMPLICIT NONE … … 129 130 INTEGER :: ioptio = 0 ! temporary integer 130 131 INTEGER :: ios ! Local integer output status for namelist read 132 !! 133 INTEGER :: ji, jj, jk, ikt ! dummy loop indices ISF 134 REAL(wp) :: znad 135 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop, zrhd ! hypothesys on isf density 136 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_isf ! density at bottom of ISF 137 REAL(wp), POINTER, DIMENSION(:,:) :: ziceload ! density at bottom of ISF 131 138 !! 132 139 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & … … 189 196 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 190 197 ! 191 ! initialisation of ice load 192 riceload(:,:)=0.0 198 ! initialisation of ice shelf load 199 IF ( .NOT. ln_isfcav ) riceload(:,:)=0.0 200 IF ( ln_isfcav ) THEN 201 CALL wrk_alloc( jpi,jpj, 2, ztstop) 202 CALL wrk_alloc( jpi,jpj,jpk, zrhd ) 203 CALL wrk_alloc( jpi,jpj, zrhdtop_isf, ziceload) 204 ! 205 IF(lwp) WRITE(numout,*) 206 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 207 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~' 208 209 ! To use density and not density anomaly 210 znad=1._wp 211 212 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 213 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 214 215 ! compute density of the water displaced by the ice shelf 216 DO jk = 1, jpk 217 CALL eos(ztstop(:,:,:),fsdept_n(:,:,jk),zrhd(:,:,jk)) 218 END DO 219 220 ! compute rhd at the ice/oce interface (ice shelf side) 221 CALL eos(ztstop,risfdep,zrhdtop_isf) 222 223 ! Surface value + ice shelf gradient 224 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 225 ! divided by 2 later 226 ziceload = 0._wp 227 DO jj = 1, jpj 228 DO ji = 1, jpi 229 ikt=mikt(ji,jj) 230 ziceload(ji,jj) = ziceload(ji,jj) + (znad + zrhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 231 DO jk=2,ikt-1 232 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhd(ji,jj,jk-1) + zrhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 233 & * (1._wp - tmask(ji,jj,jk)) 234 END DO 235 IF (ikt >= 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + zrhd(ji,jj,ikt-1)) & 236 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 237 END DO 238 END DO 239 riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 240 241 CALL wrk_dealloc( jpi,jpj, 2, ztstop) 242 CALL wrk_dealloc( jpi,jpj,jpk, zrhd ) 243 CALL wrk_dealloc( jpi,jpj, zrhdtop_isf, ziceload) 244 END IF 193 245 ! 194 246 END SUBROUTINE dyn_hpg_init … … 444 496 SUBROUTINE hpg_isf( kt ) 445 497 !!--------------------------------------------------------------------- 446 !! *** ROUTINE hpg_ sco***498 !! *** ROUTINE hpg_isf *** 447 499 !! 448 500 !! ** Method : s-coordinate case. Jacobian scheme. … … 463 515 INTEGER, INTENT(in) :: kt ! ocean time-step index 464 516 !! 465 INTEGER :: ji, jj, jk, ik u, ikv, ikt, iktp1i, iktp1j! dummy loop indices466 REAL(wp) :: zcoef0, zuap, zvap, znad , ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1! temporary scalars467 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj , zrhd517 INTEGER :: ji, jj, jk, ikt, iktp1i, iktp1j ! dummy loop indices 518 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 519 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 468 520 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 469 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 470 !!---------------------------------------------------------------------- 471 ! 472 CALL wrk_alloc( jpi,jpj, 2, ztstop ) 473 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 474 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj ) 475 ! 476 IF( kt == nit000 ) THEN 477 IF(lwp) WRITE(numout,*) 478 IF(lwp) WRITE(numout,*) 'dyn:hpg_isf : hydrostatic pressure gradient trend for ice shelf' 479 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OPA original scheme used' 480 ENDIF 481 ! 521 REAL(wp), POINTER, DIMENSION(:,:) :: zrhdtop_oce 522 !!---------------------------------------------------------------------- 523 ! 524 CALL wrk_alloc( jpi,jpj, 2, ztstop) 525 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj) 526 CALL wrk_alloc( jpi,jpj, zrhdtop_oce ) 527 ! 528 ! Local constant initialization 482 529 zcoef0 = - grav * 0.5_wp 483 IF( ln_linssh ) THEN ; znad = 0._wp ! Fixed volume: density anomaly 484 ELSE ; znad = 1._wp ! Variable volume: density 485 ENDIF 486 zhpi(:,:,:) = 0._wp 487 zhpj(:,:,:) = 0._wp 530 531 ! To use density and not density anomaly 532 znad=1._wp 533 534 ! iniitialised to 0. zhpi zhpi 535 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 536 537 ! compute rhd at the ice/oce interface (ocean side) 538 ! usefull to reduce residual current in the test case ISOMIP with no melting 539 DO ji=1,jpi 540 DO jj=1,jpj 541 ikt=mikt(ji,jj) 542 ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 543 ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 544 END DO 545 END DO 546 CALL eos( ztstop, risfdep, zrhdtop_oce ) 488 547 489 548 !================================================================================== 490 !=====Compute iceload and contribution of the half first wet layer ================= 491 !=================================================================================== 492 493 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 494 ztstop(:,:,jp_tem) = -1.9_wp 495 ztstop(:,:,jp_sal) = 34.4_wp 496 497 !!gm I have the feeling that a much simplier and faster computation can be performed... 498 !!gm ====>>>> We have to discuss ! 499 500 !!gm below, faster to compute the ISF density in zrhd and remplace rhd value where tmask=0 501 !!gm furthermore, this calculation does not depends on time : do it at the first time-step only.... 502 503 ! compute density of the water displaced by the ice shelf 504 zrhd(:,:,:) = rhd(:,:,:) ! save rhd 505 DO jk = 1, jpk 506 zdept(:,:) = gdept_1d(jk) 507 CALL eos( ztstop(:,:,:), zdept(:,:), rhd(:,:,jk) ) 508 END DO 509 WHERE( tmask(:,:,:) == 1._wp ) 510 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 511 END WHERE 512 513 ! compute rhd at the ice/oce interface (ice shelf side) 514 CALL eos( ztstop, risfdep, zrhdtop_isf ) 515 516 ! compute rhd at the ice/oce interface (ocean side) 517 DO ji = 1, jpi 518 DO jj = 1, jpj 519 ikt = mikt(ji,jj) 520 ztstop(ji,jj,jp_tem) = tsn(ji,jj,ikt,jp_tem) 521 ztstop(ji,jj,jp_sal) = tsn(ji,jj,ikt,jp_sal) 522 END DO 523 END DO 524 CALL eos( ztstop, risfdep, zrhdtop_oce ) 525 ! 526 ! Surface value + ice shelf gradient 527 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 528 ziceload = 0._wp 529 DO jj = 1, jpj 530 DO ji = 1, jpi ! vector opt. 531 ikt = mikt(ji,jj) 532 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * e3w_n(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 533 DO jk = 2, ikt-1 534 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * e3w_n(ji,jj,jk) & 535 & * (1._wp - tmask(ji,jj,jk)) 536 END DO 537 IF( ikt >= 2 ) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 538 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 539 END DO 540 END DO 541 riceload(:,:) = 0._wp ; riceload(:,:) = ziceload(:,:) ! need to be saved for diaar5 542 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 549 !===== Compute surface value ===================================================== 550 !================================================================================== 543 551 DO jj = 2, jpjm1 544 552 DO ji = fs_2, fs_jpim1 ! vector opt. … … 548 556 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 549 557 ! we assume ISF is in isostatic equilibrium 550 zhpi(ji,jj,1) = zcoef0 * ( & 551 & 0.5_wp * e3w_n(ji+1,jj,iktp1i) * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 552 & - 0.5_wp * e3w_n(ji ,jj,ikt ) * ( 2._wp * znad + rhd(ji ,jj,ikt ) + zrhdtop_oce(ji ,jj) ) & 553 & + ( ziceload(ji+1,jj) - ziceload(ji,jj) ) ) * r1_e1u(ji,jj) 554 zhpj(ji,jj,1) = zcoef0 * ( & 555 & 0.5_wp * e3w_n(ji,jj+1,iktp1j) * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 556 & - 0.5_wp * e3w_n(ji,jj ,ikt ) * ( 2._wp * znad + rhd(ji,jj ,ikt ) + zrhdtop_oce(ji,jj ) ) & 557 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) * r1_e2v(ji,jj) 558 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj,iktp1i) & 559 & * ( 2._wp * znad + rhd(ji+1,jj,iktp1i) + zrhdtop_oce(ji+1,jj) ) & 560 & - 0.5_wp * e3w_n(ji,jj,ikt) & 561 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 562 & + ( riceload(ji+1,jj) - riceload(ji,jj)) ) 563 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji,jj+1,iktp1j) & 564 & * ( 2._wp * znad + rhd(ji,jj+1,iktp1j) + zrhdtop_oce(ji,jj+1) ) & 565 & - 0.5_wp * e3w_n(ji,jj,ikt) & 566 & * ( 2._wp * znad + rhd(ji,jj,ikt) + zrhdtop_oce(ji,jj) ) & 567 & + ( riceload(ji,jj+1) - riceload(ji,jj)) ) 558 568 ! s-coordinate pressure gradient correction (=0 if z coordinate) 559 569 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & … … 567 577 END DO 568 578 !================================================================================== 569 !===== Compute partial cell contribution for the top cell =========================570 !==================================================================================571 DO jj = 2, jpjm1572 DO ji = fs_2, fs_jpim1 ! vector opt.573 iku = miku(ji,jj)574 zpshpi(ji,jj) = 0._wp575 zpshpj(ji,jj) = 0._wp576 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))577 ! u direction578 IF( iku > 1 ) THEN579 ! case iku580 zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu &581 & * ( rhd(ji+1,jj,iku) + rhd(ji,jj,iku) &582 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad )583 ! corrective term ( = 0 if z coordinate )584 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj)585 ! zhpi will be added in interior loop586 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap587 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure588 IF( mbku(ji,jj) == iku + 1 ) zpshpi(ji,jj) = zhpi(ji,jj,iku)589 590 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps)591 zhpiint = zcoef0 * r1_e1u(ji,jj) &592 & * ( e3w_n(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) &593 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) &594 & - e3w_n(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) &595 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) )596 zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint597 END IF598 599 ! v direction600 ikv = mikv(ji,jj)601 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))602 IF( ikv > 1 ) THEN603 ! case ikv604 zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv &605 & * ( rhd(ji,jj+1,ikv) + rhd(ji,jj,ikv) &606 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad )607 ! corrective term ( = 0 if z coordinate )608 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj)609 ! zhpi will be added in interior loop610 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap611 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure612 IF( mbkv(ji,jj) == ikv + 1 ) zpshpj(ji,jj) = zhpj(ji,jj,ikv)613 614 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps)615 zhpjint = zcoef0 * r1_e2v(ji,jj) &616 & * ( e3w_n(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) &617 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) &618 & - e3w_n(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) &619 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) )620 zhpj(ji,jj,ikv+1) = zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint621 ENDIF622 END DO623 END DO624 625 !==================================================================================626 579 !===== Compute interior value ===================================================== 627 580 !================================================================================== 628 629 DO j j = 2, jpjm1630 DO j i = fs_2, fs_jpim1 ! vector opt.631 DO j k = 2, jpkm1581 ! interior value (2=<jk=<jpkm1) 582 DO jk = 2, jpkm1 583 DO jj = 2, jpjm1 584 DO ji = fs_2, fs_jpim1 ! vector opt. 632 585 ! hydrostatic pressure gradient along s-surfaces 633 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 634 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 635 & + zcoef0 * r1_e1u(ji,jj) & 636 & * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 637 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 638 & - e3w_n(ji ,jj,jk) * ( (rhd(ji ,jj,jk ) + znad) & 639 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 586 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 587 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) * wmask(ji+1,jj,jk) & 588 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) * wmask(ji ,jj,jk) ) 589 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 590 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) * wmask(ji,jj+1,jk) & 591 & - e3w_n(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) * wmask(ji,jj ,jk) ) 640 592 ! s-coordinate pressure gradient correction 641 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 642 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 643 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 644 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 645 646 ! hydrostatic pressure gradient along s-surfaces 647 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 648 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 649 & + zcoef0 * r1_e2v(ji,jj) & 650 & * ( e3w_n(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 651 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 652 & - e3w_n(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 653 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 654 ! s-coordinate pressure gradient correction 655 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 656 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 657 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 593 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 594 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 595 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 596 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 658 597 ! add to the general momentum trend 659 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 598 ua(ji,jj,jk) = ua(ji,jj,jk) + (zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 599 va(ji,jj,jk) = va(ji,jj,jk) + (zhpj(ji,jj,jk) + zvap) * vmask(ji,jj,jk) 660 600 END DO 661 601 END DO 662 602 END DO 663 664 !================================================================================== 665 !===== Compute bottom cell contribution (partial cell) ============================ 666 !================================================================================== 667 668 DO jj = 2, jpjm1 669 DO ji = 2, jpim1 670 iku = mbku(ji,jj) 671 ikv = mbkv(ji,jj) 672 673 IF (iku .GT. 1) THEN 674 ! remove old value (interior case) 675 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 676 & * ( gde3w_n(ji+1,jj ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 677 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 678 ! put new value 679 ! -zpshpi to avoid double contribution of the partial step in the top layer 680 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) * r1_e1u(ji,jj) 681 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 682 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 683 END IF 684 ! v direction 685 IF (ikv .GT. 1) THEN 686 ! remove old value (interior case) 687 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 688 & * ( gde3w_n(ji ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj) 689 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 690 ! put new value 691 ! -zpshpj to avoid double contribution of the partial step in the top layer 692 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj) 693 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 694 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 695 END IF 696 END DO 697 END DO 698 699 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 700 rhd = zrhd 701 ! 702 CALL wrk_dealloc( jpi,jpj,2, ztstop) 703 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 704 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 603 ! 604 CALL wrk_dealloc( jpi,jpj,2 , ztstop) 605 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj) 606 CALL wrk_dealloc( jpi,jpj , zrhdtop_oce ) 705 607 ! 706 608 END SUBROUTINE hpg_isf
Note: See TracChangeset
for help on using the changeset viewer.
