- Timestamp:
- 2014-11-27T16:21:44+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r4688 r4899 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 !!---------------------------------------------------------------------- … … 66 66 INTEGER, INTENT(in) :: kt ! time step index 67 67 ! 68 INTEGER :: ji, jj, jk, jl, ja, jm, jbnd1, jbnd2 ! ice typesdummy loop index68 INTEGER :: ji, jj, jk, jl ! dummy loop index 69 69 ! 70 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b … … 86 86 ! Given thermodynamic growth rates, transport ice between 87 87 ! 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 88 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl, kt ) 93 89 ! 94 90 CALL lim_var_glo2eqv ! only for info … … 123 119 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_th : sm_i : ') 124 120 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_th : smv_i : ') 125 DO j a= 1, nlay_i121 DO jk = 1, nlay_i 126 122 CALL prt_ctl_info(' ') 127 CALL prt_ctl_info(' - Layer : ', ivar1=j a)123 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 128 124 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 : ')125 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : t_i : ') 126 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_th : e_i : ') 131 127 END DO 132 128 END DO … … 140 136 ! 141 137 142 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, ntyp,kt )138 SUBROUTINE lim_itd_th_rem( klbnd, kubnd, kt ) 143 139 !!------------------------------------------------------------------ 144 140 !! *** ROUTINE lim_itd_th_rem *** … … 153 149 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 154 150 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied 155 INTEGER , INTENT (in) :: ntyp ! Number of the type used156 151 INTEGER , INTENT (in) :: kt ! Ocean time step 157 152 ! … … 171 166 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 172 167 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 173 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_ o! old ice thickness168 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 174 169 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es 175 170 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume … … 189 184 CALL wrk_alloc( jpi,jpj, zremap_flag ) ! integer 190 185 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 )186 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 192 187 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 193 188 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 218 213 WRITE(numout,*) ' klbnd : ', klbnd 219 214 WRITE(numout,*) ' kubnd : ', kubnd 220 WRITE(numout,*) ' ntyp : ', ntyp221 215 ENDIF 222 216 … … 227 221 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 228 222 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * zindb 229 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)223 zindb = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes 224 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * zindb 225 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 232 226 END DO 233 227 END DO … … 274 268 ! 275 269 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 ) THEN270 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 277 271 !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) THEN272 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 273 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 274 ELSEIF ( a_i_b(ii,ij,jl) > epsi10) THEN 281 275 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 282 ELSEIF ( old_a_i(ii,ij,jl+1) > epsi10) THEN276 ELSEIF ( a_i_b(ii,ij,jl+1) > epsi10) THEN 283 277 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 284 278 ENDIF … … 321 315 DO jj = 1, jpj 322 316 DO ji = 1, jpi 323 zhb0(ji,jj) = hi_max _typ(0,ntyp) ! 0eme324 zhb1(ji,jj) = hi_max _typ(1,ntyp) ! 1er317 zhb0(ji,jj) = hi_max(0) ! 0eme 318 zhb1(ji,jj) = hi_max(1) ! 1er 325 319 326 320 zhbnew(ji,jj,klbnd-1) = 0._wp … … 343 337 !----------------------------------------------------------------------------------------------- 344 338 !- 7.1 g(h) for category 1 at start of time step 345 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_ o(:,:,klbnd), &339 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), & 346 340 & g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), & 347 341 & hR(:,:,klbnd), zremap_flag ) … … 368 362 ! Constrain new thickness <= ht_i 369 363 zdamax = a_i(ii,ij,klbnd) * & 370 (1.0 - ht_i(ii,ij,klbnd)/zht_i_ o(ii,ij,klbnd)) ! zdamax > 0364 (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0 371 365 !ice area lost due to melting of thin ice 372 366 zda0 = MIN(zda0, zdamax) … … 382 376 ELSE ! if ice accretion 383 377 ! ji, a_i > epsi10; zdh0 > 0 384 IF ( ntyp .EQ. 1 )zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))378 zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 385 379 ! zhbnew was 0, and is shifted to the right to account for thin ice 386 380 ! 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 381 ENDIF ! zdh0 391 382 … … 493 484 CALL wrk_dealloc( jpi,jpj, zremap_flag ) ! integer 494 485 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 )486 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b, dummy_es ) 496 487 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 497 488 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) … … 839 830 840 831 841 SUBROUTINE lim_itd_th_reb( klbnd, kubnd , ntyp)832 SUBROUTINE lim_itd_th_reb( klbnd, kubnd ) 842 833 !!------------------------------------------------------------------ 843 834 !! *** ROUTINE lim_itd_th_reb *** … … 849 840 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point 850 841 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 842 ! 853 843 INTEGER :: ji,jj, jl ! dummy loop indices … … 889 879 890 880 !------------------------------------------------------------------------------ 891 ! 2) Make sure thickness of cat klbnd is at least hi_max _typ(klbnd)881 ! 2) Make sure thickness of cat klbnd is at least hi_max(klbnd) 892 882 !------------------------------------------------------------------------------ 893 883 DO jj = 1, jpj 894 884 DO ji = 1, jpi 895 885 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)886 IF( ht_i(ji,jj,klbnd) <= hi_max(0) .AND. hi_max(0) > 0._wp ) THEN 887 a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max(0) 888 ht_i(ji,jj,klbnd) = hi_max(0) 899 889 ENDIF 900 890 ENDIF
Note: See TracChangeset
for help on using the changeset viewer.
