- Timestamp:
- 2014-05-27T11:28:12+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
r4634 r4649 30 30 USE wrk_nemo ! work arrays 31 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 USE limthd_ent 32 33 33 34 IMPLICIT NONE … … 85 86 REAL(wp) :: zEw ! seawater specific enthalpy (J/kg) 86 87 REAL(wp) :: zfmdt ! mass flux x time step (kg/m2, >0 towards ocean) 87 88 INTEGER , POINTER, DIMENSION(:) :: zcatac ! indexes of categories where new ice grows 88 89 REAL(wp) :: zv_newfra 90 91 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 89 92 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 90 93 … … 97 100 REAL(wp), POINTER, DIMENSION(:) :: zdv_res ! residual volume in case of excessive heat budget 98 101 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 99 REAL(wp), POINTER, DIMENSION(:) :: zat_i_ ac! total ice fraction102 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 100 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 101 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_ ac! relative ice / frazil velocity (1D vector)105 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 103 106 104 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_old ! old volume of ice in category jl 105 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_old ! old area of ice in category jl 106 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_ac ! 1-D version of a_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_ac ! 1-D version of v_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_ac ! 1-D version of oa_i 109 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_ac ! 1-D version of smv_i 110 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_ac !: 1-D version of e_i 112 113 REAL(wp), POINTER, DIMENSION(:,:,:) :: zqm0 ! old layer-system heat content 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: zthick0 ! old ice thickness 115 116 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 117 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories 118 REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories 119 REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init ! snow energy summed over categories 109 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 110 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 111 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d ! 1-D version of oa_i 112 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 113 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i 115 120 116 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 121 117 !!-----------------------------------------------------------------------! 122 118 123 CALL wrk_alloc( jpij, zcatac) ! integer119 CALL wrk_alloc( jpij, jcat ) ! integer 124 120 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 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 ) 127 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac ) 128 CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 129 CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 130 131 et_i_init(:,:) = 0._wp 132 et_s_init(:,:) = 0._wp 133 vt_i_init(:,:) = 0._wp 134 vt_s_init(:,:) = 0._wp 135 136 !------------------------------------------------------------------------------! 137 ! 1) Conservation check and changes in each ice category 138 !------------------------------------------------------------------------------! 139 IF( con_i ) THEN 140 CALL lim_column_sum ( jpl, v_i , vt_i_init) 141 CALL lim_column_sum ( jpl, v_s , vt_s_init) 142 CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 143 CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init) 144 ENDIF 121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d ) 124 CALL wrk_alloc( jpi,jpj, zvrel ) 145 125 146 126 !------------------------------------------------------------------------------| … … 204 184 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp 205 185 ! Square root of wind stress 206 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy) )186 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 207 187 208 188 !--------------------- 209 189 ! Frazil ice velocity 210 190 !--------------------- 211 zvfrx = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 212 zvfry = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 191 zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 192 zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 193 zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 213 194 214 195 !------------------- … … 305 286 IF ( nbpac > 0 ) THEN 306 287 307 CALL tab_2d_1d( nbpac, zat_i_ ac(1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) )288 CALL tab_2d_1d( nbpac, zat_i_1d (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 308 289 DO jl = 1, jpl 309 CALL tab_2d_1d( nbpac, za_i_ ac(1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) )310 CALL tab_2d_1d( nbpac, zv_i_ ac(1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) )311 CALL tab_2d_1d( nbpac, zoa_i_ ac(1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )312 CALL tab_2d_1d( nbpac, zsmv_i_ ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )290 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 291 CALL tab_2d_1d( nbpac, zv_i_1d (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 292 CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 293 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 313 294 DO jk = 1, nlay_i 314 CALL tab_2d_1d( nbpac, ze_i_ ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )295 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 315 296 END DO ! jk 316 297 END DO ! jl … … 322 303 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 323 304 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 324 CALL tab_2d_1d( nbpac, zvrel_ ac(1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )305 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 325 306 326 307 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) )308 CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac) , hfx_opw, jpi, jpj, npac(1:nbpac) ) 328 309 329 310 !------------------------------------------------------------------------------! … … 334 315 ! Keep old ice areas and volume in memory 335 316 !----------------------------------------- 336 zv_old(:,:) = zv_i_ ac(:,:)337 za_old(:,:) = za_i_ ac(:,:)317 zv_old(:,:) = zv_i_1d(:,:) 318 za_old(:,:) = za_i_1d(:,:) 338 319 339 320 !---------------------- … … 343 324 zh_newice(ji) = hiccrit(1) 344 325 END DO 345 IF( fraz_swi == 1.0 ) 326 IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:) 346 327 347 328 !---------------------- … … 370 351 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 371 352 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 372 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt) ) &353 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) ) & 373 354 & - 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 happen376 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) &377 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus378 355 END DO ! ji 379 356 !---------------- … … 405 382 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 406 383 ! Total heat flux used in this process [W.m-2] 407 hfx_ tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * zdE * r1_rdtice384 hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 408 385 ! mass flux 409 wfx_opw_1d(ji) = wfx_opw_1d(ji) +zv_newice(ji) * rhoic * r1_rdtice386 wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 410 387 ! salt flux 411 388 sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 412 389 413 390 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 414 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 391 zinda = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 392 zfrazb = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 415 393 zv_frazb(ji) = zfrazb * zv_newice(ji) 416 394 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 417 395 END DO 418 396 419 !------------------------------------420 ! Diags for energy conservation test421 !------------------------------------422 DO ji = 1, nbpac423 ii = MOD( npac(ji) - 1 , jpi ) + 1424 ij = ( npac(ji) - 1 ) / jpi + 1425 !426 zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)427 !zde = ze_newice(ji) * area(ii,ij) * zv_newice(ji)428 !429 ! clem: change that?430 vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji) ! volume431 et_i_init(ii,ij) = et_i_init(ii,ij) + zde ! Energy432 433 END DO434 397 435 398 !----------------- … … 438 401 DO ji = 1, nbpac 439 402 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 440 END DO !ji403 END DO 441 404 442 405 !------------------------------------------------------------------------------! … … 450 413 ! we keep the excessive volume in memory and attribute it later to bottom accretion 451 414 DO ji = 1, nbpac 452 IF ( za_newice(ji) > ( amax - zat_i_ ac(ji) ) ) THEN453 zda_res(ji) = za_newice(ji) - ( amax - zat_i_ ac(ji) )415 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN 416 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) ) 454 417 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 455 418 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 464 427 ! Laterally redistribute new ice volume and area 465 428 !------------------------------------------------ 466 zat_i_ ac(:) = 0._wp429 zat_i_1d(:) = 0._wp 467 430 DO jl = 1, jpl 468 431 DO ji = 1, nbpac 469 IF( hi_max (jl-1) < zh_newice(ji) .AND. & 470 & zh_newice(ji) <= hi_max (jl) ) THEN 471 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 472 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 473 zcatac (ji) = jl 432 IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 433 za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 434 zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 435 jcat (ji) = jl 474 436 ENDIF 475 zat_i_ ac(ji) = zat_i_ac(ji) + za_i_ac(ji,jl)437 zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d (ji,jl) 476 438 END DO 477 439 END DO … … 481 443 !---------------------------------- 482 444 DO ji = 1, nbpac 483 jl = zcatac(ji)! categroy in which new ice is put484 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10) ) ! 0 if old ice445 jl = jcat(ji) ! categroy in which new ice is put 446 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) ) ! 0 if old ice 485 447 END DO 486 448 487 449 DO jk = 1, nlay_i 488 450 DO ji = 1, nbpac 489 jl = zcatac(ji) 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) 451 jl = jcat(ji) 452 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 453 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 454 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) ) & 455 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 492 456 END DO 493 457 END DO … … 496 460 ! Add excessive volume of new ice at the bottom 497 461 !----------------------------------------------- 498 jm = 1 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 --- ! 511 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 462 DO jl = 1, jpl 463 h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 464 qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 465 512 466 DO jk = 1, nlay_i 513 467 DO ji = 1, nbpac 514 zthick0(ji,jk,jl) = zv_i_ac(ji,jl) / REAL( nlay_i )515 zqm0 (ji,jk,jl) = ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)468 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 469 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 516 470 END DO 517 471 END DO 518 END DO 519 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 472 520 473 DO ji = 1, nbpac 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) 524 END DO ! ji 525 END DO ! jl 526 527 ze_i_ac(:,:,:) = 0._wp 528 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 529 DO jk = 1, nlay_i 530 DO layer = 1, nlay_i + 1 531 DO ji = 1, nbpac 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 --- ! 552 DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 553 DO jk = 1, nlay_i 554 DO ji = 1, nbpac 555 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 556 ze_i_ac(ji,jk,jl) = zindb * ze_i_ac(ji,jk,jl) / MAX( zv_i_ac(ji,jl), epsi10 ) * REAL( nlay_i ) 557 END DO 558 END DO 559 END DO 560 474 zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 475 zv_newfra = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 476 za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl) 477 478 h_i_old (ji,nlay_i+1) = zv_newfra 479 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 480 481 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 482 ENDDO 483 484 ! --- Ice enthalpy remapping --- ! 485 CALL lim_thd_ent( 1, nbpac, jl, ze_i_1d(1:nbpac,:,jl) ) 486 487 ENDDO 561 488 562 489 !------------ … … 565 492 DO jl = 1, jpl 566 493 DO ji = 1, nbpac 567 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ ac(ji,jl) + epsi10 ) ) ! 0 if no ice and 1 if yes568 zoa_i_ ac(ji,jl) = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb494 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 495 zoa_i_1d(ji,jl) = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 569 496 END DO 570 497 END DO … … 576 503 DO jl = 1, jpl 577 504 DO ji = 1, nbpac 578 zdv = zv_i_ ac(ji,jl) - zv_old(ji,jl)579 zsmv_i_ ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji)505 zdv = zv_i_1d(ji,jl) - zv_old(ji,jl) 506 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 580 507 END DO 581 508 END DO … … 586 513 !------------------------------------------------------------------------------! 587 514 DO jl = 1, jpl 588 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ ac(1:nbpac,jl), jpi, jpj )589 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ ac(1:nbpac,jl), jpi, jpj )590 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_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 )515 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 516 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 517 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 518 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 592 519 DO jk = 1, nlay_i 593 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ ac(1:nbpac,jk,jl), jpi, jpj )520 CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 594 521 END DO 595 522 END DO … … 599 526 600 527 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 )528 CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 602 529 ! 603 530 ENDIF ! nbpac > 0 … … 612 539 ! heat content in Joules 613 540 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 541 END DO 616 542 END DO … … 618 544 END DO 619 545 620 !------------------------------------------------------------------------------|621 ! 10) Conservation check and changes in each ice category622 !------------------------------------------------------------------------------|623 IF( con_i ) THEN624 CALL lim_column_sum (jpl, v_i, vt_i_final)625 fieldid = 'v_i, limthd_lac'626 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)627 !628 CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)629 fieldid = 'e_i, limthd_lac'630 CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)631 !632 CALL lim_column_sum (jpl, v_s, vt_s_final)633 fieldid = 'v_s, limthd_lac'634 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)635 !636 ! CALL lim_column_sum (jpl, e_s(:,:,1,:) , et_s_init)637 ! fieldid = 'e_s, limthd_lac'638 ! CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)639 IF( ln_nicep ) THEN640 DO ji = mi0(jiindx), mi1(jiindx)641 DO jj = mj0(jjindx), mj1(jjindx)642 WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)643 WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)644 WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)645 WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)646 END DO647 END DO648 ENDIF649 !650 ENDIF651 546 ! 652 CALL wrk_dealloc( jpij, zcatac) ! integer547 CALL wrk_dealloc( jpij, jcat ) ! integer 653 548 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 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 ) 656 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac ) 657 CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 658 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 549 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 550 CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 551 CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d ) 552 CALL wrk_dealloc( jpi,jpj, zvrel ) 659 553 ! 660 554 END SUBROUTINE lim_thd_lac
Note: See TracChangeset
for help on using the changeset viewer.
