Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4688 r5208 6 6 !! History : - ! (W. H. Lipscomb and E.C. Hunke) CICE (c) original code 7 7 !! 3.0 ! 2005-12 (M. Vancoppenolle) adaptation to LIM-3 8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age and types8 !! - ! 2006-06 (M. Vancoppenolle) adaptation to include salt, age 9 9 !! - ! 2007-04 (M. Vancoppenolle) Mass conservation checked 10 10 !!---------------------------------------------------------------------- … … 46 46 PUBLIC lim_itd_shiftice 47 47 48 REAL(wp) :: epsi10 = 1.e-10_wp !49 REAL(wp) :: epsi06 = 1.e-6_wp !50 51 48 !!---------------------------------------------------------------------- 52 49 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) … … 66 63 INTEGER, INTENT(in) :: kt ! time step index 67 64 ! 68 INTEGER :: ji, jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice typesdummy loop index65 INTEGER :: ji, jj, jk, jl ! dummy loop index 69 66 ! 70 67 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b … … 86 83 ! Given thermodynamic growth rates, transport ice between 87 84 ! thickness categories. 88 DO jm = 1, jpm 89 jbnd1 = ice_cat_bounds(jm,1) 90 jbnd2 = ice_cat_bounds(jm,2) 91 IF( ice_ncat_types(jm) > 1 ) CALL lim_itd_th_rem( jbnd1, jbnd2, jm, kt ) 92 END DO 85 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 93 86 ! 94 87 CALL lim_var_glo2eqv ! only for info … … 123 116 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 124 117 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 125 DO j a= 1, nlay_i118 DO jk = 1, nlay_i 126 119 CALL prt_ctl_info(' ') 127 CALL prt_ctl_info(' - Layer : ', ivar1=j a)120 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 128 121 CALL prt_ctl_info(' ~~~~~~~') 129 CALL prt_ctl(tab2d_1=t_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : t_i : ')130 CALL prt_ctl(tab2d_1=e_i(:,:,j a,jl) , clinfo1= ' lim_itd_th : e_i : ')122 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ') 123 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ') 131 124 END DO 132 125 END DO … … 140 133 ! 141 134 142 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp,kt )135 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 143 136 !!------------------------------------------------------------------ 144 137 !! *** ROUTINE lim_itd_th_rem *** … … 153 146 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 154 147 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 155 INTEGER , INTENT (in) :: ntyp ! Number of the type used156 148 INTEGER , INTENT (in) :: kt ! Ocean time step 157 149 ! … … 161 153 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 162 154 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 163 REAL(wp) :: zx3, zareamin , zindb! - -155 REAL(wp) :: zx3, zareamin ! - - 164 156 CHARACTER (len = 15) :: fieldid 165 157 … … 171 163 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 172 164 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 173 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness165 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 174 166 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 175 167 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 189 181 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer 190 182 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) ! integer 191 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )183 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 192 184 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 193 185 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 218 210 WRITE(numout,*) ' klbnd : ', klbnd 219 211 WRITE(numout,*) ' kubnd : ', kubnd 220 WRITE(numout,*) ' ntyp : ', ntyp221 212 ENDIF 222 213 … … 225 216 DO jj = 1, jpj 226 217 DO ji = 1, jpi 227 zindb= 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes228 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb229 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - old_a_i(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes230 zht_i_ o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX( old_a_i(ji,jj,jl), epsi10 ) * zindb231 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_ o(ji,jj,jl)218 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 219 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch 220 rswitch = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 221 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch 222 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 232 223 END DO 233 224 END DO … … 274 265 ! 275 266 zhbnew(ii,ij,jl) = hi_max(jl) 276 IF ( old_a_i(ii,ij,jl) > epsi10 .AND. old_a_i(ii,ij,jl+1) > epsi10 ) THEN267 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 277 268 !interpolate between adjacent category growth rates 278 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_ o(ii,ij,jl+1) - zht_i_o(ii,ij,jl) )279 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_ o(ii,ij,jl) )280 ELSEIF ( old_a_i(ii,ij,jl) > epsi10) THEN269 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 270 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 271 ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 281 272 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 282 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN273 ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 283 274 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 284 275 ENDIF … … 321 312 DO jj = 1, jpj 322 313 DO ji = 1, jpi 323 zhb0(ji,jj) = hi_max _typ(0,ntyp) ! 0eme324 zhb1(ji,jj) = hi_max _typ(1,ntyp) ! 1er314 zhb0(ji,jj) = hi_max(0) ! 0eme 315 zhb1(ji,jj) = hi_max(1) ! 1er 325 316 326 317 zhbnew(ji,jj,klbnd-1) = 0._wp … … 343 334 !----------------------------------------------------------------------------------------------- 344 335 !- 7.1 g(h) for category 1 at start of time step 345 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_ o(:,:,klbnd), &336 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 346 337 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 347 338 & hR(:,:,klbnd), zremap_flag ) … … 368 359 ! Constrain new thickness <= ht_i 369 360 zdamax = a_i(ii,ij,klbnd) * & 370 (1.0 - ht_i(ii,ij,klbnd)/zht_i_ o(ii,ij,klbnd)) ! zdamax > 0361 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 371 362 !ice area lost due to melting of thin ice 372 363 zda0 = MIN(zda0, zdamax) … … 382 373 ELSE ! if ice accretion 383 374 ! ji, a_i > epsi10; zdh0 > 0 384 IF ( ntyp .EQ. 1 )zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))375 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 385 376 ! zhbnew was 0, and is shifted to the right to account for thin ice 386 377 ! growth in openwater (F0 = f1) 387 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0388 ! in other types there is389 ! no open water growth (F0 = 0)390 378 ENDIF ! zdh0 391 379 … … 493 481 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer 494 482 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) ! integer 495 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_ o, dummy_es )483 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 496 484 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 497 485 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 598 586 REAL(wp) :: zdo_aice ! ice age times volume transferred 599 587 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 600 REAL(wp) :: zindsn ! snow or not601 REAL(wp) :: zindb ! ice or not602 588 603 589 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions … … 726 712 727 713 jl1 = zdonor(ii,ij,jl) 728 zindb= MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )729 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb714 rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 715 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch 730 716 IF( jl1 == jl) THEN ; jl2 = jl1+1 731 717 ELSE ; jl2 = jl … … 823 809 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 824 810 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 825 zindsn= 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes811 rswitch = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 826 812 ELSE 827 813 ht_i(ji,jj,jl) = 0._wp … … 839 825 840 826 841 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)827 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 842 828 !!------------------------------------------------------------------ 843 829 !! *** ROUTINE lim_itd_th_reb *** … … 849 835 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 850 836 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 851 INTEGER , INTENT (in) :: ntyp ! number of the ice type involved in the rebinning process852 837 ! 853 838 INTEGER :: ji,jj, jl ! dummy loop indices … … 889 874 890 875 !------------------------------------------------------------------------------ 891 ! 2) Make sure thickness of cat klbnd is at least hi_max _typ(klbnd)876 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 892 877 !------------------------------------------------------------------------------ 893 878 DO jj = 1, jpj 894 879 DO ji = 1, jpi 895 880 IF( a_i(ji,jj,klbnd) > epsi10 ) THEN 896 IF( ht_i(ji,jj,klbnd) <= hi_max _typ(0,ntyp) .AND. hi_max_typ(0,ntyp) > 0._wp ) THEN897 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max _typ(0,ntyp)898 ht_i(ji,jj,klbnd) = hi_max _typ(0,ntyp)881 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 882 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 883 ht_i(ji,jj,klbnd) = hi_max(0) 899 884 ENDIF 900 885 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
