Changeset 5034 for branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 2015-01-15T14:48:42+01:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r4624 r5034 29 29 !!---------------------------------------------------------------------- 30 30 USE oce ! ocean dynamics and tracers 31 USE sbc_oce ! surface variable (only for the flag with ice shelf) 31 32 USE dom_oce ! ocean space and time domain 32 33 USE phycst ! physical constants 33 USE trdmod ! ocean dynamics trends 34 USE trdmod_oce ! ocean variables trends 34 USE trd_oce ! trends: ocean variables 35 USE trddyn ! trend manager: dynamics 36 ! 35 37 USE in_out_manager ! I/O manager 36 38 USE prtctl ! Print control 37 USE lbclnk ! lateral boundary condition 39 USE lbclnk ! lateral boundary condition 38 40 USE lib_mpp ! MPP library 41 USE eosbn2 ! compute density 39 42 USE wrk_nemo ! Memory Allocation 40 43 USE timing ! Timing … … 74 77 !! 75 78 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 76 !! - Save the trend(l_trddyn=T)79 !! - send trends to trd_dyn for futher diagnostics (l_trddyn=T) 77 80 !!---------------------------------------------------------------------- 78 81 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 99 102 ztrdu(:,:,:) = ua(:,:,:) - ztrdu(:,:,:) 100 103 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 101 CALL trd_ mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt )104 CALL trd_dyn( ztrdu, ztrdv, jpdyn_hpg, kt ) 102 105 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 103 106 ENDIF … … 175 178 IF( ln_hpg_prj ) ioptio = ioptio + 1 176 179 IF( ioptio /= 1 ) CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 180 IF( (ln_hpg_zco .OR. ln_hpg_zps .OR. ln_hpg_djc .OR. ln_hpg_prj ) .AND. nn_isf .NE. 0 ) & 181 & CALL ctl_stop( 'Only hpg_sco has been corrected to work with ice shelf cavity.' ) 177 182 ! 178 183 END SUBROUTINE dyn_hpg_init … … 315 320 316 321 ! partial steps correction at the last level (use gru & grv computed in zpshde.F90) 317 # if defined key_vectopt_loop318 jj = 1319 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)320 # else321 322 DO jj = 2, jpjm1 322 323 DO ji = 2, jpim1 323 # endif324 324 iku = mbku(ji,jj) 325 325 ikv = mbkv(ji,jj) … … 338 338 va (ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) ! add the new one to the general momentum trend 339 339 ENDIF 340 # if ! defined key_vectopt_loop 341 END DO 342 # endif 340 END DO 343 341 END DO 344 342 ! … … 368 366 INTEGER, INTENT(in) :: kt ! ocean time-step index 369 367 !! 370 INTEGER :: ji, jj, jk ! dummy loop indices 371 REAL(wp) :: zcoef0, zuap, zvap, znad ! temporary scalars 372 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 373 !!---------------------------------------------------------------------- 374 ! 375 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 376 ! 377 IF( kt == nit000 ) THEN 368 INTEGER :: ji, jj, jk, iku, ikv, ikt, iktp1i, iktp1j ! dummy loop indices 369 REAL(wp) :: zcoef0, zuap, zvap, znad, ze3wu, ze3wv, zuapint, zvapint, zhpjint, zhpiint, zdzwt, zdzwtjp1, zdzwtip1 ! temporary scalars 370 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj, zrhd 371 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztstop 372 REAL(wp), POINTER, DIMENSION(:,:) :: ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj 373 !!---------------------------------------------------------------------- 374 ! 375 CALL wrk_alloc( jpi,jpj, 2, ztstop) 376 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 377 CALL wrk_alloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 378 ! 379 IF( kt == nit000 ) THEN 378 380 IF(lwp) WRITE(numout,*) 379 381 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' … … 384 386 zcoef0 = - grav * 0.5_wp 385 387 ! To use density and not density anomaly 386 IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 387 ELSE ; znad = 0._wp ! Fixed volume 388 ENDIF 389 390 ! Surface value 388 ! IF ( lk_vvl ) THEN ; znad = 1._wp ! Variable volume 389 ! ELSE ; znad = 0._wp ! Fixed volume 390 ! ENDIF 391 znad=1._wp 392 ! iniitialised to 0. zhpi zhpi 393 zhpi(:,:,:)=0._wp ; zhpj(:,:,:)=0._wp 394 395 !================================================================================== 396 !=====Compute iceload and contribution of the half first wet layer ================= 397 !=================================================================================== 398 399 ! assume water displaced by the ice shelf is at T=-1.9 and S=34.4 (rude) 400 ztstop(:,:,1)=-1.9_wp ; ztstop(:,:,2)=34.4_wp 401 402 ! compute density of the water displaced by the ice shelf 403 zrhd = rhd ! save rhd 404 DO jk = 1, jpk 405 zdept(:,:)=gdept_1d(jk) 406 CALL eos(ztstop(:,:,:),zdept(:,:),rhd(:,:,jk)) 407 END DO 408 WHERE ( tmask(:,:,:) == 1._wp) 409 rhd(:,:,:) = zrhd(:,:,:) ! replace wet cell by the saved rhd 410 END WHERE 411 412 ! compute rhd at the ice/oce interface (ice shelf side) 413 CALL eos(ztstop,risfdep,zrhdtop_isf) 414 415 ! compute rhd at the ice/oce interface (ocean side) 416 DO ji=1,jpi 417 DO jj=1,jpj 418 ikt=mikt(ji,jj) 419 ztstop(ji,jj,1)=tsn(ji,jj,ikt,1) 420 ztstop(ji,jj,2)=tsn(ji,jj,ikt,2) 421 END DO 422 END DO 423 CALL eos(ztstop,risfdep,zrhdtop_oce) 424 ! 425 ! Surface value + ice shelf gradient 426 ! compute pressure due to ice shelf load (used to compute hpgi/j for all the level from 1 to miku/v) 427 ziceload = 0._wp 428 DO jj = 1, jpj 429 DO ji = 1, jpi ! vector opt. 430 ikt=mikt(ji,jj) 431 ziceload(ji,jj) = ziceload(ji,jj) + (znad + rhd(ji,jj,1) ) * fse3w(ji,jj,1) * (1._wp - tmask(ji,jj,1)) 432 DO jk=2,ikt-1 433 ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + rhd(ji,jj,jk-1) + rhd(ji,jj,jk)) * fse3w(ji,jj,jk) & 434 & * (1._wp - tmask(ji,jj,jk)) 435 END DO 436 IF (ikt .GE. 2) ziceload(ji,jj) = ziceload(ji,jj) + (2._wp * znad + zrhdtop_isf(ji,jj) + rhd(ji,jj,ikt-1)) & 437 & * ( risfdep(ji,jj) - gdept_1d(ikt-1) ) 438 END DO 439 END DO 440 riceload(:,:) = 0.0_wp ; riceload(:,:)=ziceload(:,:) ! need to be saved for diaar5 441 ! compute zp from z=0 to first T wet point (correction due to zps not yet applied) 391 442 DO jj = 2, jpjm1 392 443 DO ji = fs_2, fs_jpim1 ! vector opt. 393 ! hydrostatic pressure gradient along s-surfaces 394 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( fse3w(ji+1,jj ,1) * ( znad + rhd(ji+1,jj ,1) ) & 395 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 396 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( fse3w(ji ,jj+1,1) * ( znad + rhd(ji ,jj+1,1) ) & 397 & - fse3w(ji ,jj ,1) * ( znad + rhd(ji ,jj ,1) ) ) 398 ! s-coordinate pressure gradient correction 444 ikt=mikt(ji,jj) ; iktp1i=mikt(ji+1,jj); iktp1j=mikt(ji,jj+1) 445 ! hydrostatic pressure gradient along s-surfaces and ice shelf pressure 446 ! we assume ISF is in isostatic equilibrium 447 zhpi(ji,jj,1) = zcoef0 / e1u(ji,jj) * ( 0.5_wp * fse3w(ji+1,jj ,iktp1i) & 448 & * ( 2._wp * znad + rhd(ji+1,jj ,iktp1i) + zrhdtop_oce(ji+1,jj ) ) & 449 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 450 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 451 & + ( ziceload(ji+1,jj) - ziceload(ji,jj)) ) 452 zhpj(ji,jj,1) = zcoef0 / e2v(ji,jj) * ( 0.5_wp * fse3w(ji ,jj+1,iktp1j) & 453 & * ( 2._wp * znad + rhd(ji ,jj+1,iktp1j) + zrhdtop_oce(ji ,jj+1) ) & 454 & - 0.5_wp * fse3w(ji ,jj ,ikt ) & 455 & * ( 2._wp * znad + rhd(ji ,jj ,ikt ) + zrhdtop_oce(ji ,jj ) ) & 456 & + ( ziceload(ji,jj+1) - ziceload(ji,jj) ) ) 457 ! s-coordinate pressure gradient correction (=0 if z coordinate) 399 458 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) + 2._wp * znad ) & 400 459 & * ( fsde3w(ji+1,jj,1) - fsde3w(ji,jj,1) ) / e1u(ji,jj) … … 402 461 & * ( fsde3w(ji,jj+1,1) - fsde3w(ji,jj,1) ) / e2v(ji,jj) 403 462 ! add to the general momentum trend 404 ua(ji,jj,1) = ua(ji,jj,1) + zhpi(ji,jj,1) + zuap 405 va(ji,jj,1) = va(ji,jj,1) + zhpj(ji,jj,1) + zvap 406 END DO 407 END DO 408 409 ! interior value (2=<jk=<jpkm1) 410 DO jk = 2, jpkm1 411 DO jj = 2, jpjm1 412 DO ji = fs_2, fs_jpim1 ! vector opt. 463 ua(ji,jj,1) = ua(ji,jj,1) + (zhpi(ji,jj,1) + zuap) * umask(ji,jj,1) 464 va(ji,jj,1) = va(ji,jj,1) + (zhpj(ji,jj,1) + zvap) * vmask(ji,jj,1) 465 END DO 466 END DO 467 !================================================================================== 468 !===== Compute partial cell contribution for the top cell ========================= 469 !================================================================================== 470 DO jj = 2, jpjm1 471 DO ji = fs_2, fs_jpim1 ! vector opt. 472 iku = miku(ji,jj) ; 473 zpshpi(ji,jj)=0.0_wp ; zpshpj(ji,jj)=0.0_wp 474 ze3wu = (gdepw_0(ji+1,jj,iku+1) - gdept_0(ji+1,jj,iku)) - (gdepw_0(ji,jj,iku+1) - gdept_0(ji,jj,iku)) 475 ! u direction 476 IF ( iku .GT. 1 ) THEN 477 ! case iku 478 zhpi(ji,jj,iku) = zcoef0 / e1u(ji,jj) * ze3wu & 479 & * ( rhd (ji+1,jj,iku) + rhd (ji,jj,iku) & 480 & + SIGN(1._wp,ze3wu) * grui(ji,jj) + 2._wp * znad ) 481 ! corrective term ( = 0 if z coordinate ) 482 zuap = -zcoef0 * ( arui(ji,jj) + 2._wp * znad ) * gzui(ji,jj) / e1u(ji,jj) 483 ! zhpi will be added in interior loop 484 ua(ji,jj,iku) = ua(ji,jj,iku) + zuap 485 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 486 IF (mbku(ji,jj) == iku + 1) zpshpi(ji,jj) = zhpi(ji,jj,iku) 487 488 ! case iku + 1 (remove the zphi term added in the interior loop and compute the one corrected for zps) 489 zhpiint = zcoef0 / e1u(ji,jj) & 490 & * ( fse3w(ji+1,jj ,iku+1) * ( (rhd(ji+1,jj,iku+1) + znad) & 491 & + (rhd(ji+1,jj,iku ) + znad) ) * tmask(ji+1,jj,iku) & 492 & - fse3w(ji ,jj ,iku+1) * ( (rhd(ji ,jj,iku+1) + znad) & 493 & + (rhd(ji ,jj,iku ) + znad) ) * tmask(ji ,jj,iku) ) 494 zhpi(ji,jj,iku+1) = zcoef0 / e1u(ji,jj) * ge3rui(ji,jj) - zhpiint 495 END IF 496 497 ! v direction 498 ikv = mikv(ji,jj) 499 ze3wv = (gdepw_0(ji,jj+1,ikv+1) - gdept_0(ji,jj+1,ikv)) - (gdepw_0(ji,jj,ikv+1) - gdept_0(ji,jj,ikv)) 500 IF ( ikv .GT. 1 ) THEN 501 ! case ikv 502 zhpj(ji,jj,ikv) = zcoef0 / e2v(ji,jj) * ze3wv & 503 & * ( rhd(ji,jj+1,ikv) + rhd (ji,jj,ikv) & 504 & + SIGN(1._wp,ze3wv) * grvi(ji,jj) + 2._wp * znad ) 505 ! corrective term ( = 0 if z coordinate ) 506 zvap = -zcoef0 * ( arvi(ji,jj) + 2._wp * znad ) * gzvi(ji,jj) / e2v(ji,jj) 507 ! zhpi will be added in interior loop 508 va(ji,jj,ikv) = va(ji,jj,ikv) + zvap 509 ! in case of 2 cell water column, need to save the pressure gradient to compute the bottom pressure 510 IF (mbkv(ji,jj) == ikv + 1) zpshpj(ji,jj) = zhpj(ji,jj,ikv) 511 512 ! case ikv + 1 (remove the zphj term added in the interior loop and compute the one corrected for zps) 513 zhpjint = zcoef0 / e2v(ji,jj) & 514 & * ( fse3w(ji ,jj+1,ikv+1) * ( (rhd(ji,jj+1,ikv+1) + znad) & 515 & + (rhd(ji,jj+1,ikv ) + znad) ) * tmask(ji,jj+1,ikv) & 516 & - fse3w(ji ,jj ,ikv+1) * ( (rhd(ji,jj ,ikv+1) + znad) & 517 & + (rhd(ji,jj ,ikv ) + znad) ) * tmask(ji,jj ,ikv) ) 518 zhpj(ji,jj,ikv+1) = zcoef0 / e2v(ji,jj) * ge3rvi(ji,jj) - zhpjint 519 END IF 520 END DO 521 END DO 522 523 !================================================================================== 524 !===== Compute interior value ===================================================== 525 !================================================================================== 526 527 DO jj = 2, jpjm1 528 DO ji = fs_2, fs_jpim1 ! vector opt. 529 iku=miku(ji,jj); ikv=mikv(ji,jj) 530 DO jk = 2, jpkm1 413 531 ! hydrostatic pressure gradient along s-surfaces 414 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + zcoef0 / e1u(ji,jj) & 415 & * ( fse3w(ji+1,jj,jk) * ( rhd(ji+1,jj,jk) + rhd(ji+1,jj,jk-1) + 2*znad ) & 416 & - fse3w(ji ,jj,jk) * ( rhd(ji ,jj,jk) + rhd(ji ,jj,jk-1) + 2*znad ) ) 417 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) + zcoef0 / e2v(ji,jj) & 418 & * ( fse3w(ji,jj+1,jk) * ( rhd(ji,jj+1,jk) + rhd(ji,jj+1,jk-1) + 2*znad ) & 419 & - fse3w(ji,jj ,jk) * ( rhd(ji,jj, jk) + rhd(ji,jj ,jk-1) + 2*znad ) ) 532 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 533 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) + zhpi(ji,jj,jk-1) & 534 & + zcoef0 / e1u(ji,jj) & 535 & * ( fse3w(ji+1,jj ,jk) * ( (rhd(ji+1,jj,jk ) + znad) & 536 & + (rhd(ji+1,jj,jk-1) + znad) ) * tmask(ji+1,jj,jk-1) & 537 & - fse3w(ji ,jj ,jk) * ( (rhd(ji ,jj,jk ) + znad) & 538 & + (rhd(ji ,jj,jk-1) + znad) ) * tmask(ji ,jj,jk-1) ) 420 539 ! s-coordinate pressure gradient correction 421 zuap = -zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 422 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) 423 zvap = -zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 424 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) 540 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 541 zuap = - zcoef0 * ( rhd (ji+1,jj ,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 542 & * ( fsde3w(ji+1,jj ,jk) - fsde3w(ji,jj,jk) ) / e1u(ji,jj) * umask(ji,jj,jk-1) 543 ua(ji,jj,jk) = ua(ji,jj,jk) + ( zhpi(ji,jj,jk) + zuap) * umask(ji,jj,jk) 544 545 ! hydrostatic pressure gradient along s-surfaces 546 ! zhpi is masked for the first wet cell (contribution already done in the upper bloc) 547 zhpj(ji,jj,jk) = zhpj(ji,jj,jk) + zhpj(ji,jj,jk-1) & 548 & + zcoef0 / e2v(ji,jj) & 549 & * ( fse3w(ji ,jj+1,jk) * ( (rhd(ji,jj+1,jk ) + znad) & 550 & + (rhd(ji,jj+1,jk-1) + znad) ) * tmask(ji,jj+1,jk-1) & 551 & - fse3w(ji ,jj ,jk) * ( (rhd(ji,jj ,jk ) + znad) & 552 & + (rhd(ji,jj ,jk-1) + znad) ) * tmask(ji,jj ,jk-1) ) 553 ! s-coordinate pressure gradient correction 554 ! corrective term, we mask this term for the first wet level beneath the ice shelf (contribution done in the upper bloc) 555 zvap = - zcoef0 * ( rhd (ji ,jj+1,jk) + rhd (ji,jj,jk) + 2._wp * znad ) & 556 & * ( fsde3w(ji ,jj+1,jk) - fsde3w(ji,jj,jk) ) / e2v(ji,jj) * vmask(ji,jj,jk-1) 425 557 ! add to the general momentum trend 426 ua(ji,jj,jk) = ua(ji,jj,jk) + zhpi(ji,jj,jk) + zuap 427 va(ji,jj,jk) = va(ji,jj,jk) + zhpj(ji,jj,jk) + zvap 558 va(ji,jj,jk) = va(ji,jj,jk) + ( zhpj(ji,jj,jk) + zvap ) * vmask(ji,jj,jk) 428 559 END DO 429 560 END DO 430 561 END DO 431 ! 432 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 562 563 !================================================================================== 564 !===== Compute bottom cell contribution (partial cell) ============================ 565 !================================================================================== 566 567 # if defined key_vectopt_loop 568 jj = 1 569 DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling) 570 # else 571 DO jj = 2, jpjm1 572 DO ji = 2, jpim1 573 # endif 574 iku = mbku(ji,jj) 575 ikv = mbkv(ji,jj) 576 577 IF (iku .GT. 1) THEN 578 ! remove old value (interior case) 579 zuap = -zcoef0 * ( rhd (ji+1,jj ,iku) + rhd (ji,jj,iku) + 2._wp * znad ) & 580 & * ( fsde3w(ji+1,jj ,iku) - fsde3w(ji,jj,iku) ) / e1u(ji,jj) 581 ua(ji,jj,iku) = ua(ji,jj,iku) - zhpi(ji,jj,iku) - zuap 582 ! put new value 583 ! -zpshpi to avoid double contribution of the partial step in the top layer 584 zuap = -zcoef0 * ( aru(ji,jj) + 2._wp * znad ) * gzu(ji,jj) / e1u(ji,jj) 585 zhpi(ji,jj,iku) = zhpi(ji,jj,iku-1) + zcoef0 / e1u(ji,jj) * ge3ru(ji,jj) - zpshpi(ji,jj) 586 ua(ji,jj,iku) = ua(ji,jj,iku) + zhpi(ji,jj,iku) + zuap 587 END IF 588 ! v direction 589 IF (ikv .GT. 1) THEN 590 ! remove old value (interior case) 591 zvap = -zcoef0 * ( rhd (ji ,jj+1,ikv) + rhd (ji,jj,ikv) + 2._wp * znad ) & 592 & * ( fsde3w(ji ,jj+1,ikv) - fsde3w(ji,jj,ikv) ) / e2v(ji,jj) 593 va(ji,jj,ikv) = va(ji,jj,ikv) - zhpj(ji,jj,ikv) - zvap 594 ! put new value 595 ! -zpshpj to avoid double contribution of the partial step in the top layer 596 zvap = -zcoef0 * ( arv(ji,jj) + 2._wp * znad ) * gzv(ji,jj) / e2v(ji,jj) 597 zhpj(ji,jj,ikv) = zhpj(ji,jj,ikv-1) + zcoef0 / e2v(ji,jj) * ge3rv(ji,jj) - zpshpj(ji,jj) 598 va(ji,jj,ikv) = va(ji,jj,ikv) + zhpj(ji,jj,ikv) + zvap 599 END IF 600 # if ! defined key_vectopt_loop 601 END DO 602 # endif 603 END DO 604 605 ! set back to original density value into the ice shelf cell (maybe useless because it is masked) 606 rhd = zrhd 607 ! 608 CALL wrk_dealloc( jpi,jpj,2, ztstop) 609 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj, zrhd) 610 CALL wrk_dealloc( jpi,jpj, ze3w, zp, zrhdtop_isf, zrhdtop_oce, ziceload, zdept, zpshpi, zpshpj) 433 611 ! 434 612 END SUBROUTINE hpg_sco 613 435 614 436 615 SUBROUTINE hpg_djc( kt ) … … 671 850 !! 672 851 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 673 !! - Save the trend (l_trddyn=T)674 !!675 852 !!---------------------------------------------------------------------- 676 853 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear … … 724 901 725 902 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 726 DO jj = 1, jpj; DO ji = 1, jpi 727 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 728 END DO ; END DO 729 730 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 731 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 732 END DO ; END DO ; END DO 733 734 fsp(:,:,:) = zrhh(:,:,:) 903 DO jj = 1, jpj 904 DO ji = 1, jpi 905 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 906 END DO 907 END DO 908 909 DO jk = 2, jpk 910 DO jj = 1, jpj 911 DO ji = 1, jpi 912 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 913 END DO 914 END DO 915 END DO 916 917 fsp(:,:,:) = zrhh (:,:,:) 735 918 xsp(:,:,:) = zdept(:,:,:) 736 919 … … 933 1116 END SUBROUTINE hpg_prj 934 1117 1118 935 1119 SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 936 1120 !!---------------------------------------------------------------------- … … 940 1124 !! 941 1125 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 1126 !! 942 1127 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 943 !!944 1128 !!---------------------------------------------------------------------- 945 1129 IMPLICIT NONE … … 949 1133 INTEGER, INTENT(in) :: polynomial_type ! 1: cubic spline 950 1134 ! 2: Linear 951 952 ! Local Variables 1135 ! 953 1136 INTEGER :: ji, jj, jk ! dummy loop indices 954 1137 INTEGER :: jpi, jpj, jpkm1 … … 1040 1223 ENDIF 1041 1224 1042 1043 1225 END SUBROUTINE cspline 1044 1226 … … 1050 1232 !! ** Purpose : 1-d linear interpolation 1051 1233 !! 1052 !! ** Method : 1053 !! interpolation is straight forward 1234 !! ** Method : interpolation is straight forward 1054 1235 !! extrapolation is also permitted (no value limit) 1055 !!1056 1236 !!---------------------------------------------------------------------- 1057 1237 IMPLICIT NONE … … 1070 1250 END FUNCTION interp1 1071 1251 1252 1072 1253 FUNCTION interp2(x, a, b, c, d) RESULT(f) 1073 1254 !!---------------------------------------------------------------------- … … 1133 1314 END FUNCTION integ_spline 1134 1315 1135 1136 1316 !!====================================================================== 1137 1317 END MODULE dynhpg
Note: See TracChangeset
for help on using the changeset viewer.
