- Timestamp:
- 2020-09-15T12:49:18+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/temporary_r4_trunk/src/ICE/icedyn_adv_pra.F90
r13466 r13469 110 110 END WHERE 111 111 DO jl = 1, jpl 112 DO jj = 2, jpjm1 113 DO ji = fs_2, fs_jpim1 114 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 115 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 116 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 117 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 118 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 119 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 120 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 121 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 122 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 123 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 124 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 125 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 126 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 127 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 128 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 129 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 130 END DO 131 END DO 112 DO_2D_00_00 113 zhip_max(ji,jj,jl) = MAX( epsi20, ph_ip(ji,jj,jl), ph_ip(ji+1,jj ,jl), ph_ip(ji ,jj+1,jl), & 114 & ph_ip(ji-1,jj ,jl), ph_ip(ji ,jj-1,jl), & 115 & ph_ip(ji+1,jj+1,jl), ph_ip(ji-1,jj-1,jl), & 116 & ph_ip(ji+1,jj-1,jl), ph_ip(ji-1,jj+1,jl) ) 117 zhi_max (ji,jj,jl) = MAX( epsi20, ph_i (ji,jj,jl), ph_i (ji+1,jj ,jl), ph_i (ji ,jj+1,jl), & 118 & ph_i (ji-1,jj ,jl), ph_i (ji ,jj-1,jl), & 119 & ph_i (ji+1,jj+1,jl), ph_i (ji-1,jj-1,jl), & 120 & ph_i (ji+1,jj-1,jl), ph_i (ji-1,jj+1,jl) ) 121 zhs_max (ji,jj,jl) = MAX( epsi20, ph_s (ji,jj,jl), ph_s (ji+1,jj ,jl), ph_s (ji ,jj+1,jl), & 122 & ph_s (ji-1,jj ,jl), ph_s (ji ,jj-1,jl), & 123 & ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 124 & ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 125 zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj ,jl), zs_i (ji ,jj+1,jl), & 126 & zs_i (ji-1,jj ,jl), zs_i (ji ,jj-1,jl), & 127 & zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 128 & zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 129 END_2D 132 130 END DO 133 131 CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1., zhs_max, 'T', 1., zhip_max, 'T', 1., zsi_max, 'T', 1. ) … … 145 143 END DO 146 144 DO jl = 1, jpl 147 DO jk = 1, nlay_i 148 DO jj = 2, jpjm1 149 DO ji = fs_2, fs_jpim1 150 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), & 151 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 152 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 153 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 154 END DO 155 END DO 156 END DO 145 DO_3D_00_00( 1, nlay_i ) 146 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), & 147 & ze_i(ji-1,jj ,jk,jl), ze_i(ji ,jj-1,jk,jl), & 148 & ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 149 & ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 150 END_3D 157 151 END DO 158 152 DO jl = 1, jpl 159 DO jk = 1, nlay_s 160 DO jj = 2, jpjm1 161 DO ji = fs_2, fs_jpim1 162 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), & 163 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 164 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 165 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 166 END DO 167 END DO 168 END DO 153 DO_3D_00_00( 1, nlay_s ) 154 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), & 155 & ze_s(ji-1,jj ,jk,jl), ze_s(ji ,jj-1,jk,jl), & 156 & ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 157 & ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 158 END_3D 169 159 END DO 170 160 CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) … … 317 307 ! derive open water from ice concentration 318 308 zati2(:,:) = SUM( pa_i(:,:,:), dim=3 ) 319 DO jj = 2, jpjm1 320 DO ji = fs_2, fs_jpim1 321 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 322 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 323 END DO 324 END DO 309 DO_2D_00_00 310 pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) & !--- open water 311 & - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt 312 END_2D 325 313 CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T', 1. ) 326 314 ! … … 375 363 ! 376 364 ! Limitation of moments. 377 DO jj = 2, jpjm1 378 DO ji = 1, jpi 379 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 380 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 381 ! 382 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 383 zs1max = 1.5 * zslpmax 384 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 385 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 386 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 387 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 388 389 ps0 (ji,jj,jl) = zslpmax 390 psx (ji,jj,jl) = zs1new * rswitch 391 psxx(ji,jj,jl) = zs2new * rswitch 392 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 393 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 394 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 395 END DO 396 END DO 365 DO_2D_00_11 366 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 367 psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 368 ! 369 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 370 zs1max = 1.5 * zslpmax 371 zs1new = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) ) 372 zs2new = MIN( 2.0 * zslpmax - 0.3334 * ABS( zs1new ), & 373 & MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) ) ) 374 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 375 376 ps0 (ji,jj,jl) = zslpmax 377 psx (ji,jj,jl) = zs1new * rswitch 378 psxx(ji,jj,jl) = zs2new * rswitch 379 psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch 380 psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch 381 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 382 END_2D 397 383 398 384 ! Calculate fluxes and moments between boxes i<-->i+1 399 DO jj = 2, jpjm1 ! Flux from i to i+1 WHEN u GT 0 400 DO ji = 1, jpi 401 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 402 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 403 zalfq = zalf * zalf 404 zalf1 = 1.0 - zalf 405 zalf1q = zalf1 * zalf1 406 ! 407 zfm (ji,jj) = zalf * psm (ji,jj,jl) 408 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 409 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 410 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 411 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 412 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 413 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 414 415 ! Readjust moments remaining in the box. 416 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 417 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 418 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 419 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 420 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 421 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 422 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 423 END DO 424 END DO 425 426 DO jj = 2, jpjm1 ! Flux from i+1 to i when u LT 0. 427 DO ji = 1, fs_jpim1 428 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 429 zalg (ji,jj) = zalf 430 zalfq = zalf * zalf 431 zalf1 = 1.0 - zalf 432 zalg1 (ji,jj) = zalf1 433 zalf1q = zalf1 * zalf1 434 zalg1q(ji,jj) = zalf1q 435 ! 436 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 437 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 438 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 439 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 440 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 441 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 442 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 443 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 444 END DO 445 END DO 446 447 DO jj = 2, jpjm1 ! Readjust moments remaining in the box. 448 DO ji = fs_2, fs_jpim1 449 zbt = zbet(ji-1,jj) 450 zbt1 = 1.0 - zbet(ji-1,jj) 451 ! 452 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 453 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 454 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 455 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 456 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 457 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 458 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 459 END DO 460 END DO 385 DO_2D_00_11 386 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 387 zalf = MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 388 zalfq = zalf * zalf 389 zalf1 = 1.0 - zalf 390 zalf1q = zalf1 * zalf1 391 ! 392 zfm (ji,jj) = zalf * psm (ji,jj,jl) 393 zf0 (ji,jj) = zalf * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) ) 394 zfx (ji,jj) = zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) ) 395 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) * zalfq 396 zfy (ji,jj) = zalf * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 397 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 398 zfyy(ji,jj) = zalf * psyy(ji,jj,jl) 399 400 ! Readjust moments remaining in the box. 401 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 402 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 403 psx (ji,jj,jl) = zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) ) 404 psxx(ji,jj,jl) = zalf1 * zalf1q * psxx(ji,jj,jl) 405 psy (ji,jj,jl) = psy (ji,jj,jl) - zfy(ji,jj) 406 psyy(ji,jj,jl) = psyy(ji,jj,jl) - zfyy(ji,jj) 407 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 408 END_2D 409 410 DO_2D_00_10 411 zalf = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 412 zalg (ji,jj) = zalf 413 zalfq = zalf * zalf 414 zalf1 = 1.0 - zalf 415 zalg1 (ji,jj) = zalf1 416 zalf1q = zalf1 * zalf1 417 zalg1q(ji,jj) = zalf1q 418 ! 419 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji+1,jj,jl) 420 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji+1,jj,jl) & 421 & - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) ) 422 zfx (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) ) 423 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji+1,jj,jl) * zalfq 424 zfy (ji,jj) = zfy (ji,jj) + zalf * ( psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) ) 425 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji+1,jj,jl) 426 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji+1,jj,jl) 427 END_2D 428 429 DO_2D_00_00 430 zbt = zbet(ji-1,jj) 431 zbt1 = 1.0 - zbet(ji-1,jj) 432 ! 433 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) ) 434 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) ) 435 psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) ) 436 psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl) 437 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) ) 438 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) ) 439 psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl) 440 END_2D 461 441 462 442 ! Put the temporary moments into appropriate neighboring boxes. 463 DO jj = 2, jpjm1 ! Flux from i to i+1 IF u GT 0. 464 DO ji = fs_2, fs_jpim1 465 zbt = zbet(ji-1,jj) 466 zbt1 = 1.0 - zbet(ji-1,jj) 467 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 468 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 469 zalf1 = 1.0 - zalf 470 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 471 ! 472 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 473 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 474 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 475 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 476 & + zbt1 * psxx(ji,jj,jl) 477 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 478 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 479 & + zbt1 * psxy(ji,jj,jl) 480 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 481 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 482 END DO 483 END DO 484 485 DO jj = 2, jpjm1 ! Flux from i+1 to i IF u LT 0. 486 DO ji = fs_2, fs_jpim1 487 zbt = zbet(ji,jj) 488 zbt1 = 1.0 - zbet(ji,jj) 489 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 490 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 491 zalf1 = 1.0 - zalf 492 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 493 ! 494 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 495 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 496 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 497 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 498 & + ( zalf1 - zalf ) * ztemp ) ) 499 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 500 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 501 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 502 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 503 END DO 504 END DO 443 DO_2D_00_00 444 zbt = zbet(ji-1,jj) 445 zbt1 = 1.0 - zbet(ji-1,jj) 446 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl) 447 zalf = zbt * zfm(ji-1,jj) / psm(ji,jj,jl) 448 zalf1 = 1.0 - zalf 449 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj) 450 ! 451 ps0 (ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl) 452 psx (ji,jj,jl) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl) 453 psxx(ji,jj,jl) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 454 & + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp ) ) & 455 & + zbt1 * psxx(ji,jj,jl) 456 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl) & 457 & + 3.0 * (- zalf1*zfy(ji-1,jj) + zalf * psy(ji,jj,jl) ) ) & 458 & + zbt1 * psxy(ji,jj,jl) 459 psy (ji,jj,jl) = zbt * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl) 460 psyy(ji,jj,jl) = zbt * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl) 461 END_2D 462 463 DO_2D_00_00 464 zbt = zbet(ji,jj) 465 zbt1 = 1.0 - zbet(ji,jj) 466 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 467 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 468 zalf1 = 1.0 - zalf 469 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 470 ! 471 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 472 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) 473 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) & 474 & + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) ) & 475 & + ( zalf1 - zalf ) * ztemp ) ) 476 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 477 & + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) ) 478 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) ) 479 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) ) 480 END_2D 505 481 506 482 END DO … … 544 520 ! 545 521 ! Limitation of moments. 546 DO jj = 1, jpj 547 DO ji = fs_2, fs_jpim1 548 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 549 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 550 ! 551 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 552 zs1max = 1.5 * zslpmax 553 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 554 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 555 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 556 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 557 ! 558 ps0 (ji,jj,jl) = zslpmax 559 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 560 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 561 psy (ji,jj,jl) = zs1new * rswitch 562 psyy(ji,jj,jl) = zs2new * rswitch 563 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 564 END DO 565 END DO 522 DO_2D_11_00 523 ! Initialize volumes of boxes (=area if adv_x first called, =psm otherwise) 524 psm(ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 ) 525 ! 526 zslpmax = MAX( 0._wp, ps0(ji,jj,jl) ) 527 zs1max = 1.5 * zslpmax 528 zs1new = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) ) 529 zs2new = MIN( ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ), & 530 & MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) ) ) 531 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1) ! Case of empty boxes & Apply mask 532 ! 533 ps0 (ji,jj,jl) = zslpmax 534 psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch 535 psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch 536 psy (ji,jj,jl) = zs1new * rswitch 537 psyy(ji,jj,jl) = zs2new * rswitch 538 psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch 539 END_2D 566 540 567 541 ! Calculate fluxes and moments between boxes j<-->j+1 568 DO jj = 1, jpj ! Flux from j to j+1 WHEN v GT 0 569 DO ji = fs_2, fs_jpim1 570 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 571 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 572 zalfq = zalf * zalf 573 zalf1 = 1.0 - zalf 574 zalf1q = zalf1 * zalf1 575 ! 576 zfm (ji,jj) = zalf * psm(ji,jj,jl) 577 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 578 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 579 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 580 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 581 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 582 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 583 ! 584 ! Readjust moments remaining in the box. 585 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 586 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 587 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 588 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 589 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 590 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 591 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 592 END DO 593 END DO 594 ! 595 DO jj = 1, jpjm1 ! Flux from j+1 to j when v LT 0. 596 DO ji = fs_2, fs_jpim1 597 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 598 zalg (ji,jj) = zalf 599 zalfq = zalf * zalf 600 zalf1 = 1.0 - zalf 601 zalg1 (ji,jj) = zalf1 602 zalf1q = zalf1 * zalf1 603 zalg1q(ji,jj) = zalf1q 604 ! 605 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 606 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 607 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 608 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 609 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 610 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 611 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 612 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 613 END DO 614 END DO 542 DO_2D_11_00 543 zbet(ji,jj) = MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 544 zalf = MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 545 zalfq = zalf * zalf 546 zalf1 = 1.0 - zalf 547 zalf1q = zalf1 * zalf1 548 ! 549 zfm (ji,jj) = zalf * psm(ji,jj,jl) 550 zf0 (ji,jj) = zalf * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl) + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 551 zfy (ji,jj) = zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) ) 552 zfyy(ji,jj) = zalf * zalfq * psyy(ji,jj,jl) 553 zfx (ji,jj) = zalf * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) ) 554 zfxy(ji,jj) = zalfq * psxy(ji,jj,jl) 555 zfxx(ji,jj) = zalf * psxx(ji,jj,jl) 556 ! 557 ! Readjust moments remaining in the box. 558 psm (ji,jj,jl) = psm (ji,jj,jl) - zfm(ji,jj) 559 ps0 (ji,jj,jl) = ps0 (ji,jj,jl) - zf0(ji,jj) 560 psy (ji,jj,jl) = zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) ) 561 psyy(ji,jj,jl) = zalf1 * zalf1q * psyy(ji,jj,jl) 562 psx (ji,jj,jl) = psx (ji,jj,jl) - zfx(ji,jj) 563 psxx(ji,jj,jl) = psxx(ji,jj,jl) - zfxx(ji,jj) 564 psxy(ji,jj,jl) = zalf1q * psxy(ji,jj,jl) 565 END_2D 566 ! 567 DO_2D_10_00 568 zalf = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 569 zalg (ji,jj) = zalf 570 zalfq = zalf * zalf 571 zalf1 = 1.0 - zalf 572 zalg1 (ji,jj) = zalf1 573 zalf1q = zalf1 * zalf1 574 zalg1q(ji,jj) = zalf1q 575 ! 576 zfm (ji,jj) = zfm (ji,jj) + zalf * psm (ji,jj+1,jl) 577 zf0 (ji,jj) = zf0 (ji,jj) + zalf * ( ps0 (ji,jj+1,jl) & 578 & - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) ) 579 zfy (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) ) 580 zfyy (ji,jj) = zfyy(ji,jj) + zalf * psyy(ji,jj+1,jl) * zalfq 581 zfx (ji,jj) = zfx (ji,jj) + zalf * ( psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) ) 582 zfxy (ji,jj) = zfxy(ji,jj) + zalfq * psxy(ji,jj+1,jl) 583 zfxx (ji,jj) = zfxx(ji,jj) + zalf * psxx(ji,jj+1,jl) 584 END_2D 615 585 616 586 ! Readjust moments remaining in the box. 617 DO jj = 2, jpjm1 618 DO ji = fs_2, fs_jpim1 619 zbt = zbet(ji,jj-1) 620 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 621 ! 622 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 623 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 624 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 625 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 626 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 627 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 628 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 629 END DO 630 END DO 587 DO_2D_00_00 588 zbt = zbet(ji,jj-1) 589 zbt1 = ( 1.0 - zbet(ji,jj-1) ) 590 ! 591 psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) ) 592 ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) ) 593 psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) ) 594 psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl) 595 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) ) 596 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) ) 597 psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl) 598 END_2D 631 599 632 600 ! Put the temporary moments into appropriate neighboring boxes. 633 DO jj = 2, jpjm1 ! Flux from j to j+1 IF v GT 0. 634 DO ji = fs_2, fs_jpim1 635 zbt = zbet(ji,jj-1) 636 zbt1 = 1.0 - zbet(ji,jj-1) 637 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 638 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 639 zalf1 = 1.0 - zalf 640 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 641 ! 642 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 643 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 644 & + zbt1 * psy(ji,jj,jl) 645 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 646 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 647 & + zbt1 * psyy(ji,jj,jl) 648 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 649 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 650 & + zbt1 * psxy(ji,jj,jl) 651 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 652 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 653 END DO 654 END DO 655 656 DO jj = 2, jpjm1 ! Flux from j+1 to j IF v LT 0. 657 DO ji = fs_2, fs_jpim1 658 zbt = zbet(ji,jj) 659 zbt1 = 1.0 - zbet(ji,jj) 660 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 661 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 662 zalf1 = 1.0 - zalf 663 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 664 ! 665 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 666 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 667 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 668 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 669 & + ( zalf1 - zalf ) * ztemp ) ) 670 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 671 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 672 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 673 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 674 END DO 675 END DO 601 DO_2D_00_00 602 zbt = zbet(ji,jj-1) 603 zbt1 = 1.0 - zbet(ji,jj-1) 604 psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 605 zalf = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 606 zalf1 = 1.0 - zalf 607 ztemp = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1) 608 ! 609 ps0(ji,jj,jl) = zbt * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl) 610 psy(ji,jj,jl) = zbt * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) & 611 & + zbt1 * psy(ji,jj,jl) 612 psyy(ji,jj,jl) = zbt * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl) & 613 & + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 614 & + zbt1 * psyy(ji,jj,jl) 615 psxy(ji,jj,jl) = zbt * ( zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl) & 616 & + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) ) & 617 & + zbt1 * psxy(ji,jj,jl) 618 psx (ji,jj,jl) = zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl) 619 psxx(ji,jj,jl) = zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl) 620 END_2D 621 622 DO_2D_00_00 623 zbt = zbet(ji,jj) 624 zbt1 = 1.0 - zbet(ji,jj) 625 psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) ) 626 zalf = zbt1 * zfm(ji,jj) / psm(ji,jj,jl) 627 zalf1 = 1.0 - zalf 628 ztemp = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj) 629 ! 630 ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) ) 631 psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp ) 632 psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) & 633 & + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) ) & 634 & + ( zalf1 - zalf ) * ztemp ) ) 635 psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl) & 636 & + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) ) 637 psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) ) 638 psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) ) 639 END_2D 676 640 677 641 END DO … … 715 679 ! 716 680 DO jl = 1, jpl 717 DO jj = 1, jpj 718 DO ji = 1, jpi 719 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 681 DO_2D_11_11 682 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 683 ! 684 ! ! -- check h_ip -- ! 685 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 686 IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 687 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 688 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 689 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 690 ENDIF 691 ENDIF 692 ! 693 ! ! -- check h_i -- ! 694 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 695 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 696 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 697 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 698 ENDIF 699 ! 700 ! ! -- check h_s -- ! 701 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 702 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 703 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 704 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 720 705 ! 721 ! ! -- check h_ip -- ! 722 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 723 IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 724 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 725 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 726 pa_ip(ji,jj,jl) = pv_ip(ji,jj,jl) / phip_max(ji,jj,jl) 727 ENDIF 728 ENDIF 706 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 707 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 729 708 ! 730 ! ! -- check h_i -- ! 731 ! if h_i is larger than the surrounding 9 pts => reduce h_i and increase a_i 732 zhi = pv_i(ji,jj,jl) / pa_i(ji,jj,jl) 733 IF( zhi > phi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 734 pa_i(ji,jj,jl) = pv_i(ji,jj,jl) / MIN( phi_max(ji,jj,jl), hi_max(jpl) ) !-- bound h_i to hi_max (99 m) 735 ENDIF 736 ! 737 ! ! -- check h_s -- ! 738 ! if h_s is larger than the surrounding 9 pts => put the snow excess in the ocean 739 zhs = pv_s(ji,jj,jl) / pa_i(ji,jj,jl) 740 IF( pv_s(ji,jj,jl) > 0._wp .AND. zhs > phs_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 741 zfra = phs_max(ji,jj,jl) / MAX( zhs, epsi20 ) 742 ! 743 wfx_res(ji,jj) = wfx_res(ji,jj) + ( pv_s(ji,jj,jl) - pa_i(ji,jj,jl) * phs_max(ji,jj,jl) ) * rhos * z1_dt 744 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 745 ! 746 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 747 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 748 ENDIF 749 ! 750 ! ! -- check s_i -- ! 751 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 752 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 753 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 754 zfra = psi_max(ji,jj,jl) / zsi 755 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 756 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 757 ENDIF 758 ! 759 ENDIF 760 END DO 761 END DO 709 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 710 pv_s(ji,jj,jl) = pa_i(ji,jj,jl) * phs_max(ji,jj,jl) 711 ENDIF 712 ! 713 ! ! -- check s_i -- ! 714 ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 715 zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 716 IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 717 zfra = psi_max(ji,jj,jl) / zsi 718 sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 719 psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 720 ENDIF 721 ! 722 ENDIF 723 END_2D 762 724 END DO 763 725 ! 764 726 ! ! -- check e_i/v_i -- ! 765 727 DO jl = 1, jpl 766 DO jk = 1, nlay_i 767 DO jj = 1, jpj 768 DO ji = 1, jpi 769 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 770 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 771 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 772 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 773 zfra = pei_max(ji,jj,jk,jl) / zei 774 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 775 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 776 ENDIF 777 ENDIF 778 END DO 779 END DO 780 END DO 728 DO_3D_11_11( 1, nlay_i ) 729 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 730 ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 731 zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 732 IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 733 zfra = pei_max(ji,jj,jk,jl) / zei 734 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 735 pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 736 ENDIF 737 ENDIF 738 END_3D 781 739 END DO 782 740 ! ! -- check e_s/v_s -- ! 783 741 DO jl = 1, jpl 784 DO jk = 1, nlay_s 785 DO jj = 1, jpj 786 DO ji = 1, jpi 787 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 788 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 789 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 790 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 791 zfra = pes_max(ji,jj,jk,jl) / zes 792 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 793 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 794 ENDIF 795 ENDIF 796 END DO 797 END DO 798 END DO 742 DO_3D_11_11( 1, nlay_s ) 743 IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 744 ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 745 zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 746 IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 747 zfra = pes_max(ji,jj,jk,jl) / zes 748 hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 749 pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 750 ENDIF 751 ENDIF 752 END_3D 799 753 END DO 800 754 ! … … 829 783 ! -- check snow load -- ! 830 784 DO jl = 1, jpl 831 DO jj = 1, jpj 832 DO ji = 1, jpi 833 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 834 ! 835 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 836 ! 837 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 838 ! put snow excess in the ocean 839 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 840 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 841 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 842 ! correct snow volume and heat content 843 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 844 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 845 ENDIF 846 ! 847 ENDIF 848 END DO 849 END DO 785 DO_2D_11_11 786 IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 787 ! 788 zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos ) 789 ! 790 IF( zvs_excess > 0._wp ) THEN ! snow-ice interface deplets below the ocean surface 791 ! put snow excess in the ocean 792 zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 ) 793 wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt 794 hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 795 ! correct snow volume and heat content 796 pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra 797 pv_s(ji,jj,jl) = pv_s(ji,jj,jl) - zvs_excess 798 ENDIF 799 ! 800 ENDIF 801 END_2D 850 802 END DO 851 803 !
Note: See TracChangeset
for help on using the changeset viewer.
