- Timestamp:
- 2017-09-22T16:55:24+02:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
r8534 r8559 42 42 !! ice_var_eqv2glo : transform from VEQV to VGLO 43 43 !! ice_var_salprof : salinity profile in the ice 44 !! ice_var_bv : brine volume45 44 !! ice_var_salprof1d : salinity profile in the ice 1D 46 45 !! ice_var_zapsmall : remove very small area and volume 47 46 !! ice_var_itd : convert 1-cat to multiple cat 47 !! ice_var_bv : brine volume 48 48 !!---------------------------------------------------------------------- 49 49 USE dom_oce ! ocean space and time domain … … 64 64 PUBLIC ice_var_eqv2glo 65 65 PUBLIC ice_var_salprof 66 PUBLIC ice_var_bv67 66 PUBLIC ice_var_salprof1d 68 67 PUBLIC ice_var_zapsmall 69 68 PUBLIC ice_var_itd 69 PUBLIC ice_var_bv 70 70 71 71 !!---------------------------------------------------------------------- … … 170 170 END WHERE 171 171 ! 172 WHERE( v_i(:,:,:) > epsi20 ) ; z1_v_i(:,:,:) = 1._wp / v_i(:,:,:) 173 ELSEWHERE ; z1_v_i(:,:,:) = 0._wp 174 END WHERE 175 ! 172 176 ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) !--- ice thickness 173 177 … … 177 181 ht_i (:,:,jpl) = zhmax 178 182 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 179 z1_a_i(:,:,jpl) = zhmax /v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max183 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) ! NB: v_i always /=0 as ht_i > hi_max 180 184 END WHERE 181 185 … … 185 189 186 190 IF( nn_icesal == 2 ) THEN !--- salinity (with a minimum value imposed everywhere) 187 WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )191 WHERE( v_i(:,:,:) > epsi20 ) ; sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) ) 188 192 ELSEWHERE ; sm_i(:,:,:) = rn_simin 189 193 END WHERE … … 202 206 IF ( v_i(ji,jj,jl) > epsi20 ) THEN !--- icy area 203 207 ! 204 ze_i = e_i(ji,jj,jk,jl) /v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3]208 ze_i = e_i(ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i ! Energy of melting e(S,T) [J.m-3] 205 209 ztmelts = - s_i(ji,jj,jk,jl) * tmut ! Ice layer melt temperature [C] 206 210 ! Conversion q(S,T) -> T (second order equation) … … 295 299 ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) ) 296 300 ! 297 DO jk = 1, nlay_i 298 s_i(:,:,jk,:) = sm_i(:,:,:) 301 DO jl = 1, jpl 302 DO jk = 1, nlay_i 303 s_i(:,:,jk,jl) = sm_i(:,:,jl) 304 END DO 299 305 END DO 300 306 ! ! Slope of the linear profile … … 353 359 END SUBROUTINE ice_var_salprof 354 360 355 356 SUBROUTINE ice_var_bv357 !!-------------------------------------------------------------------358 !! *** ROUTINE ice_var_bv ***359 !!360 !! ** Purpose : computes mean brine volume (%) in sea ice361 !!362 !! ** Method : e = - 0.054 * S (ppt) / T (C)363 !!364 !! References : Vancoppenolle et al., JGR, 2007365 !!-------------------------------------------------------------------366 INTEGER :: ji, jj, jk, jl ! dummy loop indices367 !!-------------------------------------------------------------------368 !369 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done370 !! instead of setting everything to zero as just below371 bv_i (:,:,:) = 0._wp372 DO jl = 1, jpl373 DO jk = 1, nlay_i374 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )375 bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )376 END WHERE377 END DO378 END DO379 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)380 ELSEWHERE ; bvm_i(:,:) = 0._wp381 END WHERE382 !383 END SUBROUTINE ice_var_bv384 385 386 361 SUBROUTINE ice_var_salprof1d 387 362 !!------------------------------------------------------------------- … … 393 368 INTEGER :: ji, jk ! dummy loop indices 394 369 REAL(wp) :: zargtemp, zsal, z1_dS ! local scalars 395 REAL(wp) :: z alpha, zs, zs0 ! - -396 ! 397 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s !370 REAL(wp) :: zs, zs0 ! - - 371 ! 372 REAL(wp), ALLOCATABLE, DIMENSION(:) :: z_slope_s, zalpha ! 398 373 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 399 374 REAL(wp), PARAMETER :: zsi1 = 4.5_wp … … 405 380 CASE( 1 ) ! constant salinity in time and space ! 406 381 ! !---------------------------------------! 407 s_i_1d( :,:) = rn_icesal382 s_i_1d(1:nidx,:) = rn_icesal 408 383 ! 409 384 ! !---------------------------------------------! … … 411 386 ! !---------------------------------------------! 412 387 ! 413 ALLOCATE( z_slope_s(jpij) )414 ! 415 DO ji = 1, nidx ! Slope of the linear profile416 rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ))417 z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )418 END DO419 388 ALLOCATE( z_slope_s(jpij), zalpha(jpij) ) 389 ! 390 ! ! Slope of the linear profile 391 WHERE( ht_i_1d(1:nidx) > epsi20 ) ; z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx) 392 ELSEWHERE ; z_slope_s(1:nidx) = 0._wp 393 END WHERE 394 420 395 z1_dS = 1._wp / ( zsi1 - zsi0 ) 396 DO ji = 1, nidx 397 zalpha(ji) = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) 398 ! ! force a constant profile when SSS too low (Baltic Sea) 399 IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha(ji) = 0._wp 400 END DO 401 ! 402 ! Computation of the profile 421 403 DO jk = 1, nlay_i 422 404 DO ji = 1, nidx 423 zalpha = MAX( 0._wp , MIN( ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp ) ) 424 IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) ) zalpha = 0._wp ! cst profile when SSS too low (Baltic Sea) 425 ! 426 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i ! linear profile with 0 surface value 427 zs = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji) ! weighting the profile 428 s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) ! bounding salinity 429 END DO 430 END DO 431 ! 432 DEALLOCATE( z_slope_s ) 405 ! ! linear profile with 0 surface value 406 zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 407 zs = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji) 408 s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) ) 409 END DO 410 END DO 411 ! 412 DEALLOCATE( z_slope_s, zalpha ) 433 413 434 414 ! !-------------------------------------------! … … 436 416 ! !-------------------------------------------! (mean = 2.30) 437 417 ! 438 sm_i_1d( :) = 2.30_wp418 sm_i_1d(1:nidx) = 2.30_wp 439 419 ! 440 420 !!gm cf remark in ice_var_salprof routine, CASE( 3 ) … … 459 439 !!------------------------------------------------------------------- 460 440 INTEGER :: ji, jj, jl, jk ! dummy loop indices 461 REAL(wp) :: zsal, zvi, zvs, zei, zes, zvp441 REAL(wp), DIMENSION(jpi,jpj) :: zswitch 462 442 !!------------------------------------------------------------------- 463 443 ! … … 467 447 ! Zap ice energy and use ocean heat to melt ice 468 448 !----------------------------------------------------------------- 449 WHERE( a_i(:,:,jl) > epsi20 ) ; ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl) 450 ELSEWHERE ; ht_i(:,:,jl) = 0._wp 451 END WHERE 452 ! 453 WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 ) ; zswitch(:,:) = 0._wp 454 ELSEWHERE ; zswitch(:,:) = 1._wp 455 END WHERE 456 469 457 DO jk = 1, nlay_i 470 458 DO jj = 1 , jpj 471 459 DO ji = 1 , jpi 472 !!gm I think we can do better/faster : to be discussed473 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )474 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch475 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch &476 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch477 zei = e_i(ji,jj,jk,jl)478 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch479 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )480 460 ! update exchanges with ocean 481 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 461 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0 462 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj) 463 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 482 464 END DO 483 465 END DO … … 486 468 DO jj = 1 , jpj 487 469 DO ji = 1 , jpi 488 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 489 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch 490 rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch & 491 & / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch 492 zsal = smv_i(ji,jj, jl) 493 zvi = v_i (ji,jj, jl) 494 zvs = v_s (ji,jj, jl) 495 zes = e_s (ji,jj,1,jl) 496 IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) zvp = v_ip (ji,jj ,jl) 470 ! update exchanges with ocean 471 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl) * rhoic * r1_rdtice 472 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoic * r1_rdtice 473 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhosn * r1_rdtice 474 hfx_res(ji,jj) = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 497 475 !----------------------------------------------------------------- 498 476 ! Zap snow energy 499 477 !----------------------------------------------------------------- 500 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch)501 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch478 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) ) 479 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj) 502 480 503 481 !----------------------------------------------------------------- 504 482 ! zap ice and snow volume, add water and salt to ocean 505 483 !----------------------------------------------------------------- 506 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 507 a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch 508 v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch 509 v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch 510 t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 511 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 512 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 513 514 ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch 515 ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch 516 517 ! MV MP 2016 518 IF ( ln_pnd ) THEN 519 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch 520 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch 521 IF ( ln_pnd_fw ) wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_ip(ji,jj,jl) - zvp ) * rhofw * r1_rdtice 522 ENDIF 523 ! END MV MP 2016 524 525 ! update exchanges with ocean 526 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 527 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 528 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 529 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 530 END DO 531 END DO 484 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj) 485 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) 486 v_i (ji,jj,jl) = v_i (ji,jj,jl) * zswitch(ji,jj) 487 v_s (ji,jj,jl) = v_s (ji,jj,jl) * zswitch(ji,jj) 488 t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) ) 489 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj) 490 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj) 491 492 ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj) 493 ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj) 494 495 END DO 496 END DO 497 498 IF( ln_pnd ) THEN 499 DO jj = 1 , jpj 500 DO ji = 1 , jpi 501 IF( ln_pnd_fw ) & 502 & wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 503 a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 504 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 505 END DO 506 END DO 507 ENDIF 508 532 509 END DO 533 510 … … 688 665 END SUBROUTINE ice_var_itd 689 666 667 668 SUBROUTINE ice_var_bv 669 !!------------------------------------------------------------------- 670 !! *** ROUTINE ice_var_bv *** 671 !! 672 !! ** Purpose : computes mean brine volume (%) in sea ice 673 !! 674 !! ** Method : e = - 0.054 * S (ppt) / T (C) 675 !! 676 !! References : Vancoppenolle et al., JGR, 2007 677 !!------------------------------------------------------------------- 678 INTEGER :: ji, jj, jk, jl ! dummy loop indices 679 !!------------------------------------------------------------------- 680 ! 681 !!gm I prefere to use WHERE / ELSEWHERE to set it to zero only where needed <<<=== to be done 682 !! instead of setting everything to zero as just below 683 bv_i (:,:,:) = 0._wp 684 DO jl = 1, jpl 685 DO jk = 1, nlay_i 686 WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 ) 687 bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 ) 688 END WHERE 689 END DO 690 END DO 691 WHERE( vt_i(:,:) > epsi20 ) ; bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:) 692 ELSEWHERE ; bvm_i(:,:) = 0._wp 693 END WHERE 694 ! 695 END SUBROUTINE ice_var_bv 696 697 690 698 #else 691 699 !!----------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
