- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icedyn_rdgrft.F90 (modified) (24 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/ICE/icedyn_rdgrft.F90
r10994 r13463 75 75 REAL(wp) :: rn_fpndrft ! fractional pond loss to the ocean during rafting 76 76 ! 77 !! * Substitutions 78 # include "do_loop_substitute.h90" 77 79 !!---------------------------------------------------------------------- 78 80 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 86 88 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 89 !!------------------------------------------------------------------- 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),&90 ALLOCATE( closing_net(jpij) , opning(jpij) , closing_gross(jpij) , & 91 & apartf(jpij,0:jpl) , hrmin (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), & 92 & hrmax (jpij,jpl) , hi_hrdg(jpij,jpl) , araft(jpij,jpl) , & 91 93 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 92 94 … … 137 139 REAL(wp) :: zfac ! local scalar 138 140 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 141 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 142 ! … … 145 146 IF( ln_timing ) CALL timing_start('icedyn_rdgrft') ! timing 146 147 IF( ln_icediachk ) CALL ice_cons_hsm(0, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 148 IF( ln_icediachk ) CALL ice_cons2D (0, 'icedyn_rdgrft', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) ! conservation 147 149 148 150 IF( kt == nit000 ) THEN … … 159 161 npti = 0 ; nptidx(:) = 0 160 162 ipti = 0 ; iptidx(:) = 0 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 IF ( at_i(ji,jj) > epsi10 ) THEN 164 npti = npti + 1 165 nptidx( npti ) = (jj - 1) * jpi + ji 166 ENDIF 167 END DO 168 END DO 163 DO_2D( 1, 1, 1, 1 ) 164 IF ( at_i(ji,jj) > epsi10 ) THEN 165 npti = npti + 1 166 nptidx( npti ) = (jj - 1) * jpi + ji 167 ENDIF 168 END_2D 169 169 170 170 !-------------------------------------------------------- … … 174 174 175 175 ! just needed here 176 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i )177 176 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 178 177 ! needed here and in the iteration loop 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 179 CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d (1:npti,1:jpl), a_i ) 180 180 CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d (1:npti,1:jpl), v_i ) … … 186 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 187 187 ! 188 ! divergence given by the advection scheme 189 ! (which may not be equal to divu as computed from the velocity field) 190 IF ( ln_adv_Pra ) THEN 191 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 192 ELSEIF( ln_adv_UMx ) THEN 193 zdivu_adv(ji) = zdivu(ji) 194 ENDIF 195 ! 196 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) ! make sure the closing rate is large enough 197 ! ! 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 198 190 ! Opening rate (non-negative) that will give asum = 1.0 after ridging. 199 opning(ji) = closing_net(ji) + zdivu _adv(ji)191 opning(ji) = closing_net(ji) + zdivu(ji) 200 192 END DO 201 193 ! … … 214 206 ato_i_1d (ipti) = ato_i_1d (ji) 215 207 closing_net(ipti) = closing_net(ji) 216 zdivu _adv (ipti) = zdivu_adv(ji)208 zdivu (ipti) = zdivu (ji) 217 209 opning (ipti) = opning (ji) 218 210 ENDIF … … 258 250 ELSE 259 251 iterate_ridging = 1 260 zdivu _adv (ji) = zfac * r1_rdtice261 closing_net(ji) = MAX( 0._wp, -zdivu _adv(ji) )262 opning (ji) = MAX( 0._wp, zdivu _adv(ji) )252 zdivu (ji) = zfac * r1_Dt_ice 253 closing_net(ji) = MAX( 0._wp, -zdivu(ji) ) 254 opning (ji) = MAX( 0._wp, zdivu(ji) ) 263 255 ENDIF 264 256 END DO … … 276 268 277 269 ! controls 270 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('icedyn_rdgrft') ! prints 271 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt,-1, ' - ice dyn rdgrft - ') ! prints 278 272 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'icedyn_rdgrft', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) ! conservation 279 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 280 274 IF( ln_timing ) CALL timing_stop ('icedyn_rdgrft') ! timing 281 275 ! … … 306 300 307 301 ! ! Ice thickness needed for rafting 308 WHERE( pa_i(1:npti,:) > epsi20 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 302 ! In single precision there were floating point invalids due a sqrt of zhi which happens to have negative values 303 ! To solve that an extra check about the value of pv_i was added. 304 ! Although adding this condition is safe, the double definition (one for single other for double) has been kept to preserve the results of the sette test. 305 #if defined key_single 306 307 WHERE( pa_i(1:npti,:) > epsi10 .and. pv_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 308 #else 309 WHERE( pa_i(1:npti,:) > epsi10 ) ; zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:) 310 #endif 309 311 ELSEWHERE ; zhi(1:npti,:) = 0._wp 310 312 END WHERE … … 325 327 zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 ) 326 328 ! 327 WHERE( zasum(1:npti) > epsi 20 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti)329 WHERE( zasum(1:npti) > epsi10 ) ; z1_asum(1:npti) = 1._wp / zasum(1:npti) 328 330 ELSEWHERE ; z1_asum(1:npti) = 0._wp 329 331 END WHERE … … 451 453 ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate. 452 454 ! NOTE: 0 < aksum <= 1 453 WHERE( zaksum(1:npti) > epsi 20 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)455 WHERE( zaksum(1:npti) > epsi10 ) ; closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti) 454 456 ELSEWHERE ; closing_gross(1:npti) = 0._wp 455 457 END WHERE … … 461 463 DO jl = 1, jpl 462 464 DO ji = 1, npti 463 zfac = apartf(ji,jl) * closing_gross(ji) * r dt_ice465 zfac = apartf(ji,jl) * closing_gross(ji) * rDt_ice 464 466 IF( zfac > pa_i(ji,jl) .AND. apartf(ji,jl) /= 0._wp ) THEN 465 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_ rdtice467 closing_gross(ji) = pa_i(ji,jl) / apartf(ji,jl) * r1_Dt_ice 466 468 ENDIF 467 469 END DO … … 473 475 ! Reduce the opening rate in proportion 474 476 DO ji = 1, npti 475 zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * r dt_ice477 zfac = pato_i(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice 476 478 IF( zfac < 0._wp ) THEN ! would lead to negative ato_i 477 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_ rdtice479 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_Dt_ice 478 480 ELSEIF( zfac > zasum(ji) ) THEN ! would lead to ato_i > asum 479 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_ rdtice481 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_Dt_ice 480 482 ENDIF 481 483 END DO … … 521 523 !-------------------------------------------------------- 522 524 DO ji = 1, npti 523 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * r dt_ice )525 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rDt_ice ) 524 526 END DO 525 527 … … 534 536 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN ! only if ice is ridging 535 537 536 IF( a_i_2d(ji,jl1) > epsi 20 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)538 IF( a_i_2d(ji,jl1) > epsi10 ) THEN ; z1_ai(ji) = 1._wp / a_i_2d(ji,jl1) 537 539 ELSE ; z1_ai(ji) = 0._wp 538 540 ENDIF 539 541 540 542 ! area of ridging / rafting ice (airdg1) and of new ridge (airdg2) 541 airdg1 = aridge(ji,jl1) * closing_gross(ji) * r dt_ice542 airft1 = araft (ji,jl1) * closing_gross(ji) * r dt_ice543 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rDt_ice 544 airft1 = araft (ji,jl1) * closing_gross(ji) * rDt_ice 543 545 544 546 airdg2(ji) = airdg1 * hi_hrdg(ji,jl1) … … 581 583 582 584 ! Ice-ocean exchanges associated with ice porosity 583 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_ rdtice ! increase in ice volume due to seawater frozen in voids584 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_ rdtice585 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_ rdtice ! > 0 [W.m-2]585 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoi * r1_Dt_ice ! increase in ice volume due to seawater frozen in voids 586 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoi * r1_Dt_ice 587 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_Dt_ice ! > 0 [W.m-2] 586 588 587 589 ! Put the snow lost by ridging into the ocean 588 590 ! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 589 591 wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhos * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean 590 & + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_ rdtice592 & + rhos * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 591 593 592 594 ! virtual salt flux to keep salinity constant 593 595 IF( nn_icesal /= 2 ) THEN 594 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) 595 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_ rdtice & ! put back sss_m into the ocean596 & - s_i_1d(ji) * vsw * rhoi * r1_ rdtice ! and get s_i from the ocean596 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - s_i_1d(ji) ) ! ridge salinity = s_i 597 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_Dt_ice & ! put back sss_m into the ocean 598 & - s_i_1d(ji) * vsw * rhoi * r1_Dt_ice ! and get s_i from the ocean 597 599 ENDIF 598 600 … … 617 619 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 618 620 ! Compute ridging /rafting fractions 619 afrdg = aridge(ji,jl1) * closing_gross(ji) * r dt_ice * z1_ai(ji)620 afrft = araft (ji,jl1) * closing_gross(ji) * r dt_ice * z1_ai(ji)621 afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) 622 afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) 621 623 ! Compute ridging /rafting ice and new ridges for es 622 624 esrdg(ji,jk) = ze_s_2d (ji,jk,jl1) * afrdg … … 624 626 ! Put the snow lost by ridging into the ocean 625 627 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji,jk) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2) 626 & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_ rdtice628 & - esrft(ji,jk) * ( 1._wp - rn_fsnwrft ) ) * r1_Dt_ice 627 629 ! 628 630 ! Remove energy of new ridge to each category jl1 … … 638 640 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 639 641 ! Compute ridging /rafting fractions 640 afrdg = aridge(ji,jl1) * closing_gross(ji) * r dt_ice * z1_ai(ji)641 afrft = araft (ji,jl1) * closing_gross(ji) * r dt_ice * z1_ai(ji)642 afrdg = aridge(ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) 643 afrft = araft (ji,jl1) * closing_gross(ji) * rDt_ice * z1_ai(ji) 642 644 ! Compute ridging ice and new ridges for ei 643 645 eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i … … 772 774 ! !--------------------------------------------------! 773 775 CASE( 1 ) !--- Spatial smoothing 774 DO jj = 2, jpjm1 775 DO ji = 2, jpim1 776 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 777 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 778 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 779 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 780 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 781 ELSE 782 zworka(ji,jj) = 0._wp 783 ENDIF 784 END DO 785 END DO 776 DO_2D( 0, 0, 0, 0 ) 777 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 778 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 779 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & 780 & + strength(ji,jj-1) * tmask(ji,jj-1,1) + strength(ji,jj+1) * tmask(ji,jj+1,1) & 781 & ) / ( 4.0 + tmask(ji-1,jj,1) + tmask(ji+1,jj,1) + tmask(ji,jj-1,1) + tmask(ji,jj+1,1) ) 782 ELSE 783 zworka(ji,jj) = 0._wp 784 ENDIF 785 END_2D 786 786 787 DO jj = 2, jpjm1 788 DO ji = 2, jpim1 789 strength(ji,jj) = zworka(ji,jj) 790 END DO 791 END DO 792 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 787 DO_2D( 0, 0, 0, 0 ) 788 strength(ji,jj) = zworka(ji,jj) 789 END_2D 790 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) 793 791 ! 794 792 CASE( 2 ) !--- Temporal smoothing … … 798 796 ENDIF 799 797 ! 800 DO jj = 2, jpjm1 801 DO ji = 2, jpim1 802 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 803 itframe = 1 ! number of time steps for the running mean 804 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 805 IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 806 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe 807 zstrp2 (ji,jj) = zstrp1 (ji,jj) 808 zstrp1 (ji,jj) = strength(ji,jj) 809 strength(ji,jj) = zp 810 ENDIF 811 END DO 812 END DO 813 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1. ) 798 DO_2D( 0, 0, 0, 0 ) 799 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 800 itframe = 1 ! number of time steps for the running mean 801 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 802 IF ( zstrp2(ji,jj) > 0._wp ) itframe = itframe + 1 803 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / itframe 804 zstrp2 (ji,jj) = zstrp1 (ji,jj) 805 zstrp1 (ji,jj) = strength(ji,jj) 806 strength(ji,jj) = zp 807 ENDIF 808 END_2D 809 CALL lbc_lnk( 'icedyn_rdgrft', strength, 'T', 1.0_wp ) 814 810 ! 815 811 END SELECT … … 914 910 !!------------------------------------------------------------------- 915 911 ! 916 REWIND( numnam_ice_ref ) ! Namelist namicetdme in reference namelist : Ice mechanical ice redistribution917 912 READ ( numnam_ice_ref, namdyn_rdgrft, IOSTAT = ios, ERR = 901) 918 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist', lwp ) 919 REWIND( numnam_ice_cfg ) ! Namelist namdyn_rdgrft in configuration namelist : Ice mechanical ice redistribution 913 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namdyn_rdgrft in reference namelist' ) 920 914 READ ( numnam_ice_cfg, namdyn_rdgrft, IOSTAT = ios, ERR = 902 ) 921 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' ) 922 916 IF(lwm) WRITE ( numoni, namdyn_rdgrft ) 923 917 !
Note: See TracChangeset
for help on using the changeset viewer.
