- Timestamp:
- 2016-08-08T11:02:18+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r6746 r6853 636 636 INTEGER :: ji, jk, jl ! dummy loop indices 637 637 INTEGER :: ijpij, i_fill, jl0 638 REAL(wp) :: zarg, zV, zconv, zdh 638 REAL(wp) :: zarg, zV, zconv, zdh, zdv 639 639 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 640 640 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables 641 641 INTEGER , POINTER, DIMENSION(:) :: itest 642 642 643 C ALLwrk_alloc( 4, itest )643 Call wrk_alloc( 4, itest ) 644 644 !-------------------------------------------------------------------- 645 645 ! initialisation of variables … … 657 657 IF( zhti(ji) > 0._wp ) THEN 658 658 659 ! initialisation of tests 660 itest(:) = 0 659 ! find which category (jl0) the input ice thickness falls into 660 jl0 = jpl 661 DO jl = 1, jpl 662 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 663 jl0 = jl 664 CYCLE 665 ENDIF 666 END DO 667 668 ! initialisation of tests 669 itest(:) = 0 661 670 662 i_fill = jpl + 1 !==================================== 663 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 664 ! iteration !==================================== 665 i_fill = i_fill - 1 671 i_fill = jpl + 1 !==================================== 672 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 673 ! iteration !==================================== 674 i_fill = i_fill - 1 675 676 ! initialisation of ice variables for each try 677 zht_i(ji,1:jpl) = 0._wp 678 za_i (ji,1:jpl) = 0._wp 679 itest(:) = 0 680 681 ! *** case very thin ice: fill only category 1 682 IF ( i_fill == 1 ) THEN 683 zht_i(ji,1) = zhti(ji) 684 za_i (ji,1) = zai (ji) 685 686 ! *** case ice is thicker: fill categories >1 687 ELSE 688 689 ! Fill ice thicknesses in the (i_fill-1) cat by hmean 690 DO jl = 1, i_fill - 1 691 zht_i(ji,jl) = hi_mean(jl) 692 END DO 693 694 ! Concentrations in the (i_fill-1) categories 695 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 696 DO jl = 1, i_fill - 1 697 IF ( jl /= jl0 ) THEN 698 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 699 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 700 ENDIF 701 END DO 702 703 ! Concentration in the last (i_fill) category 704 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 705 706 ! Ice thickness in the last (i_fill) category 707 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 708 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 709 710 ! clem: correction if concentration of upper cat is greater than lower cat 711 ! (it should be a gaussian around jl0 but sometimes it is not) 712 IF ( jl0 /= jpl ) THEN 713 DO jl = jpl, jl0+1, -1 714 IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN 715 zdv = zht_i(ji,jl) * za_i(ji,jl) 716 zht_i(ji,jl ) = 0._wp 717 za_i (ji,jl ) = 0._wp 718 za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 ) 719 END IF 720 ENDDO 721 ENDIF 722 723 ENDIF ! case ice is thick or thin 724 725 !--------------------- 726 ! Compatibility tests 727 !--------------------- 728 ! Test 1: area conservation 729 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 730 IF ( zconv < epsi06 ) itest(1) = 1 666 731 667 ! initialisation of ice variables for each try 668 zht_i(ji,1:jpl) = 0._wp 669 za_i (ji,1:jpl) = 0._wp 670 671 ! *** case very thin ice: fill only category 1 672 IF ( i_fill == 1 ) THEN 673 zht_i(ji,1) = zhti(ji) 674 za_i (ji,1) = zai (ji) 675 676 ! *** case ice is thicker: fill categories >1 677 ELSE 678 679 ! Fill ice thicknesses except the last one (i_fill) by hmean 680 DO jl = 1, i_fill - 1 681 zht_i(ji,jl) = hi_mean(jl) 682 END DO 732 ! Test 2: volume conservation 733 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 734 IF ( zconv < epsi06 ) itest(2) = 1 683 735 684 ! find which category (jl0) the input ice thickness falls into 685 jl0 = i_fill 736 ! Test 3: thickness of the last category is in-bounds ? 737 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 738 739 ! Test 4: positivity of ice concentrations 740 itest(4) = 1 686 741 DO jl = 1, i_fill 687 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 688 jl0 = jl 689 CYCLE 690 ENDIF 691 END DO 692 693 ! Concentrations in the (i_fill-1) categories 694 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 695 DO jl = 1, i_fill - 1 696 IF ( jl == jl0 ) CYCLE 697 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 698 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 699 END DO 700 701 ! Concentration in the last (i_fill) category 702 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 703 704 ! Ice thickness in the last (i_fill) category 705 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 706 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 707 708 ENDIF ! case ice is thick or thin 709 710 !--------------------- 711 ! Compatibility tests 712 !--------------------- 713 ! Test 1: area conservation 714 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 715 IF ( zconv < epsi06 ) itest(1) = 1 716 717 ! Test 2: volume conservation 718 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 719 IF ( zconv < epsi06 ) itest(2) = 1 720 721 ! Test 3: thickness of the last category is in-bounds ? 722 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 723 724 ! Test 4: positivity of ice concentrations 725 itest(4) = 1 726 DO jl = 1, i_fill 727 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 728 END DO 729 !============================ 730 END DO ! end iteration on categories 731 !============================ 742 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 743 END DO 744 ! !============================ 745 END DO ! end iteration on categories 746 ! !============================ 732 747 ENDIF ! if zhti > 0 733 748 END DO ! i loop 734 749 735 750 ! ------------------------------------------------ 736 751 ! Adding Snow in each category where za_i is not 0
Note: See TracChangeset
for help on using the changeset viewer.
