- Timestamp:
- 2012-11-21T14:19:18+01:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_NOC_2012_rev3555/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3294 r3625 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 … … 143 145 ! 1) Conservation check and changes in each ice category 144 146 !------------------------------------------------------------------------------! 145 IF 146 CALL lim_column_sum (jpl, v_i, vt_i_init)147 CALL lim_column_sum (jpl, v_s, vt_s_init)148 CALL lim_column_sum_energy ( jpl, nlay_i, e_i, et_i_init)149 CALL lim_column_sum (jpl,e_s(:,:,1,:) , et_s_init)147 IF( con_i ) THEN 148 CALL lim_column_sum ( jpl, v_i , vt_i_init) 149 CALL lim_column_sum ( jpl, v_s , vt_s_init) 150 CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 151 CALL lim_column_sum ( jpl, e_s(:,:,1,:) , et_s_init) 150 152 ENDIF 151 153 … … 158 160 DO ji = 1, jpi 159 161 !Energy of melting q(S,T) [J.m-3] 160 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / & 161 MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * & 162 nlay_i 163 zindb = 1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl))) !0 if no ice and 1 if yes 164 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl)*unit_fac*zindb 162 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * nlay_i 163 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 164 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 165 165 END DO 166 166 END DO … … 182 182 ! 183 183 184 zvrel(:,:) = 0. 0184 zvrel(:,:) = 0._wp 185 185 186 186 ! Default new ice thickness 187 DO jj = 1, jpj 188 DO ji = 1, jpi 189 hicol(ji,jj) = hiccrit(1) 190 END DO 191 END DO 192 193 IF (fraz_swi.eq.1.0) THEN 187 hicol(:,:) = hiccrit(1) 188 189 IF( fraz_swi == 1._wp ) THEN 194 190 195 191 !-------------------- 196 192 ! Physical constants 197 193 !-------------------- 198 hicol(:,:) = 0. 0194 hicol(:,:) = 0._wp 199 195 200 196 zhicrit = 0.04 ! frazil ice thickness … … 211 207 !------------- 212 208 ! C-grid wind stress components 213 ztaux = ( utau_ice(ji-1,jj ) * tmu(ji-1,jj ) &214 & + utau_ice(ji ,jj ) * tmu(ji ,jj ) ) / 2.0215 ztauy = ( vtau_ice(ji ,jj-1) * tmv(ji ,jj-1) &216 & + 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 217 213 ! Square root of wind stress 218 214 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) … … 228 224 !------------------- 229 225 ! C-grid ice velocity 230 zindb = MAX( 0.0, SIGN(1.0, at_i(ji,jj) ))231 zvgx = zindb * ( u_ice(ji-1,jj ) * tmu(ji-1,jj )&232 + u_ice(ji,jj ) * tmu(ji ,jj ) ) / 2.0233 zvgy = zindb * ( v_ice(ji ,jj-1) * tmv(ji ,jj-1)&234 + 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 235 231 236 232 !----------------------------------- … … 238 234 !----------------------------------- 239 235 ! absolute relative velocity 240 zvrel2 = MAX( ( zvfrx - zvgx ) * ( zvfrx - zvgx ) + & 241 ( zvfry - zvgy ) * ( zvfry - zvgy ) & 242 , 0.15 * 0.15 ) 243 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 ) 244 239 245 240 !--------------------- … … 247 242 !--------------------- 248 243 hicol(ji,jj) = zhicrit + 0.1 249 hicol(ji,jj) = zhicrit + hicol(ji,jj) / & 250 ( hicol(ji,jj) * hicol(ji,jj) - & 251 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 252 252 253 253 iter = 1 … … 284 284 DO jj = 1, jpj 285 285 DO ji = 1, jpi 286 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 287 287 nbpac = nbpac + 1 288 288 npac( nbpac ) = (jj - 1) * jpi + ji 289 IF ( (ji.eq.jiindx).AND.(jj.eq.jjindx) ) THEN 290 jiindex_1d = nbpac 291 ENDIF 289 IF( ji == jiindx .AND. jj == jjindx ) jiindex_1d = nbpac 292 290 ENDIF 293 291 END DO 294 292 END DO 295 293 296 IF( ln_nicep ) THEN 297 WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 298 ENDIF 294 IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac 299 295 300 296 !------------------------------ … … 306 302 IF ( nbpac > 0 ) THEN 307 303 308 CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , & 309 jpi, jpj, npac(1:nbpac) ) 304 CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 310 305 DO jl = 1, jpl 311 CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl) , a_i(:,:,jl) , & 312 jpi, jpj, npac(1:nbpac) ) 313 CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl) , v_i(:,:,jl) , & 314 jpi, jpj, npac(1:nbpac) ) 315 CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , & 316 jpi, jpj, npac(1:nbpac) ) 317 CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), & 318 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) ) 309 CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 319 310 DO jk = 1, nlay_i 320 CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , & 321 jpi, jpj, npac(1:nbpac) ) 311 CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 322 312 END DO ! jk 323 313 END DO ! jl 324 314 325 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , & 326 jpi, jpj, npac(1:nbpac) ) 327 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , & 328 jpi, jpj, npac(1:nbpac) ) 329 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , & 330 jpi, jpj, npac(1:nbpac) ) 331 CALL tab_2d_1d( nbpac, fseqv_1d (1:nbpac) , fseqv , & 332 jpi, jpj, npac(1:nbpac) ) 333 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , & 334 jpi, jpj, npac(1:nbpac) ) 335 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , & 336 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) ) 319 CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac) , rdm_ice, 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) ) 337 322 338 323 !------------------------------------------------------------------------------! … … 344 329 !---------------------- 345 330 DO ji = 1, nbpac 346 zh_newice(ji) 347 END DO 348 IF ( fraz_swi .EQ. 1.0 )zh_newice(:) = hicol_b(:)331 zh_newice(ji) = hiccrit(1) 332 END DO 333 IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:) 349 334 350 335 !---------------------- … … 352 337 !---------------------- 353 338 354 IF ( num_sal .EQ. 1 ) THEN 355 zs_newice(:) = bulk_sal 356 ENDIF ! num_sal 357 358 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN 359 360 DO ji = 1, nbpac 361 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max ) 362 zji = MOD( npac(ji) - 1, jpi ) + 1 363 zjj = ( npac(ji) - 1 ) / jpi + 1 364 zs_newice(ji) = MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) ) 365 END DO ! jl 366 367 ENDIF ! num_sal 368 369 IF ( num_sal .EQ. 3 ) THEN 370 zs_newice(:) = 2.3 371 ENDIF ! num_sal 339 SELECT CASE ( num_sal ) 340 CASE ( 1 ) ! Sice = constant 341 zs_newice(:) = bulk_sal 342 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 343 DO ji = 1, nbpac 344 zji = MOD( npac(ji) - 1 , jpi ) + 1 345 zjj = ( npac(ji) - 1 ) / jpi + 1 346 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj) ) 347 END DO 348 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 349 zs_newice(:) = 2.3 350 END SELECT 351 372 352 373 353 !------------------------- … … 376 356 ! We assume that new ice is formed at the seawater freezing point 377 357 DO ji = 1, nbpac 378 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 379 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 380 + lfus * ( 1.0 - ( ztmelts - rtt ) & 381 / ( t_bo_b(ji) - rtt ) ) & 382 - rcp * ( ztmelts-rtt ) ) 383 ze_newice(ji) = MAX( ze_newice(ji) , 0.0 ) + & 384 MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) & 385 * rhoic * lfus 358 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 359 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 360 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 361 & - rcp * ( ztmelts - rtt ) ) 362 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) & 363 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus 386 364 END DO ! ji 387 365 !---------------- … … 389 367 !---------------- 390 368 DO ji = 1, nbpac 391 zo_newice(ji) = 0.0369 zo_newice(ji) = 0._wp 392 370 END DO ! ji 393 371 … … 396 374 !-------------------------- 397 375 DO ji = 1, nbpac 398 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)!<0376 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0 399 377 END DO ! ji 400 378 … … 403 381 !------------------- 404 382 DO ji = 1, nbpac 405 zv_newice(ji) 383 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 406 384 407 385 ! A fraction zfrazb of frazil ice is accreted at the ice bottom 408 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) & 409 + 1.0 ) / 2.0 * maxfrazb 410 zdh_frazb(ji) = zfrazb*zv_newice(ji) 386 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 387 zdh_frazb(ji) = zfrazb * zv_newice(ji) 411 388 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 412 389 END DO … … 415 392 ! Salt flux due to new ice growth 416 393 !--------------------------------- 417 IF ( ( num_sal .EQ. 4 ) ) THEN 418 DO ji = 1, nbpac 419 zji = MOD( npac(ji) - 1, jpi ) + 1 420 zjj = ( npac(ji) - 1 ) / jpi + 1 421 fseqv_1d(ji) = fseqv_1d(ji) + & 422 ( sss_m(zji,zjj) - bulk_sal ) * rhoic * & 423 zv_newice(ji) / rdt_ice 424 END DO 425 ELSE 426 DO ji = 1, nbpac 427 zji = MOD( npac(ji) - 1, jpi ) + 1 428 zjj = ( npac(ji) - 1 ) / jpi + 1 429 fseqv_1d(ji) = fseqv_1d(ji) + & 430 ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic * & 431 zv_newice(ji) / rdt_ice 432 END DO ! ji 433 ENDIF 394 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 395 DO ji = 1, nbpac 396 sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 397 rdm_ice_1d(ji) = rdm_ice_1d(ji) + rhoic * zv_newice(ji) 398 END DO ! ji 434 399 435 400 !------------------------------------ … … 437 402 !------------------------------------ 438 403 DO ji = 1, nbpac 439 ! Volume440 zj i = MOD( npac(ji) - 1, jpi )+ 1441 zjj = ( npac(ji) - 1 ) / jpi + 1442 vt_i_init(zji,zjj) = vt_i_init(zji,zjj) +zv_newice(ji)443 ! Energy444 zde = ze_newice(ji) / unit_fac445 zde = zde * area(zji,zjj) * zv_newice(ji)446 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 447 412 END DO 448 413 449 414 ! keep new ice volume in memory 450 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , & 451 jpi, jpj ) 415 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj ) 452 416 453 417 !----------------- … … 455 419 !----------------- 456 420 DO ji = 1, nbpac 457 za_newice(ji) = zv_newice(ji) / zh_newice(ji) 458 ! diagnostic 459 zji = MOD( npac(ji) - 1, jpi ) + 1 460 zjj = ( npac(ji) - 1 ) / jpi + 1 461 diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice 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) 424 diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 462 425 END DO !ji 463 426 … … 476 439 !------------------------------------------- 477 440 ! If lateral ice growth gives an ice concentration gt 1, then 478 ! we keep the excessive volume in memory and attribute it later 479 ! to bottom accretion 480 DO ji = 1, nbpac 481 ! vectorize 482 IF ( za_newice(ji) .GT. ( 1.0 - zat_i_ac(ji) ) ) THEN 483 zda_res(ji) = za_newice(ji) - (1.0 - zat_i_ac(ji) ) 484 zdv_res(ji) = zda_res(ji) * zh_newice(ji) 485 za_newice(ji) = za_newice(ji) - zda_res(ji) 486 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) 487 448 ELSE 488 zda_res(ji) = 0. 0489 zdv_res(ji) = 0. 0449 zda_res(ji) = 0._wp 450 zdv_res(ji) = 0._wp 490 451 ENDIF 491 452 END DO ! ji … … 497 458 DO jl = 1, jpl 498 459 DO ji = 1, nbpac 499 IF( hi_max (jl-1) < zh_newice(ji) .AND. &500 & zh_newice(ji) <= hi_max (jl) ) THEN460 IF( hi_max (jl-1) < zh_newice(ji) .AND. & 461 & zh_newice(ji) <= hi_max (jl) ) THEN 501 462 za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 502 463 zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) … … 504 465 zcatac (ji) = jl 505 466 ENDIF 506 END DO ! ji507 END DO ! jl467 END DO 468 END DO 508 469 509 470 !---------------------------------- … … 521 482 DO ji = 1, nbpac 522 483 jl = zcatac(ji) 523 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]484 zqold = ze_i_ac(ji,jk,jl) ! [ J.m-3 ] 524 485 zalphai = MIN( zhice_old(ji,jl) * jk / nlay_i , zh_newice(ji) ) & 525 486 & - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) ) … … 527 488 + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl) * zqold * zhice_old(ji,jl) / nlay_i & 528 489 + za_newice(ji) * ze_newice(ji) * zalphai & 529 + 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 ) 530 491 END DO 531 492 END DO … … 567 528 zdhicbot (ji,jl) = zdv_res(ji) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb & 568 529 & + zindb * zdh_frazb(ji) ! frazil ice may coalesce 569 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 570 531 END DO 571 532 END DO … … 628 589 ! Update salinity 629 590 !----------------- 630 IF( num_sal == 2 .OR. num_sal == 4 ) THEN591 IF( num_sal == 2 ) THEN ! Sice = F(z,t) 631 592 DO jl = 1, jpl 632 593 DO ji = 1, nbpac 633 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes594 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes 634 595 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 635 596 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb … … 645 606 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 646 607 CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 647 IF ( num_sal == 2 .OR. num_sal == 4) &608 IF ( num_sal == 2 ) & 648 609 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 649 610 DO jk = 1, nlay_i … … 651 612 END DO 652 613 END DO 653 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 ) 615 CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj ) 654 616 ! 655 617 ENDIF ! nbpac > 0 … … 660 622 DO jl = 1, jpl 661 623 DO jk = 1, nlay_i ! heat content in 10^9 Joules 662 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 663 625 END DO 664 626 END DO
Note: See TracChangeset
for help on using the changeset viewer.
