Changeset 5864
- Timestamp:
- 2015-11-05T17:38:47+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r5845 r5864 234 234 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 235 235 ! hydrostatic pressure gradient 236 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 237 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 238 !!gm zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 239 !!gm zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 236 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 237 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 240 238 ! add to the general momentum trend 241 239 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 253 251 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 254 252 & + zcoef1 * ( ( rhd(ji+1,jj,jk)+rhd(ji+1,jj,jk-1) ) & 255 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj) 256 !!gm & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 253 & - ( rhd(ji ,jj,jk)+rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 257 254 258 255 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 259 256 & + zcoef1 * ( ( rhd(ji,jj+1,jk)+rhd(ji,jj+1,jk-1) ) & 260 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj) 261 !!gm & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 257 & - ( rhd(ji,jj, jk)+rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 262 258 ! add to the general momentum trend 263 259 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 305 301 zcoef1 = zcoef0 * e3w_n(ji,jj,1) 306 302 ! hydrostatic pressure gradient 307 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) / e1u(ji,jj) 308 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) / e2v(ji,jj) 309 !!gm zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 310 !!gm zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 303 zhpi(ji,jj,1) = zcoef1 * ( rhd(ji+1,jj ,1) - rhd(ji,jj,1) ) * r1_e1u(ji,jj) 304 zhpj(ji,jj,1) = zcoef1 * ( rhd(ji ,jj+1,1) - rhd(ji,jj,1) ) * r1_e2v(ji,jj) 311 305 ! add to the general momentum trend 312 306 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 324 318 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 325 319 & + zcoef1 * ( ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) ) & 326 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) / e1u(ji,jj) 327 !!gm & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 320 & - ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 328 321 329 322 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 330 323 & + zcoef1 * ( ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) ) & 331 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) / e2v(ji,jj) 332 !!gm & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 324 & - ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 333 325 ! add to the general momentum trend 334 326 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 349 341 ua (ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) ! subtract old value 350 342 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) & ! compute the new one 351 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) / e1u(ji,jj) 352 !!gm & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 343 & + zcoef2 * ( rhd(ji+1,jj,iku-1) - rhd(ji,jj,iku-1) + gru(ji,jj) ) * r1_e1u(ji,jj) 353 344 ua (ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) ! add the new one to the general momentum trend 354 345 ENDIF … … 356 347 va (ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) ! subtract old value 357 348 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) & ! compute the new one 358 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) / e2v(ji,jj) 359 !!gm & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 349 & + zcoef3 * ( rhd(ji,jj+1,ikv-1) - rhd(ji,jj,ikv-1) + grv(ji,jj) ) * r1_e2v(ji,jj) 360 350 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 361 351 ENDIF … … 411 401 DO ji = fs_2, fs_jpim1 ! vector opt. 412 402 ! hydrostatic pressure gradient along s-surfaces 413 !!gm zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 414 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 415 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 416 !!gm zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 417 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 418 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 403 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( e3w_n(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 404 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 405 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( e3w_n(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 406 & - e3w_n(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 419 407 ! s-coordinate pressure gradient correction 420 408 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 421 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj) 422 !!gm & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 409 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 423 410 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 424 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj) 425 !!gm & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 411 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 426 412 ! add to the general momentum trend 427 413 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap … … 435 421 DO ji = fs_2, fs_jpim1 ! vector opt. 436 422 ! hydrostatic pressure gradient along s-surfaces 437 !!gm zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 438 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 423 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 * r1_e1u(ji,jj) & 439 424 & * ( e3w_n(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 440 425 & - e3w_n(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 441 !!gm zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) &442 426 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 * r1_e2v(ji,jj) & 443 427 & * ( e3w_n(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & … … 445 429 ! s-coordinate pressure gradient correction 446 430 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 447 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) 448 !!gm & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 431 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) 449 432 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 450 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) 451 !!gm & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 433 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) 452 434 ! add to the general momentum trend 453 435 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap … … 561 543 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 562 544 ! we assume ISF is in isostatic equilibrium 563 !!gm zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj ,iktp1i) & 564 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj ,iktp1i) & 565 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 566 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 567 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 568 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)) ) 569 !!gm zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji ,jj+1,iktp1j) & 570 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * e3w_n(ji ,jj+1,iktp1j) & 571 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 572 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 573 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 574 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 545 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( 0.5_wp * e3w_n(ji+1,jj ,iktp1i) & 546 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 547 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 548 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 549 & + ( ziceload(ji+1,jj) - ziceload(ji,jj) ) ) 550 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( 0.5_wp * e3w_n(ji ,jj+1,iktp1j) & 551 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 552 & - 0.5_wp * e3w_n(ji ,jj ,ikt ) & 553 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 554 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 575 555 ! s-coordinate pressure gradient correction (=0 if z coordinate) 576 556 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 577 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) / e1u(ji,jj) 578 !!gm & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 557 & * ( gde3w_n(ji+1,jj,1) - gde3w_n(ji,jj,1) ) * r1_e1u(ji,jj) 579 558 zvap = -zcoef0 * ( rhd (ji,jj+1,1) + rhd (ji,jj,1) + 2._wp * znad ) & 580 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) / e2v(ji,jj) 581 !!gm & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 559 & * ( gde3w_n(ji,jj+1,1) - gde3w_n(ji,jj,1) ) * r1_e2v(ji,jj) 582 560 ! add to the general momentum trend 583 561 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) … … 596 574 IF ( iku .GT. 1 ) THEN 597 575 ! case iku 598 !!gm zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu & 599 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu & 576 zhpi(ji,jj,iku) = zcoef0 * r1_e1u(ji,jj) * ze3wu & 600 577 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) & 601 578 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 602 579 ! corrective term ( = 0 if z coordinate ) 603 !!gm zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 604 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 580 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) * r1_e1u(ji,jj) 605 581 ! zhpi will be added in interior loop 606 582 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap … … 609 585 610 586 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 611 !!gm zhpiint = zcoef0 * r1_e1u(ji,jj) & 612 zhpiint = zcoef0 / e1u(ji,jj) & 587 zhpiint = zcoef0 * r1_e1u(ji,jj) & 613 588 & * ( e3w_n(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 614 589 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 615 590 & - e3w_n(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 616 591 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 617 !!gm zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 618 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 592 zhpi(ji,jj,iku+1) = zcoef0 * r1_e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 619 593 END IF 620 594 … … 624 598 IF ( ikv .GT. 1 ) THEN 625 599 ! case ikv 626 !!gm zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv & 627 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv & 628 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) & 629 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 600 zhpj(ji,jj,ikv) = zcoef0 * r1_e2v(ji,jj) * ze3wv & 601 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) & 602 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 630 603 ! corrective term ( = 0 if z coordinate ) 631 !!gm zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 632 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 604 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) * r1_e2v(ji,jj) 633 605 ! zhpi will be added in interior loop 634 606 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap … … 637 609 638 610 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 639 !!gm zhpjint = zcoef0 * r1_e2v(ji,jj) & 640 zhpjint = zcoef0 / e2v(ji,jj) & 611 zhpjint = zcoef0 * r1_e2v(ji,jj) & 641 612 & * ( e3w_n(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 642 613 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 643 614 & - e3w_n(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 644 615 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 645 !!gm zhpj(ji,jj,ikv+1) = zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 646 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 616 zhpj(ji,jj,ikv+1) = zcoef0 * r1_e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 647 617 END IF 648 618 END DO … … 655 625 DO jj = 2, jpjm1 656 626 DO ji = fs_2, fs_jpim1 ! vector opt. 657 !!gm useles !658 iku=miku(ji,jj); ikv=mikv(ji,jj)659 !!gm660 627 DO jk = 2, jpkm1 661 628 ! hydrostatic pressure gradient along s-surfaces 662 629 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 663 630 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 664 !!gm & + zcoef0 * r1_e1u(ji,jj) & 665 & + zcoef0 / e1u(ji,jj) & 631 & + zcoef0 * r1_e1u(ji,jj) & 666 632 & * ( e3w_n(ji+1,jj,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 667 633 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & … … 671 637 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 672 638 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 673 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 674 !!gm & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 639 & * ( gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk) ) * r1_e1u(ji,jj) * umask(ji,jj,jk-1) 675 640 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 676 641 … … 678 643 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 679 644 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 680 !!gm & + zcoef0 * r1_e2v(ji,jj) & 681 & + zcoef0 / e2v(ji,jj) & 645 & + zcoef0 * r1_e2v(ji,jj) & 682 646 & * ( e3w_n(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 683 647 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & … … 687 651 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 688 652 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 689 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 690 !!gm & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 653 & * ( gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk-1) 691 654 ! add to the general momentum trend 692 655 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) … … 707 670 ! remove old value (interior case) 708 671 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 709 & * ( gde3w_n(ji+1,jj ,iku) - gde3w_n(ji,jj,iku) ) / e1u(ji,jj) 710 !!gm & * ( gde3w_n(ji+1,jj ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 672 & * ( gde3w_n(ji+1,jj ,iku) - gde3w_n(ji,jj,iku) ) * r1_e1u(ji,jj) 711 673 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 712 674 ! put new value 713 675 ! -zpshpi to avoid double contribution of the partial step in the top layer 714 !!gm zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) * r1_e1u(ji,jj) 715 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 716 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 717 !!gm zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 676 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) * r1_e1u(ji,jj) 677 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 * r1_e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 718 678 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 719 679 END IF … … 722 682 ! remove old value (interior case) 723 683 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 724 & * ( gde3w_n(ji ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) / e2v(ji,jj) 725 !!gm & * ( gde3w_n(ji ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj) 684 & * ( gde3w_n(ji ,jj+1,ikv) - gde3w_n(ji,jj,ikv) ) * r1_e2v(ji,jj) 726 685 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 727 686 ! put new value 728 687 ! -zpshpj to avoid double contribution of the partial step in the top layer 729 !!gm zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj) 730 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 731 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 732 !!gm zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 688 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) * r1_e2v(ji,jj) 689 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 * r1_e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 733 690 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 734 691 END IF … … 938 895 DO jj = 2, jpjm1 939 896 DO ji = fs_2, fs_jpim1 ! vector opt. 940 !!gm zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 941 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) / e1u(ji,jj) 942 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) / e2v(ji,jj) 943 !!gm zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 897 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) - rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 898 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) - rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 944 899 ! add to the general momentum trend 945 900 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) … … 957 912 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 958 913 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) & 959 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) / e1u(ji,jj) 960 !!gm & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 914 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj) 961 915 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 962 916 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) & 963 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) / e2v(ji,jj) 964 !!gm & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 917 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj) 965 918 ! add to the general momentum trend 966 919 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) … … 992 945 !! 993 946 INTEGER :: ji, jj, jk, jkk ! dummy loop indices 994 REAL(wp) :: zcoef0, znad ! temporaryscalars995 ! !947 REAL(wp) :: zcoef0, znad ! local scalars 948 ! 996 949 !! The local variables for the correction term 997 950 INTEGER :: jk1, jis, jid, jjs, jjd … … 1004 957 !!---------------------------------------------------------------------- 1005 958 ! 1006 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp )1007 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh )1008 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n )959 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 960 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 961 CALL wrk_alloc( jpi,jpj, zsshu_n, zsshv_n ) 1009 962 ! 1010 963 IF( kt == nit000 ) THEN … … 1014 967 ENDIF 1015 968 1016 !!----------------------------------------------------------------------1017 969 ! Local constant initialization 1018 970 zcoef0 = - grav … … 1060 1012 ! constrained cubic spline interpolation 1061 1013 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 1062 CALL cspline( fsp,xsp,asp,bsp,csp,dsp,polynomial_type)1014 CALL cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1063 1015 1064 1016 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1065 1017 DO jj = 2, jpj 1066 1018 DO ji = 2, jpi 1067 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 1068 bsp(ji,jj,1), csp(ji,jj,1), & 1069 dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1019 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1020 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w_n(ji,jj,1) 1070 1021 1071 1022 ! assuming linear profile across the top half surface layer … … 1078 1029 DO jj = 2, jpj 1079 1030 DO ji = 2, jpi 1080 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1081 integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),&1082 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), &1083 csp(ji,jj,jk-1), dsp(ji,jj,jk-1))1031 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1032 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1033 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1034 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1084 1035 END DO 1085 1036 END DO … … 1091 1042 DO jj = 2, jpjm1 1092 1043 DO ji = 2, jpim1 1044 !!gm BUG ? if it is ssh at u- & v-point then it should be: 1045 ! zsshu_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji+1,jj) * sshn(ji+1,jj)) * & 1046 ! & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1047 ! zsshv_n(ji,jj) = (e1e2t(ji,jj) * sshn(ji,jj) + e1e2t(ji,jj+1) * sshn(ji,jj+1)) * & 1048 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1049 !!gm not this: 1093 1050 zsshu_n(ji,jj) = (e1e2u(ji,jj) * sshn(ji,jj) + e1e2u(ji+1, jj) * sshn(ji+1,jj)) * & 1094 1051 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp … … 1108 1065 DO jj = 2, jpjm1 1109 1066 DO ji = 2, jpim1 1110 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk)1111 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk)1067 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u_n(ji,jj,jk) 1068 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v_n(ji,jj,jk) 1112 1069 END DO 1113 1070 END DO … … 1126 1083 DO jj = 2, jpjm1 1127 1084 DO ji = 2, jpim1 1128 zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1129 zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk)))1130 zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1131 zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk)))1085 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1086 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1087 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1088 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1132 1089 END DO 1133 1090 END DO … … 1187 1144 ! update the momentum trends in u direction 1188 1145 1189 !!gm zdpdx1 = zcoef0 * r1_e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1190 zdpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 1146 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1191 1147 IF( lk_vvl ) THEN 1192 !!gm zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1193 zdpdx2 = zcoef0 / e1u(ji,jj) * & 1194 ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1148 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1149 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (sshn(ji+1,jj)-sshn(ji,jj)) ) 1195 1150 ELSE 1196 !!gm zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1197 zdpdx2 = zcoef0 / e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1151 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1198 1152 ENDIF 1199 1200 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * & 1201 &umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk)1153 !!gm Since umask(ji,:,:) = tmask(ji,:,:) * tmask(ji+1,:,:) by definition 1154 !!gm in the line below only * umask(ji,jj,jk) is needed !! 1155 ua(ji,jj,jk) = ua(ji,jj,jk) + (zdpdx1 + zdpdx2) * umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 1202 1156 ENDIF 1203 1157 … … 1247 1201 ! update the momentum trends in v direction 1248 1202 1249 !!gm zdpdy1 = zcoef0 * r1_e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1250 zdpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 1203 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1251 1204 IF( lk_vvl ) THEN 1252 !!gm zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1253 zdpdy2 = zcoef0 / e2v(ji,jj) * & 1205 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1254 1206 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (sshn(ji,jj+1)-sshn(ji,jj)) ) 1255 1207 ELSE 1256 !!gm zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1257 zdpdy2 = zcoef0 / e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1208 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1258 1209 ENDIF 1259 1260 va(ji,jj,jk) = va(ji,jj,jk) + (zdpdy1 + zdpdy2)*& 1261 &vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk)1210 !!gm Since vmask(:,jj,:) = tmask(:,jj,:) * tmask(:,jj+1,:) by definition 1211 !!gm in the line below only * vmask(ji,jj,jk) is needed !! 1212 va(ji,jj,jk) = va(ji,jj,jk) + ( zdpdy1 + zdpdy2 ) * vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 1262 1213 ENDIF 1263 1264 1265 END DO 1266 END DO 1267 END DO 1268 ! 1269 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1270 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1271 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1214 ! 1215 END DO 1216 END DO 1217 END DO 1218 ! 1219 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 1220 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 1221 CALL wrk_dealloc( jpi,jpj, zsshu_n, zsshv_n ) 1272 1222 ! 1273 1223 END SUBROUTINE hpg_prj 1274 1224 1275 1225 1276 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type)1226 SUBROUTINE cspline( fsp, xsp, asp, bsp, csp, dsp, polynomial_type ) 1277 1227 !!---------------------------------------------------------------------- 1278 1228 !! *** ROUTINE cspline *** … … 1297 1247 REAL(wp) :: zdf(size(fsp,3)) 1298 1248 !!---------------------------------------------------------------------- 1299 1249 ! 1250 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1300 1251 jpi = size(fsp,1) 1301 1252 jpj = size(fsp,2) 1302 1253 jpkm1 = size(fsp,3) - 1 1303 1304 1254 ! 1305 1255 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1306 1256 DO ji = 1, jpi … … 1335 1285 1336 1286 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1337 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2)1287 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1338 1288 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1339 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1340 & 0.5_wp * zdf(jpkm1 - 1) 1289 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1341 1290 1342 1291 DO jk = 1, jpkm1 - 1 … … 1361 1310 END DO 1362 1311 1363 ELSE IF (polynomial_type == 2) THEN ! Linear1312 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1364 1313 DO ji = 1, jpi 1365 1314 DO jj = 1, jpj … … 1392 1341 !! extrapolation is also permitted (no value limit) 1393 1342 !!---------------------------------------------------------------------- 1394 IMPLICIT NONE1395 1343 REAL(wp), INTENT(in) :: x, xl, xr, fl, fr 1396 1344 REAL(wp) :: f ! result of the interpolation (extrapolation) 1397 1345 REAL(wp) :: zdeltx 1398 1346 !!---------------------------------------------------------------------- 1399 1347 ! 1400 1348 zdeltx = xr - xl 1401 IF( abs(zdeltx) <= 10._wp * EPSILON(x)) THEN1402 f = 0.5_wp * (fl + fr)1349 IF( abs(zdeltx) <= 10._wp * EPSILON(x) ) THEN 1350 f = 0.5_wp * (fl + fr) 1403 1351 ELSE 1404 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx1352 f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 1405 1353 ENDIF 1406 1354 ! 1407 1355 END FUNCTION interp1 1408 1356 1409 1357 1410 FUNCTION interp2( x, a, b, c, d) RESULT(f)1358 FUNCTION interp2( x, a, b, c, d ) RESULT(f) 1411 1359 !!---------------------------------------------------------------------- 1412 1360 !! *** ROUTINE interp1 *** … … 1417 1365 !! 1418 1366 !!---------------------------------------------------------------------- 1419 IMPLICIT NONE1420 1367 REAL(wp), INTENT(in) :: x, a, b, c, d 1421 1368 REAL(wp) :: f ! value from the interpolation 1422 1369 !!---------------------------------------------------------------------- 1423 1370 ! 1424 1371 f = a + x* ( b + x * ( c + d * x ) ) 1425 1372 ! 1426 1373 END FUNCTION interp2 1427 1374 1428 1375 1429 FUNCTION interp3( x, a, b, c, d) RESULT(f)1376 FUNCTION interp3( x, a, b, c, d ) RESULT(f) 1430 1377 !!---------------------------------------------------------------------- 1431 1378 !! *** ROUTINE interp1 *** … … 1437 1384 !! 1438 1385 !!---------------------------------------------------------------------- 1439 IMPLICIT NONE1440 1386 REAL(wp), INTENT(in) :: x, a, b, c, d 1441 1387 REAL(wp) :: f ! value from the interpolation 1442 1388 !!---------------------------------------------------------------------- 1443 1389 ! 1444 1390 f = b + x * ( 2._wp * c + 3._wp * d * x) 1445 1391 ! 1446 1392 END FUNCTION interp3 1447 1393 1448 1394 1449 FUNCTION integ_spline( xl, xr, a, b, c, d) RESULT(f)1395 FUNCTION integ_spline( xl, xr, a, b, c, d ) RESULT(f) 1450 1396 !!---------------------------------------------------------------------- 1451 1397 !! *** ROUTINE interp1 *** … … 1456 1402 !! 1457 1403 !!---------------------------------------------------------------------- 1458 IMPLICIT NONE1459 1404 REAL(wp), INTENT(in) :: xl, xr, a, b, c, d 1460 1405 REAL(wp) :: za1, za2, za3 1461 1406 REAL(wp) :: f ! integration result 1462 1407 !!---------------------------------------------------------------------- 1463 1408 ! 1464 1409 za1 = 0.5_wp * b 1465 1410 za2 = c / 3.0_wp 1466 1411 za3 = 0.25_wp * d 1467 1412 ! 1468 1413 f = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 1469 1414 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1470 1415 ! 1471 1416 END FUNCTION integ_spline 1472 1417
Note: See TracChangeset
for help on using the changeset viewer.
