- Timestamp:
- 2013-11-07T11:01:27+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r3764 r4161 19 19 !! lim_itd_shiftice : 20 20 !!---------------------------------------------------------------------- 21 USE par_oce ! ocean parameters22 USE dom_oce ! ocean domain23 USE phycst ! physical constants (ocean directory)24 USE ice ! LIM-3 variables25 USE par_ice ! LIM-3 parameters26 USE dom_ice ! LIM-3 domain27 USE thd_ice ! LIM-3 thermodynamic variables28 USE limthd_lac ! LIM-3 lateral accretion29 USE limvar ! LIM-3 variables30 USE limcons ! LIM-3 conservation31 USE prtctl ! Print control32 USE in_out_manager ! I/O manager33 USE lib_mpp ! MPP library34 USE wrk_nemo ! work arrays35 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)36 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)21 USE dom_ice ! LIM-3 domain 22 USE par_oce ! ocean parameters 23 USE dom_oce ! ocean domain 24 USE phycst ! physical constants (ocean directory) 25 USE thd_ice ! LIM-3 thermodynamic variables 26 USE ice ! LIM-3 variables 27 USE par_ice ! LIM-3 parameters 28 USE limthd_lac ! LIM-3 lateral accretion 29 USE limvar ! LIM-3 variables 30 USE limcons ! LIM-3 conservation 31 USE prtctl ! Print control 32 USE in_out_manager ! I/O manager 33 USE lib_mpp ! MPP library 34 USE wrk_nemo ! work arrays 35 USE lib_fortran ! to use key_nosignedzero 36 USE timing ! Timing 37 37 38 38 IMPLICIT NONE … … 45 45 PUBLIC lim_itd_shiftice 46 46 47 REAL(wp) :: epsi20 = 1 e-20_wp ! constant values48 REAL(wp) :: epsi13 = 1 e-13_wp !49 REAL(wp) :: epsi10 = 1 e-10_wp !47 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 48 REAL(wp) :: epsi13 = 1.e-13_wp ! 49 REAL(wp) :: epsi10 = 1.e-10_wp ! 50 50 51 51 !!---------------------------------------------------------------------- 52 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2010)52 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 53 53 !! $Id$ 54 54 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 67 67 ! 68 68 INTEGER :: jl, ja, jm, jbnd1, jbnd2 ! ice types dummy loop index 69 70 !!------------------------------------------------------------------ 69 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 70 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 71 !!------------------------------------------------------------------ 72 IF( nn_timing == 1 ) CALL timing_start('limitd_th') 73 74 ! ------------------------------- 75 !- check conservation (C Rousset) 76 IF (ln_limdiahsb) THEN 77 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 78 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 79 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 80 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 81 ENDIF 82 !- check conservation (C Rousset) 83 ! ------------------------------- 71 84 72 85 IF( kt == nit000 .AND. lwp ) THEN … … 107 120 d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 108 121 d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 109 122 !?? d_oa_i_thd(:,:,:) = oa_i (:,:,:) - old_oa_i (:,:,:) 110 123 d_smv_i_thd(:,:,:) = 0._wp 111 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 124 IF( num_sal == 2 ) d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 125 126 ! diag only (clem) 127 dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 112 128 113 129 IF(ln_ctl) THEN ! Control print … … 142 158 END DO 143 159 ENDIF 144 160 ! 161 ! ------------------------------- 162 !- check conservation (C Rousset) 163 IF( ln_limdiahsb ) THEN 164 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 165 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 166 167 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 168 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 169 170 zchk_vmin = glob_min(v_i) 171 zchk_amax = glob_max(SUM(a_i,dim=3)) 172 zchk_amin = glob_min(a_i) 173 174 IF(lwp) THEN 175 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_th) = ',(zchk_v_i * rday) 176 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 177 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_th) = ',(zchk_vmin * 1.e-3) 178 IF ( zchk_amax > amax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_th) = ',zchk_amax 179 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_th) = ',zchk_amin 180 ENDIF 181 ENDIF 182 !- check conservation (C Rousset) 183 ! ------------------------------- 184 ! 145 185 !- Recover Old values 146 a_i(:,:,:) = old_a_i (:,:,:)147 v_s(:,:,:) = old_v_s (:,:,:)148 v_i(:,:,:) = old_v_i (:,:,:)149 e_s(:,:,:,:) = old_e_s (:,:,:,:)150 e_i(:,:,:,:) = old_e_i (:,:,:,:)151 ! 186 a_i(:,:,:) = old_a_i (:,:,:) 187 v_s(:,:,:) = old_v_s (:,:,:) 188 v_i(:,:,:) = old_v_i (:,:,:) 189 e_s(:,:,:,:) = old_e_s (:,:,:,:) 190 e_i(:,:,:,:) = old_e_i (:,:,:,:) 191 !?? oa_i(:,:,:) = old_oa_i(:,:,:) 152 192 IF( num_sal == 2 ) smv_i(:,:,:) = old_smv_i(:,:,:) 153 193 ! 194 IF( nn_timing == 1 ) CALL timing_stop('limitd_th') 154 195 END SUBROUTINE lim_itd_th 155 196 ! … … 172 213 ! 173 214 INTEGER :: ji, jj, jl ! dummy loop index 174 INTEGER :: zji, zjj, nd ! local integer 215 INTEGER :: ii, ij ! 2D corresponding indices to ji 216 INTEGER :: nd ! local integer 175 217 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 176 REAL(wp) :: zx2, zwk2, zda0, zetamax , zhimin! - -218 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 177 219 REAL(wp) :: zx3, zareamin, zindb ! - - 178 220 CHARACTER (len = 15) :: fieldid … … 210 252 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 211 253 212 zhimin = 0.1 !minimum ice thickness tolerated by the model213 254 zareamin = epsi10 !minimum area in thickness categories tolerated by the conceptors of the model 214 255 … … 240 281 DO jj = 1, jpj 241 282 DO ji = 1, jpi 242 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl) )) !0 if no ice and 1 if yes283 zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 243 284 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 244 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl) )) !0 if no ice and 1 if yes285 zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 245 286 zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 246 287 IF( a_i(ji,jj,jl) > 1e-6 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl) … … 285 326 DO jl = klbnd, kubnd - 1 286 327 DO ji = 1, nbrem 287 zji = nind_i(ji)288 zjj = nind_j(ji)328 ii = nind_i(ji) 329 ij = nind_j(ji) 289 330 ! 290 IF ( ( zht_i_o( zji,zjj,jl) .GT.epsi10 ) .AND. &291 ( zht_i_o( zji,zjj,jl+1).GT.epsi10 ) ) THEN331 IF ( ( zht_i_o(ii,ij,jl) .GT.epsi10 ) .AND. & 332 ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 292 333 !interpolate between adjacent category growth rates 293 zslope = ( zdhice( zji,zjj,jl+1) - zdhice(zji,zjj,jl) ) / &294 ( zht_i_o ( zji,zjj,jl+1) - zht_i_o (zji,zjj,jl) )295 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + &296 zslope * ( hi_max(jl) - zht_i_o( zji,zjj,jl) )297 ELSEIF (zht_i_o( zji,zjj,jl).gt.epsi10) THEN298 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl)299 ELSEIF (zht_i_o( zji,zjj,jl+1).gt.epsi10) THEN300 zhbnew( zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1)334 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / & 335 ( zht_i_o (ii,ij,jl+1) - zht_i_o (ii,ij,jl) ) 336 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 337 zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 338 ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 339 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 340 ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 341 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 301 342 ELSE 302 zhbnew( zji,zjj,jl) = hi_max(jl)343 zhbnew(ii,ij,jl) = hi_max(jl) 303 344 ENDIF 304 345 END DO … … 307 348 DO ji = 1, nbrem 308 349 ! jl, ji 309 zji = nind_i(ji)310 zjj = nind_j(ji)350 ii = nind_i(ji) 351 ij = nind_j(ji) 311 352 ! jl, ji 312 IF ( ( a_i( zji,zjj,jl) .GT.epsi10) .AND. &313 ( ht_i( zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) &353 IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. & 354 ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 314 355 ) THEN 315 zremap_flag( zji,zjj) = 0316 ELSEIF ( ( a_i( zji,zjj,jl+1) .GT. epsi10 ) .AND. &317 ( ht_i( zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) &356 zremap_flag(ii,ij) = 0 357 ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 358 ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 318 359 ) THEN 319 zremap_flag( zji,zjj) = 0360 zremap_flag(ii,ij) = 0 320 361 ENDIF 321 362 322 363 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 323 364 ! jl, ji 324 IF (zhbnew( zji,zjj,jl).gt.hi_max(jl+1)) THEN325 zremap_flag( zji,zjj) = 0365 IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 366 zremap_flag(ii,ij) = 0 326 367 ENDIF 327 368 ! jl, ji 328 IF (zhbnew( zji,zjj,jl).lt.hi_max(jl-1)) THEN329 zremap_flag( zji,zjj) = 0369 IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 370 zremap_flag(ii,ij) = 0 330 371 ENDIF 331 372 ! jl, ji … … 379 420 !- 7.2 Area lost due to melting of thin ice (first category, klbnd) 380 421 DO ji = 1, nbrem 381 zji = nind_i(ji)382 zjj = nind_j(ji)422 ii = nind_i(ji) 423 ij = nind_j(ji) 383 424 384 425 !ji 385 IF (a_i( zji,zjj,klbnd) .gt. epsi10) THEN386 zdh0 = zdhice( zji,zjj,klbnd) !decrease of ice thickness in the lower category426 IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 427 zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 387 428 ! ji, a_i > epsi10 388 429 IF (zdh0 .lt. 0.0) THEN !remove area from category 1 … … 391 432 392 433 !Integrate g(1) from 0 to dh0 to estimate area melted 393 zetamax = MIN(zdh0,hR( zji,zjj,klbnd)) - hL(zji,zjj,klbnd)434 zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 394 435 IF (zetamax.gt.0.0) THEN 395 436 zx1 = zetamax 396 437 zx2 = 0.5 * zetamax*zetamax 397 zda0 = g1( zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed438 zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 398 439 ! Constrain new thickness <= ht_i 399 zdamax = a_i( zji,zjj,klbnd) * &400 (1.0 - ht_i( zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0440 zdamax = a_i(ii,ij,klbnd) * & 441 (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 401 442 !ice area lost due to melting of thin ice 402 443 zda0 = MIN(zda0, zdamax) 403 444 404 445 ! Remove area, conserving volume 405 ht_i( zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &406 * a_i( zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 )407 a_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0408 v_i( zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)446 ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 447 * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 448 a_i(ii,ij,klbnd) = a_i(ii,ij,klbnd) - zda0 449 v_i(ii,ij,klbnd) = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 409 450 ENDIF ! zetamax > 0 410 451 ! ji, a_i > epsi10 … … 412 453 ELSE ! if ice accretion 413 454 ! ji, a_i > epsi10; zdh0 > 0 414 IF ( ntyp .EQ. 1 ) zhbnew( zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))455 IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 415 456 ! zhbnew was 0, and is shifted to the right to account for thin ice 416 457 ! growth in openwater (F0 = f1) 417 IF ( ntyp .NE. 1 ) zhbnew( zji,zjj,0) = 0458 IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0 418 459 ! in other types there is 419 460 ! no open water growth (F0 = 0) … … 446 487 447 488 DO ji = 1, nbrem 448 zji = nind_i(ji)449 zjj = nind_j(ji)450 451 IF (zhbnew( zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1489 ii = nind_i(ji) 490 ij = nind_j(ji) 491 492 IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 452 493 453 494 ! left and right integration limits in eta space 454 zvetamin(ji) = MAX(hi_max(jl), hL( zji,zjj,jl)) - hL(zji,zjj,jl)455 zvetamax(ji) = MIN(zhbnew( zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)456 zdonor( zji,zjj,jl) = jl495 zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 496 zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 497 zdonor(ii,ij,jl) = jl 457 498 458 499 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl … … 460 501 ! left and right integration limits in eta space 461 502 zvetamin(ji) = 0.0 462 zvetamax(ji) = MIN(hi_max(jl), hR( zji,zjj,jl+1)) - hL(zji,zjj,jl+1)463 zdonor( zji,zjj,jl) = jl + 1503 zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 504 zdonor(ii,ij,jl) = jl + 1 464 505 465 506 ENDIF ! zhbnew(jl) > hi_max(jl) … … 475 516 zwk2 = zwk2 * zetamax 476 517 zx3 = 1.0/3.0 * (zwk2 - zwk1) 477 nd = zdonor( zji,zjj,jl)478 zdaice( zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1479 zdvice( zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + &480 zdaice( zji,zjj,jl)*hL(zji,zjj,nd)518 nd = zdonor(ii,ij,jl) 519 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 520 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 521 zdaice(ii,ij,jl)*hL(ii,ij,nd) 481 522 482 523 END DO ! ji … … 493 534 494 535 DO ji = 1, nbrem 495 zji = nind_i(ji) 496 zjj = nind_j(ji) 497 IF ( ( zhimin .GT. 0.0 ) .AND. & 498 ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 499 ) THEN 500 a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin 501 ht_i(zji,zjj,1) = zhimin 502 v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 536 ii = nind_i(ji) 537 ij = nind_j(ji) 538 IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 539 a_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim 540 ht_i(ii,ij,1) = hiclim 541 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 503 542 ENDIF 504 543 END DO !ji … … 625 664 626 665 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 627 INTEGER :: zji, zjj ! indices when changing from 2D-1D is done666 INTEGER :: ii, ij ! indices when changing from 2D-1D is done 628 667 629 668 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn … … 759 798 760 799 DO ji = 1, nbrem 761 zji = nind_i(ji)762 zjj = nind_j(ji)763 764 jl1 = zdonor( zji,zjj,jl)765 zindb = MAX( 0.0 , SIGN( 1.0 , v_i( zji,zjj,jl1) - epsi10 ) )766 zworka( zji,zjj) = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),epsi10) * zindb800 ii = nind_i(ji) 801 ij = nind_j(ji) 802 803 jl1 = zdonor(ii,ij,jl) 804 zindb = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 805 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb 767 806 IF( jl1 == jl) THEN ; jl2 = jl1+1 768 807 ELSE ; jl2 = jl … … 773 812 !-------------- 774 813 775 a_i( zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl)776 a_i( zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl)814 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 815 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 777 816 778 817 !-------------- … … 780 819 !-------------- 781 820 782 v_i( zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl)783 v_i( zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl)821 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 822 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 784 823 785 824 !-------------- … … 787 826 !-------------- 788 827 789 zdvsnow = v_s( zji,zjj,jl1) * zworka(zji,zjj)790 v_s( zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow791 v_s( zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow828 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 829 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 830 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 792 831 793 832 !-------------------- … … 795 834 !-------------------- 796 835 797 zdesnow = e_s( zji,zjj,1,jl1) * zworka(zji,zjj)798 e_s( zji,zjj,1,jl1) = e_s(zji,zjj,1,jl1) - zdesnow799 e_s( zji,zjj,1,jl2) = e_s(zji,zjj,1,jl2) + zdesnow836 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 837 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 838 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow 800 839 801 840 !-------------- … … 803 842 !-------------- 804 843 805 zdo_aice = oa_i( zji,zjj,jl1) * zdaice(zji,zjj,jl)806 oa_i( zji,zjj,jl1) = oa_i(zji,zjj,jl1) - zdo_aice807 oa_i( zji,zjj,jl2) = oa_i(zji,zjj,jl2) + zdo_aice844 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 845 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 846 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice 808 847 809 848 !-------------- … … 811 850 !-------------- 812 851 813 zdsm_vice = smv_i( zji,zjj,jl1) * zworka(zji,zjj)814 smv_i( zji,zjj,jl1) = smv_i(zji,zjj,jl1) - zdsm_vice815 smv_i( zji,zjj,jl2) = smv_i(zji,zjj,jl2) + zdsm_vice852 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 853 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 854 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice 816 855 817 856 !--------------------- … … 819 858 !--------------------- 820 859 821 zdaTsf = t_su( zji,zjj,jl1) * zdaice(zji,zjj,jl)822 zaTsfn( zji,zjj,jl1) = zaTsfn(zji,zjj,jl1) - zdaTsf823 zaTsfn( zji,zjj,jl2) = zaTsfn(zji,zjj,jl2) + zdaTsf860 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 861 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf 862 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 824 863 825 864 END DO ! ji … … 832 871 !CDIR NODEP 833 872 DO ji = 1, nbrem 834 zji = nind_i(ji)835 zjj = nind_j(ji)836 837 jl1 = zdonor( zji,zjj,jl)873 ii = nind_i(ji) 874 ij = nind_j(ji) 875 876 jl1 = zdonor(ii,ij,jl) 838 877 IF (jl1 .EQ. jl) THEN 839 878 jl2 = jl+1 … … 842 881 ENDIF 843 882 844 zdeice = e_i( zji,zjj,jk,jl1) * zworka(zji,zjj)845 e_i( zji,zjj,jk,jl1) = e_i(zji,zjj,jk,jl1) - zdeice846 e_i( zji,zjj,jk,jl2) = e_i(zji,zjj,jk,jl2) + zdeice883 zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 884 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - zdeice 885 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 847 886 END DO ! ji 848 887 END DO ! jk … … 860 899 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 861 900 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 862 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl) )) !0 if no ice and 1 if yes901 zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 863 902 ELSE 864 903 ht_i(ji,jj,jl) = 0._wp … … 967 1006 zshiftflag = 1 968 1007 zdonor(ji,jj,jl) = jl 969 zdaice(ji,jj,jl) = a_i(ji,jj,jl) 970 zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1008 ! begin TECLIM change 1009 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) 1010 !zdvice(ji,jj,jl) = v_i(ji,jj,jl) 1011 zdaice(ji,jj,jl) = a_i(ji,jj,jl)/2 1012 zdvice(ji,jj,jl) = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 1013 ! end TECLIM change 971 1014 ENDIF 972 1015 END DO ! ji
Note: See TracChangeset
for help on using the changeset viewer.
