- Timestamp:
- 2017-09-06T12:01:34+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
r8486 r8500 29 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 30 USE timing ! Timing 31 USE iom ! 31 32 32 33 IMPLICIT NONE … … 68 69 ! 69 70 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 70 REAL(wp) :: zdv , zda71 REAL(wp) :: zdv 71 72 REAL(wp), DIMENSION(jpi,jpj) :: zatold, zeiold, zesold, zsmvold 72 73 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zhimax, zviold, zvsold … … 80 81 ! MV MP 2016 81 82 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z0ap , z0vp 82 REAL(wp) :: za_old83 83 ! END MV MP 2016 84 84 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: z0ei … … 111 111 !--- Thickness correction init. --- ! 112 112 zatold(:,:) = at_i 113 DO jl = 1, jpl 114 DO jj = 1, jpj 115 DO ji = 1, jpi 116 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 117 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 118 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 119 END DO 120 END DO 121 END DO 113 WHERE( a_i(:,:,:) >= epsi20 ) 114 ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:) 115 ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:) 116 ELSEWHERE 117 ht_i(:,:,:) = 0._wp 118 ht_s(:,:,:) = 0._wp 119 END WHERE 120 122 121 ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- ! 123 122 zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:) … … 126 125 DO ji = 2, jpim1 127 126 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values 128 zhimax(ji,jj,jl) = MAX VAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )127 zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) ) 129 128 END DO 130 129 END DO … … 399 398 END DO 400 399 END DO 401 402 IF( nn_limdyn == 2) THEN 400 IF( iom_use('icetrp') ) CALL iom_put( "icetrp" , diag_trp_vi * rday ) ! ice volume transport 401 IF( iom_use('snwtrp') ) CALL iom_put( "snwtrp" , diag_trp_vs * rday ) ! snw volume transport 402 IF( iom_use('saltrp') ) CALL iom_put( "saltrp" , diag_trp_smv * rday * rhoic ) ! salt content transport 403 IF( iom_use('deitrp') ) CALL iom_put( "deitrp" , diag_trp_ei ) ! advected ice enthalpy (W/m2) 404 IF( iom_use('destrp') ) CALL iom_put( "destrp" , diag_trp_es ) ! advected snw enthalpy (W/m2) 405 406 IF( nn_limdyn == 2 ) THEN 403 407 ! 404 408 CALL ice_var_zapsmall !--- zap small areas ---! … … 407 411 DO jj = 1, jpj 408 412 DO ji = 1, jpi 409 !410 413 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 411 414 ! 412 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 413 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 414 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 415 ! 416 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 415 ht_i (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 416 ht_s (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl) 417 zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl) 417 418 ! 418 419 IF ( ( zdv > 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. & 419 420 & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN 420 ! 421 rswitch = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) ) 422 a_i(ji,jj,jl) = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 ) 423 ! 424 ! small correction due to *rswitch for a_i 425 v_i (ji,jj,jl) = rswitch * v_i (ji,jj,jl) 426 v_s (ji,jj,jl) = rswitch * v_s (ji,jj,jl) 427 smv_i(ji,jj,jl) = rswitch * smv_i(ji,jj,jl) 428 e_s(ji,jj,1,jl) = rswitch * e_s(ji,jj,1,jl) 429 e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl) 430 431 ! MV MP 2016 432 IF ( nn_pnd_scheme > 0 ) THEN 433 a_ip (ji,jj,jl) = rswitch * a_ip (ji,jj,jl) 434 v_ip (ji,jj,jl) = rswitch * v_ip (ji,jj,jl) 435 ENDIF 436 ! END MV MP 2016 437 ! 421 a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl) 422 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 438 423 ENDIF 439 424 ! 440 425 ENDIF 441 !442 426 END DO 443 427 END DO 444 428 END DO 445 446 DO jj = 1, jpj !--- bound ht_i to hi_max (99 m). 447 DO ji = 1, jpi 448 ! MV MP 2016 449 za_old = a_i(ji,jj,jpl) 450 ! END MV MP 2016 451 rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) ) 452 ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) ) 453 a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch 454 ! MV MP 2016 455 IF ( nn_pnd_scheme > 0 ) THEN 456 ! correct pond fraction to avoid a_ip > a_i 457 a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch 458 ENDIF 459 ! END MP 2016 460 END DO 461 END DO 429 430 WHERE( ht_i(:,:,jpl) > hi_max(jpl) ) !--- bound ht_i to hi_max (99 m) 431 ht_i(:,:,jpl) = hi_max(jpl) 432 a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl) 433 END WHERE 434 435 IF ( nn_pnd_scheme > 0 ) THEN !--- correct pond fraction to avoid a_ip > a_i 436 WHERE( a_ip(:,:,:) > a_i(:,:,:) ) a_ip(:,:,:) = a_i(:,:,:) 437 ENDIF 462 438 ! 463 439 ENDIF … … 467 443 !------------------------------------------------------------ 468 444 ! 469 !!gm remplace the test by, l_piling a logical compute one for all in icestp.F90 (and its declaration in ice.F90 470 !!gm IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2 471 IF( l_piling ) THEN 445 IF( l_piling ) THEN ! simple conservative piling, comparable with 1-cat models 472 446 at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 473 447 DO jl = 1, jpl 474 DO jj = 1, jpj 475 DO ji = 1, jpi 476 rswitch = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) ) 477 zda = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 ) 478 a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda 479 END DO 480 END DO 481 END DO 482 !!gm better and faster coding? 483 ! DO jl = 1, jpl 484 ! WHERE( at_i(:,:) > epsi20 ) 485 ! a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) ) 486 ! END WHERE 487 ! END DO 488 !!gm end 448 WHERE( at_i(:,:) > epsi20 ) 449 a_i(:,:,jl) = a_i(:,:,jl) * ( 1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:) ) 450 END WHERE 451 END DO 489 452 ENDIF 490 453
Note: See TracChangeset
for help on using the changeset viewer.
