- Timestamp:
- 2012-10-26T12:13:21+02:00 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
r3294 r3517 143 143 ! 1) Conservation check and changes in each ice category 144 144 !------------------------------------------------------------------------------! 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)145 IF( con_i ) THEN 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) 150 150 ENDIF 151 151 … … 158 158 DO ji = 1, jpi 159 159 !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 160 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) , epsi10 ) * nlay_i 161 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) ) ) !0 if no ice and 1 if yes 162 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 165 163 END DO 166 164 END DO … … 182 180 ! 183 181 184 zvrel(:,:) = 0. 0182 zvrel(:,:) = 0._wp 185 183 186 184 ! 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 185 hicol(:,:) = hiccrit(1) 186 187 IF( fraz_swi == 1._wp ) THEN 194 188 195 189 !-------------------- 196 190 ! Physical constants 197 191 !-------------------- 198 hicol(:,:) = 0. 0192 hicol(:,:) = 0._wp 199 193 200 194 zhicrit = 0.04 ! frazil ice thickness … … 306 300 IF ( nbpac > 0 ) THEN 307 301 308 CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , & 309 jpi, jpj, npac(1:nbpac) ) 302 CALL tab_2d_1d( nbpac, zat_i_ac (1:nbpac) , at_i , jpi, jpj, npac(1:nbpac) ) 310 303 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) ) 304 CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl) , a_i(:,:,jl) , jpi, jpj, npac(1:nbpac) ) 305 CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl) , v_i(:,:,jl) , jpi, jpj, npac(1:nbpac) ) 306 CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , jpi, jpj, npac(1:nbpac) ) 307 CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 319 308 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) ) 309 CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 322 310 END DO ! jk 323 311 END DO ! jl 324 312 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) ) 313 CALL tab_2d_1d( nbpac, qldif_1d (1:nbpac) , qldif , jpi, jpj, npac(1:nbpac) ) 314 CALL tab_2d_1d( nbpac, qcmif_1d (1:nbpac) , qcmif , jpi, jpj, npac(1:nbpac) ) 315 CALL tab_2d_1d( nbpac, t_bo_b (1:nbpac) , t_bo , jpi, jpj, npac(1:nbpac) ) 316 CALL tab_2d_1d( nbpac, fseqv_1d (1:nbpac) , fseqv , jpi, jpj, npac(1:nbpac) ) 317 CALL tab_2d_1d( nbpac, hicol_b (1:nbpac) , hicol , jpi, jpj, npac(1:nbpac) ) 318 CALL tab_2d_1d( nbpac, zvrel_ac (1:nbpac) , zvrel , jpi, jpj, npac(1:nbpac) ) 337 319 338 320 !------------------------------------------------------------------------------! … … 344 326 !---------------------- 345 327 DO ji = 1, nbpac 346 zh_newice(ji) 347 END DO 348 IF ( fraz_swi .EQ. 1.0 )zh_newice(:) = hicol_b(:)328 zh_newice(ji) = hiccrit(1) 329 END DO 330 IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:) 349 331 350 332 !---------------------- … … 352 334 !---------------------- 353 335 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 336 SELECT CASE ( num_sal ) 337 CASE ( 1 ) ! Sice = constant 338 zs_newice(:) = bulk_sal 339 CASE ( 2 ) ! Sice = F(z,t) [Vancoppenolle et al (2005)] 340 DO ji = 1, nbpac 341 zji = MOD( npac(ji) - 1 , jpi ) + 1 342 zjj = ( npac(ji) - 1 ) / jpi + 1 343 zs_newice(ji) = MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj) ) 344 END DO 345 CASE ( 3 ) ! Sice = F(z) [multiyear ice] 346 zs_newice(:) = 2.3 347 END SELECT 348 372 349 373 350 !------------------------- … … 376 353 ! We assume that new ice is formed at the seawater freezing point 377 354 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 355 ztmelts = - tmut * zs_newice(ji) + rtt ! Melting point (K) 356 ze_newice(ji) = rhoic * ( cpic * ( ztmelts - t_bo_b(ji) ) & 357 & + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) ) & 358 & - rcp * ( ztmelts - rtt ) ) 359 ze_newice(ji) = MAX( ze_newice(ji) , 0._wp ) & 360 & + MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) ) * rhoic * lfus 386 361 END DO ! ji 387 362 !---------------- … … 389 364 !---------------- 390 365 DO ji = 1, nbpac 391 zo_newice(ji) = 0.0366 zo_newice(ji) = 0._wp 392 367 END DO ! ji 393 368 … … 396 371 !-------------------------- 397 372 DO ji = 1, nbpac 398 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)!<0373 zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji) !<0 399 374 END DO ! ji 400 375 … … 403 378 !------------------- 404 379 DO ji = 1, nbpac 405 zv_newice(ji) 380 zv_newice(ji) = - zqbgow(ji) / ze_newice(ji) 406 381 407 382 ! 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) 383 zfrazb = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 384 zdh_frazb(ji) = zfrazb * zv_newice(ji) 411 385 zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 412 386 END DO … … 415 389 ! Salt flux due to new ice growth 416 390 !--------------------------------- 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 391 ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 392 DO ji = 1, nbpac 393 fseqv_1d(ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 394 END DO ! ji 434 395 435 396 !------------------------------------ … … 448 409 449 410 ! keep new ice volume in memory 450 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , & 451 jpi, jpj ) 411 CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj ) 452 412 453 413 !----------------- … … 459 419 zji = MOD( npac(ji) - 1, jpi ) + 1 460 420 zjj = ( npac(ji) - 1 ) / jpi + 1 461 diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice421 diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice 462 422 END DO !ji 463 423 … … 628 588 ! Update salinity 629 589 !----------------- 630 IF( num_sal == 2 .OR. num_sal == 4 ) THEN590 IF( num_sal == 2 ) THEN ! Sice = F(z,t) 631 591 DO jl = 1, jpl 632 592 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 yes593 zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) ! 0 if no ice and 1 if yes 634 594 zdv = zv_i_ac(ji,jl) - zv_old(ji,jl) 635 595 zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb … … 645 605 CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 646 606 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) &607 IF ( num_sal == 2 ) & 648 608 CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 649 609 DO jk = 1, nlay_i
Note: See TracChangeset
for help on using the changeset viewer.
