Changeset 14852 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
- Timestamp:
- 2021-05-12T15:05:29+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/OCE/DYN/dynhpg.F90
r14789 r14852 266 266 INTEGER :: ji, jj, jk ! dummy loop indices 267 267 REAL(wp) :: zcoef0, zcoef1 ! temporary scalars 268 REAL(wp), DIMENSION(jpi,jpj) :: zhpi, zhpj 269 !!---------------------------------------------------------------------- 270 ! 271 IF( kt == nit000 ) THEN 272 IF(lwp) WRITE(numout,*) 273 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 274 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 268 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj 269 !!---------------------------------------------------------------------- 270 ! 271 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 272 IF( kt == nit000 ) THEN 273 IF(lwp) WRITE(numout,*) 274 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 275 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 276 ENDIF 275 277 ENDIF 276 278 ! … … 318 320 INTEGER :: iku, ikv ! temporary integers 319 321 REAL(wp) :: zcoef0, zcoef1, zcoef2, zcoef3 ! temporary scalars 320 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 321 REAL(wp), DIMENSION(jpi,jpj,jpts) :: zgtsu, zgtsv 322 REAL(wp), DIMENSION(jpi,jpj) :: zgru, zgrv 323 !!---------------------------------------------------------------------- 324 ! 325 IF( kt == nit000 ) THEN 326 IF(lwp) WRITE(numout,*) 327 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 328 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 322 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 323 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zgtsu, zgtsv 324 REAL(wp), DIMENSION(A2D(nn_hls) ) :: zgru, zgrv 325 !!---------------------------------------------------------------------- 326 ! 327 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 328 IF( kt == nit000 ) THEN 329 IF(lwp) WRITE(numout,*) 330 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 331 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 332 ENDIF 329 333 ENDIF 330 334 … … 410 414 REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars 411 415 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 412 REAL(wp), DIMENSION( jpi,jpj,jpk):: zhpi, zhpj416 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 413 417 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 414 418 !!---------------------------------------------------------------------- 415 419 ! 416 IF( ln_wd_il ) ALLOCATE(zcpx(jpi,jpj), zcpy(jpi,jpj)) 417 ! 418 IF( kt == nit000 ) THEN 419 IF(lwp) WRITE(numout,*) 420 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 421 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' 420 IF( ln_wd_il ) ALLOCATE(zcpx(A2D(nn_hls)), zcpy(A2D(nn_hls))) 421 ! 422 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 423 IF( kt == nit000 ) THEN 424 IF(lwp) WRITE(numout,*) 425 IF(lwp) WRITE(numout,*) 'dyn:hpg_sco : hydrostatic pressure gradient trend' 426 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, OCE original scheme used' 427 ENDIF 422 428 ENDIF 423 429 ! … … 462 468 END IF 463 469 END_2D 464 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )465 470 END IF 466 471 ! … … 548 553 REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars 549 554 REAL(wp) :: zcoef0, zuap, zvap ! - - 550 REAL(wp), DIMENSION( jpi,jpj,jpk ) :: zhpi, zhpj551 REAL(wp), DIMENSION( jpi,jpj,jpts) :: zts_top552 REAL(wp), DIMENSION( jpi,jpj) :: zrhdtop_oce555 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 556 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts_top 557 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrhdtop_oce 553 558 !!---------------------------------------------------------------------- 554 559 ! … … 560 565 ! compute rhd at the ice/oce interface (ocean side) 561 566 ! usefull to reduce residual current in the test case ISOMIP with no melting 562 DO ji = 1, jpi 563 DO jj = 1, jpj 564 ikt = mikt(ji,jj) 565 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 566 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 567 END DO 568 END DO 567 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 568 ikt = mikt(ji,jj) 569 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 570 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 571 END_2D 569 572 CALL eos( zts_top, risfdep, zrhdtop_oce ) 570 573 … … 636 639 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 637 640 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 638 REAL(wp) :: z_grav_10, z1_12 641 REAL(wp) :: z_grav_10, z1_12, z1_cff 639 642 REAL(wp) :: cffu, cffx ! " " 640 643 REAL(wp) :: cffv, cffy ! " " 641 644 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 642 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zhpi, zhpj643 644 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz')645 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz')646 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho')647 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho')648 REAL(wp), DIMENSION( jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals649 REAL(wp), DIMENSION( jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays645 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 646 647 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 648 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 649 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 650 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 651 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 652 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 650 653 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 651 654 !!---------------------------------------------------------------------- 652 655 ! 653 656 IF( ln_wd_il ) THEN 654 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )657 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 655 658 DO_2D( 0, 0, 0, 0 ) 656 659 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & … … 689 692 END IF 690 693 END_2D 691 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )692 694 END IF 693 695 694 IF( kt == nit000 ) THEN 695 IF(lwp) WRITE(numout,*) 696 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 697 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 696 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 697 IF( kt == nit000 ) THEN 698 IF(lwp) WRITE(numout,*) 699 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 700 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 701 ENDIF 698 702 ENDIF 699 703 … … 723 727 zdz_k (:,:,:) = 0._wp 724 728 725 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 726 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 727 IF( cffw > zep) THEN 728 zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 729 ENDIF 729 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 730 cffw = MAX( 2._wp * zdrhoz(ji,jj,jk) * zdrhoz(ji,jj,jk+1), 0._wp ) 731 z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) 732 zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 730 733 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 731 734 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) … … 737 740 738 741 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 739 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd (:,:,2) - rhd (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 740 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 742 DO_2D( 1, 1, 1, 1 ) 743 zdrho_k(ji,jj,1) = aco_bc_vrt * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_vrt * zdrho_k(ji,jj,2) 744 zdz_k (ji,jj,1) = aco_bc_vrt * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_vrt * zdz_k (ji,jj,2) 745 END_2D 741 746 742 747 DO_2D( 1, 1, 1, 1 ) … … 785 790 ! 5. compute and store elementary horizontal differences in provisional arrays 786 791 !---------------------------------------------------------------------------------------- 787 788 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 789 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 790 zdzx (ji,jj,jk) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 791 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 792 zdzy (ji,jj,jk) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 793 END_3D 794 795 CALL lbc_lnk( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 792 zdrhox(:,:,:) = 0._wp 793 zdzx (:,:,:) = 0._wp 794 zdrhoy(:,:,:) = 0._wp 795 zdzy (:,:,:) = 0._wp 796 797 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 798 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji ,jj ,jk) 799 zdzx (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji+1,jj ,jk) 800 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji ,jj ,jk) 801 zdzy (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji ,jj+1,jk) 802 END_3D 803 804 IF( nn_hls == 1 ) CALL lbc_lnk( 'dynhpg', zdrhox, 'U', -1._wp, zdzx, 'U', -1._wp, zdrhoy, 'V', -1._wp, zdzy, 'V', -1._wp ) 796 805 797 806 !------------------------------------------------------------------------- … … 800 809 801 810 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 802 cffu = 2._wp * zdrhox(ji-1,jj ,jk) * zdrhox(ji,jj,jk ) 803 IF( cffu > zep ) THEN 804 zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 805 ELSE 806 zdrho_i(ji,jj,jk ) = 0._wp 807 ENDIF 808 809 cffx = 2._wp * zdzx (ji-1,jj ,jk) * zdzx (ji,jj,jk ) 810 IF( cffx > zep ) THEN 811 zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 812 ELSE 813 zdz_i(ji,jj,jk) = 0._wp 814 ENDIF 815 816 cffv = 2._wp * zdrhoy(ji ,jj-1,jk) * zdrhoy(ji,jj,jk ) 817 IF( cffv > zep ) THEN 818 zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 819 ELSE 820 zdrho_j(ji,jj,jk) = 0._wp 821 ENDIF 822 823 cffy = 2._wp * zdzy (ji ,jj-1,jk) * zdzy (ji,jj,jk ) 824 IF( cffy > zep ) THEN 825 zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 826 ELSE 827 zdz_j(ji,jj,jk) = 0._wp 828 ENDIF 811 cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp ) 812 z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) 813 zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 814 815 cffx = MAX( 2._wp * zdzx(ji-1,jj,jk) * zdzx(ji,jj,jk), 0._wp ) 816 z1_cff = zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) 817 zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 818 819 cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp ) 820 z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) 821 zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 822 823 cffy = MAX( 2._wp * zdzy(ji,jj-1,jk) * zdzy(ji,jj,jk), 0._wp ) 824 z1_cff = zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) 825 zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 829 826 END_3D 830 827 … … 840 837 zz_drho_j(:,:) = zdrho_j(:,:,jk) 841 838 zz_dz_j (:,:) = zdz_j (:,:,jk) 842 DO_2D( 0, 1, 0, 1) 843 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 844 IF (ji < jpi) THEN 845 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 846 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 847 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 848 END IF 839 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 840 DO_2D( 0, 0, 0, 1 ) 841 IF ( umask(ji,jj,jk) > 0.5_wp .AND. umask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 842 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 843 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 849 844 END IF 850 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj)851 IF (ji > 2) THEN852 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN853 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk)854 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i(ji-1,jj,jk)855 END IF845 END_2D 846 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 847 DO_2D( -1, 1, 0, 1 ) 848 IF ( umask(ji,jj,jk) < 0.5_wp .AND. umask(ji-1,jj,jk) > 0.5_wp .AND. umask(ji-2,jj,jk) > 0.5_wp) THEN 849 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) 850 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) 856 851 END IF 857 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi)858 IF (jj < jpj) THEN859 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN860 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk)861 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j(ji,jj+1,jk)862 END IF863 END IF 864 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi)865 IF (jj > 2) THEN866 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN867 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk)868 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j(ji,jj-1,jk)869 END IF852 END_2D 853 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 854 DO_2D( 0, 1, 0, 0 ) 855 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. vmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 856 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 857 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) 858 END IF 859 END_2D 860 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 861 DO_2D( 0, 1, -1, 1 ) 862 IF ( vmask(ji,jj,jk) < 0.5_wp .AND. vmask(ji,jj-1,jk) > 0.5_wp .AND. vmask(ji,jj-2,jk) > 0.5_wp) THEN 863 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd (ji,jj,jk) - rhd (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 864 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w(ji,jj,jk) + gde3w(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) 870 865 END IF 871 866 END_2D … … 974 969 REAL(wp) :: zrhdt1 975 970 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 976 REAL(wp), DIMENSION( jpi,jpj) :: zpgu, zpgv ! 2D workspace977 REAL(wp), DIMENSION( jpi,jpj) :: zsshu_n, zsshv_n978 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zdept, zrhh979 REAL(wp), DIMENSION( jpi,jpj,jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp971 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpgu, zpgv ! 2D workspace 972 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsshu_n, zsshv_n 973 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdept, zrhh 974 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 980 975 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 981 976 !!---------------------------------------------------------------------- 982 977 ! 983 IF( kt == nit000 ) THEN 984 IF(lwp) WRITE(numout,*) 985 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 986 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 978 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 979 IF( kt == nit000 ) THEN 980 IF(lwp) WRITE(numout,*) 981 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 982 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 983 ENDIF 987 984 ENDIF 988 985 … … 1001 998 ! 1002 999 IF( ln_wd_il ) THEN 1003 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )1000 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 1004 1001 DO_2D( 0, 0, 0, 0 ) 1005 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1006 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1007 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) & 1008 & > rn_wdmin1 + rn_wdmin2 1009 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. ( & 1010 & MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1011 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1012 1013 IF(ll_tmp1) THEN 1014 zcpx(ji,jj) = 1.0_wp 1015 ELSE IF(ll_tmp2) THEN 1016 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1017 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1018 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1019 1020 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1021 ELSE 1022 zcpx(ji,jj) = 0._wp 1023 END IF 1024 1025 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1026 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1027 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) & 1028 & > rn_wdmin1 + rn_wdmin2 1029 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. ( & 1030 & MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1031 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1032 1033 IF(ll_tmp1) THEN 1034 zcpy(ji,jj) = 1.0_wp 1035 ELSE IF(ll_tmp2) THEN 1036 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1037 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1038 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1039 zcpy(ji,jj) = max(min( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1040 1002 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1003 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1004 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) > & 1005 & rn_wdmin1 + rn_wdmin2 1006 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. & 1007 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1008 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1009 1010 IF(ll_tmp1) THEN 1011 zcpx(ji,jj) = 1.0_wp 1012 ELSE IF(ll_tmp2) THEN 1013 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1014 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1015 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1016 zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1017 ELSE 1018 zcpx(ji,jj) = 0._wp 1019 END IF 1020 1021 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1022 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1023 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) > & 1024 & rn_wdmin1 + rn_wdmin2 1025 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. & 1026 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1027 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1028 1029 IF(ll_tmp1) THEN 1030 zcpy(ji,jj) = 1.0_wp 1031 ELSE IF(ll_tmp2) THEN 1032 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1033 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1034 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1035 zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1041 1036 ELSE 1042 1037 zcpy(ji,jj) = 0._wp 1043 1038 ENDIF 1044 1039 END_2D 1045 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )1046 1040 ENDIF 1047 1041 1048 1042 ! Clean 3-D work arrays 1049 1043 zhpi(:,:,:) = 0._wp 1050 zrhh(:,:,:) = rhd( :,:,:)1044 zrhh(:,:,:) = rhd(A2D(nn_hls),:) 1051 1045 1052 1046 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 1053 1047 DO_2D( 1, 1, 1, 1 ) 1054 jk = mbkt(ji,jj)1055 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp1056 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk)1057 ELSEIF( jk < jpkm1 ) THEN1058 DO jkk = jk+1, jpk1059 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), &1060 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2))1061 END DO1062 ENDIF1048 jk = mbkt(ji,jj) 1049 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1050 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 1051 ELSEIF( jk < jpkm1 ) THEN 1052 DO jkk = jk+1, jpk 1053 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1054 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2)) 1055 END DO 1056 ENDIF 1063 1057 END_2D 1064 1058 … … 1082 1076 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1083 1077 DO_2D( 0, 1, 0, 1 ) 1084 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), &1085 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm)1086 1087 ! assuming linear profile across the top half surface layer1088 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt11078 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1079 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1080 1081 ! assuming linear profile across the top half surface layer 1082 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1089 1083 END_2D 1090 1084 1091 1085 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 1092 1086 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 1093 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + &1094 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), &1095 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), &1096 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) )1087 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1088 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1089 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1090 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1097 1091 END_3D 1098 1092 … … 1107 1101 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1108 1102 !!gm not this: 1109 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1110 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1111 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1112 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1113 END_2D 1114 1115 CALL lbc_lnk ('dynhpg', zsshu_n, 'U', 1.0_wp, zsshv_n, 'V', 1.0_wp ) 1103 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1104 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1105 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1106 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1107 END_2D 1116 1108 1117 1109 DO_2D( 0, 0, 0, 0 ) 1118 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) )1119 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) )1110 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 1111 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 1120 1112 END_2D 1121 1113 1122 1114 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 1123 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm)1124 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm)1115 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1116 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1125 1117 END_3D 1126 1118 1127 1119 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1128 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm)1129 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm)1120 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1121 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1130 1122 END_3D 1131 1123 1132 1124 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1133 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )1134 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) )1135 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )1136 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) )1125 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1126 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1127 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1128 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1137 1129 END_3D 1138 1130 1139 1131 1140 1132 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 1141 zpwes = 0._wp; zpwed = 0._wp 1142 zpnss = 0._wp; zpnsd = 0._wp 1143 zuijk = zu(ji,jj,jk) 1144 zvijk = zv(ji,jj,jk) 1145 1146 !!!!! for u equation 1147 IF( jk <= mbku(ji,jj) ) THEN 1148 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1149 jis = ji + 1; jid = ji 1150 ELSE 1151 jis = ji; jid = ji +1 1133 zpwes = 0._wp; zpwed = 0._wp 1134 zpnss = 0._wp; zpnsd = 0._wp 1135 zuijk = zu(ji,jj,jk) 1136 zvijk = zv(ji,jj,jk) 1137 1138 !!!!! for u equation 1139 IF( jk <= mbku(ji,jj) ) THEN 1140 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1141 jis = ji + 1; jid = ji 1142 ELSE 1143 jis = ji; jid = ji +1 1144 ENDIF 1145 1146 ! integrate the pressure on the shallow side 1147 jk1 = jk 1148 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1149 IF( jk1 == mbku(ji,jj) ) THEN 1150 zuijk = -zdept(jis,jj,jk1) 1151 EXIT 1152 ENDIF 1153 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1154 zpwes = zpwes + & 1155 integ_spline(zdept(jis,jj,jk1), zdeps, & 1156 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1157 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1158 jk1 = jk1 + 1 1159 END DO 1160 1161 ! integrate the pressure on the deep side 1162 jk1 = jk 1163 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1164 IF( jk1 == 1 ) THEN 1165 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1166 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1167 bsp(jid,jj,1) , csp(jid,jj,1), & 1168 dsp(jid,jj,1)) * zdeps 1169 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1170 EXIT 1171 ENDIF 1172 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1173 zpwed = zpwed + & 1174 integ_spline(zdeps, zdept(jid,jj,jk1), & 1175 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1176 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1177 jk1 = jk1 - 1 1178 END DO 1179 1180 ! update the momentum trends in u direction 1181 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1182 IF( .NOT.ln_linssh ) THEN 1183 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1184 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1185 ELSE 1186 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1187 ENDIF 1188 IF( ln_wd_il ) THEN 1189 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1190 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1191 ENDIF 1192 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 1152 1193 ENDIF 1153 1194 1154 ! integrate the pressure on the shallow side 1155 jk1 = jk 1156 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1157 IF( jk1 == mbku(ji,jj) ) THEN 1158 zuijk = -zdept(jis,jj,jk1) 1159 EXIT 1160 ENDIF 1161 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1162 zpwes = zpwes + & 1163 integ_spline(zdept(jis,jj,jk1), zdeps, & 1164 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1165 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1166 jk1 = jk1 + 1 1167 END DO 1168 1169 ! integrate the pressure on the deep side 1170 jk1 = jk 1171 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1172 IF( jk1 == 1 ) THEN 1173 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1174 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1175 bsp(jid,jj,1), csp(jid,jj,1), & 1176 dsp(jid,jj,1)) * zdeps 1177 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1178 EXIT 1179 ENDIF 1180 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1181 zpwed = zpwed + & 1182 integ_spline(zdeps, zdept(jid,jj,jk1), & 1183 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1184 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1185 jk1 = jk1 - 1 1186 END DO 1187 1188 ! update the momentum trends in u direction 1189 1190 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1191 IF( .NOT.ln_linssh ) THEN 1192 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1193 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1194 ELSE 1195 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1195 !!!!! for v equation 1196 IF( jk <= mbkv(ji,jj) ) THEN 1197 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 1198 jjs = jj + 1; jjd = jj 1199 ELSE 1200 jjs = jj ; jjd = jj + 1 1201 ENDIF 1202 1203 ! integrate the pressure on the shallow side 1204 jk1 = jk 1205 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 1206 IF( jk1 == mbkv(ji,jj) ) THEN 1207 zvijk = -zdept(ji,jjs,jk1) 1208 EXIT 1209 ENDIF 1210 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 1211 zpnss = zpnss + & 1212 integ_spline(zdept(ji,jjs,jk1), zdeps, & 1213 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1214 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1215 jk1 = jk1 + 1 1216 END DO 1217 1218 ! integrate the pressure on the deep side 1219 jk1 = jk 1220 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1221 IF( jk1 == 1 ) THEN 1222 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1223 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1224 bsp(ji,jjd,1) , csp(ji,jjd,1), & 1225 dsp(ji,jjd,1) ) * zdeps 1226 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1227 EXIT 1228 ENDIF 1229 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 1230 zpnsd = zpnsd + & 1231 integ_spline(zdeps, zdept(ji,jjd,jk1), & 1232 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1233 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1234 jk1 = jk1 - 1 1235 END DO 1236 1237 ! update the momentum trends in v direction 1238 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1239 IF( .NOT.ln_linssh ) THEN 1240 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1241 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1242 ELSE 1243 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1244 ENDIF 1245 IF( ln_wd_il ) THEN 1246 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1247 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1248 ENDIF 1249 1250 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 1196 1251 ENDIF 1197 IF( ln_wd_il ) THEN1198 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj)1199 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj)1200 ENDIF1201 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk)1202 ENDIF1203 1204 !!!!! for v equation1205 IF( jk <= mbkv(ji,jj) ) THEN1206 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN1207 jjs = jj + 1; jjd = jj1208 ELSE1209 jjs = jj ; jjd = jj + 11210 ENDIF1211 1212 ! integrate the pressure on the shallow side1213 jk1 = jk1214 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk )1215 IF( jk1 == mbkv(ji,jj) ) THEN1216 zvijk = -zdept(ji,jjs,jk1)1217 EXIT1218 ENDIF1219 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk)1220 zpnss = zpnss + &1221 integ_spline(zdept(ji,jjs,jk1), zdeps, &1222 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), &1223 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) )1224 jk1 = jk1 + 11225 END DO1226 1227 ! integrate the pressure on the deep side1228 jk1 = jk1229 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk )1230 IF( jk1 == 1 ) THEN1231 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad)1232 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), &1233 bsp(ji,jjd,1), csp(ji,jjd,1), &1234 dsp(ji,jjd,1) ) * zdeps1235 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps1236 EXIT1237 ENDIF1238 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk)1239 zpnsd = zpnsd + &1240 integ_spline(zdeps, zdept(ji,jjd,jk1), &1241 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), &1242 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) )1243 jk1 = jk1 - 11244 END DO1245 1246 1247 ! update the momentum trends in v direction1248 1249 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) )1250 IF( .NOT.ln_linssh ) THEN1251 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * &1252 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) )1253 ELSE1254 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd )1255 ENDIF1256 IF( ln_wd_il ) THEN1257 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj)1258 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj)1259 ENDIF1260 1261 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk)1262 ENDIF1263 1252 ! 1264 1253 END_3D … … 1279 1268 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1280 1269 !!---------------------------------------------------------------------- 1281 REAL(wp), DIMENSION( :,:,:), INTENT(in ) :: fsp, xsp ! value and coordinate1282 REAL(wp), DIMENSION( :,:,:), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function1283 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear1270 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: fsp, xsp ! value and coordinate 1271 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1272 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1284 1273 ! 1285 1274 INTEGER :: ji, jj, jk ! dummy loop indices 1286 INTEGER :: jpi, jpj, jpkm11287 1275 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 1288 1276 REAL(wp) :: zdxtmp1, zdxtmp2, zalpha 1289 REAL(wp) :: zdf(size(fsp,3)) 1290 !!---------------------------------------------------------------------- 1291 ! 1292 !!gm WHAT !!!!! THIS IS VERY DANGEROUS !!!!! 1293 jpi = size(fsp,1) 1294 jpj = size(fsp,2) 1295 jpkm1 = MAX( 1, size(fsp,3) - 1 ) 1277 REAL(wp) :: zdf(jpk) 1278 !!---------------------------------------------------------------------- 1296 1279 ! 1297 1280 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 1298 DO ji = 1, jpi 1299 DO jj = 1, jpj 1300 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 1301 ! DO jk = 2, jpkm1-1 1302 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1303 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1304 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1305 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1306 ! 1307 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1308 ! 1309 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1310 ! zdf(jk) = 0._wp 1311 ! ELSE 1312 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1313 ! ENDIF 1314 ! END DO 1315 1316 !!Simply geometric average 1317 DO jk = 2, jpkm1-1 1318 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1319 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1320 1321 IF(zdf1 * zdf2 <= 0._wp) THEN 1322 zdf(jk) = 0._wp 1323 ELSE 1324 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1325 ENDIF 1326 END DO 1327 1328 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1329 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1330 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1331 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpkm1 - 1) 1332 1333 DO jk = 1, jpkm1 - 1 1334 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1335 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1336 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1337 zddf1 = -2._wp * ztmp1 + ztmp2 1338 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1339 zddf2 = 2._wp * ztmp1 - ztmp2 1340 1341 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1342 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1343 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1344 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1345 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1346 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1347 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1348 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1349 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1350 END DO 1281 DO_2D( 1, 1, 1, 1 ) 1282 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 1283 ! DO jk = 2, jpkm1-1 1284 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1285 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1286 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1287 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1288 ! 1289 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1290 ! 1291 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1292 ! zdf(jk) = 0._wp 1293 ! ELSE 1294 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1295 ! ENDIF 1296 ! END DO 1297 1298 !!Simply geometric average 1299 DO jk = 2, jpk-2 1300 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1301 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1302 1303 IF(zdf1 * zdf2 <= 0._wp) THEN 1304 zdf(jk) = 0._wp 1305 ELSE 1306 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1307 ENDIF 1351 1308 END DO 1352 END DO 1309 1310 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1311 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1312 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1313 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2) 1314 1315 DO jk = 1, jpk-2 1316 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1317 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1318 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1319 zddf1 = -2._wp * ztmp1 + ztmp2 1320 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1321 zddf2 = 2._wp * ztmp1 - ztmp2 1322 1323 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1324 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1325 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1326 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1327 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1328 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1329 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1330 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1331 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1332 END DO 1333 END_2D 1353 1334 1354 1335 ELSEIF ( polynomial_type == 2 ) THEN ! Linear 1355 DO ji = 1, jpi 1356 DO jj = 1, jpj 1357 DO jk = 1, jpkm1-1 1358 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1359 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1360 1361 dsp(ji,jj,jk) = 0._wp 1362 csp(ji,jj,jk) = 0._wp 1363 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1364 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1365 END DO 1366 END DO 1367 END DO 1336 DO_3D( 1, 1, 1, 1, 1, jpk-2 ) 1337 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1338 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1339 1340 dsp(ji,jj,jk) = 0._wp 1341 csp(ji,jj,jk) = 0._wp 1342 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1343 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1344 END_3D 1368 1345 ! 1369 1346 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
