Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icevar.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icevar.F90
r10589 r11822 32 32 !! - vt_s(jpi,jpj) 33 33 !! - at_i(jpi,jpj) 34 !! - st_i(jpi,jpj) 34 35 !! - et_s(jpi,jpj) total snow heat content 35 36 !! - et_i(jpi,jpj) total ice thermal content … … 44 45 !! ice_var_salprof1d : salinity profile in the ice 1D 45 46 !! ice_var_zapsmall : remove very small area and volume 46 !! ice_var_zapneg : remove negative ice fields (to debug the advection scheme UM3-5) 47 !! ice_var_itd : convert 1-cat to jpl-cat 48 !! ice_var_itd2 : convert N-cat to jpl-cat 47 !! ice_var_zapneg : remove negative ice fields 48 !! ice_var_roundoff : remove negative values arising from roundoff erros 49 49 !! ice_var_bv : brine volume 50 50 !! ice_var_enthalpy : compute ice and snow enthalpies from temperature 51 51 !! ice_var_sshdyn : compute equivalent ssh in lead 52 !! ice_var_itd : convert N-cat to M-cat 52 53 !!---------------------------------------------------------------------- 53 54 USE dom_oce ! ocean space and time domain … … 71 72 PUBLIC ice_var_zapsmall 72 73 PUBLIC ice_var_zapneg 73 PUBLIC ice_var_itd 74 PUBLIC ice_var_itd2 74 PUBLIC ice_var_roundoff 75 75 PUBLIC ice_var_bv 76 76 PUBLIC ice_var_enthalpy 77 77 PUBLIC ice_var_sshdyn 78 PUBLIC ice_var_itd 79 80 INTERFACE ice_var_itd 81 MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc 82 END INTERFACE 78 83 79 84 !!---------------------------------------------------------------------- … … 99 104 ! 100 105 ! ! integrated values 101 vt_i(:,:) = SUM( v_i(:,:,:) , dim=3 ) 102 vt_s(:,:) = SUM( v_s(:,:,:) , dim=3 ) 103 at_i(:,:) = SUM( a_i(:,:,:) , dim=3 ) 104 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 105 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 106 vt_i(:,:) = SUM( v_i (:,:,:) , dim=3 ) 107 vt_s(:,:) = SUM( v_s (:,:,:) , dim=3 ) 108 st_i(:,:) = SUM( sv_i(:,:,:) , dim=3 ) 109 at_i(:,:) = SUM( a_i (:,:,:) , dim=3 ) 110 et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 111 et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 106 112 ! 107 113 at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds … … 133 139 tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 134 140 om_i (:,:) = SUM( oa_i(:,:,:) , dim=3 ) * z1_at_i(:,:) 135 sm_i (:,:) = SUM( sv_i(:,:,:) , dim=3 )* z1_vt_i(:,:)141 sm_i (:,:) = st_i(:,:) * z1_vt_i(:,:) 136 142 ! 137 143 tm_i(:,:) = 0._wp … … 153 159 tm_s (:,:) = rt0 154 160 END WHERE 155 161 ! 162 ! ! mean melt pond depth 163 WHERE( at_ip(:,:) > epsi20 ) ; hm_ip(:,:) = vt_ip(:,:) / at_ip(:,:) 164 ELSEWHERE ; hm_ip(:,:) = 0._wp 165 END WHERE 166 ! 156 167 DEALLOCATE( z1_at_i , z1_vt_i , z1_vt_s ) 168 ! 157 169 ENDIF 158 170 ! … … 229 241 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 230 242 ! 231 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i 243 ze_i = e_i (ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 232 244 ztmelts = - sz_i(ji,jj,jk,jl) * rTmlt ! Ice layer melt temperature [C] 233 245 ! Conversion q(S,T) -> T (second order equation) … … 236 248 t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * 0.5_wp * r1_rcpi , ztmelts ) ) + rt0 ! [K] with bounds: -100 < t_i < ztmelts 237 249 ! 238 ELSE !--- no ice250 ELSE !--- no ice 239 251 t_i(ji,jj,jk,jl) = rt0 240 252 ENDIF … … 258 270 ! 259 271 ! integrated values 260 vt_i (:,:) = SUM( v_i , dim=3 )261 vt_s (:,:) = SUM( v_s , dim=3 )262 at_i (:,:) = SUM( a_i , dim=3 )272 vt_i (:,:) = SUM( v_i , dim=3 ) 273 vt_s (:,:) = SUM( v_s , dim=3 ) 274 at_i (:,:) = SUM( a_i , dim=3 ) 263 275 ! 264 276 END SUBROUTINE ice_var_glo2eqv … … 528 540 529 541 ! to be sure that at_i is the sum of a_i(jl) 530 at_i (:,:) = SUM( a_i(:,:,:), dim=3 ) 531 vt_i (:,:) = SUM( v_i(:,:,:), dim=3 ) 542 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) 543 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) 544 !!clem add? 545 ! vt_s (:,:) = SUM( v_s (:,:,:), dim=3 ) 546 ! st_i (:,:) = SUM( sv_i(:,:,:), dim=3 ) 547 ! et_s(:,:) = SUM( SUM( e_s (:,:,:,:), dim=4 ), dim=3 ) 548 ! et_i(:,:) = SUM( SUM( e_i (:,:,:,:), dim=4 ), dim=3 ) 549 !!clem 532 550 533 551 ! open water = 1 if at_i=0 … … 537 555 538 556 539 SUBROUTINE ice_var_zapneg( p ato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )557 SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 540 558 !!------------------------------------------------------------------- 541 559 !! *** ROUTINE ice_var_zapneg *** … … 543 561 !! ** Purpose : Remove negative sea ice fields and correct fluxes 544 562 !!------------------------------------------------------------------- 545 INTEGER :: ji, jj, jl, jk ! dummy loop indices 546 ! 563 REAL(wp) , INTENT(in ) :: pdt ! tracer time-step 547 564 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pato_i ! open water area 548 565 REAL(wp), DIMENSION(:,:,:) , INTENT(inout) :: pv_i ! ice volume … … 555 572 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_s ! snw heat content 556 573 REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) :: pe_i ! ice heat content 557 !!------------------------------------------------------------------- 558 ! 574 ! 575 INTEGER :: ji, jj, jl, jk ! dummy loop indices 576 REAL(wp) :: z1_dt 577 !!------------------------------------------------------------------- 578 ! 579 z1_dt = 1._wp / pdt 559 580 ! 560 581 DO jl = 1, jpl !== loop over the categories ==! 561 582 ! 583 ! make sure a_i=0 where v_i<=0 584 WHERE( pv_i(:,:,:) <= 0._wp ) pa_i(:,:,:) = 0._wp 585 562 586 !---------------------------------------- 563 587 ! zap ice energy and send it to the ocean … … 566 590 DO jj = 1 , jpj 567 591 DO ji = 1 , jpi 568 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN569 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * r1_rdtice! W.m-2 >0592 IF( pe_i(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 593 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * z1_dt ! W.m-2 >0 570 594 pe_i(ji,jj,jk,jl) = 0._wp 571 595 ENDIF … … 577 601 DO jj = 1 , jpj 578 602 DO ji = 1 , jpi 579 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN580 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * r1_rdtice! W.m-2 <0603 IF( pe_s(ji,jj,jk,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 604 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * z1_dt ! W.m-2 <0 581 605 pe_s(ji,jj,jk,jl) = 0._wp 582 606 ENDIF … … 590 614 DO jj = 1 , jpj 591 615 DO ji = 1 , jpi 592 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN593 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * r1_rdtice616 IF( pv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 617 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_i (ji,jj,jl) * rhoi * z1_dt 594 618 pv_i (ji,jj,jl) = 0._wp 595 619 ENDIF 596 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN597 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * r1_rdtice620 IF( pv_s(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp ) THEN 621 wfx_res(ji,jj) = wfx_res(ji,jj) + pv_s (ji,jj,jl) * rhos * z1_dt 598 622 pv_s (ji,jj,jl) = 0._wp 599 623 ENDIF 600 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) < 0._wp ) THEN601 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * r1_rdtice624 IF( psv_i(ji,jj,jl) < 0._wp .OR. pa_i(ji,jj,jl) <= 0._wp .OR. pv_i(ji,jj,jl) <= 0._wp ) THEN 625 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * rhoi * z1_dt 602 626 psv_i (ji,jj,jl) = 0._wp 603 627 ENDIF … … 616 640 END SUBROUTINE ice_var_zapneg 617 641 642 643 SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i ) 644 !!------------------------------------------------------------------- 645 !! *** ROUTINE ice_var_roundoff *** 646 !! 647 !! ** Purpose : Remove negative sea ice values arising from roundoff errors 648 !!------------------------------------------------------------------- 649 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_i ! ice concentration 650 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_i ! ice volume 651 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_s ! snw volume 652 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: psv_i ! salt content 653 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: poa_i ! age content 654 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pa_ip ! melt pond fraction 655 REAL(wp), DIMENSION(:,:) , INTENT(inout) :: pv_ip ! melt pond volume 656 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_s ! snw heat content 657 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: pe_i ! ice heat content 658 !!------------------------------------------------------------------- 659 ! 660 WHERE( pa_i (1:npti,:) < 0._wp .AND. pa_i (1:npti,:) > -epsi10 ) pa_i (1:npti,:) = 0._wp ! a_i must be >= 0 661 WHERE( pv_i (1:npti,:) < 0._wp .AND. pv_i (1:npti,:) > -epsi10 ) pv_i (1:npti,:) = 0._wp ! v_i must be >= 0 662 WHERE( pv_s (1:npti,:) < 0._wp .AND. pv_s (1:npti,:) > -epsi10 ) pv_s (1:npti,:) = 0._wp ! v_s must be >= 0 663 WHERE( psv_i(1:npti,:) < 0._wp .AND. psv_i(1:npti,:) > -epsi10 ) psv_i(1:npti,:) = 0._wp ! sv_i must be >= 0 664 WHERE( poa_i(1:npti,:) < 0._wp .AND. poa_i(1:npti,:) > -epsi10 ) poa_i(1:npti,:) = 0._wp ! oa_i must be >= 0 665 WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 666 WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 667 IF( ln_pnd_H12 ) THEN 668 WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 669 WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 670 ENDIF 671 ! 672 END SUBROUTINE ice_var_roundoff 618 673 619 SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )620 !!-------------------------------------------------------------------621 !! *** ROUTINE ice_var_itd ***622 !!623 !! ** Purpose : converting 1-cat ice to multiple ice categories624 !!625 !! ice thickness distribution follows a gaussian law626 !! around the concentration of the most likely ice thickness627 !! (similar as iceistate.F90)628 !!629 !! ** Method: Iterative procedure630 !!631 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian632 !!633 !! 2) Check whether the distribution conserves area and volume, positivity and634 !! category boundaries635 !!636 !! 3) If not (input ice is too thin), the last category is empty and637 !! the number of categories is reduced (jpl-1)638 !!639 !! 4) Iterate until ok (SUM(itest(:) = 4)640 !!641 !! ** Arguments : zhti: 1-cat ice thickness642 !! zhts: 1-cat snow depth643 !! zati: 1-cat ice concentration644 !!645 !! ** Output : jpl-cat646 !!647 !! (Example of application: BDY forcings when input are cell averaged)648 !!-------------------------------------------------------------------649 INTEGER :: ji, jk, jl ! dummy loop indices650 INTEGER :: idim, i_fill, jl0651 REAL(wp) :: zarg, zV, zconv, zdh, zdv652 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables653 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables654 INTEGER , DIMENSION(4) :: itest655 !!-------------------------------------------------------------------656 !657 ! ----------------------------------------658 ! distribution over the jpl ice categories659 ! ----------------------------------------660 ! a gaussian distribution for ice concentration is used661 ! then we check whether the distribution fullfills662 ! volume and area conservation, positivity and ice categories bounds663 idim = SIZE( zhti , 1 )664 zh_i(1:idim,1:jpl) = 0._wp665 zh_s(1:idim,1:jpl) = 0._wp666 za_i(1:idim,1:jpl) = 0._wp667 !668 DO ji = 1, idim669 !670 IF( zhti(ji) > 0._wp ) THEN671 !672 ! find which category (jl0) the input ice thickness falls into673 jl0 = jpl674 DO jl = 1, jpl675 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN676 jl0 = jl677 CYCLE678 ENDIF679 END DO680 !681 itest(:) = 0682 i_fill = jpl + 1 !------------------------------------683 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories684 ! !------------------------------------685 i_fill = i_fill - 1686 !687 zh_i(ji,1:jpl) = 0._wp688 za_i(ji,1:jpl) = 0._wp689 itest(:) = 0690 !691 IF ( i_fill == 1 ) THEN !-- case very thin ice: fill only category 1692 zh_i(ji,1) = zhti(ji)693 za_i (ji,1) = zati (ji)694 ELSE !-- case ice is thicker: fill categories >1695 ! thickness696 DO jl = 1, i_fill - 1697 zh_i(ji,jl) = hi_mean(jl)698 END DO699 !700 ! concentration701 za_i(ji,jl0) = zati(ji) / SQRT(REAL(jpl))702 DO jl = 1, i_fill - 1703 IF ( jl /= jl0 ) THEN704 zarg = ( zh_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )705 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2)706 ENDIF707 END DO708 !709 ! last category710 za_i(ji,i_fill) = zati(ji) - SUM( za_i(ji,1:i_fill-1) )711 zV = SUM( za_i(ji,1:i_fill-1) * zh_i(ji,1:i_fill-1) )712 zh_i(ji,i_fill) = ( zhti(ji) * zati(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 )713 !714 ! correction if concentration of upper cat is greater than lower cat715 ! (it should be a gaussian around jl0 but sometimes it is not)716 IF ( jl0 /= jpl ) THEN717 DO jl = jpl, jl0+1, -1718 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN719 zdv = zh_i(ji,jl) * za_i(ji,jl)720 zh_i(ji,jl ) = 0._wp721 za_i (ji,jl ) = 0._wp722 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )723 END IF724 END DO725 ENDIF726 !727 ENDIF728 !729 ! Compatibility tests730 zconv = ABS( zati(ji) - SUM( za_i(ji,1:jpl) ) )731 IF ( zconv < epsi06 ) itest(1) = 1 ! Test 1: area conservation732 !733 zconv = ABS( zhti(ji)*zati(ji) - SUM( za_i(ji,1:jpl)*zh_i(ji,1:jpl) ) )734 IF ( zconv < epsi06 ) itest(2) = 1 ! Test 2: volume conservation735 !736 IF ( zh_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 ! Test 3: thickness of the last category is in-bounds ?737 !738 itest(4) = 1739 DO jl = 1, i_fill740 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 ! Test 4: positivity of ice concentrations741 END DO742 ! !----------------------------743 END DO ! end iteration on categories744 ! !----------------------------745 ENDIF746 END DO747 748 ! Add Snow in each category where za_i is not 0749 DO jl = 1, jpl750 DO ji = 1, idim751 IF( za_i(ji,jl) > 0._wp ) THEN752 zh_s(ji,jl) = zh_i(ji,jl) * ( zhts(ji) / zhti(ji) )753 ! In case snow load is in excess that would lead to transformation from snow to ice754 ! Then, transfer the snow excess into the ice (different from icethd_dh)755 zdh = MAX( 0._wp, ( rhos * zh_s(ji,jl) + ( rhoi - rau0 ) * zh_i(ji,jl) ) * r1_rau0 )756 ! recompute h_i, h_s avoiding out of bounds values757 zh_i(ji,jl) = MIN( hi_max(jl), zh_i(ji,jl) + zdh )758 zh_s(ji,jl) = MAX( 0._wp, zh_s(ji,jl) - zdh * rhoi * r1_rhos )759 ENDIF760 END DO761 END DO762 !763 END SUBROUTINE ice_var_itd764 765 766 SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )767 !!-------------------------------------------------------------------768 !! *** ROUTINE ice_var_itd2 ***769 !!770 !! ** Purpose : converting N-cat ice to jpl ice categories771 !!772 !! ice thickness distribution follows a gaussian law773 !! around the concentration of the most likely ice thickness774 !! (similar as iceistate.F90)775 !!776 !! ** Method: Iterative procedure777 !!778 !! 1) Fill ice cat that correspond to input thicknesses779 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled780 !!781 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1782 !! by removing 25% ice area from jlmin and jlmax (resp.)783 !!784 !! 3) Expand the filling to the empty cat between jlmin and jlmax785 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax)786 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin)787 !!788 !! ** Arguments : zhti: N-cat ice thickness789 !! zhts: N-cat snow depth790 !! zati: N-cat ice concentration791 !!792 !! ** Output : jpl-cat793 !!794 !! (Example of application: BDY forcings when inputs have N-cat /= jpl)795 !!-------------------------------------------------------------------796 INTEGER :: ji, jl, jl1, jl2 ! dummy loop indices797 INTEGER :: idim, icat798 REAL(wp), PARAMETER :: ztrans = 0.25_wp799 REAL(wp), DIMENSION(:,:), INTENT(in) :: zhti, zhts, zati ! input ice/snow variables800 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zh_i, zh_s, za_i ! output ice/snow variables801 INTEGER , DIMENSION(:,:), ALLOCATABLE :: jlfil, jlfil2802 INTEGER , DIMENSION(:) , ALLOCATABLE :: jlmax, jlmin803 !!-------------------------------------------------------------------804 !805 idim = SIZE( zhti, 1 )806 icat = SIZE( zhti, 2 )807 !808 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays809 ALLOCATE( jlmin(idim), jlmax(idim) )810 811 ! --- initialize output fields to 0 --- !812 zh_i(1:idim,1:jpl) = 0._wp813 zh_s(1:idim,1:jpl) = 0._wp814 za_i(1:idim,1:jpl) = 0._wp815 !816 ! --- fill the categories --- !817 ! find where cat-input = cat-output and fill cat-output fields818 jlmax(:) = 0819 jlmin(:) = 999820 jlfil(:,:) = 0821 DO jl1 = 1, jpl822 DO jl2 = 1, icat823 DO ji = 1, idim824 IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN825 ! fill the right category826 zh_i(ji,jl1) = zhti(ji,jl2)827 zh_s(ji,jl1) = zhts(ji,jl2)828 za_i(ji,jl1) = zati(ji,jl2)829 ! record categories that are filled830 jlmax(ji) = MAX( jlmax(ji), jl1 )831 jlmin(ji) = MIN( jlmin(ji), jl1 )832 jlfil(ji,jl1) = jl1833 ENDIF834 END DO835 END DO836 END DO837 !838 ! --- fill the gaps between categories --- !839 ! transfer from categories filled at the previous step to the empty ones in between840 DO ji = 1, idim841 jl1 = jlmin(ji)842 jl2 = jlmax(ji)843 IF( jl1 > 1 ) THEN844 ! fill the lower cat (jl1-1)845 za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)846 zh_i(ji,jl1-1) = hi_mean(jl1-1)847 ! remove from cat jl1848 za_i(ji,jl1 ) = ( 1._wp - ztrans ) * za_i(ji,jl1)849 ENDIF850 IF( jl2 < jpl ) THEN851 ! fill the upper cat (jl2+1)852 za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)853 zh_i(ji,jl2+1) = hi_mean(jl2+1)854 ! remove from cat jl2855 za_i(ji,jl2 ) = ( 1._wp - ztrans ) * za_i(ji,jl2)856 ENDIF857 END DO858 !859 jlfil2(:,:) = jlfil(:,:)860 ! fill categories from low to high861 DO jl = 2, jpl-1862 DO ji = 1, idim863 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN864 ! fill high865 za_i(ji,jl) = ztrans * za_i(ji,jl-1)866 zh_i(ji,jl) = hi_mean(jl)867 jlfil(ji,jl) = jl868 ! remove low869 za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)870 ENDIF871 END DO872 END DO873 !874 ! fill categories from high to low875 DO jl = jpl-1, 2, -1876 DO ji = 1, idim877 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN878 ! fill low879 za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)880 zh_i(ji,jl) = hi_mean(jl)881 jlfil2(ji,jl) = jl882 ! remove high883 za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)884 ENDIF885 END DO886 END DO887 !888 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays889 DEALLOCATE( jlmin, jlmax )890 !891 END SUBROUTINE ice_var_itd2892 893 674 894 675 SUBROUTINE ice_var_bv … … 952 733 END SUBROUTINE ice_var_enthalpy 953 734 735 954 736 FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b) 955 737 !!--------------------------------------------------------------------- … … 998 780 END FUNCTION ice_var_sshdyn 999 781 782 783 !!------------------------------------------------------------------- 784 !! *** INTERFACE ice_var_itd *** 785 !! 786 !! ** Purpose : converting N-cat ice to jpl ice categories 787 !!------------------------------------------------------------------- 788 SUBROUTINE ice_var_itd_1c1c( phti, phts, pati , ph_i, ph_s, pa_i, & 789 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 790 !!------------------------------------------------------------------- 791 !! ** Purpose : converting 1-cat ice to 1 ice category 792 !!------------------------------------------------------------------- 793 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 794 REAL(wp), DIMENSION(:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 795 REAL(wp), DIMENSION(:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 796 REAL(wp), DIMENSION(:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 797 !!------------------------------------------------------------------- 798 ! == thickness and concentration == ! 799 ph_i(:) = phti(:) 800 ph_s(:) = phts(:) 801 pa_i(:) = pati(:) 802 ! 803 ! == temperature and salinity and ponds == ! 804 pt_i (:) = ptmi (:) 805 pt_s (:) = ptms (:) 806 pt_su(:) = ptmsu(:) 807 ps_i (:) = psmi (:) 808 pa_ip(:) = patip(:) 809 ph_ip(:) = phtip(:) 810 811 END SUBROUTINE ice_var_itd_1c1c 812 813 SUBROUTINE ice_var_itd_Nc1c( phti, phts, pati , ph_i, ph_s, pa_i, & 814 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 815 !!------------------------------------------------------------------- 816 !! ** Purpose : converting N-cat ice to 1 ice category 817 !!------------------------------------------------------------------- 818 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 819 REAL(wp), DIMENSION(:) , INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 820 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 821 REAL(wp), DIMENSION(:) , INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 822 ! 823 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs 824 ! 825 INTEGER :: idim 826 !!------------------------------------------------------------------- 827 ! 828 idim = SIZE( phti, 1 ) 829 ! 830 ! == thickness and concentration == ! 831 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim) ) 832 ! 833 pa_i(:) = SUM( pati(:,:), dim=2 ) 834 835 WHERE( ( pa_i(:) ) /= 0._wp ) ; z1_ai(:) = 1._wp / pa_i(:) 836 ELSEWHERE ; z1_ai(:) = 0._wp 837 END WHERE 838 839 ph_i(:) = SUM( phti(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 840 ph_s(:) = SUM( phts(:,:) * pati(:,:), dim=2 ) * z1_ai(:) 841 ! 842 ! == temperature and salinity == ! 843 WHERE( ( pa_i(:) * ph_i(:) ) /= 0._wp ) ; z1_vi(:) = 1._wp / ( pa_i(:) * ph_i(:) ) 844 ELSEWHERE ; z1_vi(:) = 0._wp 845 END WHERE 846 WHERE( ( pa_i(:) * ph_s(:) ) /= 0._wp ) ; z1_vs(:) = 1._wp / ( pa_i(:) * ph_s(:) ) 847 ELSEWHERE ; z1_vs(:) = 0._wp 848 END WHERE 849 pt_i (:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 850 pt_s (:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 851 pt_su(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 852 ps_i (:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 853 854 ! == ponds == ! 855 pa_ip(:) = SUM( patip(:,:), dim=2 ) 856 WHERE( pa_ip(:) /= 0._wp ) ; ph_ip(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / pa_ip(:) 857 ELSEWHERE ; ph_ip(:) = 0._wp 858 END WHERE 859 ! 860 DEALLOCATE( z1_ai, z1_vi, z1_vs ) 861 ! 862 END SUBROUTINE ice_var_itd_Nc1c 863 864 SUBROUTINE ice_var_itd_1cMc( phti, phts, pati , ph_i, ph_s, pa_i, & 865 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 866 !!------------------------------------------------------------------- 867 !! 868 !! ** Purpose : converting 1-cat ice to jpl ice categories 869 !! 870 !! 871 !! ** Method: ice thickness distribution follows a gamma function from Abraham et al. (2015) 872 !! it has the property of conserving total concentration and volume 873 !! 874 !! 875 !! ** Arguments : phti: 1-cat ice thickness 876 !! phts: 1-cat snow depth 877 !! pati: 1-cat ice concentration 878 !! 879 !! ** Output : jpl-cat 880 !! 881 !! Abraham, C., Steiner, N., Monahan, A. and Michel, C., 2015. 882 !! Effects of subgrid‐scale snow thickness variability on radiative transfer in sea ice. 883 !! Journal of Geophysical Research: Oceans, 120(8), pp.5597-5614 884 !!------------------------------------------------------------------- 885 REAL(wp), DIMENSION(:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 886 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 887 REAL(wp), DIMENSION(:) , INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 888 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 889 ! 890 REAL(wp), ALLOCATABLE, DIMENSION(:) :: zfra, z1_hti 891 INTEGER :: ji, jk, jl 892 INTEGER :: idim 893 REAL(wp) :: zv, zdh 894 !!------------------------------------------------------------------- 895 ! 896 idim = SIZE( phti , 1 ) 897 ! 898 ph_i(1:idim,1:jpl) = 0._wp 899 ph_s(1:idim,1:jpl) = 0._wp 900 pa_i(1:idim,1:jpl) = 0._wp 901 ! 902 ALLOCATE( z1_hti(idim) ) 903 WHERE( phti(:) /= 0._wp ) ; z1_hti(:) = 1._wp / phti(:) 904 ELSEWHERE ; z1_hti(:) = 0._wp 905 END WHERE 906 ! 907 ! == thickness and concentration == ! 908 ! for categories 1:jpl-1, integrate the gamma function from hi_max(jl-1) to hi_max(jl) 909 DO jl = 1, jpl-1 910 DO ji = 1, idim 911 ! 912 IF( phti(ji) > 0._wp ) THEN 913 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 914 pa_i(ji,jl) = pati(ji) * z1_hti(ji) * ( ( phti(ji) + 2.*hi_max(jl-1) ) * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 915 & - ( phti(ji) + 2.*hi_max(jl ) ) * EXP( -2.*hi_max(jl )*z1_hti(ji) ) ) 916 ! 917 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jl-1) to hi_max(jl) 918 zv = pati(ji) * z1_hti(ji) * ( ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl-1) + 2.*hi_max(jl-1)*hi_max(jl-1) ) & 919 & * EXP( -2.*hi_max(jl-1)*z1_hti(ji) ) & 920 & - ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jl) + 2.*hi_max(jl)*hi_max(jl) ) & 921 & * EXP(-2.*hi_max(jl)*z1_hti(ji)) ) 922 ! thickness 923 IF( pa_i(ji,jl) > epsi06 ) THEN 924 ph_i(ji,jl) = zv / pa_i(ji,jl) 925 ELSE 926 ph_i(ji,jl) = 0. 927 pa_i(ji,jl) = 0. 928 ENDIF 929 ENDIF 930 ! 931 ENDDO 932 ENDDO 933 ! 934 ! for the last category (jpl), integrate the gamma function from hi_max(jpl-1) to infinity 935 DO ji = 1, idim 936 ! 937 IF( phti(ji) > 0._wp ) THEN 938 ! concentration : integrate ((4A/H^2)xexp(-2x/H))dx from x=hi_max(jpl-1) to infinity 939 pa_i(ji,jpl) = pati(ji) * z1_hti(ji) * ( phti(ji) + 2.*hi_max(jpl-1) ) * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 940 941 ! volume : integrate ((4A/H^2)x^2exp(-2x/H))dx from x=hi_max(jpl-1) to infinity 942 zv = pati(ji) * z1_hti(ji) * ( phti(ji)*phti(ji) + 2.*phti(ji)*hi_max(jpl-1) + 2.*hi_max(jpl-1)*hi_max(jpl-1) ) & 943 & * EXP( -2.*hi_max(jpl-1)*z1_hti(ji) ) 944 ! thickness 945 IF( pa_i(ji,jpl) > epsi06 ) THEN 946 ph_i(ji,jpl) = zv / pa_i(ji,jpl) 947 else 948 ph_i(ji,jpl) = 0. 949 pa_i(ji,jpl) = 0. 950 ENDIF 951 ENDIF 952 ! 953 ENDDO 954 ! 955 ! Add Snow in each category where pa_i is not 0 956 DO jl = 1, jpl 957 DO ji = 1, idim 958 IF( pa_i(ji,jl) > 0._wp ) THEN 959 ph_s(ji,jl) = ph_i(ji,jl) * phts(ji) * z1_hti(ji) 960 ! In case snow load is in excess that would lead to transformation from snow to ice 961 ! Then, transfer the snow excess into the ice (different from icethd_dh) 962 zdh = MAX( 0._wp, ( rhos * ph_s(ji,jl) + ( rhoi - rau0 ) * ph_i(ji,jl) ) * r1_rau0 ) 963 ! recompute h_i, h_s avoiding out of bounds values 964 ph_i(ji,jl) = MIN( hi_max(jl), ph_i(ji,jl) + zdh ) 965 ph_s(ji,jl) = MAX( 0._wp, ph_s(ji,jl) - zdh * rhoi * r1_rhos ) 966 ENDIF 967 END DO 968 END DO 969 ! 970 DEALLOCATE( z1_hti ) 971 ! 972 ! == temperature and salinity == ! 973 DO jl = 1, jpl 974 pt_i (:,jl) = ptmi (:) 975 pt_s (:,jl) = ptms (:) 976 pt_su(:,jl) = ptmsu(:) 977 ps_i (:,jl) = psmi (:) 978 ps_i (:,jl) = psmi (:) 979 END DO 980 ! 981 ! == ponds == ! 982 ALLOCATE( zfra(idim) ) 983 ! keep the same pond fraction atip/ati for each category 984 WHERE( pati(:) /= 0._wp ) ; zfra(:) = patip(:) / pati(:) 985 ELSEWHERE ; zfra(:) = 0._wp 986 END WHERE 987 DO jl = 1, jpl 988 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 989 END DO 990 ! keep the same v_ip/v_i ratio for each category 991 WHERE( ( phti(:) * pati(:) ) /= 0._wp ) ; zfra(:) = ( phtip(:) * patip(:) ) / ( phti(:) * pati(:) ) 992 ELSEWHERE ; zfra(:) = 0._wp 993 END WHERE 994 DO jl = 1, jpl 995 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 996 ELSEWHERE ; ph_ip(:,jl) = 0._wp 997 END WHERE 998 END DO 999 DEALLOCATE( zfra ) 1000 ! 1001 END SUBROUTINE ice_var_itd_1cMc 1002 1003 SUBROUTINE ice_var_itd_NcMc( phti, phts, pati , ph_i, ph_s, pa_i, & 1004 & ptmi, ptms, ptmsu, psmi, patip, phtip, pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ) 1005 !!------------------------------------------------------------------- 1006 !! 1007 !! ** Purpose : converting N-cat ice to jpl ice categories 1008 !! 1009 !! ice thickness distribution follows a gaussian law 1010 !! around the concentration of the most likely ice thickness 1011 !! (similar as iceistate.F90) 1012 !! 1013 !! ** Method: Iterative procedure 1014 !! 1015 !! 1) Fill ice cat that correspond to input thicknesses 1016 !! Find the lowest(jlmin) and highest(jlmax) cat that are filled 1017 !! 1018 !! 2) Expand the filling to the cat jlmin-1 and jlmax+1 1019 !! by removing 25% ice area from jlmin and jlmax (resp.) 1020 !! 1021 !! 3) Expand the filling to the empty cat between jlmin and jlmax 1022 !! by a) removing 25% ice area from the lower cat (ascendant loop jlmin=>jlmax) 1023 !! b) removing 25% ice area from the higher cat (descendant loop jlmax=>jlmin) 1024 !! 1025 !! ** Arguments : phti: N-cat ice thickness 1026 !! phts: N-cat snow depth 1027 !! pati: N-cat ice concentration 1028 !! 1029 !! ** Output : jpl-cat 1030 !! 1031 !! (Example of application: BDY forcings when inputs have N-cat /= jpl) 1032 !!------------------------------------------------------------------- 1033 REAL(wp), DIMENSION(:,:), INTENT(in) :: phti, phts, pati ! input ice/snow variables 1034 REAL(wp), DIMENSION(:,:), INTENT(inout) :: ph_i, ph_s, pa_i ! output ice/snow variables 1035 REAL(wp), DIMENSION(:,:), INTENT(in) :: ptmi, ptms, ptmsu, psmi, patip, phtip ! input ice/snow temp & sal & ponds 1036 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pt_i, pt_s, pt_su, ps_i, pa_ip, ph_ip ! output ice/snow temp & sal & ponds 1037 ! 1038 INTEGER , ALLOCATABLE, DIMENSION(:,:) :: jlfil, jlfil2 1039 INTEGER , ALLOCATABLE, DIMENSION(:) :: jlmax, jlmin 1040 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z1_ai, z1_vi, z1_vs, ztmp, zfra 1041 ! 1042 REAL(wp), PARAMETER :: ztrans = 0.25_wp 1043 INTEGER :: ji, jl, jl1, jl2 1044 INTEGER :: idim, icat 1045 !!------------------------------------------------------------------- 1046 ! 1047 idim = SIZE( phti, 1 ) 1048 icat = SIZE( phti, 2 ) 1049 ! 1050 ! == thickness and concentration == ! 1051 ! ! ---------------------- ! 1052 IF( icat == jpl ) THEN ! input cat = output cat ! 1053 ! ! ---------------------- ! 1054 ph_i(:,:) = phti(:,:) 1055 ph_s(:,:) = phts(:,:) 1056 pa_i(:,:) = pati(:,:) 1057 ! 1058 ! == temperature and salinity and ponds == ! 1059 pt_i (:,:) = ptmi (:,:) 1060 pt_s (:,:) = ptms (:,:) 1061 pt_su(:,:) = ptmsu(:,:) 1062 ps_i (:,:) = psmi (:,:) 1063 pa_ip(:,:) = patip(:,:) 1064 ph_ip(:,:) = phtip(:,:) 1065 ! ! ---------------------- ! 1066 ELSEIF( icat == 1 ) THEN ! input cat = 1 ! 1067 ! ! ---------------------- ! 1068 CALL ice_var_itd_1cMc( phti(:,1), phts(:,1), pati (:,1), & 1069 & ph_i(:,:), ph_s(:,:), pa_i (:,:), & 1070 & ptmi(:,1), ptms(:,1), ptmsu(:,1), psmi(:,1), patip(:,1), phtip(:,1), & 1071 & pt_i(:,:), pt_s(:,:), pt_su(:,:), ps_i(:,:), pa_ip(:,:), ph_ip(:,:) ) 1072 ! ! ---------------------- ! 1073 ELSEIF( jpl == 1 ) THEN ! output cat = 1 ! 1074 ! ! ---------------------- ! 1075 CALL ice_var_itd_Nc1c( phti(:,:), phts(:,:), pati (:,:), & 1076 & ph_i(:,1), ph_s(:,1), pa_i (:,1), & 1077 & ptmi(:,:), ptms(:,:), ptmsu(:,:), psmi(:,:), patip(:,:), phtip(:,:), & 1078 & pt_i(:,1), pt_s(:,1), pt_su(:,1), ps_i(:,1), pa_ip(:,1), ph_ip(:,1) ) 1079 ! ! ----------------------- ! 1080 ELSE ! input cat /= output cat ! 1081 ! ! ----------------------- ! 1082 1083 ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) ) ! allocate arrays 1084 ALLOCATE( jlmin(idim), jlmax(idim) ) 1085 1086 ! --- initialize output fields to 0 --- ! 1087 ph_i(1:idim,1:jpl) = 0._wp 1088 ph_s(1:idim,1:jpl) = 0._wp 1089 pa_i(1:idim,1:jpl) = 0._wp 1090 ! 1091 ! --- fill the categories --- ! 1092 ! find where cat-input = cat-output and fill cat-output fields 1093 jlmax(:) = 0 1094 jlmin(:) = 999 1095 jlfil(:,:) = 0 1096 DO jl1 = 1, jpl 1097 DO jl2 = 1, icat 1098 DO ji = 1, idim 1099 IF( hi_max(jl1-1) <= phti(ji,jl2) .AND. hi_max(jl1) > phti(ji,jl2) ) THEN 1100 ! fill the right category 1101 ph_i(ji,jl1) = phti(ji,jl2) 1102 ph_s(ji,jl1) = phts(ji,jl2) 1103 pa_i(ji,jl1) = pati(ji,jl2) 1104 ! record categories that are filled 1105 jlmax(ji) = MAX( jlmax(ji), jl1 ) 1106 jlmin(ji) = MIN( jlmin(ji), jl1 ) 1107 jlfil(ji,jl1) = jl1 1108 ENDIF 1109 END DO 1110 END DO 1111 END DO 1112 ! 1113 ! --- fill the gaps between categories --- ! 1114 ! transfer from categories filled at the previous step to the empty ones in between 1115 DO ji = 1, idim 1116 jl1 = jlmin(ji) 1117 jl2 = jlmax(ji) 1118 IF( jl1 > 1 ) THEN 1119 ! fill the lower cat (jl1-1) 1120 pa_i(ji,jl1-1) = ztrans * pa_i(ji,jl1) 1121 ph_i(ji,jl1-1) = hi_mean(jl1-1) 1122 ! remove from cat jl1 1123 pa_i(ji,jl1 ) = ( 1._wp - ztrans ) * pa_i(ji,jl1) 1124 ENDIF 1125 IF( jl2 < jpl ) THEN 1126 ! fill the upper cat (jl2+1) 1127 pa_i(ji,jl2+1) = ztrans * pa_i(ji,jl2) 1128 ph_i(ji,jl2+1) = hi_mean(jl2+1) 1129 ! remove from cat jl2 1130 pa_i(ji,jl2 ) = ( 1._wp - ztrans ) * pa_i(ji,jl2) 1131 ENDIF 1132 END DO 1133 ! 1134 jlfil2(:,:) = jlfil(:,:) 1135 ! fill categories from low to high 1136 DO jl = 2, jpl-1 1137 DO ji = 1, idim 1138 IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN 1139 ! fill high 1140 pa_i(ji,jl) = ztrans * pa_i(ji,jl-1) 1141 ph_i(ji,jl) = hi_mean(jl) 1142 jlfil(ji,jl) = jl 1143 ! remove low 1144 pa_i(ji,jl-1) = ( 1._wp - ztrans ) * pa_i(ji,jl-1) 1145 ENDIF 1146 END DO 1147 END DO 1148 ! 1149 ! fill categories from high to low 1150 DO jl = jpl-1, 2, -1 1151 DO ji = 1, idim 1152 IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN 1153 ! fill low 1154 pa_i(ji,jl) = pa_i(ji,jl) + ztrans * pa_i(ji,jl+1) 1155 ph_i(ji,jl) = hi_mean(jl) 1156 jlfil2(ji,jl) = jl 1157 ! remove high 1158 pa_i(ji,jl+1) = ( 1._wp - ztrans ) * pa_i(ji,jl+1) 1159 ENDIF 1160 END DO 1161 END DO 1162 ! 1163 DEALLOCATE( jlfil, jlfil2 ) ! deallocate arrays 1164 DEALLOCATE( jlmin, jlmax ) 1165 ! 1166 ! == temperature and salinity == ! 1167 ! 1168 ALLOCATE( z1_ai(idim), z1_vi(idim), z1_vs(idim), ztmp(idim) ) 1169 ! 1170 WHERE( SUM( pa_i(:,:), dim=2 ) /= 0._wp ) ; z1_ai(:) = 1._wp / SUM( pa_i(:,:), dim=2 ) 1171 ELSEWHERE ; z1_ai(:) = 0._wp 1172 END WHERE 1173 WHERE( SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) /= 0._wp ) ; z1_vi(:) = 1._wp / SUM( pa_i(:,:) * ph_i(:,:), dim=2 ) 1174 ELSEWHERE ; z1_vi(:) = 0._wp 1175 END WHERE 1176 WHERE( SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) /= 0._wp ) ; z1_vs(:) = 1._wp / SUM( pa_i(:,:) * ph_s(:,:), dim=2 ) 1177 ELSEWHERE ; z1_vs(:) = 0._wp 1178 END WHERE 1179 ! 1180 ! fill all the categories with the same value 1181 ztmp(:) = SUM( ptmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1182 DO jl = 1, jpl 1183 pt_i (:,jl) = ztmp(:) 1184 END DO 1185 ztmp(:) = SUM( ptms (:,:) * pati(:,:) * phts(:,:), dim=2 ) * z1_vs(:) 1186 DO jl = 1, jpl 1187 pt_s (:,jl) = ztmp(:) 1188 END DO 1189 ztmp(:) = SUM( ptmsu(:,:) * pati(:,:) , dim=2 ) * z1_ai(:) 1190 DO jl = 1, jpl 1191 pt_su(:,jl) = ztmp(:) 1192 END DO 1193 ztmp(:) = SUM( psmi (:,:) * pati(:,:) * phti(:,:), dim=2 ) * z1_vi(:) 1194 DO jl = 1, jpl 1195 ps_i (:,jl) = ztmp(:) 1196 END DO 1197 ! 1198 DEALLOCATE( z1_ai, z1_vi, z1_vs, ztmp ) 1199 ! 1200 ! == ponds == ! 1201 ALLOCATE( zfra(idim) ) 1202 ! keep the same pond fraction atip/ati for each category 1203 WHERE( SUM( pati(:,:), dim=2 ) /= 0._wp ) ; zfra(:) = SUM( patip(:,:), dim=2 ) / SUM( pati(:,:), dim=2 ) 1204 ELSEWHERE ; zfra(:) = 0._wp 1205 END WHERE 1206 DO jl = 1, jpl 1207 pa_ip(:,jl) = zfra(:) * pa_i(:,jl) 1208 END DO 1209 ! keep the same v_ip/v_i ratio for each category 1210 WHERE( SUM( phti(:,:) * pati(:,:), dim=2 ) /= 0._wp ) 1211 zfra(:) = SUM( phtip(:,:) * patip(:,:), dim=2 ) / SUM( phti(:,:) * pati(:,:), dim=2 ) 1212 ELSEWHERE 1213 zfra(:) = 0._wp 1214 END WHERE 1215 DO jl = 1, jpl 1216 WHERE( pa_ip(:,jl) /= 0._wp ) ; ph_ip(:,jl) = zfra(:) * ( ph_i(:,jl) * pa_i(:,jl) ) / pa_ip(:,jl) 1217 ELSEWHERE ; ph_ip(:,jl) = 0._wp 1218 END WHERE 1219 END DO 1220 DEALLOCATE( zfra ) 1221 ! 1222 ENDIF 1223 ! 1224 END SUBROUTINE ice_var_itd_NcMc 1000 1225 1001 1226 #else
Note: See TracChangeset
for help on using the changeset viewer.
