Changeset 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
- Timestamp:
- 2015-03-04T17:06:03+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r4990 r5123 30 30 !!====================================================================== 31 31 !! History : - ! 2006-01 (M. Vancoppenolle) Original code 32 !! 4.0! 2011-02 (G. Madec) dynamical allocation32 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 33 33 !!---------------------------------------------------------------------- 34 34 #if defined key_lim3 … … 36 36 !! 'key_lim3' LIM3 sea-ice model 37 37 !!---------------------------------------------------------------------- 38 !! lim_var_agg :39 !! lim_var_glo2eqv :40 !! lim_var_eqv2glo :41 !! lim_var_salprof :42 !! lim_var_salprof1d :43 !! lim_var_bv :44 !!----------------------------------------------------------------------45 38 USE par_oce ! ocean parameters 46 39 USE phycst ! physical constants (ocean directory) 47 40 USE sbc_oce ! Surface boundary condition: ocean fields 48 41 USE ice ! ice variables 49 USE par_ice ! ice parameters50 42 USE thd_ice ! ice variables (thermodynamics) 51 43 USE dom_ice ! ice domain … … 58 50 PRIVATE 59 51 60 PUBLIC lim_var_agg ! 61 PUBLIC lim_var_glo2eqv ! 62 PUBLIC lim_var_eqv2glo ! 63 PUBLIC lim_var_salprof ! 64 PUBLIC lim_var_icetm ! 65 PUBLIC lim_var_bv ! 66 PUBLIC lim_var_salprof1d ! 52 PUBLIC lim_var_agg 53 PUBLIC lim_var_glo2eqv 54 PUBLIC lim_var_eqv2glo 55 PUBLIC lim_var_salprof 56 PUBLIC lim_var_icetm 57 PUBLIC lim_var_bv 58 PUBLIC lim_var_salprof1d 59 PUBLIC lim_var_zapsmall 60 PUBLIC lim_var_itd 67 61 68 62 !!---------------------------------------------------------------------- 69 !! NEMO/LIM3 4.0, UCL - NEMO Consortium (2011)63 !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011) 70 64 !! $Id$ 71 65 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 129 123 DO jj = 1, jpj 130 124 DO ji = 1, jpi 131 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content125 et_s(ji,jj) = et_s(ji,jj) + e_s(ji,jj,1,jl) ! snow heat content 132 126 rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 133 127 smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi10 ) * rswitch ! ice salinity … … 175 169 END DO 176 170 177 IF( n um_sal == 2 )THEN171 IF( nn_icesal == 2 )THEN 178 172 DO jl = 1, jpl 179 173 DO jj = 1, jpj … … 191 185 ! Ice temperatures 192 186 !------------------- 193 !CDIR NOVERRCHK 194 DO jl = 1, jpl 195 !CDIR NOVERRCHK 187 DO jl = 1, jpl 196 188 DO jk = 1, nlay_i 197 !CDIR NOVERRCHK 198 DO jj = 1, jpj 199 !CDIR NOVERRCHK 189 DO jj = 1, jpj 200 190 DO ji = 1, jpi 201 191 ! ! Energy of melting q(S,T) [J.m-3] 202 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 203 zq_i = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 204 zq_i = zq_i * unit_fac !convert units 205 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt ! Ice layer melt temperature 192 rswitch = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 193 zq_i = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 194 ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0 ! Ice layer melt temperature 206 195 ! 207 196 zaaa = cpic ! Conversion q(S,T) -> T (second order equation) 208 zbbb = ( rcp - cpic ) * ( ztmelts - rt t ) + zq_i /rhoic - lfus209 zccc = lfus * (ztmelts-rt t)197 zbbb = ( rcp - cpic ) * ( ztmelts - rt0 ) + zq_i * r1_rhoic - lfus 198 zccc = lfus * (ztmelts-rt0) 210 199 zdiscrim = SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) ) 211 t_i(ji,jj,jk,jl) = rt t+ rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )212 t_i(ji,jj,jk,jl) = MIN( rt t, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt200 t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa ) 201 t_i(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15_wp, t_i(ji,jj,jk,jl) ) ) ! 100-rt0 < t_i < rt0 213 202 END DO 214 203 END DO … … 226 215 DO ji = 1, jpi 227 216 !Energy of melting q(S,T) [J.m-3] 228 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) ) ! rswitch = 0 if no ice and 1 if yes 229 zq_s = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp) 230 zq_s = zq_s * unit_fac ! convert units 217 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) ) ! rswitch = 0 if no ice and 1 if yes 218 zq_s = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp) 231 219 ! 232 t_s(ji,jj,jk,jl) = rt t+ rswitch * ( - zfac1 * zq_s + zfac2 )233 t_s(ji,jj,jk,jl) = MIN( rt t, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rtt < t_i < rtt220 t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * zq_s + zfac2 ) 221 t_s(ji,jj,jk,jl) = MIN( rt0, MAX( 173.15, t_s(ji,jj,jk,jl) ) ) ! 100-rt0 < t_i < rt0 234 222 END DO 235 223 END DO … … 281 269 !! ** Purpose : computes salinity profile in function of bulk salinity 282 270 !! 283 !! ** Method : If bulk salinity greater than s_i_1,271 !! ** Method : If bulk salinity greater than zsi1, 284 272 !! the profile is assumed to be constant (S_inf) 285 !! If bulk salinity lower than s_i_0,273 !! If bulk salinity lower than zsi0, 286 274 !! the profile is linear with 0 at the surface (S_zero) 287 !! If it is between s_i_0 and s_i_1, it is a275 !! If it is between zsi0 and zsi1, it is a 288 276 !! alpha-weighted linear combination of s_inf and s_zero 289 277 !! 290 !! ** References : Vancoppenolle et al., 2007 (in preparation)278 !! ** References : Vancoppenolle et al., 2007 291 279 !!------------------------------------------------------------------ 292 280 INTEGER :: ji, jj, jk, jl ! dummy loop index 293 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac, zsal ! local scalar 294 REAL(wp) :: zswi0, zswi01, zswibal, zargtemp , zs_zero ! - - 295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha ! 3D pointer 281 REAL(wp) :: zfac0, zfac1, zsal 282 REAL(wp) :: zswi0, zswi01, zargtemp , zs_zero 283 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_slope_s, zalpha 284 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 285 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 296 286 !!------------------------------------------------------------------ 297 287 … … 301 291 ! Vertically constant, constant in time 302 292 !--------------------------------------- 303 IF( n um_sal == 1 ) s_i(:,:,:,:) = bulk_sal293 IF( nn_icesal == 1 ) s_i(:,:,:,:) = rn_icesal 304 294 305 295 !----------------------------------- 306 296 ! Salinity profile, varying in time 307 297 !----------------------------------- 308 IF( n um_sal == 2 ) THEN298 IF( nn_icesal == 2 ) THEN 309 299 ! 310 300 DO jk = 1, nlay_i … … 320 310 END DO 321 311 ! 322 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) ! Weighting factor between zs_zero and zs_inf323 dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )312 zfac0 = 1._wp / ( zsi0 - zsi1 ) ! Weighting factor between zs_zero and zs_inf 313 zfac1 = zsi1 / ( zsi1 - zsi0 ) 324 314 ! 325 315 zalpha(:,:,:) = 0._wp … … 327 317 DO jj = 1, jpj 328 318 DO ji = 1, jpi 329 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise330 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i(ji,jj,jl) ) )331 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws332 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i(ji,jj,jl) ) )333 ! If 2.sm_i GE sss_m then zswibal= 1319 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 320 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i(ji,jj,jl) ) ) 321 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 322 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i(ji,jj,jl) ) ) 323 ! If 2.sm_i GE sss_m then rswitch = 1 334 324 ! this is to force a constant salinity profile in the Baltic Sea 335 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )336 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )337 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal)338 END DO 339 END DO 340 END DO 341 342 dummy_fac = 1._wp / REAL( nlay_i )! Computation of the profile325 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) ) 326 zalpha(ji,jj,jl) = zswi0 + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 ) 327 zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch ) 328 END DO 329 END DO 330 END DO 331 332 ! Computation of the profile 343 333 DO jl = 1, jpl 344 334 DO jk = 1, nlay_i … … 346 336 DO ji = 1, jpi 347 337 ! ! linear profile with 0 at the surface 348 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac338 zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i 349 339 ! ! weighting the profile 350 340 s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl) … … 354 344 END DO ! jl 355 345 ! 356 ENDIF ! n um_sal346 ENDIF ! nn_icesal 357 347 358 348 !------------------------------------------------------- … … 360 350 !------------------------------------------------------- 361 351 362 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)352 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 363 353 ! 364 354 sm_i(:,:,:) = 2.30_wp 365 355 ! 366 356 DO jl = 1, jpl 367 !CDIR NOVERRCHK368 357 DO jk = 1, nlay_i 369 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)358 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 370 359 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ) ) 371 360 s_i(:,:,jk,jl) = zsal … … 373 362 END DO 374 363 ! 375 ENDIF ! n um_sal364 ENDIF ! nn_icesal 376 365 ! 377 366 CALL wrk_dealloc( jpi, jpj, jpl, z_slope_s, zalpha ) … … 397 386 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 398 387 tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl) & 399 & / ( REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 ))388 & * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 ) 400 389 END DO 401 390 END DO … … 425 414 DO jj = 1, jpj 426 415 DO ji = 1, jpi 427 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt t) + epsi10 ) ) )428 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt t, - epsi10 ) &429 & * v_i(ji,jj,jl) / REAL(nlay_i,wp)416 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) ) ) 417 zbvi = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) & 418 & * v_i(ji,jj,jl) * r1_nlay_i 430 419 rswitch = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) ) ) 431 420 bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi / MAX( vt_i(ji,jj) , epsi10 ) … … 448 437 ! 449 438 INTEGER :: ji, jk ! dummy loop indices 450 INTEGER :: ii, ij ! local integers451 REAL(wp) :: dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal ! local scalars452 REAL(wp) :: zalpha, zswi0, zswi01, zs wibal, zs_zero ! - -439 INTEGER :: ii, ij ! local integers 440 REAL(wp) :: zfac0, zfac1, zargtemp, zsal ! local scalars 441 REAL(wp) :: zalpha, zswi0, zswi01, zs_zero ! - - 453 442 ! 454 443 REAL(wp), POINTER, DIMENSION(:) :: z_slope_s 444 REAL(wp), PARAMETER :: zsi0 = 3.5_wp 445 REAL(wp), PARAMETER :: zsi1 = 4.5_wp 455 446 !!--------------------------------------------------------------------- 456 447 … … 460 451 ! Vertically constant, constant in time 461 452 !--------------------------------------- 462 IF( n um_sal == 1 ) s_i_1d(:,:) = bulk_sal453 IF( nn_icesal == 1 ) s_i_1d(:,:) = rn_icesal 463 454 464 455 !------------------------------------------------------ … … 466 457 !------------------------------------------------------ 467 458 468 IF( n um_sal == 2 ) THEN459 IF( nn_icesal == 2 ) THEN 469 460 ! 470 461 DO ji = kideb, kiut ! Slope of the linear profile zs_zero … … 474 465 ! Weighting factor between zs_zero and zs_inf 475 466 !--------------------------------------------- 476 dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 ) 477 dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 ) 478 dummy_fac2 = 1._wp / REAL(nlay_i,wp) 479 480 !CDIR NOVERRCHK 467 zfac0 = 1._wp / ( zsi0 - zsi1 ) 468 zfac1 = zsi1 / ( zsi1 - zsi0 ) 481 469 DO jk = 1, nlay_i 482 !CDIR NOVERRCHK483 470 DO ji = kideb, kiut 484 471 ii = MOD( npb(ji) - 1 , jpi ) + 1 485 472 ij = ( npb(ji) - 1 ) / jpi + 1 486 ! zswi0 = 1 if sm_i le s_i_0 and 0 otherwise487 zswi0 = MAX( 0._wp , SIGN( 1._wp , s_i_0 - sm_i_1d(ji) ) )488 ! zswi01 = 1 if sm_i is between s_i_0 and s_i_1 and 0 othws489 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , s_i_1 - sm_i_1d(ji) ) )490 ! if 2.sm_i GE sss_m then zswibal= 1473 ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise 474 zswi0 = MAX( 0._wp , SIGN( 1._wp , zsi0 - sm_i_1d(ji) ) ) 475 ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 476 zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 477 ! if 2.sm_i GE sss_m then rswitch = 1 491 478 ! this is to force a constant salinity profile in the Baltic Sea 492 zswibal= MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )479 rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) ) 493 480 ! 494 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 ) ) * ( 1.0 - zswibal)481 zalpha = ( zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 ) ) * ( 1._wp - rswitch ) 495 482 ! 496 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2483 zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i 497 484 ! weighting the profile 498 485 s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji) 499 END DO ! ji500 END DO ! jk501 502 ENDIF ! num_sal486 END DO 487 END DO 488 489 ENDIF 503 490 504 491 !------------------------------------------------------- … … 506 493 !------------------------------------------------------- 507 494 508 IF( n um_sal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)495 IF( nn_icesal == 3 ) THEN ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30) 509 496 ! 510 497 sm_i_1d(:) = 2.30_wp 511 498 ! 512 !CDIR NOVERRCHK513 499 DO jk = 1, nlay_i 514 zargtemp = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)515 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) ))500 zargtemp = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i 501 zsal = 1.6_wp * ( 1._wp - COS( rpi * zargtemp**( 0.407_wp / ( 0.573_wp + zargtemp ) ) ) ) 516 502 DO ji = kideb, kiut 517 503 s_i_1d(ji,jk) = zsal … … 524 510 ! 525 511 END SUBROUTINE lim_var_salprof1d 512 513 SUBROUTINE lim_var_zapsmall 514 !!------------------------------------------------------------------- 515 !! *** ROUTINE lim_var_zapsmall *** 516 !! 517 !! ** Purpose : Remove too small sea ice areas and correct fluxes 518 !! 519 !! history : LIM3.5 - 01-2014 (C. Rousset) original code 520 !!------------------------------------------------------------------- 521 INTEGER :: ji, jj, jl, jk ! dummy loop indices 522 REAL(wp) :: zsal, zvi, zvs, zei, zes 523 !!------------------------------------------------------------------- 524 at_i (:,:) = 0._wp 525 DO jl = 1, jpl 526 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 527 END DO 528 529 DO jl = 1, jpl 530 531 !----------------------------------------------------------------- 532 ! Zap ice energy and use ocean heat to melt ice 533 !----------------------------------------------------------------- 534 DO jk = 1, nlay_i 535 DO jj = 1 , jpj 536 DO ji = 1 , jpi 537 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 538 rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch 539 zei = e_i(ji,jj,jk,jl) 540 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch 541 t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 542 ! update exchanges with ocean 543 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0 544 END DO 545 END DO 546 END DO 547 548 DO jj = 1 , jpj 549 DO ji = 1 , jpi 550 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 551 rswitch = MAX( 0._wp , SIGN( 1._wp, at_i(ji,jj ) - epsi10 ) ) * rswitch 552 553 zsal = smv_i(ji,jj, jl) 554 zvi = v_i (ji,jj, jl) 555 zvs = v_s (ji,jj, jl) 556 zes = e_s (ji,jj,1,jl) 557 !----------------------------------------------------------------- 558 ! Zap snow energy 559 !----------------------------------------------------------------- 560 t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch ) 561 e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch 562 563 !----------------------------------------------------------------- 564 ! zap ice and snow volume, add water and salt to ocean 565 !----------------------------------------------------------------- 566 ato_i(ji,jj) = a_i (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj) 567 a_i (ji,jj,jl) = a_i (ji,jj,jl) * rswitch 568 v_i (ji,jj,jl) = v_i (ji,jj,jl) * rswitch 569 v_s (ji,jj,jl) = v_s (ji,jj,jl) * rswitch 570 t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch ) 571 oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch 572 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch 573 574 ! ice salinity must stay in bounds 575 IF( nn_icesal == 2 ) THEN 576 smv_i(ji,jj,jl) = MAX( MIN( rn_simax * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), rn_simin * v_i(ji,jj,jl) ) 577 ENDIF 578 ! update exchanges with ocean 579 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 580 wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 581 wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 582 hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0 583 END DO 584 END DO 585 END DO 586 587 ! to be sure that at_i is the sum of a_i(jl) 588 at_i (:,:) = 0._wp 589 DO jl = 1, jpl 590 at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 591 END DO 592 593 ! open water = 1 if at_i=0 594 DO jj = 1, jpj 595 DO ji = 1, jpi 596 rswitch = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) ) 597 ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj) 598 END DO 599 END DO 600 601 ! 602 END SUBROUTINE lim_var_zapsmall 603 604 SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i ) 605 !!------------------------------------------------------------------ 606 !! *** ROUTINE lim_var_itd *** 607 !! 608 !! ** Purpose : converting 1-cat ice to multiple ice categories 609 !! 610 !! ice thickness distribution follows a gaussian law 611 !! around the concentration of the most likely ice thickness 612 !! (similar as limistate.F90) 613 !! 614 !! ** Method: Iterative procedure 615 !! 616 !! 1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian 617 !! 618 !! 2) Check whether the distribution conserves area and volume, positivity and 619 !! category boundaries 620 !! 621 !! 3) If not (input ice is too thin), the last category is empty and 622 !! the number of categories is reduced (jpl-1) 623 !! 624 !! 4) Iterate until ok (SUM(itest(:) = 4) 625 !! 626 !! ** Arguments : zhti: 1-cat ice thickness 627 !! zhts: 1-cat snow depth 628 !! zai : 1-cat ice concentration 629 !! 630 !! ** Output : jpl-cat 631 !! 632 !! (Example of application: BDY forcings when input are cell averaged) 633 !! 634 !!------------------------------------------------------------------- 635 !! History : LIM3.5 - 2012 (M. Vancoppenolle) Original code 636 !! 2014 (C. Rousset) Rewriting 637 !!------------------------------------------------------------------- 638 !! Local variables 639 INTEGER :: ji, jk, jl ! dummy loop indices 640 INTEGER :: ijpij, i_fill, jl0 641 REAL(wp) :: zarg, zV, zconv, zdh 642 REAL(wp), DIMENSION(:), INTENT(in) :: zhti, zhts, zai ! input ice/snow variables 643 REAL(wp), DIMENSION(:,:), INTENT(inout) :: zht_i, zht_s, za_i ! output ice/snow variables 644 INTEGER , POINTER, DIMENSION(:) :: itest 645 646 CALL wrk_alloc( 4, itest ) 647 !-------------------------------------------------------------------- 648 ! initialisation of variables 649 !-------------------------------------------------------------------- 650 ijpij = SIZE(zhti,1) 651 zht_i(1:ijpij,1:jpl) = 0._wp 652 zht_s(1:ijpij,1:jpl) = 0._wp 653 za_i (1:ijpij,1:jpl) = 0._wp 654 655 ! ---------------------------------------- 656 ! distribution over the jpl ice categories 657 ! ---------------------------------------- 658 DO ji = 1, ijpij 659 660 IF( zhti(ji) > 0._wp ) THEN 661 662 ! initialisation of tests 663 itest(:) = 0 664 665 i_fill = jpl + 1 !==================================== 666 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 667 ! iteration !==================================== 668 i_fill = i_fill - 1 669 670 ! initialisation of ice variables for each try 671 zht_i(ji,1:jpl) = 0._wp 672 za_i (ji,1:jpl) = 0._wp 673 674 ! *** case very thin ice: fill only category 1 675 IF ( i_fill == 1 ) THEN 676 zht_i(ji,1) = zhti(ji) 677 za_i (ji,1) = zai (ji) 678 679 ! *** case ice is thicker: fill categories >1 680 ELSE 681 682 ! Fill ice thicknesses except the last one (i_fill) by hmean 683 DO jl = 1, i_fill - 1 684 zht_i(ji,jl) = hi_mean(jl) 685 END DO 686 687 ! find which category (jl0) the input ice thickness falls into 688 jl0 = i_fill 689 DO jl = 1, i_fill 690 IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN 691 jl0 = jl 692 CYCLE 693 ENDIF 694 END DO 695 696 ! Concentrations in the (i_fill-1) categories 697 za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl)) 698 DO jl = 1, i_fill - 1 699 IF ( jl == jl0 ) CYCLE 700 zarg = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp ) 701 za_i(ji,jl) = za_i (ji,jl0) * EXP(-zarg**2) 702 END DO 703 704 ! Concentration in the last (i_fill) category 705 za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) ) 706 707 ! Ice thickness in the last (i_fill) category 708 zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) ) 709 zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / za_i(ji,i_fill) 710 711 ENDIF ! case ice is thick or thin 712 713 !--------------------- 714 ! Compatibility tests 715 !--------------------- 716 ! Test 1: area conservation 717 zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) ) 718 IF ( zconv < epsi06 ) itest(1) = 1 719 720 ! Test 2: volume conservation 721 zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) ) 722 IF ( zconv < epsi06 ) itest(2) = 1 723 724 ! Test 3: thickness of the last category is in-bounds ? 725 IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 726 727 ! Test 4: positivity of ice concentrations 728 itest(4) = 1 729 DO jl = 1, i_fill 730 IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0 731 END DO 732 !============================ 733 END DO ! end iteration on categories 734 !============================ 735 ENDIF ! if zhti > 0 736 END DO ! i loop 737 738 ! ------------------------------------------------ 739 ! Adding Snow in each category where za_i is not 0 740 ! ------------------------------------------------ 741 DO jl = 1, jpl 742 DO ji = 1, ijpij 743 IF( za_i(ji,jl) > 0._wp ) THEN 744 zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) ) 745 ! In case snow load is in excess that would lead to transformation from snow to ice 746 ! Then, transfer the snow excess into the ice (different from limthd_dh) 747 zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 748 ! recompute ht_i, ht_s avoiding out of bounds values 749 zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh ) 750 zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn ) 751 ENDIF 752 ENDDO 753 ENDDO 754 755 CALL wrk_dealloc( 4, itest ) 756 ! 757 END SUBROUTINE lim_var_itd 758 526 759 527 760 #else … … 542 775 SUBROUTINE lim_var_salprof1d ! Emtpy routines 543 776 END SUBROUTINE lim_var_salprof1d 777 SUBROUTINE lim_var_zapsmall 778 END SUBROUTINE lim_var_zapsmall 779 SUBROUTINE lim_var_itd 780 END SUBROUTINE lim_var_itd 544 781 #endif 545 782
Note: See TracChangeset
for help on using the changeset viewer.
