Changeset 8369 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM
- Timestamp:
- 2017-07-25T16:38:38+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 8 edited
-
limitd_th.F90 (modified) (14 diffs)
-
limmp.F90 (modified) (3 diffs)
-
limtab.F90 (modified) (3 diffs)
-
limthd.F90 (modified) (3 diffs)
-
limthd_da.F90 (modified) (1 diff)
-
limthd_dh.F90 (modified) (2 diffs)
-
limthd_dif.F90 (modified) (1 diff)
-
thd_ice.F90 (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r8361 r8369 25 25 USE limvar ! LIM-3 variables 26 26 USE limcons ! conservation tests 27 USE limtab 27 28 ! 28 29 USE prtctl ! Print control … … 35 36 PRIVATE 36 37 37 PUBLIC lim_itd_th_rem 38 PUBLIC lim_itd_th_reb 38 PUBLIC lim_itd_th_rem ! called in limthd 39 PUBLIC lim_itd_th_reb ! called in limupdate 39 40 40 41 !!---------------------------------------------------------------------- … … 58 59 INTEGER , INTENT (in) :: kt ! Ocean time step 59 60 ! 60 INTEGER :: ji, jj, jl ! dummy loop index 61 INTEGER :: ii, ij ! 2D corresponding indices to ji 62 INTEGER :: nd ! local integer 61 INTEGER :: ji, jj, jl, jcat ! dummy loop index 62 INTEGER :: nidx2 ! local integer 63 63 REAL(wp) :: zx1, zwk1, zdh0, zetamin, zdamax ! local scalars 64 64 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 65 65 REAL(wp) :: zx3 66 67 INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index 68 69 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdhice ! ice thickness increment 70 REAL(wp), POINTER, DIMENSION(:,:,:) :: g0 ! coefficients for fitting the line of the ITD 71 REAL(wp), POINTER, DIMENSION(:,:,:) :: g1 ! coefficients for fitting the line of the ITD 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: hL ! left boundary for the ITD for each thickness 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume 75 INTEGER , DIMENSION(jpij) :: nind_i, nind_j ! compressed indices for i/j directions 76 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 77 REAL(wp), POINTER, DIMENSION(:,:) :: zhb0, zhb1 ! category boundaries for thinnes categories 78 INTEGER , POINTER, DIMENSION(:,:) :: zremap_flag ! compute remapping or not ???? 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhbnew ! new boundaries of ice categories 80 !!------------------------------------------------------------------ 81 82 CALL wrk_alloc( jpi,jpj, zremap_flag ) 83 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 84 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR ) 85 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 86 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 87 CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 66 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 67 68 INTEGER , DIMENSION(jpij) :: idxice2 ! compute remapping or not 69 INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index 70 REAL(wp), DIMENSION(jpij,jpl) :: zdhice ! ice thickness increment 71 REAL(wp), DIMENSION(jpij,jpl) :: g0, g1 ! coefficients for fitting the line of the ITD 72 REAL(wp), DIMENSION(jpij,jpl) :: hL, hR ! left and right boundary for the ITD for each thickness 73 REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! local increment of ice area and volume 74 REAL(wp), DIMENSION(jpij) :: zhb0, zhb1 ! category boundaries for thinnes categories 75 REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories 76 !!------------------------------------------------------------------ 88 77 89 78 IF( kt == nit000 .AND. lwp) THEN … … 94 83 95 84 !----------------------------------------------------------------------------------------------- 96 ! 3) Identify grid cells with ice85 ! 1) Identify grid cells with ice 97 86 !----------------------------------------------------------------------------------------------- 98 87 nidx = 0 ; idxice(:) = 0 … … 101 90 IF ( at_i(ji,jj) > epsi10 ) THEN 102 91 nidx = nidx + 1 103 nind_i(nidx) = ji104 nind_j(nidx) = jj105 92 idxice( nidx ) = (jj - 1) * jpi + ji 106 zremap_flag(ji,jj) = 1107 ELSE108 zremap_flag(ji,jj) = 0109 93 ENDIF 110 94 END DO 111 95 END DO 112 96 113 97 !----------------------------------------------------------------------------------------------- 114 ! 4) Compute new category boundaries: 1:jpl-198 ! 2) Compute new category boundaries 115 99 !----------------------------------------------------------------------------------------------- 116 zdhice(:,:,:) = 0._wp117 zhbnew(:,:,:) = 0._wp118 zhb0(:,:) = hi_max(0)119 zhb1(:,:) = hi_max(1)120 121 100 IF( nidx > 0 ) THEN 122 101 123 ! Compute thickness change in each ice category 102 zdhice(:,:) = 0._wp 103 zhbnew(:,:) = 0._wp 104 105 CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i ) 106 CALL tab_3d_2d( nidx, idxice(1:nidx), ht_ib_2d(1:nidx,1:jpl), ht_i_b ) 107 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 108 CALL tab_3d_2d( nidx, idxice(1:nidx), a_ib_2d (1:nidx,1:jpl), a_i_b ) 109 124 110 DO jl = 1, jpl 111 ! Compute thickness change in each ice category 125 112 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) 113 zdhice(ji,jl) = ht_i_2d(ji,jl) - ht_ib_2d(ji,jl) 129 114 END DO 130 115 END DO 131 116 117 ! --- New boundaries for category 1:jpl-1 --- ! 132 118 DO jl = 1, jpl - 1 119 133 120 DO ji = 1, nidx 134 ii = nind_i(ji)135 ij = nind_j(ji)136 121 ! 137 ! --- New boundary categoriesHn* = Hn + Fn*dt --- !122 ! --- New boundary: Hn* = Hn + Fn*dt --- ! 138 123 ! Fn*dt = ( fn + (fn+1 - fn)/(hn+1 - hn) * (Hn - hn) ) * dt = zdhice + zslope * (Hmax - ht_i_b) 139 124 ! 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) /= 0141 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*dt144 zhbnew( ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl)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*dt146 zhbnew( ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1)147 ELSE ! a_i_b(jl+1) & a_i_b(jl) = 0148 zhbnew( ii,ij,jl) = hi_max(jl)125 IF ( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl+1) & a(jl) /= 0 126 zslope = ( zdhice(ji,jl+1) - zdhice(ji,jl) ) / ( ht_ib_2d(ji,jl+1) - ht_ib_2d(ji,jl) ) 127 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) + zslope * ( hi_max(jl) - ht_ib_2d(ji,jl) ) 128 ELSEIF( a_ib_2d(ji,jl) > epsi10 .AND. a_ib_2d(ji,jl+1) <= epsi10 ) THEN ! a(jl+1)=0 => Hn* = Hn + fn*dt 129 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl) 130 ELSEIF( a_ib_2d(ji,jl) <= epsi10 .AND. a_ib_2d(ji,jl+1) > epsi10 ) THEN ! a(jl)=0 => Hn* = Hn + fn+1*dt 131 zhbnew(ji,jl) = hi_max(jl) + zdhice(ji,jl+1) 132 ELSE ! a(jl+1) & a(jl) = 0 133 zhbnew(ji,jl) = hi_max(jl) 149 134 ENDIF 150 135 … … 153 138 ! Note: hn(t+1) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 154 139 ! in lim_itd_glinear in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 155 IF( a_i (ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) zremap_flag(ii,ij) = 0156 IF( a_i (ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) zremap_flag(ii,ij) = 0140 IF( a_i_2d(ji,jl ) > epsi10 .AND. ht_i_2d(ji,jl ) > ( zhbnew(ji,jl) - epsi10 ) ) idxice(ji) = 0 141 IF( a_i_2d(ji,jl+1) > epsi10 .AND. ht_i_2d(ji,jl+1) < ( zhbnew(ji,jl) + epsi10 ) ) idxice(ji) = 0 157 142 158 143 ! 2) Hn-1 < Hn* < Hn+1 159 IF( zhbnew( ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0160 IF( zhbnew( ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0144 IF( zhbnew(ji,jl) < hi_max(jl-1) ) idxice(ji) = 0 145 IF( zhbnew(ji,jl) > hi_max(jl+1) ) idxice(ji) = 0 161 146 162 147 END DO 163 148 END DO 164 149 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(jpl-1), 3._wp * ht_i(ii,ij,jpl) - 2._wp * zhbnew(ii,ij,jpl-1) ) 150 ! --- New boundaries for category jpl --- ! 151 DO ji = 1, nidx 152 IF( a_i_2d(ji,jpl) > epsi10 ) THEN 153 zhbnew(ji,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i_2d(ji,jpl) - 2._wp * zhbnew(ji,jpl-1) ) 171 154 ELSE 172 zhbnew( ii,ij,jpl) = hi_max(jpl)173 ENDIF 174 175 ! 1 condition for remapping for the 1st category155 zhbnew(ji,jpl) = hi_max(jpl) 156 ENDIF 157 158 ! --- 1 additional condition for remapping (1st category) --- ! 176 159 ! H0+epsi < h1(t) < H1-epsi 177 160 ! h1(t) must not be too close to either HR or HL otherwise a division by nearly 0 is possible 178 161 ! in lim_itd_glinear in the case (HR-HL) = 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 162 IF( ht_ib_2d(ji,1) < ( hi_max(0) + epsi10 ) ) idxice(ji) = 0 163 IF( ht_ib_2d(ji,1) > ( hi_max(1) - epsi10 ) ) idxice(ji) = 0 164 END DO 165 166 !----------------------------------------------------------------------------------------------- 167 ! 3) Identify cells where remapping 168 !----------------------------------------------------------------------------------------------- 169 nidx2 = 0 ; idxice2(:) = 0 170 DO ji = 1, nidx 171 IF( idxice(ji) /= 0 ) THEN 172 nidx2 = nidx2 + 1 173 idxice2(nidx2) = idxice(ji) 174 zhbnew(nidx2,:) = zhbnew(ji,:) ! adjust zhbnew to new indices 175 ENDIF 176 END DO 177 idxice(:) = idxice2(:) 178 nidx = nidx2 182 179 183 180 ENDIF 184 181 182 185 183 !----------------------------------------------------------------------------------------------- 186 ! 5) Identify cells where ITD is to be remapped 187 !----------------------------------------------------------------------------------------------- 188 nidx = 0 189 DO jj = 1, jpj 190 DO ji = 1, jpi 191 IF( zremap_flag(ji,jj) == 1 ) THEN 192 nidx = nidx + 1 193 nind_i(nidx) = ji 194 nind_j(nidx) = jj 195 ENDIF 196 END DO 197 END DO 198 199 !----------------------------------------------------------------------------------------------- 200 ! 7) Compute g(h) 184 ! 4) Compute g(h) 201 185 !----------------------------------------------------------------------------------------------- 202 186 IF( nidx > 0 ) THEN 203 187 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) 188 zhb0(:) = hi_max(0) ; zhb1(:) = hi_max(1) 189 g0(:,:) = 0._wp ; g1(:,:) = 0._wp 190 hL(:,:) = 0._wp ; hR(:,:) = 0._wp 191 192 DO jl = 1, jpl 193 194 CALL tab_2d_1d( nidx, idxice(1:nidx), ht_ib_1d(1:nidx), ht_i_b(:,:,jl) ) 195 CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) ) 196 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 197 CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) 198 199 IF( jl == 1 ) THEN 200 201 ! --- g(h) for category 1 --- ! 202 CALL lim_itd_glinear( zhb0(1:nidx), zhb1(1:nidx), ht_ib_1d(1:nidx), a_i_1d(1:nidx), & ! in 203 & g0(1:nidx,1), g1(1:nidx,1), hL(1:nidx,1) , hR(1:nidx,1) ) ! out 204 205 ! Area lost due to melting of thin ice 206 DO ji = 1, nidx 221 207 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) ! clem-useless ? 208 IF( a_i_1d(ji) > epsi10 ) THEN 209 210 zdh0 = ht_i_1d(ji) - ht_ib_1d(ji) 211 IF( zdh0 < 0.0 ) THEN !remove area from category 1 212 zdh0 = MIN( -zdh0, hi_max(1) ) 213 !Integrate g(1) from 0 to dh0 to estimate area melted 214 zetamax = MIN( zdh0, hR(ji,1) ) - hL(ji,1) 215 216 IF( zetamax > 0.0 ) THEN 217 zx1 = zetamax 218 zx2 = 0.5 * zetamax * zetamax 219 zda0 = g1(ji,1) * zx2 + g0(ji,1) * zx1 ! ice area removed 220 zdamax = a_i_1d(ji) * (1.0 - ht_i_1d(ji) / ht_ib_1d(ji) ) ! Constrain new thickness <= ht_i 221 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 222 ! of thin ice (zdamax > 0) 223 ! Remove area, conserving volume 224 ht_i_1d(ji) = ht_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 ) 225 a_i_1d(ji) = a_i_1d(ji) - zda0 226 v_i_1d(ji) = a_i_1d(ji) * ht_i_1d(ji) ! clem-useless ? 227 ENDIF 228 229 ELSE ! if ice accretion zdh0 > 0 230 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 231 zhbnew(ji,0) = MIN( zdh0, hi_max(1) ) 232 ENDIF 233 233 234 ENDIF 234 235 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 240 ENDIF 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(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 247 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) 236 END DO 237 238 CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) ) 239 CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 240 CALL tab_1d_2d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) 241 242 ENDIF ! jl=1 243 244 ! --- g(h) for each thickness category --- ! 245 CALL lim_itd_glinear( zhbnew(1:nidx,jl-1), zhbnew(1:nidx,jl), ht_i_1d(1:nidx), a_i_1d(1:nidx), & ! in 246 & g0(1:nidx,jl) , g1(1:nidx,jl) , hL(1:nidx,jl) , hR(1:nidx,jl) ) ! out 247 248 248 END DO 249 249 250 250 !----------------------------------------------------------------------------------------------- 251 ! 8) Compute area and volume to be shifted across each boundary (Eq. 18)251 ! 5) Compute area and volume to be shifted across each boundary (Eq. 18) 252 252 !----------------------------------------------------------------------------------------------- 253 254 253 DO jl = 1, jpl - 1 255 DO jj = 1, jpj256 DO ji = 1, jpi257 zdonor(ji,jj,jl) = 0258 zdaice(ji,jj,jl) = 0.0259 zdvice(ji,jj,jl) = 0.0260 END DO261 END DO262 254 263 255 DO ji = 1, nidx 264 ii = nind_i(ji)265 ij = nind_j(ji)266 256 267 257 ! 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+1269 zetamin = MAX( hi_max(jl) , hL(ii,ij,jl) ) - hL(ii,ij,jl)! hi_max(jl) - hL270 zetamax = MIN( zhbnew( ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)! hR - hL271 zdonor(ii,ij,jl) = jl272 ELSE 258 IF (zhbnew(ji,jl) > hi_max(jl)) THEN ! Hn* > Hn => transfer from jl to jl+1 259 zetamin = MAX( hi_max(jl) , hL(ji,jl) ) - hL(ji,jl) ! hi_max(jl) - hL 260 zetamax = MIN( zhbnew(ji,jl), hR(ji,jl) ) - hL(ji,jl) ! hR - hL 261 jdonor(ji,jl) = jl 262 ELSE ! Hn* <= Hn => transfer from jl+1 to jl 273 263 zetamin = 0.0 274 zetamax = MIN( hi_max(jl), hR( ii,ij,jl+1) ) - hL(ii,ij,jl+1)! hi_max(jl) - hL275 zdonor(ii,ij,jl) = jl + 1264 zetamax = MIN( hi_max(jl), hR(ji,jl+1) ) - hL(ji,jl+1) ! hi_max(jl) - hL 265 jdonor(ji,jl) = jl + 1 276 266 ENDIF 277 267 zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin … … 284 274 zwk2 = zwk2 * zetamax 285 275 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)*zx1288 zdvice( ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)276 jcat = jdonor(ji,jl) 277 zdaice(ji,jl) = g1(ji,jcat)*zx2 + g0(ji,jcat)*zx1 278 zdvice(ji,jl) = g1(ji,jcat)*zx3 + g0(ji,jcat)*zx2 + zdaice(ji,jl)*hL(ji,jcat) 289 279 290 280 END DO 291 281 END DO 292 282 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 303 ii = nind_i(ji) 304 ij = nind_j(ji) 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 283 !---------------------------------------------------------------------------------------------- 284 ! 6) Shift ice between categories 285 !---------------------------------------------------------------------------------------------- 286 CALL lim_itd_shiftice ( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) 287 288 !---------------------------------------------------------------------------------------------- 289 ! 7) Make sure ht_i >= minimum ice thickness hi_min 290 !---------------------------------------------------------------------------------------------- 291 CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1) ) 292 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) ) 293 CALL tab_2d_1d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) ) 294 295 DO ji = 1, nidx 296 IF ( a_i_1d(ji) > epsi10 .AND. ht_i_1d(ji) < rn_himin ) THEN 297 a_i_1d (ji) = a_i_1d(ji) * ht_i_1d(ji) / rn_himin 307 298 ! MV MP 2016 308 299 IF ( nn_pnd_scheme > 0 ) THEN 309 a_ip (ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin300 a_ip_1d(ji) = a_ip_1d(ji) * ht_i_1d(ji) / rn_himin 310 301 ENDIF 311 302 ! END MV MP 2016 312 ht_i(ii,ij,1) = rn_himin 313 ENDIF 314 END DO 315 303 ht_i_1d(ji) = rn_himin 304 ENDIF 305 END DO 306 307 CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,1) ) 308 CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,1) ) 309 CALL tab_1d_2d( nidx, idxice(1:nidx), a_ip_1d (1:nidx), a_ip(:,:,1) ) 310 316 311 ENDIF 317 312 318 CALL wrk_dealloc( jpi,jpj, zremap_flag )319 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor )320 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR )321 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice )322 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 )323 CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 )324 325 313 END SUBROUTINE lim_itd_th_rem 326 314 327 315 328 SUBROUTINE lim_itd_glinear( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag)316 SUBROUTINE lim_itd_glinear( HbL, Hbr, phice, paice, pg0, pg1, phL, phR ) 329 317 !!------------------------------------------------------------------ 330 318 !! *** ROUTINE lim_itd_glinear *** … … 336 324 !! 337 325 !!------------------------------------------------------------------ 338 INTEGER , INTENT(in ) :: num_cat ! category index 339 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: HbL, HbR ! left and right category boundaries 340 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: hice ! ice thickness 341 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: g0, g1 ! coefficients in linear equation for g(eta) 342 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hL ! min value of range over which g(h) > 0 343 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: hR ! max value of range over which g(h) > 0 344 INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: zremap_flag ! 326 REAL(wp), DIMENSION(:), INTENT(in ) :: HbL, HbR ! left and right category boundaries 327 REAL(wp), DIMENSION(:), INTENT(in ) :: phice, paice ! ice thickness and concentration 328 REAL(wp), DIMENSION(:), INTENT(inout) :: pg0, pg1 ! coefficients in linear equation for g(eta) 329 REAL(wp), DIMENSION(:), INTENT(inout) :: phL, phR ! min and max value of range over which g(h) > 0 345 330 ! 346 INTEGER :: ji ,jj! horizontal indices331 INTEGER :: ji ! horizontal indices 347 332 REAL(wp) :: zh13 ! HbL + 1/3 * (HbR - HbL) 348 333 REAL(wp) :: zh23 ! HbL + 2/3 * (HbR - HbL) … … 351 336 !!------------------------------------------------------------------ 352 337 ! 353 DO jj = 1, jpj 354 DO ji = 1, jpi 338 DO ji = 1, nidx 355 339 ! 356 IF( zremap_flag(ji,jj) == 1 .AND. a_i(ji,jj,num_cat) > epsi10 & 357 & .AND. hice(ji,jj) > 0._wp ) THEN 340 IF( paice(ji) > epsi10 .AND. phice(ji) > 0._wp ) THEN 358 341 359 342 ! Initialize hL and hR 360 hL(ji,jj) = HbL(ji,jj)361 hR(ji,jj) = HbR(ji,jj)343 phL(ji) = HbL(ji) 344 phR(ji) = HbR(ji) 362 345 363 346 ! Change hL or hR if hice falls outside central third of range, 364 347 ! so that hice is in the central third of the range [HL HR] 365 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )366 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )367 368 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left369 ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right348 zh13 = 1.0 / 3.0 * ( 2.0 * phL(ji) + phR(ji) ) 349 zh23 = 1.0 / 3.0 * ( phL(ji) + 2.0 * phR(ji) ) 350 351 IF ( phice(ji) < zh13 ) THEN ; phR(ji) = 3._wp * phice(ji) - 2._wp * phL(ji) ! move HR to the left 352 ELSEIF( phice(ji) > zh23 ) THEN ; phL(ji) = 3._wp * phice(ji) - 2._wp * phR(ji) ! move HL to the right 370 353 ENDIF 371 354 372 355 ! Compute coefficients of g(eta) = g0 + g1*eta 373 zdhr = 1._wp / ( hR(ji,jj) - hL(ji,jj))374 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr375 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr376 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) ! Eq. 14377 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14356 zdhr = 1._wp / (phR(ji) - phL(ji)) 357 zwk1 = 6._wp * paice(ji) * zdhr 358 zwk2 = ( phice(ji) - phL(ji) ) * zdhr 359 pg0(ji) = zwk1 * ( 2._wp / 3._wp - zwk2 ) ! Eq. 14 360 pg1(ji) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 378 361 ! 379 362 ELSE ! remap_flag = .false. or a_i < epsi10 380 hL(ji,jj) = 0._wp381 hR(ji,jj) = 0._wp382 g0(ji,jj) = 0._wp383 g1(ji,jj) = 0._wp363 phL(ji) = 0._wp 364 phR(ji) = 0._wp 365 pg0(ji) = 0._wp 366 pg1(ji) = 0._wp 384 367 ENDIF 385 368 ! 386 369 END DO 387 END DO388 370 ! 389 371 END SUBROUTINE lim_itd_glinear 390 372 391 373 392 SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice )374 SUBROUTINE lim_itd_shiftice( kdonor, pdaice, pdvice ) 393 375 !!------------------------------------------------------------------ 394 376 !! *** ROUTINE lim_itd_shiftice *** … … 399 381 !! ** Method : 400 382 !!------------------------------------------------------------------ 401 INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in ) :: zdonor ! donor category index 402 REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdaice ! ice area transferred across boundary 403 REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdvice ! ice volume transferred across boundary 404 405 INTEGER :: ji, jj, jl, jl2, jl1, jk ! dummy loop indices 406 407 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zaTsfn 408 REAL(wp), DIMENSION(jpi,jpj) :: zworka ! temporary array used here 383 INTEGER , DIMENSION(:,:), INTENT(in) :: kdonor ! donor category index 384 REAL(wp), DIMENSION(:,:), INTENT(in) :: pdaice ! ice area transferred across boundary 385 REAL(wp), DIMENSION(:,:), INTENT(in) :: pdvice ! ice volume transferred across boundary 386 387 INTEGER :: ii,ij, ji, jj, jl, jl2, jl1, jk ! dummy loop indices 388 389 REAL(wp), DIMENSION(jpij,jpl) :: zaTsfn 390 REAL(wp), DIMENSION(jpij) :: zworka ! temporary array used here 391 REAL(wp), DIMENSION(jpij) :: zworkv ! temporary array used here 409 392 410 393 REAL(wp) :: ztrans ! ice/snow transferred 411 394 !!------------------------------------------------------------------ 395 396 CALL tab_3d_2d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i ) 397 CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 398 CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i ) 399 CALL tab_3d_2d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s ) 400 CALL tab_3d_2d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i ) 401 CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i ) 402 CALL tab_3d_2d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip ) 403 CALL tab_3d_2d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip ) 404 CALL tab_3d_2d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su ) 412 405 413 406 !---------------------------------------------------------------------------------------------- … … 415 408 !---------------------------------------------------------------------------------------------- 416 409 DO jl = 1, jpl 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 410 DO ji = 1, nidx 411 zaTsfn(ji,jl) = a_i_2d(ji,jl) * t_su_2d(ji,jl) 421 412 END DO 422 413 END DO … … 426 417 !------------------------------------------------------------------------------- 427 418 DO jl = 1, jpl - 1 428 DO jj = 1, jpj 429 DO ji = 1, jpi 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 419 DO ji = 1, nidx 420 421 jl1 = kdonor(ji,jl) 422 IF ( jl1 == jl ) THEN ; jl2 = jl1+1 423 ELSE ; jl2 = jl 424 ENDIF 425 426 rswitch = MAX( 0._wp , SIGN( 1._wp , v_i_2d(ji,jl1) - epsi10 ) ) 427 zworkv(ji) = pdvice(ji,jl) / MAX( v_i_2d(ji,jl1), epsi10 ) * rswitch 428 429 rswitch = MAX( 0._wp , SIGN( 1._wp , a_i_2d(ji,jl1) - epsi10 ) ) 430 zworka(ji) = pdaice(ji,jl) / MAX( a_i_2d(ji,jl1), epsi10 ) * rswitch 431 432 ! Ice areas 433 a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl) 434 a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl) 435 436 ! Ice volumes 437 v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl) 438 v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl) 439 440 ! Snow volumes 441 ztrans = v_s_2d(ji,jl1) * zworkv(ji) 442 v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans 443 v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 444 445 ! Ice age 446 ztrans = oa_i_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 447 oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - ztrans 448 oa_i_2d(ji,jl2) = oa_i_2d(ji,jl2) + ztrans 449 450 ! Ice salinity 451 ztrans = smv_i_2d(ji,jl1) * zworkv(ji) 452 smv_i_2d(ji,jl1) = smv_i_2d(ji,jl1) - ztrans 453 smv_i_2d(ji,jl2) = smv_i_2d(ji,jl2) + ztrans 454 455 ! Surface temperature 456 ztrans = t_su_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 457 zaTsfn(ji,jl1) = zaTsfn(ji,jl1) - ztrans 458 zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans 459 460 ! MV MP 2016 461 IF ( nn_pnd_scheme > 0 ) THEN 462 ! Pond fraction 463 ztrans = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 464 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 465 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 466 467 ! Pond volume (also proportional to da/a) 468 ztrans = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 469 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 470 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 471 ENDIF 472 ! END MV MP 2016 473 474 END DO 475 476 ! Snow heat content 477 DO jk = 1, nlay_s 478 479 DO ji = 1, nidx 480 ii = MOD( idxice(ji) - 1, jpi ) + 1 481 ij = ( idxice(ji) - 1 ) / jpi + 1 482 483 jl1 = kdonor(ji,jl) 484 IF(jl1 == jl) THEN ; jl2 = jl+1 485 ELSE ; jl2 = jl 486 ENDIF 487 488 ztrans = e_s(ii,ij,jk,jl1) * zworkv(ji) 489 e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans 490 e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 486 491 END DO 487 492 END DO 493 488 494 489 495 ! Ice heat content 490 496 DO jk = 1, nlay_i 491 DO j j = 1, jpj492 DO ji = 1, jpi493 494 jl1 = zdonor(ji,jj,jl)495 IF(jl1 == jl) THEN ; jl2 = jl+1496 ELSE ; jl2 = jl497 ENDIF498 499 ztrans = e_i(ji,jj,jk,jl1) * zworka(ji,jj)500 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - ztrans501 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + ztrans502 END DO497 DO ji = 1, nidx 498 ii = MOD( idxice(ji) - 1, jpi ) + 1 499 ij = ( idxice(ji) - 1 ) / jpi + 1 500 501 jl1 = kdonor(ji,jl) 502 IF(jl1 == jl) THEN ; jl2 = jl+1 503 ELSE ; jl2 = jl 504 ENDIF 505 506 ztrans = e_i(ii,ij,jk,jl1) * zworkv(ji) 507 e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans 508 e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 503 509 END DO 504 510 END DO … … 506 512 END DO ! boundaries, 1 to jpl-1 507 513 508 ! Update ice thickness and temperature 514 !------------------------------------------------------------------------------- 515 ! 3) Update ice thickness and temperature 516 !------------------------------------------------------------------------------- 509 517 DO jl = 1, jpl 510 DO jj = 1, jpj 511 DO ji = 1, jpi 512 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 513 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl) 514 t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 515 ELSE 516 ht_i(ji,jj,jl) = 0._wp 517 t_su(ji,jj,jl) = rt0 518 ENDIF 519 END DO 520 END DO 521 END DO 518 DO ji = 1, nidx 519 IF ( a_i_2d(ji,jl) > epsi10 ) THEN 520 ht_i_2d(ji,jl) = v_i_2d (ji,jl) / a_i_2d(ji,jl) 521 t_su_2d(ji,jl) = zaTsfn(ji,jl) / a_i_2d(ji,jl) 522 ELSE 523 ht_i_2d(ji,jl) = 0._wp 524 t_su_2d(ji,jl) = rt0 525 ENDIF 526 END DO 527 END DO 528 529 CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i ) 530 CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d (1:nidx,1:jpl), a_i ) 531 CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d (1:nidx,1:jpl), v_i ) 532 CALL tab_2d_3d( nidx, idxice(1:nidx), v_s_2d (1:nidx,1:jpl), v_s ) 533 CALL tab_2d_3d( nidx, idxice(1:nidx), oa_i_2d (1:nidx,1:jpl), oa_i ) 534 CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i ) 535 CALL tab_2d_3d( nidx, idxice(1:nidx), a_ip_2d (1:nidx,1:jpl), a_ip ) 536 CALL tab_2d_3d( nidx, idxice(1:nidx), v_ip_2d (1:nidx,1:jpl), v_ip ) 537 CALL tab_2d_3d( nidx, idxice(1:nidx), t_su_2d (1:nidx,1:jpl), t_su ) 538 522 539 ! 523 540 END SUBROUTINE lim_itd_shiftice … … 530 547 !! ** Purpose : rebin - rebins thicknesses into defined categories 531 548 !! 532 !! ** Method : 549 !! ** Method : If a category thickness is out of bounds, shift part (for down to top) 550 !! or entire (for top to down) area, volume, and energy 551 !! to the neighboring category 533 552 !!------------------------------------------------------------------ 534 553 INTEGER :: ji, jj, jl ! dummy loop indices 535 INTEGER :: zshiftflag ! = .true. if ice must be shifted 536 537 INTEGER , DIMENSION(jpi,jpj,jpl) :: zdonor ! donor category index 538 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zdaice, zdvice ! ice area and volume transferred 554 ! 555 INTEGER , DIMENSION(jpij,jpl-1) :: jdonor ! donor category index 556 REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred 539 557 !!------------------------------------------------------------------ 540 558 541 !------------------------------------------------------------------------------ 542 ! 1) Compute ice thickness. 543 !------------------------------------------------------------------------------ 544 DO jl = 1, jpl 559 DO jl = 1, jpl-1 560 jdonor(:,jl) = 0 561 zdaice(:,jl) = 0._wp 562 zdvice(:,jl) = 0._wp 563 END DO 564 565 !--------------------------------------- 566 ! identify thicknesses that are too big 567 !--------------------------------------- 568 DO jl = 1, jpl - 1 569 570 nidx = 0 ; idxice(:) = 0 545 571 DO jj = 1, jpj 546 DO ji = 1, jpi 547 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) ) 548 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * rswitch 549 END DO 550 END DO 551 END DO 552 553 !------------------------------------------------------------------------------ 554 ! 2) If a category thickness is not in bounds, shift the 555 ! entire area, volume, and energy to the neighboring category 556 !------------------------------------------------------------------------------ 557 !------------------------- 558 ! Initialize shift arrays 559 !------------------------- 560 DO jl = 1, jpl 561 zdonor(:,:,jl) = 0 562 zdaice(:,:,jl) = 0._wp 563 zdvice(:,:,jl) = 0._wp 564 END DO 565 566 !------------------------- 567 ! Move thin categories up 568 !------------------------- 569 570 DO jl = 1, jpl - 1 ! loop over category boundaries 571 572 !--------------------------------------- 573 ! identify thicknesses that are too big 574 !--------------------------------------- 575 zshiftflag = 0 576 577 DO jj = 1, jpj 578 DO ji = 1, jpi 579 IF( a_i(ji,jj,jl) > epsi10 .AND. ht_i(ji,jj,jl) > hi_max(jl) ) THEN 580 zshiftflag = 1 581 zdonor(ji,jj,jl) = jl 582 ! clem: how much of a_i you send in cat sup is somewhat arbitrary 583 zdaice(ji,jj,jl) = a_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - hi_max(jl) + epsi20 ) / ht_i(ji,jj,jl) 584 zdvice(ji,jj,jl) = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi20 ) 585 ENDIF 586 END DO 587 END DO 588 IF(lk_mpp) CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 589 590 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 591 CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 572 DO ji = 1, jpi 573 IF( a_i(ji,jj,jl) > epsi10 .AND. v_i(ji,jj,jl) > (a_i(ji,jj,jl) * hi_max(jl)) ) THEN 574 nidx = nidx + 1 575 idxice( nidx ) = (jj - 1) * jpi + ji 576 ENDIF 577 ENDDO 578 ENDDO 579 580 !!clem CALL tab_2d_1d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) ) 581 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl) ) 582 CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl) ) 583 584 DO ji = 1, nidx 585 jdonor(ji,jl) = jl 586 ! how much of a_i you send in cat sup is somewhat arbitrary 587 !!clem: these do not work properly after a restart (I do not know why) 588 !! zdaice(ji,jl) = a_i_1d(ji) * ( ht_i_1d(ji) - hi_max(jl) + epsi10 ) / ht_i_1d(ji) 589 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 590 !!clem: these do not work properly after a restart (I do not know why) 591 !! zdaice(ji,jl) = a_i_1d(ji) 592 !! zdvice(ji,jl) = v_i_1d(ji) 593 !!clem: these are from UCL and work ok 594 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 595 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 596 597 END DO 598 599 IF( nidx > 0 ) THEN 600 CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl=>jl+1 592 601 ! Reset shift parameters 593 zdonor(:,:,jl) = 0594 zdaice( :,:,jl) = 0._wp595 zdvice( :,:,jl) = 0._wp602 jdonor(1:nidx,jl) = 0 603 zdaice(1:nidx,jl) = 0._wp 604 zdvice(1:nidx,jl) = 0._wp 596 605 ENDIF 597 606 ! 598 607 END DO 599 608 600 !---------------------------- 601 ! Move thick categories down 602 !---------------------------- 603 604 DO jl = jpl - 1, 1, -1 ! loop over category boundaries 605 606 !----------------------------------------- 607 ! Identify thicknesses that are too small 608 !----------------------------------------- 609 zshiftflag = 0 610 609 !----------------------------------------- 610 ! Identify thicknesses that are too small 611 !----------------------------------------- 612 DO jl = jpl - 1, 1, -1 613 614 nidx = 0 ; idxice(:) = 0 611 615 DO jj = 1, jpj 612 616 DO ji = 1, jpi 613 IF( a_i(ji,jj,jl+1) > epsi10 .AND. ht_i(ji,jj,jl+1) <= hi_max(jl) ) THEN 614 ! 615 zshiftflag = 1 616 zdonor(ji,jj,jl) = jl + 1 617 zdaice(ji,jj,jl) = a_i(ji,jj,jl+1) 618 zdvice(ji,jj,jl) = v_i(ji,jj,jl+1) 619 ENDIF 620 END DO 621 END DO 622 623 IF(lk_mpp) CALL mpp_max( zshiftflag ) ! clem: for reproducibility ??? 624 625 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 626 CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 617 IF( a_i(ji,jj,jl+1) > epsi10 .AND. v_i(ji,jj,jl+1) <= (a_i(ji,jj,jl+1) * hi_max(jl)) ) THEN 618 nidx = nidx + 1 619 idxice( nidx ) = (jj - 1) * jpi + ji 620 ENDIF 621 ENDDO 622 ENDDO 623 624 CALL tab_2d_1d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i(:,:,jl+1) ) ! jl+1 is ok 625 CALL tab_2d_1d( nidx, idxice(1:nidx), v_i_1d (1:nidx), v_i(:,:,jl+1) ) ! jl+1 is ok 626 DO ji = 1, nidx 627 jdonor(ji,jl) = jl + 1 628 zdaice(ji,jl) = a_i_1d(ji) 629 zdvice(ji,jl) = v_i_1d(ji) 630 END DO 631 632 IF( nidx > 0 ) THEN 633 CALL lim_itd_shiftice( jdonor(1:nidx,:), zdaice(1:nidx,:), zdvice(1:nidx,:) ) ! Shift jl+1=>jl 627 634 ! Reset shift parameters 628 zdonor(:,:,jl) = 0629 zdaice( :,:,jl) = 0._wp630 zdvice( :,:,jl) = 0._wp635 jdonor(1:nidx,jl) = 0 636 zdaice(1:nidx,jl) = 0._wp 637 zdvice(1:nidx,jl) = 0._wp 631 638 ENDIF 632 639 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
r8321 r8369 110 110 111 111 IF ( .NOT. ln_pnd ) THEN 112 WRITE(numout,*)113 WRITE(numout,*) ' Melt ponds are not activated '114 WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '115 WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '112 IF(lwp) WRITE(numout,*) 113 IF(lwp) WRITE(numout,*) ' Melt ponds are not activated ' 114 IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. ' 115 IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero ' 116 116 ln_pnd_rad = .FALSE. 117 117 ln_pnd_fw = .FALSE. … … 130 130 131 131 IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN 132 WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false '132 IF(lwp) WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false ' 133 133 ln_pnd_fw = .FALSE. 134 134 IF(lwp) THEN ! control print … … 138 138 139 139 IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN 140 WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '141 WRITE(numout,*) ' Run aborted '140 IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 ' 141 IF(lwp) WRITE(numout,*) ' Run aborted ' 142 142 CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' ) 143 143 ENDIF -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90
r8342 r8369 13 13 USE par_kind 14 14 USE par_oce 15 USE ice, ONLY : jpl 15 16 16 17 IMPLICIT NONE 17 18 PRIVATE 18 19 19 PUBLIC tab_2d_1d ! called by limthd 20 PUBLIC tab_1d_2d ! called by limthd 20 PUBLIC tab_3d_2d 21 PUBLIC tab_2d_1d 22 PUBLIC tab_2d_3d 23 PUBLIC tab_1d_2d 21 24 22 25 !!---------------------------------------------------------------------- … … 26 29 !!---------------------------------------------------------------------- 27 30 CONTAINS 31 32 33 SUBROUTINE tab_3d_2d( ndim1d, tab_ind, tab1d, tab2d ) 34 !!---------------------------------------------------------------------- 35 !! *** ROUTINE tab_2d_1d *** 36 !!---------------------------------------------------------------------- 37 INTEGER , INTENT(in ) :: ndim1d ! 1d size 38 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 39 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in ) :: tab2d ! input 2D field 40 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT( out) :: tab1d ! output 1D field 41 ! 42 INTEGER :: jl, jn, jid, jjd 43 !!---------------------------------------------------------------------- 44 DO jl = 1, jpl 45 DO jn = 1, ndim1d 46 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 47 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 48 tab1d(jn,jl) = tab2d(jid,jjd,jl) 49 END DO 50 END DO 51 END SUBROUTINE tab_3d_2d 28 52 29 53 SUBROUTINE tab_2d_1d( ndim1d, tab_ind, tab1d, tab2d ) … … 45 69 END SUBROUTINE tab_2d_1d 46 70 71 SUBROUTINE tab_2d_3d( ndim1d, tab_ind, tab1d, tab2d ) 72 !!---------------------------------------------------------------------- 73 !! *** ROUTINE tab_2d_1d *** 74 !!---------------------------------------------------------------------- 75 INTEGER , INTENT(in ) :: ndim1d ! 1D size 76 INTEGER , DIMENSION(ndim1d) , INTENT(in ) :: tab_ind ! input index 77 REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in ) :: tab1d ! input 1D field 78 REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT( out) :: tab2d ! output 2D field 79 ! 80 INTEGER :: jl, jn, jid, jjd 81 !!---------------------------------------------------------------------- 82 DO jl = 1, jpl 83 DO jn = 1, ndim1d 84 jid = MOD( tab_ind(jn) - 1 , jpi ) + 1 85 jjd = ( tab_ind(jn) - 1 ) / jpi + 1 86 tab2d(jid,jjd,jl) = tab1d(jn,jl) 87 END DO 88 END DO 89 END SUBROUTINE tab_2d_3d 47 90 48 91 SUBROUTINE tab_1d_2d( ndim1d, tab_ind, tab1d, tab2d ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r8355 r8369 237 237 IF( nidx > 0 ) THEN ! If there is no ice, do nothing. 238 238 ! 239 s_i_new (:) = 0._wp ; dh_s_tot (:) = 0._wp 239 s_i_new (:) = 0._wp ; dh_s_tot (:) = 0._wp ! --- some init --- ! 240 240 dh_i_surf (:) = 0._wp ; dh_i_bott(:) = 0._wp 241 241 dh_snowice(:) = 0._wp ; dh_i_sub (:) = 0._wp 242 242 243 CALL lim_thd_1d2d( jl, 1 ) 244 ! 245 DO jk = 1, nlay_i 246 WHERE( ht_i_1d( :)>0._wp ) e_i_1d(:,jk) = e_i_1d(:,jk) / (ht_i_1d(:) * a_i_1d(:)) * nlay_i243 CALL lim_thd_1d2d( jl, 1 ) ! --- Move to 1D arrays --- ! 244 ! 245 DO jk = 1, nlay_i ! --- Change units from J/m2 to J/m3 --- ! 246 WHERE( ht_i_1d(1:nidx)>0._wp ) e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) / (ht_i_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_i 247 247 ENDDO 248 248 DO jk = 1, nlay_s 249 WHERE( ht_s_1d( :)>0._wp ) e_s_1d(:,jk) = e_s_1d(:,jk) / (ht_s_1d(:) * a_i_1d(:)) * nlay_s249 WHERE( ht_s_1d(1:nidx)>0._wp ) e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) / (ht_s_1d(1:nidx) * a_i_1d(1:nidx)) * nlay_s 250 250 ENDDO 251 251 ! 252 IF( ln_limdH ) CALL lim_thd_dif ! --- Ice/Snow Temperature profile --- !253 ! 254 IF( ln_limdH ) CALL lim_thd_dh ! --- Ice/Snow thickness --- !252 IF( ln_limdH ) CALL lim_thd_dif ! --- Ice/Snow Temperature profile --- ! 253 ! 254 IF( ln_limdH ) CALL lim_thd_dh ! --- Ice/Snow thickness --- ! 255 255 ! 256 256 IF( ln_limdH ) CALL lim_thd_ent( e_i_1d(1:nidx,:) ) ! --- Ice enthalpy remapping --- ! 257 257 ! 258 CALL lim_thd_sal ! --- Ice salinity --- !259 ! 260 CALL lim_thd_temp ! --- temperature update --- !258 CALL lim_thd_sal ! --- Ice salinity --- ! 259 ! 260 CALL lim_thd_temp ! --- temperature update --- ! 261 261 ! 262 262 IF( ln_limdH ) THEN 263 263 IF ( ( nn_monocat == 1 .OR. nn_monocat == 4 ) .AND. jpl == 1 ) THEN 264 CALL lim_thd_lam ! --- extra lateral melting if monocat --- !264 CALL lim_thd_lam ! --- extra lateral melting if monocat --- ! 265 265 END IF 266 266 END IF 267 267 ! 268 DO jk = 1, nlay_i 269 e_i_1d( :,jk) = e_i_1d(:,jk) * ht_i_1d(:) * a_i_1d(:) * r1_nlay_i268 DO jk = 1, nlay_i ! --- Change units from J/m3 to J/m2 --- ! 269 e_i_1d(1:nidx,jk) = e_i_1d(1:nidx,jk) * ht_i_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_i 270 270 ENDDO 271 271 DO jk = 1, nlay_s 272 e_s_1d( :,jk) = e_s_1d(:,jk) * ht_s_1d(:) * a_i_1d(:) * r1_nlay_s272 e_s_1d(1:nidx,jk) = e_s_1d(1:nidx,jk) * ht_s_1d(1:nidx) * a_i_1d(1:nidx) * r1_nlay_s 273 273 ENDDO 274 274 ! 275 CALL lim_thd_1d2d( jl, 2 ) 275 CALL lim_thd_1d2d( jl, 2 ) ! --- Move to 2D arrays --- ! 276 276 ! 277 277 IF( lk_mpp ) CALL mpp_comm_free( ncomm_ice ) !RB necessary ?? … … 288 288 289 289 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_thd_da', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 290 IF( ln_limdA) CALL lim_thd_da 290 IF( ln_limdA) CALL lim_thd_da ! --- lateral melting --- ! 291 291 292 292 ! Change thickness to volume … … 441 441 ! 442 442 CALL tab_2d_1d( nidx, idxice(1:nidx), qprec_ice_1d(1:nidx), qprec_ice ) 443 CALL tab_2d_1d( nidx, idxice(1:nidx), qevap_ice_1d(1:nidx), qevap_ice(:,:,jl) )444 443 CALL tab_2d_1d( nidx, idxice(1:nidx), qsr_ice_1d (1:nidx), qsr_ice(:,:,jl) ) 445 444 CALL tab_2d_1d( nidx, idxice(1:nidx), fr1_i0_1d (1:nidx), fr1_i0 ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_da.F90
r8360 r8369 163 163 164 164 DO ji = 1, nidx 165 ! decrease of concentration for the category jl 166 ! each category contributes to melting in proportion to its concentration 167 zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) 168 169 ! Contribution to salt flux 170 sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic * ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice 171 172 ! Contribution to heat flux into the ocean [W.m-2], (<0) 173 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 174 175 ! Contribution to mass flux 176 wfx_lam_1d(ji) = wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) ) 177 178 !! adjust e_i ??? 179 e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 180 e_s_1d(ji,1) = e_s_1d(ji,1) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 181 182 ! new concentration 183 a_i_1d(ji) = a_i_1d(ji) + zda 184 185 ! ensure that ht_i = 0 where a_i = 0 186 IF( a_i_1d(ji) == 0._wp ) THEN 187 ht_i_1d(ji) = 0._wp 188 ht_s_1d(ji) = 0._wp 189 ENDIF 190 END DO 191 165 IF( a_i_1d(ji) > 0._wp ) THEN 166 ! decrease of concentration for the category jl 167 ! each category contributes to melting in proportion to its concentration 168 zda = zda_tot(ji) * a_i_1d(ji) / at_i_1d(ji) 169 170 ! Contribution to salt flux 171 sfx_lam_1d(ji) = sfx_lam_1d(ji) - rhoic * ht_i_1d(ji) * zda * sm_i_1d(ji) * r1_rdtice 172 173 ! Contribution to heat flux into the ocean [W.m-2], (<0) 174 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zda_tot(ji) / at_i_1d(ji) * SUM( e_i_1d(ji,1:nlay_i) + e_s_1d(ji,1) ) * r1_rdtice 175 176 ! Contribution to mass flux 177 wfx_lam_1d(ji) = wfx_lam_1d(ji) - zda * r1_rdtice * ( rhoic * ht_i_1d(ji) + rhosn * ht_s_1d(ji) ) 178 179 !! adjust e_i ??? 180 e_i_1d(ji,1:nlay_i) = e_i_1d(ji,1:nlay_i) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 181 e_s_1d(ji,1) = e_s_1d(ji,1) * ( 1._wp + zda_tot(ji) / at_i_1d(ji) ) 182 183 ! new concentration 184 a_i_1d(ji) = a_i_1d(ji) + zda 185 186 ! ensure that ht_i = 0 where a_i = 0 187 IF( a_i_1d(ji) == 0._wp ) THEN 188 ht_i_1d(ji) = 0._wp 189 ht_s_1d(ji) = 0._wp 190 ENDIF 191 ENDIF 192 END DO 193 192 194 CALL tab_1d_2d( nidx, idxice(1:nidx), a_i_1d (1:nidx), a_i (:,:,jl) ) 193 195 CALL tab_1d_2d( nidx, idxice(1:nidx), ht_i_1d (1:nidx), ht_i(:,:,jl) ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
r8342 r8369 631 631 IF( nn_icesal == 1 .OR. nn_icesal == 3 ) THEN 632 632 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_rdtice & ! put back sss_m into the ocean 633 & - sm_i_1d(ji) 633 & - sm_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice ! and get rn_icesal from the ocean 634 634 ENDIF 635 635 … … 665 665 666 666 ! --- ensure that a_i = 0 where ht_i = 0 --- 667 WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp 668 667 DO ji = 1, nidx 668 IF( ht_i_1d(ji) == 0._wp ) a_i_1d(ji) = 0._wp 669 END DO 670 669 671 CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema ) 670 672 CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zeh_i ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
r8342 r8369 735 735 DO ji = 1, nidx 736 736 zfac = 1. / MAX( epsi10 , rn_cdsn * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) ) 737 t_si_1d(ji) = ( rn_cdsn * zh_i(ji) * t_s_1d(ji,1) + & 738 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac 739 END DO 740 WHERE( ( zh_s < 1.0e-3 ) ) ; t_si_1d(:) = t_su_1d(:) ; END WHERE 737 IF( zh_s(ji) >= 1.e-3 ) THEN 738 t_si_1d(ji) = ( rn_cdsn * zh_i(ji) * t_s_1d(ji,1) + & 739 & ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) * zfac 740 ELSE 741 t_si_1d(ji) = t_su_1d(ji) 742 ENDIF 743 END DO 741 744 ! END MV SIMIP 2016 742 745 -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
r8342 r8369 12 12 USE in_out_manager ! I/O manager 13 13 USE lib_mpp ! MPP library 14 USE ice, ONLY : nlay_i, nlay_s 14 USE ice, ONLY : nlay_i, nlay_s, jpl 15 15 16 16 IMPLICIT NONE … … 19 19 PUBLIC thd_ice_alloc ! Routine called by nemogcm.F90 20 20 21 !!---------------------- -------22 !! * Share1D Module variables23 !!---------------------- -------21 !!---------------------- 22 !! * 1D Module variables 23 !!---------------------- 24 24 !: In ice thermodynamics, to spare memory, the vectors are folded 25 25 !: from 1D to 2D vectors. The following variables, with ending _1d … … 90 90 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: evap_ice_1d !: <==> the 2D evap_ice 91 91 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qprec_ice_1d !: <==> the 2D qprec_ice 92 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: qevap_ice_1d !: <==> the 3D qevap_ice93 92 ! ! to reintegrate longwave flux inside the ice thermodynamics 94 93 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: i0 !: fraction of radiation transmitted to the ice … … 98 97 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: t_si_1d !: <==> the 2D t_si 99 98 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_i_1d !: <==> the 2D a_i 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_1d !: <==> the 2D ht_s 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_1d !: <==> the 2D ht_i 99 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ib_1d !: <==> the 2D a_i_b 100 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_i_1d !: <==> the 2D ht_i 101 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_ib_1d !: <==> the 2D ht_i_b 102 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: ht_s_1d !: <==> the 2D ht_s 102 103 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_su !: Surface Conduction flux 103 104 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: fc_bo_i !: Bottom Conduction flux … … 109 110 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sm_i_1d !: Ice bulk salinity [ppt] 110 111 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: s_i_new !: Salinity of new ice at the bottom 112 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: a_ip_1d !: 113 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_ip_1d !: 114 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: v_i_1d !: 111 115 112 116 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_s_1d !: corresponding to the 2D var t_s … … 127 131 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sss_1d 128 132 133 ! 134 !!---------------------- 135 !! * 2D Module variables 136 !!---------------------- 137 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_i_2d 138 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_i_2d 139 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_s_2d 140 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: oa_i_2d 141 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smv_i_2d 142 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ip_2d 143 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: v_ip_2d 144 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_su_2d 145 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_i_2d 146 147 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: a_ib_2d 148 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_ib_2d 149 129 150 INTEGER , PUBLIC :: jiindex_1d ! 1D index of debugging point 130 151 … … 141 162 !!---------------------------------------------------------------------! 142 163 INTEGER :: thd_ice_alloc ! return value 143 INTEGER :: ierr( 6), ii164 INTEGER :: ierr(7), ii 144 165 !!---------------------------------------------------------------------! 145 166 ierr(:) = 0 … … 164 185 & wfx_snw_sub_1d(jpij), wfx_ice_sub_1d(jpij), wfx_err_sub_1d(jpij) , & 165 186 & wfx_lam_1d(jpij) , dqns_ice_1d(jpij) , evap_ice_1d (jpij), & 166 & qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0 (jpij) ,&187 & qprec_ice_1d(jpij), i0 (jpij) , & 167 188 & sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij), & 168 189 & sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij), & … … 170 191 ! 171 192 ii = ii + 1 172 ALLOCATE( t_su_1d (jpij) , t_si_1d (jpij) , a_i_1d (jpij) , ht_i_1d (jpij) ,&173 & ht_ s_1d (jpij) , fc_su (jpij) , fc_bo_i (jpij) ,&174 & dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub (jpij) ,&175 & dh_i_bott (jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new (jpij) ,&176 & STAT=ierr(ii) )193 ALLOCATE( t_su_1d (jpij) , t_si_1d (jpij) , a_i_1d (jpij) , a_ib_1d(jpij) , & 194 & ht_i_1d (jpij) , ht_ib_1d (jpij) , ht_s_1d (jpij) , fc_su (jpij) , fc_bo_i(jpij) , & 195 & dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) , & 196 & dh_i_bott(jpij) , dh_snowice(jpij) , sm_i_1d (jpij) , s_i_new(jpij) , & 197 & a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d (jpij) , STAT=ierr(ii) ) 177 198 ! 178 199 ii = ii + 1 … … 186 207 ii = ii + 1 187 208 ALLOCATE( sst_1d(jpij) , sss_1d(jpij) , STAT=ierr(ii) ) 209 ! 210 ii = ii + 1 211 ALLOCATE( a_i_2d(jpij,jpl) , a_ib_2d(jpij,jpl) , ht_i_2d(jpij,jpl) , ht_ib_2d(jpij,jpl) , & 212 & v_i_2d(jpij,jpl) ,v_s_2d(jpij,jpl) ,oa_i_2d(jpij,jpl) ,smv_i_2d(jpij,jpl) , & 213 & a_ip_2d(jpij,jpl) ,v_ip_2d(jpij,jpl) ,t_su_2d(jpij,jpl) , & 214 & STAT=ierr(ii) ) 188 215 189 216 thd_ice_alloc = MAXVAL( ierr(:) )
Note: See TracChangeset
for help on using the changeset viewer.
