Changeset 5208 for branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2015-04-13T15:08:59+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO11_restart_functionality/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4688 r5208 29 29 USE lib_mpp ! MPP library 30 30 USE wrk_nemo ! work arrays 31 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 31 32 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 32 33 USE limthd_ent … … 36 37 37 38 PUBLIC lim_thd_lac ! called by lim_thd 38 39 REAL(wp) :: epsi10 = 1.e-10_wp !40 REAL(wp) :: epsi20 = 1.e-20_wp !41 39 42 40 !!---------------------------------------------------------------------- … … 71 69 !! - Computation of variation of ice volume and mass 72 70 !! - Computation of frldb after lateral accretion and 73 !! update ht_s_ b, ht_i_band tbif_1d(:,:)71 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 74 72 !!------------------------------------------------------------------------ 75 INTEGER :: ji,jj,jk,jl ,jm! dummy loop indices76 INTEGER :: layer, nbpac! local integers77 INTEGER :: ii, ij, iter ! - -78 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, z indb, zinda, zde! local scalars73 INTEGER :: ji,jj,jk,jl ! dummy loop indices 74 INTEGER :: nbpac ! local integers 75 INTEGER :: ii, ij, iter ! - - 76 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zde ! local scalars 79 77 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - 80 78 REAL(wp) :: ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit ! - - … … 89 87 REAL(wp) :: zv_newfra 90 88 91 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows89 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 92 90 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 93 91 … … 101 99 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 102 100 REAL(wp), POINTER, DIMENSION(:) :: zat_i_1d ! total ice fraction 103 REAL(wp), POINTER, DIMENSION(:) :: zat_i_lev ! total ice fraction for level ice only (type 1) 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 101 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 105 102 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 103 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_ old! old volume of ice in category jl108 REAL(wp), POINTER, DIMENSION(:,:) :: za_ old! old area of ice in category jl109 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d 110 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d 111 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d 112 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d 104 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 105 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl 106 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_i_1d ! 1-D version of v_i 108 REAL(wp), POINTER, DIMENSION(:,:) :: zoa_i_1d ! 1-D version of oa_i 109 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 113 110 114 111 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d !: 1-D version of e_i … … 119 116 CALL wrk_alloc( jpij, jcat ) ! integer 120 117 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, z at_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 )118 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 119 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 120 CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 124 121 CALL wrk_alloc( jpi,jpj, zvrel ) 125 122 … … 132 129 DO ji = 1, jpi 133 130 !Energy of melting q(S,T) [J.m-3] 134 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 135 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 ) 131 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes 132 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) & 133 & / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i, wp ) 136 134 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 137 135 END DO … … 171 169 zgamafr = 0.03 172 170 173 DO jj = 1, jpj 174 DO ji = 1, jpi 175 171 DO jj = 2, jpj 172 DO ji = 2, jpi 176 173 IF ( qlead(ji,jj) < 0._wp ) THEN 177 174 !------------- … … 189 186 ! Frazil ice velocity 190 187 !--------------------- 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 )188 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 189 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 190 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 194 191 195 192 !------------------- … … 197 194 !------------------- 198 195 ! C-grid ice velocity 199 zindb = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 200 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) & 201 & + u_ice(ji,jj ) * tmu(ji ,jj ) ) * 0.5_wp 202 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) & 203 & + v_ice(ji,jj ) * tmv(ji ,jj ) ) * 0.5_wp 196 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) ) ) 197 zvgx = rswitch * ( u_ice(ji-1,jj ) * tmu(ji-1,jj ) + u_ice(ji,jj) * tmu(ji,jj) ) * 0.5_wp 198 zvgy = rswitch * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp 204 199 205 200 !----------------------------------- … … 243 238 END DO ! loop on ji ends 244 239 END DO ! loop on jj ends 240 ! 241 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 242 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 245 243 246 244 ENDIF ! End of computation of frazil ice collection thickness … … 255 253 ! This occurs if open water energy budget is negative 256 254 nbpac = 0 255 npac(:) = 0 256 ! 257 257 DO jj = 1, jpj 258 258 DO ji = 1, jpi … … 298 298 299 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) 300 CALL tab_2d_1d( nbpac, t_bo_ b(1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) )300 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 301 301 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 302 302 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac) , wfx_opw, jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 303 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 305 304 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 306 305 … … 315 314 ! Keep old ice areas and volume in memory 316 315 !----------------------------------------- 317 zv_old(:,:) = zv_i_1d(:,:) 318 za_old(:,:) = za_i_1d(:,:) 319 316 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 317 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 320 318 !---------------------- 321 319 ! Thickness of new ice … … 324 322 zh_newice(ji) = hiccrit 325 323 END DO 326 IF( fraz_swi == 1 ) zh_newice( :) = hicol_b(:)324 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 327 325 328 326 !---------------------- … … 331 329 SELECT CASE ( num_sal ) 332 330 CASE ( 1 ) ! Sice = constant 333 zs_newice( :) = bulk_sal331 zs_newice(1:nbpac) = bulk_sal 334 332 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 335 333 DO ji = 1, nbpac … … 339 337 END DO 340 338 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 341 zs_newice( :) = 2.3339 zs_newice(1:nbpac) = 2.3 342 340 END SELECT 343 341 … … 348 346 DO ji = 1, nbpac 349 347 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 350 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_ b(ji) ) &351 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_ b(ji) - rtt, -epsi10 ) ) &348 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 349 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) & 352 350 & - rcp * ( ztmelts - rtt ) ) 353 351 END DO ! ji … … 367 365 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 368 366 369 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b[J/kg]367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 370 368 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 371 369 … … 388 386 389 387 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 390 zinda= 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )391 zfrazb = zinda* ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb388 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 389 zfrazb = rswitch * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 392 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 393 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 438 436 DO ji = 1, nbpac 439 437 jl = jcat(ji) ! categroy in which new ice is put 440 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_ old(ji,jl) ) ) ! 0 if old ice438 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 441 439 END DO 442 440 … … 444 442 DO ji = 1, nbpac 445 443 jl = jcat(ji) 446 zinda= MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )444 rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 447 445 ze_i_1d(ji,jk,jl) = zswinew(ji) * ze_newice(ji) + & 448 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_ old(ji,jl) ) &449 & * zinda/ MAX( zv_i_1d(ji,jl), epsi20 )446 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 447 & * rswitch / MAX( zv_i_1d(ji,jl), epsi20 ) 450 448 END DO 451 449 END DO … … 468 466 ! new volumes including lateral/bottom accretion + residual 469 467 DO ji = 1, nbpac 470 zinda= MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )471 zv_newfra = zinda* ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )472 za_i_1d(ji,jl) = zinda* za_i_1d(ji,jl)468 rswitch = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 469 zv_newfra = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 470 za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl) 473 471 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 474 475 472 ! for remapping 476 473 h_i_old (ji,nlay_i+1) = zv_newfra 477 474 qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 478 475 ENDDO 479 480 476 ! --- Ice enthalpy remapping --- ! 481 IF( zv_newfra > 0._wp ) THEN 482 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 483 ENDIF 484 477 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 485 478 ENDDO 486 479 … … 490 483 DO jl = 1, jpl 491 484 DO ji = 1, nbpac 492 zindb= 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes493 zoa_i_1d(ji,jl) = za_ old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb485 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 486 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch 494 487 END DO 495 488 END DO … … 500 493 DO jl = 1, jpl 501 494 DO ji = 1, nbpac 502 zdv = zv_i_1d(ji,jl) - zv_ old(ji,jl)495 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 503 496 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 504 497 END DO … … 519 512 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 520 513 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj ) 521 CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )522 514 523 515 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) … … 534 526 DO ji = 1, jpi 535 527 ! heat content in Joules 536 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 )528 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 537 529 END DO 538 530 END DO … … 543 535 CALL wrk_dealloc( jpij, jcat ) ! integer 544 536 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 545 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, z at_i_lev, zv_frazb, zvrel_1d )546 CALL wrk_dealloc( jpij,jpl, zv_ old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d )547 CALL wrk_dealloc( jpij, jkmax,jpl, ze_i_1d )537 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 538 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 539 CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 548 540 CALL wrk_dealloc( jpi,jpj, zvrel ) 549 541 !
Note: See TracChangeset
for help on using the changeset viewer.