- Timestamp:
- 2013-11-07T11:01:27+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r4147 r4161 5 5 !!====================================================================== 6 6 !! History : LIM ! 2006-02 (M. Vancoppenolle) Original code 7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo7 !! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 8 8 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 9 9 !!---------------------------------------------------------------------- … … 12 12 !! 'key_lim3' LIM-3 sea-ice model 13 13 !!---------------------------------------------------------------------- 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 31 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 14 USE par_oce ! ocean parameters 15 USE dom_oce ! ocean domain 16 USE phycst ! physical constants (ocean directory) 17 USE sbc_oce ! surface boundary condition: ocean fields 18 USE thd_ice ! LIM thermodynamics 19 USE ice ! LIM variables 20 USE par_ice ! LIM parameters 21 USE dom_ice ! LIM domain 22 USE limthd_lac ! LIM 23 USE limvar ! LIM 24 USE limcons ! LIM 25 USE in_out_manager ! I/O manager 26 USE lbclnk ! lateral boundary condition - MPP exchanges 27 USE lib_mpp ! MPP library 28 USE wrk_nemo ! work arrays 29 USE prtctl ! Print control 30 ! Check budget (Rousset) 31 USE iom ! I/O manager 32 USE lib_fortran ! glob_sum 33 USE limdiahsb 34 USE timing ! Timing 32 35 33 36 IMPLICIT NONE … … 62 65 REAL(wp), PARAMETER :: krdgmin = 1.1_wp ! min ridge thickness multiplier 63 66 REAL(wp), PARAMETER :: kraft = 2.0_wp ! rafting multipliyer 67 REAL(wp), PARAMETER :: kamax = 1.0 64 68 65 69 REAL(wp) :: Cp ! … … 141 145 REAL(wp), POINTER, DIMENSION(:,:) :: esnow_mlt ! energy needed to melt snow in ocean (J m-2) 142 146 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories 147 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 148 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 149 ! mass and salt flux (clem) 150 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold, zvsold, zsmvold ! old ice volume... 143 151 !!----------------------------------------------------------------------------- 152 IF( nn_timing == 1 ) CALL timing_start('limitd_me') 144 153 145 154 CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 155 156 CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 146 157 147 158 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) … … 151 162 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 152 163 ENDIF 164 165 IF( ln_limdyn ) THEN ! Start ridging and rafting ! 166 ! ------------------------------- 167 !- check conservation (C Rousset) 168 IF (ln_limdiahsb) THEN 169 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 170 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 171 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 172 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 173 ENDIF 174 !- check conservation (C Rousset) 175 ! ------------------------------- 176 177 ! mass and salt flux init (clem) 178 zviold(:,:,:) = v_i(:,:,:) 179 zvsold(:,:,:) = v_s(:,:,:) 180 zsmvold(:,:,:) = smv_i(:,:,:) 153 181 154 182 !-----------------------------------------------------------------------------! … … 204 232 ! to give asum = 1.0 after ridging. 205 233 206 divu_adv(ji,jj) = ( 1._wp- asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep234 divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep 207 235 208 236 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) … … 288 316 DO jj = 1, jpj 289 317 DO ji = 1, jpi 290 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN318 IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 291 319 closing_net(ji,jj) = 0._wp 292 320 opning (ji,jj) = 0._wp 293 321 ELSE 294 322 iterate_ridging = 1 295 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj)) * r1_rdtice323 divu_adv (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 296 324 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 297 325 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) … … 330 358 DO ji = 1, jpi 331 359 332 IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )asum_error = .true.360 IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 333 361 334 362 dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice … … 349 377 DO jj = 1, jpj 350 378 DO ji = 1, jpi 351 IF( ABS( asum(ji,jj) - 1._wp) > epsi11 ) THEN ! there is a bug379 IF( ABS( asum(ji,jj) - kamax) > epsi11 ) THEN ! there is a bug 352 380 WRITE(numout,*) ' ' 353 381 WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) … … 377 405 CALL lim_var_glo2eqv 378 406 CALL lim_itd_me_zapsmall 407 408 !-------------------------------- 409 ! Update mass/salt fluxes (clem) 410 !-------------------------------- 411 DO jl = 1, jpl 412 DO jj = 1, jpj 413 DO ji = 1, jpi 414 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 415 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic 416 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn 417 sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice 418 END DO 419 END DO 420 END DO 379 421 380 422 !----------------- … … 425 467 ENDIF 426 468 469 ! ------------------------------- 470 !- check conservation (C Rousset) 471 IF (ln_limdiahsb) THEN 472 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 473 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 474 475 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 476 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 477 478 zchk_vmin = glob_min(v_i) 479 zchk_amax = glob_max(SUM(a_i,dim=3)) 480 zchk_amin = glob_min(a_i) 481 482 IF(lwp) THEN 483 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limitd_me) = ',(zchk_v_i * rday) 484 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 485 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limitd_me) = ',(zchk_vmin * 1.e-3) 486 IF ( zchk_amax > kamax+epsi10 ) WRITE(numout,*) 'violation a_i>amax (limitd_me) = ',zchk_amax 487 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limitd_me) = ',zchk_amin 488 ENDIF 489 ENDIF 490 !- check conservation (C Rousset) 491 ! ------------------------------- 492 427 493 !-------------------------! 428 494 ! Back to initial values … … 448 514 449 515 ! heat content has to be corrected before ice volume 450 DO jl = 1, jpl 451 DO jk = 1, nlay_i 452 DO jj = 1, jpj 453 DO ji = 1, jpi 454 IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 455 ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 456 old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 457 d_e_i_trp(ji,jj,jk,jl) = 0._wp 458 ENDIF 459 END DO 460 END DO 461 END DO 462 END DO 463 464 DO jl = 1, jpl 465 DO jj = 1, jpj 466 DO ji = 1, jpi 467 IF( old_v_i (ji,jj,jl) < epsi06 .AND. & 468 d_v_i_trp(ji,jj,jl) > epsi06 ) THEN 469 old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl) 470 d_v_i_trp (ji,jj,jl) = 0._wp 471 old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl) 472 d_a_i_trp (ji,jj,jl) = 0._wp 473 old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl) 474 d_v_s_trp (ji,jj,jl) = 0._wp 475 old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 476 d_e_s_trp (ji,jj,1,jl) = 0._wp 477 old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 478 d_oa_i_trp(ji,jj,jl) = 0._wp 479 IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 480 d_smv_i_trp(ji,jj,jl) = 0._wp 481 ENDIF 482 END DO 483 END DO 484 END DO 516 !clem@order 517 ! DO jl = 1, jpl 518 ! DO jk = 1, nlay_i 519 ! DO jj = 1, jpj 520 ! DO ji = 1, jpi 521 ! IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 522 ! ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 523 ! old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 524 ! d_e_i_trp(ji,jj,jk,jl) = 0._wp 525 ! ENDIF 526 ! END DO 527 ! END DO 528 ! END DO 529 ! END DO 530 ! 531 ! DO jl = 1, jpl 532 ! DO jj = 1, jpj 533 ! DO ji = 1, jpi 534 ! IF( old_v_i (ji,jj,jl) < epsi06 .AND. & 535 ! d_v_i_trp(ji,jj,jl) > epsi06 ) THEN 536 ! old_v_i (ji,jj,jl) = d_v_i_trp(ji,jj,jl) 537 ! d_v_i_trp (ji,jj,jl) = 0._wp 538 ! old_a_i (ji,jj,jl) = d_a_i_trp(ji,jj,jl) 539 ! d_a_i_trp (ji,jj,jl) = 0._wp 540 ! old_v_s (ji,jj,jl) = d_v_s_trp(ji,jj,jl) 541 ! d_v_s_trp (ji,jj,jl) = 0._wp 542 ! old_e_s (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 543 ! d_e_s_trp (ji,jj,1,jl) = 0._wp 544 ! old_oa_i (ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 545 ! d_oa_i_trp(ji,jj,jl) = 0._wp 546 ! IF( num_sal == 2 ) old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 547 ! d_smv_i_trp(ji,jj,jl) = 0._wp 548 ! ENDIF 549 ! END DO 550 ! END DO 551 ! END DO 552 !clem@order 553 ENDIF ! ln_limdyn=.true. 485 554 ! 486 555 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 487 556 ! 557 CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold ) ! clem 558 ! 559 IF( nn_timing == 1 ) CALL timing_stop('limitd_me') 488 560 END SUBROUTINE lim_itd_me 489 561 … … 1086 1158 afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 1087 1159 1088 IF (afrac(ji,jj) > 1.0+ epsi11) THEN !riging1160 IF (afrac(ji,jj) > kamax + epsi11) THEN !riging 1089 1161 large_afrac = .true. 1090 ELSEIF (afrac(ji,jj) > 1.0) THEN ! roundoff error1091 afrac(ji,jj) = 1.01162 ELSEIF (afrac(ji,jj) > kamax) THEN ! roundoff error 1163 afrac(ji,jj) = kamax 1092 1164 ENDIF 1093 IF (afrft(ji,jj) > 1.0+ epsi11) THEN !rafting1165 IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 1094 1166 large_afrft = .true. 1095 ELSEIF (afrft(ji,jj) > 1.0) THEN ! roundoff error1096 afrft(ji,jj) = 1.01167 ELSEIF (afrft(ji,jj) > kamax) THEN ! roundoff error 1168 afrft(ji,jj) = kamax 1097 1169 ENDIF 1098 1170 … … 1137 1209 1138 1210 ! ! excess of salt is flushed into the ocean 1139 sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice1140 1141 rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0 ! increase in ice volume du to seawater frozen in voids1142 1211 !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 1212 1213 !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic ! gurvan: increase in ice volume du to seawater frozen in voids 1214 1143 1215 !------------------------------------ 1144 1216 ! 3.6 Increment ridging diagnostics … … 1150 1222 dardg1dt (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 1151 1223 dardg2dt (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 1152 diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice1224 !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 1153 1225 opening (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 1154 1226 … … 1217 1289 1218 1290 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1219 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i1291 ersw (ji,jj,jk) = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 1220 1292 1221 1293 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) … … 1240 1312 ji = indxi(ij) 1241 1313 jj = indxj(ij) 1242 IF( afrac(ji,jj) > 1.0+ epsi11 ) THEN1314 IF( afrac(ji,jj) > kamax + epsi11 ) THEN 1243 1315 WRITE(numout,*) '' 1244 1316 WRITE(numout,*) ' ardg > a_i' … … 1252 1324 ji = indxi(ij) 1253 1325 jj = indxj(ij) 1254 IF( afrft(ji,jj) > 1.0+ epsi11 ) THEN1326 IF( afrft(ji,jj) > kamax + epsi11 ) THEN 1255 1327 WRITE(numout,*) '' 1256 1328 WRITE(numout,*) ' arft > a_i'
Note: See TracChangeset
for help on using the changeset viewer.
