Changeset 13618
- Timestamp:
- 2020-10-16T10:07:55+02:00 (4 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 4 edited
-
icedyn_adv_pra.F90 (modified) (9 diffs)
-
icedyn_adv_umx.F90 (modified) (39 diffs)
-
icedyn_rdgrft.F90 (modified) (1 diff)
-
iceitd.F90 (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/icedyn_adv_pra.F90
r13613 r13618 111 111 ELSEWHERE ; zs_i(:,:,:) = 0._wp 112 112 END WHERE 113 DO jl = 1, jpl 114 DO_2D( 0, 0, 0, 0 ) 115 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 116 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 117 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 118 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 119 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 120 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 121 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 122 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 123 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 124 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 125 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 126 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 127 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 128 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 129 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 130 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 131 END_2D 132 END DO 113 CALL icemax3D( ph_i , zhi_max ) 114 CALL icemax3D( ph_s , zhs_max ) 115 CALL icemax3D( ph_ip, zhip_max) 116 CALL icemax3D( zs_i , zsi_max ) 133 117 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 134 118 ! … … 143 127 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 144 128 END WHERE 145 END DO 146 DO jl = 1, jpl 147 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 148 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 149 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 150 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 151 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 152 END_3D 153 END DO 154 DO jl = 1, jpl 155 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 156 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 157 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 158 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 159 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 160 END_3D 161 END DO 162 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 163 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 129 END DO 130 CALL icemax4D( ze_i , zei_max ) 131 CALL icemax4D( ze_s , zes_max ) 132 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1._wp ) 133 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1._wp ) 164 134 ! 165 135 ! … … 399 369 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 400 370 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 371 REAL(wp) :: zpsm, zps0 372 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 401 373 REAL(wp), DIMENSION(jpi,jpj) :: zf0 , zfx , zfy , zbet ! 2D workspace 402 374 REAL(wp), DIMENSION(jpi,jpj) :: zfm , zfxx , zfyy , zfxy ! - - … … 417 389 DO ji = Nis0 - 1, Nie0 + 1 418 390 391 zpsm = psm (ji,jj,jl) ! optimization 392 zps0 = ps0 (ji,jj,jl) 393 zpsx = psx (ji,jj,jl) 394 zpsxx = psxx(ji,jj,jl) 395 zpsy = psy (ji,jj,jl) 396 zpsyy = psyy(ji,jj,jl) 397 zpsxy = psxy(ji,jj,jl) 398 419 399 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 420 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl), epsi20 )421 ! 422 zslpmax = MAX( 0._wp, ps0(ji,jj,jl))400 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 401 ! 402 zslpmax = MAX( 0._wp, zps0 ) 423 403 zs1max = 1.5 * zslpmax 424 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl)) )425 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl)) )404 zs1new = MIN( zs1max, MAX( -zs1max, zpsx ) ) 405 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), MAX( ABS( zs1new ) - zslpmax, zpsxx ) ) 426 406 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 427 428 ps0 (ji,jj,jl)= zslpmax429 psx (ji,jj,jl) = zs1new* rswitch430 psxx(ji,jj,jl) = zs2new* rswitch431 psy (ji,jj,jl) = psy (ji,jj,jl)* rswitch432 psyy(ji,jj,jl) = psyy(ji,jj,jl)* rswitch433 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl)) ) * rswitch434 407 408 zps0 = zslpmax 409 zpsx = zs1new * rswitch 410 zpsxx = zs2new * rswitch 411 zpsy = zpsy * rswitch 412 zpsyy = zpsyy * rswitch 413 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 414 435 415 ! Calculate fluxes and moments between boxes i<-->i+1 436 416 ! ! Flux from i to i+1 WHEN u GT 0 437 417 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 438 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)418 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / zpsm 439 419 zalfq = zalf * zalf 440 420 zalf1 = 1.0 - zalf 441 421 zalf1q = zalf1 * zalf1 442 422 ! 443 zfm (ji,jj) = zalf * psm (ji,jj,jl)444 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl)) )445 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl))446 zfxx(ji,jj) = zalf * psxx(ji,jj,jl)* zalfq447 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl))448 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl)449 zfyy(ji,jj) = zalf * psyy(ji,jj,jl)450 423 zfm (ji,jj) = zalf * zpsm 424 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsx + (zalf1 - zalf) * zpsxx ) ) 425 zfx (ji,jj) = zalfq * ( zpsx + 3.0 * zalf1 * zpsxx ) 426 zfxx(ji,jj) = zalf * zpsxx * zalfq 427 zfy (ji,jj) = zalf * ( zpsy + zalf1 * zpsxy ) 428 zfxy(ji,jj) = zalfq * zpsxy 429 zfyy(ji,jj) = zalf * zpsyy 430 451 431 ! ! Readjust moments remaining in the box. 452 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 453 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 454 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 455 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 456 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 457 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 458 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 432 zpsm = zpsm - zfm(ji,jj) 433 zps0 = zps0 - zf0(ji,jj) 434 zpsx = zalf1q * ( zpsx - 3.0 * zalf * zpsxx ) 435 zpsxx = zalf1 * zalf1q * zpsxx 436 zpsy = zpsy - zfy (ji,jj) 437 zpsyy = zpsyy - zfyy(ji,jj) 438 zpsxy = zalf1q * zpsxy 439 ! 440 psm (ji,jj,jl) = zpsm ! optimization 441 ps0 (ji,jj,jl) = zps0 442 psx (ji,jj,jl) = zpsx 443 psxx(ji,jj,jl) = zpsxx 444 psy (ji,jj,jl) = zpsy 445 psyy(ji,jj,jl) = zpsyy 446 psxy(ji,jj,jl) = zpsxy 447 ! 459 448 END DO 460 449 … … 480 469 481 470 DO ji = Nis0, Nie0 471 ! 472 zpsm = psm (ji,jj,jl) ! optimization 473 zps0 = ps0 (ji,jj,jl) 474 zpsx = psx (ji,jj,jl) 475 zpsxx = psxx(ji,jj,jl) 476 zpsy = psy (ji,jj,jl) 477 zpsyy = psyy(ji,jj,jl) 478 zpsxy = psxy(ji,jj,jl) 482 479 ! ! Readjust moments remaining in the box. 483 480 zbt = zbet(ji-1,jj) 484 481 zbt1 = 1.0 - zbet(ji-1,jj) 485 482 ! 486 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl)- zfm(ji-1,jj) )487 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl)- zf0(ji-1,jj) )488 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl))489 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)490 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl)- zfy (ji-1,jj) )491 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl)- zfyy(ji-1,jj) )492 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)493 483 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji-1,jj) ) 484 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji-1,jj) ) 485 zpsx = zalg1q(ji-1,jj) * ( zpsx + 3.0 * zalg(ji-1,jj) * zpsxx ) 486 zpsxx = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * zpsxx 487 zpsy = zbt * zpsy + zbt1 * ( zpsy - zfy (ji-1,jj) ) 488 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy - zfyy(ji-1,jj) ) 489 zpsxy = zalg1q(ji-1,jj) * zpsxy 490 494 491 ! Put the temporary moments into appropriate neighboring boxes. 495 492 ! ! Flux from i to i+1 IF u GT 0. 496 zbt = zbet(ji-1,jj)497 zbt1 = 1.0 - zbet(ji-1,jj)498 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)499 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)500 zalf1 501 ztemp = zalf * ps0(ji,jj,jl)- zalf1 * zf0(ji-1,jj)502 ! 503 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)504 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)505 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)&506 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )) &507 & + zbt1 * psxx(ji,jj,jl)508 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)&509 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) )&510 & + zbt1 * psxy(ji,jj,jl)511 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)512 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)513 493 zbt = zbet(ji-1,jj) 494 zbt1 = 1.0 - zbet(ji-1,jj) 495 zpsm = zbt * ( zpsm + zfm(ji-1,jj) ) + zbt1 * zpsm 496 zalf = zbt * zfm(ji-1,jj) / zpsm 497 zalf1 = 1.0 - zalf 498 ztemp = zalf * zps0 - zalf1 * zf0(ji-1,jj) 499 ! 500 zps0 = zbt * ( zps0 + zf0(ji-1,jj) ) + zbt1 * zps0 501 zpsx = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * zpsx + 3.0 * ztemp ) + zbt1 * zpsx 502 zpsxx = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * zpsxx & 503 & + 5.0 * ( zalf * zalf1 * ( zpsx - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 504 & + zbt1 * zpsxx 505 zpsxy = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * zpsxy & 506 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * zpsy ) ) & 507 & + zbt1 * zpsxy 508 zpsy = zbt * ( zpsy + zfy (ji-1,jj) ) + zbt1 * zpsy 509 zpsyy = zbt * ( zpsyy + zfyy(ji-1,jj) ) + zbt1 * zpsyy 510 514 511 ! ! Flux from i+1 to i IF u LT 0. 515 zbt = zbet(ji,jj) 516 zbt1 = 1.0 - zbet(ji,jj) 517 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 518 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 519 zalf1 = 1.0 - zalf 520 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 521 ! 522 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 523 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 524 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 525 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 526 & + ( zalf1 - zalf ) * ztemp ) ) 527 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 528 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 529 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 530 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 512 zbt = zbet(ji,jj) 513 zbt1 = 1.0 - zbet(ji,jj) 514 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 515 zalf = zbt1 * zfm(ji,jj) / zpsm 516 zalf1 = 1.0 - zalf 517 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 518 ! 519 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 520 zpsx = zbt * zpsx + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * zpsx + 3.0 * ztemp ) 521 zpsxx = zbt * zpsxx + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * zpsxx & 522 & + 5.0 * ( zalf * zalf1 * ( - zpsx + zfx(ji,jj) ) & 523 & + ( zalf1 - zalf ) * ztemp ) ) 524 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 525 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * zpsy ) ) 526 zpsy = zbt * zpsy + zbt1 * ( zpsy + zfy (ji,jj) ) 527 zpsyy = zbt * zpsyy + zbt1 * ( zpsyy + zfyy(ji,jj) ) 528 ! 529 psm (ji,jj,jl) = zpsm ! optimization 530 ps0 (ji,jj,jl) = zps0 531 psx (ji,jj,jl) = zpsx 532 psxx(ji,jj,jl) = zpsxx 533 psy (ji,jj,jl) = zpsy 534 psyy(ji,jj,jl) = zpsyy 535 psxy(ji,jj,jl) = zpsxy 531 536 END DO 532 537 ! … … 559 564 REAL(wp) :: zs1new, zalf , zalfq , zbt ! - - 560 565 REAL(wp) :: zs2new, zalf1, zalf1q, zbt1 ! - - 566 REAL(wp) :: zpsm, zps0 567 REAL(wp) :: zpsx, zpsy, zpsxx, zpsyy, zpsxy 561 568 REAL(wp), DIMENSION(jpi,jpj) :: zf0, zfx , zfy , zbet ! 2D workspace 562 569 REAL(wp), DIMENSION(jpi,jpj) :: zfm, zfxx, zfyy, zfxy ! - - … … 574 581 ! Limitation of moments. 575 582 DO_2D( 1, 1, ji0, ji0 ) 576 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 577 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 578 ! 579 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 583 ! 584 zpsm = psm (ji,jj,jl) ! optimization 585 zps0 = ps0 (ji,jj,jl) 586 zpsx = psx (ji,jj,jl) 587 zpsxx = psxx(ji,jj,jl) 588 zpsy = psy (ji,jj,jl) 589 zpsyy = psyy(ji,jj,jl) 590 zpsxy = psxy(ji,jj,jl) 591 ! 592 ! Initialize volumes of boxes (=area if adv_y first called, =psm otherwise) 593 zpsm = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * zpsm , epsi20 ) 594 ! 595 zslpmax = MAX( 0._wp, zps0 ) 580 596 zs1max = 1.5 * zslpmax 581 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl)) )582 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl)) )597 zs1new = MIN( zs1max, MAX( -zs1max, zpsy ) ) 598 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), MAX( ABS( zs1new )-zslpmax, zpsyy ) ) 583 599 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 584 600 ! 585 ps0 (ji,jj,jl)= zslpmax586 psx (ji,jj,jl) = psx (ji,jj,jl)* rswitch587 psxx(ji,jj,jl) = psxx(ji,jj,jl)* rswitch588 psy (ji,jj,jl)= zs1new * rswitch589 psyy(ji,jj,jl)= zs2new * rswitch590 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl)) ) * rswitch591 601 zps0 = zslpmax 602 zpsx = zpsx * rswitch 603 zpsxx = zpsxx * rswitch 604 zpsy = zs1new * rswitch 605 zpsyy = zs2new * rswitch 606 zpsxy = MIN( zslpmax, MAX( -zslpmax, zpsxy ) ) * rswitch 607 592 608 ! Calculate fluxes and moments between boxes j<-->j+1 593 609 ! ! Flux from j to j+1 WHEN v GT 0 594 610 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 595 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)611 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / zpsm 596 612 zalfq = zalf * zalf 597 613 zalf1 = 1.0 - zalf 598 614 zalf1q = zalf1 * zalf1 599 615 ! 600 zfm (ji,jj) = zalf * psm(ji,jj,jl)601 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl)) )602 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl))603 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl)604 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl))605 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl)606 zfxx(ji,jj) = zalf * psxx(ji,jj,jl)616 zfm (ji,jj) = zalf * zpsm 617 zf0 (ji,jj) = zalf * ( zps0 + zalf1 * ( zpsy + (zalf1-zalf) * zpsyy ) ) 618 zfy (ji,jj) = zalfq *( zpsy + 3.0*zalf1*zpsyy ) 619 zfyy(ji,jj) = zalf * zalfq * zpsyy 620 zfx (ji,jj) = zalf * ( zpsx + zalf1 * zpsxy ) 621 zfxy(ji,jj) = zalfq * zpsxy 622 zfxx(ji,jj) = zalf * zpsxx 607 623 ! 608 624 ! ! Readjust moments remaining in the box. 609 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 610 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 611 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 612 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 613 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 614 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 615 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 625 zpsm = zpsm - zfm(ji,jj) 626 zps0 = zps0 - zf0(ji,jj) 627 zpsy = zalf1q * ( zpsy -3.0 * zalf * zpsyy ) 628 zpsyy = zalf1 * zalf1q * zpsyy 629 zpsx = zpsx - zfx(ji,jj) 630 zpsxx = zpsxx - zfxx(ji,jj) 631 zpsxy = zalf1q * zpsxy 632 ! 633 psm (ji,jj,jl) = zpsm ! optimization 634 ps0 (ji,jj,jl) = zps0 635 psx (ji,jj,jl) = zpsx 636 psxx(ji,jj,jl) = zpsxx 637 psy (ji,jj,jl) = zpsy 638 psyy(ji,jj,jl) = zpsyy 639 psxy(ji,jj,jl) = zpsxy 616 640 END_2D 617 641 ! … … 641 665 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 642 666 ! 643 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 644 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 645 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 646 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 647 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 648 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 649 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 667 zpsm = psm (ji,jj,jl) ! optimization 668 zps0 = ps0 (ji,jj,jl) 669 zpsx = psx (ji,jj,jl) 670 zpsxx = psxx(ji,jj,jl) 671 zpsy = psy (ji,jj,jl) 672 zpsyy = psyy(ji,jj,jl) 673 zpsxy = psxy(ji,jj,jl) 674 ! 675 zpsm = zbt * zpsm + zbt1 * ( zpsm - zfm(ji,jj-1) ) 676 zps0 = zbt * zps0 + zbt1 * ( zps0 - zf0(ji,jj-1) ) 677 zpsy = zalg1q(ji,jj-1) * ( zpsy + 3.0 * zalg(ji,jj-1) * zpsyy ) 678 zpsyy = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * zpsyy 679 zpsx = zbt * zpsx + zbt1 * ( zpsx - zfx (ji,jj-1) ) 680 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx - zfxx(ji,jj-1) ) 681 zpsxy = zalg1q(ji,jj-1) * zpsxy 650 682 651 683 ! Put the temporary moments into appropriate neighboring boxes. 652 684 ! ! Flux from j to j+1 IF v GT 0. 653 zbt = zbet(ji,jj-1)654 zbt1 = 1.0 - zbet(ji,jj-1)655 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl)656 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl)657 zalf1 658 ztemp = zalf * ps0(ji,jj,jl)- zalf1 * zf0(ji,jj-1)659 ! 660 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)661 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl)+ 3.0 * ztemp ) &662 & + zbt1 * psy(ji,jj,jl)663 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)&664 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl)- zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) &665 & + zbt1 * psyy(ji,jj,jl)666 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)&667 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl)) ) &668 & + zbt1 * psxy(ji,jj,jl)669 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)670 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)685 zbt = zbet(ji,jj-1) 686 zbt1 = 1.0 - zbet(ji,jj-1) 687 zpsm = zbt * ( zpsm + zfm(ji,jj-1) ) + zbt1 * zpsm 688 zalf = zbt * zfm(ji,jj-1) / zpsm 689 zalf1 = 1.0 - zalf 690 ztemp = zalf * zps0 - zalf1 * zf0(ji,jj-1) 691 ! 692 zps0 = zbt * ( zps0 + zf0(ji,jj-1) ) + zbt1 * zps0 693 zpsy = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * zpsy + 3.0 * ztemp ) & 694 & + zbt1 * zpsy 695 zpsyy = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * zpsyy & 696 & + 5.0 * ( zalf * zalf1 * ( zpsy - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 697 & + zbt1 * zpsyy 698 zpsxy = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * zpsxy & 699 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * zpsx ) ) & 700 & + zbt1 * zpsxy 701 zpsx = zbt * ( zpsx + zfx (ji,jj-1) ) + zbt1 * zpsx 702 zpsxx = zbt * ( zpsxx + zfxx(ji,jj-1) ) + zbt1 * zpsxx 671 703 672 704 ! ! Flux from j+1 to j IF v LT 0. 673 zbt = zbet(ji,jj) 674 zbt1 = 1.0 - zbet(ji,jj) 675 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 676 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 677 zalf1 = 1.0 - zalf 678 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 679 ! 680 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 681 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 682 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 683 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 684 & + ( zalf1 - zalf ) * ztemp ) ) 685 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 686 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 687 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 688 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 705 zbt = zbet(ji,jj) 706 zbt1 = 1.0 - zbet(ji,jj) 707 zpsm = zbt * zpsm + zbt1 * ( zpsm + zfm(ji,jj) ) 708 zalf = zbt1 * zfm(ji,jj) / zpsm 709 zalf1 = 1.0 - zalf 710 ztemp = - zalf * zps0 + zalf1 * zf0(ji,jj) 711 ! 712 zps0 = zbt * zps0 + zbt1 * ( zps0 + zf0(ji,jj) ) 713 zpsy = zbt * zpsy + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * zpsy + 3.0 * ztemp ) 714 zpsyy = zbt * zpsyy + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * zpsyy & 715 & + 5.0 * ( zalf * zalf1 * ( - zpsy + zfy(ji,jj) ) & 716 & + ( zalf1 - zalf ) * ztemp ) ) 717 zpsxy = zbt * zpsxy + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * zpsxy & 718 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * zpsx ) ) 719 zpsx = zbt * zpsx + zbt1 * ( zpsx + zfx (ji,jj) ) 720 zpsxx = zbt * zpsxx + zbt1 * ( zpsxx + zfxx(ji,jj) ) 721 ! 722 psm (ji,jj,jl) = zpsm ! optimization 723 ps0 (ji,jj,jl) = zps0 724 psx (ji,jj,jl) = zpsx 725 psxx(ji,jj,jl) = zpsxx 726 psy (ji,jj,jl) = zpsy 727 psyy(ji,jj,jl) = zpsyy 728 psxy(ji,jj,jl) = zpsxy 689 729 END_2D 690 730 ! … … 1101 1141 END SUBROUTINE adv_pra_rst 1102 1142 1143 SUBROUTINE icemax3D( pice , pmax ) 1144 !!--------------------------------------------------------------------- 1145 !! *** ROUTINE icemax3D *** 1146 !! ** Purpose : compute the max of the 9 points around 1147 !!---------------------------------------------------------------------- 1148 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1149 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1150 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1151 INTEGER :: ji, jj, jl ! dummy loop indices 1152 !!---------------------------------------------------------------------- 1153 DO jl = 1, jpl 1154 DO jj = Njs0-1, Nje0+1 1155 DO ji = Nis0, Nie0 1156 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1157 END DO 1158 END DO 1159 DO jj = Njs0, Nje0 1160 DO ji = Nis0, Nie0 1161 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1162 END DO 1163 END DO 1164 END DO 1165 END SUBROUTINE icemax3D 1166 1167 SUBROUTINE icemax4D( pice , pmax ) 1168 !!--------------------------------------------------------------------- 1169 !! *** ROUTINE icemax4D *** 1170 !! ** Purpose : compute the max of the 9 points around 1171 !!---------------------------------------------------------------------- 1172 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1173 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1174 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1175 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1176 !!---------------------------------------------------------------------- 1177 jlay = SIZE( pice , 3 ) ! size of input arrays 1178 DO jl = 1, jpl 1179 DO jk = 1, jlay 1180 DO jj = Njs0-1, Nje0+1 1181 DO ji = Nis0, Nie0 1182 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1183 END DO 1184 END DO 1185 DO jj = Njs0, Nje0 1186 DO ji = Nis0, Nie0 1187 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1188 END DO 1189 END DO 1190 END DO 1191 END DO 1192 END SUBROUTINE icemax4D 1193 1103 1194 #else 1104 1195 !!---------------------------------------------------------------------- -
NEMO/trunk/src/ICE/icedyn_adv_umx.F90
r13609 r13618 115 115 ELSEWHERE ; zs_i(:,:,:) = 0._wp 116 116 END WHERE 117 DO jl = 1, jpl 118 DO_2D( 0, 0, 0, 0 ) 119 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 120 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 121 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 122 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 123 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 124 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 125 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 126 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 127 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 128 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 129 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 130 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 131 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 132 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 133 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 134 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 135 END_2D 136 END DO 117 CALL icemax3D( ph_i , zhi_max ) 118 CALL icemax3D( ph_s , zhs_max ) 119 CALL icemax3D( ph_ip, zhip_max) 120 CALL icemax3D( zs_i , zsi_max ) 137 121 CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 138 122 ! … … 147 131 ELSEWHERE ; ze_s(:,:,jk,:) = 0._wp 148 132 END WHERE 149 END DO 150 DO jl = 1, jpl 151 DO_3D( 0, 0, 0, 0, 1, nlay_i ) 152 zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj ,jk,jl), ze_i(ji ,jj+1,jk,jl), & 153 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 154 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 155 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 156 END_3D 157 END DO 158 DO jl = 1, jpl 159 DO_3D( 0, 0, 0, 0, 1, nlay_s ) 160 zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj ,jk,jl), ze_s(ji ,jj+1,jk,jl), & 161 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 162 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 163 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 164 END_3D 165 END DO 166 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 167 CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 133 END DO 134 CALL icemax4D( ze_i , zei_max ) 135 CALL icemax4D( ze_s , zes_max ) 136 CALL lbc_lnk( 'icedyn_adv_umx', zei_max, 'T', 1._wp ) 137 CALL lbc_lnk( 'icedyn_adv_umx', zes_max, 'T', 1._wp ) 168 138 ! 169 139 ! … … 386 356 ENDIF 387 357 ENDIF 358 359 ! --- Lateral boundary conditions --- ! 360 IF ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 361 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 362 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 363 ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 364 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 365 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 366 ELSE 367 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp ) 368 ENDIF 369 CALL lbc_lnk( 'icedyn_adv_umx', pe_i, 'T', 1._wp ) 370 CALL lbc_lnk( 'icedyn_adv_umx', pe_s, 'T', 1._wp ) 388 371 ! 389 372 !== Open water area ==! … … 393 376 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 394 377 END_2D 395 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1. 0_wp )378 CALL lbc_lnk( 'icedyn_adv_umx', pato_i, 'T', 1._wp ) 396 379 ! 397 380 ! --- diagnostics --- ! … … 462 445 !! work on H (and not V). It is partly related to the multi-category approach 463 446 !! Therefore, after advection we limit the thickness to the largest value of the 9-points around (only if ice 464 !! concentration is small). Since we do not limit S and T, large values can occur at the edge but it does not really matter 465 !! since sv_i and e_i are still good. 447 !! concentration is small). We also limit S and T. 466 448 !!---------------------------------------------------------------------- 467 449 REAL(wp) , INTENT(in ) :: pamsk ! advection of concentration (1) or other tracers (0) … … 507 489 IF( pamsk == 0._wp ) THEN 508 490 DO jl = 1, jpl 509 DO_2D( 1, 0, 1, 0 )491 DO_2D( 0, 0, 1, 0 ) 510 492 IF( ABS( pu(ji,jj) ) > epsi10 ) THEN 511 493 zfu_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) * puc (ji,jj,jl) / pu(ji,jj) … … 516 498 ENDIF 517 499 ! 500 END_2D 501 DO_2D( 1, 0, 0, 0 ) 518 502 IF( ABS( pv(ji,jj) ) > epsi10 ) THEN 519 503 zfv_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) * pvc (ji,jj,jl) / pv(ji,jj) … … 550 534 IF( PRESENT( pua_ho ) ) THEN 551 535 DO jl = 1, jpl 552 DO_2D( 1, 0, 1, 0 ) 553 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) ; pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 554 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) ; pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 536 DO_2D( 0, 0, 1, 0 ) 537 pua_ho (ji,jj,jl) = zfu_ho (ji,jj,jl) 538 pua_ups(ji,jj,jl) = zfu_ups(ji,jj,jl) 539 END_2D 540 DO_2D( 1, 0, 0, 0 ) 541 pva_ho (ji,jj,jl) = zfv_ho (ji,jj,jl) 542 pva_ups(ji,jj,jl) = zfv_ups(ji,jj,jl) 555 543 END_2D 556 544 END DO … … 566 554 END_2D 567 555 END DO 568 CALL lbc_lnk( 'icedyn_adv_umx', ptc, 'T', 1.0_wp )569 556 ! 570 557 END SUBROUTINE adv_umx … … 605 592 ! 606 593 DO jl = 1, jpl !-- flux in x-direction 607 DO_2D( 1, 0, 1, 0 )594 DO_2D( 1, 1, 1, 0 ) 608 595 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * pt(ji+1,jj,jl) 609 596 END_2D … … 611 598 ! 612 599 DO jl = 1, jpl !-- first guess of tracer from u-flux 613 DO_2D( 0, 0, 0, 0 )600 DO_2D( 1, 1, 0, 0 ) 614 601 ztra = - ( pfu_ups(ji,jj,jl) - pfu_ups(ji-1,jj,jl) ) & 615 602 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 618 605 END_2D 619 606 END DO 620 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )621 607 ! 622 608 DO jl = 1, jpl !-- flux in y-direction 623 DO_2D( 1, 0, 1, 0 )609 DO_2D( 1, 0, 0, 0 ) 624 610 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * zpt(ji,jj+1,jl) 625 611 END_2D … … 629 615 ! 630 616 DO jl = 1, jpl !-- flux in y-direction 631 DO_2D( 1, 0, 1, 0)617 DO_2D( 1, 0, 1, 1 ) 632 618 pfv_ups(ji,jj,jl) = MAX( pv(ji,jj), 0._wp ) * pt(ji,jj,jl) + MIN( pv(ji,jj), 0._wp ) * pt(ji,jj+1,jl) 633 619 END_2D … … 635 621 ! 636 622 DO jl = 1, jpl !-- first guess of tracer from v-flux 637 DO_2D( 0, 0, 0, 0)623 DO_2D( 0, 0, 1, 1 ) 638 624 ztra = - ( pfv_ups(ji,jj,jl) - pfv_ups(ji,jj-1,jl) ) & 639 625 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 642 628 END_2D 643 629 END DO 644 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )645 630 ! 646 631 DO jl = 1, jpl !-- flux in x-direction 647 DO_2D( 1, 0, 1, 0 )632 DO_2D( 0, 0, 1, 0 ) 648 633 pfu_ups(ji,jj,jl) = MAX( pu(ji,jj), 0._wp ) * zpt(ji,jj,jl) + MIN( pu(ji,jj), 0._wp ) * zpt(ji+1,jj,jl) 649 634 END_2D … … 694 679 ! 695 680 DO jl = 1, jpl 696 DO_2D( 1, 0, 1, 0 )681 DO_2D( 1, 1, 1, 0 ) 697 682 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj ,jl) ) 683 END_2D 684 DO_2D( 1, 0, 1, 1 ) 698 685 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji ,jj+1,jl) ) 699 686 END_2D … … 712 699 ! 713 700 DO jl = 1, jpl !-- flux in x-direction 714 DO_2D( 1, 0, 1, 0 )701 DO_2D( 1, 1, 1, 0 ) 715 702 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( pt(ji,jj,jl) + pt(ji+1,jj,jl) ) 716 703 END_2D … … 719 706 720 707 DO jl = 1, jpl !-- first guess of tracer from u-flux 721 DO_2D( 0, 0, 0, 0 )708 DO_2D( 1, 1, 0, 0 ) 722 709 ztra = - ( pfu_ho(ji,jj,jl) - pfu_ho(ji-1,jj,jl) ) & 723 710 & + ( pu (ji,jj ) - pu (ji-1,jj ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 726 713 END_2D 727 714 END DO 728 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )729 715 730 716 DO jl = 1, jpl !-- flux in y-direction 731 DO_2D( 1, 0, 1, 0 )717 DO_2D( 1, 0, 0, 0 ) 732 718 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji,jj+1,jl) ) 733 719 END_2D … … 738 724 ! 739 725 DO jl = 1, jpl !-- flux in y-direction 740 DO_2D( 1, 0, 1, 0)726 DO_2D( 1, 0, 1, 1 ) 741 727 pfv_ho(ji,jj,jl) = 0.5_wp * pv(ji,jj) * ( pt(ji,jj,jl) + pt(ji,jj+1,jl) ) 742 728 END_2D … … 745 731 ! 746 732 DO jl = 1, jpl !-- first guess of tracer from v-flux 747 DO_2D( 0, 0, 0, 0)733 DO_2D( 0, 0, 1, 1 ) 748 734 ztra = - ( pfv_ho(ji,jj,jl) - pfv_ho(ji,jj-1,jl) ) & 749 735 & + ( pv (ji,jj ) - pv (ji,jj-1 ) ) * pt(ji,jj,jl) * (1.-pamsk) … … 752 738 END_2D 753 739 END DO 754 CALL lbc_lnk( 'icedyn_adv_umx', zpt, 'T', 1.0_wp )755 740 ! 756 741 DO jl = 1, jpl !-- flux in x-direction 757 DO_2D( 1, 0, 1, 0 )742 DO_2D( 0, 0, 1, 0 ) 758 743 pfu_ho(ji,jj,jl) = 0.5_wp * pu(ji,jj) * ( zpt(ji,jj,jl) + zpt(ji+1,jj,jl) ) 759 744 END_2D … … 912 897 ! 913 898 DO jl = 1, jpl 914 DO_2D( 1, 0, 1, 0 )899 DO_2D( 0, 0, 1, 0 ) 915 900 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & 916 901 & - SIGN( 1._wp, pu(ji,jj) ) * ( pt(ji+1,jj,jl) - pt(ji,jj,jl) ) ) … … 921 906 ! 922 907 DO jl = 1, jpl 923 DO_2D( 1, 0, 1, 0 )908 DO_2D( 0, 0, 1, 0 ) 924 909 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 925 910 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 931 916 ! 932 917 DO jl = 1, jpl 933 DO_2D( 1, 0, 1, 0 )918 DO_2D( 0, 0, 1, 0 ) 934 919 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 935 920 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 945 930 ! 946 931 DO jl = 1, jpl 947 DO_2D( 1, 0, 1, 0 )932 DO_2D( 0, 0, 1, 0 ) 948 933 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 949 934 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 959 944 ! 960 945 DO jl = 1, jpl 961 DO_2D( 1, 0, 1, 0 )946 DO_2D( 0, 0, 1, 0 ) 962 947 zcu = pu(ji,jj) * r1_e2u(ji,jj) * pdt * r1_e1u(ji,jj) 963 948 zdx2 = e1u(ji,jj) * e1u(ji,jj) … … 980 965 IF( ll_neg ) THEN 981 966 DO jl = 1, jpl 982 DO_2D( 1, 0, 1, 0 )967 DO_2D( 0, 0, 1, 0 ) 983 968 IF( pt_u(ji,jj,jl) < 0._wp .OR. ( imsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 984 969 pt_u(ji,jj,jl) = 0.5_wp * umask(ji,jj,1) * ( pt(ji+1,jj,jl) + pt(ji,jj,jl) & … … 1048 1033 CASE( 1 ) !== 1st order central TIM ==! (Eq. 21) 1049 1034 DO jl = 1, jpl 1050 DO_2D( 1, 0, 1, 0 )1035 DO_2D( 1, 0, 0, 0 ) 1051 1036 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & 1052 1037 & - SIGN( 1._wp, pv(ji,jj) ) * ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) ) … … 1056 1041 CASE( 2 ) !== 2nd order central TIM ==! (Eq. 23) 1057 1042 DO jl = 1, jpl 1058 DO_2D( 1, 0, 1, 0 )1043 DO_2D( 1, 0, 0, 0 ) 1059 1044 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1060 1045 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( pt(ji,jj+1,jl) + pt(ji,jj,jl) & … … 1065 1050 CASE( 3 ) !== 3rd order central TIM ==! (Eq. 24) 1066 1051 DO jl = 1, jpl 1067 DO_2D( 1, 0, 1, 0 )1052 DO_2D( 1, 0, 0, 0 ) 1068 1053 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1069 1054 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1078 1063 CASE( 4 ) !== 4th order central TIM ==! (Eq. 27) 1079 1064 DO jl = 1, jpl 1080 DO_2D( 1, 0, 1, 0 )1065 DO_2D( 1, 0, 0, 0 ) 1081 1066 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1082 1067 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1091 1076 CASE( 5 ) !== 5th order central TIM ==! (Eq. 29) 1092 1077 DO jl = 1, jpl 1093 DO_2D( 1, 0, 1, 0 )1078 DO_2D( 1, 0, 0, 0 ) 1094 1079 zcv = pv(ji,jj) * r1_e1v(ji,jj) * pdt * r1_e2v(ji,jj) 1095 1080 zdy2 = e2v(ji,jj) * e2v(ji,jj) … … 1112 1097 IF( ll_neg ) THEN 1113 1098 DO jl = 1, jpl 1114 DO_2D( 1, 0, 1, 0 )1099 DO_2D( 1, 0, 0, 0 ) 1115 1100 IF( pt_v(ji,jj,jl) < 0._wp .OR. ( jmsk_small(ji,jj,jl) == 0 .AND. pamsk == 0. ) ) THEN 1116 1101 pt_v(ji,jj,jl) = 0.5_wp * vmask(ji,jj,1) * ( ( pt(ji,jj+1,jl) + pt(ji,jj,jl) ) & … … 1122 1107 ! !-- High order flux in j-direction --! 1123 1108 DO jl = 1, jpl 1124 DO_2D( 1, 0, 1, 0 )1109 DO_2D( 1, 0, 0, 0 ) 1125 1110 pfv_ho(ji,jj,jl) = pv(ji,jj) * pt_v(ji,jj,jl) 1126 1111 END_2D … … 1158 1143 ! -------------------------------------------------- 1159 1144 DO jl = 1, jpl 1160 DO_2D( 1, 0, 1, 0 )1145 DO_2D( 0, 0, 1, 0 ) 1161 1146 pfu_ho(ji,jj,jl) = pfu_ho(ji,jj,jl) - pfu_ups(ji,jj,jl) 1147 END_2D 1148 DO_2D( 1, 0, 0, 0 ) 1162 1149 pfv_ho(ji,jj,jl) = pfv_ho(ji,jj,jl) - pfv_ups(ji,jj,jl) 1163 1150 END_2D … … 1265 1252 ! --------------------------------- 1266 1253 DO jl = 1, jpl 1267 DO_2D( 1, 0, 1, 0 )1254 DO_2D( 0, 0, 1, 0 ) 1268 1255 zau = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji+1,jj,jl) ) 1269 1256 zbu = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji+1,jj,jl) ) … … 1276 1263 END_2D 1277 1264 1278 DO_2D( 1, 0, 1, 0 )1265 DO_2D( 1, 0, 0, 0 ) 1279 1266 zav = MIN( 1._wp , zbetdo(ji,jj,jl) , zbetup(ji,jj+1,jl) ) 1280 1267 zbv = MIN( 1._wp , zbetup(ji,jj,jl) , zbetdo(ji,jj+1,jl) ) … … 1633 1620 END SUBROUTINE Hsnow 1634 1621 1622 SUBROUTINE icemax3D( pice , pmax ) 1623 !!--------------------------------------------------------------------- 1624 !! *** ROUTINE icemax3D *** 1625 !! ** Purpose : compute the max of the 9 points around 1626 !!---------------------------------------------------------------------- 1627 REAL(wp), DIMENSION(:,:,:) , INTENT(in ) :: pice ! input 1628 REAL(wp), DIMENSION(:,:,:) , INTENT(out) :: pmax ! output 1629 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1630 INTEGER :: ji, jj, jl ! dummy loop indices 1631 !!---------------------------------------------------------------------- 1632 DO jl = 1, jpl 1633 DO jj = Njs0-1, Nje0+1 1634 DO ji = Nis0, Nie0 1635 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jl), pice(ji-1,jj,jl), pice(ji+1,jj,jl) ) 1636 END DO 1637 END DO 1638 DO jj = Njs0, Nje0 1639 DO ji = Nis0, Nie0 1640 pmax(ji,jj,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1641 END DO 1642 END DO 1643 END DO 1644 END SUBROUTINE icemax3D 1645 1646 SUBROUTINE icemax4D( pice , pmax ) 1647 !!--------------------------------------------------------------------- 1648 !! *** ROUTINE icemax4D *** 1649 !! ** Purpose : compute the max of the 9 points around 1650 !!---------------------------------------------------------------------- 1651 REAL(wp), DIMENSION(:,:,:,:) , INTENT(in ) :: pice ! input 1652 REAL(wp), DIMENSION(:,:,:,:) , INTENT(out) :: pmax ! output 1653 REAL(wp), DIMENSION(2:jpim1,jpj) :: zmax ! temporary array 1654 INTEGER :: jlay, ji, jj, jk, jl ! dummy loop indices 1655 !!---------------------------------------------------------------------- 1656 jlay = SIZE( pice , 3 ) ! size of input arrays 1657 DO jl = 1, jpl 1658 DO jk = 1, jlay 1659 DO jj = Njs0-1, Nje0+1 1660 DO ji = Nis0, Nie0 1661 zmax(ji,jj) = MAX( epsi20, pice(ji,jj,jk,jl), pice(ji-1,jj,jk,jl), pice(ji+1,jj,jk,jl) ) 1662 END DO 1663 END DO 1664 DO jj = Njs0, Nje0 1665 DO ji = Nis0, Nie0 1666 pmax(ji,jj,jk,jl) = MAX( epsi20, zmax(ji,jj), zmax(ji,jj-1), zmax(ji,jj+1) ) 1667 END DO 1668 END DO 1669 END DO 1670 END DO 1671 END SUBROUTINE icemax4D 1635 1672 1636 1673 #else -
NEMO/trunk/src/ICE/icedyn_rdgrft.F90
r13495 r13618 349 349 ELSEIF( zGsum(ji,jl-1) < rn_gstar ) THEN 350 350 apartf(ji,jl) = z1_gstar * ( rn_gstar - zGsum(ji,jl-1) ) * & 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar 351 & ( 2._wp - ( zGsum(ji,jl-1) + rn_gstar ) * z1_gstar ) 352 352 ELSE 353 353 apartf(ji,jl) = 0._wp -
NEMO/trunk/src/ICE/iceitd.F90
r13472 r13618 627 627 END_2D 628 628 ! 629 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 630 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) )631 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) )632 !633 DO ji = 1, npti634 jdonor(ji,jl) = jl635 ! how much of a_i you send in cat sup is somewhat arbitrary636 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 637 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 638 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 639 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 640 !! zdaice(ji,jl) = a_i_1d(ji) 641 !! zdvice(ji,jl) = v_i_1d(ji)642 !!clem: these are from UCL and work ok 643 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp644 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1)) * 0.5_wp645 END DO646 !647 IF( npti > 0 ) THEN629 IF( npti > 0 ) THEN 630 !!clem CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,jl) ) 631 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl) ) 632 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl) ) 633 ! 634 DO ji = 1, npti 635 jdonor(ji,jl) = jl 636 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 646 END DO 647 ! 648 648 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl=>jl+1 649 649 ! Reset shift parameters … … 666 666 END_2D 667 667 ! 668 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok669 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok670 DO ji = 1, npti671 jdonor(ji,jl) = jl + 1672 zdaice(ji,jl) = a_i_1d(ji)673 zdvice(ji,jl) = v_i_1d(ji)674 END DO675 !676 668 IF( npti > 0 ) THEN 669 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d(1:npti), a_i(:,:,jl+1) ) ! jl+1 is ok 670 CALL tab_2d_1d( npti, nptidx(1:npti), v_i_1d(1:npti), v_i(:,:,jl+1) ) ! jl+1 is ok 671 DO ji = 1, npti 672 jdonor(ji,jl) = jl + 1 673 zdaice(ji,jl) = a_i_1d(ji) 674 zdvice(ji,jl) = v_i_1d(ji) 675 END DO 676 ! 677 677 CALL itd_shiftice( jdonor(1:npti,:), zdaice(1:npti,:), zdvice(1:npti,:) ) ! Shift jl+1=>jl 678 678 ! Reset shift parameters
Note: See TracChangeset
for help on using the changeset viewer.
