Changeset 8355 for branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO
- Timestamp:
- 2017-07-20T10:20:49+02:00 (7 years ago)
- Location:
- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3
- Files:
-
- 5 edited
-
limcons.F90 (modified) (2 diffs)
-
limitd_th.F90 (modified) (35 diffs)
-
limthd.F90 (modified) (1 diff)
-
limupdate1.F90 (modified) (1 diff)
-
limupdate2.F90 (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r8341 r8355 28 28 PRIVATE 29 29 30 PUBLIC lim_column_sum31 PUBLIC lim_column_sum_energy32 PUBLIC lim_cons_check33 30 PUBLIC lim_cons_hsm 34 31 PUBLIC lim_cons_final … … 40 37 !!---------------------------------------------------------------------- 41 38 CONTAINS 42 43 SUBROUTINE lim_column_sum( ksum, pin, pout )44 !!-------------------------------------------------------------------45 !! *** ROUTINE lim_column_sum ***46 !!47 !! ** Purpose : Compute the sum of xin over nsum categories48 !!49 !! ** Method : Arithmetics50 !!51 !! ** Action : Gets xin(ji,jj,jl) and computes xout(ji,jj)52 !!---------------------------------------------------------------------53 INTEGER , INTENT(in ) :: ksum ! number of categories/layers54 REAL(wp), DIMENSION(:,:,:), INTENT(in ) :: pin ! input field55 REAL(wp), DIMENSION(:,:) , INTENT( out) :: pout ! output field56 !57 INTEGER :: jl ! dummy loop indices58 !!---------------------------------------------------------------------59 !60 pout(:,:) = pin(:,:,1)61 DO jl = 2, ksum62 pout(:,:) = pout(:,:) + pin(:,:,jl)63 END DO64 !65 END SUBROUTINE lim_column_sum66 67 68 SUBROUTINE lim_column_sum_energy( ksum, klay, pin, pout)69 !!-------------------------------------------------------------------70 !! *** ROUTINE lim_column_sum_energy ***71 !!72 !! ** Purpose : Compute the sum of xin over nsum categories73 !! and nlay layers74 !!75 !! ** Method : Arithmetics76 !!---------------------------------------------------------------------77 INTEGER , INTENT(in ) :: ksum !: number of categories78 INTEGER , INTENT(in ) :: klay !: number of vertical layers79 REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl), INTENT(in ) :: pin !: input field80 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pout !: output field81 !82 INTEGER :: jk, jl ! dummy loop indices83 !!---------------------------------------------------------------------84 !85 pout(:,:) = 0._wp86 DO jl = 1, ksum87 DO jk = 2, klay88 pout(:,:) = pout(:,:) + pin(:,:,jk,jl)89 END DO90 END DO91 !92 END SUBROUTINE lim_column_sum_energy93 94 95 SUBROUTINE lim_cons_check( px1, px2, pmax_err, cd_fieldid )96 !!-------------------------------------------------------------------97 !! *** ROUTINE lim_cons_check ***98 !!99 !! ** Purpose : Test the conservation of a certain variable100 !! For each physical grid cell, check that initial101 !! and final values102 !! of a conserved field are equal to within a small value.103 !!104 !! ** Method :105 !!---------------------------------------------------------------------106 REAL(wp), DIMENSION(:,:), INTENT(in ) :: px1 !: initial field107 REAL(wp), DIMENSION(:,:), INTENT(in ) :: px2 !: final field108 REAL(wp) , INTENT(in ) :: pmax_err !: max allowed error109 CHARACTER(len=15) , INTENT(in ) :: cd_fieldid !: field identifyer110 !111 INTEGER :: ji, jj ! dummy loop indices112 INTEGER :: inb_error ! number of g.c where there is a cons. error113 LOGICAL :: llconserv_err ! = .true. if conservation check failed114 REAL(wp) :: zmean_error ! mean error on error points115 !!---------------------------------------------------------------------116 !117 IF(lwp) WRITE(numout,*) ' lim_cons_check '118 IF(lwp) WRITE(numout,*) ' ~~~~~~~~~~~~~~ '119 120 llconserv_err = .FALSE.121 inb_error = 0122 zmean_error = 0._wp123 IF( MAXVAL( px2(:,:) - px1(:,:) ) > pmax_err ) llconserv_err = .TRUE.124 125 IF( llconserv_err ) THEN126 DO jj = 1, jpj127 DO ji = 1, jpi128 IF( ABS( px2(ji,jj) - px1(ji,jj) ) > pmax_err ) THEN129 inb_error = inb_error + 1130 zmean_error = zmean_error + ABS( px2(ji,jj) - px1(ji,jj) )131 !132 IF(lwp) THEN133 WRITE (numout,*) ' ALERTE 99 '134 WRITE (numout,*) ' Conservation error: ', cd_fieldid135 WRITE (numout,*) ' Point : ', ji, jj136 WRITE (numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj)137 WRITE (numout,*) ' Initial value : ', px1(ji,jj)138 WRITE (numout,*) ' Final value : ', px2(ji,jj)139 WRITE (numout,*) ' Difference : ', px2(ji,jj) - px1(ji,jj)140 ENDIF141 ENDIF142 END DO143 END DO144 !145 ENDIF146 IF(lk_mpp) CALL mpp_sum( inb_error )147 IF(lk_mpp) CALL mpp_sum( zmean_error )148 !149 IF( inb_error > 0 .AND. lwp ) THEN150 zmean_error = zmean_error / REAL( inb_error, wp )151 WRITE(numout,*) ' Conservation check for : ', cd_fieldid152 WRITE(numout,*) ' Number of error points : ', inb_error153 WRITE(numout,*) ' Mean error on these pts: ', zmean_error154 ENDIF155 !156 END SUBROUTINE lim_cons_check157 158 39 159 40 SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r8233 r8355 44 44 CONTAINS 45 45 46 SUBROUTINE lim_itd_th_rem( k lbnd, kubnd, kt )46 SUBROUTINE lim_itd_th_rem( kt ) 47 47 !!------------------------------------------------------------------ 48 48 !! *** ROUTINE lim_itd_th_rem *** … … 55 55 !! References : W.H. Lipscomb, JGR 2001 56 56 !!------------------------------------------------------------------ 57 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point58 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied59 57 INTEGER , INTENT (in) :: kt ! Ocean time step 60 58 ! … … 65 63 REAL(wp) :: zx2, zwk2, zda0, zetamax ! - - 66 64 REAL(wp) :: zx3 67 CHARACTER (len = 15) :: fieldid68 65 69 66 INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index … … 75 72 REAL(wp), POINTER, DIMENSION(:,:,:) :: hR ! left boundary for the ITD for each thickness 76 73 REAL(wp), POINTER, DIMENSION(:,:,:) :: zht_i_b ! old ice thickness 77 REAL(wp), POINTER, DIMENSION(:,:,:) :: dummy_es78 74 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! local increment of ice area and volume 79 REAL(wp), POINTER, DIMENSION(:) :: zvetamin, zvetamax ! maximum values for etas80 75 INTEGER , POINTER, DIMENSION(:) :: nind_i, nind_j ! compressed indices for i/j directions 81 76 INTEGER :: nbrem ! number of cells with ice to transfer 82 77 REAL(wp) :: zslope ! used to compute local thermodynamic "speeds" 83 78 REAL(wp), POINTER, DIMENSION(:,:) :: zhb0, zhb1 ! category boundaries for thinnes categories 84 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories85 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories86 REAL(wp), POINTER, DIMENSION(:,:) :: et_i_init, et_i_final ! ice energy summed over categories87 REAL(wp), POINTER, DIMENSION(:,:) :: et_s_init, et_s_final ! snow energy summed over categories88 79 INTEGER , POINTER, DIMENSION(:,:) :: zremap_flag ! compute remapping or not ???? 89 80 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhbnew ! new boundaries of ice categories … … 92 83 CALL wrk_alloc( jpi,jpj, zremap_flag ) 93 84 CALL wrk_alloc( jpi,jpj,jpl-1, zdonor ) 94 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b , dummy_es)85 CALL wrk_alloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 95 86 CALL wrk_alloc( jpi,jpj,jpl-1, zdaice, zdvice ) 96 87 CALL wrk_alloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 97 CALL wrk_alloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )98 88 CALL wrk_alloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 99 CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 100 101 !!---------------------------------------------------------------------------------------------- 102 !! 0) Conservation checkand changes in each ice category 103 !!---------------------------------------------------------------------------------------------- 104 IF( con_i ) THEN 105 CALL lim_column_sum (jpl, v_i, vt_i_init) 106 CALL lim_column_sum (jpl, v_s, vt_s_init) 107 CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init) 108 dummy_es(:,:,:) = e_s(:,:,1,:) 109 CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init) 110 ENDIF 89 CALL wrk_alloc( jpi,jpj, zhb0, zhb1 ) 111 90 112 91 !!---------------------------------------------------------------------------------------------- … … 117 96 WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution' 118 97 WRITE(numout,*) '~~~~~~~~~~~~~~~' 119 WRITE(numout,*) ' klbnd : ', klbnd120 WRITE(numout,*) ' kubnd : ', kubnd121 98 ENDIF 122 99 123 100 zdhice(:,:,:) = 0._wp 124 DO jl = klbnd, kubnd101 DO jl = 1, jpl 125 102 DO jj = 1, jpj 126 103 DO ji = 1, jpi … … 134 111 END DO 135 112 136 !-----------------------------------------------------------------------------------------------137 ! 2) Compute fractional ice area in each grid cell138 !-----------------------------------------------------------------------------------------------139 at_i(:,:) = 0._wp140 DO jl = klbnd, kubnd 141 at_i(:,:) = at_i(:,:) + a_i(:,:,jl)142 END DO113 !! !----------------------------------------------------------------------------------------------- 114 !! ! 2) Compute fractional ice area in each grid cell 115 !! !----------------------------------------------------------------------------------------------- 116 !! at_i(:,:) = 0._wp 117 !! DO jl = 1, jpl 118 !! at_i(:,:) = at_i(:,:) + a_i(:,:,jl) 119 !! END DO 143 120 144 121 !----------------------------------------------------------------------------------------------- … … 160 137 161 138 !----------------------------------------------------------------------------------------------- 162 ! 4) Compute new category boundaries 163 !----------------------------------------------------------------------------------------------- 164 !- 4.1 Compute category boundaries 139 ! 4) Compute new category boundaries: 1:jpl-1 140 !----------------------------------------------------------------------------------------------- 165 141 zhbnew(:,:,:) = 0._wp 166 142 167 DO jl = klbnd, kubnd - 1 168 DO ji = 1, nbrem 169 ii = nind_i(ji) 170 ij = nind_j(ji) 171 ! 172 zhbnew(ii,ij,jl) = hi_max(jl) 173 IF ( a_i_b(ii,ij,jl) > epsi10 .AND. a_i_b(ii,ij,jl+1) > epsi10 ) THEN 174 !interpolate between adjacent category growth rates 175 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 176 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + zslope * ( hi_max(jl) - zht_i_b(ii,ij,jl) ) 177 ELSEIF( a_i_b(ii,ij,jl) > epsi10) THEN 178 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 179 ELSEIF( a_i_b(ii,ij,jl+1) > epsi10) THEN 180 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 181 ENDIF 182 END DO 183 184 !- 4.2 Check that each zhbnew lies between adjacent values of ice thickness 185 DO ji = 1, nbrem 186 ii = nind_i(ji) 187 ij = nind_j(ji) 188 189 ! clem: we do not want ht_i to be too close to either HR or HL otherwise a division by nearly 0 is possible 190 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 191 IF ( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) THEN 192 zremap_flag(ii,ij) = 0 193 ELSEIF( a_i(ii,ij,jl+1) > epsi10 .AND. ht_i(ii,ij,jl+1) < ( zhbnew(ii,ij,jl) + epsi10 ) ) THEN 194 zremap_flag(ii,ij) = 0 195 ENDIF 196 197 !- 4.3 Check that each zhbnew does not exceed maximal values hi_max 198 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 199 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 200 ! clem bug: why is not the following instead? 201 !!IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 202 !!IF( zhbnew(ii,ij,jl) > hi_max(jl ) ) zremap_flag(ii,ij) = 0 203 204 END DO 205 206 END DO 207 143 IF( nbrem > 0 ) THEN 144 DO jl = 1, jpl - 1 145 DO ji = 1, nbrem 146 ii = nind_i(ji) 147 ij = nind_j(ji) 148 ! 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) /= 0 152 zslope = ( zdhice(ii,ij,jl+1) - zdhice(ii,ij,jl) ) / ( zht_i_b(ii,ij,jl+1) - zht_i_b(ii,ij,jl) ) 153 zhbnew(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*dt 156 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 159 zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 160 ELSE !! a_i_b(jl+1) & a_i_b(jl) = 0 161 zhbnew(ii,ij,jl) = hi_max(jl) 162 ENDIF 163 164 ! --- 2 conditions for remapping --- ! 165 ! 1) hn(t+1)+espi < Hn* < hn+1(t+1)-epsi 166 ! 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_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 168 IF( a_i(ii,ij,jl ) > epsi10 .AND. ht_i(ii,ij,jl ) > ( zhbnew(ii,ij,jl) - epsi10 ) ) zremap_flag(ii,ij) = 0 169 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 170 171 ! 2) Hn-1 < Hn* < Hn+1 172 IF( zhbnew(ii,ij,jl) < hi_max(jl-1) ) zremap_flag(ii,ij) = 0 173 IF( zhbnew(ii,ij,jl) > hi_max(jl+1) ) zremap_flag(ii,ij) = 0 174 175 END DO 176 END DO 177 ENDIF 178 208 179 !----------------------------------------------------------------------------------------------- 209 180 ! 5) Identify cells where ITD is to be remapped … … 221 192 222 193 !----------------------------------------------------------------------------------------------- 223 ! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories194 ! 6) Compute new category boundaries: jpl 224 195 !----------------------------------------------------------------------------------------------- 225 196 DO jj = 1, jpj … … 228 199 zhb1(ji,jj) = hi_max(1) 229 200 230 IF( a_i(ji,jj,kubnd) > epsi10 ) THEN 231 zhbnew(ji,jj,kubnd) = MAX( hi_max(kubnd-1), 3._wp * ht_i(ji,jj,kubnd) - 2._wp * zhbnew(ji,jj,kubnd-1) ) 201 ! boundary for jpl 202 IF( a_i(ji,jj,jpl) > epsi10 ) THEN 203 zhbnew(ji,jj,jpl) = MAX( hi_max(jpl-1), 3._wp * ht_i(ji,jj,jpl) - 2._wp * zhbnew(ji,jj,jpl-1) ) 232 204 ELSE 233 !clem bug zhbnew(ji,jj,kubnd) = hi_max(kubnd) 234 zhbnew(ji,jj,kubnd) = hi_max(kubnd-1) ! not used anyway 235 ENDIF 236 237 ! clem: we do not want ht_i_b to be too close to either HR or HL otherwise a division by nearly 0 is possible 238 ! in lim_itd_fitline in the case (HR-HL) = 3(Hice - HL) or = 3(HR - Hice) 239 IF ( zht_i_b(ji,jj,klbnd) < ( zhb0(ji,jj) + epsi10 ) ) THEN 240 zremap_flag(ji,jj) = 0 241 ELSEIF( zht_i_b(ji,jj,klbnd) > ( zhb1(ji,jj) - epsi10 ) ) THEN 242 zremap_flag(ji,jj) = 0 243 ENDIF 205 !clem bug zhbnew(ji,jj,jpl) = hi_max(jpl) 206 zhbnew(ji,jj,jpl) = hi_max(jpl-1) ! not used anyway 207 ENDIF 208 209 ! 1 condition for remapping for the 1st category 210 ! H0+epsi < h1(t) < H1-epsi 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 (HR-HL) = 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 244 215 245 216 END DO … … 250 221 !----------------------------------------------------------------------------------------------- 251 222 !- 7.1 g(h) for category 1 at start of time step 252 CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &253 & hR(:,:, klbnd), zremap_flag )254 255 !- 7.2 Area lost due to melting of thin ice (first category, klbnd)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) 256 227 DO ji = 1, nbrem 257 228 ii = nind_i(ji) 258 229 ij = nind_j(ji) 259 230 260 IF( a_i(ii,ij, klbnd) > epsi10 ) THEN261 262 zdh0 = zdhice(ii,ij, klbnd) !decrease of ice thickness in the lower category231 IF( a_i(ii,ij,1) > epsi10 ) THEN 232 233 zdh0 = zdhice(ii,ij,1) !decrease of ice thickness in the lower category 263 234 264 235 IF( zdh0 < 0.0 ) THEN !remove area from category 1 265 zdh0 = MIN( -zdh0, hi_max( klbnd) )236 zdh0 = MIN( -zdh0, hi_max(1) ) 266 237 !Integrate g(1) from 0 to dh0 to estimate area melted 267 zetamax = MIN( zdh0, hR(ii,ij, klbnd) ) - hL(ii,ij,klbnd)238 zetamax = MIN( zdh0, hR(ii,ij,1) ) - hL(ii,ij,1) 268 239 269 240 IF( zetamax > 0.0 ) THEN 270 241 zx1 = zetamax 271 242 zx2 = 0.5 * zetamax * zetamax 272 zda0 = g1(ii,ij, klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 ! ice area removed273 zdamax = a_i(ii,ij, klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! Constrain new thickness <= ht_i243 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 274 245 zda0 = MIN( zda0, zdamax ) ! ice area lost due to melting 275 246 ! of thin ice (zdamax > 0) 276 247 ! Remove area, conserving volume 277 ht_i(ii,ij, klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )278 a_i(ii,ij, klbnd) = a_i(ii,ij,klbnd) - zda0279 v_i(ii,ij, klbnd) = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?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) ! clem-useless ? 280 251 ENDIF 281 252 282 253 ELSE ! if ice accretion zdh0 > 0 283 254 ! zhbnew was 0, and is shifted to the right to account for thin ice growth in openwater (F0 = f1) 284 zhbnew(ii,ij, klbnd-1) = MIN( zdh0, hi_max(klbnd) )255 zhbnew(ii,ij,0) = MIN( zdh0, hi_max(1) ) 285 256 ENDIF 286 257 … … 290 261 291 262 !- 7.3 g(h) for each thickness category 292 DO jl = klbnd, kubnd263 DO jl = 1, jpl 293 264 CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), & 294 265 & g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag ) … … 296 267 297 268 !----------------------------------------------------------------------------------------------- 298 ! 8) Compute area and volume to be shifted across each boundary 299 !----------------------------------------------------------------------------------------------- 300 301 DO jl = klbnd, kubnd- 1269 ! 8) Compute area and volume to be shifted across each boundary (Eq. 18) 270 !----------------------------------------------------------------------------------------------- 271 272 DO jl = 1, jpl - 1 302 273 DO jj = 1, jpj 303 274 DO ji = 1, jpi … … 312 283 ij = nind_j(ji) 313 284 314 IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1315 ! left and right integration limits in eta space316 z vetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)317 z vetamax(ji) = MIN( zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)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 318 289 zdonor(ii,ij,jl) = jl 319 320 ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 321 ! left and right integration limits in eta space 322 zvetamin(ji) = 0.0 323 zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1) 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 324 293 zdonor(ii,ij,jl) = jl + 1 325 326 ENDIF 327 328 zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin 329 zetamin = zvetamin(ji) 294 ENDIF 295 zetamax = MAX( zetamax, zetamin ) ! no transfer if etamax < etamin 330 296 331 297 zx1 = zetamax - zetamin … … 346 312 !! 9) Shift ice between categories 347 313 !!---------------------------------------------------------------------------------------------- 348 CALL lim_itd_shiftice ( klbnd, kubnd,zdonor, zdaice, zdvice )314 CALL lim_itd_shiftice ( zdonor, zdaice, zdvice ) 349 315 350 316 !!---------------------------------------------------------------------------------------------- … … 366 332 END DO 367 333 368 !!----------------------------------------------------------------------------------------------369 !! 11) Conservation check370 !!----------------------------------------------------------------------------------------------371 IF ( con_i ) THEN372 CALL lim_column_sum (jpl, v_i, vt_i_final)373 fieldid = ' v_i : limitd_th '374 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)375 376 CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_final)377 fieldid = ' e_i : limitd_th '378 CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)379 380 CALL lim_column_sum (jpl, v_s, vt_s_final)381 fieldid = ' v_s : limitd_th '382 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)383 384 dummy_es(:,:,:) = e_s(:,:,1,:)385 CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)386 fieldid = ' e_s : limitd_th '387 CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)388 ENDIF389 390 334 CALL wrk_dealloc( jpi,jpj, zremap_flag ) 391 335 CALL wrk_dealloc( jpi,jpj,jpl-1, zdonor ) 392 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b , dummy_es)336 CALL wrk_dealloc( jpi,jpj,jpl, zdhice, g0, g1, hL, hR, zht_i_b ) 393 337 CALL wrk_dealloc( jpi,jpj,jpl-1, zdaice, zdvice ) 394 338 CALL wrk_dealloc( jpi,jpj,jpl+1, zhbnew, kkstart = 0 ) 395 CALL wrk_dealloc( (jpi+1)*(jpj+1), zvetamin, zvetamax )396 339 CALL wrk_dealloc( (jpi+1)*(jpj+1), nind_i, nind_j ) 397 CALL wrk_dealloc( jpi,jpj, zhb0, zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final)340 CALL wrk_dealloc( jpi,jpj, zhb0, zhb1 ) 398 341 399 342 END SUBROUTINE lim_itd_th_rem … … 404 347 !! *** ROUTINE lim_itd_fitline *** 405 348 !! 406 !! ** Purpose : fit g(h) with a line using area, volume constraints 407 !! 408 !! ** Method : Fit g(h) with a line, satisfying area and volume constraints. 409 !! To reduce roundoff errors caused by large values of g0 and g1, 410 !! we actually compute g(eta), where eta = h - hL, and hL is the 411 !! left boundary. 349 !! ** Purpose : build g(h) satisfying area and volume constraints (Eq. 6 and 9) 350 !! 351 !! ** Method : g(h) is linear and written as: g(eta) = g1(eta) + g0 352 !! with eta = h - HL 353 !! 412 354 !!------------------------------------------------------------------ 413 355 INTEGER , INTENT(in ) :: num_cat ! category index … … 436 378 hR(ji,jj) = HbR(ji,jj) 437 379 438 ! Change hL or hR if hice falls outside central third of range 439 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 440 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 441 442 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) 443 ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) 380 ! Change hL or hR if hice falls outside central third of range, 381 ! so that hice is in the central third of the range [HL HR] 382 zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) ) 383 zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) ) 384 385 IF ( hice(ji,jj) < zh13 ) THEN ; hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj) ! move HR to the left 386 ELSEIF( hice(ji,jj) > zh23 ) THEN ; hL(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hR(ji,jj) ! move HL to the right 444 387 ENDIF 445 388 … … 448 391 zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr 449 392 zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr 450 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) 451 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) 393 g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 ) ! Eq. 14 394 g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 ) ! Eq. 14 452 395 ! 453 396 ELSE ! remap_flag = .false. or a_i < epsi10 … … 464 407 465 408 466 SUBROUTINE lim_itd_shiftice( klbnd, kubnd,zdonor, zdaice, zdvice )409 SUBROUTINE lim_itd_shiftice( zdonor, zdaice, zdvice ) 467 410 !!------------------------------------------------------------------ 468 411 !! *** ROUTINE lim_itd_shiftice *** … … 473 416 !! ** Method : 474 417 !!------------------------------------------------------------------ 475 INTEGER , INTENT(in ) :: klbnd ! Start thickness category index point476 INTEGER , INTENT(in ) :: kubnd ! End point on which the the computation is applied477 418 INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(in ) :: zdonor ! donor category index 478 419 REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(inout) :: zdaice ! ice area transferred across boundary … … 508 449 !---------------------------------------------------------------------------------------------- 509 450 510 DO jl = klbnd, kubnd451 DO jl = 1, jpl 511 452 zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl) 512 453 END DO … … 516 457 !------------------------------------------------------------------------------- 517 458 518 DO jl = klbnd, kubnd- 1459 DO jl = 1, jpl - 1 519 460 nbrem = 0 520 461 DO jj = 1, jpj … … 539 480 ENDIF 540 481 541 !--------------542 482 ! Ice areas 543 !--------------544 483 a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 545 484 a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 546 485 547 !--------------548 486 ! Ice volumes 549 !--------------550 487 v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl) 551 488 v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 552 489 553 !--------------554 490 ! Snow volumes 555 !--------------556 491 zdvsnow = v_s(ii,ij,jl1) * zworka(ii,ij) 557 492 v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 558 493 v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow 559 494 560 !--------------------561 495 ! Snow heat content 562 !--------------------563 496 zdesnow = e_s(ii,ij,1,jl1) * zworka(ii,ij) 564 497 e_s(ii,ij,1,jl1) = e_s(ii,ij,1,jl1) - zdesnow 565 498 e_s(ii,ij,1,jl2) = e_s(ii,ij,1,jl2) + zdesnow 566 499 567 !--------------568 500 ! Ice age 569 !--------------570 501 zdo_aice = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 571 502 oa_i(ii,ij,jl1) = oa_i(ii,ij,jl1) - zdo_aice 572 503 oa_i(ii,ij,jl2) = oa_i(ii,ij,jl2) + zdo_aice 573 504 574 !--------------575 505 ! Ice salinity 576 !--------------577 506 zdsm_vice = smv_i(ii,ij,jl1) * zworka(ii,ij) 578 507 smv_i(ii,ij,jl1) = smv_i(ii,ij,jl1) - zdsm_vice 579 508 smv_i(ii,ij,jl2) = smv_i(ii,ij,jl2) + zdsm_vice 580 509 581 !---------------------582 510 ! Surface temperature 583 !---------------------584 511 zdaTsf = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 585 512 zaTsfn(ii,ij,jl1) = zaTsfn(ii,ij,jl1) - zdaTsf … … 588 515 ! MV MP 2016 589 516 IF ( nn_pnd_scheme > 0 ) THEN 590 !---------------------591 517 ! Pond fraction 592 !---------------------593 518 zdapnd = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 594 519 a_ip(ii,ij,jl1) = a_ip(ii,ij,jl1) - zdapnd 595 520 a_ip(ii,ij,jl2) = a_ip(ii,ij,jl2) + zdapnd 596 521 597 !---------------------598 522 ! Pond volume 599 !---------------------600 523 zdvpnd = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl) 601 524 v_ip(ii,ij,jl1) = v_ip(ii,ij,jl1) - zdvpnd … … 607 530 END DO 608 531 609 !------------------610 532 ! Ice heat content 611 !------------------612 613 533 DO jk = 1, nlay_i 614 534 DO ji = 1, nbrem … … 635 555 !----------------------------------------------------------------- 636 556 637 DO jl = klbnd, kubnd557 DO jl = 1, jpl 638 558 DO jj = 1, jpj 639 559 DO ji = 1, jpi … … 656 576 657 577 658 SUBROUTINE lim_itd_th_reb ( klbnd, kubnd )578 SUBROUTINE lim_itd_th_reb 659 579 !!------------------------------------------------------------------ 660 580 !! *** ROUTINE lim_itd_th_reb *** … … 664 584 !! ** Method : 665 585 !!------------------------------------------------------------------ 666 INTEGER , INTENT (in) :: klbnd ! Start thickness category index point667 INTEGER , INTENT (in) :: kubnd ! End point on which the the computation is applied668 !669 586 INTEGER :: ji,jj, jl ! dummy loop indices 670 587 INTEGER :: zshiftflag ! = .true. if ice must be shifted 671 CHARACTER (len = 15) :: fieldid672 588 673 589 INTEGER , POINTER, DIMENSION(:,:,:) :: zdonor ! donor category index 674 590 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdaice, zdvice ! ice area and volume transferred 675 676 REAL(wp), POINTER, DIMENSION(:,:) :: vt_i_init, vt_i_final ! ice volume summed over categories677 REAL(wp), POINTER, DIMENSION(:,:) :: vt_s_init, vt_s_final ! snow volume summed over categories678 591 !!------------------------------------------------------------------ 679 592 680 593 CALL wrk_alloc( jpi,jpj,jpl, zdonor ) ! interger 681 594 CALL wrk_alloc( jpi,jpj,jpl, zdaice, zdvice ) 682 CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )683 595 ! 684 IF( con_i ) THEN ! conservation check685 CALL lim_column_sum (jpl, v_i, vt_i_init)686 CALL lim_column_sum (jpl, v_s, vt_s_init)687 ENDIF688 689 !690 596 !------------------------------------------------------------------------------ 691 597 ! 1) Compute ice thickness. 692 598 !------------------------------------------------------------------------------ 693 DO jl = klbnd, kubnd599 DO jl = 1, jpl 694 600 DO jj = 1, jpj 695 601 DO ji = 1, jpi … … 707 613 ! Initialize shift arrays 708 614 !------------------------- 709 DO jl = klbnd, kubnd615 DO jl = 1, jpl 710 616 zdonor(:,:,jl) = 0 711 617 zdaice(:,:,jl) = 0._wp … … 717 623 !------------------------- 718 624 719 DO jl = klbnd, kubnd- 1 ! loop over category boundaries625 DO jl = 1, jpl - 1 ! loop over category boundaries 720 626 721 627 !--------------------------------------- … … 742 648 743 649 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 744 CALL lim_itd_shiftice( klbnd, kubnd,zdonor, zdaice, zdvice )650 CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 745 651 ! Reset shift parameters 746 652 zdonor(:,:,jl) = 0 … … 755 661 !---------------------------- 756 662 757 DO jl = kubnd- 1, 1, -1 ! loop over category boundaries663 DO jl = jpl - 1, 1, -1 ! loop over category boundaries 758 664 759 665 !----------------------------------------- … … 777 683 778 684 IF( zshiftflag == 1 ) THEN ! Shift ice between categories 779 CALL lim_itd_shiftice( klbnd, kubnd,zdonor, zdaice, zdvice )685 CALL lim_itd_shiftice( zdonor, zdaice, zdvice ) 780 686 ! Reset shift parameters 781 687 zdonor(:,:,jl) = 0 … … 785 691 786 692 END DO 787 788 !------------------------------------------------------------------------------789 ! 3) Conservation check790 !------------------------------------------------------------------------------791 792 IF( con_i ) THEN793 CALL lim_column_sum (jpl, v_i, vt_i_final)794 fieldid = ' v_i : limitd_reb '795 CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)796 797 CALL lim_column_sum (jpl, v_s, vt_s_final)798 fieldid = ' v_s : limitd_reb '799 CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)800 ENDIF801 693 ! 802 694 CALL wrk_dealloc( jpi,jpj,jpl, zdonor ) 803 695 CALL wrk_dealloc( jpi,jpj,jpl, zdaice, zdvice ) 804 CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final )805 696 806 697 END SUBROUTINE lim_itd_th_reb -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
r8342 r8355 319 319 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 320 320 321 IF( jpl > 1 ) CALL lim_itd_th_rem( 1, jpl,kt )321 IF( jpl > 1 ) CALL lim_itd_th_rem( kt ) 322 322 323 323 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
r8316 r8355 104 104 ! Rebin categories with thickness out of bounds 105 105 !---------------------------------------------------- 106 IF ( jpl > 1 ) CALL lim_itd_th_reb (1, jpl)106 IF ( jpl > 1 ) CALL lim_itd_th_reb 107 107 108 108 !----------------- -
branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
r8316 r8355 118 118 ! Rebin categories with thickness out of bounds 119 119 !---------------------------------------------------- 120 IF ( jpl > 1 ) CALL lim_itd_th_reb ( 1, jpl )120 IF ( jpl > 1 ) CALL lim_itd_th_reb 121 121 122 122 !-----------------
Note: See TracChangeset
for help on using the changeset viewer.
