- Timestamp:
- 2012-11-02T07:13:40+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3523 r3524 13 13 !! 'key_lim3' LIM3 sea-ice model 14 14 !!---------------------------------------------------------------------- 15 !! lim_lat_acr : lateral accretion of ice 16 !!---------------------------------------------------------------------- 17 USE par_oce ! ocean parameters 18 USE dom_oce ! domain variables 19 USE phycst ! physical constants 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE sbc_ice ! Surface boundary condition: ice fields 22 USE thd_ice ! LIM thermodynamics 23 USE dom_ice ! LIM domain 24 USE par_ice ! LIM parameters 25 USE ice ! LIM variables 26 USE limtab ! LIM 2D <==> 1D 27 USE limcons ! LIM conservation 28 USE in_out_manager ! I/O manager 29 USE lib_mpp ! MPP library 30 USE wrk_nemo ! work arrays 15 !! lim_lat_acr : lateral accretion of ice 16 !!---------------------------------------------------------------------- 17 USE par_oce ! ocean parameters 18 USE dom_oce ! domain variables 19 USE phycst ! physical constants 20 USE sbc_oce ! Surface boundary condition: ocean fields 21 USE sbc_ice ! Surface boundary condition: ice fields 22 USE thd_ice ! LIM thermodynamics 23 USE dom_ice ! LIM domain 24 USE par_ice ! LIM parameters 25 USE ice ! LIM variables 26 USE limtab ! LIM 2D <==> 1D 27 USE limcons ! LIM conservation 28 USE in_out_manager ! I/O manager 29 USE lib_mpp ! MPP library 30 USE wrk_nemo ! work arrays 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 32 32 33 IMPLICIT NONE … … 45 46 46 47 !!---------------------------------------------------------------------- 47 !! NEMO/LIM3 4.0, UCL - NEMO Consortium (2011)48 !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 48 49 !! $Id$ 49 50 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 77 78 !! update ht_s_b, ht_i_b and tbif_1d(:,:) 78 79 !!------------------------------------------------------------------------ 79 INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices80 INTEGER :: layer, nbpac ! local integers81 INTEGER :: zji, zjj, iter ! - -82 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde! local scalars80 INTEGER :: ji,jj,jk,jl,jm ! dummy loop indices 81 INTEGER :: layer, nbpac ! local integers 82 INTEGER :: zji, zjj, iter ! - - 83 REAL(wp) :: ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde ! local scalars 83 84 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 84 85 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - 86 REAL(wp) :: zcoef ! - - 85 87 LOGICAL :: iterate_frazil ! iterate frazil ice collection thickness 86 88 CHARACTER (len = 15) :: fieldid … … 205 207 !------------- 206 208 ! C-grid wind stress components 207 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) &208 & + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) / 2.0209 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &210 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) / 2.0209 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) & 210 & + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) * 0.5_wp 211 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) & 212 & + vtau_ice(ji ,jj ) * tmv(ji ,jj ) ) * 0.5_wp 211 213 ! Square root of wind stress 212 214 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) … … 222 224 !------------------- 223 225 ! C-grid ice velocity 224 zindb = MAX( 0.0, SIGN(1.0, at_i(ji,jj) ))225 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj )&226 + u_ice(ji,jj ) * tmu(ji ,jj ) ) / 2.0227 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1)&228 + v_ice(ji,jj ) * tmv(ji ,jj ) ) / 2.0226 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 227 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 228 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 229 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 230 & + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp 229 231 230 232 !----------------------------------- … … 232 234 !----------------------------------- 233 235 ! absolute relative velocity 234 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + & 235 ( zvfry - zvgy ) * ( zvfry - zvgy ) & 236 , 0.15 * 0.15 ) 237 zvrel(ji,jj) = SQRT(zvrel2) 236 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) & 237 & + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 ) 238 zvrel(ji,jj) = SQRT( zvrel2 ) 238 239 239 240 !--------------------- … … 241 242 !--------------------- 242 243 hicol(ji,jj) = zhicrit + 0.1 243 hicol(ji,jj) = zhicrit + hicol(ji,jj) / & 244 ( hicol(ji,jj) * hicol(ji,jj) - & 245 zhicrit * zhicrit ) * ztwogp * zvrel2 244 hicol(ji,jj) = zhicrit + hicol(ji,jj) & 245 & / ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) * ztwogp * zvrel2 246 247 !!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1) 248 !!gm = zhicrit**2 + 0.2*zhicrit +0.01 249 !!gm therefore the 2 lines with hicol can be replaced by 1 line: 250 !!gm hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2 251 !!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop 246 252 247 253 iter = 1 … … 278 284 DO jj = 1, jpj 279 285 DO ji = 1, jpi 280 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0) THEN286 IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0._wp ) THEN 281 287 nbpac = nbpac + 1 282 288 npac( nbpac ) = (jj - 1) * jpi + ji 283 IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 284 jiindex_1d = nbpac 285 ENDIF 289 IF( ji == jiindx .AND. jj == jjindx ) jiindex_1d = nbpac 286 290 ENDIF 287 291 END DO 288 292 END DO 289 293 290 IF( ln_nicep ) THEN 291 WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 292 ENDIF 294 IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 293 295 294 296 !------------------------------ … … 302 304 CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 303 305 DO jl = 1, jpl 304 CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl) , a_i(:,:,jl), jpi, jpj, npac(1:nbpac) )305 CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl) , v_i(:,:,jl), jpi, jpj, npac(1:nbpac) )306 CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl) , oa_i(:,:,jl), jpi, jpj, npac(1:nbpac) )306 CALL tab_2d_1d( nbpac, za_i_ac (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, zv_i_ac (1:nbpac,jl), v_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 308 CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 307 309 CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 308 310 DO jk = 1, nlay_i … … 311 313 END DO ! jl 312 314 313 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) )314 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) )315 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )316 CALL tab_2d_1d( nbpac, fseqv_1d (1:nbpac) , fseqv, jpi, jpj, npac(1:nbpac) )315 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) 316 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) 317 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 318 CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac) , sfx_thd, jpi, jpj, npac(1:nbpac) ) 317 319 CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, jpi, jpj, npac(1:nbpac) ) 318 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) )319 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) )320 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 321 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 320 322 321 323 !------------------------------------------------------------------------------! … … 392 394 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 393 395 DO ji = 1, nbpac 394 fseqv_1d (ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice396 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 395 397 rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * zv_newice(ji) 396 398 END DO ! ji … … 400 402 !------------------------------------ 401 403 DO ji = 1, nbpac 402 ! Volume403 zj i = MOD( npac(ji) - 1, jpi )+ 1404 zjj = ( npac(ji) - 1 ) / jpi + 1405 vt_i_init(zji,zjj) = vt_i_init(zji,zjj) +zv_newice(ji)406 ! Energy407 zde = ze_newice(ji) / unit_fac408 zde = zde * area(zji,zjj) * zv_newice(ji)409 et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde 404 zji = MOD( npac(ji) - 1 , jpi ) + 1 405 zjj = ( npac(ji) - 1 ) / jpi + 1 406 ! 407 zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji) 408 ! 409 vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji) ! volume 410 et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde ! Energy 411 410 412 END DO 411 413 … … 417 419 !----------------- 418 420 DO ji = 1, nbpac 419 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 420 ! diagnostic 421 zji = MOD( npac(ji) - 1, jpi ) + 1 422 zjj = ( npac(ji) - 1 ) / jpi + 1 421 zji = MOD( npac(ji) - 1 , jpi ) + 1 422 zjj = ( npac(ji) - 1 ) / jpi + 1 423 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 423 424 diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 424 425 END DO !ji … … 438 439 !------------------------------------------- 439 440 ! If lateral ice growth gives an ice concentration gt 1, then 440 ! we keep the excessive volume in memory and attribute it later 441 ! to bottom accretion 442 DO ji = 1, nbpac 443 ! vectorize 444 IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN 445 zda_res(ji) = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 446 zdv_res(ji) = zda_res(ji) * zh_newice(ji) 447 za_newice(ji) = za_newice(ji) - zda_res(ji) 448 zv_newice(ji) = zv_newice(ji) - zdv_res(ji) 441 ! we keep the excessive volume in memory and attribute it later to bottom accretion 442 DO ji = 1, nbpac 443 IF ( za_newice(ji) > ( 1._wp - zat_i_ac(ji) ) ) THEN 444 zda_res(ji) = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 445 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 446 za_newice(ji) = za_newice(ji) - zda_res (ji) 447 zv_newice(ji) = zv_newice(ji) - zdv_res (ji) 449 448 ELSE 450 zda_res(ji) = 0. 0451 zdv_res(ji) = 0. 0449 zda_res(ji) = 0._wp 450 zdv_res(ji) = 0._wp 452 451 ENDIF 453 452 END DO ! ji … … 459 458 DO jl = 1, jpl 460 459 DO ji = 1, nbpac 461 IF( hi_max (jl-1) < zh_newice(ji) .AND. &462 & zh_newice(ji) <= hi_max (jl) ) THEN460 IF( hi_max (jl-1) < zh_newice(ji) .AND. & 461 & zh_newice(ji) <= hi_max (jl) ) THEN 463 462 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 464 463 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) … … 466 465 zcatac (ji) = jl 467 466 ENDIF 468 END DO ! ji469 END DO ! jl467 END DO 468 END DO 470 469 471 470 !---------------------------------- … … 483 482 DO ji = 1, nbpac 484 483 jl = zcatac(ji) 485 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]484 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 486 485 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i , zh_newice(ji) ) & 487 486 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) ) … … 489 488 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i & 490 489 + za_newice(ji) * ze_newice(ji) * zalphai & 491 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( ( zv_i_ac(ji,jl)) / nlay_i )490 + za_newice(ji) * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i ) 492 491 END DO 493 492 END DO … … 529 528 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & 530 529 & + zindb * zdh_frazb(ji) ! frazil ice may coalesce 531 zdummy(ji,jl) = zv_i_ac(ji,jl) /MAX(za_i_ac(ji,jl),epsi10)*zindb ! thickness of residual ice530 zdummy(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb ! thickness of residual ice 532 531 END DO 533 532 END DO … … 613 612 END DO 614 613 END DO 615 CALL tab_1d_2d( nbpac, fseqv , npac(1:nbpac), fseqv_1d(1:nbpac), jpi, jpj )614 CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj ) 616 615 CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 617 616 ! … … 623 622 DO jl = 1, jpl 624 623 DO jk = 1, nlay_i ! heat content in 10^9 Joules 625 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i 624 e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac 626 625 END DO 627 626 END DO
Note: See TracChangeset
for help on using the changeset viewer.