- Timestamp:
- 2014-11-28T18:24:01+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO2_ice_shelves/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r4333 r4924 10 10 !! ! 2006-11 (X. Fettweis) Vectorized 11 11 !! 3.0 ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 12 !! 3.4 ! 2011-02 (G. Madec) dynamical allocation 13 !! - ! 2014-05 (C. Rousset) complete rewriting 13 14 !!---------------------------------------------------------------------- 14 15 #if defined key_lim3 … … 22 23 USE domain ! 23 24 USE phycst ! physical constants 25 USE sbc_oce ! Surface boundary condition: ocean fields 24 26 USE ice ! LIM variables 25 27 USE par_ice ! LIM parameters … … 34 36 PRIVATE 35 37 36 PUBLIC lim_thd_ent ! called by lim _thd38 PUBLIC lim_thd_ent ! called by limthd and limthd_lac 37 39 38 REAL(wp) :: epsi20 = 1.e-20_wp ! constant values 39 REAL(wp) :: epsi10 = 1.e-10_wp ! 40 REAL(wp) :: zzero = 0._wp ! 41 REAL(wp) :: zone = 1._wp ! 40 REAL(wp) :: epsi20 = 1.e-20 ! constant values 41 REAL(wp) :: epsi10 = 1.e-10 ! constant values 42 42 43 43 !!---------------------------------------------------------------------- … … 48 48 CONTAINS 49 49 50 SUBROUTINE lim_thd_ent( kideb, kiut, jl)50 SUBROUTINE lim_thd_ent( kideb, kiut, qnew ) 51 51 !!------------------------------------------------------------------- 52 52 !! *** ROUTINE lim_thd_ent *** 53 53 !! 54 54 !! ** Purpose : 55 !! This routine computes new vertical grids 56 !! in the ice and in the snow, and consistently redistributes 57 !! temperatures in the snow / ice. 55 !! This routine computes new vertical grids in the ice, 56 !! and consistently redistributes temperatures. 58 57 !! Redistribution is made so as to ensure to energy conservation 59 58 !! … … 61 60 !! ** Method : linear conservative remapping 62 61 !! 63 !! ** Steps : 1) Grid 64 !! 2) Switches 65 !! 3) Snow redistribution 66 !! 4) Ice enthalpy redistribution 67 !! 5) Ice salinity, recover temperature 62 !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 63 !! 2) linear remapping on the new layers 64 !! 65 !! ------------ cum0(0) ------------- cum1(0) 66 !! NEW ------------- 67 !! ------------ cum0(1) ==> ------------- 68 !! ... ------------- 69 !! ------------ ------------- 70 !! ------------ cum0(nlay_i+2) ------------- cum1(nlay_i) 71 !! 68 72 !! 69 73 !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 70 74 !!------------------------------------------------------------------- 71 75 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 72 INTEGER , INTENT(in) :: jl ! Thickness cateogry number73 76 74 INTEGER :: ji,jk ! dummy loop indices 75 INTEGER :: ii, ij , & ! dummy indices 76 ntop0 , & ! old layer top index 77 nbot1 , & ! new layer bottom index 78 ntop1 , & ! new layer top index 79 limsum , & ! temporary loop index 80 nlayi0,nlays0 , & ! old number of layers 81 maxnbot0 , & ! old layer bottom index 82 layer0, layer1 ! old/new layer indexes 77 REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew ! new enthlapies (J.m-3, remapped) 83 78 84 85 REAL(wp) :: & 86 ztmelts , & ! ice melting point 87 zqsnic , & ! enthalpy of snow ice layer 88 zhsnow , & ! temporary snow thickness variable 89 zswitch , & ! dummy switch argument 90 zfac1 , & ! dummy factor 91 zfac2 , & ! dummy factor 92 ztform , & !: bottom formation temperature 93 zaaa , & !: dummy factor 94 zbbb , & !: dummy factor 95 zccc , & !: dummy factor 96 zdiscrim !: dummy factor 97 98 INTEGER, POINTER, DIMENSION(:) :: snswi ! snow switch 99 INTEGER, POINTER, DIMENSION(:) :: nbot0 ! old layer bottom index 100 INTEGER, POINTER, DIMENSION(:) :: icsuind ! ice surface index 101 INTEGER, POINTER, DIMENSION(:) :: icsuswi ! ice surface switch 102 INTEGER, POINTER, DIMENSION(:) :: icboind ! ice bottom index 103 INTEGER, POINTER, DIMENSION(:) :: icboswi ! ice bottom switch 104 INTEGER, POINTER, DIMENSION(:) :: snicind ! snow ice index 105 INTEGER, POINTER, DIMENSION(:) :: snicswi ! snow ice switch 106 INTEGER, POINTER, DIMENSION(:) :: snind ! snow index 79 INTEGER :: ji ! dummy loop indices 80 INTEGER :: jk0, jk1 ! old/new layer indices 81 REAL(wp) :: zswitch 107 82 ! 108 REAL(wp), POINTER, DIMENSION(:) :: zh_i ! thickness of an ice layer 109 REAL(wp), POINTER, DIMENSION(:) :: zh_s ! thickness of a snow layer 110 REAL(wp), POINTER, DIMENSION(:) :: zqsnow ! enthalpy of the snow put in snow ice 111 REAL(wp), POINTER, DIMENSION(:) :: zdeltah ! temporary variable 112 REAL(wp), POINTER, DIMENSION(:) :: zqti_in, zqts_in 113 REAL(wp), POINTER, DIMENSION(:) :: zqti_fin, zqts_fin 114 115 REAL(wp), POINTER, DIMENSION(:,:) :: zm0 ! old layer-system vertical cotes 116 REAL(wp), POINTER, DIMENSION(:,:) :: qm0 ! old layer-system heat content 117 REAL(wp), POINTER, DIMENSION(:,:) :: z_s ! new snow system vertical cotes 118 REAL(wp), POINTER, DIMENSION(:,:) :: z_i ! new ice system vertical cotes 119 REAL(wp), POINTER, DIMENSION(:,:) :: zthick0 ! old ice thickness 120 REAL(wp), POINTER, DIMENSION(:,:) :: zhl0 ! old and new layer thicknesses 121 REAL(wp), POINTER, DIMENSION(:,:) :: zrl01 122 123 REAL(wp) :: zinda 83 REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces 84 REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces 85 REAL(wp), POINTER, DIMENSION(:) :: zhnew ! new layers thicknesses 124 86 !!------------------------------------------------------------------- 125 87 126 CALL wrk_alloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer 127 CALL wrk_alloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real 128 CALL wrk_alloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) 129 CALL wrk_alloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) 88 CALL wrk_alloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 89 CALL wrk_alloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 90 CALL wrk_alloc( jpij, zhnew ) 130 91 131 zthick0(:,:) = 0._wp 132 zm0 (:,:) = 0._wp 133 qm0 (:,:) = 0._wp 134 zrl01 (:,:) = 0._wp 135 zhl0 (:,:) = 0._wp 136 z_i (:,:) = 0._wp 137 z_s (:,:) = 0._wp 138 139 ! 140 !------------------------------------------------------------------------------| 141 ! 1) Grid | 142 !------------------------------------------------------------------------------| 143 nlays0 = nlay_s 144 nlayi0 = nlay_i 145 146 DO ji = kideb, kiut 147 zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i ) 148 zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 149 END DO 150 151 ! 152 !------------------------------------------------------------------------------| 153 ! 2) Switches | 154 !------------------------------------------------------------------------------| 155 ! 2.1 snind(ji), snswi(ji) 156 ! snow surface behaviour : computation of snind(ji)-snswi(ji) 157 ! snind(ji) : index which equals 158 ! 0 if snow is accumulating 159 ! 1 if 1st layer is melting 160 ! 2 if 2nd layer is melting ... 161 DO ji = kideb, kiut 162 snind (ji) = 0 163 zdeltah(ji) = 0._wp 164 ENDDO !ji 165 166 DO jk = 1, nlays0 92 !-------------------------------------------------------------------------- 93 ! 1) Cumulative integral of old enthalpy * thicnkess and layers interfaces 94 !-------------------------------------------------------------------------- 95 zqh_cum0(:,0:nlay_i+2) = 0._wp 96 zh_cum0 (:,0:nlay_i+2) = 0._wp 97 DO jk0 = 1, nlay_i+2 167 98 DO ji = kideb, kiut 168 snind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 169 + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 170 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 171 END DO ! ji 172 END DO ! jk 173 174 ! snswi(ji) : switch which value equals 1 if snow melts 175 ! 0 if not 176 DO ji = kideb, kiut 177 snswi(ji) = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 178 END DO ! ji 179 180 ! 2.2 icsuind(ji), icsuswi(ji) 181 ! ice surface behaviour : computation of icsuind(ji)-icsuswi(ji) 182 ! icsuind(ji) : index which equals 183 ! 0 if nothing happens at the surface 184 ! 1 if first layer is melting 185 ! 2 if 2nd layer is reached by melt ... 186 DO ji = kideb, kiut 187 icsuind(ji) = 0 188 zdeltah(ji) = 0._wp 189 END DO !ji 190 DO jk = 1, nlayi0 191 DO ji = kideb, kiut 192 icsuind(ji) = jk * NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 193 + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 194 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 195 END DO ! ji 196 ENDDO !jk 197 198 ! icsuswi(ji) : switch which equals 199 ! 1 if ice melts at the surface 200 ! 0 if not 201 DO ji = kideb, kiut 202 icsuswi(ji) = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 99 zqh_cum0(ji,jk0) = zqh_cum0(ji,jk0-1) + qh_i_old(ji,jk0-1) 100 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + h_i_old (ji,jk0-1) 101 ENDDO 203 102 ENDDO 204 103 205 ! 2.3 icboind(ji), icboswi(ji) 206 ! ice bottom behaviour : computation of icboind(ji)-icboswi(ji) 207 ! icboind(ji) : index which equals 208 ! 0 if accretion is on the way 209 ! 1 if last layer has started to melt 210 ! 2 if penultiem layer is melting ... and so on 211 ! N+1 if all layers melt and that snow transforms into ice 212 DO ji = kideb, kiut 213 icboind(ji) = 0 214 zdeltah(ji) = 0._wp 215 END DO 216 DO jk = nlayi0, 1, -1 217 DO ji = kideb, kiut 218 icboind(ji) = (nlayi0+1-jk) * NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 219 & + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji))))) 220 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 221 END DO 222 END DO 223 104 !------------------------------------ 105 ! 2) Interpolation on the new layers 106 !------------------------------------ 107 ! new layer thickesses 224 108 DO ji = kideb, kiut 225 ! case of total ablation with remaining snow 226 IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. & 227 ( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1 228 END DO 229 230 ! icboswi(ji) : switch which equals 231 ! 1 if ice accretion is on the way 232 ! 0 if ablation is on the way 233 DO ji = kideb, kiut 234 icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 235 END DO 236 237 ! 2.4 snicind(ji), snicswi(ji) 238 ! snow ice formation : calcul de snicind(ji)-snicswi(ji) 239 ! snicind(ji) : index which equals 240 ! 0 if no snow-ice forms 241 ! 1 if last layer of snow has started to melt 242 ! 2 if penultiem layer ... 243 DO ji = kideb, kiut 244 snicind(ji) = 0 245 zdeltah(ji) = 0._wp 246 END DO 247 DO jk = nlays0, 1, -1 248 DO ji = kideb, kiut 249 snicind(ji) = (nlays0+1-jk) & 250 * NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji) & 251 * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 252 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 253 END DO 254 END DO 255 256 ! snicswi(ji) : switch which equals 257 ! 1 if snow-ice forms 258 ! 0 if not 259 DO ji = kideb, kiut 260 snicswi(ji) = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 109 zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i ) 261 110 ENDDO 262 111 263 ! 264 !------------------------------------------------------------------------------| 265 ! 3) Snow redistribution | 266 !------------------------------------------------------------------------------| 267 ! 268 !------------- 269 ! Old profile 270 !------------- 271 272 ! by 'old', it is meant that layers coming from accretion are included, 273 ! and that interfacial layers which were partly melted are reduced 274 275 ! indexes of the vectors 276 !------------------------ 277 ntop0 = 1 278 maxnbot0 = 0 279 280 DO ji = kideb, kiut 281 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 282 ! cotes of the top of the layers 283 zm0(ji,0) = 0._wp 284 maxnbot0 = MAX ( maxnbot0 , nbot0(ji) ) 285 END DO 286 IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice ) 287 288 DO jk = 1, maxnbot0 112 ! new layers interfaces 113 zh_cum1(:,0:nlay_i) = 0._wp 114 DO jk1 = 1, nlay_i 289 115 DO ji = kideb, kiut 290 !change 291 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 292 limsum = MIN( limsum , nlay_s ) 293 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 294 END DO 295 END DO 296 297 DO ji = kideb, kiut 298 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 299 zm0(ji,1) = dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 300 END DO 301 302 DO jk = ntop0, maxnbot0 303 DO ji = kideb, kiut 304 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) ! layer thickness 305 END DO 306 END DO 307 308 zqts_in(:) = 0._wp 309 310 DO ji = kideb, kiut ! layer heat content 311 qm0 (ji,1) = rhosn * ( cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji) & 312 & - REAL( snswi(ji) ) * t_s_b (ji,1) ) & 313 & + lfus ) * zthick0(ji,1) 314 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) 315 END DO 316 317 DO jk = 2, maxnbot0 318 DO ji = kideb, kiut 319 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 320 limsum = MIN( limsum , nlay_s ) 321 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 322 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 323 zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 324 END DO ! jk 325 END DO ! ji 326 327 !------------------------------------------------ 328 ! Energy given by the snow in snow-ice formation 329 !------------------------------------------------ 330 ! zqsnow, enthalpy of the flooded snow 331 DO ji = kideb, kiut 332 zqsnow (ji) = rhosn * lfus 333 zdeltah(ji) = 0._wp 334 END DO 335 336 DO jk = nlays0, 1, -1 337 DO ji = kideb, kiut 338 zhsnow = MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) ) 339 zqsnow (ji) = zqsnow (ji) + rhosn*cpic*(rtt-t_s_b(ji,jk)) 340 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 341 END DO 342 END DO 343 344 DO ji = kideb, kiut 345 zqsnow(ji) = zqsnow(ji) * dh_snowice(ji) 346 END DO 347 348 !------------------ 349 ! new snow profile 350 !------------------ 351 352 !-------------- 353 ! Vector index 354 !-------------- 355 ntop1 = 1 356 nbot1 = nlay_s 357 358 !------------------- 359 ! Layer coordinates 360 !------------------- 361 DO ji = kideb, kiut 362 zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 363 z_s(ji,0) = 0._wp 116 zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 117 ENDDO 364 118 ENDDO 365 119 366 DO jk = 1, nlay_s 120 zqh_cum1(:,0:nlay_i) = 0._wp 121 ! new cumulative q*h => linear interpolation 122 DO jk0 = 1, nlay_i+1 123 DO jk1 = 1, nlay_i-1 124 DO ji = kideb, kiut 125 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 126 zqh_cum1(ji,jk1) = ( zqh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & 127 & zqh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) & 128 & / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 129 ENDIF 130 ENDDO 131 ENDDO 132 ENDDO 133 ! to ensure that total heat content is strictly conserved, set: 134 zqh_cum1(:,nlay_i) = zqh_cum0(:,nlay_i+2) 135 136 ! new enthalpies 137 DO jk1 = 1, nlay_i 367 138 DO ji = kideb, kiut 368 z_s(ji,jk) = zh_s(ji) * REAL( jk ) 369 END DO 370 END DO 371 372 !----------------- 373 ! Layer thickness 374 !----------------- 375 DO layer0 = ntop0, maxnbot0 376 DO ji = kideb, kiut 377 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) 378 END DO 379 END DO 380 381 DO layer1 = ntop1, nbot1 382 DO ji = kideb, kiut 383 q_s_b(ji,layer1) = 0._wp 384 END DO 385 END DO 386 387 !---------------- 388 ! Weight factors 389 !---------------- 390 DO layer0 = ntop0, maxnbot0 391 DO layer1 = ntop1, nbot1 392 DO ji = kideb, kiut 393 zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 394 zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & 395 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 396 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 397 & * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 398 END DO 399 END DO 400 END DO 401 402 ! Heat conservation 403 zqts_fin(:) = 0._wp 404 DO jk = 1, nlay_s 405 DO ji = kideb, kiut 406 zqts_fin(ji) = zqts_fin(ji) + q_s_b(ji,jk) 407 END DO 408 END DO 409 410 IF ( con_i .AND. jiindex_1d > 0 ) THEN 411 DO ji = kideb, kiut 412 IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 413 ii = MOD( npb(ji) - 1, jpi ) + 1 414 ij = ( npb(ji) - 1 ) / jpi + 1 415 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 416 WRITE(numout,*) ' ji, jj : ', ii, ij 417 WRITE(numout,*) ' ht_s_b : ', ht_s_b(ji) 418 WRITE(numout,*) ' zqts_in : ', zqts_in (ji) * r1_rdtice 419 WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice 420 WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji) 421 WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji) 422 WRITE(numout,*) ' snswi : ', snswi(ji) 423 ENDIF 424 END DO 425 ENDIF 426 427 !--------------------- 428 ! Recover heat content 429 !--------------------- 430 DO jk = 1, nlay_s 431 DO ji = kideb, kiut 432 zinda = MAX( 0._wp, SIGN( 1._wp , zh_s(ji) - epsi10 ) ) 433 q_s_b(ji,jk) = zinda * q_s_b(ji,jk) / MAX( zh_s(ji) , epsi10 ) 434 END DO !ji 435 END DO !jk 436 437 !--------------------- 438 ! Recover temperature 439 !--------------------- 440 zfac1 = 1. / ( rhosn * cpic ) 441 zfac2 = lfus / cpic 442 DO jk = 1, nlay_s 443 DO ji = kideb, kiut 444 zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 445 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 446 END DO 447 END DO 448 ! 449 !------------------------------------------------------------------------------| 450 ! 4) Ice redistribution | 451 !------------------------------------------------------------------------------| 452 ! 453 !------------- 454 ! OLD PROFILE 455 !------------- 456 457 !---------------- 458 ! Vector indexes 459 !---------------- 460 ntop0 = 1 461 maxnbot0 = 0 462 463 DO ji = kideb, kiut 464 ! reference number of the bottommost layer 465 nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + & 466 & ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) ) 467 ! maximum reference number of the bottommost layer over all domain 468 maxnbot0 = MAX( maxnbot0 , nbot0(ji) ) 469 END DO 470 471 !------------------------- 472 ! Cotes of old ice layers 473 !------------------------- 474 zm0(:,0) = 0._wp 475 476 DO jk = 1, maxnbot0 477 DO ji = kideb, kiut 478 ! jk goes from 1 to nbot0 479 ! the ice layer number goes from 1 to nlay_i 480 ! limsum is the real ice layer number corresponding to present jk 481 limsum = ( (icsuswi(ji)*(icsuind(ji)+jk-1) + & 482 (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 483 zm0(ji,jk)= REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 484 + REAL(limsum) * zh_i(ji) 485 END DO 486 END DO 487 488 DO ji = kideb, kiut 489 zm0(ji,nbot0(ji)) = REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 490 + zh_i(ji) * REAL(nlayi0) 491 zm0(ji,1) = REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 492 END DO 493 494 !----------------------------- 495 ! Thickness of old ice layers 496 !----------------------------- 497 DO jk = ntop0, maxnbot0 498 DO ji = kideb, kiut 499 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) 500 END DO 501 END DO 502 503 !--------------------------- 504 ! Inner layers heat content 505 !--------------------------- 506 qm0(:,:) = 0.0 507 zqti_in(:) = 0.0 508 509 DO jk = ntop0, maxnbot0 510 DO ji = kideb, kiut 511 limsum = MAX(1,MIN(snicswi(ji)*(jk-1) + icsuswi(ji)*(jk-1+icsuind(ji)) + & 512 (1-icsuswi(ji))*(1-snicswi(ji))*jk,nlay_i)) 513 ztmelts = -tmut * s_i_b(ji,limsum) + rtt 514 qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ & 515 MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) & 516 * zthick0(ji,jk) 517 END DO 518 END DO 519 520 !---------------------------- 521 ! Bottom layers heat content 522 !---------------------------- 523 DO ji = kideb, kiut 524 ztmelts = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice 525 & + REAL( icboswi(ji) ) * (-tmut * s_i_new(ji) ) & ! case of forming ice 526 & + rtt ! in Kelvin 527 528 ! bottom formation temperature 529 ztform = t_i_b(ji,nlay_i) 530 IF( num_sal == 2 ) ztform = t_bo_b(ji) 531 qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 532 & + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice 533 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 534 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) 535 END DO 536 537 !----------------------------- 538 ! Snow ice layer heat content 539 !----------------------------- 540 DO ji = kideb, kiut 541 ! energy of the flooding seawater 542 zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 543 (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 544 ! Heat conservation diagnostic 545 qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic 546 547 qldif_1d(ji) = qldif_1d(ji) + zqsnic * a_i_b(ji) 548 549 ! enthalpy of the newly formed snow-ice layer 550 ! = enthalpy of snow + enthalpy of frozen water 551 zqsnic = zqsnow(ji) + zqsnic 552 qm0(ji,1) = REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 553 554 END DO ! ji 555 556 DO jk = ntop0, maxnbot0 557 DO ji = kideb, kiut 558 ! Heat conservation 559 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi10) ) & 560 & * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 561 END DO 562 END DO 563 564 !------------- 565 ! NEW PROFILE 566 !------------- 567 568 !--------------- 569 ! Vectors index 570 !--------------- 571 ntop1 = 1 572 nbot1 = nlay_i 573 574 !------------------ 575 ! Layers thickness 576 !------------------ 577 DO ji = kideb, kiut 578 zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 139 zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 140 qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 ) 141 ENDDO 579 142 ENDDO 580 143 581 !------------- 582 ! Layer cotes 583 !------------- 584 z_i(:,0) = 0._wp 585 DO jk = 1, nlay_i 586 DO ji = kideb, kiut 587 z_i(ji,jk) = zh_i(ji) * jk 588 END DO 144 ! --- diag error on heat remapping --- ! 145 ! comment: if input h_i_old and qh_i_old are already multiplied by a_i (as in limthd_lac), 146 ! then we should not (* a_i) again but not important since this is just to check that remap error is ~0 147 DO ji = kideb, kiut 148 hfx_err_rem_1d(ji) = hfx_err_rem_1d(ji) + a_i_1d(ji) * r1_rdtice * & 149 & ( SUM( qnew(ji,1:nlay_i) ) * zhnew(ji) - SUM( qh_i_old(ji,0:nlay_i+1) ) ) 589 150 END DO 590 591 !--thicknesses of the layers 592 DO layer0 = ntop0, maxnbot0 593 DO ji = kideb, kiut 594 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) ! thicknesses of the layers 595 END DO 596 END DO 597 598 !------------------------ 599 ! Weights for relayering 600 !------------------------ 601 q_i_b(:,:) = 0._wp 602 DO layer0 = ntop0, maxnbot0 603 DO layer1 = ntop1, nbot1 604 DO ji = kideb, kiut 605 zinda = MAX( 0._wp, SIGN( 1._wp , zhl0(ji,layer0) - epsi10 ) ) 606 zrl01(layer1,layer0) = zinda * MAX(0.0,( MIN(zm0(ji,layer0),z_i(ji,layer1)) & 607 - MAX(zm0(ji,layer0-1), z_i(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10)) 608 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 609 + zrl01(layer1,layer0)*qm0(ji,layer0) & 610 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi10)) & 611 * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 612 END DO 613 END DO 614 END DO 615 616 !------------------------- 617 ! Heat conservation check 618 !------------------------- 619 zqti_fin(:) = 0._wp 620 DO jk = 1, nlay_i 621 DO ji = kideb, kiut 622 zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) 623 END DO 624 END DO 151 625 152 ! 626 IF ( con_i .AND. jiindex_1d > 0 ) THEN 627 DO ji = kideb, kiut 628 IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice > 1.0e-6 ) THEN 629 ii = MOD( npb(ji) - 1, jpi ) + 1 630 ij = ( npb(ji) - 1 ) / jpi + 1 631 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 632 WRITE(numout,*) ' ji, jj : ', ii, ij 633 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) 634 WRITE(numout,*) ' zqti_in : ', zqti_in (ji) * r1_rdtice 635 WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 636 WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 637 WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 638 WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 639 WRITE(numout,*) ' icsuswi : ', icsuswi(ji) 640 WRITE(numout,*) ' icboswi : ', icboswi(ji) 641 WRITE(numout,*) ' snicswi : ', snicswi(ji) 642 ENDIF 643 END DO 644 ENDIF 645 646 !---------------------- 647 ! Recover heat content 648 !---------------------- 649 DO jk = 1, nlay_i 650 DO ji = kideb, kiut 651 zinda = MAX( 0._wp, SIGN( 1._wp , zh_i(ji) - epsi10 ) ) 652 q_i_b(ji,jk) = zinda * q_i_b(ji,jk) / MAX( zh_i(ji) , epsi10 ) 653 END DO !ji 654 END DO !jk 655 656 ! Heat conservation 657 zqti_fin(:) = 0.0 658 DO jk = 1, nlay_i 659 DO ji = kideb, kiut 660 zqti_fin(ji) = zqti_fin(ji) + q_i_b(ji,jk) * zh_i(ji) 661 END DO 662 END DO 663 664 ! 665 !------------------------------------------------------------------------------| 666 ! 5) Update salinity and recover temperature | 667 !------------------------------------------------------------------------------| 668 ! 669 ! Update salinity (basal entrapment, snow ice formation) 670 DO ji = kideb, kiut 671 sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 672 END DO !ji 673 674 ! Recover temperature 675 DO jk = 1, nlay_i 676 DO ji = kideb, kiut 677 ztmelts = -tmut*s_i_b(ji,jk) + rtt 678 !Conversion q(S,T) -> T (second order equation) 679 zaaa = cpic 680 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 681 zccc = lfus * ( ztmelts - rtt ) 682 zdiscrim = SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) ) 683 t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) 684 END DO !ji 685 686 END DO !jk 687 ! 688 CALL wrk_dealloc( jpij, snswi, nbot0, icsuind, icsuswi, icboind, icboswi, snicind, snicswi, snind ) ! integer 689 CALL wrk_dealloc( jpij, zh_i, zh_s, zqsnow, zdeltah, zqti_in, zqts_in, zqti_fin, zqts_fin ) ! real 690 CALL wrk_dealloc( jpij,jkmax+4, zm0, qm0, z_s, z_i, zthick0, zhl0, kjstart = 0 ) 691 CALL wrk_dealloc( jkmax+4,jkmax+4, zrl01, kistart = 0, kjstart = 0 ) 153 CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 154 CALL wrk_dealloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 155 CALL wrk_dealloc( jpij, zhnew ) 692 156 ! 693 157 END SUBROUTINE lim_thd_ent
Note: See TracChangeset
for help on using the changeset viewer.