- Timestamp:
- 2015-07-16T13:55:15+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5003_MERCATOR6_CRS/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r4990 r5602 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 … … 32 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 33 32 USE limthd_ent 33 USE limvar 34 34 35 35 IMPLICIT NONE … … 106 106 REAL(wp), POINTER, DIMENSION(:,:) :: za_i_1d ! 1-D version of a_i 107 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_i109 108 REAL(wp), POINTER, DIMENSION(:,:) :: zsmv_i_1d ! 1-D version of smv_i 110 109 … … 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 … … 117 118 CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 118 119 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, z oa_i_1d, zsmv_i_1d )120 CALL wrk_alloc( jpij,nlay_i +1,jpl, ze_i_1d )120 CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 121 CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d ) 121 122 CALL wrk_alloc( jpi,jpj, zvrel ) 122 123 124 CALL lim_var_agg(1) 125 CALL lim_var_glo2eqv 123 126 !------------------------------------------------------------------------------| 124 127 ! 2) Convert units for ice internal energy … … 129 132 DO ji = 1, jpi 130 133 !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 134 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl) - epsi20 ) ) !0 if no ice 135 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 136 END DO 136 137 END DO … … 155 156 156 157 ! Default new ice thickness 157 hicol(:,:) = hiccrit158 159 IF( fraz_swi == 1) THEN158 hicol(:,:) = rn_hnewice 159 160 IF( ln_frazil ) THEN 160 161 161 162 !-------------------- … … 166 167 zhicrit = 0.04 ! frazil ice thickness 167 168 ztwogp = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav 168 zsqcd = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)169 zsqcd = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag) 169 170 zgamafr = 0.03 170 171 … … 176 177 !------------- 177 178 ! 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_wp179 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) & 180 & + utau_ice(ji ,jj ) * umask(ji ,jj ,1) ) * 0.5_wp 181 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) & 182 & + vtau_ice(ji ,jj ) * vmask(ji ,jj ,1) ) * 0.5_wp 182 183 ! Square root of wind stress 183 184 ztenagm = SQRT( SQRT( ztaux**2 + ztauy**2 ) ) … … 195 196 ! C-grid ice velocity 196 197 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_wp198 zvgx = rswitch * ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 199 zvgy = rswitch * ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 199 200 200 201 !----------------------------------- … … 222 223 iterate_frazil = .true. 223 224 224 DO WHILE ( iter .LT.100 .AND. iterate_frazil )225 DO WHILE ( iter < 100 .AND. iterate_frazil ) 225 226 zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) & 226 227 - hicol(ji,jj) * zhicrit * ztwogp * zvrel2 … … 266 267 ! debug point to follow 267 268 jiindex_1d = 0 268 IF( ln_ nicep) THEN269 DO ji = mi0( jiindx), mi1(jiindx)270 DO jj = mj0(j jindx), mj1(jjindx)269 IF( ln_icectl ) THEN 270 DO ji = mi0(iiceprt), mi1(iiceprt) 271 DO jj = mj0(jiceprt), mj1(jiceprt) 271 272 IF ( qlead(ji,jj) < 0._wp ) THEN 272 273 jiindex_1d = (jj - 1) * jpi + ji … … 276 277 ENDIF 277 278 278 IF( ln_ nicep) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac279 IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 279 280 280 281 !------------------------------ … … 290 291 CALL tab_2d_1d( nbpac, za_i_1d (1:nbpac,jl), a_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 291 292 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 293 CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 294 294 DO jk = 1, nlay_i 295 295 CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 296 END DO ! jk297 END DO ! jl296 END DO 297 END DO 298 298 299 299 CALL tab_2d_1d( nbpac, qlead_1d (1:nbpac) , qlead , jpi, jpj, npac(1:nbpac) ) … … 320 320 !---------------------- 321 321 DO ji = 1, nbpac 322 zh_newice(ji) = hiccrit323 END DO 324 IF( fraz_swi == 1) zh_newice(1:nbpac) = hicol_1d(1:nbpac)322 zh_newice(ji) = rn_hnewice 323 END DO 324 IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac) 325 325 326 326 !---------------------- 327 327 ! Salinity of new ice 328 328 !---------------------- 329 SELECT CASE ( n um_sal )329 SELECT CASE ( nn_icesal ) 330 330 CASE ( 1 ) ! Sice = constant 331 zs_newice(1:nbpac) = bulk_sal331 zs_newice(1:nbpac) = rn_icesal 332 332 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 333 333 DO ji = 1, nbpac 334 334 ii = MOD( npac(ji) - 1 , jpi ) + 1 335 335 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) )336 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij) ) 337 337 END DO 338 338 CASE ( 3 ) ! Sice = F(z) [multiyear ice] … … 345 345 ! We assume that new ice is formed at the seawater freezing point 346 346 DO ji = 1, nbpac 347 ztmelts = - tmut * zs_newice(ji) + rt t! Melting point (K)347 ztmelts = - tmut * zs_newice(ji) + rt0 ! Melting point (K) 348 348 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 ! ji349 & + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) ) & 350 & - rcp * ( ztmelts - rt0 ) ) 351 END DO 352 352 353 353 !---------------- … … 356 356 DO ji = 1, nbpac 357 357 zo_newice(ji) = 0._wp 358 END DO ! ji358 END DO 359 359 360 360 !------------------- … … 363 363 DO ji = 1, nbpac 364 364 365 zEi = - ze_newice(ji) / rhoic! specific enthalpy of forming ice [J/kg]366 367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) 365 zEi = - ze_newice(ji) * r1_rhoic ! specific enthalpy of forming ice [J/kg] 366 367 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! specific enthalpy of seawater at t_bo_1d [J/kg] 368 368 ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 369 369 … … 372 372 zfmdt = - qlead_1d(ji) / zdE ! Fm.dt [kg/m2] (<0) 373 373 ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point 374 zv_newice(ji) = - zfmdt /rhoic374 zv_newice(ji) = - zfmdt * r1_rhoic 375 375 376 376 zQm = zfmdt * zEw ! heat to the ocean >0 associated with mass flux … … 387 387 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 388 388 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 *maxfrazb389 zfrazb = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb 390 390 zv_frazb(ji) = zfrazb * zv_newice(ji) 391 391 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) … … 409 409 ! we keep the excessive volume in memory and attribute it later to bottom accretion 410 410 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) )411 IF ( za_newice(ji) > ( rn_amax - zat_i_1d(ji) ) ) THEN 412 zda_res(ji) = za_newice(ji) - ( rn_amax - zat_i_1d(ji) ) 413 413 zdv_res(ji) = zda_res (ji) * zh_newice(ji) 414 414 za_newice(ji) = za_newice(ji) - zda_res (ji) … … 459 459 DO jk = 1, nlay_i 460 460 DO ji = 1, nbpac 461 h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )461 h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i 462 462 qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 463 463 END DO … … 478 478 ENDDO 479 479 480 !------------481 ! Update age482 !------------483 DO jl = 1, jpl484 DO ji = 1, nbpac485 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) ) ! 0 if no ice and 1 if yes486 zoa_i_1d(ji,jl) = za_b(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * rswitch487 END DO488 END DO489 490 480 !----------------- 491 481 ! Update salinity … … 504 494 CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 505 495 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 506 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj )507 496 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 508 497 DO jk = 1, nlay_i … … 525 514 DO jj = 1, jpj 526 515 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 )516 ! heat content in J/m2 517 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 529 518 END DO 530 519 END DO … … 536 525 CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 537 526 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, z oa_i_1d, zsmv_i_1d )539 CALL wrk_dealloc( jpij,nlay_i +1,jpl, ze_i_1d )527 CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d ) 528 CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d ) 540 529 CALL wrk_dealloc( jpi,jpj, zvrel ) 541 530 !
Note: See TracChangeset
for help on using the changeset viewer.