Changeset 8554
- Timestamp:
- 2017-09-21T17:25:29+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 3 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
r8534 r8554 47 47 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_opw_1d 48 48 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_snw_1d 49 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_dyn_1d 49 50 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_rem_1d 50 51 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hfx_err_dif_1d … … 61 62 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_sni_1d 62 63 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_sum_1d 64 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_dyn_1d 63 65 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_sub_1d 64 66 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_snw_sub_1d … … 74 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_res_1d 75 77 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_spr_1d 78 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_dyn_1d 79 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: wfx_pnd_1d 76 80 77 81 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_bri_1d … … 84 88 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_sub_1d 85 89 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_lam_1d 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sfx_dyn_1d 86 91 87 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sprecip_1d !: <==> the 2D sprecip 88 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: at_i_1d !: <==> the 2D at_i 94 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ato_i_1d !: <==> the 2D ato_i 89 95 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhtur_1d !: <==> the 2D fhtur 90 96 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fhld_1d !: <==> the 2D fhld … … 100 106 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_1d !: <==> the 2D a_i 101 107 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ib_1d !: <==> the 2D a_i_b 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_1d !: <==> the 2D ht_i103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_ib_1d !: <==> the 2D ht_i_b104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_1d !: <==> the 2D ht_s108 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_1d !: 109 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_ib_1d !: 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_1d !: 105 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux 106 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux … … 117 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_s_1d !: 118 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: smv_i_1d !: 125 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: oa_i_1d !: 119 126 120 127 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 174 181 & t_bo_1d (jpij) , & 175 182 & hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) , & 176 & hfx_dif_1d(jpij) , hfx_opw_1d(jpij) , 183 & hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,hfx_dyn_1d(jpij) , & 177 184 & rn_amax_1d(jpij) , & 178 185 & hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , & … … 181 188 ! 182 189 ii = ii + 1 183 ALLOCATE( sprecip_1d (jpij) , at_i_1d (jpij) , &190 ALLOCATE( sprecip_1d (jpij) , at_i_1d (jpij) , ato_i_1d(jpij) , & 184 191 & fhtur_1d (jpij) , wfx_snw_sni_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , & 185 192 & fhld_1d (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) , & 186 193 & wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) , & 187 & wfx_snw_sub_1d(jpij), wfx_ ice_sub_1d(jpij), wfx_err_sub_1d(jpij) ,&188 & wfx_lam_1d(jpij) , dqns_ice_1d(jpij) , evap_ice_1d (jpij),&194 & wfx_snw_sub_1d(jpij), wfx_snw_dyn_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) , & 195 & wfx_lam_1d(jpij) , wfx_dyn_1d(jpij), wfx_pnd_1d(jpij),dqns_ice_1d(jpij) , evap_ice_1d (jpij), & 189 196 & qprec_ice_1d(jpij), i0 (jpij) , & 190 197 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 191 198 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij), & 192 & sfx_lam_1d (jpij) , hicol_1d (jpij), STAT=ierr(ii) )199 & sfx_lam_1d (jpij) , sfx_dyn_1d(jpij), hicol_1d (jpij) , STAT=ierr(ii) ) 193 200 ! 194 201 ii = ii + 1 … … 197 204 & dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) , & 198 205 & dh_i_bott(jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new(jpij) , & 199 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , smv_i_1d(jpij) , STAT=ierr(ii) ) 206 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , v_s_1d (jpij) , & 207 & smv_i_1d (jpij) , oa_i_1d (jpij) , STAT=ierr(ii) ) 200 208 ! 201 209 ii = ii + 1 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
r8534 r8554 21 21 USE ice1D ! sea-ice: thermodynamics 22 22 USE ice ! sea-ice: variables 23 USE icetab ! sea-ice: 1D <==> 2D transformation 23 24 USE icevar ! sea-ice: operations 24 25 USE icectl ! sea-ice: control prints … … 39 40 40 41 ! Variables shared among ridging subroutines 41 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: asum ! sum of total ice and open water area 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aksum ! ratio of area removed to area ridged 43 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmin ! minimum ridge thickness 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hrmax ! maximum ridge thickness 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hraft ! thickness of rafted ice 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: aridge ! participating ice ridging 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: araft ! participating ice rafting 42 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_net ! net rate at which area is removed (1/s) 43 ! ! (ridging ice area - area of new ridges) / dt 44 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: opning ! rate of opening due to divergence/shear 45 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: closing_gross ! rate at which area removed, not counting area of new ridges 46 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: apartf ! participation function; fraction of ridging/closing associated w/ category n 47 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmin ! minimum ridge thickness 48 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hrmax ! maximum ridge thickness 49 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hraft ! thickness of rafted ice 50 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: hi_hrdg ! thickness of ridging ice / mean ridge thickness 51 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: aridge ! participating ice ridging 52 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: araft ! participating ice rafting 53 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_i_2d 54 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ze_s_2d 50 55 ! 51 56 REAL(wp), PARAMETER :: hrdg_hi_min = 1.1_wp ! min ridge thickness multiplier: min(hrdg/hi) 52 57 REAL(wp), PARAMETER :: hi_hrft = 0.5_wp ! rafting multipliyer: (hi/hraft) 53 REAL(wp) :: zdrho !54 58 ! 55 59 ! ** namelist (namdyn_rdgrft) ** … … 83 87 84 88 INTEGER FUNCTION ice_dyn_rdgrft_alloc() 85 !!------------------------------------------------------------------- --!89 !!------------------------------------------------------------------- 86 90 !! *** ROUTINE ice_dyn_rdgrft_alloc *** 87 !!---------------------------------------------------------------------! 88 ALLOCATE( asum (jpi,jpj) , apartf (jpi,jpj,0:jpl) , aksum (jpi,jpj) , & 89 & hrmin(jpi,jpj,jpl) , hraft (jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 90 & hrmax(jpi,jpj,jpl) , hi_hrdg(jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=ice_dyn_rdgrft_alloc ) 91 !!------------------------------------------------------------------- 92 ALLOCATE( closing_net(jpij), opning(jpij) , closing_gross(jpij), & 93 & apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl) , aridge(jpij,jpl), & 94 & hrmax(jpij,jpl), hi_hrdg(jpij,jpl) , araft (jpij,jpl), & 95 & ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc ) 91 96 92 97 IF( lk_mpp ) CALL mpp_sum ( ice_dyn_rdgrft_alloc ) … … 97 102 98 103 SUBROUTINE ice_dyn_rdgrft( kt ) 99 !!------------------------------------------------------------------- --!104 !!------------------------------------------------------------------- 100 105 !! *** ROUTINE ice_dyn_rdgrft *** 101 106 !! … … 121 126 !! This routine is based on CICE code and authors William H. Lipscomb, 122 127 !! and Elizabeth C. Hunke, LANL are gratefully acknowledged 123 !!------------------------------------------------------------------- -!128 !!------------------------------------------------------------------- 124 129 INTEGER, INTENT(in) :: kt ! number of iteration 125 130 !! 126 INTEGER :: ji, jj, jk, jl ! dummy loop index 127 INTEGER :: niter ! local integer 128 INTEGER :: iterate_ridging ! if =1, repeat the ridging 129 REAL(wp) :: zfac ! local scalar 130 REAL(wp), DIMENSION(jpi,jpj) :: closing_net ! net rate at which area is removed (1/s) 131 ! ! (ridging ice area - area of new ridges) / dt 132 REAL(wp), DIMENSION(jpi,jpj) :: divu_adv ! divu as implied by transport scheme (1/s) 133 REAL(wp), DIMENSION(jpi,jpj) :: opning ! rate of opening due to divergence/shear 134 REAL(wp), DIMENSION(jpi,jpj) :: closing_gross ! rate at which area removed, not counting area of new ridges 131 INTEGER :: ji, jj, jk, jl ! dummy loop index 132 INTEGER :: niter, iterate_ridging ! local integer 133 INTEGER :: nidx2 ! local integer 134 REAL(wp) :: zfac ! local scalar 135 INTEGER , DIMENSION(jpij) :: idxice2 ! compute ridge/raft or not 136 REAL(wp), DIMENSION(jpij) :: zdivu_adv ! divu as implied by transport scheme (1/s) 137 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 135 138 ! 136 139 INTEGER, PARAMETER :: nitermax = 20 137 !!------------------------------------------------------------------- ----------140 !!------------------------------------------------------------------- 138 141 ! controls 139 142 IF( nn_timing == 1 ) CALL timing_start('icedyn_rdgrft') ! timing … … 143 146 IF(lwp) WRITE(numout,*) 144 147 IF(lwp) WRITE(numout,*)'ice_dyn_rdgrft: ice ridging and rafting' 145 IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~' 146 ! 147 zdrho = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0 ! proport const for PE 148 ! 148 IF(lwp) WRITE(numout,*)'~~~~~~~~~~~~~~' 149 149 ENDIF 150 150 151 !-----------------------------------------------------------------------------!152 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 153 !-------------------------------- ---------------------------------------------!154 ! 155 CALL rdgrft_prep ! prepare ridging156 !157 DO jj = 1, jpj ! Initialize arrays.151 CALL ice_var_zapsmall ! Zero out categories with very small areas 152 153 !-------------------------------- 154 ! 0) Identify grid cells with ice 155 !-------------------------------- 156 nidx = 0 ; idxice(:) = 0 157 DO jj = 1, jpj 158 158 DO ji = 1, jpi 159 IF ( SUM(a_i(ji,jj,:)) > 0._wp ) THEN 160 nidx = nidx + 1 161 idxice( nidx ) = (jj - 1) * jpi + ji 162 ENDIF 163 END DO 164 END DO 165 166 IF( nidx > 0 ) THEN 167 168 ! just needed here 169 CALL tab_2d_1d( nidx, idxice(1:nidx), zdivu(1:nidx), divu_i(:,:) ) 170 CALL tab_2d_1d( nidx, idxice(1:nidx), zdelt(1:nidx), delta_i(:,:) ) 171 ! needed here and in the iteration loop 172 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i(:,:,:) ) 173 CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i(:,:,:) ) 174 CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx) , ato_i(:,:) ) 175 176 DO ji = 1, nidx 159 177 !-----------------------------------------------------------------------------! 160 178 ! 2) Dynamical inputs (closing rate, divu_adv, opning) … … 176 194 ! (thick, newly ridged ice). 177 195 178 closing_net(ji ,jj) = rn_csrdg * 0.5_wp * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )196 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 179 197 180 198 ! 2.2 divu_adv … … 186 204 ! If divu_adv < 0, make sure the closing rate is large enough 187 205 ! to give asum = 1.0 after ridging. 188 189 divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice ! asum found in ridgeprep 190 191 IF( divu_adv(ji,jj) < 0._wp ) closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 206 zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice 207 208 IF( zdivu_adv(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) ) 192 209 193 210 ! 2.3 opning 194 211 !------------ 195 212 ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging. 196 opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 197 END DO 198 END DO 199 200 !-----------------------------------------------------------------------------! 201 ! 3) Ridging iteration 202 !-----------------------------------------------------------------------------! 203 niter = 1 ! iteration counter 204 iterate_ridging = 1 205 206 DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 207 208 ! 3.1 closing_gross 213 opning(ji) = closing_net(ji) + zdivu_adv(ji) 214 END DO 215 216 ! 217 !------------------------------------------------ 218 ! 2.1) Identify grid cells with nonzero ridging 219 !------------------------------------------------ 220 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 221 nidx2 = 0 ; idxice2(:) = 0 222 DO ji = 1, nidx 223 IF( SUM( ABS(apartf(ji,1:jpl)) ) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 224 nidx2 = nidx2 + 1 225 idxice2(nidx2) = idxice(ji) 226 a_i_2d(nidx2,:) = a_i_2d(ji,:) ! adjust to new indices 227 v_i_2d(nidx2,:) = v_i_2d(ji,:) ! adjust to new indices 228 ato_i_1d(nidx2) = ato_i_1d(ji) ! adjust to new indices 229 ENDIF 230 END DO 231 idxice(:) = idxice2(:) 232 nidx = nidx2 233 234 ENDIF 235 236 !-------------- 237 ! Start ridging 238 !-------------- 239 IF( nidx > 0 ) THEN 240 241 !----------- 242 ! 2D <==> 1D 243 !----------- 244 ! fields used but not modified 245 CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d(1:nidx), sss_m(:,:) ) 246 CALL tab_2d_1d( nidx, idxice(1:nidx), sst_1d(1:nidx), sst_m(:,:) ) 247 ! the following fields are modified in this routine 248 !!CALL tab_2d_1d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) 249 !!CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) 250 !!CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) 251 CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) ) 252 CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 253 CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) 254 IF ( nn_pnd_scheme > 0 ) THEN 255 CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) 256 CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) 257 ENDIF 258 DO jl = 1, jpl 259 DO jk = 1, nlay_s 260 CALL tab_2d_1d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) 261 END DO 262 DO jk = 1, nlay_i 263 CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 264 END DO 265 END DO 266 CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) ) 267 CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) ) 268 CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) ) 269 CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) ) 270 CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) 271 IF ( nn_pnd_scheme > 0 ) THEN 272 CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) 273 ENDIF 274 209 275 !-----------------------------------------------------------------------------! 210 ! Based on the ITD of ridging and ridged ice, convert the net 211 ! closing rate to a gross closing rate. 212 ! NOTE: 0 < aksum <= 1 213 closing_gross(:,:) = closing_net(:,:) / aksum(:,:) 214 215 ! correction to closing rate and opening if closing rate is excessive 216 !--------------------------------------------------------------------- 217 ! Reduce the closing rate if more than 100% of the open water 218 ! would be removed. Reduce the opening rate proportionately. 219 DO jj = 1, jpj 220 DO ji = 1, jpi 221 zfac = ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 222 IF ( zfac < 0._wp .AND. zfac > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i 223 opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 224 ELSEIF( zfac > 0._wp .AND. zfac > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum 225 opning(ji,jj) = apartf(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 276 ! 3) Ridging iteration 277 !-----------------------------------------------------------------------------! 278 niter = 1 279 iterate_ridging = 1 280 DO WHILE( iterate_ridging > 0 .AND. niter < nitermax ) 281 282 CALL rdgrft_prep( a_i_2d, v_i_2d, ato_i_1d ) 283 284 ! 3.2 Redistribute area, volume, and energy. 285 !-----------------------------------------------------------------------------! 286 CALL rdgrft_shift 287 288 ! 3.4 Do we keep on iterating? 289 !-----------------------------------------------------------------------------! 290 ! Check whether asum = 1. If not (because the closing and opening 291 ! rates were reduced above), ridge again with new rates. 292 iterate_ridging = 0 293 DO ji = 1, nidx 294 zfac = 1._wp - ( ato_i_1d(ji) + SUM( a_i_2d(ji,:) ) ) 295 IF( ABS( zfac ) < epsi10 ) THEN 296 closing_net(ji) = 0._wp 297 opning (ji) = 0._wp 298 ato_i_1d (ji) = MAX( 0._wp, 1._wp - SUM( a_i_2d(ji,:) ) ) 299 ELSE 300 iterate_ridging = 1 301 zdivu_adv (ji) = zfac * r1_rdtice 302 closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) ) 303 opning (ji) = MAX( 0._wp, zdivu_adv(ji) ) 226 304 ENDIF 227 305 END DO 228 END DO 229 230 ! correction to closing rate / opening if excessive ice removal 231 !--------------------------------------------------------------- 232 ! Reduce the closing rate if more than 100% of any ice category 233 ! would be removed. Reduce the opening rate proportionately. 306 ! 307 niter = niter + 1 308 IF( niter > nitermax ) CALL ctl_warn( 'icedyn_rdgrft: non-converging ridging scheme' ) 309 ! 310 END DO 311 312 !----------- 313 ! 1D <==> 2D 314 !----------- 315 CALL tab_1d_2d( nidx, idxice(1:nidx), ato_i_1d(1:nidx), ato_i(:,:) ) 316 CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d(1:nidx,1:jpl), a_i(:,:,:) ) 317 CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i (:,:,:) ) 318 CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s (:,:,:) ) 319 CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) ) 320 CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i (:,:,:) ) 321 IF ( nn_pnd_scheme > 0 ) THEN 322 CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d(1:nidx,1:jpl), a_ip(:,:,:) ) 323 CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d(1:nidx,1:jpl), v_ip(:,:,:) ) 324 ENDIF 234 325 DO jl = 1, jpl 235 DO jj = 1, jpj 236 DO ji = 1, jpi 237 zfac = apartf(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 238 IF( zfac > a_i(ji,jj,jl) ) THEN 239 closing_gross(ji,jj) = closing_gross(ji,jj) * a_i(ji,jj,jl) / zfac 240 ENDIF 241 END DO 242 END DO 243 END DO 244 245 ! 3.2 Redistribute area, volume, and energy. 246 !-----------------------------------------------------------------------------! 247 CALL rdgrft_shift( opning, closing_gross ) 248 249 ! 3.3 Compute total area of ice plus open water after ridging. 250 !-----------------------------------------------------------------------------! 251 ! This is in general not equal to one because of divergence during transport 252 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 253 254 ! 3.4 Do we keep on iterating? 255 !-----------------------------------------------------------------------------! 256 ! Check whether asum = 1. If not (because the closing and opening 257 ! rates were reduced above), ridge again with new rates. 258 iterate_ridging = 0 259 DO jj = 1, jpj 260 DO ji = 1, jpi 261 IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN 262 closing_net(ji,jj) = 0._wp 263 opning (ji,jj) = 0._wp 264 ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 265 ELSE 266 iterate_ridging = 1 267 divu_adv (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice 268 closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 269 opning (ji,jj) = MAX( 0._wp, divu_adv(ji,jj) ) 270 ENDIF 271 END DO 272 END DO 273 IF( lk_mpp ) CALL mpp_max( iterate_ridging ) 274 275 ! Repeat if necessary. 276 ! NOTE: If strength smoothing is turned on, the ridging must be 277 ! iterated globally because of the boundary update in the smoothing. 278 niter = niter + 1 279 ! 280 IF( iterate_ridging == 1 ) THEN 281 CALL rdgrft_prep 282 IF( niter > nitermax ) THEN 283 WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 284 WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 285 ENDIF 326 DO jk = 1, nlay_s 327 CALL tab_1d_2d( nidx, idxice(1:nidx), ze_s_2d(1:nidx,jk,jl), e_s(:,:,jk,jl) ) 328 END DO 329 DO jk = 1, nlay_i 330 CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) ) 331 END DO 332 END DO 333 CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_dyn_1d (1:nidx), sfx_dyn (:,:) ) 334 CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_bri_1d (1:nidx), sfx_bri (:,:) ) 335 CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_dyn_1d (1:nidx), wfx_dyn (:,:) ) 336 CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_dyn_1d (1:nidx), hfx_dyn (:,:) ) 337 CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_snw_dyn_1d(1:nidx), wfx_snw_dyn(:,:) ) 338 IF ( nn_pnd_scheme > 0 ) THEN 339 CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_pnd_1d(1:nidx), wfx_pnd(:,:) ) 286 340 ENDIF 287 ! 288 END DO !! on the do while over iter289 341 342 ENDIF ! nidx > 0 343 290 344 CALL ice_var_agg( 1 ) 291 345 … … 298 352 299 353 300 SUBROUTINE rdgrft_prep 301 !!------------------------------------------------------------------- --!354 SUBROUTINE rdgrft_prep( pa_i, pv_i, pato_i ) 355 !!------------------------------------------------------------------- 302 356 !! *** ROUTINE rdgrft_prep *** 303 357 !! 304 !! ** Purpose : preparation for ridging and strengthcalculations358 !! ** Purpose : preparation for ridging calculations 305 359 !! 306 360 !! ** Method : Compute the thickness distribution of the ice and open water 307 361 !! participating in ridging and of the resulting ridges. 308 !!---------------------------------------------------------------------! 309 INTEGER :: ji, jj, jl ! dummy loop indices 310 REAL(wp) :: z1_gstar, z1_astar, zhmean, zdummy ! local scalar 311 REAL(wp), DIMENSION(jpi,jpj,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 312 !------------------------------------------------------------------------------! 313 314 z1_gstar = 1._wp / rn_gstar 315 z1_astar = 1._wp / rn_astar 316 317 CALL ice_var_zapsmall ! Zero out categories with very small areas 362 !!------------------------------------------------------------------- 363 REAL(wp), DIMENSION(:) , INTENT(in) :: pato_i 364 REAL(wp), DIMENSION(:,:), INTENT(in) :: pa_i, pv_i 365 !! 366 INTEGER :: ji, jl ! dummy loop indices 367 REAL(wp) :: z1_gstar, z1_astar, zhmean, zfac ! local scalar 368 REAL(wp), DIMENSION(jpij) :: zasum, z1_asum, zaksum ! sum of a_i+ato_i and reverse 369 REAL(wp), DIMENSION(jpij,jpl) :: zhi ! ice thickness 370 REAL(wp), DIMENSION(jpij,-1:jpl) :: zGsum ! zGsum(n) = sum of areas in categories 0 to n 371 !-------------------------------------------------------------------- 372 373 z1_gstar = 1._wp / rn_gstar 374 z1_astar = 1._wp / rn_astar 318 375 319 376 ! ! Ice thickness needed for rafting 320 WHERE( a_i(:,:,:) >= epsi20 ) ; ht_i(:,:,:) = v_i (:,:,:) / a_i(:,:,:)321 ELSEWHERE ; ht_i(:,:,:) = 0._wp377 WHERE( pa_i(1:nidx,:) > 0._wp ) ; zhi(1:nidx,:) = pv_i(1:nidx,:) / pa_i(1:nidx,:) 378 ELSEWHERE ; zhi(1:nidx,:) = 0._wp 322 379 END WHERE 323 380 … … 338 395 ! Compute total area of ice plus open water. 339 396 ! This is in general not equal to one because of divergence during transport 340 asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 ) 341 ! 397 zasum(1:nidx) = pato_i(1:nidx) + SUM( pa_i(1:nidx,:), dim=2 ) 398 ! 399 WHERE( zasum(1:nidx) > 0._wp ) ; z1_asum(1:nidx) = 1._wp / zasum(1:nidx) 400 ELSEWHERE ; z1_asum(1:nidx) = 0._wp 401 END WHERE 342 402 ! Compute cumulative thickness distribution function 343 403 ! Compute the cumulative thickness distribution function zGsum, 344 404 ! where zGsum(n) is the fractional area in categories 0 to n. 345 405 ! initial value (in h = 0) equals open water area 346 zGsum( :,:,-1) = 0._wp347 zGsum( :,:,0 ) = ato_i(:,:) / asum(:,:)406 zGsum(1:nidx,-1) = 0._wp 407 zGsum(1:nidx,0 ) = pato_i(1:nidx) * z1_asum(1:nidx) 348 408 DO jl = 1, jpl 349 zGsum( :,:,jl) = ( ato_i(:,:) + SUM( a_i(:,:,1:jl), dim=3 ) ) / asum(:,:)! sum(1:jl) is ok (and not jpl)409 zGsum(1:nidx,jl) = ( pato_i(1:nidx) + SUM( pa_i(1:nidx,1:jl), dim=2 ) ) * z1_asum(1:nidx) ! sum(1:jl) is ok (and not jpl) 350 410 END DO 351 352 411 ! 353 412 IF( ln_partf_lin ) THEN !--- Linear formulation (Thorndike et al., 1975) 354 413 DO jl = 0, jpl 355 DO jj = 1, jpj 356 DO ji = 1, jpi 357 IF ( zGsum(ji,jj,jl) < rn_gstar ) THEN 358 apartf(ji,jj,jl) = z1_gstar * ( zGsum(ji,jj,jl) - zGsum(ji,jj,jl-1) ) * & 359 & ( 2._wp - ( zGsum(ji,jj,jl-1) + zGsum(ji,jj,jl) ) * z1_gstar ) 360 ELSEIF( zGsum(ji,jj,jl-1) < rn_gstar ) THEN 361 apartf(ji,jj,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jj,jl-1) ) * & 362 & ( 2._wp - ( zGsum(ji,jj,jl-1) + rn_gstar ) * z1_gstar ) 363 ELSE 364 apartf(ji,jj,jl) = 0._wp 365 ENDIF 366 END DO 414 DO ji = 1, nidx 415 IF ( zGsum(ji,jl) < rn_gstar ) THEN 416 apartf(ji,jl) = z1_gstar * ( zGsum(ji,jl) - zGsum(ji,jl-1) ) * & 417 & ( 2._wp - ( zGsum(ji,jl-1) + zGsum(ji,jl) ) * z1_gstar ) 418 ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 419 apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & 420 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) 421 ELSE 422 apartf(ji,jl) = 0._wp 423 ENDIF 367 424 END DO 368 425 END DO … … 370 427 ELSEIF( ln_partf_exp ) THEN !--- Exponential, more stable formulation (Lipscomb et al, 2007) 371 428 ! 372 z dummy = 1._wp / ( 1._wp - EXP(-z1_astar) ) ! precompute exponential terms using zGsum as a work array429 zfac = 1._wp / ( 1._wp - EXP(-z1_astar) ) 373 430 DO jl = -1, jpl 374 zGsum(:,:,jl) = EXP( -zGsum(:,:,jl) * z1_astar ) * zdummy 431 DO ji = 1, nidx 432 zGsum(ji,jl) = EXP( -zGsum(ji,jl) * z1_astar ) * zfac 433 END DO 375 434 END DO 376 435 DO jl = 0, jpl 377 apartf(:,:,jl) = zGsum(:,:,jl-1) - zGsum(:,:,jl) 436 DO ji = 1, nidx 437 apartf(ji,jl) = zGsum(ji,jl-1) - zGsum(ji,jl) 438 END DO 378 439 END DO 379 440 ! … … 383 444 IF( ln_rafting .AND. ln_ridging ) THEN !- ridging & rafting 384 445 DO jl = 1, jpl 385 DO jj = 1, jpj 386 DO ji = 1, jpi 387 aridge(ji,jj,jl) = ( 1._wp + TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jj,jl) 388 araft (ji,jj,jl) = apartf(ji,jj,jl) - aridge(ji,jj,jl) 389 END DO 446 DO ji = 1, nidx 447 aridge(ji,jl) = ( 1._wp + TANH ( rn_craft * ( zhi(ji,jl) - rn_hraft ) ) ) * 0.5_wp * apartf(ji,jl) 448 araft (ji,jl) = apartf(ji,jl) - aridge(ji,jl) 390 449 END DO 391 450 END DO 392 451 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN !- ridging alone 393 aridge(:,:,:) = apartf(:,:,1:jpl) 394 araft (:,:,:) = 0._wp 452 DO jl = 1, jpl 453 DO ji = 1, nidx 454 aridge(ji,jl) = apartf(ji,jl) 455 araft (ji,jl) = 0._wp 456 END DO 457 END DO 395 458 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN !- rafting alone 396 aridge(:,:,:) = 0._wp 397 araft (:,:,:) = apartf(:,:,1:jpl) 459 DO jl = 1, jpl 460 DO ji = 1, nidx 461 aridge(ji,jl) = 0._wp 462 araft (ji,jl) = apartf(ji,jl) 463 END DO 464 END DO 398 465 ELSE !- no ridging & no rafting 399 aridge(:,:,:) = 0._wp 400 araft (:,:,:) = 0._wp 466 DO jl = 1, jpl 467 DO ji = 1, nidx 468 aridge(ji,jl) = 0._wp 469 araft (ji,jl) = 0._wp 470 END DO 471 END DO 401 472 ENDIF 402 473 … … 422 493 ! ridging. 423 494 ! 424 ! aksum = net area removed/ total area removed495 ! zaksum = net area removed/ total area removed 425 496 ! where total area removed = area of ice that ridges 426 497 ! net area removed = total area removed - area of new ridges 427 498 !----------------------------------------------------------------- 428 499 429 aksum(:,:) = apartf(:,:,0)430 z dummy= 1._wp / hi_hrft500 zaksum(:) = apartf(:,0) 501 zfac = 1._wp / hi_hrft 431 502 ! Transfer function 432 503 DO jl = 1, jpl !all categories have a specific transfer function 433 DO jj = 1, jpj 434 DO ji = 1, jpi 435 IF ( apartf(ji,jj,jl) > 0._wp ) THEN 436 zhmean = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * hrdg_hi_min ) 437 hrmin (ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( zhmean + ht_i(ji,jj,jl) ) ) 438 hrmax (ji,jj,jl) = 2._wp * zhmean - hrmin(ji,jj,jl) 439 hraft (ji,jj,jl) = ht_i(ji,jj,jl) * zdummy 440 hi_hrdg(ji,jj,jl) = ht_i(ji,jj,jl) / MAX( zhmean, epsi20 ) 441 ! 442 ! Normalization factor : aksum, ensures mass conservation 443 aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - hi_hrdg(ji,jj,jl) ) & 444 & + araft (ji,jj,jl) * ( 1._wp - hi_hrft ) 445 ELSE 446 hrmin (ji,jj,jl) = 0._wp 447 hrmax (ji,jj,jl) = 0._wp 448 hraft (ji,jj,jl) = 0._wp 449 hi_hrdg(ji,jj,jl) = 1._wp 450 ENDIF 451 END DO 504 DO ji = 1, nidx 505 IF ( apartf(ji,jl) > 0._wp ) THEN 506 zhmean = MAX( SQRT( rn_hstar * zhi(ji,jl) ), zhi(ji,jl) * hrdg_hi_min ) 507 hrmin (ji,jl) = MIN( 2._wp * zhi(ji,jl), 0.5_wp * ( zhmean + zhi(ji,jl) ) ) 508 hrmax (ji,jl) = 2._wp * zhmean - hrmin(ji,jl) 509 hraft (ji,jl) = zhi(ji,jl) * zfac 510 hi_hrdg(ji,jl) = zhi(ji,jl) / MAX( zhmean, epsi20 ) 511 ! 512 ! Normalization factor : zaksum, ensures mass conservation 513 zaksum(ji) = zaksum(ji) + aridge(ji,jl) * ( 1._wp - hi_hrdg(ji,jl) ) & 514 & + araft (ji,jl) * ( 1._wp - hi_hrft ) 515 ELSE 516 hrmin (ji,jl) = 0._wp 517 hrmax (ji,jl) = 0._wp 518 hraft (ji,jl) = 0._wp 519 hi_hrdg(ji,jl) = 1._wp 520 ENDIF 452 521 END DO 453 522 END DO 454 523 ! 524 ! 3.1 closing_gross 525 !-----------------------------------------------------------------------------! 526 ! Based on the ITD of ridging and ridged ice, convert the net 527 ! closing rate to a gross closing rate. 528 ! NOTE: 0 < aksum <= 1 529 WHERE( zaksum(1:nidx) > 0._wp ) ; closing_gross(1:nidx) = closing_net(1:nidx) / zaksum(1:nidx) 530 ELSEWHERE ; closing_gross(1:nidx) = 0._wp 531 END WHERE 532 533 DO ji = 1, nidx 534 535 ! correction to closing rate and opening if closing rate is excessive 536 !--------------------------------------------------------------------- 537 ! Reduce the closing rate if more than 100% of the open water 538 ! would be removed. Reduce the opening rate proportionately. 539 zfac = ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice 540 IF ( zfac < 0._wp .AND. zfac > - pato_i(ji) ) THEN ! would lead to negative ato_i 541 opning(ji) = apartf(ji,0) * closing_gross(ji) - pato_i(ji) * r1_rdtice 542 ELSEIF( zfac > 0._wp .AND. zfac > ( zasum(ji) - pato_i(ji) ) ) THEN ! would lead to ato_i > asum 543 opning(ji) = apartf(ji,0) * closing_gross(ji) + ( zasum(ji) - pato_i(ji) ) * r1_rdtice 544 ENDIF 545 END DO 546 547 ! correction to closing rate / opening if excessive ice removal 548 !--------------------------------------------------------------- 549 ! Reduce the closing rate if more than 100% of any ice category 550 ! would be removed. Reduce the opening rate proportionately. 551 DO jl = 1, jpl 552 DO ji = 1, nidx 553 zfac = apartf(ji,jl) * closing_gross(ji) * rdt_ice 554 !! IF( zfac > pa_i(ji,jl) .AND. zfac > 0._wp ) THEN 555 IF( zfac > pa_i(ji,jl) ) THEN 556 closing_gross(ji) = closing_gross(ji) * pa_i(ji,jl) / zfac 557 ENDIF 558 END DO 559 END DO 560 ! 455 561 END SUBROUTINE rdgrft_prep 456 562 457 563 458 SUBROUTINE rdgrft_shift ( opning, closing_gross )459 !!------------------------------------------------------------------- ---564 SUBROUTINE rdgrft_shift 565 !!------------------------------------------------------------------- 460 566 !! *** ROUTINE rdgrft_shift *** 461 567 !! … … 464 570 !! ** Method : Remove area, volume, and energy from each ridging category 465 571 !! and add to thicker ice categories. 466 !!---------------------------------------------------------------------- 467 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: opning ! rate of opening due to divergence/shear 468 REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: closing_gross ! rate at which area retreats, excluding area of new ridges 572 !!------------------------------------------------------------------- 469 573 ! 470 574 INTEGER :: ji, jj, jl, jl1, jl2, jk ! dummy loop indices 471 INTEGER :: ij ! horizontal index, combines i and j loops472 INTEGER :: icells ! number of cells with a_i > puny473 575 REAL(wp) :: hL, hR, farea ! left and right limits of integration and new area going to jl2 474 475 INTEGER , DIMENSION(jpij) :: indxi, indxj ! compressed indices 576 REAL(wp) :: vsw ! vol of water trapped into ridges 577 REAL(wp) :: afrdg, afrft ! fraction of category area ridged/rafted 578 REAL(wp) :: airdg1, oirdg1, aprdg1, virdg1, sirdg1 579 REAL(wp) :: airft1, oirft1, aprft1 580 REAL(wp), DIMENSION(jpij) :: airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, esrdg ! area etc of new ridges 581 REAL(wp), DIMENSION(jpij) :: airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, esrft ! area etc of rafted ice 582 ! 583 REAL(wp), DIMENSION(jpij) :: ersw ! enth of water trapped into ridges 476 584 REAL(wp), DIMENSION(jpij) :: zswitch, fvol ! new ridge volume going to jl2 477 478 REAL(wp), DIMENSION(jpij) :: afrac ! fraction of category area ridged 479 REAL(wp), DIMENSION(jpij) :: ardg1 , ardg2 ! area of ice ridged & new ridges 480 REAL(wp), DIMENSION(jpij) :: vsrdg , esrdg ! snow volume & energy of ridging ice 481 ! MV MP 2016 482 REAL(wp), DIMENSION(jpij) :: vprdg ! pond volume of ridging ice 483 REAL(wp), DIMENSION(jpij) :: aprdg1 ! pond area of ridging ice 484 REAL(wp), DIMENSION(jpij) :: aprdg2 ! pond area of ridging ice 485 ! END MV MP 2016 486 REAL(wp), DIMENSION(jpij) :: dhr, dhr2 ! hrmax - hrmin & hrmax^2 - hrmin^2 487 488 REAL(wp), DIMENSION(jpij) :: vrdg1 ! volume of ice ridged 489 REAL(wp), DIMENSION(jpij) :: vrdg2 ! volume of new ridges 490 REAL(wp), DIMENSION(jpij) :: vsw ! volume of seawater trapped into ridges 491 REAL(wp), DIMENSION(jpij) :: srdg1 ! sal*volume of ice ridged 492 REAL(wp), DIMENSION(jpij) :: srdg2 ! sal*volume of new ridges 493 REAL(wp), DIMENSION(jpij) :: smsw ! sal*volume of water trapped into ridges 494 REAL(wp), DIMENSION(jpij) :: oirdg1, oirdg2 ! ice age of ice ridged 495 496 REAL(wp), DIMENSION(jpij) :: afrft ! fraction of category area rafted 497 REAL(wp), DIMENSION(jpij) :: arft1 , arft2 ! area of ice rafted and new rafted zone 498 REAL(wp), DIMENSION(jpij) :: virft , vsrft ! ice & snow volume of rafting ice 499 ! MV MP 2016 500 REAL(wp), DIMENSION(jpij) :: vprft ! pond volume of rafting ice 501 REAL(wp), DIMENSION(jpij) :: aprft1 ! pond area of rafted ice 502 REAL(wp), DIMENSION(jpij) :: aprft2 ! pond area of new rafted ice 503 ! END MV MP 2016 504 REAL(wp), DIMENSION(jpij) :: esrft , smrft ! snow energy & salinity of rafting ice 505 REAL(wp), DIMENSION(jpij) :: oirft1, oirft2 ! ice age of ice rafted 506 585 ! 586 REAL(wp), DIMENSION(jpij,jpl) :: z1_ai ! 1 / a 587 ! 507 588 REAL(wp), DIMENSION(jpij,nlay_i) :: eirft ! ice energy of rafting ice 508 REAL(wp), DIMENSION(jpij,nlay_i) :: erdg1 ! enth*volume of ice ridged 509 REAL(wp), DIMENSION(jpij,nlay_i) :: erdg2 ! enth*volume of new ridges 510 REAL(wp), DIMENSION(jpij,nlay_i) :: ersw ! enth of water trapped into ridges 511 !!---------------------------------------------------------------------- 589 REAL(wp), DIMENSION(jpij,nlay_i) :: eirdg ! enth*volume of new ridges 590 !!------------------------------------------------------------------- 512 591 513 592 !------------------------------------------------------------------------------- 514 593 ! 1) Compute change in open water area due to closing and opening. 515 594 !------------------------------------------------------------------------------- 516 DO jj = 1, jpj 517 DO ji = 1, jpi 518 ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) + & 519 & ( opning(ji,jj) - apartf(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice ) 520 END DO 595 DO ji = 1, nidx 596 ato_i_1d(ji) = MAX( 0._wp, ato_i_1d(ji) + ( opning(ji) - apartf(ji,0) * closing_gross(ji) ) * rdt_ice ) 521 597 END DO 522 598 523 599 !----------------------------------------------------------------- 524 ! 2) Pump everything from ice which is being ridged / rafted600 ! 2) compute categories participating to ridging/rafting (jl1) 525 601 !----------------------------------------------------------------- 526 ! Compute the area, volume, and energy of ice ridging in each 527 ! category, along with the area of the resulting ridge. 528 529 DO jl1 = 1, jpl !jl1 describes the ridging category 530 531 !------------------------------------------------ 532 ! 2.1) Identify grid cells with nonzero ridging 533 !------------------------------------------------ 534 icells = 0 535 DO jj = 1, jpj 536 DO ji = 1, jpi 537 IF( apartf(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN 538 icells = icells + 1 539 indxi(icells) = ji 540 indxj(icells) = jj 602 DO jl1 = 1, jpl 603 604 CALL tab_2d_1d( nidx, idxice(1:nidx), sm_i_1d(1:nidx), sm_i(:,:,jl1) ) 605 606 DO ji = 1, nidx 607 608 !------------------------------------------------ 609 ! 2.1) Identify grid cells with nonzero ridging 610 !------------------------------------------------ 611 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 612 613 z1_ai(ji,jl1) = 1._wp / a_i_2d(ji,jl1) 614 615 !-------------------------------------------------------------------- 616 ! 2.2) Compute area of ridging ice (airdg1) and of new ridge (airdg2) 617 !-------------------------------------------------------------------- 618 airdg1 = aridge(ji,jl1) * closing_gross(ji) * rdt_ice 619 airft1 = araft (ji,jl1) * closing_gross(ji) * rdt_ice 620 621 airdg2(ji) = airdg1 * hi_hrdg(ji,jl1) 622 airft2(ji) = airft1 * hi_hrft 623 624 !--------------------------------------------------------------- 625 ! 2.3) Compute ridging /rafting fractions, make sure afrdg <=1 626 !--------------------------------------------------------------- 627 afrdg = airdg1 * z1_ai(ji,jl1) 628 afrft = airft1 * z1_ai(ji,jl1) 629 630 ! volume and enthalpy (J/m2, >0) of seawater trapped into ridges 631 vsw = v_i_2d(ji,jl1) * afrdg * rn_porordg 632 ersw(ji) = -rhoic * vsw * rcp * sst_1d(ji) ! clem: if sst>0, then ersw <0 (is that possible?) 633 634 !--------------------------------------------------------------- 635 ! 2.4) Compute ridging ice and new ridges (vi, vs, sm, oi, es, ei) 636 !--------------------------------------------------------------- 637 virdg1 = v_i_2d (ji,jl1) * afrdg 638 virdg2(ji) = v_i_2d (ji,jl1) * afrdg * ( 1. + rn_porordg ) 639 vsrdg(ji) = v_s_2d (ji,jl1) * afrdg 640 sirdg1 = smv_i_2d(ji,jl1) * afrdg 641 sirdg2(ji) = smv_i_2d(ji,jl1) * afrdg + vsw * sss_1d(ji) 642 oirdg1 = oa_i_2d (ji,jl1) * afrdg 643 oirdg2(ji) = oa_i_2d (ji,jl1) * afrdg * hi_hrdg(ji,jl1) 644 esrdg(ji) = ze_s_2d (ji,1,jl1) * afrdg 645 646 virft(ji) = v_i_2d (ji,jl1) * afrft 647 vsrft(ji) = v_s_2d (ji,jl1) * afrft 648 sirft(ji) = smv_i_2d(ji,jl1) * afrft 649 oirft1 = oa_i_2d (ji,jl1) * afrft 650 oirft2(ji) = oa_i_2d (ji,jl1) * afrft * hi_hrft 651 esrft(ji) = ze_s_2d (ji,1,jl1) * afrft 652 653 !MV MP 2016 654 IF ( nn_pnd_scheme > 0 ) THEN 655 aprdg1 = a_ip_2d(ji,jl1) * afrdg 656 aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 657 vprdg (ji) = v_ip_2d(ji,jl1) * afrdg 658 aprft1 = a_ip_2d(ji,jl1) * afrft 659 aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 660 vprft (ji) = v_ip_2d(ji,jl1) * afrft 541 661 ENDIF 542 END DO 543 END DO 544 545 DO ij = 1, icells 546 ji = indxi(ij) ; jj = indxj(ij) 547 548 !-------------------------------------------------------------------- 549 ! 2.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 550 !-------------------------------------------------------------------- 551 ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 552 arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice 553 554 !--------------------------------------------------------------- 555 ! 2.3) Compute ridging /rafting fractions, make sure afrac <=1 556 !--------------------------------------------------------------- 557 afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging 558 afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting 559 ardg2(ij) = ardg1(ij) * hi_hrdg(ji,jj,jl1) 560 arft2(ij) = arft1(ij) * hi_hrft 561 562 !-------------------------------------------------------------------------- 563 ! 2.4) Substract area, volume, and energy from ridging 564 ! / rafting category n1. 565 !-------------------------------------------------------------------------- 566 vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij) 567 vrdg2(ij) = vrdg1(ij) * ( 1. + rn_porordg ) 568 vsw (ij) = vrdg1(ij) * rn_porordg 569 570 vsrdg (ij) = v_s (ji,jj, jl1) * afrac(ij) 571 esrdg (ij) = e_s (ji,jj,1,jl1) * afrac(ij) 572 !MV MP 2016 573 IF ( nn_pnd_scheme > 0 ) THEN 574 aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij) 575 aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1) 576 vprdg(ij) = v_ip(ji,jj, jl1) * afrac(ij) 662 ! END MV MP 2016 663 664 !----------------------------------------------------------------- 665 ! 2.5) Ice-ocean exchanges 666 !----------------------------------------------------------------- 667 wfx_dyn_1d(ji) = wfx_dyn_1d(ji) - vsw * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids 668 sfx_dyn_1d(ji) = sfx_dyn_1d(ji) - vsw * sss_1d(ji) * rhoic * r1_rdtice 669 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ersw(ji) * r1_rdtice ! > 0 [W.m-2] 670 671 ! Put the snow lost by ridging into the ocean 672 !------------------------------------------------ 673 ! Note that esrdg > 0; the ocean must cool to melt snow. If the ocean temp = Tf already, new ice must grow. 674 wfx_snw_dyn_1d(ji) = wfx_snw_dyn_1d(ji) + ( rhosn * vsrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! fresh water source for ocean 675 & + rhosn * vsrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 676 677 hfx_dyn_1d(ji) = hfx_dyn_1d(ji) + ( - esrdg(ji) * ( 1._wp - rn_fsnwrdg ) & ! heat sink for ocean (<0, W.m-2) 678 & - esrft(ji) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice 679 680 ! Put the melt pond water into the ocean 681 !------------------------------------------ 682 IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 683 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg ) & ! fresh water source for ocean 684 & + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 685 ENDIF 686 687 ! virtual salt flux to keep salinity constant 688 !------------------------------------------ 689 IF( nn_icesal /= 2 ) THEN 690 sirdg2(ji) = sirdg2(ji) - vsw * ( sss_1d(ji) - sm_i_1d(ji) ) ! ridge salinity = sm_i 691 sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoic * r1_rdtice & ! put back sss_m into the ocean 692 & - sm_i_1d(ji) * vsw * rhoic * r1_rdtice ! and get sm_i from the ocean 693 ENDIF 694 695 696 !------------------------------------------ 697 ! 2.6) update jl1 (removing ridged/rafted area) 698 !------------------------------------------ 699 a_i_2d (ji,jl1) = a_i_2d (ji,jl1) - airdg1 - airft1 700 v_i_2d (ji,jl1) = v_i_2d (ji,jl1) - virdg1 - virft(ji) 701 v_s_2d (ji,jl1) = v_s_2d (ji,jl1) - vsrdg(ji) - vsrft(ji) 702 smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - sirdg1 - sirft(ji) 703 oa_i_2d (ji,jl1) = oa_i_2d (ji,jl1) - oirdg1 - oirft1 704 ! MV MP 2016 705 IF ( nn_pnd_scheme > 0 ) THEN 706 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 707 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 708 ENDIF 709 ! END MV MP 2016 710 ze_s_2d(ji,1,jl1) = ze_s_2d(ji,1,jl1) - esrdg (ji) - esrft (ji) 577 711 ENDIF 578 ! END MV MP 2016 579 srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij) 580 oirdg1(ij) = oa_i (ji,jj, jl1) * afrac(ij) 581 oirdg2(ij) = oa_i (ji,jj, jl1) * afrac(ij) * hi_hrdg(ji,jj,jl1) 582 583 ! rafting volumes, heat contents ... 584 virft (ij) = v_i (ji,jj, jl1) * afrft(ij) 585 vsrft (ij) = v_s (ji,jj, jl1) * afrft(ij) 586 !MV MP 2016 587 IF ( nn_pnd_scheme > 0 ) THEN 588 aprft1(ij) = a_ip (ji,jj, jl1) * afrft(ij) 589 aprft2(ij) = a_ip (ji,jj, jl1) * afrft(ij) * hi_hrft 590 vprft(ij) = v_ip(ji,jj,jl1) * afrft(ij) 591 ENDIF 592 ! END MV MP 2016 593 srdg1 (ij) = smv_i(ji,jj, jl1) * afrac(ij) 594 esrft (ij) = e_s (ji,jj,1,jl1) * afrft(ij) 595 smrft (ij) = smv_i(ji,jj, jl1) * afrft(ij) 596 oirft1(ij) = oa_i (ji,jj, jl1) * afrft(ij) 597 oirft2(ij) = oa_i (ji,jj, jl1) * afrft(ij) * hi_hrft 598 599 !----------------------------------------------------------------- 600 ! 2.5) Compute properties of new ridges 601 !----------------------------------------------------------------- 602 smsw(ij) = vsw(ij) * sss_m(ji,jj) ! salt content of seawater frozen in voids 603 srdg2(ij) = srdg1(ij) + smsw(ij) ! salt content of new ridge 604 605 sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice 606 wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice ! increase in ice volume due to seawater frozen in voids 607 608 ! virtual salt flux to keep salinity constant 609 IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN 610 srdg2(ij) = srdg2(ij) - vsw(ij) * ( sss_m(ji,jj) - sm_i(ji,jj,jl1) ) ! ridge salinity = sm_i 611 sfx_bri(ji,jj) = sfx_bri(ji,jj) + sss_m(ji,jj) * vsw(ij) * rhoic * r1_rdtice & ! put back sss_m into the ocean 612 & - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean 613 ENDIF 614 615 !------------------------------------------ 616 ! 2.6 Put the snow somewhere in the ocean 617 !------------------------------------------ 618 ! Place part of the snow lost by ridging into the ocean. 619 ! Note that esrdg > 0; the ocean must cool to melt snow. 620 ! If the ocean temp = Tf already, new ice must grow. 621 ! During the next time step, thermo_rates will determine whether 622 ! the ocean cools or new ice grows. 623 wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnwrdg ) & 624 & + rhosn * vsrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! fresh water source for ocean 625 626 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnwrdg ) & 627 & - esrft(ij) * ( 1._wp - rn_fsnwrft ) ) * r1_rdtice ! heat sink for ocean (<0, W.m-2) 628 629 ! MV MP 2016 630 !------------------------------------------ 631 ! 2.7 Put the melt pond water in the ocean 632 !------------------------------------------ 633 ! Place part of the melt pond volume into the ocean. 634 IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 635 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpndrdg ) & 636 & + rhofw * vprft(ij) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice ! fresh water source for ocean 637 ENDIF 638 ! END MV MP 2016 639 640 !----------------------------------------------------------------- 641 ! 2.8 Compute quantities used to apportion ice among categories 642 ! in the n2 loop below 643 !----------------------------------------------------------------- 644 dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) ) 645 dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) ) 646 647 648 ! update jl1 (removing ridged/rafted area) 649 a_i (ji,jj, jl1) = a_i (ji,jj, jl1) - ardg1 (ij) - arft1 (ij) 650 v_i (ji,jj, jl1) = v_i (ji,jj, jl1) - vrdg1 (ij) - virft (ij) 651 v_s (ji,jj, jl1) = v_s (ji,jj, jl1) - vsrdg (ij) - vsrft (ij) 652 e_s (ji,jj,1,jl1) = e_s (ji,jj,1,jl1) - esrdg (ij) - esrft (ij) 653 smv_i(ji,jj, jl1) = smv_i(ji,jj, jl1) - srdg1 (ij) - smrft (ij) 654 oa_i (ji,jj, jl1) = oa_i (ji,jj, jl1) - oirdg1(ij) - oirft1(ij) 655 656 ! MV MP 2016 657 IF ( nn_pnd_scheme > 0 ) THEN 658 v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij) 659 a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij) 660 ENDIF 661 ! END MV MP 2016 662 663 END DO 664 665 !-------------------------------------------------------------------- 666 ! 2.9 Compute ridging ice enthalpy, remove it from ridging ice and 667 ! compute ridged ice enthalpy 668 !-------------------------------------------------------------------- 712 713 END DO ! ji 714 715 ! special loop for e_i because of layers jk 669 716 DO jk = 1, nlay_i 670 DO ij = 1, icells 671 ji = indxi(ij) ; jj = indxj(ij) 672 ! heat content of ridged ice 673 erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 674 eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij) 675 676 ! enthalpy of the trapped seawater (J/m2, >0) 677 ! clem: if sst>0, then ersw <0 (is that possible?) 678 ersw(ij,jk) = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i 679 680 ! heat flux to the ocean 681 hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice ! > 0 [W.m-2] ocean->ice flux 682 683 ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 684 erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk) 685 686 ! update jl1 687 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk) 688 689 END DO 690 END DO 717 DO ji = 1, nidx 718 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 719 ! Compute ridging /rafting fractions 720 afrdg = aridge(ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) 721 afrft = araft (ji,jl1) * closing_gross(ji) * rdt_ice * z1_ai(ji,jl1) 722 ! Compute ridging ice and new ridges for ei 723 eirdg(ji,jk) = ze_i_2d (ji,jk,jl1) * afrdg + ersw(ji) * r1_nlay_i 724 eirft(ji,jk) = ze_i_2d (ji,jk,jl1) * afrft 725 ! Update jl1 726 ze_i_2d(ji,jk,jl1) = ze_i_2d(ji,jk,jl1) * ( 1._wp - afrdg - afrft ) 727 ENDIF 728 END DO 729 END DO 730 691 731 692 732 !------------------------------------------------------------------------------- … … 694 734 !------------------------------------------------------------------------------- 695 735 DO jl2 = 1, jpl 696 ! over categories to which ridged/rafted ice is transferred 697 DO ij = 1, icells 698 ji = indxi(ij) ; jj = indxj(ij) 699 700 ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 701 IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN 702 hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) ) 703 hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2) ) 704 farea = ( hR - hL ) * dhr(ij) 705 fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij) 706 ELSE 707 farea = 0._wp 708 fvol(ij) = 0._wp 736 ! 737 DO ji = 1, nidx 738 739 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN 740 741 ! Compute the fraction of ridged ice area and volume going to thickness category jl2. 742 IF( hrmin(ji,jl1) <= hi_max(jl2) .AND. hrmax(ji,jl1) > hi_max(jl2-1) ) THEN 743 !----------------------------------------------------------------- 744 ! 2.8 Compute quantities used to apportion ice among categories 745 ! in the n2 loop below 746 !----------------------------------------------------------------- 747 hL = MAX( hrmin(ji,jl1), hi_max(jl2-1) ) 748 hR = MIN( hrmax(ji,jl1), hi_max(jl2) ) 749 farea = ( hR - hL ) / ( hrmax(ji,jl1) - hrmin(ji,jl1) ) 750 fvol(ji) = ( hR * hR - hL * hL ) / ( hrmax(ji,jl1) * hrmax(ji,jl1) - hrmin(ji,jl1) * hrmin(ji,jl1) ) 751 ELSE 752 farea = 0._wp 753 fvol(ji) = 0._wp 754 ENDIF 755 756 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 757 !!gm see above IF( hraft(ji) <= hi_max(jl2) .AND. hraft(ji) > hi_max(jl2-1) ) THEN 758 IF( hi_max(jl2-1) < hraft(ji,jl1) .AND. hraft(ji,jl1) <= hi_max(jl2) ) THEN ; zswitch(ji) = 1._wp 759 ELSE ; zswitch(ji) = 0._wp 760 ENDIF 761 ! 762 a_i_2d (ji,jl2) = a_i_2d (ji,jl2) + ( airdg2(ji) * farea + airft2(ji) * zswitch(ji) ) 763 oa_i_2d (ji,jl2) = oa_i_2d (ji,jl2) + ( oirdg2(ji) * farea + oirft2(ji) * zswitch(ji) ) 764 v_i_2d (ji,jl2) = v_i_2d (ji,jl2) + ( virdg2(ji) * fvol(ji) + virft (ji) * zswitch(ji) ) 765 smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ( sirdg2(ji) * fvol(ji) + sirft (ji) * zswitch(ji) ) 766 v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + & 767 & vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 768 ! MV MP 2016 769 IF ( nn_pnd_scheme > 0 ) THEN 770 v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) & 771 & + vprft (ji) * rn_fpndrft * zswitch(ji) ) 772 a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + ( aprdg2(ji) * rn_fpndrdg * farea & 773 & + aprft2(ji) * rn_fpndrft * zswitch(ji) ) 774 ENDIF 775 ! END MV MP 2016 776 ze_s_2d(ji,1,jl2) = ze_s_2d(ji,1,jl2) + ( esrdg (ji) * rn_fsnwrdg * fvol(ji) + & 777 & esrft (ji) * rn_fsnwrft * zswitch(ji) ) 778 709 779 ENDIF 710 780 711 ! Compute the fraction of rafted ice area and volume going to thickness category jl2 712 !!gm see above IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1) ) THEN 713 IF( hi_max(jl2-1) < hraft(ji,jj,jl1) .AND. hraft(ji,jj,jl1) <= hi_max(jl2) ) THEN ; zswitch(ij) = 1._wp 714 ELSE ; zswitch(ij) = 0._wp 715 ENDIF 716 ! 717 a_i (ji,jj ,jl2) = a_i (ji,jj ,jl2) + ( ardg2 (ij) * farea + arft2 (ij) * zswitch(ij) ) 718 oa_i (ji,jj ,jl2) = oa_i (ji,jj ,jl2) + ( oirdg2(ij) * farea + oirft2(ij) * zswitch(ij) ) 719 v_i (ji,jj ,jl2) = v_i (ji,jj ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) ) 720 smv_i(ji,jj ,jl2) = smv_i(ji,jj ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) ) 721 v_s (ji,jj ,jl2) = v_s (ji,jj ,jl2) + ( vsrdg (ij) * rn_fsnwrdg * fvol(ij) + & 722 & vsrft (ij) * rn_fsnwrft * zswitch(ij) ) 723 e_s (ji,jj,1,jl2) = e_s (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnwrdg * fvol(ij) + & 724 & esrft (ij) * rn_fsnwrft * zswitch(ij) ) 725 ! MV MP 2016 726 IF ( nn_pnd_scheme > 0 ) THEN 727 v_ip (ji,jj,jl2) = v_ip(ji,jj,jl2) + ( vprdg (ij) * rn_fpndrdg * fvol (ij) & 728 & + vprft (ij) * rn_fpndrft * zswitch(ij) ) 729 a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2) + ( aprdg2(ij) * rn_fpndrdg * farea & 730 & + aprft2(ij) * rn_fpndrft * zswitch(ji) ) 731 ENDIF 732 ! END MV MP 2016 733 END DO 734 735 ! Transfer ice energy to category jl2 by ridging 781 END DO 782 736 783 DO jk = 1, nlay_i 737 DO ij = 1, icells738 ji = indxi(ij) ; jj = indxj(ij)739 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)784 DO ji = 1, nidx 785 IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) & 786 & ze_i_2d(ji,jk,jl2) = ze_i_2d(ji,jk,jl2) + eirdg(ji,jk) * fvol(ji) + eirft(ji,jk) * zswitch(ji) 740 787 END DO 741 788 END DO … … 743 790 END DO ! jl2 744 791 ! 745 END DO ! jl1 (deforming categories) 792 END DO ! jl1 793 ! 794 ! In case ridging/rafting lead to very small negative values (sometimes it happens) 795 WHERE( a_i_2d(1:nidx,:) < 0._wp ) a_i_2d(1:nidx,:) = 0._wp 796 WHERE( v_i_2d(1:nidx,:) < 0._wp ) v_i_2d(1:nidx,:) = 0._wp 746 797 ! 747 798 END SUBROUTINE rdgrft_shift … … 760 811 !! ice strength has to be smoothed 761 812 !!---------------------------------------------------------------------- 762 INTEGER :: ji, jj, jl! dummy loop indices813 INTEGER :: ji, jj, jl ! dummy loop indices 763 814 INTEGER :: ismooth ! smoothing the resistance to deformation 764 815 INTEGER :: itframe ! number of time steps for the P smoothing … … 767 818 REAL(wp), DIMENSION(jpi,jpj) :: zstrp1, zstrp2 ! strength at previous time steps 768 819 !!---------------------------------------------------------------------- 769 770 ! !--------------------------------------------------!771 CALL rdgrft_prep ! Thickness distribution of ridging and ridged ice !772 ! !--------------------------------------------------!773 774 820 ! !--------------------------------------------------! 775 821 IF( ln_str_R75 ) THEN ! Ice strength => Rothrock (1975) method ! 776 822 ! !--------------------------------------------------! 777 z1_3 = 1._wp / 3._wp778 DO jl = 1, jpl779 WHERE( apartf(:,:,jl) > 0._wp )780 strength(:,:) = - apartf(:,:,jl) * ht_i (:,:,jl) * ht_i(:,:,jl) & ! PE loss from deforming ice781 & + 2._wp * araft (:,:,jl) * ht_i (:,:,jl) * ht_i(:,:,jl) & ! PE gain from rafting ice782 & + aridge(:,:,jl) * hi_hrdg(:,:,jl) * z1_3 * & ! PE gain from ridging ice783 & ( hrmax(:,:,jl) * hrmax (:,:,jl) + &784 & hrmin(:,:,jl) * hrmin (:,:,jl) + &785 & hrmax(:,:,jl) * hrmin (:,:,jl) )786 ELSEWHERE787 strength(:,:) = 0._wp788 END WHERE789 END DO790 strength(:,:) = rn_perdg * zdrho * strength(:,:) / aksum(:,:) * tmask(:,:,1)791 ! where zdrho = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_perdg accounts for frictional dissipation792 ismooth = 1793 823 ! !--------------------------------------------------! 794 824 ELSEIF( ln_str_H79 ) THEN ! Ice strength => Hibler (1979) method ! … … 805 835 DO jj = 2, jpjm1 806 836 DO ji = 2, jpim1 807 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN837 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 808 838 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 809 839 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & … … 831 861 DO jj = 2, jpjm1 832 862 DO ji = 2, jpim1 833 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN863 IF ( SUM( a_i(ji,jj,:) ) > 0._wp ) THEN 834 864 itframe = 1 ! number of time steps for the running mean 835 865 IF ( zstrp1(ji,jj) > 0._wp ) itframe = itframe + 1 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icetab.F90
r8534 r8554 13 13 !! tab_1d_2d : 1-D <==> 2-D 14 14 !!---------------------------------------------------------------------- 15 USE par_kind16 15 USE par_oce 17 16 USE ice, ONLY : jpl
Note: See TracChangeset
for help on using the changeset viewer.