- Timestamp:
- 2014-11-27T16:21:44+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4743_NOC2_ZTS/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4688 r4899 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 … … 71 72 !! - Computation of variation of ice volume and mass 72 73 !! - Computation of frldb after lateral accretion and 73 !! update ht_s_ b, ht_i_band tbif_1d(:,:)74 !! update ht_s_1d, ht_i_1d and tbif_1d(:,:) 74 75 !!------------------------------------------------------------------------ 75 INTEGER :: ji,jj,jk,jl ,jm! dummy loop indices76 INTEGER :: layer, nbpac! local integers77 INTEGER :: ii, ij, iter ! - -76 INTEGER :: ji,jj,jk,jl ! dummy loop indices 77 INTEGER :: nbpac ! local integers 78 INTEGER :: ii, ij, iter ! - - 78 79 REAL(wp) :: ztmelts, zdv, zfrazb, zweight, zindb, zinda, zde ! local scalars 79 80 REAL(wp) :: zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new ! - - … … 89 90 REAL(wp) :: zv_newfra 90 91 91 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows92 INTEGER , POINTER, DIMENSION(:) :: jcat ! indexes of categories where new ice grows 92 93 REAL(wp), POINTER, DIMENSION(:) :: zswinew ! switch for new ice or not 93 94 … … 101 102 REAL(wp), POINTER, DIMENSION(:) :: zda_res ! residual area in case of excessive heat budget 102 103 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 104 REAL(wp), POINTER, DIMENSION(:) :: zv_frazb ! accretion of frazil ice at the ice bottom 105 105 REAL(wp), POINTER, DIMENSION(:) :: zvrel_1d ! relative ice / frazil velocity (1D vector) 106 106 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 113 114 REAL(wp), POINTER, DIMENSION(:,:,:) :: ze_i_1d 107 REAL(wp), POINTER, DIMENSION(:,:) :: zv_b ! old volume of ice in category jl 108 REAL(wp), POINTER, DIMENSION(:,:) :: za_b ! old area of ice in category jl 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 115 116 116 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity … … 119 119 CALL wrk_alloc( jpij, jcat ) ! integer 120 120 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 )121 CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 122 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 123 CALL wrk_alloc( jpij,nlay_i+1,jpl, ze_i_1d ) 124 124 CALL wrk_alloc( jpi,jpj, zvrel ) 125 125 … … 133 133 !Energy of melting q(S,T) [J.m-3] 134 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 )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, wp ) 136 136 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 137 137 END DO … … 171 171 zgamafr = 0.03 172 172 173 DO jj = 1, jpj 174 DO ji = 1, jpi 175 173 DO jj = 2, jpj 174 DO ji = 2, jpi 176 175 IF ( qlead(ji,jj) < 0._wp ) THEN 177 176 !------------- … … 243 242 END DO ! loop on ji ends 244 243 END DO ! loop on jj ends 244 ! 245 CALL lbc_lnk( zvrel(:,:), 'T', 1. ) 246 CALL lbc_lnk( hicol(:,:), 'T', 1. ) 245 247 246 248 ENDIF ! End of computation of frazil ice collection thickness … … 255 257 ! This occurs if open water energy budget is negative 256 258 nbpac = 0 259 npac(:) = 0 260 ! 257 261 DO jj = 1, jpj 258 262 DO ji = 1, jpi … … 298 302 299 303 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) )304 CALL tab_2d_1d( nbpac, t_bo_1d (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 301 305 CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac) , sfx_opw, jpi, jpj, npac(1:nbpac) ) 302 306 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) ) 307 CALL tab_2d_1d( nbpac, hicol_1d (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 305 308 CALL tab_2d_1d( nbpac, zvrel_1d (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 306 309 … … 315 318 ! Keep old ice areas and volume in memory 316 319 !----------------------------------------- 317 zv_old(:,:) = zv_i_1d(:,:) 318 za_old(:,:) = za_i_1d(:,:) 319 320 zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 321 za_b(1:nbpac,:) = za_i_1d(1:nbpac,:) 320 322 !---------------------- 321 323 ! Thickness of new ice … … 324 326 zh_newice(ji) = hiccrit 325 327 END DO 326 IF( fraz_swi == 1 ) zh_newice( :) = hicol_b(:)328 IF( fraz_swi == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 327 329 328 330 !---------------------- … … 331 333 SELECT CASE ( num_sal ) 332 334 CASE ( 1 ) ! Sice = constant 333 zs_newice( :) = bulk_sal335 zs_newice(1:nbpac) = bulk_sal 334 336 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 335 337 DO ji = 1, nbpac … … 339 341 END DO 340 342 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 341 zs_newice( :) = 2.3343 zs_newice(1:nbpac) = 2.3 342 344 END SELECT 343 345 … … 348 350 DO ji = 1, nbpac 349 351 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 ) ) &352 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 353 & + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) & 352 354 & - rcp * ( ztmelts - rtt ) ) 353 355 END DO ! ji … … 367 369 zEi = - ze_newice(ji) / rhoic ! specific enthalpy of forming ice [J/kg] 368 370 369 zEw = rcp * ( t_bo_ b(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_b[J/kg]371 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 370 372 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 371 373 … … 438 440 DO ji = 1, nbpac 439 441 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 ice442 zswinew (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) ) ! 0 if old ice 441 443 END DO 442 444 … … 446 448 zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 447 449 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) ) &450 & ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) ) & 449 451 & * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 450 452 END DO … … 472 474 za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl) 473 475 zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 474 475 476 ! for remapping 476 477 h_i_old (ji,nlay_i+1) = zv_newfra … … 479 480 480 481 ! --- 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 482 CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 485 483 ENDDO 486 484 … … 491 489 DO ji = 1, nbpac 492 490 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes 493 zoa_i_1d(ji,jl) = za_ old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb491 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb 494 492 END DO 495 493 END DO … … 500 498 DO jl = 1, jpl 501 499 DO ji = 1, nbpac 502 zdv = zv_i_1d(ji,jl) - zv_ old(ji,jl)500 zdv = zv_i_1d(ji,jl) - zv_b(ji,jl) 503 501 zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 504 502 END DO … … 519 517 CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj ) 520 518 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 519 523 520 CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) … … 534 531 DO ji = 1, jpi 535 532 ! 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 )533 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 534 END DO 538 535 END DO … … 543 540 CALL wrk_dealloc( jpij, jcat ) ! integer 544 541 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 )542 CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d ) 543 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 544 CALL wrk_dealloc( jpij,nlay_i+1,jpl, ze_i_1d ) 548 545 CALL wrk_dealloc( jpi,jpj, zvrel ) 549 546 !
Note: See TracChangeset
for help on using the changeset viewer.
