Changeset 11822 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rdgrft.F90
- Timestamp:
- 2019-10-29T11:41:36+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rdgrft.F90
r10531 r11822 86 86 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 87 !!------------------------------------------------------------------- 88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij),&89 & apartf(jpij,0:jpl) , hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl),&90 & hrmax (jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl),&88 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 89 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 90 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 91 91 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 92 … … 137 137 REAL(wp) :: zfac ! local scalar 138 138 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 139 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s)140 139 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 140 ! 142 141 INTEGER, PARAMETER :: jp_itermax = 20 143 142 !!------------------------------------------------------------------- 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 143 ! controls 149 144 IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing 150 145 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 146 IF( ln_icediachk ) CALL ice_cons2D (0, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 151 147 152 148 IF( kt == nit000 ) THEN … … 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), zdelt (1:npti) , delta_i ) 182 177 ! 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(:,:) ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_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 … … 190 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 191 187 ! 192 ! divergence given by the advection scheme 193 ! (which may not be equal to divu as computed from the velocity field) 194 IF ( ln_adv_Pra ) THEN 195 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 196 ELSEIF( ln_adv_UMx ) THEN 197 zdivu_adv(ji) = zdivu(ji) 198 ENDIF 199 ! 200 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough 201 ! ! to give asum = 1.0 after ridging 188 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough 189 ! ! to give asum = 1.0 after ridging 202 190 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 203 opning(ji) = closing_net(ji) + zdivu _adv(ji)191 opning(ji) = closing_net(ji) + zdivu(ji) 204 192 END DO 205 193 ! … … 218 206 ato_i_1d (ipti) = ato_i_1d (ji) 219 207 closing_net(ipti) = closing_net(ji) 220 zdivu _adv (ipti) = zdivu_adv(ji)208 zdivu (ipti) = zdivu (ji) 221 209 opning (ipti) = opning (ji) 222 210 ENDIF … … 262 250 ELSE 263 251 iterate_ridging = 1 264 zdivu _adv(ji) = zfac * r1_rdtice265 closing_net(ji) = MAX( 0._wp, -zdivu _adv(ji) )266 opning (ji) = MAX( 0._wp, zdivu _adv(ji) )252 zdivu (ji) = zfac * r1_rdtice 253 closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 254 opning (ji) = MAX( 0._wp, zdivu(ji) ) 267 255 ENDIF 268 256 END DO … … 280 268 281 269 ! controls 270 IF( ln_ctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints 271 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints 282 272 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') ! prints273 IF( ln_icediachk ) CALL ice_cons2D (1, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 284 274 IF( ln_timing ) CALL timing_stop ('icedyn_rdgrft') ! timing 285 275 ! … … 310 300 311 301 ! ! 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._wp302 WHERE( pa_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 303 ELSEWHERE ; zhi(1:npti,:) = 0._wp 314 304 END WHERE 315 305 … … 329 319 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 330 320 ! 331 WHERE( zasum(1:npti) > 0._wp) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)332 ELSEWHERE ; z1_asum(1:npti) = 0._wp321 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 322 ELSEWHERE ; z1_asum(1:npti) = 0._wp 333 323 END WHERE 334 324 ! … … 455 445 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 456 446 ! 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._wp447 WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 448 ELSEWHERE ; closing_gross(1:npti) = 0._wp 459 449 END WHERE 460 450 … … 466 456 DO ji = 1, npti 467 457 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 468 IF( zfac > pa_i(ji,jl) ) THEN458 IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 469 459 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_rdtice 470 460 ENDIF … … 510 500 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 511 501 REAL(wp), DIMENSION(jpij) :: z1_ai ! 1 / a 502 REAL(wp), DIMENSION(jpij) :: zvti ! sum(v_i) 512 503 ! 513 504 REAL(wp), DIMENSION(jpij,nlay_s) :: esrft ! snow energy of rafting ice … … 518 509 INTEGER , DIMENSION(jpij) :: itest_rdg, itest_rft ! test for conservation 519 510 !!------------------------------------------------------------------- 520 511 ! 512 zvti(1:npti) = SUM( v_i_2d(1:npti,:), dim=2 ) ! total ice volume 513 ! 521 514 ! 1) Change in open water area due to closing and opening 522 515 !-------------------------------------------------------- … … 535 528 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 536 529 537 z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 538 530 IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 531 ELSE ; z1_ai(ji) = 0._wp 532 ENDIF 533 539 534 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 540 535 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice … … 549 544 550 545 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 551 vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 546 IF ( zvti(ji) <= 10. ) THEN ; vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg ! v <= 10m then porosity = rn_porordg 547 ELSEIF( zvti(ji) >= 20. ) THEN ; vsw = 0._wp ! v >= 20m then porosity = 0 548 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 549 ENDIF 552 550 ersw(ji) = -rhoi * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 553 551 554 552 ! volume etc of ridging / rafting ice and new ridges (vi, vs, sm, oi, es, ei) 555 553 virdg1 = v_i_2d (ji,jl1) * afrdg 556 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg )554 virdg2(ji) = v_i_2d (ji,jl1) * afrdg + vsw 557 555 vsrdg(ji) = v_s_2d (ji,jl1) * afrdg 558 556 sirdg1 = sv_i_2d(ji,jl1) * afrdg … … 588 586 ! virtual salt flux to keep salinity constant 589 587 IF( nn_icesal /= 2 ) THEN 590 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) 588 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 591 589 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice & ! put back sss_m into the ocean 592 590 & - s_i_1d(ji) * vsw * rhoi * r1_rdtice ! and get s_i from the ocean … … 726 724 END DO ! jl1 727 725 ! 726 ! roundoff errors 727 !---------------- 728 728 ! 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 729 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 730 ! 732 731 END SUBROUTINE rdgrft_shift … … 854 853 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 855 854 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd (:,:) ) 856 855 ! 857 856 ! !---------------------! 858 857 CASE( 2 ) !== from 1D to 2D ==! … … 911 910 REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution 912 911 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)912 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' ) 914 913 REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution 915 914 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)915 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in configuration namelist' ) 917 916 IF(lwm) WRITE ( numoni, namdyn_rdgrft ) 918 917 ! … … 945 944 CALL ctl_stop( 'ice_dyn_rdgrft_init: choose one and only one participation function (ln_partf_lin or ln_partf_exp)' ) 946 945 ENDIF 947 ! ! allocate tke arrays 946 ! 947 IF( .NOT. ln_icethd ) THEN 948 rn_porordg = 0._wp 949 rn_fsnwrdg = 1._wp ; rn_fsnwrft = 1._wp 950 rn_fpndrdg = 1._wp ; rn_fpndrft = 1._wp 951 IF( lwp ) THEN 952 WRITE(numout,*) ' ==> only ice dynamics is activated, thus some parameters must be changed' 953 WRITE(numout,*) ' rn_porordg = ', rn_porordg 954 WRITE(numout,*) ' rn_fsnwrdg = ', rn_fsnwrdg 955 WRITE(numout,*) ' rn_fpndrdg = ', rn_fpndrdg 956 WRITE(numout,*) ' rn_fsnwrft = ', rn_fsnwrft 957 WRITE(numout,*) ' rn_fpndrft = ', rn_fpndrft 958 ENDIF 959 ENDIF 960 ! ! allocate arrays 948 961 IF( ice_dyn_rdgrft_alloc() /= 0 ) CALL ctl_stop( 'STOP', 'ice_dyn_rdgrft_init: unable to allocate arrays' ) 949 962 !
Note: See TracChangeset
for help on using the changeset viewer.