- Timestamp:
- 2019-06-06T16:29:54+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0_GO6_mixing/src/ICE/icedyn_rdgrft.F90
r10888 r11083 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 … … 156 152 ENDIF 157 153 158 CALL ice_var_zapsmall ! Zero out categories with very small areas159 160 154 !-------------------------------- 161 155 ! 0) Identify grid cells with ice 162 156 !-------------------------------- 157 at_i(:,:) = SUM( a_i, dim=3 ) 158 ! 163 159 npti = 0 ; nptidx(:) = 0 164 160 ipti = 0 ; iptidx(:) = 0 165 161 DO jj = 1, jpj 166 162 DO ji = 1, jpi 167 IF ( at_i(ji,jj) > 0._wp) THEN163 IF ( at_i(ji,jj) > epsi10 ) THEN 168 164 npti = npti + 1 169 165 nptidx( npti ) = (jj - 1) * jpi + ji … … 178 174 179 175 ! 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(:,:))176 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i ) 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 182 178 ! 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 (:,:))179 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 180 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) 181 CALL tab_2d_1d( npti, nptidx(1:npti), ato_i_1d(1:npti) , ato_i ) 186 182 187 183 DO ji = 1, npti … … 310 306 311 307 ! ! 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._wp308 WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 309 ELSEWHERE ; zhi(1:npti,:) = 0._wp 314 310 END WHERE 315 311 … … 329 325 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 330 326 ! 331 WHERE( zasum(1:npti) > 0._wp) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)332 ELSEWHERE ; z1_asum(1:npti) = 0._wp327 WHERE( zasum(1:npti) > epsi20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 328 ELSEWHERE ; z1_asum(1:npti) = 0._wp 333 329 END WHERE 334 330 ! … … 455 451 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 456 452 ! 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._wp453 WHERE( zaksum(1:npti) > epsi20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 454 ELSEWHERE ; closing_gross(1:npti) = 0._wp 459 455 END WHERE 460 456 … … 466 462 DO ji = 1, npti 467 463 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 468 IF( zfac > pa_i(ji,jl) ) THEN464 IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 469 465 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 470 466 ENDIF … … 510 506 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 511 507 REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a 508 REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i) 512 509 ! 513 510 REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice … … 518 515 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 519 516 !!------------------------------------------------------------------- 520 517 ! 518 zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume 519 ! 521 520 ! 1) Change in open water area due to closing and opening 522 521 !-------------------------------------------------------- … … 535 534 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 536 535 537 z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 538 536 IF( a_i_2d(ji,jl1) > epsi20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 537 ELSE ; z1_ai(ji) = 0._wp 538 ENDIF 539 539 540 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 540 541 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice … … 549 550 550 551 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 551 vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 552 IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg 553 ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0 554 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 555 ENDIF 552 556 ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 553 557 554 558 ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 555 559 virdg1 = v_i_2d (ji,jl1) * afrdg 556 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg )560 virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw 557 561 vsrdg(ji) = v_s_2d (ji,jl1) * afrdg 558 562 sirdg1 = sv_i_2d(ji,jl1) * afrdg … … 726 730 END DO ! jl1 727 731 ! 732 ! roundoff errors 733 !---------------- 728 734 ! 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 735 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 736 ! 732 737 END SUBROUTINE rdgrft_shift … … 854 859 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 855 860 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 856 861 ! 857 862 ! !---------------------! 858 863 CASE( 2 ) !== from 1D to 2D ==! … … 945 950 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 946 951 ENDIF 947 ! ! allocate tke arrays 952 ! 953 IF( .NOT. ln_icethd ) THEN 954 rn_porordg = 0._wp 955 rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 956 rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 957 IF( lwp ) THEN 958 WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed' 959 WRITE(numout,*) ' rn_porordg = ', rn_porordg 960 WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg 961 WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg 962 WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft 963 WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft 964 ENDIF 965 ENDIF 966 ! ! allocate arrays 948 967 IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 949 968 !
Note: See TracChangeset
for help on using the changeset viewer.