- Timestamp:
- 2019-11-22T15:29:17+01:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/SBC/sbcisf.F90
r11536 r11949 75 75 CONTAINS 76 76 77 SUBROUTINE sbc_isf( kt )77 SUBROUTINE sbc_isf( kt, Kmm ) 78 78 !!--------------------------------------------------------------------- 79 79 !! *** ROUTINE sbc_isf *** … … 89 89 !!---------------------------------------------------------------------- 90 90 INTEGER, INTENT(in) :: kt ! ocean time step 91 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 91 92 ! 92 93 INTEGER :: ji, jj, jk ! loop index … … 102 103 CASE ( 1 ) ! realistic ice shelf formulation 103 104 ! compute T/S/U/V for the top boundary layer 104 CALL sbc_isf_tbl(ts n(:,:,:,jp_tem),ttbl(:,:),'T')105 CALL sbc_isf_tbl(ts n(:,:,:,jp_sal),stbl(:,:),'T')106 CALL sbc_isf_tbl(u n(:,:,:) ,utbl(:,:),'U')107 CALL sbc_isf_tbl(v n(:,:,:) ,vtbl(:,:),'V')105 CALL sbc_isf_tbl(ts(:,:,:,jp_tem,Kmm),ttbl(:,:),'T',Kmm) 106 CALL sbc_isf_tbl(ts(:,:,:,jp_sal,Kmm),stbl(:,:),'T',Kmm) 107 CALL sbc_isf_tbl(uu(:,:,:,Kmm) ,utbl(:,:),'U',Kmm) 108 CALL sbc_isf_tbl(vv(:,:,:,Kmm) ,vtbl(:,:),'V',Kmm) 108 109 ! iom print 109 110 CALL iom_put('ttbl',ttbl(:,:)) … … 113 114 ! compute fwf and heat flux 114 115 ! compute fwf and heat flux 115 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt )116 IF( .NOT.l_isfcpl ) THEN ; CALL sbc_isf_cav (kt, Kmm) 116 117 ELSE ; qisf(:,:) = fwfisf(:,:) * rLfusisf ! heat flux 117 118 ENDIF … … 119 120 CASE ( 2 ) ! Beckmann and Goosse parametrisation 120 121 stbl(:,:) = soce 121 CALL sbc_isf_bg03(kt )122 CALL sbc_isf_bg03(kt, Kmm) 122 123 ! 123 124 CASE ( 3 ) ! specified runoff in depth (Mathiot et al., XXXX in preparation) … … 148 149 DO jj = 1,jpj 149 150 DO ji = 1,jpi 150 zdep(ji,jj)=gdepw _n(ji,jj,misfkt(ji,jj))151 zdep(ji,jj)=gdepw(ji,jj,misfkt(ji,jj),Kmm) 151 152 END DO 152 153 END DO … … 179 180 ikb = misfkb(ji,jj) 180 181 DO jk = ikt, ikb - 1 181 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t _n(ji,jj,jk)182 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t _n(ji,jj,jk)183 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t _n(ji,jj,jk)182 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 183 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 184 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) * e3t(ji,jj,jk,Kmm) 184 185 END DO 185 186 zfwfisf3d (ji,jj,jk) = zfwfisf3d (ji,jj,jk) + fwfisf (ji,jj) * r1_hisf_tbl(ji,jj) & 186 & * ralpha(ji,jj) * e3t _n(ji,jj,jk)187 & * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 187 188 zqhcisf3d (ji,jj,jk) = zqhcisf3d (ji,jj,jk) + zqhcisf2d(ji,jj) * r1_hisf_tbl(ji,jj) & 188 & * ralpha(ji,jj) * e3t _n(ji,jj,jk)189 & * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 189 190 zqlatisf3d(ji,jj,jk) = zqlatisf3d(ji,jj,jk) + qisf (ji,jj) * r1_hisf_tbl(ji,jj) & 190 & * ralpha(ji,jj) * e3t _n(ji,jj,jk)191 & * ralpha(ji,jj) * e3t(ji,jj,jk,Kmm) 191 192 END DO 192 193 END DO … … 251 252 252 253 253 SUBROUTINE sbc_isf_init 254 SUBROUTINE sbc_isf_init( Kmm ) 254 255 !!--------------------------------------------------------------------- 255 256 !! *** ROUTINE sbc_isf_init *** … … 263 264 !! 4 : specified fwf and heat flux forcing beneath the ice shelf 264 265 !!---------------------------------------------------------------------- 266 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 265 267 INTEGER :: ji, jj, jk ! loop index 266 268 INTEGER :: ik ! current level index … … 355 357 ik = 2 356 358 !!gm potential bug: use gdepw_0 not _n 357 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw _n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO359 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw(ji,jj,ik,Kmm) < rzisf_tbl(ji,jj) ) ; ik = ik + 1 ; END DO 358 360 misfkt(ji,jj) = ik-1 359 361 END DO … … 386 388 ikb = misfkt(ji,jj) 387 389 ! thickness of boundary layer at least the top level thickness 388 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t _n(ji,jj,ikt))390 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm)) 389 391 390 392 ! determine the deepest level influenced by the boundary layer 391 393 DO jk = ikt+1, mbkt(ji,jj) 392 IF( (SUM(e3t _n(ji,jj,ikt:jk-1)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk393 END DO 394 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t _n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.394 IF( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) < rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 395 END DO 396 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness. 395 397 misfkb(ji,jj) = ikb ! last wet level of the tbl 396 398 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 397 399 398 zhk = SUM( e3t _n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1399 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t _n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer400 zhk = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 401 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm) ! proportion of bottom cell influenced by boundary layer 400 402 END DO 401 403 END DO … … 411 413 412 414 413 SUBROUTINE sbc_isf_bg03( kt)415 SUBROUTINE sbc_isf_bg03( kt, Kmm ) 414 416 !!--------------------------------------------------------------------- 415 417 !! *** ROUTINE sbc_isf_bg03 *** … … 426 428 !!---------------------------------------------------------------------- 427 429 INTEGER, INTENT ( in ) :: kt 430 INTEGER, INTENT ( in ) :: Kmm ! ocean time level indices 428 431 ! 429 432 INTEGER :: ji, jj, jk ! dummy loop index … … 444 447 DO jk = misfkt(ji,jj),misfkb(ji,jj) 445 448 ! Calculate freezing temperature 446 zpress = grav*rau0*gdept _n(ji,jj,ik)*1.e-04449 zpress = grav*rau0*gdept(ji,jj,ik,Kmm)*1.e-04 447 450 CALL eos_fzp(stbl(ji,jj), zt_frz, zpress) 448 zt_sum = zt_sum + (ts n(ji,jj,jk,jp_tem)-zt_frz) * e3t_n(ji,jj,jk) * tmask(ji,jj,jk) ! sum temp451 zt_sum = zt_sum + (ts(ji,jj,jk,jp_tem,Kmm)-zt_frz) * e3t(ji,jj,jk,Kmm) * tmask(ji,jj,jk) ! sum temp 449 452 END DO 450 453 zt_ave = zt_sum/rhisf_tbl(ji,jj) ! calcul mean value … … 466 469 467 470 468 SUBROUTINE sbc_isf_cav( kt )471 SUBROUTINE sbc_isf_cav( kt, Kmm ) 469 472 !!--------------------------------------------------------------------- 470 473 !! *** ROUTINE sbc_isf_cav *** … … 480 483 !!--------------------------------------------------------------------- 481 484 INTEGER, INTENT(in) :: kt ! ocean time step 485 INTEGER, INTENT(in) :: Kmm ! ocean time level index 482 486 ! 483 487 INTEGER :: ji, jj ! dummy loop indices … … 520 524 521 525 ! compute gammat every where (2d) 522 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx )526 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm) 523 527 524 528 ! compute upward heat flux zhtflx and upward water flux zwflx … … 536 540 CASE ( 2 ) ! ISOMIP+ formulation (3 equations) for volume flux (Asay-Davis et al., 2015) 537 541 ! compute gammat every where (2d) 538 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx )542 CALL sbc_isf_gammats(zgammat, zgammas, zhtflx, zfwflx, Kmm) 539 543 540 544 ! compute upward heat flux zhtflx and upward water flux zwflx … … 600 604 601 605 602 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf )606 SUBROUTINE sbc_isf_gammats(pgt, pgs, pqhisf, pqwisf, Kmm ) 603 607 !!---------------------------------------------------------------------- 604 608 !! ** Purpose : compute the coefficient echange for heat flux … … 611 615 REAL(wp), DIMENSION(:,:), INTENT( out) :: pgt , pgs ! 612 616 REAL(wp), DIMENSION(:,:), INTENT(in ) :: pqhisf, pqwisf ! 617 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 613 618 ! 614 619 INTEGER :: ji, jj ! loop index … … 679 684 !!gm better to do it like in the new zdfric.F90 i.e. avm weighted Ri computation 680 685 !!gm moreover, use Max(rn2,0) to take care of static instabilities.... 681 zcoef = 0.5_wp / e3w _n(ji,jj,ikt+1)686 zcoef = 0.5_wp / e3w(ji,jj,ikt+1,Kmm) 682 687 ! ! shear of horizontal velocity 683 zdku = zcoef * ( u n(ji-1,jj ,ikt ) + un(ji,jj,ikt) &684 & -u n(ji-1,jj ,ikt+1) - un(ji,jj,ikt+1) )685 zdkv = zcoef * ( v n(ji ,jj-1,ikt ) + vn(ji,jj,ikt) &686 & -v n(ji ,jj-1,ikt+1) - vn(ji,jj,ikt+1) )688 zdku = zcoef * ( uu(ji-1,jj ,ikt ,Kmm) + uu(ji,jj,ikt ,Kmm) & 689 & -uu(ji-1,jj ,ikt+1,Kmm) - uu(ji,jj,ikt+1,Kmm) ) 690 zdkv = zcoef * ( vv(ji ,jj-1,ikt ,Kmm) + vv(ji,jj,ikt ,Kmm) & 691 & -vv(ji ,jj-1,ikt+1,Kmm) - vv(ji,jj,ikt+1,Kmm) ) 687 692 ! ! richardson number (minimum value set to zero) 688 693 zRc = rn2(ji,jj,ikt+1) / MAX( zdku*zdku + zdkv*zdkv, zeps ) … … 691 696 zts(jp_tem) = ttbl(ji,jj) 692 697 zts(jp_sal) = stbl(ji,jj) 693 zdep = gdepw _n(ji,jj,ikt)698 zdep = gdepw(ji,jj,ikt,Kmm) 694 699 ! 695 CALL eos_rab( zts, zdep, zab )700 CALL eos_rab( zts, zdep, zab, Kmm ) 696 701 ! 697 702 !! compute length scale … … 700 705 !! compute Monin Obukov Length 701 706 ! Maximum boundary layer depth 702 zhmax = gdept _n(ji,jj,mbkt(ji,jj)) - gdepw_n(ji,jj,mikt(ji,jj)) -0.001_wp707 zhmax = gdept(ji,jj,mbkt(ji,jj),Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) -0.001_wp 703 708 ! Compute Monin obukhov length scale at the surface and Ekman depth: 704 709 zmob = zustar(ji,jj) ** 3 / (vkarmn * (zbuofdep + zeps)) … … 727 732 728 733 729 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin )734 SUBROUTINE sbc_isf_tbl( pvarin, pvarout, cd_ptin, Kmm ) 730 735 !!---------------------------------------------------------------------- 731 736 !! *** SUBROUTINE sbc_isf_tbl *** … … 737 742 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pvarout 738 743 CHARACTER(len=1), INTENT(in ) :: cd_ptin ! point of variable in/out 744 INTEGER , INTENT(in ) :: Kmm ! ocean time level indices 739 745 ! 740 746 INTEGER :: ji, jj, jk ! loop index … … 753 759 ikt = miku(ji,jj) ; ikb = miku(ji,jj) 754 760 ! thickness of boundary layer at least the top level thickness 755 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u _n(ji,jj,ikt) )761 zhisf_tbl(ji,jj) = MAX( rhisf_tbl_0(ji,jj) , e3u(ji,jj,ikt,Kmm) ) 756 762 757 763 ! determine the deepest level influenced by the boundary layer 758 764 DO jk = ikt+1, mbku(ji,jj) 759 IF ( (SUM(e3u _n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk760 END DO 761 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u _n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.765 IF ( (SUM(e3u(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 766 END DO 767 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness. 762 768 763 769 ! level fully include in the ice shelf boundary layer 764 770 DO jk = ikt, ikb - 1 765 ze3 = e3u _n(ji,jj,jk)771 ze3 = e3u(ji,jj,jk,Kmm) 766 772 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 767 773 END DO 768 774 769 775 ! level partially include in ice shelf boundary layer 770 zhk = SUM( e3u _n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)776 zhk = SUM( e3u(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj) 771 777 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 772 778 END DO … … 785 791 ikt = mikv(ji,jj) ; ikb = mikv(ji,jj) 786 792 ! thickness of boundary layer at least the top level thickness 787 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v _n(ji,jj,ikt))793 zhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3v(ji,jj,ikt,Kmm)) 788 794 789 795 ! determine the deepest level influenced by the boundary layer 790 796 DO jk = ikt+1, mbkv(ji,jj) 791 IF ( (SUM(e3v _n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk792 END DO 793 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v _n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.797 IF ( (SUM(e3v(ji,jj,ikt:jk-1,Kmm)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 798 END DO 799 zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness. 794 800 795 801 ! level fully include in the ice shelf boundary layer 796 802 DO jk = ikt, ikb - 1 797 ze3 = e3v _n(ji,jj,jk)803 ze3 = e3v(ji,jj,jk,Kmm) 798 804 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) / zhisf_tbl(ji,jj) * ze3 799 805 END DO 800 806 801 807 ! level partially include in ice shelf boundary layer 802 zhk = SUM( e3v _n(ji, jj, ikt:ikb - 1)) / zhisf_tbl(ji,jj)808 zhk = SUM( e3v(ji, jj, ikt:ikb - 1,Kmm)) / zhisf_tbl(ji,jj) 803 809 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 804 810 END DO … … 820 826 ! level fully include in the ice shelf boundary layer 821 827 DO jk = ikt, ikb - 1 822 ze3 = e3t _n(ji,jj,jk)828 ze3 = e3t(ji,jj,jk,Kmm) 823 829 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,jk) * r1_hisf_tbl(ji,jj) * ze3 824 830 END DO 825 831 826 832 ! level partially include in ice shelf boundary layer 827 zhk = SUM( e3t _n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj)833 zhk = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) 828 834 pvarout(ji,jj) = pvarout(ji,jj) + pvarin(ji,jj,ikb) * (1._wp - zhk) 829 835 END DO … … 837 843 838 844 839 SUBROUTINE sbc_isf_div( phdivn )845 SUBROUTINE sbc_isf_div( phdivn, Kmm ) 840 846 !!---------------------------------------------------------------------- 841 847 !! *** SUBROUTINE sbc_isf_div *** … … 850 856 !!---------------------------------------------------------------------- 851 857 REAL(wp), DIMENSION(:,:,:), INTENT( inout ) :: phdivn ! horizontal divergence 858 INTEGER , INTENT( in ) :: Kmm ! ocean time level indices 852 859 ! 853 860 INTEGER :: ji, jj, jk ! dummy loop indices … … 865 872 ikb = misfkt(ji,jj) 866 873 ! thickness of boundary layer at least the top level thickness 867 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t _n(ji,jj,ikt))874 rhisf_tbl(ji,jj) = MAX(rhisf_tbl_0(ji,jj), e3t(ji,jj,ikt,Kmm)) 868 875 869 876 ! determine the deepest level influenced by the boundary layer 870 877 DO jk = ikt, mbkt(ji,jj) 871 IF ( (SUM(e3t _n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk872 END DO 873 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t _n(ji,jj,ikt:ikb))) ! limit the tbl to water thickness.878 IF ( (SUM(e3t(ji,jj,ikt:jk-1,Kmm)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 879 END DO 880 rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t(ji,jj,ikt:ikb,Kmm))) ! limit the tbl to water thickness. 874 881 misfkb(ji,jj) = ikb ! last wet level of the tbl 875 882 r1_hisf_tbl(ji,jj) = 1._wp / rhisf_tbl(ji,jj) 876 883 877 zhk = SUM( e3t _n(ji, jj, ikt:ikb - 1)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1878 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t _n(ji,jj,ikb) ! proportion of bottom cell influenced by boundary layer884 zhk = SUM( e3t(ji, jj, ikt:ikb - 1,Kmm)) * r1_hisf_tbl(ji,jj) ! proportion of tbl cover by cell from ikt to ikb - 1 885 ralpha(ji,jj) = rhisf_tbl(ji,jj) * (1._wp - zhk ) / e3t(ji,jj,ikb,Kmm) ! proportion of bottom cell influenced by boundary layer 879 886 END DO 880 887 END DO
Note: See TracChangeset
for help on using the changeset viewer.