Changeset 11564 for NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/ICE/icedyn_rdgrft.F90
- Timestamp:
- 2019-09-18T16:11:52+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10973_AGRIF-01_jchanut_small_jpi_jpj/src/ICE/icedyn_rdgrft.F90
r10945 r11564 142 142 INTEGER, PARAMETER :: jp_itermax = 20 143 143 !!------------------------------------------------------------------- 144 ! clem: The redistribution of ice between categories can lead to small negative values (as for the remapping in ice_itd_rem)145 ! likely due to truncation error ( i.e. 1. - 1. /= 0 )146 ! I do not think it should be a concern since small areas and volumes are erased (in ice_var_zapsmall.F90)147 148 144 ! controls 149 145 IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing 150 146 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 147 IF( ln_icediachk ) CALL ice_cons2D (0, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 151 148 152 149 IF( kt == nit000 ) THEN … … 156 153 ENDIF 157 154 158 CALL ice_var_zapsmall ! Zero out categories with very small areas159 160 155 !-------------------------------- 161 156 ! 0) Identify grid cells with ice 162 157 !-------------------------------- 158 at_i(:,:) = SUM( a_i, dim=3 ) 159 ! 163 160 npti = 0 ; nptidx(:) = 0 164 161 ipti = 0 ; iptidx(:) = 0 165 162 DO jj = 1, jpj 166 163 DO ji = 1, jpi 167 IF ( at_i(ji,jj) > 0._wp) THEN164 IF ( at_i(ji,jj) > epsi10 ) THEN 168 165 npti = npti + 1 169 166 nptidx( npti ) = (jj - 1) * jpi + ji … … 178 175 179 176 ! just needed here 180 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti), divu_i(:,:))181 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti), delta_i(:,:))177 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 182 179 ! needed here and in the iteration loop 183 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i (:,:,:))184 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i (:,:,:))185 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i (:,:))180 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 181 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 182 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 186 183 187 184 DO ji = 1, npti … … 280 277 281 278 ! controls 279 IF( ln_ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints 280 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints 282 281 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 283 IF( ln_ ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints282 IF( ln_icediachk ) CALL ice_cons2D (1, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 284 283 IF( ln_timing ) CALL timing_stop ('icedyn_rdgrft') ! timing 285 284 ! … … 310 309 311 310 ! ! Ice thickness needed for rafting 312 WHERE( pa_i(1:npti,:) > 0._wp) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)313 ELSEWHERE ; zhi(1:npti,:) = 0._wp311 WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 312 ELSEWHERE ; zhi(1:npti,:) = 0._wp 314 313 END WHERE 315 314 … … 329 328 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 330 329 ! 331 WHERE( zasum(1:npti) > 0._wp) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)332 ELSEWHERE ; z1_asum(1:npti) = 0._wp330 WHERE( zasum(1:npti) > epsi20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 331 ELSEWHERE ; z1_asum(1:npti) = 0._wp 333 332 END WHERE 334 333 ! … … 455 454 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 456 455 ! NOTE: 0 < aksum <= 1 457 WHERE( zaksum(1:npti) > 0._wp) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)458 ELSEWHERE ; closing_gross(1:npti) = 0._wp456 WHERE( zaksum(1:npti) > epsi20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 457 ELSEWHERE ; closing_gross(1:npti) = 0._wp 459 458 END WHERE 460 459 … … 466 465 DO ji = 1, npti 467 466 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 468 IF( zfac > pa_i(ji,jl) ) THEN467 IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 469 468 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 470 469 ENDIF … … 510 509 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 511 510 REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a 511 REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i) 512 512 ! 513 513 REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice … … 518 518 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 519 519 !!------------------------------------------------------------------- 520 520 ! 521 zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume 522 ! 521 523 ! 1) Change in open water area due to closing and opening 522 524 !-------------------------------------------------------- … … 535 537 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 536 538 537 z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 538 539 IF( a_i_2d(ji,jl1) > epsi20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 540 ELSE ; z1_ai(ji) = 0._wp 541 ENDIF 542 539 543 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 540 544 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice … … 549 553 550 554 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 551 vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 555 IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg 556 ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0 557 ELSE ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg * MAX( 0._wp, 2._wp - 0.1_wp * zvti(ji) ) ! v > 10m and v < 20m then porosity = linear transition to 0 558 ENDIF 552 559 ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 553 560 554 561 ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 555 562 virdg1 = v_i_2d (ji,jl1) * afrdg 556 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg )563 virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw 557 564 vsrdg(ji) = v_s_2d (ji,jl1) * afrdg 558 565 sirdg1 = sv_i_2d(ji,jl1) * afrdg … … 726 733 END DO ! jl1 727 734 ! 735 ! roundoff errors 736 !---------------- 728 737 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 729 WHERE( a_i_2d(1:npti,:) < 0._wp ) a_i_2d(1:npti,:) = 0._wp 730 WHERE( v_i_2d(1:npti,:) < 0._wp ) v_i_2d(1:npti,:) = 0._wp 738 CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 731 739 ! 732 740 END SUBROUTINE rdgrft_shift … … 854 862 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 855 863 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 856 864 ! 857 865 ! !---------------------! 858 866 CASE( 2 ) !== from 1D to 2D ==! … … 911 919 REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution 912 920 READ ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901) 913 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' , lwp)921 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' ) 914 922 REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution 915 923 READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 ) 916 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' , lwp)924 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' ) 917 925 IF(lwm) WRITE ( numoni, namdyn_rdgrft ) 918 926 !
Note: See TracChangeset
for help on using the changeset viewer.
