- Timestamp:
- 2014-05-12T22:46:18+02:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4332 r4634 37 37 38 38 REAL(wp) :: epsi10 = 1.e-10_wp ! 39 REAL(wp) :: zzero = 0._wp ! 40 REAL(wp) :: zone = 1._wp ! 39 REAL(wp) :: epsi20 = 1.e-20_wp ! 41 40 42 41 !!---------------------------------------------------------------------- … … 76 75 INTEGER :: layer, nbpac ! local integers 77 76 INTEGER :: ii, ij, iter ! - - 78 REAL(wp) :: ztmelts, zdv, z qold, zfrazb, zweight, zalphai, zindb, zinda, zde ! local scalars77 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde ! local scalars 79 78 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 79 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 81 80 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 82 81 CHARACTER (len = 15) :: fieldid 83 ! 82 83 REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean) 84 REAL(wp) :: zEi ! sea ice specific enthalpy (J/kg) 85 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 86 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 87 84 88 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 85 89 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not … … 95 99 REAL(wp), POINTER, DIMENSION(:) :: zat_i_ac ! total ice fraction 96 100 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 97 REAL(wp), POINTER, DIMENSION(:) :: z dh_frazb ! accretion of frazil ice at the ice bottom101 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 98 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_ac ! relative ice / frazil velocity (1D vector) 99 103 100 REAL(wp), POINTER, DIMENSION(:,:) :: zhice_old ! previous ice thickness101 REAL(wp), POINTER, DIMENSION(:,:) :: zdummy ! dummy thickness of new ice102 REAL(wp), POINTER, DIMENSION(:,:) :: zdhicbot ! thickness of new ice which is accreted vertically103 104 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 104 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl … … 110 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i 111 112 112 REAL(wp), POINTER, DIMENSION(:) :: zqbgow ! heat budget of the open water (negative)113 REAL(wp), POINTER, DIMENSION(:) :: zdhex ! excessively thick accreted sea ice (hlead-hice)114 115 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content 116 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness … … 125 123 CALL wrk_alloc( jpij, zcatac ) ! integer 126 124 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 127 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, z dh_frazb, zvrel_ac, zqbgow, zdhex)128 CALL wrk_alloc( jpij,jpl, z hice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )125 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac ) 126 CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 129 127 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac ) 130 128 CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) … … 154 152 DO ji = 1, jpi 155 153 !Energy of melting q(S,T) [J.m-3] 156 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * REAL( nlay_i )157 154 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 158 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 155 e_i(ji,jj,jk,jl) = zindb * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i ) 156 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 159 157 END DO 160 158 END DO … … 196 194 DO ji = 1, jpi 197 195 198 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0) THEN196 IF ( qlead(ji,jj) < 0._wp ) THEN 199 197 !------------- 200 198 ! Wind stress … … 278 276 DO jj = 1, jpj 279 277 DO ji = 1, jpi 280 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN278 IF ( qlead(ji,jj) < 0._wp ) THEN 281 279 nbpac = nbpac + 1 282 280 npac( nbpac ) = (jj - 1) * jpi + ji … … 290 288 DO ji = mi0(jiindx), mi1(jiindx) 291 289 DO jj = mj0(jjindx), mj1(jjindx) 292 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) >0._wp ) THEN290 IF ( qlead(ji,jj) < 0._wp ) THEN 293 291 jiindex_1d = (jj - 1) * jpi + ji 294 292 ENDIF … … 318 316 END DO ! jl 319 317 320 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) 321 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) 318 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 322 319 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 323 CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) ) 324 CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) ) 320 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 321 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 322 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 325 323 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 326 324 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 325 326 CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac) , hfx_thd, jpi, jpj, npac(1:nbpac) ) 327 CALL tab_2d_1d( nbpac, hfx_tot_1d(1:nbpac) , hfx_tot, jpi, jpj, npac(1:nbpac) ) 327 328 328 329 !------------------------------------------------------------------------------! 329 330 ! 5) Compute thickness, salinity, enthalpy, age, area and volume of new ice 330 331 !------------------------------------------------------------------------------! 332 333 !----------------------------------------- 334 ! Keep old ice areas and volume in memory 335 !----------------------------------------- 336 zv_old(:,:) = zv_i_ac(:,:) 337 za_old(:,:) = za_i_ac(:,:) 331 338 332 339 !---------------------- … … 365 372 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 366 373 & - rcp * ( ztmelts - rtt ) ) 374 ! MV HC 2014 comment I dont see why this line below is here... ? 375 ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen 367 376 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) & 368 377 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus … … 375 384 END DO ! ji 376 385 377 !--------------------------378 ! Open water energy budget379 !--------------------------380 DO ji = 1, nbpac381 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0382 END DO ! ji383 384 386 !------------------- 385 387 ! Volume of new ice 386 388 !------------------- 387 389 DO ji = 1, nbpac 388 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 390 391 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 392 393 zEw = rcp * ( t_bo_b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b [J/kg] 394 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 395 396 zdE = zEi - zEw ! specific enthalpy difference [J/kg] 397 398 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 399 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 400 zv_newice(ji) = - zfmdt / rhoic 401 402 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux 403 404 ! Contribution to heat flux to the ocean [W.m-2], >0 405 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 406 ! Total heat flux used in this process [W.m-2] 407 hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * zdE * r1_rdtice 408 ! mass flux 409 wfx_opw_1d(ji) = wfx_opw_1d(ji) + zv_newice(ji) * rhoic * r1_rdtice 410 ! salt flux 411 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 389 412 390 413 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 391 414 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 392 z dh_frazb(ji)= zfrazb * zv_newice(ji)415 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 416 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 394 417 END DO … … 402 425 ! 403 426 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 427 !zde = ze_newice(ji) * area(ii,ij) * zv_newice(ji) 404 428 ! 429 ! clem: change that? 405 430 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume 406 431 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy 407 432 408 433 END DO 409 410 ! keep new ice volume in memory411 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )412 434 413 435 !----------------- … … 415 437 !----------------- 416 438 DO ji = 1, nbpac 417 ii = MOD( npac(ji) - 1 , jpi ) + 1418 ij = ( npac(ji) - 1 ) / jpi + 1419 439 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 420 diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem421 440 END DO !ji 422 441 … … 424 443 ! 6) Redistribute new ice area and volume into ice categories ! 425 444 !------------------------------------------------------------------------------! 426 427 !-----------------------------------------428 ! Keep old ice areas and volume in memory429 !-----------------------------------------430 zv_old(:,:) = zv_i_ac(:,:)431 za_old(:,:) = za_i_ac(:,:)432 445 433 446 !------------------------------------------- … … 458 471 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 459 472 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 460 zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl)461 473 zcatac (ji) = jl 462 474 ENDIF 475 zat_i_ac(ji) = zat_i_ac(ji) + za_i_ac (ji,jl) 463 476 END DO 464 477 END DO … … 469 482 DO ji = 1, nbpac 470 483 jl = zcatac(ji) ! categroy in which new ice is put 471 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 472 zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb ! old ice thickness 473 zdhex (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) ) ! difference in thickness 474 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! ice totally new in jl category 484 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) ) ! 0 if old ice 475 485 END DO 476 486 … … 478 488 DO ji = 1, nbpac 479 489 jl = zcatac(ji) 480 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 481 zalphai = MIN( zhice_old(ji,jl) * REAL( jk ) / REAL( nlay_i ), zh_newice(ji) ) & 482 & - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) ) 483 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 484 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / REAL( nlay_i ) & 485 + za_newice(ji) * ze_newice(ji) * zalphai & 486 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) ) 490 ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji) & 491 & + ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_ac(ji,jk,jl) * zv_old(ji,jl) ) / zv_i_ac(ji,jl) 487 492 END DO 488 493 END DO … … 491 496 ! Add excessive volume of new ice at the bottom 492 497 !----------------------------------------------- 493 ! If the ice concentration exceeds 1, the remaining volume of new ice494 ! is equally redistributed among all ice categories in which there is495 ! ice496 497 ! Fraction of level ice498 498 jm = 1 499 zat_i_lev(:) = 0._wp 500 501 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 502 DO ji = 1, nbpac 503 zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 504 END DO 505 END DO 506 507 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 508 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 509 DO ji = 1, nbpac 510 zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) ) 511 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) ) ! clem 512 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 ) 513 END DO 514 END DO 515 IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl) 516 517 !--------------------------------- 518 ! Heat content - bottom accretion 519 !--------------------------------- 520 jm = 1 521 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 522 DO ji = 1, nbpac 523 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) ) ! zindb=1 if ice =0 otherwise 524 zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb 525 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & 526 & + zindb * zdh_frazb(ji) ! frazil ice may coalesce 527 zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice 528 END DO 529 END DO 530 531 ! old layers thicknesses and enthalpies 499 500 ! --- Redistributing energy on the new grid (energy is equally distributed in every layer) --- ! 501 ! DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 502 ! DO jk = 1, nlay_i 503 ! DO ji = 1, nbpac 504 ! ze_i_ac(ji,jk,jl) = ( ze_i_ac(ji,jk,jl) * zv_i_ac(ji,jl) + ze_newice(ji) * ( zdv_res(ji) + zv_frazb(ji) ) ) / & 505 ! & ( zv_i_ac(ji,jl) + ( zdv_res(ji) + zv_frazb(ji) ) ) 506 ! END DO 507 ! END DO 508 ! END DO 509 510 ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 1 --- ! 532 511 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 533 512 DO jk = 1, nlay_i 534 513 DO ji = 1, nbpac 535 zthick0(ji,jk,jl) = z hice_old(ji,jl) / REAL( nlay_i )514 zthick0(ji,jk,jl) = zv_i_ac(ji,jl) / REAL( nlay_i ) 536 515 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 537 516 END DO 538 517 END DO 539 518 END DO 540 !!gm ??? why the previous do loop if ocerwriten by the following one ?541 519 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 542 520 DO ji = 1, nbpac 543 zthick0(ji,nlay_i+1,jl) = zdhicbot(ji,jl) 544 zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zdhicbot(ji,jl) 521 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) ) 522 zthick0(ji,nlay_i+1,jl) = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 ) 523 zqm0 (ji,nlay_i+1,jl) = ze_newice(ji) * zthick0(ji,nlay_i+1,jl) 545 524 END DO ! ji 546 525 END DO ! jl 547 526 548 ! Redistributing energy on the new grid549 527 ze_i_ac(:,:,:) = 0._wp 550 528 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) … … 552 530 DO layer = 1, nlay_i + 1 553 531 DO ji = 1, nbpac 554 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 555 ! Redistributing energy on the new grid 556 zweight = MAX ( MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) ) & 557 & - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp ) & 558 & /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb 559 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 560 END DO ! ji 561 END DO ! layer 562 END DO ! jk 563 END DO ! jl 564 532 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zthick0(ji,layer,jl) + epsi10 ) ) 533 zweight = zindb * MAX( 0._wp, & 534 & MIN( zv_i_ac(ji,jl) * REAL( layer ), ( zv_i_ac(ji,jl) + zthick0(ji,nlay_i+1,jl) ) * REAL( jk ) ) & 535 & - MAX( zv_i_ac(ji,jl) * REAL( layer - 1 ), ( zv_i_ac(ji,jl) + zthick0(ji,nlay_i+1,jl) ) * REAL( jk - 1 ) ) ) & 536 & / ( REAL( nlay_i ) * MAX( zthick0(ji,layer,jl), epsi10 ) ) 537 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl) 538 END DO 539 END DO 540 END DO 541 END DO 542 543 ! --- new volumes and layer thickness --- 544 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 545 DO ji = 1, nbpac 546 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) ) 547 zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 ) 548 END DO 549 END DO 550 551 ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 2 --- ! 565 552 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 566 553 DO jk = 1, nlay_i 567 554 DO ji = 1, nbpac 568 555 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 569 ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) & 570 & / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb 556 ze_i_ac(ji,jk,jl) = zindb * ze_i_ac(ji,jk,jl) / MAX( zv_i_ac(ji,jl), epsi10 ) * REAL( nlay_i ) 571 557 END DO 572 558 END DO 573 559 END DO 560 574 561 575 562 !------------ … … 589 576 DO jl = 1, jpl 590 577 DO ji = 1, nbpac 591 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes592 578 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 593 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif579 zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) 594 580 END DO 595 581 END DO 596 582 !clem ENDIF 597 598 !--------------------------------599 ! Update mass/salt fluxes (clem)600 !--------------------------------601 DO jl = 1, jpl602 DO ji = 1, nbpac603 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes604 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl)605 rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb606 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb607 END DO608 END DO609 583 610 584 !------------------------------------------------------------------------------! … … 615 589 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 616 590 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 617 !clem IF ( num_sal == 2 ) & 618 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 591 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 619 592 DO jk = 1, nlay_i 620 593 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj ) 621 594 END DO 622 595 END DO 623 CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 624 CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 596 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 597 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 598 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 599 600 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 601 CALL tab_1d_2d( nbpac, hfx_tot, npac(1:nbpac), hfx_tot_1d(1:nbpac), jpi, jpj ) 625 602 ! 626 603 ENDIF ! nbpac > 0 … … 630 607 !------------------------------------------------------------------------------! 631 608 DO jl = 1, jpl 632 DO jk = 1, nlay_i ! heat content in 10^9 Joules 633 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i ) / unit_fac 609 DO jk = 1, nlay_i 610 DO jj = 1, jpj 611 DO ji = 1, jpi 612 ! heat content in Joules 613 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac ) 614 !e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) ) 615 END DO 616 END DO 634 617 END DO 635 618 END DO … … 669 652 CALL wrk_dealloc( jpij, zcatac ) ! integer 670 653 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 671 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, z dh_frazb, zvrel_ac, zqbgow, zdhex)672 CALL wrk_dealloc( jpij,jpl, z hice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )654 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac ) 655 CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 673 656 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac ) 674 657 CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
Note: See TracChangeset
for help on using the changeset viewer.
