Changeset 8361
 Timestamp:
 20170723T12:51:10+02:00 (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r8355 r8361 15 15 !! lim_itd_th_rem : 16 16 !! lim_itd_th_reb : 17 !! lim_itd_ fitline:17 !! lim_itd_glinear : 18 18 !! lim_itd_shiftice : 19 19 !! … … 24 24 USE ice ! LIM3 variables 25 25 USE limvar ! LIM3 variables 26 USE limcons ! conservation tests 27 ! 26 28 USE prtctl ! Print control 27 29 USE in_out_manager ! I/O manager … … 29 31 USE wrk_nemo ! work arrays 30 32 USE lib_fortran ! to use key_nosignedzero 31 USE limcons ! conservation tests32 33 33 34 IMPLICIT NONE … … 71 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 72 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness74 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume 75 INTEGER , POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 76 INTEGER :: nbrem ! number of cells with ice to transfer 75 INTEGER , DIMENSION(jpij) :: nind_i, nind_j ! compressed indices for i/j directions 77 76 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 78 77 REAL(wp), POINTER, DIMENSION(:,:) :: zhb0, zhb1 ! category boundaries for thinnes categories … … 83 82 CALL wrk_alloc( jpi,jpj, zremap_flag ) 84 83 CALL wrk_alloc( jpi,jpj,jpl1, zdonor ) 85 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR , zht_i_b)84 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 86 85 CALL wrk_alloc( jpi,jpj,jpl1, zdaice, zdvice ) 87 86 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 88 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j )89 87 CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 90 88 91 !!92 !! 1) Compute thickness and changes in each ice category93 !!94 89 IF( kt == nit000 .AND. lwp) THEN 95 90 WRITE(numout,*) … … 98 93 ENDIF 99 94 100 zdhice(:,:,:) = 0._wp101 DO jl = 1, jpl102 DO jj = 1, jpj103 DO ji = 1, jpi104 rswitch = MAX( 0.0, SIGN( 1.0, a_i(ji,jj,jl)  epsi10 ) ) !0 if no ice and 1 if yes105 ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch106 rswitch = MAX( 0.0, SIGN( 1.0, a_i_b(ji,jj,jl)  epsi10) )107 zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch108 IF( a_i(ji,jj,jl) > epsi10 ) zdhice(ji,jj,jl) = ht_i(ji,jj,jl)  zht_i_b(ji,jj,jl) ! clem: useless IF statement?109 END DO110 END DO111 END DO112 113 !! !114 !! ! 2) Compute fractional ice area in each grid cell115 !! !116 !! at_i(:,:) = 0._wp117 !! DO jl = 1, jpl118 !! at_i(:,:) = at_i(:,:) + a_i(:,:,jl)119 !! END DO120 121 95 ! 122 96 ! 3) Identify grid cells with ice 123 97 ! 124 n brem= 098 nidx = 0 ; idxice(:) = 0 125 99 DO jj = 1, jpj 126 100 DO ji = 1, jpi 127 101 IF ( at_i(ji,jj) > epsi10 ) THEN 128 nbrem = nbrem + 1 129 nind_i(nbrem) = ji 130 nind_j(nbrem) = jj 102 nidx = nidx + 1 103 nind_i(nidx) = ji 104 nind_j(nidx) = jj 105 idxice( nidx ) = (jj  1) * jpi + ji 131 106 zremap_flag(ji,jj) = 1 132 107 ELSE … … 139 114 ! 4) Compute new category boundaries: 1:jpl1 140 115 ! 116 zdhice(:,:,:) = 0._wp 141 117 zhbnew(:,:,:) = 0._wp 142 143 IF( nbrem > 0 ) THEN 118 zhb0(:,:) = hi_max(0) 119 zhb1(:,:) = hi_max(1) 120 121 IF( nidx > 0 ) THEN 122 123 ! Compute thickness change in each ice category 124 DO jl = 1, jpl 125 DO ji = 1, nidx 126 ii = nind_i(ji) 127 ij = nind_j(ji) 128 zdhice(ii,ij,jl) = ht_i(ii,ij,jl)  ht_i_b(ii,ij,jl) 129 END DO 130 END DO 131 144 132 DO jl = 1, jpl  1 145 DO ji = 1, n brem133 DO ji = 1, nidx 146 134 ii = nind_i(ji) 147 135 ij = nind_j(ji) 148 136 ! 149 !  New boundary categories Hn* = Hn + Fn*dt  150 ! !Fn*dt = ( fn + (fn+1  fn)/(hn+1  hn) * (Hn  hn) ) * dt = zdhice + zslope * (Hmax  ht_i_b)151 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl+1) & a_i_b(jl) /= 0152 zslope = ( zdhice(ii,ij,jl+1)  zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1)  zht_i_b(ii,ij,jl) )153 z hbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl)  zht_i_b(ii,ij,jl) )154 !! ELSEIF( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,jj,jl+1) <= epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 155 ELSEIF( a_i_b(ii,ij,jl) > epsi10 ) THEN !! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt137 !  New boundary categories Hn* = Hn + Fn*dt  ! 138 ! Fn*dt = ( fn + (fn+1  fn)/(hn+1  hn) * (Hn  hn) ) * dt = zdhice + zslope * (Hmax  ht_i_b) 139 ! 140 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN ! a_i_b(jl+1) & a_i_b(jl) /= 0 141 zslope = ( zdhice(ii,ij,jl+1)  zdhice(ii,ij,jl) ) / ( ht_i_b(ii,ij,jl+1)  ht_i_b(ii,ij,jl) ) 142 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl)  ht_i_b(ii,ij,jl) ) 143 ELSEIF( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) <= epsi10 ) THEN ! a_i_b(jl+1)=0 => Hn* = Hn + fn*dt 156 144 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 157 !! ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 158 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10 ) THEN !! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 145 ELSEIF( a_i_b(ii,ij,jl) <= epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN ! a_i_b(jl)=0 => Hn* = Hn + fn+1*dt 159 146 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 160 ELSE !! a_i_b(jl+1) & a_i_b(jl) = 0147 ELSE ! a_i_b(jl+1) & a_i_b(jl) = 0 161 148 zhbnew(ii,ij,jl) = hi_max(jl) 162 149 ENDIF … … 165 152 ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)epsi 166 153 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 167 ! in lim_itd_ fitlinein the case (HRHL) = 3(Hice  HL) or = 3(HR  Hice)154 ! in lim_itd_glinear in the case (HRHL) = 3(Hice  HL) or = 3(HR  Hice) 168 155 IF( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl)  epsi10 ) ) zremap_flag(ii,ij) = 0 169 156 IF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) zremap_flag(ii,ij) = 0 … … 175 162 END DO 176 163 END DO 164 165 !  New boundary for category jpl  ! 166 DO ji = 1, nidx 167 ii = nind_i(ji) 168 ij = nind_j(ji) 169 IF( a_i(ii,ij,jpl) > epsi10 ) THEN 170 zhbnew(ii,ij,jpl) = MAX( hi_max(jpl1), 3._wp * ht_i(ii,ij,jpl)  2._wp * zhbnew(ii,ij,jpl1) ) 171 ELSE 172 zhbnew(ii,ij,jpl) = hi_max(jpl) 173 ENDIF 174 175 ! 1 condition for remapping for the 1st category 176 ! H0+epsi < h1(t) < H1epsi 177 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 178 ! in lim_itd_glinear in the case (HRHL) = 3(Hice  HL) or = 3(HR  Hice) 179 IF( ht_i_b(ii,ij,1) < ( hi_max(0) + epsi10 ) ) zremap_flag(ii,ij) = 0 180 IF( ht_i_b(ii,ij,1) > ( hi_max(1)  epsi10 ) ) zremap_flag(ii,ij) = 0 181 END DO 182 177 183 ENDIF 178 184 … … 180 186 ! 5) Identify cells where ITD is to be remapped 181 187 ! 182 n brem= 0188 nidx = 0 183 189 DO jj = 1, jpj 184 190 DO ji = 1, jpi 185 191 IF( zremap_flag(ji,jj) == 1 ) THEN 186 n brem = nbrem+ 1187 nind_i(n brem) = ji188 nind_j(n brem) = jj192 nidx = nidx + 1 193 nind_i(nidx) = ji 194 nind_j(nidx) = jj 189 195 ENDIF 190 196 END DO … … 192 198 193 199 ! 194 ! 6) Compute new category boundaries: jpl 195 ! 196 DO jj = 1, jpj 197 DO ji = 1, jpi 198 zhb0(ji,jj) = hi_max(0) 199 zhb1(ji,jj) = hi_max(1) 200 201 ! boundary for jpl 202 IF( a_i(ji,jj,jpl) > epsi10 ) THEN 203 zhbnew(ji,jj,jpl) = MAX( hi_max(jpl1), 3._wp * ht_i(ji,jj,jpl)  2._wp * zhbnew(ji,jj,jpl1) ) 204 ELSE 205 !clem bug zhbnew(ji,jj,jpl) = hi_max(jpl) 206 zhbnew(ji,jj,jpl) = hi_max(jpl1) ! not used anyway 200 ! 7) Compute g(h) 201 ! 202 IF( nidx > 0 ) THEN 203 204 ! 7.1 g(h) for category 1 at start of time step 205 CALL lim_itd_glinear( 1, zhb0, zhb1, ht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1), & 206 & hR(:,:,1), zremap_flag ) 207 208 ! 7.2 Area lost due to melting of thin ice (first category, 1) 209 DO ji = 1, nidx 210 ii = nind_i(ji) 211 ij = nind_j(ji) 212 213 IF( a_i(ii,ij,1) > epsi10 ) THEN 214 215 zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 216 217 IF( zdh0 < 0.0 ) THEN !remove area from category 1 218 zdh0 = MIN( zdh0, hi_max(1) ) 219 !Integrate g(1) from 0 to dh0 to estimate area melted 220 zetamax = MIN( zdh0, hR(ii,ij,1) )  hL(ii,ij,1) 221 222 IF( zetamax > 0.0 ) THEN 223 zx1 = zetamax 224 zx2 = 0.5 * zetamax * zetamax 225 zda0 = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1 ! ice area removed 226 zdamax = a_i(ii,ij,1) * (1.0  ht_i(ii,ij,1) / ht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i 227 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 228 ! of thin ice (zdamax > 0) 229 ! Remove area, conserving volume 230 ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1)  zda0 ) 231 a_i(ii,ij,1) = a_i(ii,ij,1)  zda0 232 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clemuseless ? 233 ENDIF 234 235 ELSE ! if ice accretion zdh0 > 0 236 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 237 zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) ) 238 ENDIF 239 207 240 ENDIF 208 209 ! 1 condition for remapping for the 1st category 210 ! H0+epsi < h1(t) < H1epsi 211 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 212 ! in lim_itd_fitline in the case (HRHL) = 3(Hice  HL) or = 3(HR  Hice) 213 IF( zht_i_b(ji,jj,1) < ( zhb0(ji,jj) + epsi10 ) ) zremap_flag(ji,jj) = 0 214 IF( zht_i_b(ji,jj,1) > ( zhb1(ji,jj)  epsi10 ) ) zremap_flag(ji,jj) = 0 215 216 END DO 217 END DO 218 219 ! 220 ! 7) Compute g(h) 221 ! 222 ! 7.1 g(h) for category 1 at start of time step 223 CALL lim_itd_fitline( 1, zhb0, zhb1, zht_i_b(:,:,1), g0(:,:,1), g1(:,:,1), hL(:,:,1), & 224 & hR(:,:,1), zremap_flag ) 225 226 ! 7.2 Area lost due to melting of thin ice (first category, 1) 227 DO ji = 1, nbrem 228 ii = nind_i(ji) 229 ij = nind_j(ji) 230 231 IF( a_i(ii,ij,1) > epsi10 ) THEN 232 233 zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 234 235 IF( zdh0 < 0.0 ) THEN !remove area from category 1 236 zdh0 = MIN( zdh0, hi_max(1) ) 237 !Integrate g(1) from 0 to dh0 to estimate area melted 238 zetamax = MIN( zdh0, hR(ii,ij,1) )  hL(ii,ij,1) 239 240 IF( zetamax > 0.0 ) THEN 241 zx1 = zetamax 242 zx2 = 0.5 * zetamax * zetamax 243 zda0 = g1(ii,ij,1) * zx2 + g0(ii,ij,1) * zx1 ! ice area removed 244 zdamax = a_i(ii,ij,1) * (1.0  ht_i(ii,ij,1) / zht_i_b(ii,ij,1) ) ! Constrain new thickness <= ht_i 245 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 246 ! of thin ice (zdamax > 0) 247 ! Remove area, conserving volume 248 ht_i(ii,ij,1) = ht_i(ii,ij,1) * a_i(ii,ij,1) / ( a_i(ii,ij,1)  zda0 ) 249 a_i(ii,ij,1) = a_i(ii,ij,1)  zda0 250 v_i(ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) ! clemuseless ? 251 ENDIF 252 253 ELSE ! if ice accretion zdh0 > 0 254 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 255 zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) ) 256 ENDIF 257 258 ENDIF 259 260 END DO 261 262 ! 7.3 g(h) for each thickness category 263 DO jl = 1, jpl 264 CALL lim_itd_fitline( jl, zhbnew(:,:,jl1), zhbnew(:,:,jl), ht_i(:,:,jl), & 265 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 266 END DO 267 268 ! 269 ! 8) Compute area and volume to be shifted across each boundary (Eq. 18) 270 ! 271 272 DO jl = 1, jpl  1 273 DO jj = 1, jpj 274 DO ji = 1, jpi 275 zdonor(ji,jj,jl) = 0 276 zdaice(ji,jj,jl) = 0.0 277 zdvice(ji,jj,jl) = 0.0 278 END DO 279 END DO 280 281 DO ji = 1, nbrem 241 242 END DO 243 244 ! 7.3 g(h) for each thickness category 245 DO jl = 1, jpl 246 CALL lim_itd_glinear( jl, zhbnew(:,:,jl1), zhbnew(:,:,jl), ht_i(:,:,jl), & 247 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 248 END DO 249 250 ! 251 ! 8) Compute area and volume to be shifted across each boundary (Eq. 18) 252 ! 253 254 DO jl = 1, jpl  1 255 DO jj = 1, jpj 256 DO ji = 1, jpi 257 zdonor(ji,jj,jl) = 0 258 zdaice(ji,jj,jl) = 0.0 259 zdvice(ji,jj,jl) = 0.0 260 END DO 261 END DO 262 263 DO ji = 1, nidx 264 ii = nind_i(ji) 265 ij = nind_j(ji) 266 267 ! left and right integration limits in eta space 268 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 269 zetamin = MAX( hi_max(jl), hL(ii,ij,jl) )  hL(ii,ij,jl) ! hi_max(jl)  hL 270 zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) )  hL(ii,ij,jl) ! hR  hL 271 zdonor(ii,ij,jl) = jl 272 ELSE ! Hn* <= Hn => transfer from jl+1 to jl 273 zetamin = 0.0 274 zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) )  hL(ii,ij,jl+1) ! hi_max(jl)  hL 275 zdonor(ii,ij,jl) = jl + 1 276 ENDIF 277 zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 278 279 zx1 = zetamax  zetamin 280 zwk1 = zetamin * zetamin 281 zwk2 = zetamax * zetamax 282 zx2 = 0.5 * ( zwk2  zwk1 ) 283 zwk1 = zwk1 * zetamin 284 zwk2 = zwk2 * zetamax 285 zx3 = 1.0 / 3.0 * ( zwk2  zwk1 ) 286 nd = zdonor(ii,ij,jl) 287 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 288 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 289 290 END DO 291 END DO 292 293 !! 294 !! 9) Shift ice between categories 295 !! 296 CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 297 298 !! 299 !! 10) Make sure ht_i >= minimum ice thickness hi_min 300 !! 301 302 DO ji = 1, nidx 282 303 ii = nind_i(ji) 283 304 ij = nind_j(ji) 284 285 ! left and right integration limits in eta space 286 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 287 zetamin = MAX( hi_max(jl), hL(ii,ij,jl) )  hL(ii,ij,jl) ! hi_max(jl)  hL 288 zetamax = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) )  hL(ii,ij,jl) ! hR  hL 289 zdonor(ii,ij,jl) = jl 290 ELSE ! Hn* <= Hn => transfer from jl+1 to jl 291 zetamin = 0.0 292 zetamax = MIN( hi_max(jl), hR(ii,ij,jl+1) )  hL(ii,ij,jl+1) ! hi_max(jl)  hL 293 zdonor(ii,ij,jl) = jl + 1 305 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 306 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 307 ! MV MP 2016 308 IF ( nn_pnd_scheme > 0 ) THEN 309 a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 310 ENDIF 311 ! END MV MP 2016 312 ht_i(ii,ij,1) = rn_himin 294 313 ENDIF 295 zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 296 297 zx1 = zetamax  zetamin 298 zwk1 = zetamin * zetamin 299 zwk2 = zetamax * zetamax 300 zx2 = 0.5 * ( zwk2  zwk1 ) 301 zwk1 = zwk1 * zetamin 302 zwk2 = zwk2 * zetamax 303 zx3 = 1.0 / 3.0 * ( zwk2  zwk1 ) 304 nd = zdonor(ii,ij,jl) 305 zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 306 zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd) 307 308 END DO 309 END DO 310 311 !! 312 !! 9) Shift ice between categories 313 !! 314 CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 315 316 !! 317 !! 10) Make sure ht_i >= minimum ice thickness hi_min 318 !! 319 320 DO ji = 1, nbrem 321 ii = nind_i(ji) 322 ij = nind_j(ji) 323 IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN 324 a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 325 ! MV MP 2016 326 IF ( nn_pnd_scheme > 0 ) THEN 327 a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 328 ENDIF 329 ! END MV MP 2016 330 ht_i(ii,ij,1) = rn_himin 331 ENDIF 332 END DO 314 END DO 315 316 ENDIF 333 317 334 318 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 335 319 CALL wrk_dealloc( jpi,jpj,jpl1, zdonor ) 336 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR , zht_i_b)320 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 337 321 CALL wrk_dealloc( jpi,jpj,jpl1, zdaice, zdvice ) 338 322 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 339 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j )340 323 CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) 341 324 … … 343 326 344 327 345 SUBROUTINE lim_itd_ fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )346 !! 347 !! *** ROUTINE lim_itd_ fitline***328 SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag ) 329 !! 330 !! *** ROUTINE lim_itd_glinear *** 348 331 !! 349 332 !! ** Purpose : build g(h) satisfying area and volume constraints (Eq. 6 and 9) … … 404 387 END DO 405 388 ! 406 END SUBROUTINE lim_itd_ fitline389 END SUBROUTINE lim_itd_glinear 407 390 408 391 … … 421 404 422 405 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 423 INTEGER :: ii, ij ! indices when changing from 2D1D is done 424 425 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaTsfn 426 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 427 428 REAL(wp) :: zdvsnow, zdesnow ! snow volume and energy transferred 429 REAL(wp) :: zdeice ! ice energy transferred 430 REAL(wp) :: zdsm_vice ! ice salinity times volume transferred 431 REAL(wp) :: zdo_aice ! ice age times volume transferred 432 REAL(wp) :: zdaTsf ! aicen*Tsfcn transferred 433 ! MV MP 2016 434 REAL(wp) :: zdapnd ! pond fraction transferred 435 REAL(wp) :: zdvpnd ! pond volume transferred 436 ! END MV MP 2016 437 438 INTEGER, POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 439 440 INTEGER :: nbrem ! number of cells with ice to transfer 441 !! 442 443 CALL wrk_alloc( jpi,jpj,jpl, zaTsfn ) 444 CALL wrk_alloc( jpi,jpj, zworka ) 445 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 406 407 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zaTsfn 408 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 409 410 REAL(wp) :: ztrans ! ice/snow transferred 411 !! 446 412 447 413 ! 448 414 ! 1) Define a variable equal to a_i*T_su 449 415 ! 450 451 416 DO jl = 1, jpl 452 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 417 DO jj = 1, jpj 418 DO ji = 1, jpi 419 zaTsfn(ji,jj,jl) = a_i(ji,jj,jl) * t_su(ji,jj,jl) 420 END DO 421 END DO 453 422 END DO 454 423 455 424 ! 456 425 ! 2) Transfer volume and energy between categories 457 426 ! 458 459 427 DO jl = 1, jpl  1 460 nbrem = 0461 428 DO jj = 1, jpj 462 429 DO ji = 1, jpi 463 IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny 464 nbrem = nbrem + 1 465 nind_i(nbrem) = ji 466 nind_j(nbrem) = jj 467 ENDIF 468 END DO 469 END DO 470 471 DO ji = 1, nbrem 472 ii = nind_i(ji) 473 ij = nind_j(ji) 474 475 jl1 = zdonor(ii,ij,jl) 476 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ii,ij,jl1)  epsi10 ) ) 477 zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch 478 IF( jl1 == jl) THEN ; jl2 = jl1+1 479 ELSE ; jl2 = jl 480 ENDIF 481 482 ! Ice areas 483 a_i(ii,ij,jl1) = a_i(ii,ij,jl1)  zdaice(ii,ij,jl) 484 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 485 486 ! Ice volumes 487 v_i(ii,ij,jl1) = v_i(ii,ij,jl1)  zdvice(ii,ij,jl) 488 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 489 490 ! Snow volumes 491 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 492 v_s(ii,ij,jl1) = v_s(ii,ij,jl1)  zdvsnow 493 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 494 495 ! Snow heat content 496 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 497 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1)  zdesnow 498 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow 499 500 ! Ice age 501 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 502 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1)  zdo_aice 503 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice 504 505 ! Ice salinity 506 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 507 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1)  zdsm_vice 508 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice 509 510 ! Surface temperature 511 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 512 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1)  zdaTsf 513 zaTsfn(ii,ij,jl2) = zaTsfn(ii,ij,jl2) + zdaTsf 514 515 ! MV MP 2016 516 IF ( nn_pnd_scheme > 0 ) THEN 517 ! Pond fraction 518 zdapnd = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 519 a_ip(ii,ij,jl1) = a_ip(ii,ij,jl1)  zdapnd 520 a_ip(ii,ij,jl2) = a_ip(ii,ij,jl2) + zdapnd 521 522 ! Pond volume 523 zdvpnd = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 524 v_ip(ii,ij,jl1) = v_ip(ii,ij,jl1)  zdvpnd 525 v_ip(ii,ij,jl2) = v_ip(ii,ij,jl2) + zdvpnd 526 527 ENDIF 528 ! END MV MP 2016 529 530 END DO 531 430 431 jl1 = zdonor(ji,jj,jl) 432 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i(ji,jj,jl1)  epsi10 ) ) 433 zworka(ji,jj) = zdvice(ji,jj,jl) / MAX( v_i(ji,jj,jl1), epsi10 ) * rswitch 434 435 IF( jl1 == jl) THEN ; jl2 = jl1+1 436 ELSE ; jl2 = jl 437 ENDIF 438 439 ! Ice areas 440 a_i(ji,jj,jl1) = a_i(ji,jj,jl1)  zdaice(ji,jj,jl) 441 a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + zdaice(ji,jj,jl) 442 443 ! Ice volumes 444 v_i(ji,jj,jl1) = v_i(ji,jj,jl1)  zdvice(ji,jj,jl) 445 v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + zdvice(ji,jj,jl) 446 447 ! Snow volumes 448 ztrans = v_s(ji,jj,jl1) * zworka(ji,jj) 449 v_s(ji,jj,jl1) = v_s(ji,jj,jl1)  ztrans 450 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + ztrans 451 452 ! Snow heat content 453 ztrans = e_s(ji,jj,1,jl1) * zworka(ji,jj) 454 e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1)  ztrans 455 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + ztrans 456 457 ! Ice age 458 ztrans = oa_i(ji,jj,jl1) * zdaice(ji,jj,jl) 459 oa_i(ji,jj,jl1) = oa_i(ji,jj,jl1)  ztrans 460 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + ztrans 461 462 ! Ice salinity 463 ztrans = smv_i(ji,jj,jl1) * zworka(ji,jj) 464 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1)  ztrans 465 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + ztrans 466 467 ! Surface temperature 468 ztrans = t_su(ji,jj,jl1) * zdaice(ji,jj,jl) 469 zaTsfn(ji,jj,jl1) = zaTsfn(ji,jj,jl1)  ztrans 470 zaTsfn(ji,jj,jl2) = zaTsfn(ji,jj,jl2) + ztrans 471 472 ! MV MP 2016 473 IF ( nn_pnd_scheme > 0 ) THEN 474 ! Pond fraction 475 ztrans = a_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 476 a_ip(ji,jj,jl1) = a_ip(ji,jj,jl1)  ztrans 477 a_ip(ji,jj,jl2) = a_ip(ji,jj,jl2) + ztrans 478 479 ! Pond volume 480 ztrans = v_ip(ji,jj,jl1) * zdaice(ji,jj,jl) 481 v_ip(ji,jj,jl1) = v_ip(ji,jj,jl1)  ztrans 482 v_ip(ji,jj,jl2) = v_ip(ji,jj,jl2) + ztrans 483 ENDIF 484 ! END MV MP 2016 485 486 END DO 487 END DO 488 532 489 ! Ice heat content 533 490 DO jk = 1, nlay_i 534 DO ji = 1, nbrem 535 ii = nind_i(ji) 536 ij = nind_j(ji) 537 538 jl1 = zdonor(ii,ij,jl) 539 IF (jl1 == jl) THEN 540 jl2 = jl+1 541 ELSE ! n1 = n+1 542 jl2 = jl 543 ENDIF 544 545 zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 546 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1)  zdeice 547 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + zdeice 548 END DO 549 END DO 550 551 END DO ! boundaries, 1 to ncat1 552 553 ! 491 DO jj = 1, jpj 492 DO ji = 1, jpi 493 494 jl1 = zdonor(ji,jj,jl) 495 IF(jl1 == jl) THEN ; jl2 = jl+1 496 ELSE ; jl2 = jl 497 ENDIF 498 499 ztrans = e_i(ji,jj,jk,jl1) * zworka(ji,jj) 500 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1)  ztrans 501 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + ztrans 502 END DO 503 END DO 504 END DO 505 506 END DO ! boundaries, 1 to jpl1 507 554 508 ! Update ice thickness and temperature 555 !556 557 509 DO jl = 1, jpl 558 510 DO jj = 1, jpj 559 DO ji = 1, jpi 511 DO ji = 1, jpi 560 512 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 561 513 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) … … 567 519 END DO 568 520 END DO 569 END DO 570 ! 571 CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn ) 572 CALL wrk_dealloc( jpi,jpj, zworka ) 573 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 521 END DO 574 522 ! 575 523 END SUBROUTINE lim_itd_shiftice … … 584 532 !! ** Method : 585 533 !! 586 INTEGER :: ji, jj, jl ! dummy loop indices534 INTEGER :: ji, jj, jl ! dummy loop indices 587 535 INTEGER :: zshiftflag ! = .true. if ice must be shifted 588 536 589 INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index590 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! ice area and volume transferred537 INTEGER , DIMENSION(jpi,jpj,jpl) :: zdonor ! donor category index 538 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zdaice, zdvice ! ice area and volume transferred 591 539 !! 592 540 593 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger594 CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice )595 !596 541 ! 597 542 ! 1) Compute ice thickness. … … 635 580 zshiftflag = 1 636 581 zdonor(ji,jj,jl) = jl 637 ! begin TECLIM change638 !zdaice(ji,jj,jl) = a_i(ji,jj,jl) * 0.5_wp639 !zdvice(ji,jj,jl) = v_i(ji,jj,jl)zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl1)) * 0.5_wp640 ! end TECLIM change641 582 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 642 583 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl)  hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) … … 645 586 END DO 646 587 END DO 647 IF(lk_mpp) CALL mpp_max( zshiftflag ) 588 IF(lk_mpp) CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 648 589 649 590 IF( zshiftflag == 1 ) THEN ! Shift ice between categories … … 680 621 END DO 681 622 682 IF(lk_mpp) CALL mpp_max( zshiftflag ) 623 IF(lk_mpp) CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 683 624 684 625 IF( zshiftflag == 1 ) THEN ! Shift ice between categories … … 692 633 END DO 693 634 ! 694 CALL wrk_dealloc( jpi,jpj,jpl, zdonor )695 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice )696 697 635 END SUBROUTINE lim_itd_th_reb 698 636 … … 704 642 SUBROUTINE lim_itd_th_rem 705 643 END SUBROUTINE lim_itd_th_rem 706 SUBROUTINE lim_itd_ fitline707 END SUBROUTINE lim_itd_ fitline644 SUBROUTINE lim_itd_glinear 645 END SUBROUTINE lim_itd_glinear 708 646 SUBROUTINE lim_itd_shiftice 709 647 END SUBROUTINE lim_itd_shiftice
Note: See TracChangeset
for help on using the changeset viewer.