Changeset 5260 for branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2015-05-12T12:37:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4624 r5260 16 16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 18 19 !!---------------------------------------------------------------------- 19 20 … … 25 26 !! hpg_zps : z-coordinate plus partial steps (interpolation) 26 27 !! hpg_sco : s-coordinate (standard jacobian formulation) 28 !! hpg_isf : s-coordinate (sco formulation) adapted to ice shelf 27 29 !! hpg_djc : s-coordinate (Density Jacobian with Cubic polynomial) 28 30 !! hpg_prj : s-coordinate (Pressure Jacobian with Cubic polynomial) 29 31 !!---------------------------------------------------------------------- 30 32 USE oce ! ocean dynamics and tracers 33 USE sbc_oce ! surface variable (only for the flag with ice shelf) 31 34 USE dom_oce ! ocean space and time domain 32 35 USE phycst ! physical constants 33 USE trdmod ! ocean dynamics trends 34 USE trdmod_oce ! ocean variables trends 36 USE trd_oce ! trends: ocean variables 37 USE trddyn ! trend manager: dynamics 38 ! 35 39 USE in_out_manager ! I/O manager 36 40 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 41 USE lbclnk ! lateral boundary condition 38 42 USE lib_mpp ! MPP library 43 USE eosbn2 ! compute density 39 44 USE wrk_nemo ! Memory Allocation 40 45 USE timing ! Timing … … 52 57 LOGICAL , PUBLIC :: ln_hpg_djc !: s-coordinate (Density Jacobian with Cubic polynomial) 53 58 LOGICAL , PUBLIC :: ln_hpg_prj !: s-coordinate (Pressure Jacobian scheme) 59 LOGICAL , PUBLIC :: ln_hpg_isf !: s-coordinate similar to sco modify for isf 54 60 LOGICAL , PUBLIC :: ln_dynhpg_imp !: semi-implicite hpg flag 55 61 … … 74 80 !! 75 81 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 !! - Save the trend(l_trddyn=T)82 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 77 83 !!---------------------------------------------------------------------- 78 84 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 94 100 CASE ( 3 ) ; CALL hpg_djc ( kt ) ! s-coordinate (Density Jacobian with Cubic polynomial) 95 101 CASE ( 4 ) ; CALL hpg_prj ( kt ) ! s-coordinate (Pressure Jacobian scheme) 102 CASE ( 5 ) ; CALL hpg_isf ( kt ) ! s-coordinate similar to sco modify for ice shelf 96 103 END SELECT 97 104 ! … … 99 106 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 107 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )108 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 102 109 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 110 ENDIF … … 125 132 !! 126 133 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 127 & ln_hpg_djc, ln_hpg_prj, ln_ dynhpg_imp134 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, ln_dynhpg_imp 128 135 !!---------------------------------------------------------------------- 129 136 ! … … 145 152 WRITE(numout,*) ' z-coord. - partial steps (interpolation) ln_hpg_zps = ', ln_hpg_zps 146 153 WRITE(numout,*) ' s-coord. (standard jacobian formulation) ln_hpg_sco = ', ln_hpg_sco 154 WRITE(numout,*) ' s-coord. (standard jacobian formulation) for isf ln_hpg_isf = ', ln_hpg_isf 147 155 WRITE(numout,*) ' s-coord. (Density Jacobian: Cubic polynomial) ln_hpg_djc = ', ln_hpg_djc 148 156 WRITE(numout,*) ' s-coord. (Pressure Jacobian: Cubic polynomial) ln_hpg_prj = ', ln_hpg_prj … … 155 163 & either ln_hpg_sco or ln_hpg_prj instead') 156 164 ! 157 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj ) ) &165 IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 158 166 & CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 159 167 & the standard jacobian formulation hpg_sco or & 160 168 & the pressure jacobian formulation hpg_prj') 169 170 IF( ln_hpg_isf .AND. .NOT. ln_isfcav ) & 171 & CALL ctl_stop( ' hpg_isf not available if ln_isfcav = false ' ) 172 IF( .NOT. ln_hpg_isf .AND. ln_isfcav ) & 173 & CALL ctl_stop( 'Only hpg_isf has been corrected to work with ice shelf cavity.' ) 161 174 ! 162 175 ! ! Set nhpg from ln_hpg_... flags … … 166 179 IF( ln_hpg_djc ) nhpg = 3 167 180 IF( ln_hpg_prj ) nhpg = 4 181 IF( ln_hpg_isf ) nhpg = 5 168 182 ! 169 183 ! ! Consistency check … … 174 188 IF( ln_hpg_djc ) ioptio = ioptio + 1 175 189 IF( ln_hpg_prj ) ioptio = ioptio + 1 190 IF( ln_hpg_isf ) ioptio = ioptio + 1 176 191 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 192 ! 193 ! initialisation of ice load 194 riceload(:,:)=0.0 177 195 ! 178 196 END SUBROUTINE dyn_hpg_init … … 315 333 316 334 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 317 # if defined key_vectopt_loop318 jj = 1319 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)320 # else321 335 DO jj = 2, jpjm1 322 336 DO ji = 2, jpim1 323 # endif324 337 iku = mbku(ji,jj) 325 338 ikv = mbkv(ji,jj) … … 338 351 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 339 352 ENDIF 340 # if ! defined key_vectopt_loop 341 END DO 342 # endif 353 END DO 343 354 END DO 344 355 ! … … 346 357 ! 347 358 END SUBROUTINE hpg_zps 348 349 359 350 360 SUBROUTINE hpg_sco( kt ) … … 434 444 END SUBROUTINE hpg_sco 435 445 446 SUBROUTINE hpg_isf( kt ) 447 !!--------------------------------------------------------------------- 448 !! *** ROUTINE hpg_sco *** 449 !! 450 !! ** Method : s-coordinate case. Jacobian scheme. 451 !! The now hydrostatic pressure gradient at a given level, jk, 452 !! is computed by taking the vertical integral of the in-situ 453 !! density gradient along the model level from the suface to that 454 !! level. s-coordinates (ln_sco): a corrective term is added 455 !! to the horizontal pressure gradient : 456 !! zhpi = grav ..... + 1/e1u mi(rhd) di[ grav dep3w ] 457 !! zhpj = grav ..... + 1/e2v mj(rhd) dj[ grav dep3w ] 458 !! add it to the general momentum trend (ua,va). 459 !! ua = ua - 1/e1u * zhpi 460 !! va = va - 1/e2v * zhpj 461 !! iceload is added and partial cell case are added to the top and bottom 462 !! 463 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 464 !!---------------------------------------------------------------------- 465 INTEGER, INTENT(in) :: kt ! ocean time-step index 466 !! 467 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 468 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars 469 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zrhd 470 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 484 ! Local constant initialization 485 zcoef0 = - grav * 0.5_wp 486 ! To use density and not density anomaly 487 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 488 ! ELSE ; znad = 0._wp ! Fixed volume 489 ! ENDIF 490 znad=1._wp 491 ! iniitialised to 0. zhpi zhpi 492 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 493 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) 541 DO jj = 2, jpjm1 542 DO ji = fs_2, fs_jpim1 ! vector opt. 543 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 544 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 545 ! 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) ) ) 556 ! s-coordinate pressure gradient correction (=0 if z coordinate) 557 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 558 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) 559 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 560 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 561 ! add to the general momentum trend 562 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 563 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 564 END DO 565 END DO 566 !================================================================================== 567 !===== Compute partial cell contribution for the top cell ========================= 568 !================================================================================== 569 DO jj = 2, jpjm1 570 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_wp 573 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 direction 575 IF ( iku .GT. 1 ) THEN 576 ! case iku 577 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 loop 583 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 584 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 585 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) - zhpiint 594 END IF 595 596 ! v direction 597 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 ) THEN 600 ! case ikv 601 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 loop 607 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 608 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 609 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) - zhpjint 618 END IF 619 END DO 620 END DO 621 622 !================================================================================== 623 !===== Compute interior value ===================================================== 624 !================================================================================== 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 630 ! 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) ) 638 ! 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) 656 ! add to the general momentum trend 657 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 658 END DO 659 END DO 660 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) 703 ! 704 END SUBROUTINE hpg_isf 705 706 436 707 SUBROUTINE hpg_djc( kt ) 437 708 !!--------------------------------------------------------------------- … … 671 942 !! 672 943 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 673 !! - Save the trend (l_trddyn=T)674 !!675 944 !!---------------------------------------------------------------------- 676 945 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear … … 687 956 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 688 957 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 958 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_n, zsshv_n 689 959 !!---------------------------------------------------------------------- 690 960 ! 691 961 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 692 962 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 963 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 693 964 ! 694 965 IF( kt == nit000 ) THEN … … 724 995 725 996 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 726 DO jj = 1, jpj; DO ji = 1, jpi 727 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 728 END DO ; END DO 729 730 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 731 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 732 END DO ; END DO ; END DO 733 734 fsp(:,:,:) = zrhh(:,:,:) 997 DO jj = 1, jpj 998 DO ji = 1, jpi 999 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 1000 END DO 1001 END DO 1002 1003 DO jk = 2, jpk 1004 DO jj = 1, jpj 1005 DO ji = 1, jpi 1006 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 1007 END DO 1008 END DO 1009 END DO 1010 1011 fsp(:,:,:) = zrhh (:,:,:) 735 1012 xsp(:,:,:) = zdept(:,:,:) 736 1013 … … 765 1042 766 1043 ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 1044 1045 ! Prepare zsshu_n and zsshv_n 767 1046 DO jj = 2, jpjm1 768 1047 DO ji = 2, jpim1 769 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshu_n for ztilde compilation 770 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshn(ji,jj) * znad) ! probable bug: changed from sshv_n for ztilde compilation 1048 zsshu_n(ji,jj) = (e12u(ji,jj) * sshn(ji,jj) + e12u(ji+1, jj) * sshn(ji+1,jj)) * & 1049 & r1_e12u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1050 zsshv_n(ji,jj) = (e12v(ji,jj) * sshn(ji,jj) + e12v(ji+1, jj) * sshn(ji,jj+1)) * & 1051 & r1_e12v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1052 END DO 1053 END DO 1054 1055 DO jj = 2, jpjm1 1056 DO ji = 2, jpim1 1057 zu(ji,jj,1) = - ( fse3u(ji,jj,1) - zsshu_n(ji,jj) * znad) 1058 zv(ji,jj,1) = - ( fse3v(ji,jj,1) - zsshv_n(ji,jj) * znad) 771 1059 END DO 772 1060 END DO … … 930 1218 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 931 1219 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1220 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 932 1221 ! 933 1222 END SUBROUTINE hpg_prj 934 1223 1224 935 1225 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 936 1226 !!---------------------------------------------------------------------- … … 940 1230 !! 941 1231 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1232 !! 942 1233 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 943 !!944 1234 !!---------------------------------------------------------------------- 945 1235 IMPLICIT NONE … … 949 1239 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 950 1240 ! 2: Linear 951 952 ! Local Variables 1241 ! 953 1242 INTEGER :: ji, jj, jk ! dummy loop indices 954 1243 INTEGER :: jpi, jpj, jpkm1 … … 1040 1329 ENDIF 1041 1330 1042 1043 1331 END SUBROUTINE cspline 1044 1332 … … 1050 1338 !! ** Purpose : 1-d linear interpolation 1051 1339 !! 1052 !! ** Method : 1053 !! interpolation is straight forward 1340 !! ** Method : interpolation is straight forward 1054 1341 !! extrapolation is also permitted (no value limit) 1055 !!1056 1342 !!---------------------------------------------------------------------- 1057 1343 IMPLICIT NONE … … 1070 1356 END FUNCTION interp1 1071 1357 1358 1072 1359 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1073 1360 !!---------------------------------------------------------------------- … … 1133 1420 END FUNCTION integ_spline 1134 1421 1135 1136 1422 !!====================================================================== 1137 1423 END MODULE dynhpg
Note: See TracChangeset
for help on using the changeset viewer.
