Changeset 5123 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
- Timestamp:
- 2015-03-04T17:06:03+01:00 (9 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4990 r5123 22 22 USE thd_ice ! LIM thermodynamics 23 23 USE dom_ice ! LIM domain 24 USE par_ice ! LIM parameters25 24 USE ice ! LIM variables 26 25 USE limtab ! LIM 2D <==> 1D … … 112 111 113 112 REAL(wp), POINTER, DIMENSION(:,:) :: zvrel ! relative ice / frazil velocity 113 114 REAL(wp) :: zcai = 1.4e-3_wp 114 115 !!-----------------------------------------------------------------------! 115 116 … … 129 130 DO ji = 1, jpi 130 131 !Energy of melting q(S,T) [J.m-3] 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 ) 134 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 132 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 ) ) !0 if no ice 133 e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp ) 135 134 END DO 136 135 END DO … … 155 154 156 155 ! Default new ice thickness 157 hicol(:,:) = hiccrit158 159 IF( fraz_swi == 1) THEN156 hicol(:,:) = rn_hnewice 157 158 IF( ln_frazil ) THEN 160 159 161 160 !-------------------- … … 166 165 zhicrit = 0.04 ! frazil ice thickness 167 166 ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 168 zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)167 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 169 168 zgamafr = 0.03 170 169 … … 176 175 !------------- 177 176 ! C-grid wind stress components 178 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj) &179 & + utau_ice(ji ,jj ) * tmu(ji ,jj) ) * 0.5_wp180 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &181 & + vtau_ice(ji ,jj ) * tmv(ji ,jj) ) * 0.5_wp177 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 178 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 179 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 180 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 182 181 ! Square root of wind stress 183 182 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) … … 195 194 ! C-grid ice velocity 196 195 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_wp198 zvgy = rswitch * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1) + v_ice(ji,jj) * tmv(ji,jj) ) * 0.5_wp196 zvgx = rswitch * ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 197 zvgy = rswitch * ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 199 198 200 199 !----------------------------------- … … 222 221 iterate_frazil = .true. 223 222 224 DO WHILE ( iter .LT.100 .AND. iterate_frazil )223 DO WHILE ( iter < 100 .AND. iterate_frazil ) 225 224 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 226 225 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 … … 320 319 !---------------------- 321 320 DO ji = 1, nbpac 322 zh_newice(ji) = hiccrit323 END DO 324 IF( fraz_swi == 1) zh_newice(1:nbpac) = hicol_1d(1:nbpac)321 zh_newice(ji) = rn_hnewice 322 END DO 323 IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 325 324 326 325 !---------------------- 327 326 ! Salinity of new ice 328 327 !---------------------- 329 SELECT CASE ( n um_sal )328 SELECT CASE ( nn_icesal ) 330 329 CASE ( 1 ) ! Sice = constant 331 zs_newice(1:nbpac) = bulk_sal330 zs_newice(1:nbpac) = rn_icesal 332 331 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 333 332 DO ji = 1, nbpac 334 333 ii = MOD( npac(ji) - 1 , jpi ) + 1 335 334 ij = ( npac(ji) - 1 ) / jpi + 1 336 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij) )335 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) ) 337 336 END DO 338 337 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 345 344 ! We assume that new ice is formed at the seawater freezing point 346 345 DO ji = 1, nbpac 347 ztmelts = - tmut * zs_newice(ji) + rt t! Melting point (K)346 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K) 348 347 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_1d(ji) ) & 349 & + lfus * ( 1.0 - ( ztmelts - rt t ) / MIN( t_bo_1d(ji) - rtt, -epsi10 ) ) &350 & - rcp * ( ztmelts - rt t) )351 END DO ! ji348 & + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) & 349 & - rcp * ( ztmelts - rt0 ) ) 350 END DO 352 351 353 352 !---------------- … … 363 362 DO ji = 1, nbpac 364 363 365 zEi = - ze_newice(ji) / rhoic! specific enthalpy of forming ice [J/kg]366 367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) 364 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] 365 366 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 368 367 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 369 368 … … 372 371 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 373 372 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 374 zv_newice(ji) = - zfmdt /rhoic373 zv_newice(ji) = - zfmdt * r1_rhoic 375 374 376 375 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux … … 387 386 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 388 387 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 *maxfrazb388 zfrazb = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 390 389 zv_frazb(ji) = zfrazb * zv_newice(ji) 391 390 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 409 408 ! we keep the excessive volume in memory and attribute it later to bottom accretion 410 409 DO ji = 1, nbpac 411 IF ( za_newice(ji) > ( amax - zat_i_1d(ji) ) ) THEN412 zda_res(ji) = za_newice(ji) - ( amax - zat_i_1d(ji) )410 IF ( za_newice(ji) > ( rn_amax - zat_i_1d(ji) ) ) THEN 411 zda_res(ji) = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 413 412 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 414 413 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 459 458 DO jk = 1, nlay_i 460 459 DO ji = 1, nbpac 461 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )460 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 462 461 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 463 462 END DO … … 525 524 DO jj = 1, jpj 526 525 DO ji = 1, jpi 527 ! heat content in J oules528 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 )526 ! heat content in J/m2 527 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 529 528 END DO 530 529 END DO
Note: See TracChangeset
for help on using the changeset viewer.
