Changeset 7698 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
- Timestamp:
- 2017-02-18T10:02:03+01:00 (7 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90 (modified) (26 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r7646 r7698 135 135 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 136 136 CALL dom_vvl_rst( nit000, 'READ' ) 137 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 137 !$OMP PARALLEL DO schedule(static) private(jj, ji) 138 DO jj = 1, jpj 139 DO ji = 1, jpi 140 e3t_a(ji,jj,jpk) = e3t_0(ji,jj,jpk) ! last level always inside the sea floor set one for all 141 END DO 142 END DO 138 143 ! 139 144 ! !== Set of all other vertical scale factors ==! (now and before) … … 153 158 ! 154 159 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 155 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) ! reference to the ocean surface (used for MLD and light penetration) 156 gdepw_n(:,:,1) = 0.0_wp 157 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) ! reference to a common level z=0 for hpg 158 gdept_b(:,:,1) = 0.5_wp * e3w_b(:,:,1) 159 gdepw_b(:,:,1) = 0.0_wp 160 !$OMP PARALLEL 161 !$OMP DO schedule(static) private(jj,ji) 162 DO jj = 1, jpj 163 DO ji = 1, jpi 164 gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) ! reference to the ocean surface (used for MLD and light penetration) 165 gdepw_n(ji,jj,1) = 0.0_wp 166 gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj) ! reference to a common level z=0 for hpg 167 gdept_b(ji,jj,1) = 0.5_wp * e3w_b(ji,jj,1) 168 gdepw_b(ji,jj,1) = 0.0_wp 169 END DO 170 END DO 160 171 DO jk = 2, jpk ! vertical sum 172 !$OMP DO schedule(static) private(jj,ji,zcoef) 161 173 DO jj = 1,jpj 162 174 DO ji = 1,jpi … … 178 190 ! 179 191 ! !== thickness of the water column !! (ocean portion only) 180 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 181 hu_b(:,:) = e3u_b(:,:,1) * umask(:,:,1) 182 hu_n(:,:) = e3u_n(:,:,1) * umask(:,:,1) 183 hv_b(:,:) = e3v_b(:,:,1) * vmask(:,:,1) 184 hv_n(:,:) = e3v_n(:,:,1) * vmask(:,:,1) 192 !$OMP DO schedule(static) private(jj,ji) 193 DO jj = 1, jpj 194 DO ji = 1, jpi 195 ht_n(ji,jj) = e3t_n(ji,jj,1) * tmask(ji,jj,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 196 hu_b(ji,jj) = e3u_b(ji,jj,1) * umask(ji,jj,1) 197 hu_n(ji,jj) = e3u_n(ji,jj,1) * umask(ji,jj,1) 198 hv_b(ji,jj) = e3v_b(ji,jj,1) * vmask(ji,jj,1) 199 hv_n(ji,jj) = e3v_n(ji,jj,1) * vmask(ji,jj,1) 200 END DO 201 END DO 185 202 DO jk = 2, jpkm1 186 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 187 hu_b(:,:) = hu_b(:,:) + e3u_b(:,:,jk) * umask(:,:,jk) 188 hu_n(:,:) = hu_n(:,:) + e3u_n(:,:,jk) * umask(:,:,jk) 189 hv_b(:,:) = hv_b(:,:) + e3v_b(:,:,jk) * vmask(:,:,jk) 190 hv_n(:,:) = hv_n(:,:) + e3v_n(:,:,jk) * vmask(:,:,jk) 203 !$OMP DO schedule(static) private(jj,ji) 204 DO jj = 1, jpj 205 DO ji = 1, jpi 206 ht_n(ji,jj) = ht_n(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 207 hu_b(ji,jj) = hu_b(ji,jj) + e3u_b(ji,jj,jk) * umask(ji,jj,jk) 208 hu_n(ji,jj) = hu_n(ji,jj) + e3u_n(ji,jj,jk) * umask(ji,jj,jk) 209 hv_b(ji,jj) = hv_b(ji,jj) + e3v_b(ji,jj,jk) * vmask(ji,jj,jk) 210 hv_n(ji,jj) = hv_n(ji,jj) + e3v_n(ji,jj,jk) * vmask(ji,jj,jk) 211 END DO 212 END DO 191 213 END DO 192 214 ! 193 215 ! !== inverse of water column thickness ==! (u- and v- points) 194 r1_hu_b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 195 r1_hu_n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) 196 r1_hv_b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) ) 197 r1_hv_n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) 198 216 !$OMP DO schedule(static) private(jj,ji) 217 DO jj = 1, jpj 218 DO ji = 1, jpi 219 r1_hu_b(ji,jj) = ssumask(ji,jj) / ( hu_b(ji,jj) + 1._wp - ssumask(ji,jj) ) ! _i mask due to ISF 220 r1_hu_n(ji,jj) = ssumask(ji,jj) / ( hu_n(ji,jj) + 1._wp - ssumask(ji,jj) ) 221 r1_hv_b(ji,jj) = ssvmask(ji,jj) / ( hv_b(ji,jj) + 1._wp - ssvmask(ji,jj) ) 222 r1_hv_n(ji,jj) = ssvmask(ji,jj) / ( hv_n(ji,jj) + 1._wp - ssvmask(ji,jj) ) 223 END DO 224 END DO 225 !$OMP END PARALLEL 199 226 ! !== z_tilde coordinate case ==! (Restoring frequencies) 200 227 IF( ln_vvl_ztilde ) THEN … … 202 229 ! ! Values in days provided via the namelist 203 230 ! ! use rsmall to avoid possible division by zero errors with faulty settings 204 frq_rst_e3t(:,:) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 205 frq_rst_hdv(:,:) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 231 !$OMP PARALLEL DO schedule(static) private(jj,ji) 232 DO jj = 1, jpj 233 DO ji = 1, jpi 234 frq_rst_e3t(ji,jj) = 2._wp * rpi / ( MAX( rn_rst_e3t , rsmall ) * 86400.0_wp ) 235 frq_rst_hdv(ji,jj) = 2._wp * rpi / ( MAX( rn_lf_cutoff, rsmall ) * 86400.0_wp ) 236 END DO 237 END DO 206 238 ! 207 239 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 208 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 209 frq_rst_hdv(:,:) = 1._wp / rdt 240 !$OMP PARALLEL DO schedule(static) private(jj,ji) 241 DO jj = 1, jpj 242 DO ji = 1, jpi 243 frq_rst_e3t(ji,jj) = 0._wp !Ignore namelist settings 244 frq_rst_hdv(ji,jj) = 1._wp / rdt 245 END DO 246 END DO 210 247 ENDIF 211 248 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 249 !$OMP PARALLEL DO schedule(static) private(jj,ji) 212 250 DO jj = 1, jpj 213 251 DO ji = 1, jpi … … 305 343 ! ! --------------------------------------------- ! 306 344 ! 307 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 345 !$OMP PARALLEL 346 !$OMP DO schedule(static) private(jj,ji) 347 DO jj = 1, jpj 348 DO ji = 1, jpi 349 z_scale(ji,jj) = ( ssha(ji,jj) - sshb(ji,jj) ) * ssmask(ji,jj) / ( ht_0(ji,jj) + sshn(ji,jj) + 1. - ssmask(ji,jj) ) 350 END DO 351 END DO 352 !$OMP DO schedule(static) private(jk,jj,ji) 308 353 DO jk = 1, jpkm1 309 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 310 e3t_a(:,:,jk) = e3t_b(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 311 END DO 354 DO jj = 1, jpj 355 DO ji = 1, jpi 356 ! formally this is the same as e3t_a = e3t_0*(1+ssha/ht_0) 357 e3t_a(ji,jj,jk) = e3t_b(ji,jj,jk) + e3t_n(ji,jj,jk) * z_scale(ji,jj) * tmask(ji,jj,jk) 358 END DO 359 END DO 360 END DO 361 !$OMP END PARALLEL 312 362 ! 313 363 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! … … 318 368 ! 1 - barotropic divergence 319 369 ! ------------------------- 320 zhdiv(:,:) = 0._wp 321 zht(:,:) = 0._wp 370 !$OMP PARALLEL 371 !$OMP DO schedule(static) private(jj,ji) 372 DO jj = 1, jpj 373 DO ji = 1, jpi 374 zhdiv(ji,jj) = 0._wp 375 zht(ji,jj) = 0._wp 376 END DO 377 END DO 322 378 DO jk = 1, jpkm1 323 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 324 zht (:,:) = zht (:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 325 END DO 326 zhdiv(:,:) = zhdiv(:,:) / ( zht(:,:) + 1. - tmask_i(:,:) ) 379 !$OMP DO schedule(static) private(jj,ji) 380 DO jj = 1, jpj 381 DO ji = 1, jpi 382 zhdiv(ji,jj) = zhdiv(ji,jj) + e3t_n(ji,jj,jk) * hdivn(ji,jj,jk) 383 zht (ji,jj) = zht (ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 384 END DO 385 END DO 386 END DO 387 !$OMP DO schedule(static) private(jj,ji) 388 DO jj = 1, jpj 389 DO ji = 1, jpi 390 zhdiv(ji,jj) = zhdiv(ji,jj) / ( zht(ji,jj) + 1. - tmask_i(ji,jj) ) 391 END DO 392 END DO 393 !$OMP END PARALLEL 327 394 328 395 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 330 397 IF( ln_vvl_ztilde ) THEN 331 398 IF( kt > nit000 ) THEN 399 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 332 400 DO jk = 1, jpkm1 333 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 334 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 401 DO jj = 1, jpj 402 DO ji = 1, jpi 403 hdiv_lf(ji,jj,jk) = hdiv_lf(ji,jj,jk) - rdt * frq_rst_hdv(ji,jj) & 404 & * ( hdiv_lf(ji,jj,jk) - e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) ) 405 END DO 406 END DO 335 407 END DO 336 408 ENDIF … … 339 411 ! II - after z_tilde increments of vertical scale factors 340 412 ! ======================================================= 341 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 413 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 414 DO jk = 1, jpk 415 DO jj = 1, jpj 416 DO ji = 1, jpi 417 tilde_e3t_a(ji,jj,jk) = 0._wp ! tilde_e3t_a used to store tendency terms 418 END DO 419 END DO 420 END DO 342 421 343 422 ! 1 - High frequency divergence term 344 423 ! ---------------------------------- 345 424 IF( ln_vvl_ztilde ) THEN ! z_tilde case 425 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 346 426 DO jk = 1, jpkm1 347 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 427 DO jj = 1, jpj 428 DO ji = 1, jpi 429 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - ( e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) - hdiv_lf(ji,jj,jk) ) 430 END DO 431 END DO 348 432 END DO 349 433 ELSE ! layer case 434 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 350 435 DO jk = 1, jpkm1 351 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 436 DO jj = 1, jpj 437 DO ji = 1, jpi 438 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - e3t_n(ji,jj,jk) * ( hdivn(ji,jj,jk) - zhdiv(ji,jj) ) * tmask(ji,jj,jk) 439 END DO 440 END DO 352 441 END DO 353 442 ENDIF … … 356 445 ! ------------------ 357 446 IF( ln_vvl_ztilde ) THEN 447 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 358 448 DO jk = 1, jpk 359 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 449 DO jj = 1, jpj 450 DO ji = 1, jpi 451 tilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - frq_rst_e3t(ji,jj) * tilde_e3t_b(ji,jj,jk) 452 END DO 453 END DO 360 454 END DO 361 455 ENDIF … … 363 457 ! 3 - Thickness diffusion term 364 458 ! ---------------------------- 365 zwu(:,:) = 0._wp 366 zwv(:,:) = 0._wp 459 !$OMP PARALLEL 460 !$OMP DO schedule(static) private(jj,ji) 461 DO jj = 1, jpj 462 DO ji = 1, jpi 463 zwu(ji,jj) = 0._wp 464 zwv(ji,jj) = 0._wp 465 END DO 466 END DO 367 467 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 468 !$OMP DO schedule(static) private(jj,ji) 368 469 DO jj = 1, jpjm1 369 470 DO ji = 1, fs_jpim1 ! vector opt. … … 377 478 END DO 378 479 END DO 480 !$OMP DO schedule(static) private(jj,ji) 379 481 DO jj = 1, jpj ! b - correction for last oceanic u-v points 380 482 DO ji = 1, jpi … … 383 485 END DO 384 486 END DO 487 !$OMP DO schedule(static) private(jk,jj,ji) 385 488 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 386 489 DO jj = 2, jpjm1 … … 392 495 END DO 393 496 END DO 497 !$OMP END PARALLEL 394 498 ! ! d - thickness diffusion transport: boundary conditions 395 499 ! (stored for tracer advction and continuity equation) … … 407 511 ENDIF 408 512 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 409 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 513 !$OMP PARALLEL 514 !$OMP DO schedule(static) private(jk,jj,ji) 515 DO jk = 1, jpk 516 DO jj = 1, jpj 517 DO ji = 1, jpi 518 tilde_e3t_a(ji,jj,jk) = tilde_e3t_b(ji,jj,jk) + z2dt * tmask(ji,jj,jk) * tilde_e3t_a(ji,jj,jk) 519 END DO 520 END DO 521 END DO 410 522 411 523 ! Maximum deformation control 412 524 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 413 ze3t(:,:,jpk) = 0._wp 525 !$OMP DO schedule(static) private(jj,ji) 526 DO jj = 1, jpj 527 DO ji = 1, jpi 528 ze3t(ji,jj,jpk) = 0._wp 529 END DO 530 END DO 531 !$OMP DO schedule(static) private(jk,jj,ji) 414 532 DO jk = 1, jpkm1 415 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 416 END DO 533 DO jj = 1, jpj 534 DO ji = 1, jpi 535 ze3t(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) / e3t_0(ji,jj,jk) * tmask(ji,jj,jk) * tmask_i(ji,jj) 536 END DO 537 END DO 538 END DO 539 !$OMP END PARALLEL 417 540 z_tmax = MAXVAL( ze3t(:,:,:) ) 418 541 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain … … 442 565 ! - ML - end test 443 566 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 444 tilde_e3t_a(:,:,:) = MIN( tilde_e3t_a(:,:,:), rn_zdef_max * e3t_0(:,:,:) ) 445 tilde_e3t_a(:,:,:) = MAX( tilde_e3t_a(:,:,:), - rn_zdef_max * e3t_0(:,:,:) ) 567 !$OMP PARALLEL 568 !$OMP DO schedule(static) private(jk,jj,ji) 569 DO jk = 1, jpk 570 DO jj = 1, jpj 571 DO ji = 1, jpi 572 tilde_e3t_a(ji,jj,jk) = MIN( tilde_e3t_a(ji,jj,jk), rn_zdef_max * e3t_0(ji,jj,jk) ) 573 tilde_e3t_a(ji,jj,jk) = MAX( tilde_e3t_a(ji,jj,jk), - rn_zdef_max * e3t_0(ji,jj,jk) ) 574 END DO 575 END DO 576 END DO 446 577 447 578 ! 448 579 ! "tilda" change in the after scale factor 449 580 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 581 !$OMP DO schedule(static) private(jk,jj,ji) 450 582 DO jk = 1, jpkm1 451 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 583 DO jj = 1, jpj 584 DO ji = 1, jpi 585 dtilde_e3t_a(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) - tilde_e3t_b(ji,jj,jk) 586 END DO 587 END DO 452 588 END DO 453 589 ! III - Barotropic repartition of the sea surface height over the baroclinic profile … … 457 593 ! i.e. locally and not spread over the water column. 458 594 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 459 zht(:,:) = 0. 595 !$OMP DO schedule(static) private(jj,ji) 596 DO jj = 1, jpj 597 DO ji = 1, jpi 598 zht(ji,jj) = 0. 599 END DO 600 END DO 460 601 DO jk = 1, jpkm1 461 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 462 END DO 463 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 602 !$OMP DO schedule(static) private(jj,ji) 603 DO jj = 1, jpj 604 DO ji = 1, jpi 605 zht(ji,jj) = zht(ji,jj) + tilde_e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 606 END DO 607 END DO 608 END DO 609 !$OMP DO schedule(static) private(jj,ji) 610 DO jj = 1, jpj 611 DO ji = 1, jpi 612 z_scale(ji,jj) = - zht(ji,jj) / ( ht_0(ji,jj) + sshn(ji,jj) + 1. - ssmask(ji,jj) ) 613 END DO 614 END DO 615 !$OMP DO schedule(static) private(jk,jj,ji) 464 616 DO jk = 1, jpkm1 465 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 466 END DO 467 617 DO jj = 1, jpj 618 DO ji = 1, jpi 619 dtilde_e3t_a(ji,jj,jk) = dtilde_e3t_a(ji,jj,jk) + e3t_n(ji,jj,jk) * z_scale(ji,jj) * tmask(ji,jj,jk) 620 END DO 621 END DO 622 END DO 623 !$OMP END PARALLEL 468 624 ENDIF 469 625 470 626 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 471 627 ! ! ---baroclinic part--------- ! 628 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 472 629 DO jk = 1, jpkm1 473 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 630 DO jj = 1, jpj 631 DO ji = 1, jpi 632 e3t_a(ji,jj,jk) = e3t_a(ji,jj,jk) + dtilde_e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 633 END DO 634 END DO 474 635 END DO 475 636 ENDIF … … 484 645 END IF 485 646 ! 486 zht(:,:) = 0.0_wp 647 !$OMP PARALLEL 648 !$OMP DO schedule(static) private(jj,ji) 649 DO jj = 1, jpj 650 DO ji = 1, jpi 651 zht(ji,jj) = 0.0_wp 652 END DO 653 END DO 487 654 DO jk = 1, jpkm1 488 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 489 END DO 655 !$OMP DO schedule(static) private(jj,ji) 656 DO jj = 1, jpj 657 DO ji = 1, jpi 658 zht(ji,jj) = zht(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 659 END DO 660 END DO 661 END DO 662 !$OMP END PARALLEL 490 663 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 491 664 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 492 665 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 493 666 ! 494 zht(:,:) = 0.0_wp 667 !$OMP PARALLEL 668 !$OMP DO schedule(static) private(jj,ji) 669 DO jj = 1, jpj 670 DO ji = 1, jpi 671 zht(ji,jj) = 0.0_wp 672 END DO 673 END DO 495 674 DO jk = 1, jpkm1 496 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 497 END DO 675 !$OMP DO schedule(static) private(jj,ji) 676 DO jj = 1, jpj 677 DO ji = 1, jpi 678 zht(ji,jj) = zht(ji,jj) + e3t_a(ji,jj,jk) * tmask(ji,jj,jk) 679 END DO 680 END DO 681 END DO 682 !$OMP END PARALLEL 498 683 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 499 684 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 500 685 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 501 686 ! 502 zht(:,:) = 0.0_wp 687 !$OMP PARALLEL 688 !$OMP DO schedule(static) private(jj,ji) 689 DO jj = 1, jpj 690 DO ji = 1, jpi 691 zht(ji,jj) = 0.0_wp 692 END DO 693 END DO 503 694 DO jk = 1, jpkm1 504 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 505 END DO 695 !$OMP DO schedule(static) private(jj,ji) 696 DO jj = 1, jpj 697 DO ji = 1, jpi 698 zht(ji,jj) = zht(ji,jj) + e3t_b(ji,jj,jk) * tmask(ji,jj,jk) 699 END DO 700 END DO 701 END DO 702 !$OMP END PARALLEL 506 703 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 507 704 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain … … 532 729 ! *********************************** ! 533 730 534 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 535 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 731 !$OMP PARALLEL 732 !$OMP DO schedule(static) private(jj,ji) 733 DO jj = 1, jpj 734 DO ji = 1, jpi 735 hu_a(ji,jj) = e3u_a(ji,jj,1) * umask(ji,jj,1) 736 hv_a(ji,jj) = e3v_a(ji,jj,1) * vmask(ji,jj,1) 737 END DO 738 END DO 536 739 DO jk = 2, jpkm1 537 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 538 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 740 !$OMP DO schedule(static) private(jj,ji) 741 DO jj = 1, jpj 742 DO ji = 1, jpi 743 hu_a(ji,jj) = hu_a(ji,jj) + e3u_a(ji,jj,jk) * umask(ji,jj,jk) 744 hv_a(ji,jj) = hv_a(ji,jj) + e3v_a(ji,jj,jk) * vmask(ji,jj,jk) 745 END DO 746 END DO 539 747 END DO 540 748 ! ! Inverse of the local depth 541 749 !!gm BUG ? don't understand the use of umask_i here ..... 542 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 543 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 750 !$OMP DO schedule(static) private(jj,ji) 751 DO jj = 1, jpj 752 DO ji = 1, jpi 753 r1_hu_a(ji,jj) = ssumask(ji,jj) / ( hu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 754 r1_hv_a(ji,jj) = ssvmask(ji,jj) / ( hv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 755 END DO 756 END DO 757 !$OMP END PARALLEL 544 758 ! 545 759 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) … … 596 810 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 597 811 IF( neuler == 0 .AND. kt == nit000 ) THEN 598 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 812 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 813 DO jk = 1, jpk 814 DO jj = 1, jpj 815 DO ji = 1, jpi 816 tilde_e3t_b(ji,jj,jk) = tilde_e3t_n(ji,jj,jk) 817 END DO 818 END DO 819 END DO 599 820 ELSE 600 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 601 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 821 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 822 DO jk = 1, jpk 823 DO jj = 1, jpj 824 DO ji = 1, jpi 825 tilde_e3t_b(ji,jj,jk) = tilde_e3t_n(ji,jj,jk) & 826 & + atfp * ( tilde_e3t_b(ji,jj,jk) - 2.0_wp * tilde_e3t_n(ji,jj,jk) + tilde_e3t_a(ji,jj,jk) ) 827 END DO 828 END DO 829 END DO 602 830 ENDIF 603 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 831 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 832 DO jk = 1, jpk 833 DO jj = 1, jpj 834 DO ji = 1, jpi 835 tilde_e3t_n(ji,jj,jk) = tilde_e3t_a(ji,jj,jk) 836 END DO 837 END DO 838 END DO 604 839 ENDIF 605 gdept_b(:,:,:) = gdept_n(:,:,:) 606 gdepw_b(:,:,:) = gdepw_n(:,:,:) 607 608 e3t_n(:,:,:) = e3t_a(:,:,:) 609 e3u_n(:,:,:) = e3u_a(:,:,:) 610 e3v_n(:,:,:) = e3v_a(:,:,:) 840 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 841 DO jk = 1, jpk 842 DO jj = 1, jpj 843 DO ji = 1, jpi 844 gdept_b(ji,jj,jk) = gdept_n(ji,jj,jk) 845 gdepw_b(ji,jj,jk) = gdepw_n(ji,jj,jk) 846 847 e3t_n(ji,jj,jk) = e3t_a(ji,jj,jk) 848 e3u_n(ji,jj,jk) = e3u_a(ji,jj,jk) 849 e3v_n(ji,jj,jk) = e3v_a(ji,jj,jk) 850 END DO 851 END DO 852 END DO 611 853 612 854 ! Compute all missing vertical scale factor and depths … … 628 870 629 871 ! t- and w- points depth (set the isf depth as it is in the initial step) 630 gdept_n(:,:,1) = 0.5_wp * e3w_n(:,:,1) 631 gdepw_n(:,:,1) = 0.0_wp 632 gde3w_n(:,:,1) = gdept_n(:,:,1) - sshn(:,:) 872 ! !$OMP PARALLEL 873 ! !$OMP DO schedule(static) private(jj,ji) 874 DO jj = 1, jpj 875 DO ji = 1, jpi 876 gdept_n(ji,jj,1) = 0.5_wp * e3w_n(ji,jj,1) 877 gdepw_n(ji,jj,1) = 0.0_wp 878 gde3w_n(ji,jj,1) = gdept_n(ji,jj,1) - sshn(ji,jj) 879 END DO 880 END DO 633 881 DO jk = 2, jpk 882 ! !$OMP DO schedule(static) private(jj,ji,zcoef) 634 883 DO jj = 1,jpj 635 884 DO ji = 1,jpi … … 647 896 ! Local depth and Inverse of the local depth of the water 648 897 ! ------------------------------------------------------- 649 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 650 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 651 ! 652 ht_n(:,:) = e3t_n(:,:,1) * tmask(:,:,1) 898 !$OMP PARALLEL 899 !$OMP DO schedule(static) private(jj,ji) 900 DO jj = 1, jpj 901 DO ji = 1, jpi 902 hu_n(ji,jj) = hu_a(ji,jj) ; r1_hu_n(ji,jj) = r1_hu_a(ji,jj) 903 hv_n(ji,jj) = hv_a(ji,jj) ; r1_hv_n(ji,jj) = r1_hv_a(ji,jj) 904 ! 905 ht_n(ji,jj) = e3t_n(ji,jj,1) * tmask(ji,jj,1) 906 END DO 907 END DO 653 908 DO jk = 2, jpkm1 654 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 655 END DO 656 909 !$OMP DO schedule(static) private(jj,ji) 910 DO jj = 1, jpj 911 DO ji = 1, jpi 912 ht_n(ji,jj) = ht_n(ji,jj) + e3t_n(ji,jj,jk) * tmask(ji,jj,jk) 913 END DO 914 END DO 915 END DO 916 !$OMP END PARALLEL 657 917 ! write restart file 658 918 ! ================== … … 694 954 ! 695 955 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 956 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 696 957 DO jk = 1, jpk 697 958 DO jj = 1, jpjm1 … … 704 965 END DO 705 966 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 706 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 967 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 968 DO jk = 1, jpk 969 DO jj = 1, jpj 970 DO ji = 1, jpi 971 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3u_0(ji,jj,jk) 972 END DO 973 END DO 974 END DO 707 975 ! 708 976 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 977 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 709 978 DO jk = 1, jpk 710 979 DO jj = 1, jpjm1 … … 717 986 END DO 718 987 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 988 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 989 DO jk = 1, jpk 990 DO jj = 1, jpj 991 DO ji = 1, jpi 992 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3v_0(ji,jj,jk) 993 END DO 994 END DO 995 END DO 720 996 ! 721 997 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 998 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 722 999 DO jk = 1, jpk 723 1000 DO jj = 1, jpjm1 … … 731 1008 END DO 732 1009 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 733 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 1010 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji) 1011 DO jk = 1, jpk 1012 DO jj = 1, jpj 1013 DO ji = 1, jpi 1014 pe3_out(ji,jj,jk) = pe3_out(ji,jj,jk) + e3f_0(ji,jj,jk) 1015 END DO 1016 END DO 1017 END DO 734 1018 ! 735 1019 CASE( 'W' ) !* from T- to W-point : vertical simple mean 736 1020 ! 737 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 1021 !$OMP PARALLEL 1022 !$OMP DO schedule(static) private(jj,ji) 1023 DO jj = 1, jpj 1024 DO ji = 1, jpi 1025 pe3_out(ji,jj,1) = e3w_0(ji,jj,1) + pe3_in(ji,jj,1) - e3t_0(ji,jj,1) 1026 END DO 1027 END DO 738 1028 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 739 1029 !!gm BUG? use here wmask in case of ISF ? to be checked 1030 !$OMP DO schedule(static) private(jk,jj,ji) 740 1031 DO jk = 2, jpk 741 pe3_out(:,:,jk) = e3w_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 742 & * ( pe3_in(:,:,jk-1) - e3t_0(:,:,jk-1) ) & 743 & + 0.5_wp * ( tmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 744 & * ( pe3_in(:,:,jk ) - e3t_0(:,:,jk ) ) 745 END DO 1032 DO jj = 1, jpj 1033 DO ji = 1, jpi 1034 pe3_out(ji,jj,jk) = e3w_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( tmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 1035 & * ( pe3_in(ji,jj,jk-1) - e3t_0(ji,jj,jk-1) ) & 1036 & + 0.5_wp * ( tmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & 1037 & * ( pe3_in(ji,jj,jk ) - e3t_0(ji,jj,jk ) ) 1038 END DO 1039 END DO 1040 END DO 1041 !$OMP END PARALLEL 746 1042 ! 747 1043 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean 748 1044 ! 749 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 1045 !$OMP PARALLEL 1046 !$OMP DO schedule(static) private(jj,ji) 1047 DO jj = 1, jpj 1048 DO ji = 1, jpi 1049 pe3_out(ji,jj,1) = e3uw_0(ji,jj,1) + pe3_in(ji,jj,1) - e3u_0(ji,jj,1) 1050 END DO 1051 END DO 750 1052 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 751 1053 !!gm BUG? use here wumask in case of ISF ? to be checked 1054 !$OMP DO schedule(static) private(jk,jj,ji) 752 1055 DO jk = 2, jpk 753 pe3_out(:,:,jk) = e3uw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 754 & * ( pe3_in(:,:,jk-1) - e3u_0(:,:,jk-1) ) & 755 & + 0.5_wp * ( umask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 756 & * ( pe3_in(:,:,jk ) - e3u_0(:,:,jk ) ) 757 END DO 1056 DO jj = 1, jpj 1057 DO ji = 1, jpi 1058 pe3_out(ji,jj,jk) = e3uw_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 1059 & * ( pe3_in(ji,jj,jk-1) - e3u_0(ji,jj,jk-1) ) & 1060 & + 0.5_wp * ( umask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & 1061 & * ( pe3_in(ji,jj,jk ) - e3u_0(ji,jj,jk ) ) 1062 END DO 1063 END DO 1064 END DO 1065 !$OMP END PARALLEL 758 1066 ! 759 1067 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean 760 1068 ! 761 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 1069 !$OMP PARALLEL 1070 !$OMP DO schedule(static) private(jj,ji) 1071 DO jj = 1, jpj 1072 DO ji = 1, jpi 1073 pe3_out(ji,jj,1) = e3vw_0(ji,jj,1) + pe3_in(ji,jj,1) - e3v_0(ji,jj,1) 1074 END DO 1075 END DO 762 1076 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 763 1077 !!gm BUG? use here wvmask in case of ISF ? to be checked 1078 !$OMP DO schedule(static) private(jk,jj,ji) 764 1079 DO jk = 2, jpk 765 pe3_out(:,:,jk) = e3vw_0(:,:,jk) + ( 1.0_wp - 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 766 & * ( pe3_in(:,:,jk-1) - e3v_0(:,:,jk-1) ) & 767 & + 0.5_wp * ( vmask(:,:,jk) * (1.0_wp - zlnwd) + zlnwd ) & 768 & * ( pe3_in(:,:,jk ) - e3v_0(:,:,jk ) ) 769 END DO 1080 DO jj = 1, jpj 1081 DO ji = 1, jpi 1082 pe3_out(ji,jj,jk) = e3vw_0(ji,jj,jk) + ( 1.0_wp - 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) ) & 1083 & * ( pe3_in(ji,jj,jk-1) - e3v_0(ji,jj,jk-1) ) & 1084 & + 0.5_wp * ( vmask(ji,jj,jk) * (1.0_wp - zlnwd) + zlnwd ) & 1085 & * ( pe3_in(ji,jj,jk ) - e3v_0(ji,jj,jk ) ) 1086 END DO 1087 END DO 1088 END DO 1089 !$OMP END PARALLEL 770 1090 END SELECT 771 1091 ! … … 905 1225 sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) !!gm I don't understand that ! 906 1226 sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 907 ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 1227 ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 908 1228 ENDIF 909 1229 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
