Changeset 921 for trunk/NEMO/LIM_SRC_3/limitd_me.F90
- Timestamp:
- 2008-05-13T10:28:52+02:00 (16 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/LIM_SRC_3/limitd_me.F90 (modified) (62 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/LIM_SRC_3/limitd_me.F90
r903 r921 32 32 USE prtctl ! Print control 33 33 USE lib_mpp 34 34 35 35 IMPLICIT NONE 36 36 PRIVATE … … 53 53 zone = 1.e0 54 54 55 !-----------------------------------------------------------------------56 ! Variables shared among ridging subroutines57 !-----------------------------------------------------------------------58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 !81 !-----------------------------------------------------------------------82 ! Ridging diagnostic arrays for history files83 !-----------------------------------------------------------------------84 !85 86 87 88 89 90 55 !----------------------------------------------------------------------- 56 ! Variables shared among ridging subroutines 57 !----------------------------------------------------------------------- 58 REAL(wp), DIMENSION (jpi,jpj) :: & 59 asum , & ! sum of total ice and open water area 60 aksum ! ratio of area removed to area ridged 61 62 REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: & 63 athorn ! participation function; fraction of ridging/ 64 ! closing associated w/ category n 65 66 REAL(wp), DIMENSION(jpi,jpj,jpl) :: & 67 hrmin , & ! minimum ridge thickness 68 hrmax , & ! maximum ridge thickness 69 hraft , & ! thickness of rafted ice 70 krdg , & ! mean ridge thickness/thickness of ridging ice 71 aridge , & ! participating ice ridging 72 araft ! participating ice rafting 73 74 REAL(wp), PARAMETER :: & 75 krdgmin = 1.1, & ! min ridge thickness multiplier 76 kraft = 2.0 ! rafting multipliyer 77 78 REAL(wp) :: & 79 Cp 80 ! 81 !----------------------------------------------------------------------- 82 ! Ridging diagnostic arrays for history files 83 !----------------------------------------------------------------------- 84 ! 85 REAL (wp), DIMENSION(jpi,jpj) :: & 86 dardg1dt , & ! rate of fractional area loss by ridging ice (1/s) 87 dardg2dt , & ! rate of fractional area gain by new ridges (1/s) 88 dvirdgdt , & ! rate of ice volume ridged (m/s) 89 opening ! rate of opening due to divergence/shear (1/s) 90 91 91 92 92 !!---------------------------------------------------------------------- … … 97 97 CONTAINS 98 98 99 !!-----------------------------------------------------------------------------!100 !!-----------------------------------------------------------------------------!99 !!-----------------------------------------------------------------------------! 100 !!-----------------------------------------------------------------------------! 101 101 102 102 SUBROUTINE lim_itd_me ! (subroutine 1/6) 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 !!-- End of declarations204 !-----------------------------------------------------------------------------!103 !!---------------------------------------------------------------------! 104 !! *** ROUTINE lim_itd_me *** 105 !! ** Purpose : 106 !! This routine computes the mechanical redistribution 107 !! of ice thickness 108 !! 109 !! ** Method : a very simple method :-) 110 !! 111 !! ** Arguments : 112 !! kideb , kiut : Starting and ending points on which the 113 !! the computation is applied 114 !! 115 !! ** Inputs / Ouputs : (global commons) 116 !! 117 !! ** External : 118 !! 119 !! ** Steps : 120 !! 1) Thickness categories boundaries, ice / o.w. concentrations 121 !! Ridge preparation 122 !! 2) Dynamical inputs (closing rate, divu_adv, opning) 123 !! 3) Ridging iteration 124 !! 4) Ridging diagnostics 125 !! 5) Heat, salt and freshwater fluxes 126 !! 6) Compute increments of tate variables and come back to old values 127 !! 128 !! ** References : There are a lot of references and can be difficult / 129 !! boring to read 130 !! 131 !! Flato, G. M., and W. D. Hibler III, 1995: Ridging and strength 132 !! in modeling the thickness distribution of Arctic sea ice, 133 !! J. Geophys. Res., 100, 18,611-18,626. 134 !! 135 !! Hibler, W. D. III, 1980: Modeling a variable thickness sea ice 136 !! cover, Mon. Wea. Rev., 108, 1943-1973, 1980. 137 !! 138 !! Rothrock, D. A., 1975: The energetics of the plastic deformation of 139 !! pack ice by ridging, J. Geophys. Res., 80, 4514-4519. 140 !! 141 !! Thorndike, A. S., D. A. Rothrock, G. A. Maykut, and R. Colony, 142 !! 1975: The thickness distribution of sea ice, J. Geophys. Res., 143 !! 80, 4501-4513. 144 !! 145 !! Bitz et al., JGR 2001 146 !! 147 !! Amundrud and Melling, JGR 2005 148 !! 149 !! Babko et al., JGR 2002 150 !! 151 !! ** History : 152 !! This routine is based on CICE code 153 !! and authors William H. Lipscomb, 154 !! and Elizabeth C. Hunke, LANL 155 !! are gratefully acknowledged 156 !! 157 !! (02-2006) Martin Vancoppenolle, UCL-ASTR 158 !! 159 !!--------------------------------------------------------------------! 160 !! * Arguments 161 162 !! * Local variables 163 INTEGER :: ji, & ! spatial dummy loop index 164 jj, & ! spatial dummy loop index 165 jk, & ! vertical layering dummy loop index 166 jl, & ! ice category dummy loop index 167 niter, & ! iteration counter 168 nitermax = 20 ! max number of ridging iterations 169 170 REAL(wp) :: & ! constant values 171 zeps = 1.0e-10, & 172 epsi10 = 1.0e-10, & 173 epsi06 = 1.0e-6 174 175 REAL(wp), DIMENSION(jpi,jpj) :: & 176 closing_net, & ! net rate at which area is removed (1/s) 177 ! (ridging ice area - area of new ridges) / dt 178 divu_adv , & ! divu as implied by transport scheme (1/s) 179 opning , & ! rate of opening due to divergence/shear 180 closing_gross, & ! rate at which area removed, not counting 181 ! area of new ridges 182 msnow_mlt , & ! mass of snow added to ocean (kg m-2) 183 esnow_mlt ! energy needed to melt snow in ocean (J m-2) 184 185 REAL(wp) :: & 186 w1, & ! temporary variable 187 tmpfac, & ! factor by which opening/closing rates are cut 188 dti ! 1 / dt 189 190 LOGICAL :: & 191 asum_error ! flag for asum .ne. 1 192 193 INTEGER :: iterate_ridging ! if true, repeat the ridging 194 195 REAL(wp) :: & 196 big = 1.0e8 197 198 REAL (wp), DIMENSION(jpi,jpj) :: & ! 199 vt_i_init, vt_i_final ! ice volume summed over categories 200 201 CHARACTER (len = 15) :: fieldid 202 203 !!-- End of declarations 204 !-----------------------------------------------------------------------------! 205 205 206 206 IF( numit == nstart ) CALL lim_itd_me_init ! Initialization (first time-step only) … … 211 211 ENDIF 212 212 213 !-----------------------------------------------------------------------------!214 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons215 !-----------------------------------------------------------------------------!216 ! Set hi_max(ncat) to a big value to ensure that all ridged ice217 ! is thinner than hi_max(ncat).213 !-----------------------------------------------------------------------------! 214 ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons 215 !-----------------------------------------------------------------------------! 216 ! Set hi_max(ncat) to a big value to ensure that all ridged ice 217 ! is thinner than hi_max(ncat). 218 218 219 219 hi_max(jpl) = 999.99 … … 225 225 IF ( con_i) CALL lim_column_sum (jpl, v_i, vt_i_init) 226 226 227 ! Initialize arrays.227 ! Initialize arrays. 228 228 DO jj = 1, jpj 229 229 DO ji = 1, jpi 230 230 231 msnow_mlt(ji,jj) = 0.0232 esnow_mlt(ji,jj) = 0.0233 dardg1dt(ji,jj) = 0.0234 dardg2dt(ji,jj) = 0.0235 dvirdgdt(ji,jj) = 0.0236 opening (ji,jj) = 0.0237 238 !-----------------------------------------------------------------------------!239 ! 2) Dynamical inputs (closing rate, divu_adv, opning)240 !-----------------------------------------------------------------------------!241 !242 ! 2.1 closing_net243 !-----------------244 ! Compute the net rate of closing due to convergence245 ! and shear, based on Flato and Hibler (1995).246 !247 ! The energy dissipation rate is equal to the net closing rate248 ! times the ice strength.249 !250 ! NOTE: The NET closing rate is equal to the rate that open water251 ! area is removed, plus the rate at which ice area is removed by252 ! ridging, minus the rate at which area is added in new ridges.253 ! The GROSS closing rate is equal to the first two terms (open254 ! water closing and thin ice ridging) without the third term255 ! (thick, newly ridged ice).256 257 closing_net(ji,jj) = &258 Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0)259 260 ! 2.2 divu_adv261 !--------------262 ! Compute divu_adv, the divergence rate given by the transport/263 ! advection scheme, which may not be equal to divu as computed264 ! from the velocity field.265 !266 ! If divu_adv < 0, make sure the closing rate is large enough267 ! to give asum = 1.0 after ridging.268 269 divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep270 271 IF (divu_adv(ji,jj) .LT. 0.0) &272 closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj))273 274 ! 2.3 opning275 !------------276 ! Compute the (non-negative) opening rate that will give277 ! asum = 1.0 after ridging.278 opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)231 msnow_mlt(ji,jj) = 0.0 232 esnow_mlt(ji,jj) = 0.0 233 dardg1dt(ji,jj) = 0.0 234 dardg2dt(ji,jj) = 0.0 235 dvirdgdt(ji,jj) = 0.0 236 opening (ji,jj) = 0.0 237 238 !-----------------------------------------------------------------------------! 239 ! 2) Dynamical inputs (closing rate, divu_adv, opning) 240 !-----------------------------------------------------------------------------! 241 ! 242 ! 2.1 closing_net 243 !----------------- 244 ! Compute the net rate of closing due to convergence 245 ! and shear, based on Flato and Hibler (1995). 246 ! 247 ! The energy dissipation rate is equal to the net closing rate 248 ! times the ice strength. 249 ! 250 ! NOTE: The NET closing rate is equal to the rate that open water 251 ! area is removed, plus the rate at which ice area is removed by 252 ! ridging, minus the rate at which area is added in new ridges. 253 ! The GROSS closing rate is equal to the first two terms (open 254 ! water closing and thin ice ridging) without the third term 255 ! (thick, newly ridged ice). 256 257 closing_net(ji,jj) = & 258 Cs*0.5*(Delta_i(ji,jj)-ABS(divu_i(ji,jj))) - MIN(divu_i(ji,jj),0.0) 259 260 ! 2.2 divu_adv 261 !-------------- 262 ! Compute divu_adv, the divergence rate given by the transport/ 263 ! advection scheme, which may not be equal to divu as computed 264 ! from the velocity field. 265 ! 266 ! If divu_adv < 0, make sure the closing rate is large enough 267 ! to give asum = 1.0 after ridging. 268 269 divu_adv(ji,jj) = (1.0-asum(ji,jj)) / rdt_ice ! asum found in ridgeprep 270 271 IF (divu_adv(ji,jj) .LT. 0.0) & 272 closing_net(ji,jj) = max(closing_net(ji,jj), -divu_adv(ji,jj)) 273 274 ! 2.3 opning 275 !------------ 276 ! Compute the (non-negative) opening rate that will give 277 ! asum = 1.0 after ridging. 278 opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj) 279 279 280 280 END DO 281 281 END DO 282 282 283 !-----------------------------------------------------------------------------!284 ! 3) Ridging iteration285 !-----------------------------------------------------------------------------!283 !-----------------------------------------------------------------------------! 284 ! 3) Ridging iteration 285 !-----------------------------------------------------------------------------! 286 286 niter = 1 ! iteration counter 287 287 iterate_ridging = 1 … … 290 290 DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax ) 291 291 292 DO jj = 1, jpj 293 DO ji = 1, jpi 294 295 ! 3.2 closing_gross 296 !-----------------------------------------------------------------------------! 297 ! Based on the ITD of ridging and ridged ice, convert the net 298 ! closing rate to a gross closing rate. 299 ! NOTE: 0 < aksum <= 1 300 closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 301 302 ! correction to closing rate and opening if closing rate is excessive 303 !--------------------------------------------------------------------- 304 ! Reduce the closing rate if more than 100% of the open water 305 ! would be removed. Reduce the opening rate proportionately. 306 IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 307 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 308 IF ( w1 .GT. ato_i(ji,jj)) THEN 309 tmpfac = ato_i(ji,jj) / w1 310 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 311 opning(ji,jj) = opning(ji,jj) * tmpfac 312 ENDIF !w1 313 ENDIF !at0i and athorn 314 315 END DO ! ji 316 END DO ! jj 317 318 ! correction to closing rate / opening if excessive ice removal 319 !--------------------------------------------------------------- 320 ! Reduce the closing rate if more than 100% of any ice category 321 ! would be removed. Reduce the opening rate proportionately. 322 323 DO jl = 1, jpl 324 DO jj = 1, jpj 325 DO ji = 1, jpi 326 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 327 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 328 IF ( w1 .GT. a_i(ji,jj,jl) ) THEN 329 tmpfac = a_i(ji,jj,jl) / w1 330 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 331 opning(ji,jj) = opning(ji,jj) * tmpfac 332 ENDIF 333 ENDIF 334 END DO !ji 335 END DO ! jj 336 END DO !jl 337 338 ! 3.3 Redistribute area, volume, and energy. 339 !-----------------------------------------------------------------------------! 340 341 CALL lim_itd_me_ridgeshift (opning, closing_gross, & 342 msnow_mlt, esnow_mlt) 343 344 ! 3.4 Compute total area of ice plus open water after ridging. 345 !-----------------------------------------------------------------------------! 346 347 CALL lim_itd_me_asumr 348 349 ! 3.5 Do we keep on iterating ??? 350 !-----------------------------------------------------------------------------! 351 ! Check whether asum = 1. If not (because the closing and opening 352 ! rates were reduced above), ridge again with new rates. 353 354 iterate_ridging = 0 355 356 DO jj = 1, jpj 357 DO ji = 1, jpi 358 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 359 closing_net(ji,jj) = 0.0 360 opning(ji,jj) = 0.0 361 ELSE 362 iterate_ridging = 1 363 divu_adv(ji,jj) = (1.0 - asum(ji,jj)) / rdt_ice 364 closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj)) 365 opning(ji,jj) = MAX(0.0, divu_adv(ji,jj)) 366 ENDIF 367 END DO 368 END DO 369 370 IF( lk_mpp ) CALL mpp_max(iterate_ridging) 371 372 ! Repeat if necessary. 373 ! NOTE: If strength smoothing is turned on, the ridging must be 374 ! iterated globally because of the boundary update in the 375 ! smoothing. 376 377 niter = niter + 1 378 379 IF (iterate_ridging == 1) THEN 380 IF (niter .GT. nitermax) THEN 381 WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 382 WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 383 ENDIF 384 CALL lim_itd_me_ridgeprep 385 ENDIF 386 387 END DO !! on the do while over iter 388 389 !-----------------------------------------------------------------------------! 390 ! 4) Ridging diagnostics 391 !-----------------------------------------------------------------------------! 392 ! Convert ridging rate diagnostics to correct units. 393 ! Update fresh water and heat fluxes due to snow melt. 394 395 dti = 1.0/rdt_ice 396 397 asum_error = .false. 398 292 399 DO jj = 1, jpj 293 400 DO ji = 1, jpi 294 401 295 ! 3.2 closing_gross 296 !-----------------------------------------------------------------------------! 297 ! Based on the ITD of ridging and ridged ice, convert the net 298 ! closing rate to a gross closing rate. 299 ! NOTE: 0 < aksum <= 1 300 closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj) 301 302 ! correction to closing rate and opening if closing rate is excessive 303 !--------------------------------------------------------------------- 304 ! Reduce the closing rate if more than 100% of the open water 305 ! would be removed. Reduce the opening rate proportionately. 306 IF ( ato_i(ji,jj) .GT. epsi11 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN 307 w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice 308 IF ( w1 .GT. ato_i(ji,jj)) THEN 309 tmpfac = ato_i(ji,jj) / w1 310 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 311 opning(ji,jj) = opning(ji,jj) * tmpfac 312 ENDIF !w1 313 ENDIF !at0i and athorn 314 315 END DO ! ji 316 END DO ! jj 317 318 ! correction to closing rate / opening if excessive ice removal 319 !--------------------------------------------------------------- 320 ! Reduce the closing rate if more than 100% of any ice category 321 ! would be removed. Reduce the opening rate proportionately. 322 323 DO jl = 1, jpl 324 DO jj = 1, jpj 325 DO ji = 1, jpi 326 IF ( a_i(ji,jj,jl) .GT. epsi11 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN 327 w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice 328 IF ( w1 .GT. a_i(ji,jj,jl) ) THEN 329 tmpfac = a_i(ji,jj,jl) / w1 330 closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac 331 opning(ji,jj) = opning(ji,jj) * tmpfac 332 ENDIF 333 ENDIF 334 END DO !ji 335 END DO ! jj 336 END DO !jl 337 338 ! 3.3 Redistribute area, volume, and energy. 339 !-----------------------------------------------------------------------------! 340 341 CALL lim_itd_me_ridgeshift (opning, closing_gross, & 342 msnow_mlt, esnow_mlt) 343 344 ! 3.4 Compute total area of ice plus open water after ridging. 345 !-----------------------------------------------------------------------------! 346 347 CALL lim_itd_me_asumr 348 349 ! 3.5 Do we keep on iterating ??? 350 !-----------------------------------------------------------------------------! 351 ! Check whether asum = 1. If not (because the closing and opening 352 ! rates were reduced above), ridge again with new rates. 353 354 iterate_ridging = 0 355 356 DO jj = 1, jpj 357 DO ji = 1, jpi 358 IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 359 closing_net(ji,jj) = 0.0 360 opning(ji,jj) = 0.0 361 ELSE 362 iterate_ridging = 1 363 divu_adv(ji,jj) = (1.0 - asum(ji,jj)) / rdt_ice 364 closing_net(ji,jj) = MAX(0.0, -divu_adv(ji,jj)) 365 opning(ji,jj) = MAX(0.0, divu_adv(ji,jj)) 366 ENDIF 367 END DO 368 END DO 369 370 IF( lk_mpp ) CALL mpp_max(iterate_ridging) 371 372 ! Repeat if necessary. 373 ! NOTE: If strength smoothing is turned on, the ridging must be 374 ! iterated globally because of the boundary update in the 375 ! smoothing. 376 377 niter = niter + 1 378 379 IF (iterate_ridging == 1) THEN 380 IF (niter .GT. nitermax) THEN 381 WRITE(numout,*) ' ALERTE : non-converging ridging scheme ' 382 WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging 383 ENDIF 384 CALL lim_itd_me_ridgeprep 385 ENDIF 386 387 END DO !! on the do while over iter 388 389 !-----------------------------------------------------------------------------! 390 ! 4) Ridging diagnostics 391 !-----------------------------------------------------------------------------! 392 ! Convert ridging rate diagnostics to correct units. 393 ! Update fresh water and heat fluxes due to snow melt. 394 395 dti = 1.0/rdt_ice 396 397 asum_error = .false. 398 399 DO jj = 1, jpj 400 DO ji = 1, jpi 401 402 IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 403 404 dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 405 dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 406 dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 407 opening (ji,jj) = opening (ji,jj) * dti 408 409 !-----------------------------------------------------------------------------! 410 ! 5) Heat, salt and freshwater fluxes 411 !-----------------------------------------------------------------------------! 412 ! fresh water source for ocean 413 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti 414 415 ! heat sink for ocean 416 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti 402 IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true. 403 404 dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti 405 dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti 406 dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti 407 opening (ji,jj) = opening (ji,jj) * dti 408 409 !-----------------------------------------------------------------------------! 410 ! 5) Heat, salt and freshwater fluxes 411 !-----------------------------------------------------------------------------! 412 ! fresh water source for ocean 413 fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj)*dti 414 415 ! heat sink for ocean 416 fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj)*dti 417 417 418 418 END DO … … 444 444 ENDIF 445 445 446 !-----------------------------------------------------------------------------!447 ! 6) Updating state variables and trend terms448 !-----------------------------------------------------------------------------!446 !-----------------------------------------------------------------------------! 447 ! 6) Updating state variables and trend terms 448 !-----------------------------------------------------------------------------! 449 449 450 450 CALL lim_var_glo2eqv … … 465 465 d_smv_i_trp(:,:,:) = 0.0 466 466 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 467 d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)467 d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 468 468 469 469 IF(ln_ctl) THEN ! Control print … … 513 513 oa_i(:,:,:) = old_oa_i(:,:,:) 514 514 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 515 smv_i(:,:,:) = old_smv_i(:,:,:)515 smv_i(:,:,:) = old_smv_i(:,:,:) 516 516 517 517 !----------------------------------------------------! … … 528 528 DO ji = 1, jpi 529 529 IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & 530 531 532 530 ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 531 old_e_i(ji,jj,jk,jl) = d_e_i_trp(ji,jj,jk,jl) 532 d_e_i_trp(ji,jj,jk,jl) = 0.0 533 533 ENDIF 534 534 END DO … … 541 541 DO ji = 1, jpi 542 542 IF ( ( old_v_i(ji,jj,jl) .LT. epsi06 ) .AND. & 543 544 545 546 547 548 549 550 551 552 553 554 555 old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)556 543 ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN 544 old_v_i(ji,jj,jl) = d_v_i_trp(ji,jj,jl) 545 d_v_i_trp(ji,jj,jl) = 0.0 546 old_a_i(ji,jj,jl) = d_a_i_trp(ji,jj,jl) 547 d_a_i_trp(ji,jj,jl) = 0.0 548 old_v_s(ji,jj,jl) = d_v_s_trp(ji,jj,jl) 549 d_v_s_trp(ji,jj,jl) = 0.0 550 old_e_s(ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 551 d_e_s_trp(ji,jj,1,jl) = 0.0 552 old_oa_i(ji,jj,jl) = d_oa_i_trp(ji,jj,jl) 553 d_oa_i_trp(ji,jj,jl) = 0.0 554 IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 555 old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 556 d_smv_i_trp(ji,jj,jl) = 0.0 557 557 ENDIF 558 558 END DO 559 559 END DO 560 560 END DO 561 561 562 562 END SUBROUTINE lim_itd_me 563 563 564 !===============================================================================564 !=============================================================================== 565 565 566 566 SUBROUTINE lim_itd_me_icestrength (kstrngth) ! (subroutine 2/6) 567 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 568 !!---------------------------------------------------------------------- 569 !! *** ROUTINE lim_itd_me_icestrength *** 570 !! ** Purpose : 571 !! This routine computes ice strength used in dynamics routines 572 !! of ice thickness 573 !! 574 !! ** Method : 575 !! Compute the strength of the ice pack, defined as the energy (J m-2) 576 !! dissipated per unit area removed from the ice pack under compression, 577 !! and assumed proportional to the change in potential energy caused 578 !! by ridging. Note that only Hibler's formulation is stable and that 579 !! ice strength has to be smoothed 580 !! 581 !! ** Inputs / Ouputs : kstrngth (what kind of ice strength we are using) 582 !! 583 !! ** External : 584 !! 585 !! ** References : 586 !! 587 !!---------------------------------------------------------------------- 588 !! * Arguments 589 590 590 INTEGER, INTENT(in) :: & 591 591 kstrngth ! = 1 for Rothrock formulation, 0 for Hibler (1979) … … 606 606 zworka !: temporary array used here 607 607 608 !------------------------------------------------------------------------------!609 ! 1) Initialize610 !------------------------------------------------------------------------------!608 !------------------------------------------------------------------------------! 609 ! 1) Initialize 610 !------------------------------------------------------------------------------! 611 611 strength(:,:) = 0.0 612 612 613 !------------------------------------------------------------------------------!614 ! 2) Compute thickness distribution of ridging and ridged ice615 !------------------------------------------------------------------------------!613 !------------------------------------------------------------------------------! 614 ! 2) Compute thickness distribution of ridging and ridged ice 615 !------------------------------------------------------------------------------! 616 616 CALL lim_itd_me_ridgeprep 617 617 618 !------------------------------------------------------------------------------!619 ! 3) Rothrock(1975)'s method620 !------------------------------------------------------------------------------!618 !------------------------------------------------------------------------------! 619 ! 3) Rothrock(1975)'s method 620 !------------------------------------------------------------------------------! 621 621 IF (kstrngth == 1) then 622 622 … … 626 626 627 627 IF( ( a_i(ji,jj,jl) .GT. epsi11 ) & 628 .AND. ( athorn(ji,jj,jl) .GT. 0.0 ) ) THEN628 .AND. ( athorn(ji,jj,jl) .GT. 0.0 ) ) THEN 629 629 hi = v_i(ji,jj,jl) / a_i(ji,jj,jl) 630 630 !---------------------------- … … 632 632 !---------------------------- 633 633 strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * & 634 hi * hi634 hi * hi 635 635 636 636 !-------------------------- … … 638 638 !-------------------------- 639 639 strength(ji,jj) = strength(ji,jj) + 2.0 * araft(ji,jj,jl) & 640 * hi * hi640 * hi * hi 641 641 642 642 !---------------------------- … … 644 644 !---------------------------- 645 645 strength(ji,jj) = strength(ji,jj) & 646 + aridge(ji,jj,jl)/krdg(ji,jj,jl) &647 * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) &648 / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl))646 + aridge(ji,jj,jl)/krdg(ji,jj,jl) & 647 * 1.0/3.0 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) & 648 / (hrmax(ji,jj,jl)-hrmin(ji,jj,jl)) 649 649 ENDIF ! aicen > epsi11 650 650 … … 656 656 DO ji = 1, jpi 657 657 strength(ji,jj) = Cf * Cp * strength(ji,jj) / aksum(ji,jj) 658 659 660 658 ! Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) 659 ! Cf accounts for frictional dissipation 660 661 661 END DO ! j 662 662 END DO ! i … … 664 664 ksmooth = 1 665 665 666 !------------------------------------------------------------------------------!667 ! 4) Hibler (1979)' method668 !------------------------------------------------------------------------------!666 !------------------------------------------------------------------------------! 667 ! 4) Hibler (1979)' method 668 !------------------------------------------------------------------------------! 669 669 ELSE ! kstrngth ne 1: Hibler (1979) form 670 670 … … 679 679 ENDIF ! kstrngth 680 680 681 !682 !------------------------------------------------------------------------------!683 ! 5) Impact of brine volume684 !------------------------------------------------------------------------------!685 ! CAN BE REMOVED686 !681 ! 682 !------------------------------------------------------------------------------! 683 ! 5) Impact of brine volume 684 !------------------------------------------------------------------------------! 685 ! CAN BE REMOVED 686 ! 687 687 IF ( brinstren_swi .EQ. 1 ) THEN 688 688 … … 700 700 ENDIF 701 701 702 !703 !------------------------------------------------------------------------------!704 ! 6) Smoothing ice strength705 !------------------------------------------------------------------------------!706 !702 ! 703 !------------------------------------------------------------------------------! 704 ! 6) Smoothing ice strength 705 !------------------------------------------------------------------------------! 706 ! 707 707 !------------------- 708 708 ! Spatial smoothing … … 715 715 DO ji = 2, jpi - 1 716 716 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 717 717 ! present 718 718 zworka(ji,jj) = 4.0 * strength(ji,jj) & 719 720 721 722 719 + strength(ji-1,jj) * tms(ji-1,jj) & 720 + strength(ji+1,jj) * tms(ji+1,jj) & 721 + strength(ji,jj-1) * tms(ji,jj-1) & 722 + strength(ji,jj+1) * tms(ji,jj+1) 723 723 724 724 zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) & 725 725 + tms(ji,jj-1) + tms(ji,jj+1) 726 726 zworka(ji,jj) = zworka(ji,jj) / zw1 727 727 ELSE … … 749 749 750 750 IF ( ksmooth .EQ. 2 ) THEN 751 752 751 752 753 753 CALL lbc_lnk( strength, 'T', 1. ) 754 754 755 755 DO jj = 1, jpj - 1 756 756 DO ji = 1, jpi - 1 757 757 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi11) THEN ! ice is 758 758 ! present 759 759 numts_rm = 1 ! number of time steps for the running mean 760 760 IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 761 761 IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1 762 762 zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / & 763 763 numts_rm 764 764 strp2(ji,jj) = strp1(ji,jj) 765 765 strp1(ji,jj) = strength(ji,jj) … … 771 771 772 772 ENDIF ! ksmooth 773 773 774 774 ! Boundary conditions 775 775 CALL lbc_lnk( strength, 'T', 1. ) 776 776 777 778 779 !===============================================================================780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 777 END SUBROUTINE lim_itd_me_icestrength 778 779 !=============================================================================== 780 781 SUBROUTINE lim_itd_me_ridgeprep !(subroutine 3/6) 782 783 !!---------------------------------------------------------------------! 784 !! *** ROUTINE lim_itd_me_ridgeprep *** 785 !! ** Purpose : 786 !! preparation for ridging and strength calculations 787 !! 788 !! ** Method : 789 !! Compute the thickness distribution of the ice and open water 790 !! participating in ridging and of the resulting ridges. 791 !! 792 !! ** Arguments : 793 !! 794 !! ** External : 795 !! 796 !!---------------------------------------------------------------------! 797 !! * Arguments 798 799 799 INTEGER :: & 800 800 ji,jj, & ! horizontal indices … … 820 820 epsi06 = 1.0e-6 821 821 822 !------------------------------------------------------------------------------!822 !------------------------------------------------------------------------------! 823 823 824 824 Gstari = 1.0/Gstar … … 833 833 krdg (:,:,:) = 1.0 834 834 835 ! ! Zero out categories with very small areas835 ! ! Zero out categories with very small areas 836 836 CALL lim_itd_me_zapsmall 837 837 838 !------------------------------------------------------------------------------!839 ! 1) Participation function840 !------------------------------------------------------------------------------!838 !------------------------------------------------------------------------------! 839 ! 1) Participation function 840 !------------------------------------------------------------------------------! 841 841 842 842 ! Compute total area of ice plus open water. … … 886 886 DO ji = 1, jpi 887 887 Gsum(ji,jj,jl) = Gsum(ji,jj,jl) * zworka(ji,jj) 888 END DO 888 END DO 889 889 END DO 890 890 END DO 891 891 892 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)893 !--------------------------------------------------------------------------------------------------894 ! Compute the participation function athorn; this is analogous to895 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).896 ! area lost from category n due to ridging/closing897 ! athorn(n) = total area lost due to ridging/closing898 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar).899 !900 ! The expressions for athorn are found by integrating b(h)g(h) between901 ! the category boundaries.902 !-----------------------------------------------------------------892 ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn) 893 !-------------------------------------------------------------------------------------------------- 894 ! Compute the participation function athorn; this is analogous to 895 ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975). 896 ! area lost from category n due to ridging/closing 897 ! athorn(n) = total area lost due to ridging/closing 898 ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 899 ! 900 ! The expressions for athorn are found by integrating b(h)g(h) between 901 ! the category boundaries. 902 !----------------------------------------------------------------- 903 903 904 904 krdg_index = 1 … … 906 906 IF ( krdg_index .EQ. 0 ) THEN 907 907 908 !--- Linear formulation (Thorndike et al., 1975)909 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates910 DO jj = 1, jpj911 DO ji = 1, jpi912 IF (Gsum(ji,jj,jl) < Gstar) THEN913 athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * &914 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari)915 ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN916 athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * &917 (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari)918 ELSE919 athorn(ji,jj,jl) = 0.0920 ENDIF921 END DO ! ji922 END DO ! jj923 END DO ! jl908 !--- Linear formulation (Thorndike et al., 1975) 909 DO jl = 0, ice_cat_bounds(1,2) ! only undeformed ice participates 910 DO jj = 1, jpj 911 DO ji = 1, jpi 912 IF (Gsum(ji,jj,jl) < Gstar) THEN 913 athorn(ji,jj,jl) = Gstari * (Gsum(ji,jj,jl)-Gsum(ji,jj,jl-1)) * & 914 (2.0 - (Gsum(ji,jj,jl-1)+Gsum(ji,jj,jl))*Gstari) 915 ELSEIF (Gsum(ji,jj,jl-1) < Gstar) THEN 916 athorn(ji,jj,jl) = Gstari * (Gstar-Gsum(ji,jj,jl-1)) * & 917 (2.0 - (Gsum(ji,jj,jl-1)+Gstar)*Gstari) 918 ELSE 919 athorn(ji,jj,jl) = 0.0 920 ENDIF 921 END DO ! ji 922 END DO ! jj 923 END DO ! jl 924 924 925 925 ELSE ! krdg_index = 1 926 927 !--- Exponential, more stable formulation (Lipscomb et al, 2007)928 ! precompute exponential terms using Gsum as a work array929 zdummy = 1.0 / (1.0-EXP(-astari))930 931 DO jl = -1, jpl932 DO jj = 1, jpj933 DO ji = 1, jpi934 Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy935 END DO !ji936 END DO !jj937 END DO !jl938 939 ! compute athorn940 DO jl = 0, ice_cat_bounds(1,2)941 DO jj = 1, jpj942 DO ji = 1, jpi943 athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl)944 END DO !ji945 END DO ! jj946 END DO !jl926 927 !--- Exponential, more stable formulation (Lipscomb et al, 2007) 928 ! precompute exponential terms using Gsum as a work array 929 zdummy = 1.0 / (1.0-EXP(-astari)) 930 931 DO jl = -1, jpl 932 DO jj = 1, jpj 933 DO ji = 1, jpi 934 Gsum(ji,jj,jl) = EXP(-Gsum(ji,jj,jl)*astari)*zdummy 935 END DO !ji 936 END DO !jj 937 END DO !jl 938 939 ! compute athorn 940 DO jl = 0, ice_cat_bounds(1,2) 941 DO jj = 1, jpj 942 DO ji = 1, jpi 943 athorn(ji,jj,jl) = Gsum(ji,jj,jl-1) - Gsum(ji,jj,jl) 944 END DO !ji 945 END DO ! jj 946 END DO !jl 947 947 948 948 ENDIF ! krdg_index … … 956 956 IF ( athorn(ji,jj,jl) .GT. 0.0 ) THEN 957 957 aridge(ji,jj,jl) = ( TANH ( Craft * ( ht_i(ji,jj,jl) - & 958 959 958 hparmeter ) ) + 1.0 ) / 2.0 * & 959 athorn(ji,jj,jl) 960 960 araft (ji,jj,jl) = ( TANH ( - Craft * ( ht_i(ji,jj,jl) - & 961 962 961 hparmeter ) ) + 1.0 ) / 2.0 * & 962 athorn(ji,jj,jl) 963 963 IF ( araft(ji,jj,jl) .LT. epsi06 ) araft(ji,jj,jl) = 0.0 964 964 aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0) … … 982 982 IF ( raftswi .EQ. 1 ) THEN 983 983 984 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 985 DO jl = 1, jpl 986 DO jj = 1, jpj 987 DO ji = 1, jpi 988 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 989 epsi11 ) THEN 990 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 991 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 992 WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj) 993 WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl) 994 WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl) 995 WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl) 996 ENDIF 984 IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi11 ) THEN 985 DO jl = 1, jpl 986 DO jj = 1, jpj 987 DO ji = 1, jpi 988 IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. & 989 epsi11 ) THEN 990 WRITE(numout,*) ' ALERTE 96 : wrong participation function ... ' 991 WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl 992 WRITE(numout,*) ' lat, lon : ', gphit(ji,jj), glamt(ji,jj) 993 WRITE(numout,*) ' aridge : ', aridge(ji,jj,1:jpl) 994 WRITE(numout,*) ' araft : ', araft(ji,jj,1:jpl) 995 WRITE(numout,*) ' athorn : ', athorn(ji,jj,1:jpl) 996 ENDIF 997 END DO 997 998 END DO 998 999 END DO 999 END DO 1000 ENDIF 1001 1000 1002 ENDIF 1001 1003 1002 ENDIF 1003 1004 !----------------------------------------------------------------- 1005 ! 2) Transfer function 1006 !----------------------------------------------------------------- 1007 ! Compute max and min ridged ice thickness for each ridging category. 1008 ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 1009 ! 1010 ! This parameterization is a modified version of Hibler (1980). 1011 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 1012 ! and for very thick ridging ice must be >= krdgmin*hi 1013 ! 1014 ! The minimum ridging thickness, hrmin, is equal to 2*hi 1015 ! (i.e., rafting) and for very thick ridging ice is 1016 ! constrained by hrmin <= (hrmean + hi)/2. 1017 ! 1018 ! The maximum ridging thickness, hrmax, is determined by 1019 ! hrmean and hrmin. 1020 ! 1021 ! These modifications have the effect of reducing the ice strength 1022 ! (relative to the Hibler formulation) when very thick ice is 1023 ! ridging. 1024 ! 1025 ! aksum = net area removed/ total area removed 1026 ! where total area removed = area of ice that ridges 1027 ! net area removed = total area removed - area of new ridges 1028 !----------------------------------------------------------------- 1004 !----------------------------------------------------------------- 1005 ! 2) Transfer function 1006 !----------------------------------------------------------------- 1007 ! Compute max and min ridged ice thickness for each ridging category. 1008 ! Assume ridged ice is uniformly distributed between hrmin and hrmax. 1009 ! 1010 ! This parameterization is a modified version of Hibler (1980). 1011 ! The mean ridging thickness, hrmean, is proportional to hi^(0.5) 1012 ! and for very thick ridging ice must be >= krdgmin*hi 1013 ! 1014 ! The minimum ridging thickness, hrmin, is equal to 2*hi 1015 ! (i.e., rafting) and for very thick ridging ice is 1016 ! constrained by hrmin <= (hrmean + hi)/2. 1017 ! 1018 ! The maximum ridging thickness, hrmax, is determined by 1019 ! hrmean and hrmin. 1020 ! 1021 ! These modifications have the effect of reducing the ice strength 1022 ! (relative to the Hibler formulation) when very thick ice is 1023 ! ridging. 1024 ! 1025 ! aksum = net area removed/ total area removed 1026 ! where total area removed = area of ice that ridges 1027 ! net area removed = total area removed - area of new ridges 1028 !----------------------------------------------------------------- 1029 1029 1030 1030 ! Transfer function … … 1062 1062 DO ji = 1, jpi 1063 1063 aksum(ji,jj) = aksum(ji,jj) & 1064 1065 1064 + aridge(ji,jj,jl) * (1.0 - 1.0/krdg(ji,jj,jl)) & 1065 + araft (ji,jj,jl) * (1.0 - 1.0/kraft) 1066 1066 END DO 1067 1067 END DO 1068 1068 END DO 1069 1069 1070 1071 1072 !===============================================================================1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1070 END SUBROUTINE lim_itd_me_ridgeprep 1071 1072 !=============================================================================== 1073 1074 SUBROUTINE lim_itd_me_ridgeshift(opning, closing_gross, & 1075 msnow_mlt, esnow_mlt) ! (subroutine 4/6) 1076 1077 !!----------------------------------------------------------------------------- 1078 !! *** ROUTINE lim_itd_me_icestrength *** 1079 !! ** Purpose : 1080 !! This routine shift ridging ice among thickness categories 1081 !! of ice thickness 1082 !! 1083 !! ** Method : 1084 !! Remove area, volume, and energy from each ridging category 1085 !! and add to thicker ice categories. 1086 !! 1087 !! ** Arguments : 1088 !! 1089 !! ** Inputs / Ouputs : 1090 !! 1091 !! ** External : 1092 !! 1093 1093 1094 1094 REAL (wp), DIMENSION(jpi,jpj), INTENT(IN) :: & 1095 1095 opning, & ! rate of opening due to divergence/shear 1096 1096 closing_gross ! rate at which area removed, not counting 1097 1097 ! area of new ridges 1098 1098 1099 1099 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: & … … 1176 1176 LOGICAL, PARAMETER :: & 1177 1177 l_conservation_check = .true. ! if true, check conservation 1178 1178 ! (useful for debugging) 1179 1179 LOGICAL :: & 1180 1180 neg_ato_i , & ! flag for ato_i(i,j) < -puny … … 1187 1187 zindb ! switch for the presence of ridge poros or not 1188 1188 1189 !----------------------------------------------------------------------------1189 !---------------------------------------------------------------------------- 1190 1190 1191 1191 ! Conservation check … … 1202 1202 epsi10 = 1.0d-10 1203 1203 1204 !-------------------------------------------------------------------------------1205 ! 1) Compute change in open water area due to closing and opening.1206 !-------------------------------------------------------------------------------1204 !------------------------------------------------------------------------------- 1205 ! 1) Compute change in open water area due to closing and opening. 1206 !------------------------------------------------------------------------------- 1207 1207 1208 1208 neg_ato_i = .false. … … 1211 1211 DO ji = 1, jpi 1212 1212 ato_i(ji,jj) = ato_i(ji,jj) & 1213 1214 1213 - athorn(ji,jj,0)*closing_gross(ji,jj)*rdt_ice & 1214 + opning(ji,jj)*rdt_ice 1215 1215 IF (ato_i(ji,jj) .LT. -epsi11) THEN 1216 1216 neg_ato_i = .true. … … 1234 1234 ENDIF ! neg_ato_i 1235 1235 1236 !-----------------------------------------------------------------1237 ! 2) Save initial state variables1238 !-----------------------------------------------------------------1236 !----------------------------------------------------------------- 1237 ! 2) Save initial state variables 1238 !----------------------------------------------------------------- 1239 1239 1240 1240 DO jl = 1, jpl … … 1252 1252 1253 1253 esnon_init(:,:,:) = e_s(:,:,1,:) 1254 1254 1255 1255 DO jl = 1, jpl 1256 1256 DO jk = 1, nlay_i … … 1263 1263 END DO !jl 1264 1264 1265 !1266 !-----------------------------------------------------------------1267 ! 3) Pump everything from ice which is being ridged / rafted1268 !-----------------------------------------------------------------1269 ! Compute the area, volume, and energy of ice ridging in each1270 ! category, along with the area of the resulting ridge.1265 ! 1266 !----------------------------------------------------------------- 1267 ! 3) Pump everything from ice which is being ridged / rafted 1268 !----------------------------------------------------------------- 1269 ! Compute the area, volume, and energy of ice ridging in each 1270 ! category, along with the area of the resulting ridge. 1271 1271 1272 1272 DO jl1 = 1, jpl !jl1 describes the ridging category 1273 1273 1274 !------------------------------------------------1275 ! 3.1) Identify grid cells with nonzero ridging1276 !------------------------------------------------1274 !------------------------------------------------ 1275 ! 3.1) Identify grid cells with nonzero ridging 1276 !------------------------------------------------ 1277 1277 1278 1278 icells = 0 … … 1280 1280 DO ji = 1, jpi 1281 1281 IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0 & 1282 1282 .AND. closing_gross(ji,jj) > 0.0) THEN 1283 1283 icells = icells + 1 1284 1284 indxi(icells) = ji … … 1296 1296 jj = indxj(ij) 1297 1297 1298 !--------------------------------------------------------------------1299 ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)1300 !--------------------------------------------------------------------1298 !-------------------------------------------------------------------- 1299 ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2) 1300 !-------------------------------------------------------------------- 1301 1301 1302 1302 ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice … … 1310 1310 oirft2(ji,jj)= oirft1(ji,jj) / kraft 1311 1311 1312 !---------------------------------------------------------------1313 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=11314 !---------------------------------------------------------------1312 !--------------------------------------------------------------- 1313 ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 1314 !--------------------------------------------------------------- 1315 1315 1316 1316 afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging … … 1328 1328 ENDIF 1329 1329 1330 !--------------------------------------------------------------------------1331 ! 3.4) Subtract area, volume, and energy from ridging1332 ! / rafting category n1.1333 !--------------------------------------------------------------------------1330 !-------------------------------------------------------------------------- 1331 ! 3.4) Subtract area, volume, and energy from ridging 1332 ! / rafting category n1. 1333 !-------------------------------------------------------------------------- 1334 1334 vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj) / & 1335 1335 ( 1.0 + ridge_por ) 1336 1336 vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + ridge_por ) 1337 1337 vsw (ji,jj) = vrdg1(ji,jj) * ridge_por … … 1340 1340 esrdg(ji,jj) = esnon_init(ji,jj,jl1) * afrac(ji,jj) 1341 1341 srdg1(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) / & 1342 1342 ( 1. + ridge_por ) 1343 1343 srdg2(ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj) 1344 1344 … … 1357 1357 smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1(ji,jj) - smrft(ji,jj) 1358 1358 1359 !-----------------------------------------------------------------1360 ! 3.5) Compute properties of new ridges1361 !-----------------------------------------------------------------1359 !----------------------------------------------------------------- 1360 ! 3.5) Compute properties of new ridges 1361 !----------------------------------------------------------------- 1362 1362 !------------- 1363 1363 ! Salinity … … 1373 1373 ! salt flux due to ridge creation 1374 1374 fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + & 1375 MAX ( zdummy - srdg2(ji,jj) , 0.0 ) &1376 * rhoic / rdt_ice1375 MAX ( zdummy - srdg2(ji,jj) , 0.0 ) & 1376 * rhoic / rdt_ice 1377 1377 1378 1378 ! sal times volume for new ridges 1379 1379 srdg2(ji,jj) = sm_newridge * vrdg2(ji,jj) 1380 1380 1381 !------------------------------------1382 ! 3.6 Increment ridging diagnostics1383 !------------------------------------1384 1385 ! jl1 looping 1-jpl1386 ! ij looping 1-icells1381 !------------------------------------ 1382 ! 3.6 Increment ridging diagnostics 1383 !------------------------------------ 1384 1385 ! jl1 looping 1-jpl 1386 ! ij looping 1-icells 1387 1387 1388 1388 dardg1dt(ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) … … 1393 1393 IF (con_i) vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj) 1394 1394 1395 !------------------------------------------1396 ! 3.7 Put the snow somewhere in the ocean1397 !------------------------------------------1398 1399 ! Place part of the snow lost by ridging into the ocean.1400 ! Note that esnow_mlt < 0; the ocean must cool to melt snow.1401 ! If the ocean temp = Tf already, new ice must grow.1402 ! During the next time step, thermo_rates will determine whether1403 ! the ocean cools or new ice grows.1404 ! jl1 looping 1-jpl1405 ! ij looping 1-icells1406 1395 !------------------------------------------ 1396 ! 3.7 Put the snow somewhere in the ocean 1397 !------------------------------------------ 1398 1399 ! Place part of the snow lost by ridging into the ocean. 1400 ! Note that esnow_mlt < 0; the ocean must cool to melt snow. 1401 ! If the ocean temp = Tf already, new ice must grow. 1402 ! During the next time step, thermo_rates will determine whether 1403 ! the ocean cools or new ice grows. 1404 ! jl1 looping 1-jpl 1405 ! ij looping 1-icells 1406 1407 1407 msnow_mlt(ji,jj) = msnow_mlt(ji,jj) & 1408 1409 !rafting included1410 1408 + rhosn*vsrdg(ji,jj)*(1.0-fsnowrdg) & 1409 !rafting included 1410 + rhosn*vsrft(ji,jj)*(1.0-fsnowrft) 1411 1411 1412 1412 esnow_mlt(ji,jj) = esnow_mlt(ji,jj) & 1413 1414 !rafting included1415 1416 1417 !-----------------------------------------------------------------1418 ! 3.8 Compute quantities used to apportion ice among categories1419 ! in the n2 loop below1420 !-----------------------------------------------------------------1421 1422 ! jl1 looping 1-jpl1423 ! ij looping 1-icells1413 + esrdg(ji,jj)*(1.0-fsnowrdg) & 1414 !rafting included 1415 + esrft(ji,jj)*(1.0-fsnowrft) 1416 1417 !----------------------------------------------------------------- 1418 ! 3.8 Compute quantities used to apportion ice among categories 1419 ! in the n2 loop below 1420 !----------------------------------------------------------------- 1421 1422 ! jl1 looping 1-jpl 1423 ! ij looping 1-icells 1424 1424 1425 1425 dhr(ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) 1426 1426 dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) & 1427 1427 - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) 1428 1428 1429 1429 1430 1430 END DO ! ij 1431 1431 1432 !--------------------------------------------------------------------1433 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and1434 ! compute ridged ice enthalpy1435 !--------------------------------------------------------------------1432 !-------------------------------------------------------------------- 1433 ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and 1434 ! compute ridged ice enthalpy 1435 !-------------------------------------------------------------------- 1436 1436 DO jk = 1, nlay_i 1437 1437 !CDIR NODEP 1438 1438 DO ij = 1, icells 1439 ji = indxi(ij)1440 jj = indxj(ij)1441 ! heat content of ridged ice1442 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / &1443 1444 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)1445 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) &1446 1447 1448 ! sea water heat content1449 ztmelts = - tmut * sss_m(ji,jj) + rtt1450 ! heat content per unit volume1451 zdummy0 = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj)1452 1453 ! corrected sea water salinity1454 zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) )1455 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / &1456 1457 1458 ztmelts = - tmut * zdummy + rtt1459 ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj)1460 1461 ! heat flux1462 fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / &1463 1464 1465 ! Correct dimensions to avoid big values1466 ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+091467 1468 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J1469 ersw(ji,jj,jk) = ersw(ji,jj,jk) * &1470 1471 1472 1473 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk)1439 ji = indxi(ij) 1440 jj = indxj(ij) 1441 ! heat content of ridged ice 1442 erdg1(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) / & 1443 ( 1.0 + ridge_por ) 1444 eirft(ji,jj,jk) = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj) 1445 e_i(ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) & 1446 - erdg1(ji,jj,jk) & 1447 - eirft(ji,jj,jk) 1448 ! sea water heat content 1449 ztmelts = - tmut * sss_m(ji,jj) + rtt 1450 ! heat content per unit volume 1451 zdummy0 = - rcp * ( sst_m(ji,jj) + rt0 - rtt ) * vsw(ji,jj) 1452 1453 ! corrected sea water salinity 1454 zindb = MAX( 0.0, SIGN( 1.0, vsw(ji,jj) - zeps ) ) 1455 zdummy = zindb * ( srdg1(ji,jj) - srdg2(ji,jj) ) / & 1456 MAX( ridge_por * vsw(ji,jj), zeps ) 1457 1458 ztmelts = - tmut * zdummy + rtt 1459 ersw(ji,jj,jk) = - rcp * ( ztmelts - rtt ) * vsw(ji,jj) 1460 1461 ! heat flux 1462 fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / & 1463 rdt_ice 1464 1465 ! Correct dimensions to avoid big values 1466 ersw(ji,jj,jk) = ersw(ji,jj,jk) / 1.0d+09 1467 1468 ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 1469 ersw(ji,jj,jk) = ersw(ji,jj,jk) * & 1470 area(ji,jj) * vsw(ji,jj) / & 1471 nlay_i 1472 1473 erdg2(ji,jj,jk) = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 1474 1474 END DO ! ij 1475 1475 END DO !jk … … 1483 1483 jj = indxj(ij) 1484 1484 eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - & 1485 erdg1(ji,jj,jk)1485 erdg1(ji,jj,jk) 1486 1486 END DO ! ij 1487 1487 END DO !jk … … 1497 1497 WRITE(numout,*) ' ardg > a_i' 1498 1498 WRITE(numout,*) ' ardg, aicen_init : ', & 1499 1499 ardg1(ji,jj), aicen_init(ji,jj,jl1) 1500 1500 ENDIF ! afrac > 1 + puny 1501 1501 ENDDO ! if … … 1510 1510 WRITE(numout,*) ' arft > a_i' 1511 1511 WRITE(numout,*) ' arft, aicen_init : ', & 1512 1512 arft1(ji,jj), aicen_init(ji,jj,jl1) 1513 1513 ENDIF ! afrft > 1 + puny 1514 1514 ENDDO ! if 1515 1515 ENDIF ! large_afrft 1516 1516 1517 !-------------------------------------------------------------------------------1518 ! 4) Add area, volume, and energy of new ridge to each category jl21519 !-------------------------------------------------------------------------------1520 ! jl1 looping 1-jpl1517 !------------------------------------------------------------------------------- 1518 ! 4) Add area, volume, and energy of new ridge to each category jl2 1519 !------------------------------------------------------------------------------- 1520 ! jl1 looping 1-jpl 1521 1521 DO jl2 = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 1522 ! over categories to which ridged ice is transferred1522 ! over categories to which ridged ice is transferred 1523 1523 !CDIR NODEP 1524 1524 DO ij = 1, icells … … 1531 1531 1532 1532 IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR. & 1533 1533 hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN 1534 1534 hL = 0.0 1535 1535 hR = 0.0 … … 1546 1546 v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + fvol(ji,jj) * vrdg2(ji,jj) 1547 1547 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & 1548 1548 + fvol(ji,jj) * vsrdg(ji,jj) * fsnowrdg 1549 1549 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & 1550 1550 + fvol(ji,jj) * esrdg(ji,jj) * fsnowrdg 1551 1551 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) + fvol(ji,jj) * srdg2(ji,jj) 1552 1552 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) + farea * oirdg2(ji,jj) … … 1561 1561 jj = indxj(ij) 1562 1562 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & 1563 1563 + fvol(ji,jj)*erdg2(ji,jj,jk) 1564 1564 END DO ! ij 1565 1565 END DO !jk … … 1576 1576 ! Compute the fraction of rafted ice area and volume going to 1577 1577 ! thickness category jl2, transfer area, volume, and energy accordingly. 1578 1578 1579 1579 IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. & 1580 1580 hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 1581 1581 a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj) 1582 1582 v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj) 1583 1583 v_s(ji,jj,jl2) = v_s(ji,jj,jl2) & 1584 1584 + vsrft(ji,jj)*fsnowrft 1585 1585 e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) & 1586 1586 + esrft(ji,jj)*fsnowrft 1587 1587 smv_i(ji,jj,jl2) = smv_i(ji,jj,jl2) & 1588 1588 + smrft(ji,jj) 1589 1589 oa_i(ji,jj,jl2) = oa_i(ji,jj,jl2) & 1590 1590 + oirft2(ji,jj) 1591 1591 ENDIF ! hraft 1592 1592 … … 1600 1600 jj = indxj(ij) 1601 1601 IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND. & 1602 1602 hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN 1603 1603 e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) & 1604 1604 + eirft(ji,jj,jk) 1605 1605 ENDIF 1606 1606 END DO ! ij … … 1628 1628 END SUBROUTINE lim_itd_me_ridgeshift 1629 1629 1630 !==============================================================================1630 !============================================================================== 1631 1631 1632 1632 SUBROUTINE lim_itd_me_asumr !(subroutine 5/6) 1633 1633 1634 1635 1636 1637 1638 1639 1640 1641 1642 1643 1644 1645 1646 1634 !!----------------------------------------------------------------------------- 1635 !! *** ROUTINE lim_itd_me_asumr *** 1636 !! ** Purpose : 1637 !! This routine finds total fractional area 1638 !! 1639 !! ** Method : 1640 !! Find the total area of ice plus open water in each grid cell. 1641 !! 1642 !! This is similar to the aggregate_area subroutine except that the 1643 !! total area can be greater than 1, so the open water area is 1644 !! included in the sum instead of being computed as a residual. 1645 !! 1646 !! ** Arguments : 1647 1647 1648 1648 INTEGER :: ji, jj, jl … … 1672 1672 END SUBROUTINE lim_itd_me_asumr 1673 1673 1674 !==============================================================================1674 !============================================================================== 1675 1675 1676 1676 SUBROUTINE lim_itd_me_init ! (subroutine 6/6) … … 1691 1691 !!------------------------------------------------------------------- 1692 1692 NAMELIST/namiceitdme/ ridge_scheme_swi, Cs, Cf, fsnowrdg, fsnowrft,& 1693 1694 1695 1696 1693 Gstar, astar, & 1694 Hstar, raftswi, hparmeter, Craft, ridge_por, & 1695 sal_max_ridge, partfun_swi, transfun_swi, & 1696 brinstren_swi 1697 1697 !!------------------------------------------------------------------- 1698 1698 … … 1725 1725 END SUBROUTINE lim_itd_me_init 1726 1726 1727 !==============================================================================1727 !============================================================================== 1728 1728 1729 1729 SUBROUTINE lim_itd_me_zapsmall … … 1743 1743 1744 1744 INTEGER :: & 1745 1746 1747 1748 ! ij, & ! combined i/j horizontal index1749 1750 1751 ! INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &1752 ! indxi, & ! compressed indices for i/j directions1753 ! indxj1745 ji,jj, & ! horizontal indices 1746 jl, & ! ice category index 1747 jk, & ! ice layer index 1748 ! ij, & ! combined i/j horizontal index 1749 icells ! number of cells with ice to zap 1750 1751 ! INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: & 1752 ! indxi, & ! compressed indices for i/j directions 1753 ! indxj 1754 1754 1755 1755 INTEGER, DIMENSION(jpi,jpj) :: zmask … … 1757 1757 1758 1758 REAL(wp) :: & 1759 1759 xtmp ! temporary variable 1760 1760 1761 1761 DO jl = 1, jpl 1762 1762 1763 !-----------------------------------------------------------------1764 ! Count categories to be zapped.1765 ! Abort model in case of negative area.1766 !-----------------------------------------------------------------1767 IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 ) THEN1763 !----------------------------------------------------------------- 1764 ! Count categories to be zapped. 1765 ! Abort model in case of negative area. 1766 !----------------------------------------------------------------- 1767 IF( MINVAL(a_i(:,:,jl)) .LT. -epsi11 .AND. ln_nicep ) THEN 1768 1768 DO jj = 1, jpj 1769 1769 DO ji = 1, jpi … … 1774 1774 ENDIF 1775 1775 END DO 1776 END DO 1776 END DO 1777 1777 ENDIF 1778 1779 icells = 0 1780 zmask = 0.e0 1781 DO jj = 1, jpj 1782 DO ji = 1, jpi 1783 IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) & 1784 .OR. & 1785 ( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) & 1786 .OR. & 1787 !new line 1788 ( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) & 1789 .OR. & 1790 ( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN 1791 zmask(ji,jj) = 1 1792 ENDIF 1778 1779 icells = 0 1780 zmask = 0.e0 1781 DO jj = 1, jpj 1782 DO ji = 1, jpi 1783 IF ( ( a_i(ji,jj,jl) .GE. -epsi11 .AND. a_i(ji,jj,jl) .LT. 0.0) & 1784 .OR. & 1785 ( a_i(ji,jj,jl) .GT. 0.0 .AND. a_i(ji,jj,jl) .LE. 1.0e-11 ) & 1786 .OR. & 1787 !new line 1788 ( v_i(ji,jj,jl) .EQ. 0.0 .AND. a_i(ji,jj,jl) .GT. 0.0 ) & 1789 .OR. & 1790 ( v_i(ji,jj,jl) .GT. 0.0 .AND. v_i(ji,jj,jl) .LT. 1.e-12 ) ) THEN 1791 zmask(ji,jj) = 1 1792 ENDIF 1793 END DO 1793 1794 END DO 1794 END DO 1795 WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 1796 1797 !----------------------------------------------------------------- 1798 ! Zap ice energy and use ocean heat to melt ice 1799 !----------------------------------------------------------------- 1795 IF( ln_nicep ) WRITE(numout,*) SUM(zmask), ' cells of ice zapped in the ocean ' 1796 1797 !----------------------------------------------------------------- 1798 ! Zap ice energy and use ocean heat to melt ice 1799 !----------------------------------------------------------------- 1800 1800 1801 1801 DO jk = 1, nlay_i … … 1803 1803 DO ji = 1 , jpi 1804 1804 1805 xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice1806 xtmp = xtmp * unit_fac1807 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp1808 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) )1805 xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice 1806 xtmp = xtmp * unit_fac 1807 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1808 e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1 - zmask(ji,jj) ) 1809 1809 END DO ! ji 1810 1810 END DO ! jj … … 1814 1814 DO ji = 1 , jpi 1815 1815 1816 !-----------------------------------------------------------------1817 ! Zap snow energy and use ocean heat to melt snow1818 !-----------------------------------------------------------------1819 1820 ! xtmp = esnon(i,j,n) / dt ! < 01821 ! fhnet(i,j) = fhnet(i,j) + xtmp1822 ! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp1823 ! xtmp is greater than 01824 ! fluxes are positive to the ocean1825 ! here the flux has to be negative for the ocean1826 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice1827 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp1828 1829 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB ???????1830 1831 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) )1832 1833 !-----------------------------------------------------------------1834 ! zap ice and snow volume, add water and salt to ocean1835 !-----------------------------------------------------------------1836 1837 ! xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt1838 ! fresh(i,j) = fresh(i,j) + xtmp1839 ! fresh_hist(i,j) = fresh_hist(i,j) + xtmp1840 1841 ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) ) * &1842 ! rhosn * v_s(ji,jj,jl) / rdt_ice1843 1844 ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * &1845 ! rhoic * v_i(ji,jj,jl) / rdt_ice1846 1847 ! emps(i,j) = emps(i,j) + xtmp1848 ! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp1849 1850 ato_i(ji,jj) = a_i(ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj)1851 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1852 v_i(ji,jj,jl) = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1853 v_s(ji,jj,jl) = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) )1854 t_su(ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj)1855 oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1856 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) )1816 !----------------------------------------------------------------- 1817 ! Zap snow energy and use ocean heat to melt snow 1818 !----------------------------------------------------------------- 1819 1820 ! xtmp = esnon(i,j,n) / dt ! < 0 1821 ! fhnet(i,j) = fhnet(i,j) + xtmp 1822 ! fhnet_hist(i,j) = fhnet_hist(i,j) + xtmp 1823 ! xtmp is greater than 0 1824 ! fluxes are positive to the ocean 1825 ! here the flux has to be negative for the ocean 1826 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice 1827 ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp 1828 1829 xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB ??????? 1830 1831 t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) ) 1832 1833 !----------------------------------------------------------------- 1834 ! zap ice and snow volume, add water and salt to ocean 1835 !----------------------------------------------------------------- 1836 1837 ! xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 1838 ! fresh(i,j) = fresh(i,j) + xtmp 1839 ! fresh_hist(i,j) = fresh_hist(i,j) + xtmp 1840 1841 ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) ) * & 1842 ! rhosn * v_s(ji,jj,jl) / rdt_ice 1843 1844 ! fsalt_res(ji,jj) = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * & 1845 ! rhoic * v_i(ji,jj,jl) / rdt_ice 1846 1847 ! emps(i,j) = emps(i,j) + xtmp 1848 ! fsalt_hist(i,j) = fsalt_hist(i,j) + xtmp 1849 1850 ato_i(ji,jj) = a_i(ji,jj,jl) * zmask(ji,jj) + ato_i(ji,jj) 1851 a_i(ji,jj,jl) = a_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1852 v_i(ji,jj,jl) = v_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1853 v_s(ji,jj,jl) = v_s(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1854 t_su(ji,jj,jl) = t_su(ji,jj,jl) * (1 -zmask(ji,jj) ) + t_bo(ji,jj) * zmask(ji,jj) 1855 oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1856 smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * ( 1 - zmask(ji,jj) ) 1857 1857 1858 1858 END DO ! ji
Note: See TracChangeset
for help on using the changeset viewer.
