Changeset 7646 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
- Timestamp:
- 2017-02-06T10:25:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r6470 r7646 27 27 !! - et_i(jpi,jpj) !total ice thermal content 28 28 !! - smt_i(jpi,jpj) !mean ice salinity 29 !! - ot_i(jpi,jpj) !average ice age29 !! - tm_i (jpi,jpj) !mean ice temperature 30 30 !!====================================================================== 31 31 !! History : - ! 2006-01 (M. Vancoppenolle) Original code … … 41 41 USE ice ! ice variables 42 42 USE thd_ice ! ice variables (thermodynamics) 43 USE dom_ice ! ice domain44 43 USE in_out_manager ! I/O manager 45 44 USE lib_mpp ! MPP library … … 54 53 PUBLIC lim_var_eqv2glo 55 54 PUBLIC lim_var_salprof 56 PUBLIC lim_var_icetm57 55 PUBLIC lim_var_bv 58 56 PUBLIC lim_var_salprof1d … … 86 84 !!------------------------------------------------------------------ 87 85 88 !-------------------- 89 ! Compute variables 90 !-------------------- 91 vt_i (:,:) = 0._wp 92 vt_s (:,:) = 0._wp 93 at_i (:,:) = 0._wp 94 ato_i(:,:) = 1._wp 95 ! 96 DO jl = 1, jpl 86 ! integrated values 87 vt_i (:,:) = SUM( v_i, dim=3 ) 88 vt_s (:,:) = SUM( v_s, dim=3 ) 89 at_i (:,:) = SUM( a_i, dim=3 ) 90 et_s(:,:) = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 91 et_i(:,:) = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 92 93 ! open water fraction 94 DO jj = 1, jpj 95 DO ji = 1, jpi 96 ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp ) 97 END DO 98 END DO 99 100 IF( kn > 1 ) THEN 101 102 ! mean ice/snow thickness 97 103 DO jj = 1, jpj 98 104 DO ji = 1, jpi 99 ! 100 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 101 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 102 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 103 ! 104 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 105 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch ! ice thickness 106 END DO 107 END DO 108 END DO 109 110 DO jj = 1, jpj 111 DO ji = 1, jpi 112 ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp ) ! open water fraction 113 END DO 114 END DO 115 116 IF( kn > 1 ) THEN 117 et_s (:,:) = 0._wp 118 ot_i (:,:) = 0._wp 105 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 106 htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 107 htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 108 ENDDO 109 ENDDO 110 111 ! mean temperature (K), salinity and age 119 112 smt_i(:,:) = 0._wp 120 et_i (:,:) = 0._wp 121 ! 113 tm_i(:,:) = 0._wp 114 tm_su(:,:) = 0._wp 115 om_i (:,:) = 0._wp 122 116 DO jl = 1, jpl 117 123 118 DO jj = 1, jpj 124 119 DO ji = 1, jpi 125 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) ) 127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch ! ice salinity 128 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) ) 129 ot_i(ji,jj) = ot_i(ji,jj) + oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi20 ) * rswitch ! ice age 130 END DO 131 END DO 132 END DO 133 ! 134 DO jl = 1, jpl 120 rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 121 tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 122 om_i (ji,jj) = om_i (ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 123 END DO 124 END DO 125 135 126 DO jk = 1, nlay_i 136 et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl) ! ice heat content 137 END DO 138 END DO 127 DO jj = 1, jpj 128 DO ji = 1, jpi 129 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 130 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 131 & / MAX( vt_i(ji,jj) , epsi10 ) 132 smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 133 & / MAX( vt_i(ji,jj) , epsi10 ) 134 END DO 135 END DO 136 END DO 137 END DO 138 tm_i = tm_i + rt0 139 tm_su = tm_su + rt0 139 140 ! 140 141 ENDIF … … 243 244 END DO 244 245 245 !------------------- 246 ! Mean temperature 247 !------------------- 248 vt_i (:,:) = 0._wp 249 DO jl = 1, jpl 250 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 251 END DO 252 253 tm_i(:,:) = 0._wp 254 DO jl = 1, jpl 255 DO jk = 1, nlay_i 256 DO jj = 1, jpj 257 DO ji = 1, jpi 258 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 259 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 260 & / MAX( vt_i(ji,jj) , epsi10 ) 261 END DO 262 END DO 263 END DO 264 END DO 265 tm_i = tm_i + rt0 246 ! integrated values 247 vt_i (:,:) = SUM( v_i, dim=3 ) 248 vt_s (:,:) = SUM( v_s, dim=3 ) 249 at_i (:,:) = SUM( a_i, dim=3 ) 250 266 251 ! 267 252 END SUBROUTINE lim_var_glo2eqv … … 398 383 399 384 400 SUBROUTINE lim_var_icetm 401 !!------------------------------------------------------------------ 402 !! *** ROUTINE lim_var_icetm *** 403 !! 404 !! ** Purpose : computes mean sea ice temperature 385 SUBROUTINE lim_var_bv 386 !!------------------------------------------------------------------ 387 !! *** ROUTINE lim_var_bv *** 388 !! 389 !! ** Purpose : computes mean brine volume (%) in sea ice 390 !! 391 !! ** Method : e = - 0.054 * S (ppt) / T (C) 392 !! 393 !! References : Vancoppenolle et al., JGR, 2007 405 394 !!------------------------------------------------------------------ 406 395 INTEGER :: ji, jj, jk, jl ! dummy loop indices 407 396 !!------------------------------------------------------------------ 408 409 ! Mean sea ice temperature 410 vt_i (:,:) = 0._wp 411 DO jl = 1, jpl 412 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 413 END DO 414 415 tm_i(:,:) = 0._wp 397 ! 398 bvm_i(:,:) = 0._wp 399 bv_i (:,:,:) = 0._wp 416 400 DO jl = 1, jpl 417 401 DO jk = 1, nlay_i 418 402 DO jj = 1, jpj 419 403 DO ji = 1, jpi 420 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 421 tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl) & 422 & / MAX( vt_i(ji,jj) , epsi10 ) 423 END DO 424 END DO 425 END DO 426 END DO 427 tm_i = tm_i + rt0 428 429 END SUBROUTINE lim_var_icetm 430 431 432 SUBROUTINE lim_var_bv 433 !!------------------------------------------------------------------ 434 !! *** ROUTINE lim_var_bv *** 435 !! 436 !! ** Purpose : computes mean brine volume (%) in sea ice 437 !! 438 !! ** Method : e = - 0.054 * S (ppt) / T (C) 439 !! 440 !! References : Vancoppenolle et al., JGR, 2007 441 !!------------------------------------------------------------------ 442 INTEGER :: ji, jj, jk, jl ! dummy loop indices 443 REAL(wp) :: zbvi ! local scalars 444 !!------------------------------------------------------------------ 445 ! 446 vt_i (:,:) = 0._wp 447 DO jl = 1, jpl 448 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 449 END DO 450 451 bv_i(:,:) = 0._wp 452 DO jl = 1, jpl 453 DO jk = 1, nlay_i 454 DO jj = 1, jpj 455 DO ji = 1, jpi 456 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 457 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) & 458 & * v_i(ji,jj,jl) * r1_nlay_i 459 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) ) ) 460 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi20 ) 461 END DO 404 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 405 bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i & 406 & / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) 407 END DO 408 END DO 409 END DO 410 411 DO jj = 1, jpj 412 DO ji = 1, jpi 413 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 414 bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 ) 462 415 END DO 463 416 END DO … … 683 636 INTEGER :: ji, jk, jl ! dummy loop indices 684 637 INTEGER :: ijpij, i_fill, jl0 685 REAL(wp) :: zarg, zV, zconv, zdh 638 REAL(wp) :: zarg, zV, zconv, zdh, zdv 686 639 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 687 640 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables … … 704 657 IF( zhti(ji) > 0._wp ) THEN 705 658 706 ! initialisation of tests 707 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 708 670 709 i_fill = jpl + 1 !==================================== 710 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 711 ! iteration !==================================== 712 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 713 731 714 ! initialisation of ice variables for each try 715 zht_i(ji,1:jpl) = 0._wp 716 za_i (ji,1:jpl) = 0._wp 717 718 ! *** case very thin ice: fill only category 1 719 IF ( i_fill == 1 ) THEN 720 zht_i(ji,1) = zhti(ji) 721 za_i (ji,1) = zai (ji) 722 723 ! *** case ice is thicker: fill categories >1 724 ELSE 725 726 ! Fill ice thicknesses except the last one (i_fill) by hmean 727 DO jl = 1, i_fill - 1 728 zht_i(ji,jl) = hi_mean(jl) 729 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 730 735 731 ! find which category (jl0) the input ice thickness falls into 732 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 733 741 DO jl = 1, i_fill 734 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 735 jl0 = jl 736 CYCLE 737 ENDIF 738 END DO 739 740 ! Concentrations in the (i_fill-1) categories 741 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 742 DO jl = 1, i_fill - 1 743 IF ( jl == jl0 ) CYCLE 744 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 745 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 746 END DO 747 748 ! Concentration in the last (i_fill) category 749 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 750 751 ! Ice thickness in the last (i_fill) category 752 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 753 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 754 755 ENDIF ! case ice is thick or thin 756 757 !--------------------- 758 ! Compatibility tests 759 !--------------------- 760 ! Test 1: area conservation 761 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 762 IF ( zconv < epsi06 ) itest(1) = 1 763 764 ! Test 2: volume conservation 765 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 766 IF ( zconv < epsi06 ) itest(2) = 1 767 768 ! Test 3: thickness of the last category is in-bounds ? 769 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 770 771 ! Test 4: positivity of ice concentrations 772 itest(4) = 1 773 DO jl = 1, i_fill 774 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 775 END DO 776 !============================ 777 END DO ! end iteration on categories 778 !============================ 742 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 743 END DO 744 ! !============================ 745 END DO ! end iteration on categories 746 ! !============================ 779 747 ENDIF ! if zhti > 0 780 748 END DO ! i loop 781 749 782 750 ! ------------------------------------------------ 783 751 ! Adding Snow in each category where za_i is not 0
Note: See TracChangeset
for help on using the changeset viewer.
