- Timestamp:
- 2013-09-25T16:38:24+02:00 (11 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
r3625 r4045 28 28 USE prtctl ! Print control 29 29 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 30 USE limvar ! clem for ice thickness correction 30 31 31 32 IMPLICIT NONE … … 36 37 REAL(wp) :: epsi06 = 1.e-06_wp ! constant values 37 38 REAL(wp) :: epsi03 = 1.e-03_wp 38 REAL(wp) :: zeps10 = 1.e-10_wp39 REAL(wp) :: epsi10 = 1.e-10_wp 39 40 REAL(wp) :: epsi16 = 1.e-16_wp 41 REAL(wp) :: epsi20 = 1.e-20_wp 40 42 REAL(wp) :: rzero = 0._wp 41 43 REAL(wp) :: rone = 1._wp … … 46 48 # include "vectopt_loop_substitute.h90" 47 49 !!---------------------------------------------------------------------- 48 !! NEMO/LIM3 3.4, UCL - NEMO Consortium (2011)50 !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 49 51 !! $Id$ 50 52 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) … … 69 71 INTEGER :: initad ! number of sub-timestep for the advection 70 72 INTEGER :: ierr ! error status 71 REAL(wp) :: zindb , zindsn , zindic ! local scalar73 REAL(wp) :: zindb , zindsn , zindic, zindh, zinda ! local scalar 72 74 REAL(wp) :: zusvosn, zusvoic, zbigval ! - - 73 75 REAL(wp) :: zcfl , zusnit , zrtt ! - - … … 77 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 78 80 REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zs0e 81 REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 82 REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 83 ! mass and salt flux (clem) 84 REAL(wp), POINTER, DIMENSION(:,:,:) :: zviold ! old ice volume... 85 ! correct ice thickness (clem) 86 REAL(wp), POINTER, DIMENSION(:,:,:) :: zaiold, zhimax ! old ice concentration and thickness 87 REAL(wp) :: zdv, zda, zvi, zvs, zsmv 79 88 !!--------------------------------------------------------------------- 80 89 … … 82 91 CALL wrk_alloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 83 92 CALL wrk_alloc( jpi, jpj, jkmax, jpl, zs0e ) 93 94 CALL wrk_alloc( jpi,jpj,jpl,zviold ) ! clem 95 CALL wrk_alloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 96 97 ! ------------------------------- 98 !- check conservation (C Rousset) 99 IF (ln_limdiahsb) THEN 100 zchk_v_i_b = glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 101 zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 102 zchk_fw_b = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 103 zchk_fs_b = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 104 ENDIF 105 !- check conservation (C Rousset) 106 ! ------------------------------- 84 107 85 108 IF( numit == nstart .AND. lwp ) THEN … … 96 119 IF( ln_limdyn ) THEN ! Advection of sea ice properties ! 97 120 ! !-------------------------------------! 98 ! 99 121 ! mass and salt flux init (clem) 122 zviold(:,:,:) = v_i(:,:,:) 123 124 !--- Thickness correction init. (clem) ------------------------------- 125 CALL lim_var_glo2eqv 126 zaiold(:,:,:) = a_i(:,:,:) 127 !--------------------------------------------------------------------- 128 ! Record max of the surrounding ice thicknesses for correction in limupdate 129 ! in case advection creates ice too thick. 130 !--------------------------------------------------------------------- 131 zhimax(:,:,:) = ht_i(:,:,:) 132 DO jl = 1, jpl 133 DO jj = 2, jpjm1 134 DO ji = 2, jpim1 135 zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) ) 136 !zhimax(ji,jj,jl) = ( ht_i(ji ,jj ,jl) * tmask(ji, jj ,1) + ht_i(ji-1,jj-1,jl) * tmask(ji-1,jj-1,1) + ht_i(ji+1,jj+1,jl) * tmask(ji+1,jj+1,1) & 137 ! & + ht_i(ji-1,jj ,jl) * tmask(ji-1,jj ,1) + ht_i(ji ,jj-1,jl) * tmask(ji ,jj-1,1) & 138 ! & + ht_i(ji+1,jj ,jl) * tmask(ji+1,jj ,1) + ht_i(ji ,jj+1,jl) * tmask(ji ,jj+1,1) & 139 ! & + ht_i(ji-1,jj+1,jl) * tmask(ji-1,jj+1,1) + ht_i(ji+1,jj-1,jl) * tmask(ji+1,jj-1,1) ) 140 END DO 141 END DO 142 CALL lbc_lnk(zhimax(:,:,jl),'T',1.) 143 END DO 144 100 145 !------------------------- 101 146 ! transported fields … … 126 171 ! ENDIF 127 172 !!gm end 128 initad = 1 + INT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) )173 initad = 1 + NINT( MAX( rzero, SIGN( rone, zcfl-0.5 ) ) ) 129 174 zusnit = 1.0 / REAL( initad ) 130 175 IF( zcfl > 0.5 .AND. lwp ) & … … 282 327 END DO 283 328 284 !-----------------------------------------285 ! Remultiply everything by ice area286 !-----------------------------------------287 zs0ow(:,:) = MAX( rzero, zs0ow(:,:) * area(:,:) )288 DO jl = 1, jpl289 zs0ice(:,:,jl) = MAX( rzero, zs0ice(:,:,jl) * area(:,:) ) !!bug: est-ce utile290 zs0sn (:,:,jl) = MAX( rzero, zs0sn (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres291 zs0sm (:,:,jl) = MAX( rzero, zs0sm (:,:,jl) * area(:,:) ) !!bug: cf /area juste apres292 zs0oi (:,:,jl) = MAX( rzero, zs0oi (:,:,jl) * area(:,:) )293 zs0a (:,:,jl) = MAX( rzero, zs0a (:,:,jl) * area(:,:) ) !! suppress both change le resultat294 zs0c0 (:,:,jl) = MAX( rzero, zs0c0 (:,:,jl) * area(:,:) )295 DO jk = 1, nlay_i296 zs0e(:,:,jk,jl) = MAX( rzero, zs0e (:,:,jk,jl) * area(:,:) )297 END DO ! jk298 END DO ! jl299 300 329 !------------------------------------------------------------------------------! 301 330 ! 5) Update and limit ice properties after transport … … 305 334 ! 5.1) Recover mean values over the grid squares. 306 335 !-------------------------------------------------- 307 308 DO jl = 1, jpl309 DO jk = 1, nlay_i310 DO jj = 1, jpj311 DO ji = 1, jpi312 zs0e(ji,jj,jk,jl) = MAX( rzero, zs0e(ji,jj,jk,jl) / area(ji,jj) )313 END DO314 END DO315 END DO316 END DO317 318 DO jj = 1, jpj319 DO ji = 1, jpi320 zs0ow(ji,jj) = MAX( rzero, zs0ow (ji,jj) / area(ji,jj) )321 END DO322 END DO323 324 336 zs0at(:,:) = 0._wp 325 337 DO jl = 1, jpl 326 338 DO jj = 1, jpj 327 339 DO ji = 1, jpi 328 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) /area(ji,jj))329 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) /area(ji,jj))330 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) /area(ji,jj))331 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) /area(ji,jj))332 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) /area(ji,jj))333 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) /area(ji,jj))340 zs0sn (ji,jj,jl) = MAX( rzero, zs0sn (ji,jj,jl) ) 341 zs0ice(ji,jj,jl) = MAX( rzero, zs0ice(ji,jj,jl) ) 342 zs0sm (ji,jj,jl) = MAX( rzero, zs0sm (ji,jj,jl) ) 343 zs0oi (ji,jj,jl) = MAX( rzero, zs0oi (ji,jj,jl) ) 344 zs0a (ji,jj,jl) = MAX( rzero, zs0a (ji,jj,jl) ) 345 zs0c0 (ji,jj,jl) = MAX( rzero, zs0c0 (ji,jj,jl) ) 334 346 zs0at (ji,jj) = zs0at(ji,jj) + zs0a(ji,jj,jl) 335 347 END DO … … 342 354 DO jj = 1, jpj 343 355 DO ji = 1, jpi 344 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - zeps10) )356 zindb = MAX( 0._wp , SIGN( 1.0, zs0at(ji,jj) - epsi10) ) 345 357 zs0ow(ji,jj) = ( 1._wp - zindb ) + zindb * MAX( zs0ow(ji,jj), 0._wp ) 346 358 ato_i(ji,jj) = zs0ow(ji,jj) … … 351 363 DO jj = 1, jpj 352 364 DO ji = 1, jpi 353 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - zeps10) ) 365 zvi = zs0ice(ji,jj,jl) 366 zvs = zs0sn(ji,jj,jl) 354 367 ! 355 zs0a(ji,jj,jl) = zindb * MIN( zs0a(ji,jj,jl), 0.99 ) 368 zindb = MAX( 0.0 , SIGN( 1.0, zs0a(ji,jj,jl) - epsi10) ) 369 ! 356 370 v_s(ji,jj,jl) = zindb * zs0sn (ji,jj,jl) 357 371 v_i(ji,jj,jl) = zindb * zs0ice(ji,jj,jl) 358 372 ! 359 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )360 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )373 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 374 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 361 375 zindb = MAX( zindsn, zindic ) 376 ! 362 377 zs0a(ji,jj,jl) = zindb * zs0a(ji,jj,jl) !ice concentration 363 378 a_i (ji,jj,jl) = zs0a(ji,jj,jl) 364 379 v_s (ji,jj,jl) = zindsn * v_s(ji,jj,jl) 365 380 v_i (ji,jj,jl) = zindic * v_i(ji,jj,jl) 381 ! 382 ! Update mass fluxes (clem) 383 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 384 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 385 END DO 386 END DO 387 END DO 388 389 !--- Thickness correction in case too high (clem) -------------------------------------------------------- 390 CALL lim_var_glo2eqv 391 DO jl = 1, jpl 392 DO jj = 1, jpj 393 DO ji = 1, jpi 394 395 IF ( v_i(ji,jj,jl) > 0._wp ) THEN 396 zvi = v_i(ji,jj,jl) 397 zvs = v_s(ji,jj,jl) 398 zdv = v_i(ji,jj,jl) - zviold(ji,jj,jl) 399 !zda = a_i(ji,jj,jl) - zaiold(ji,jj,jl) 400 401 zindh = 1._wp 402 IF ( ( zdv > 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) .AND. SUM( zaiold(ji,jj,1:jpl) ) < 0.80 ) .OR. & 403 & ( zdv < 0.0 .AND. ht_i(ji,jj,jl) > zhimax(ji,jj,jl) ) ) THEN 404 ht_i(ji,jj,jl) = MIN( zhimax(ji,jj,jl), hi_max(jl) ) 405 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 406 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 407 ELSE 408 ht_i(ji,jj,jl) = MAX( MIN( ht_i(ji,jj,jl), hi_max(jl) ), hi_max(jl-1) ) 409 zindh = MAX( rzero, SIGN( rone, ht_i(ji,jj,jl) - epsi10 ) ) 410 a_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) / MAX( ht_i(ji,jj,jl), epsi10 ) 411 ENDIF 412 413 ! small correction due to *zindh for a_i 414 v_i(ji,jj,jl) = zindh * v_i(ji,jj,jl) 415 v_s(ji,jj,jl) = zindh * v_s(ji,jj,jl) 416 417 ! Update mass fluxes 418 rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic 419 rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn 420 421 ENDIF 422 423 diag_trp_vi(ji,jj) = diag_trp_vi(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 424 366 425 END DO 367 426 END DO 368 427 END DO 369 428 429 ! --- 370 430 DO jj = 1, jpj 371 431 DO ji = 1, jpi 372 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) 432 zs0at(ji,jj) = SUM( zs0a(ji,jj,1:jpl) ) ! clem@useless?? 373 433 END DO 374 434 END DO … … 378 438 !---------------------- 379 439 380 zbigval = 1. d+13440 zbigval = 1.e+13 381 441 382 442 DO jl = 1, jpl 383 443 DO jj = 1, jpj 384 444 DO ji = 1, jpi 445 zsmv = zs0sm(ji,jj,jl) 385 446 386 447 ! Switches and dummy variables … … 388 449 zusvoic = 1.0/MAX( v_i(ji,jj,jl) , epsi16 ) 389 450 zrtt = 173.15 * rone 390 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - zeps10 ) )391 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )451 zindsn = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 452 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 392 453 zindb = MAX( zindsn, zindic ) 393 454 … … 395 456 zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj) , & 396 457 & zusvoic * zs0sm(ji,jj,jl) ) , s_i_min ) * v_i(ji,jj,jl) 397 IF( num_sal == 2 ) smv_i(ji,jj,jl) = zindic * zsal + (1.0-zindic) * 0._wp458 IF( num_sal == 2 ) smv_i(ji,jj,jl) = zindic * zsal 398 459 399 460 zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp ) * a_i(ji,jj,jl) … … 402 463 ! Snow heat content 403 464 ze = MIN( MAX( 0.0, zs0c0(ji,jj,jl)*area(ji,jj) ), zbigval ) 404 e_s(ji,jj,1,jl) = zindsn * ze + (1.0 - zindsn) * 0.0 405 465 e_s(ji,jj,1,jl) = zindsn * ze 466 467 ! Update salt fluxes (clem) 468 sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 406 469 END DO !ji 407 470 END DO !jj … … 413 476 DO ji = 1, jpi 414 477 ! Ice heat content 415 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - zeps10 ) )478 zindic = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 416 479 ze = MIN( MAX( 0.0, zs0e(ji,jj,jk,jl)*area(ji,jj) ), zbigval ) 417 e_i(ji,jj,jk,jl) = zindic * ze + ( 1.0 - zindic ) * 0.0480 e_i(ji,jj,jk,jl) = zindic * ze 418 481 END DO !ji 419 482 END DO ! jj 420 483 END DO ! jk 421 484 END DO ! jl 485 486 487 ! --- agglomerate variables (clem) ----------------- 488 vt_i (:,:) = 0._wp 489 vt_s (:,:) = 0._wp 490 at_i (:,:) = 0._wp 491 ! 492 DO jl = 1, jpl 493 DO jj = 1, jpj 494 DO ji = 1, jpi 495 ! 496 vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 497 vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 498 at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 499 ! 500 zinda = MAX( rzero , SIGN( rone , at_i(ji,jj) - epsi16 ) ) 501 icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi16 ) * zinda ! ice thickness 502 END DO 503 END DO 504 END DO 505 ! ------------------------------------------------- 506 507 422 508 423 509 ENDIF … … 454 540 END DO 455 541 ENDIF 542 ! ------------------------------- 543 !- check conservation (C Rousset) 544 IF (ln_limdiahsb) THEN 545 zchk_fs = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 546 zchk_fw = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 547 548 zchk_v_i = ( glob_sum( SUM( v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 549 zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 550 551 zchk_vmin = glob_min(v_i) 552 zchk_amax = glob_max(SUM(a_i,dim=3)) 553 zchk_amin = glob_min(a_i) 554 555 IF(lwp) THEN 556 IF ( ABS( zchk_v_i ) > 1.e-5 ) WRITE(numout,*) 'violation volume [m3/day] (limtrp) = ',(zchk_v_i * rday) 557 IF ( ABS( zchk_smv ) > 1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday) 558 IF ( zchk_vmin < 0. ) WRITE(numout,*) 'violation v_i<0 [mm] (limtrp) = ',(zchk_vmin * 1.e-3) 559 IF ( zchk_amin < 0. ) WRITE(numout,*) 'violation a_i<0 (limtrp) = ',zchk_amin 560 ENDIF 561 ENDIF 562 !- check conservation (C Rousset) 563 ! ------------------------------- 456 564 ! 457 565 CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow ) 458 566 CALL wrk_dealloc( jpi, jpj, jpl, zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi ) 459 567 CALL wrk_dealloc( jpi, jpj, jkmax, jpl, zs0e ) 568 569 CALL wrk_dealloc( jpi,jpj,jpl,zaiold, zhimax ) ! clem 460 570 ! 461 571 END SUBROUTINE lim_trp
Note: See TracChangeset
for help on using the changeset viewer.
