Changeset 7753 for trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
- Timestamp:
- 2017-03-03T12:46:59+01:00 (7 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DOM/domvvl.F90
r7698 r7753 135 135 ! ! Read or initialize e3t_(b/n), tilde_e3t_(b/n) and hdiv_lf 136 136 CALL dom_vvl_rst( nit000, 'READ' ) 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 137 e3t_a(:,:,jpk) = e3t_0(:,:,jpk) ! last level always inside the sea floor set one for all 143 138 ! 144 139 ! !== Set of all other vertical scale factors ==! (now and before) … … 158 153 ! 159 154 ! !== depth of t and w-point ==! (set the isf depth as it is in the initial timestep) 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 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 171 160 DO jk = 2, jpk ! vertical sum 172 !$OMP DO schedule(static) private(jj,ji,zcoef)173 161 DO jj = 1,jpj 174 162 DO ji = 1,jpi … … 190 178 ! 191 179 ! !== thickness of the water column !! (ocean portion only) 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 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) 185 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) 201 191 END DO 202 DO jk = 2, jpkm1203 !$OMP DO schedule(static) private(jj,ji)204 DO jj = 1, jpj205 DO ji = 1, jpi206 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 DO212 END DO213 END DO214 192 ! 215 193 ! !== inverse of water column thickness ==! (u- and v- points) 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 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 226 199 ! !== z_tilde coordinate case ==! (Restoring frequencies) 227 200 IF( ln_vvl_ztilde ) THEN … … 229 202 ! ! Values in days provided via the namelist 230 203 ! ! use rsmall to avoid possible division by zero errors with faulty settings 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 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 ) 238 206 ! 239 207 IF( ln_vvl_ztilde_as_zstar ) THEN ! z-star emulation using z-tile 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 208 frq_rst_e3t(:,:) = 0._wp !Ignore namelist settings 209 frq_rst_hdv(:,:) = 1._wp / rdt 247 210 ENDIF 248 211 IF ( ln_vvl_zstar_at_eqtor ) THEN ! use z-star in vicinity of the Equator 249 !$OMP PARALLEL DO schedule(static) private(jj,ji)250 212 DO jj = 1, jpj 251 213 DO ji = 1, jpi … … 343 305 ! ! --------------------------------------------- ! 344 306 ! 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 307 z_scale(:,:) = ( ssha(:,:) - sshb(:,:) ) * ssmask(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 308 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) 351 311 END DO 352 !$OMP DO schedule(static) private(jk,jj,ji)353 DO jk = 1, jpkm1354 DO jj = 1, jpj355 DO ji = 1, jpi356 ! 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 DO359 END DO360 END DO361 !$OMP END PARALLEL362 312 ! 363 313 IF( ln_vvl_ztilde .OR. ln_vvl_layer .AND. ll_do_bclinic ) THEN ! z_tilde or layer coordinate ! … … 368 318 ! 1 - barotropic divergence 369 319 ! ------------------------- 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 320 zhdiv(:,:) = 0._wp 321 zht(:,:) = 0._wp 378 322 DO jk = 1, jpkm1 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 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(:,:) ) 394 327 395 328 ! 2 - Low frequency baroclinic horizontal divergence (z-tilde case only) … … 397 330 IF( ln_vvl_ztilde ) THEN 398 331 IF( kt > nit000 ) THEN 399 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)400 332 DO jk = 1, jpkm1 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 333 hdiv_lf(:,:,jk) = hdiv_lf(:,:,jk) - rdt * frq_rst_hdv(:,:) & 334 & * ( hdiv_lf(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) ) 407 335 END DO 408 336 ENDIF … … 411 339 ! II - after z_tilde increments of vertical scale factors 412 340 ! ======================================================= 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 341 tilde_e3t_a(:,:,:) = 0._wp ! tilde_e3t_a used to store tendency terms 421 342 422 343 ! 1 - High frequency divergence term 423 344 ! ---------------------------------- 424 345 IF( ln_vvl_ztilde ) THEN ! z_tilde case 425 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)426 346 DO jk = 1, jpkm1 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 347 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - ( e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) - hdiv_lf(:,:,jk) ) 432 348 END DO 433 349 ELSE ! layer case 434 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)435 350 DO jk = 1, jpkm1 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 351 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - e3t_n(:,:,jk) * ( hdivn(:,:,jk) - zhdiv(:,:) ) * tmask(:,:,jk) 441 352 END DO 442 353 ENDIF … … 445 356 ! ------------------ 446 357 IF( ln_vvl_ztilde ) THEN 447 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)448 358 DO jk = 1, jpk 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 359 tilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - frq_rst_e3t(:,:) * tilde_e3t_b(:,:,jk) 454 360 END DO 455 361 ENDIF … … 457 363 ! 3 - Thickness diffusion term 458 364 ! ---------------------------- 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 365 zwu(:,:) = 0._wp 366 zwv(:,:) = 0._wp 467 367 DO jk = 1, jpkm1 ! a - first derivative: diffusive fluxes 468 !$OMP DO schedule(static) private(jj,ji)469 368 DO jj = 1, jpjm1 470 369 DO ji = 1, fs_jpim1 ! vector opt. … … 478 377 END DO 479 378 END DO 480 !$OMP DO schedule(static) private(jj,ji)481 379 DO jj = 1, jpj ! b - correction for last oceanic u-v points 482 380 DO ji = 1, jpi … … 485 383 END DO 486 384 END DO 487 !$OMP DO schedule(static) private(jk,jj,ji)488 385 DO jk = 1, jpkm1 ! c - second derivative: divergence of diffusive fluxes 489 386 DO jj = 2, jpjm1 … … 495 392 END DO 496 393 END DO 497 !$OMP END PARALLEL498 394 ! ! d - thickness diffusion transport: boundary conditions 499 395 ! (stored for tracer advction and continuity equation) … … 511 407 ENDIF 512 408 CALL lbc_lnk( tilde_e3t_a(:,:,:), 'T', 1._wp ) 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 409 tilde_e3t_a(:,:,:) = tilde_e3t_b(:,:,:) + z2dt * tmask(:,:,:) * tilde_e3t_a(:,:,:) 522 410 523 411 ! Maximum deformation control 524 412 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~ 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) 413 ze3t(:,:,jpk) = 0._wp 532 414 DO jk = 1, jpkm1 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 415 ze3t(:,:,jk) = tilde_e3t_a(:,:,jk) / e3t_0(:,:,jk) * tmask(:,:,jk) * tmask_i(:,:) 416 END DO 540 417 z_tmax = MAXVAL( ze3t(:,:,:) ) 541 418 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain … … 565 442 ! - ML - end test 566 443 ! - ML - Imposing these limits will cause a baroclinicity error which is corrected for below 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 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(:,:,:) ) 577 446 578 447 ! 579 448 ! "tilda" change in the after scale factor 580 449 ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 581 !$OMP DO schedule(static) private(jk,jj,ji)582 450 DO jk = 1, jpkm1 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 451 dtilde_e3t_a(:,:,jk) = tilde_e3t_a(:,:,jk) - tilde_e3t_b(:,:,jk) 588 452 END DO 589 453 ! III - Barotropic repartition of the sea surface height over the baroclinic profile … … 593 457 ! i.e. locally and not spread over the water column. 594 458 ! (keep in mind that the idea is to reduce Eulerian velocity as much as possible) 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 459 zht(:,:) = 0. 601 460 DO jk = 1, jpkm1 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) 461 zht(:,:) = zht(:,:) + tilde_e3t_a(:,:,jk) * tmask(:,:,jk) 462 END DO 463 z_scale(:,:) = - zht(:,:) / ( ht_0(:,:) + sshn(:,:) + 1. - ssmask(:,:) ) 616 464 DO jk = 1, jpkm1 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 465 dtilde_e3t_a(:,:,jk) = dtilde_e3t_a(:,:,jk) + e3t_n(:,:,jk) * z_scale(:,:) * tmask(:,:,jk) 466 END DO 467 624 468 ENDIF 625 469 626 470 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN ! z_tilde or layer coordinate ! 627 471 ! ! ---baroclinic part--------- ! 628 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)629 472 DO jk = 1, jpkm1 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 473 e3t_a(:,:,jk) = e3t_a(:,:,jk) + dtilde_e3t_a(:,:,jk) * tmask(:,:,jk) 635 474 END DO 636 475 ENDIF … … 645 484 END IF 646 485 ! 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 486 zht(:,:) = 0.0_wp 654 487 DO jk = 1, jpkm1 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 488 zht(:,:) = zht(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 489 END DO 663 490 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshn(:,:) - zht(:,:) ) ) 664 491 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 665 492 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+sshn-SUM(e3t_n))) =', z_tmax 666 493 ! 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 494 zht(:,:) = 0.0_wp 674 495 DO jk = 1, jpkm1 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 496 zht(:,:) = zht(:,:) + e3t_a(:,:,jk) * tmask(:,:,jk) 497 END DO 683 498 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + ssha(:,:) - zht(:,:) ) ) 684 499 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain 685 500 IF( lwp ) WRITE(numout, *) kt,' MAXVAL(abs(ht_0+ssha-SUM(e3t_a))) =', z_tmax 686 501 ! 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 502 zht(:,:) = 0.0_wp 694 503 DO jk = 1, jpkm1 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 504 zht(:,:) = zht(:,:) + e3t_b(:,:,jk) * tmask(:,:,jk) 505 END DO 703 506 z_tmax = MAXVAL( tmask(:,:,1) * tmask_i(:,:) * ABS( ht_0(:,:) + sshb(:,:) - zht(:,:) ) ) 704 507 IF( lk_mpp ) CALL mpp_max( z_tmax ) ! max over the global domain … … 729 532 ! *********************************** ! 730 533 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 534 hu_a(:,:) = e3u_a(:,:,1) * umask(:,:,1) 535 hv_a(:,:) = e3v_a(:,:,1) * vmask(:,:,1) 739 536 DO jk = 2, jpkm1 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 537 hu_a(:,:) = hu_a(:,:) + e3u_a(:,:,jk) * umask(:,:,jk) 538 hv_a(:,:) = hv_a(:,:) + e3v_a(:,:,jk) * vmask(:,:,jk) 747 539 END DO 748 540 ! ! Inverse of the local depth 749 541 !!gm BUG ? don't understand the use of umask_i here ..... 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 542 r1_hu_a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) ) 543 r1_hv_a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) ) 758 544 ! 759 545 CALL wrk_dealloc( jpi,jpj, zht, z_scale, zwu, zwv, zhdiv ) … … 810 596 IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN 811 597 IF( neuler == 0 .AND. kt == nit000 ) THEN 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 598 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) 820 599 ELSE 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 600 tilde_e3t_b(:,:,:) = tilde_e3t_n(:,:,:) & 601 & + atfp * ( tilde_e3t_b(:,:,:) - 2.0_wp * tilde_e3t_n(:,:,:) + tilde_e3t_a(:,:,:) ) 830 602 ENDIF 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 839 ENDIF 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 603 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 604 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(:,:,:) 853 611 854 612 ! Compute all missing vertical scale factor and depths … … 870 628 871 629 ! t- and w- points depth (set the isf depth as it is in the initial step) 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 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(:,:) 881 633 DO jk = 2, jpk 882 ! !$OMP DO schedule(static) private(jj,ji,zcoef)883 634 DO jj = 1,jpj 884 635 DO ji = 1,jpi … … 896 647 ! Local depth and Inverse of the local depth of the water 897 648 ! ------------------------------------------------------- 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 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) 653 DO jk = 2, jpkm1 654 ht_n(:,:) = ht_n(:,:) + e3t_n(:,:,jk) * tmask(:,:,jk) 907 655 END DO 908 DO jk = 2, jpkm1 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 656 917 657 ! write restart file 918 658 ! ================== … … 954 694 ! 955 695 CASE( 'U' ) !* from T- to U-point : hor. surface weighted mean 956 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)957 696 DO jk = 1, jpk 958 697 DO jj = 1, jpjm1 … … 965 704 END DO 966 705 CALL lbc_lnk( pe3_out(:,:,:), 'U', 1._wp ) 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 706 pe3_out(:,:,:) = pe3_out(:,:,:) + e3u_0(:,:,:) 975 707 ! 976 708 CASE( 'V' ) !* from T- to V-point : hor. surface weighted mean 977 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)978 709 DO jk = 1, jpk 979 710 DO jj = 1, jpjm1 … … 986 717 END DO 987 718 CALL lbc_lnk( pe3_out(:,:,:), 'V', 1._wp ) 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 719 pe3_out(:,:,:) = pe3_out(:,:,:) + e3v_0(:,:,:) 996 720 ! 997 721 CASE( 'F' ) !* from U-point to F-point : hor. surface weighted mean 998 !$OMP PARALLEL DO schedule(static) private(jk,jj,ji)999 722 DO jk = 1, jpk 1000 723 DO jj = 1, jpjm1 … … 1008 731 END DO 1009 732 CALL lbc_lnk( pe3_out(:,:,:), 'F', 1._wp ) 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 733 pe3_out(:,:,:) = pe3_out(:,:,:) + e3f_0(:,:,:) 1018 734 ! 1019 735 CASE( 'W' ) !* from T- to W-point : vertical simple mean 1020 736 ! 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 737 pe3_out(:,:,1) = e3w_0(:,:,1) + pe3_in(:,:,1) - e3t_0(:,:,1) 1028 738 ! - ML - The use of mask in this formulea enables the special treatment of the last w-point without indirect adressing 1029 739 !!gm BUG? use here wmask in case of ISF ? to be checked 1030 !$OMP DO schedule(static) private(jk,jj,ji)1031 740 DO jk = 2, jpk 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 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 1042 746 ! 1043 747 CASE( 'UW' ) !* from U- to UW-point : vertical simple mean 1044 748 ! 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 749 pe3_out(:,:,1) = e3uw_0(:,:,1) + pe3_in(:,:,1) - e3u_0(:,:,1) 1052 750 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 1053 751 !!gm BUG? use here wumask in case of ISF ? to be checked 1054 !$OMP DO schedule(static) private(jk,jj,ji)1055 752 DO jk = 2, jpk 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 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 1066 758 ! 1067 759 CASE( 'VW' ) !* from V- to VW-point : vertical simple mean 1068 760 ! 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 761 pe3_out(:,:,1) = e3vw_0(:,:,1) + pe3_in(:,:,1) - e3v_0(:,:,1) 1076 762 ! - ML - The use of mask in this formaula enables the special treatment of the last w- point without indirect adressing 1077 763 !!gm BUG? use here wvmask in case of ISF ? to be checked 1078 !$OMP DO schedule(static) private(jk,jj,ji)1079 764 DO jk = 2, jpk 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 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 1090 770 END SELECT 1091 771 ! … … 1225 905 sshb(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) !!gm I don't understand that ! 1226 906 sshn(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 1227 ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 907 ssha(ji,jj) = rn_wdmin1 - ht_wd(ji,jj) 1228 908 ENDIF 1229 909 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.