Changeset 2715 for trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
- Timestamp:
- 2011-03-30T17:58:35+02:00 (13 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90 (modified) (39 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
r2528 r2715 6 6 !! after vertical growth/decay 7 7 !!====================================================================== 8 !! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D 9 !! ! 2005-07 (M. Vancoppenolle) 3D version 10 !! ! 2006-11 (X. Fettweis) Vectorized 11 !! 3.0 ! 2008-03 (M. Vancoppenolle) Energy conservation and clean code 12 !! 4.0 ! 2011-02 (G. Madec) dynamical allocation 13 !!---------------------------------------------------------------------- 8 14 #if defined key_lim3 9 15 !!---------------------------------------------------------------------- … … 13 19 !!---------------------------------------------------------------------- 14 20 USE par_oce ! ocean parameters 15 USE dom_oce 16 USE domain 17 USE in_out_manager 18 USE phycst 19 USE thd_ice 20 USE ice 21 USE limvar 22 USE par_ice 23 USE lib_mpp 21 USE dom_oce ! domain variables 22 USE domain ! 23 USE phycst ! physical constants 24 USE ice ! LIM variables 25 USE par_ice ! LIM parameters 26 USE thd_ice ! LIM thermodynamics 27 USE limvar ! LIM variables 28 USE in_out_manager ! I/O manager 29 USE wrk_nemo ! workspace manager 30 USE lib_mpp ! MPP library 24 31 25 32 IMPLICIT NONE … … 28 35 PUBLIC lim_thd_ent ! called by lim_thd 29 36 30 REAL(wp) :: & ! constant values 31 epsi20 = 1.e-20 , & 32 epsi13 = 1.e-13 , & 33 zzero = 0.e0 , & 34 zone = 1.e0 , & 35 epsi10 = 1.0e-10 37 REAL(wp) :: epsi20 = 1e-20_wp ! constant values 38 REAL(wp) :: epsi13 = 1e-13_wp ! 39 REAL(wp) :: epsi10 = 1e-10_wp ! 40 REAL(wp) :: epsi06 = 1e-06_wp ! 41 REAL(wp) :: zzero = 0._wp ! 42 REAL(wp) :: zone = 1._wp ! 43 36 44 !!---------------------------------------------------------------------- 37 !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)45 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 38 46 !! $Id$ 39 47 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 41 49 CONTAINS 42 50 43 SUBROUTINE lim_thd_ent( kideb,kiut,jl)51 SUBROUTINE lim_thd_ent( kideb, kiut, jl ) 44 52 !!------------------------------------------------------------------- 45 53 !! *** ROUTINE lim_thd_ent *** … … 60 68 !! 5) Ice salinity, recover temperature 61 69 !! 62 !! ** Arguments 63 !! 64 !! ** Inputs / Outputs 65 !! 66 !! ** External 67 !! 68 !! ** References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 69 !! 70 !! ** History : (05-2003) Martin V. UCL-Astr 71 !! (07-2005) Martin for 3d adapatation 72 !! (11-2006) Vectorized by Xavier Fettweis (ASTR) 73 !! (03-2008) Energy conservation and clean code 74 !! * Arguments 75 76 INTEGER , INTENT(IN):: & 77 kideb , & ! start point on which the the computation is applied 78 kiut , & ! end point on which the the computation is applied 79 jl ! thickness category number 80 81 INTEGER :: & 82 ji,jk , & ! dummy loop indices 83 zji, zjj , & ! dummy indices 70 !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 71 !!------------------------------------------------------------------- 72 INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied 73 INTEGER , INTENT(in) :: jl ! Thickness cateogry number 74 75 INTEGER :: ji,jk ! dummy loop indices 76 INTEGER :: zji, zjj , & ! dummy indices 84 77 ntop0 , & ! old layer top index 85 78 nbot1 , & ! new layer bottom index … … 90 83 layer0, layer1 ! old/new layer indexes 91 84 92 INTEGER, DIMENSION(jpij) :: &93 snswi , & ! snow switch94 nbot0 , & ! old layer bottom index95 icsuind , & ! ice surface index96 icsuswi , & ! ice surface switch97 icboind , & ! ice bottom index98 icboswi , & ! ice bottom switch99 snicind , & ! snow ice index100 snicswi , & ! snow ice switch101 snind ! snow index102 85 103 86 REAL(wp) :: & 104 zeps, zeps6 , & ! numerical constant very small105 87 ztmelts , & ! ice melting point 106 88 zqsnic , & ! enthalpy of snow ice layer … … 115 97 zdiscrim !: dummy factor 116 98 117 REAL(wp), DIMENSION(jpij) :: & 118 zh_i , & ! thickness of an ice layer 119 zh_s , & ! thickness of a snow layer 120 zqsnow , & ! enthalpy of the snow put in snow ice 121 zdeltah ! temporary variable 122 123 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: & 124 zm0 , & ! old layer-system vertical cotes 125 qm0 , & ! old layer-system heat content 126 z_s , & ! new snow system vertical cotes 127 z_i , & ! new ice system vertical cotes 128 zthick0 ! old ice thickness 129 130 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: & 131 zhl0 ! old and new layer thicknesses 132 133 REAL(wp), DIMENSION(0:jkmax+3,0:jkmax+3) :: & 134 zrl01 135 136 ! Energy conservation 137 REAL(wp), DIMENSION(jpij) :: & 138 zqti_in, zqts_in, & 139 zqti_fin, zqts_fin 140 141 !------------------------------------------------------------------------------| 142 143 zeps = 1.0d-20 144 zeps6 = 1.0d-06 145 zthick0(:,:) = 0.0 146 zm0(:,:) = 0.0 147 qm0(:,:) = 0.0 148 zrl01(:,:) = 0.0 149 zhl0(:,:) = 0.0 150 z_i(:,:) = 0.0 151 z_s(:,:) = 0.0 99 INTEGER, DIMENSION(jpij) :: & 100 snswi , & ! snow switch 101 nbot0 , & ! old layer bottom index 102 icsuind , & ! ice surface index 103 icsuswi , & ! ice surface switch 104 icboind , & ! ice bottom index 105 icboswi , & ! ice bottom switch 106 snicind , & ! snow ice index 107 snicswi , & ! snow ice switch 108 snind ! snow index 109 ! 110 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: zm0 ! old layer-system vertical cotes 111 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: qm0 ! old layer-system heat content 112 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: z_s ! new snow system vertical cotes 113 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: z_i ! new ice system vertical cotes 114 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: zthick0 ! old ice thickness 115 REAL(wp), DIMENSION(jpij,0:jkmax+3) :: zhl0 ! old and new layer thicknesses 116 ! 117 REAL(wp), DIMENSION(0:jkmax+3,0:jkmax+3) :: zrl01 118 ! 119 REAL(wp), POINTER, DIMENSION(:) :: zh_i, zqsnow , zqti_in, zqti_fin 120 REAL(wp), POINTER, DIMENSION(:) :: zh_s, zdeltah, zqts_in, zqts_fin 121 !!------------------------------------------------------------------- 122 123 IF( wrk_in_use(1, 1,2,3,4,5,6,7,8) ) THEN 124 CALL ctl_stop('lim_thd_dh : requestead workspace arrays unavailable') ; RETURN 125 END IF 126 127 ! Set-up pointers to sub-arrays of workspace arrays 128 zh_i => wrk_1d_1 (1:jpij) ! thickness of an ice layer 129 zh_s => wrk_1d_2 (1:jpij) ! thickness of a snow layer 130 zqsnow => wrk_1d_3 (1:jpij) ! enthalpy of the snow put in snow ice 131 zdeltah => wrk_1d_4 (1:jpij) ! temporary variable 132 zqti_in => wrk_1d_5 (1:jpij) ! Energy conservation 133 zqts_in => wrk_1d_6 (1:jpij) ! - - 134 zqti_fin => wrk_1d_7 (1:jpij) ! - - 135 zqts_fin => wrk_1d_8 (1:jpij) ! - - 136 137 zthick0(:,:) = 0._wp 138 zm0 (:,:) = 0._wp 139 qm0 (:,:) = 0._wp 140 zrl01 (:,:) = 0._wp 141 zhl0 (:,:) = 0._wp 142 z_i (:,:) = 0._wp 143 z_s (:,:) = 0._wp 152 144 153 145 ! … … 155 147 ! 1) Grid | 156 148 !------------------------------------------------------------------------------| 157 ! 158 nlays0 = nlay_s 159 nlayi0 = nlay_i 160 161 DO ji = kideb, kiut 162 zh_i(ji) = old_ht_i_b(ji) / nlay_i 163 zh_s(ji) = old_ht_s_b(ji) / nlay_s 164 ENDDO 149 nlays0 = nlay_s 150 nlayi0 = nlay_i 151 152 DO ji = kideb, kiut 153 zh_i(ji) = old_ht_i_b(ji) / nlay_i 154 zh_s(ji) = old_ht_s_b(ji) / nlay_s 155 END DO 165 156 166 157 ! … … 168 159 ! 2) Switches | 169 160 !------------------------------------------------------------------------------| 170 !171 161 ! 2.1 snind(ji), snswi(ji) 172 162 ! snow surface behaviour : computation of snind(ji)-snswi(ji) … … 176 166 ! 2 if 2nd layer is melting ... 177 167 DO ji = kideb, kiut 178 snind (ji)= 0179 zdeltah(ji) = 0.0168 snind (ji) = 0 169 zdeltah(ji) = 0._wp 180 170 ENDDO !ji 181 171 182 172 DO jk = 1, nlays0 183 173 DO ji = kideb, kiut 184 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)- zeps))) &185 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)- zeps))))174 snind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) & 175 + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20)))) 186 176 zdeltah(ji)= zdeltah(ji) + zh_s(ji) 187 177 END DO ! ji 188 END DO ! jk178 END DO ! jk 189 179 190 180 ! snswi(ji) : switch which value equals 1 if snow melts 191 181 ! 0 if not 192 182 DO ji = kideb, kiut 193 snswi(ji) = MAX(0,INT(-dh_s_tot(ji)/MAX( zeps,ABS(dh_s_tot(ji)))))194 END DO ! ji183 snswi(ji) = MAX(0,INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 184 END DO ! ji 195 185 196 186 ! 2.2 icsuind(ji), icsuswi(ji) … … 201 191 ! 2 if 2nd layer is reached by melt ... 202 192 DO ji = kideb, kiut 203 icsuind(ji) 204 zdeltah(ji) = 0.0205 END DO !ji193 icsuind(ji) = 0 194 zdeltah(ji) = 0._wp 195 END DO !ji 206 196 DO jk = 1, nlayi0 207 197 DO ji = kideb, kiut 208 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)- zeps))) &209 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)- zeps))))198 icsuind(ji) = jk * INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) & 199 + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20)))) 210 200 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 211 201 END DO ! ji … … 216 206 ! 0 if not 217 207 DO ji = kideb, kiut 218 icsuswi(ji) = MAX(0,INT(-dh_i_surf(ji)/MAX( zeps, ABS(dh_i_surf(ji)) ) ) )208 icsuswi(ji) = MAX(0,INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 219 209 ENDDO 220 210 … … 227 217 ! N+1 if all layers melt and that snow transforms into ice 228 218 DO ji = kideb, kiut 229 icboind(ji) 230 zdeltah(ji) = 0.0231 END DO219 icboind(ji) = 0 220 zdeltah(ji) = 0._wp 221 END DO 232 222 DO jk = nlayi0, 1, -1 233 223 DO ji = kideb, kiut 234 icboind(ji) = (nlayi0+1-jk) & 235 * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps))) & 236 + icboind(ji) & 237 * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-zeps)))) 224 icboind(ji) = (nlayi0+1-jk) * INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) & 225 & + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20)))) 238 226 zdeltah(ji) = zdeltah(ji) + zh_i(ji) 239 227 END DO 240 END DO228 END DO 241 229 242 230 DO ji = kideb, kiut 243 231 ! case of total ablation with remaining snow 244 IF ( ( ht_i_b(ji) .GT. zeps) .AND. &245 ( ht_i_b(ji) - dh_snowice(ji) .LT. zeps) ) icboind(ji) = nlay_i + 1232 IF ( ( ht_i_b(ji) .GT. epsi20 ) .AND. & 233 ( ht_i_b(ji) - dh_snowice(ji) .LT. epsi20 ) ) icboind(ji) = nlay_i + 1 246 234 END DO 247 235 … … 250 238 ! 0 if ablation is on the way 251 239 DO ji = kideb, kiut 252 icboswi(ji) = MAX(0,INT(dh_i_bott(ji) / MAX(zeps,ABS(dh_i_bott(ji)))))253 END DO240 icboswi(ji) = MAX(0,INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 241 END DO 254 242 255 243 ! 2.4 snicind(ji), snicswi(ji) … … 260 248 ! 2 if penultiem layer ... 261 249 DO ji = kideb, kiut 262 snicind(ji) 263 zdeltah(ji) = 0.0264 END DO250 snicind(ji) = 0 251 zdeltah(ji) = 0._wp 252 END DO 265 253 DO jk = nlays0, 1, -1 266 254 DO ji = kideb, kiut 267 255 snicind(ji) = (nlays0+1-jk) & 268 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps))) & 269 + snicind(ji) & 270 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-zeps)))) 256 * INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji) & 257 * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20)))) 271 258 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 272 259 END DO 273 END DO260 END DO 274 261 275 262 ! snicswi(ji) : switch which equals … … 277 264 ! 0 if not 278 265 DO ji = kideb, kiut 279 snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX( zeps,ABS(dh_snowice(ji)))))266 snicswi(ji) = MAX(0,INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 280 267 ENDDO 281 268 … … 294 281 ! indexes of the vectors 295 282 !------------------------ 296 ntop0 = 1 297 maxnbot0 = 0 298 299 DO ji = kideb, kiut 300 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1. - snicind(ji) ) * & 301 snicswi(ji) 283 ntop0 = 1 284 maxnbot0 = 0 285 286 DO ji = kideb, kiut 287 nbot0(ji) = nlays0 + 1 - snind(ji) + ( 1. - snicind(ji) ) * snicswi(ji) 302 288 ! cotes of the top of the layers 303 zm0(ji,0) = 0.0304 maxnbot0 305 END DO306 IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice )289 zm0(ji,0) = 0._wp 290 maxnbot0 = MAX ( maxnbot0 , nbot0(ji) ) 291 END DO 292 IF( lk_mpp ) CALL mpp_max( maxnbot0, kcom=ncomm_ice ) 307 293 308 294 DO jk = 1, maxnbot0 309 295 DO ji = kideb, kiut 310 296 !change 311 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + & 312 snswi(ji) * ( jk + snind(ji) - 1 ) 297 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 298 limsum = MIN( limsum , nlay_s ) 299 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum 300 END DO 301 END DO 302 303 DO ji = kideb, kiut 304 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0 305 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1) 306 END DO 307 308 DO jk = ntop0, maxnbot0 309 DO ji = kideb, kiut 310 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) ! layer thickness 311 END DO 312 END DO 313 314 zqts_in(:) = 0._wp 315 316 DO ji = kideb, kiut ! layer heat content 317 qm0 (ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * tatm_ice_1d(ji) & 318 & - snswi(ji) * t_s_b (ji,1) ) & 319 & + lfus ) * zthick0(ji,1) 320 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) 321 END DO 322 323 DO jk = 2, maxnbot0 324 DO ji = kideb, kiut 325 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 313 326 limsum = MIN( limsum , nlay_s ) 314 zm0(ji,jk) = dh_s_tot(ji) + zh_s(ji) * limsum 315 END DO 316 ENDDO 317 318 DO ji = kideb, kiut 319 zm0(ji,nbot0(ji)) = dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + & 320 zh_s(ji) * nlays0 321 zm0(ji,1) = dh_s_tot(ji) * (1 -snswi(ji) ) + & 322 snswi(ji) * zm0(ji,1) 323 ENDDO 324 325 DO jk = ntop0, maxnbot0 326 DO ji = kideb, kiut 327 ! layer thickness 328 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) 329 END DO 330 ENDDO 331 332 zqts_in(:) = 0.0 333 334 DO ji = kideb, kiut 335 ! layer heat content 336 qm0(ji,1) = rhosn * ( cpic * ( rtt - ( 1. - snswi(ji) ) * ( tatm_ice_1d(ji) ) & 337 - snswi(ji) * t_s_b(ji,1) ) & 338 + lfus ) * zthick0(ji,1) 339 zqts_in(ji) = zqts_in(ji) + qm0(ji,1) 340 ENDDO 341 342 DO jk = 2, maxnbot0 343 DO ji = kideb, kiut 344 limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + & 345 snswi(ji) * ( jk + snind(ji) - 1 ) 346 limsum = MIN( limsum , nlay_s ) 347 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) & 348 * zthick0(ji,jk) 349 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, zeps - ht_s_b(ji) ) ) 327 qm0(ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 328 zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20 - ht_s_b(ji) ) ) 350 329 zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch 351 330 END DO ! jk 352 END DO ! ji331 END DO ! ji 353 332 354 333 !------------------------------------------------ … … 357 336 ! zqsnow, enthalpy of the flooded snow 358 337 DO ji = kideb, kiut 359 zqsnow (ji) = rhosn*lfus360 zdeltah(ji) = 0.0361 END DO338 zqsnow (ji) = rhosn * lfus 339 zdeltah(ji) = 0._wp 340 END DO 362 341 363 342 DO jk = nlays0, 1, -1 364 343 DO ji = kideb, kiut 365 zhsnow = MAX(0.0,dh_snowice(ji)-zdeltah(ji)) 366 zqsnow(ji) = zqsnow(ji) + & 367 rhosn*cpic*(rtt-t_s_b(ji,jk)) 344 zhsnow = MAX( 0._wp , dh_snowice(ji)-zdeltah(ji) ) 345 zqsnow (ji) = zqsnow (ji) + rhosn*cpic*(rtt-t_s_b(ji,jk)) 368 346 zdeltah(ji) = zdeltah(ji) + zh_s(ji) 369 347 END DO 370 END DO348 END DO 371 349 372 350 DO ji = kideb, kiut … … 381 359 ! Vector index 382 360 !-------------- 383 ntop1 384 nbot1 361 ntop1 = 1 362 nbot1 = nlay_s 385 363 386 364 !------------------- … … 389 367 DO ji = kideb, kiut 390 368 zh_s(ji) = ht_s_b(ji) / nlay_s 391 z_s(ji,0) = 0. 0369 z_s(ji,0) = 0._wp 392 370 ENDDO 393 371 … … 396 374 z_s(ji,jk) = zh_s(ji) * jk 397 375 END DO 398 END DO376 END DO 399 377 400 378 !----------------- … … 405 383 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) 406 384 END DO 407 END DO385 END DO 408 386 409 387 DO layer1 = ntop1, nbot1 410 388 DO ji = kideb, kiut 411 q_s_b(ji,layer1) = 0.0412 END DO 413 END DO389 q_s_b(ji,layer1) = 0._wp 390 END DO 391 END DO 414 392 415 393 !---------------- … … 419 397 DO layer1 = ntop1, nbot1 420 398 DO ji = kideb, kiut 421 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) &422 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1)))/MAX(zhl0(ji,layer0),epsi10))423 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) &424 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+zeps))399 zrl01(layer1,layer0) = MAX(0.0,( MIN(zm0(ji,layer0),z_s(ji,layer1)) & 400 & - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10)) 401 q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0) & 402 & * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20)) 425 403 END DO 426 404 END DO 427 END DO405 END DO 428 406 429 407 ! Heat conservation 430 zqts_fin(:) = 0. 0408 zqts_fin(:) = 0._wp 431 409 DO jk = 1, nlay_s 432 410 DO ji = kideb, kiut … … 458 436 DO jk = 1, nlay_s 459 437 DO ji = kideb, kiut 460 q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , zeps)438 q_s_b(ji,jk) = q_s_b(ji,jk) / MAX( zh_s(ji) , epsi20 ) 461 439 END DO !ji 462 END DO !jk440 END DO !jk 463 441 464 442 !--------------------- … … 469 447 DO jk = 1, nlay_s 470 448 DO ji = kideb, kiut 471 zswitch = MAX ( 0.0 , SIGN ( 1.0, zeps - ht_s_b(ji) ) ) 472 t_s_b(ji,jk) = rtt & 473 + ( 1.0 - zswitch ) * & 474 ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 475 END DO 476 ENDDO 449 zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20 - ht_s_b(ji) ) ) 450 t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 451 END DO 452 END DO 477 453 ! 478 454 !------------------------------------------------------------------------------| … … 487 463 ! Vector indexes 488 464 !---------------- 489 ntop0 490 maxnbot0 465 ntop0 = 1 466 maxnbot0 = 0 491 467 492 468 DO ji = kideb, kiut 493 469 ! reference number of the bottommost layer 494 nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + & 495 ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , & 496 nlay_i + 2 ) ) 470 nbot0(ji) = MAX( 1 , MIN( nlayi0 + ( 1 - icboind(ji) ) + & 471 & ( 1 - icsuind(ji) ) * icsuswi(ji) + snicswi(ji) , nlay_i + 2 ) ) 497 472 ! maximum reference number of the bottommost layer over all domain 498 maxnbot0 499 END DO473 maxnbot0 = MAX( maxnbot0 , nbot0(ji) ) 474 END DO 500 475 501 476 !------------------------- 502 477 ! Cotes of old ice layers 503 478 !------------------------- 504 zm0(:,0) = 0.0479 zm0(:,0) = 0.-wp 505 480 506 481 DO jk = 1, maxnbot0 … … 514 489 + limsum * zh_i(ji) 515 490 END DO 516 END DO491 END DO 517 492 518 493 DO ji = kideb, kiut … … 520 495 + zh_i(ji) * nlayi0 521 496 zm0(ji,1) = snicswi(ji)*dh_snowice(ji) + (1-snicswi(ji))*zm0(ji,1) 522 END DO497 END DO 523 498 524 499 !----------------------------- … … 529 504 zthick0(ji,jk) = zm0(ji,jk) - zm0(ji,jk-1) 530 505 END DO 531 END DO506 END DO 532 507 533 508 !--------------------------- … … 543 518 ztmelts = -tmut * s_i_b(ji,limsum) + rtt 544 519 qm0(ji,jk) = rhoic * ( cpic * (ztmelts-t_i_b(ji,limsum)) + lfus * ( 1.0-(ztmelts-rtt)/ & 545 MIN((t_i_b(ji,limsum)-rtt),- zeps) ) - rcp*(ztmelts-rtt) ) &520 MIN((t_i_b(ji,limsum)-rtt),-epsi20) ) - rcp*(ztmelts-rtt) ) & 546 521 * zthick0(ji,jk) 547 522 END DO 548 END DO523 END DO 549 524 550 525 !---------------------------- … … 552 527 !---------------------------- 553 528 DO ji = kideb, kiut 554 ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0))& ! case of melting ice555 + icboswi(ji) * (-tmut * s_i_new(ji))& ! case of forming ice556 + rtt ! this temperature is in Celsius529 ztmelts = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b (ji,nlayi0) ) & ! case of melting ice 530 & + icboswi(ji) * (-tmut * s_i_new(ji) ) & ! case of forming ice 531 & + rtt ! in Kelvin 557 532 558 533 ! bottom formation temperature 559 534 ztform = t_i_b(ji,nlay_i) 560 535 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji) 561 qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 562 + icboswi(ji) * & ! case of forming ice 563 rhoic*( cpic*(ztmelts-ztform) & 564 + lfus *( 1.0-(ztmelts-rtt)/ & 565 MIN ( (ztform-rtt) , - epsi10 ) ) & 566 - rcp*(ztmelts-rtt) ) & 567 *zthick0(ji,nbot0(ji)) 568 ENDDO 536 qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji)) & ! case of melting ice 537 & + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform) & ! case of forming ice 538 + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) ) & 539 - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji) ) 540 END DO 569 541 570 542 !----------------------------- … … 585 557 qm0(ji,1) = snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1) 586 558 587 END DO ! ji559 END DO ! ji 588 560 589 561 DO jk = ntop0, maxnbot0 590 562 DO ji = kideb, kiut 591 563 ! Heat conservation 592 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) & 593 * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-zeps6+zeps) ) & 594 * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + zeps ) ) 595 END DO 596 ENDDO 564 zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06+epsi20) ) & 565 & * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20 ) ) 566 END DO 567 END DO 597 568 598 569 !------------- … … 603 574 ! Vectors index 604 575 !--------------- 605 606 ntop1 = 1 607 nbot1 = nlay_i 576 ntop1 = 1 577 nbot1 = nlay_i 608 578 609 579 !------------------ … … 611 581 !------------------ 612 582 DO ji = kideb, kiut 613 zh_i(ji) 583 zh_i(ji) = ht_i_b(ji) / nlay_i 614 584 ENDDO 615 585 … … 617 587 ! Layer cotes 618 588 !------------- 619 z_i(:,0) = 0. 0589 z_i(:,0) = 0._wp 620 590 DO jk = 1, nlay_i 621 591 DO ji = kideb, kiut 622 592 z_i(ji,jk) = zh_i(ji) * jk 623 593 END DO 624 END DO594 END DO 625 595 626 596 !--thicknesses of the layers 627 597 DO layer0 = ntop0, maxnbot0 628 598 DO ji = kideb, kiut 629 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) !thicknesses of the layers630 END DO 631 END DO599 zhl0(ji,layer0) = zm0(ji,layer0) - zm0(ji,layer0-1) ! thicknesses of the layers 600 END DO 601 END DO 632 602 633 603 !------------------------ 634 604 ! Weights for relayering 635 605 !------------------------ 636 637 q_i_b(:,:) = 0.0 606 q_i_b(:,:) = 0._wp 638 607 DO layer0 = ntop0, maxnbot0 639 608 DO layer1 = ntop1, nbot1 … … 643 612 q_i_b(ji,layer1) = q_i_b(ji,layer1) & 644 613 + zrl01(layer1,layer0)*qm0(ji,layer0) & 645 * MAX(0.0,SIGN(1.0,ht_i_b(ji)- zeps6+zeps)) &646 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+ zeps))614 * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06+epsi20)) & 615 * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20)) 647 616 END DO 648 617 END DO 649 END DO618 END DO 650 619 651 620 !------------------------- 652 621 ! Heat conservation check 653 622 !------------------------- 654 zqti_fin(:) = 0. 0623 zqti_fin(:) = 0._wp 655 624 DO jk = 1, nlay_i 656 625 DO ji = kideb, kiut … … 663 632 zji = MOD( npb(ji) - 1, jpi ) + 1 664 633 zjj = ( npb(ji) - 1 ) / jpi + 1 665 WRITE(numout,*) ' violation of heat conservation : ', & 666 ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice 634 WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice 667 635 WRITE(numout,*) ' ji, jj : ', zji, zjj 668 636 WRITE(numout,*) ' ht_i_b : ', ht_i_b(ji) … … 683 651 DO jk = 1, nlay_i 684 652 DO ji = kideb, kiut 685 q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , zeps)653 q_i_b(ji,jk) = q_i_b(ji,jk) / MAX( zh_i(ji) , epsi20 ) 686 654 END DO !ji 687 END DO !jk655 END DO !jk 688 656 689 657 ! Heat conservation … … 702 670 ! Update salinity (basal entrapment, snow ice formation) 703 671 DO ji = kideb, kiut 704 sm_i_b(ji) = sm_i_b(ji) & 705 + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 672 sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 706 673 END DO !ji 707 674 708 675 ! Recover temperature 709 676 DO jk = 1, nlay_i 710 711 DO ji = kideb, kiut 712 677 DO ji = kideb, kiut 713 678 ztmelts = -tmut*s_i_b(ji,jk) + rtt 714 679 !Conversion q(S,T) -> T (second order equation) 715 680 zaaa = cpic 716 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + & 717 q_i_b(ji,jk) / rhoic - lfus 681 zbbb = ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 718 682 zccc = lfus * ( ztmelts - rtt ) 719 683 zdiscrim = SQRT( MAX(zbbb*zbbb - 4.0*zaaa*zccc,0.0) ) 720 t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / & 721 ( 2.0 *zaaa ) 684 t_i_b(ji,jk) = rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) 722 685 END DO !ji 723 686 724 687 END DO !jk 725 688 ! 689 IF( wrk_not_released(1, 1,2,3,4,5,6,7,8) ) CALL ctl_stop( 'lim_thd_ent : failed to release workspace arrays' ) 690 ! 726 691 END SUBROUTINE lim_thd_ent 727 692 728 693 #else 729 !!====================================================================== 730 !! *** MODULE limthd_ent *** 731 !! no sea ice model 732 !!====================================================================== 694 !!---------------------------------------------------------------------- 695 !! Default option NO LIM3 sea-ice model 696 !!---------------------------------------------------------------------- 733 697 CONTAINS 734 698 SUBROUTINE lim_thd_ent ! Empty routine 735 699 END SUBROUTINE lim_thd_ent 736 700 #endif 701 702 !!====================================================================== 737 703 END MODULE limthd_ent
Note: See TracChangeset
for help on using the changeset viewer.
