- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynhpg.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/OCE/DYN/dynhpg.F90
r14433 r15548 189 189 & CALL ctl_stop( 'dyn_hpg_init : ln_hpg_isf=T requires ln_isfcav=T and vice versa' ) 190 190 ! 191 #if defined key_qco192 IF( ln_hpg_isf ) THEN193 CALL ctl_stop( 'dyn_hpg_init : key_qco and ln_hpg_isf not yet compatible' )194 ENDIF195 #endif196 191 ! 197 192 ! ! Set nhpg from ln_hpg_... flags & consistency check … … 266 261 INTEGER :: ji, jj, jk ! dummy loop indices 267 262 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 ' 263 REAL(wp), DIMENSION(A2D(nn_hls)) :: zhpi, zhpj 264 !!---------------------------------------------------------------------- 265 ! 266 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 267 IF( kt == nit000 ) THEN 268 IF(lwp) WRITE(numout,*) 269 IF(lwp) WRITE(numout,*) 'dyn:hpg_zco : hydrostatic pressure gradient trend' 270 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate case ' 271 ENDIF 275 272 ENDIF 276 273 ! … … 318 315 INTEGER :: iku, ikv ! temporary integers 319 316 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' 317 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 318 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zgtsu, zgtsv 319 REAL(wp), DIMENSION(A2D(nn_hls) ) :: zgru, zgrv 320 !!---------------------------------------------------------------------- 321 ! 322 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 323 IF( kt == nit000 ) THEN 324 IF(lwp) WRITE(numout,*) 325 IF(lwp) WRITE(numout,*) 'dyn:hpg_zps : hydrostatic pressure gradient trend' 326 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ z-coordinate with partial steps - vector optimization' 327 ENDIF 329 328 ENDIF 330 329 … … 410 409 REAL(wp) :: zcoef0, zuap, zvap, ztmp ! local scalars 411 410 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 412 REAL(wp), DIMENSION( jpi,jpj,jpk):: zhpi, zhpj411 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 413 412 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 414 413 !!---------------------------------------------------------------------- 415 414 ! 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' 415 IF( ln_wd_il ) ALLOCATE(zcpx(A2D(nn_hls)), zcpy(A2D(nn_hls))) 416 ! 417 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 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' 422 ENDIF 422 423 ENDIF 423 424 ! … … 462 463 END IF 463 464 END_2D 464 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )465 465 END IF 466 466 ! … … 548 548 REAL(wp) :: ze3w, ze3wi1, ze3wj1 ! local scalars 549 549 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_oce550 REAL(wp), DIMENSION(A2D(nn_hls),jpk ) :: zhpi, zhpj 551 REAL(wp), DIMENSION(A2D(nn_hls),jpts) :: zts_top 552 REAL(wp), DIMENSION(A2D(nn_hls)) :: zrhd_top, zdep_top 553 553 !!---------------------------------------------------------------------- 554 554 ! … … 560 560 ! compute rhd at the ice/oce interface (ocean side) 561 561 ! 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 569 CALL eos( zts_top, risfdep, zrhdtop_oce ) 562 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 563 ikt = mikt(ji,jj) 564 zts_top(ji,jj,1) = ts(ji,jj,ikt,1,Kmm) 565 zts_top(ji,jj,2) = ts(ji,jj,ikt,2,Kmm) 566 zdep_top(ji,jj) = MAX( risfdep(ji,jj) , gdept(ji,jj,1,Kmm) ) 567 END_2D 568 CALL eos( zts_top, zdep_top, zrhd_top ) 570 569 571 570 ! !===========================! … … 579 578 ! ! we assume ISF is in isostatic equilibrium 580 579 zhpi(ji,jj,1) = zcoef0 * r1_e1u(ji,jj) * ( risfload(ji+1,jj) - risfload(ji,jj) & 581 & + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhd top_oce(ji+1,jj) ) &582 & - ze3w * ( rhd(ji ,jj,ikt ) + zrhd top_oce(ji ,jj) ) ) )580 & + 0.5_wp * ( ze3wi1 * ( rhd(ji+1,jj,ikti1) + zrhd_top(ji+1,jj) ) & 581 & - ze3w * ( rhd(ji ,jj,ikt ) + zrhd_top(ji ,jj) ) ) ) 583 582 zhpj(ji,jj,1) = zcoef0 * r1_e2v(ji,jj) * ( risfload(ji,jj+1) - risfload(ji,jj) & 584 & + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhd top_oce(ji,jj+1) ) &585 & - ze3w * ( rhd(ji,jj ,ikt ) + zrhd top_oce(ji,jj ) ) ) )583 & + 0.5_wp * ( ze3wj1 * ( rhd(ji,jj+1,iktj1) + zrhd_top(ji,jj+1) ) & 584 & - ze3w * ( rhd(ji,jj ,ikt ) + zrhd_top(ji,jj ) ) ) ) 586 585 ! ! s-coordinate pressure gradient correction (=0 if z coordinate) 587 586 zuap = -zcoef0 * ( rhd (ji+1,jj,1) + rhd (ji,jj,1) ) & … … 636 635 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 637 636 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 638 REAL(wp) :: z_grav_10, z1_12 637 REAL(wp) :: z_grav_10, z1_12, z1_cff 639 638 REAL(wp) :: cffu, cffx ! " " 640 639 REAL(wp) :: cffv, cffy ! " " 641 640 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 arrays641 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zhpj 642 643 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 644 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 645 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 646 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 647 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 648 REAL(wp), DIMENSION(A2D(nn_hls)) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 650 649 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 651 650 !!---------------------------------------------------------------------- 652 651 ! 653 652 IF( ln_wd_il ) THEN 654 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )653 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 655 654 DO_2D( 0, 0, 0, 0 ) 656 655 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & … … 689 688 END IF 690 689 END_2D 691 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )692 690 END IF 693 691 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' 692 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 693 IF( kt == nit000 ) THEN 694 IF(lwp) WRITE(numout,*) 695 IF(lwp) WRITE(numout,*) 'dyn:hpg_djc : hydrostatic pressure gradient trend' 696 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, density Jacobian with cubic polynomial scheme' 697 ENDIF 698 698 ENDIF 699 699 … … 723 723 zdz_k (:,:,:) = 0._wp 724 724 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 725 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 726 cffw = MAX( 2._wp * zdrhoz(ji,jj,jk) * zdrhoz(ji,jj,jk+1), 0._wp ) 727 z1_cff = zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) 728 zdrho_k(ji,jj,jk) = cffw / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 730 729 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 731 730 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) … … 737 736 738 737 ! 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) 738 DO_2D( 1, 1, 1, 1 ) 739 zdrho_k(ji,jj,1) = aco_bc_vrt * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) - bco_bc_vrt * zdrho_k(ji,jj,2) 740 zdz_k (ji,jj,1) = aco_bc_vrt * (-gde3w(ji,jj,2) + gde3w(ji,jj,1) ) - bco_bc_vrt * zdz_k (ji,jj,2) 741 END_2D 741 742 742 743 DO_2D( 1, 1, 1, 1 ) … … 785 786 ! 5. compute and store elementary horizontal differences in provisional arrays 786 787 !---------------------------------------------------------------------------------------- 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. ) 788 zdrhox(:,:,:) = 0._wp 789 zdzx (:,:,:) = 0._wp 790 zdrhoy(:,:,:) = 0._wp 791 zdzy (:,:,:) = 0._wp 792 793 DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpkm1 ) 794 zdrhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji ,jj ,jk) 795 zdzx (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji+1,jj ,jk) 796 zdrhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji ,jj ,jk) 797 zdzy (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji ,jj+1,jk) 798 END_3D 799 800 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 801 797 802 !------------------------------------------------------------------------- … … 800 805 801 806 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 807 cffu = MAX( 2._wp * zdrhox(ji-1,jj,jk) * zdrhox(ji,jj,jk), 0._wp ) 808 z1_cff = zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) 809 zdrho_i(ji,jj,jk) = cffu / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 810 811 cffx = MAX( 2._wp * zdzx(ji-1,jj,jk) * zdzx(ji,jj,jk), 0._wp ) 812 z1_cff = zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) 813 zdz_i(ji,jj,jk) = cffx / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 814 815 cffv = MAX( 2._wp * zdrhoy(ji,jj-1,jk) * zdrhoy(ji,jj,jk), 0._wp ) 816 z1_cff = zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) 817 zdrho_j(ji,jj,jk) = cffv / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 818 819 cffy = MAX( 2._wp * zdzy(ji,jj-1,jk) * zdzy(ji,jj,jk), 0._wp ) 820 z1_cff = zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) 821 zdz_j(ji,jj,jk) = cffy / SIGN( MAX( ABS(z1_cff), zep ), z1_cff ) 829 822 END_3D 830 823 … … 840 833 zz_drho_j(:,:) = zdrho_j(:,:,jk) 841 834 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 835 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 836 DO_2D( 0, 0, 0, 1 ) 837 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 838 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) 839 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 840 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 IF841 END_2D 842 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 843 DO_2D( -1, 1, 0, 1 ) 844 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 845 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) 846 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 847 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 IF848 END_2D 849 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 850 DO_2D( 0, 1, 0, 0 ) 851 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 852 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) 853 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) 854 END IF 855 END_2D 856 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 857 DO_2D( 0, 1, -1, 1 ) 858 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 859 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) 860 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 861 END IF 871 862 END_2D … … 974 965 REAL(wp) :: zrhdt1 975 966 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, dsp967 REAL(wp), DIMENSION(A2D(nn_hls)) :: zpgu, zpgv ! 2D workspace 968 REAL(wp), DIMENSION(A2D(nn_hls)) :: zsshu_n, zsshv_n 969 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zdept, zrhh 970 REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 980 971 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 981 972 !!---------------------------------------------------------------------- 982 973 ! 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' 974 IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile 975 IF( kt == nit000 ) THEN 976 IF(lwp) WRITE(numout,*) 977 IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 978 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate case, cubic spline pressure Jacobian' 979 ENDIF 987 980 ENDIF 988 981 … … 1001 994 ! 1002 995 IF( ln_wd_il ) THEN 1003 ALLOCATE( zcpx( jpi,jpj) , zcpy(jpi,jpj) )996 ALLOCATE( zcpx(A2D(nn_hls)) , zcpy(A2D(nn_hls)) ) 1004 997 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 998 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 999 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1000 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) ) > & 1001 & rn_wdmin1 + rn_wdmin2 1002 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji+1,jj,Kmm) ) > 1.E-12 ) .AND. & 1003 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji+1,jj,Kmm) ) > & 1004 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1005 1006 IF(ll_tmp1) THEN 1007 zcpx(ji,jj) = 1.0_wp 1008 ELSE IF(ll_tmp2) THEN 1009 ! no worries about ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm) = 0, it won't happen ! here 1010 zcpx(ji,jj) = ABS( (ssh(ji+1,jj,Kmm) + ht_0(ji+1,jj) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1011 & / (ssh(ji+1,jj,Kmm) - ssh(ji ,jj,Kmm)) ) 1012 zcpx(ji,jj) = MAX(MIN( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1013 ELSE 1014 zcpx(ji,jj) = 0._wp 1015 END IF 1016 1017 ll_tmp1 = MIN( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1018 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1019 & MAX( ssh(ji,jj,Kmm) + ht_0(ji,jj), ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) ) > & 1020 & rn_wdmin1 + rn_wdmin2 1021 ll_tmp2 = ( ABS( ssh(ji,jj,Kmm) - ssh(ji,jj+1,Kmm) ) > 1.E-12 ) .AND. & 1022 & ( MAX( ssh(ji,jj,Kmm) , ssh(ji,jj+1,Kmm) ) > & 1023 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1024 1025 IF(ll_tmp1) THEN 1026 zcpy(ji,jj) = 1.0_wp 1027 ELSE IF(ll_tmp2) THEN 1028 ! no worries about ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm) = 0, it won't happen ! here 1029 zcpy(ji,jj) = ABS( (ssh(ji,jj+1,Kmm) + ht_0(ji,jj+1) - ssh(ji,jj,Kmm) - ht_0(ji,jj)) & 1030 & / (ssh(ji,jj+1,Kmm) - ssh(ji,jj ,Kmm)) ) 1031 zcpy(ji,jj) = MAX(MIN( zcpy(ji,jj) , 1.0_wp),0.0_wp) 1041 1032 ELSE 1042 1033 zcpy(ji,jj) = 0._wp 1043 1034 ENDIF 1044 1035 END_2D 1045 CALL lbc_lnk( 'dynhpg', zcpx, 'U', 1.0_wp, zcpy, 'V', 1.0_wp )1046 1036 ENDIF 1047 1037 1048 1038 ! Clean 3-D work arrays 1049 1039 zhpi(:,:,:) = 0._wp 1050 zrhh(:,:,:) = rhd( :,:,:)1040 zrhh(:,:,:) = rhd(A2D(nn_hls),:) 1051 1041 1052 1042 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 1053 1043 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 ENDIF1044 jk = mbkt(ji,jj) 1045 IF( jk <= 1 ) THEN ; zrhh(ji,jj, : ) = 0._wp 1046 ELSEIF( jk == 2 ) THEN ; zrhh(ji,jj,jk+1:jpk) = rhd(ji,jj,jk) 1047 ELSEIF( jk < jpkm1 ) THEN 1048 DO jkk = jk+1, jpk 1049 zrhh(ji,jj,jkk) = interp1(gde3w(ji,jj,jkk ), gde3w(ji,jj,jkk-1), & 1050 & gde3w(ji,jj,jkk-2), zrhh (ji,jj,jkk-1), zrhh(ji,jj,jkk-2)) 1051 END DO 1052 ENDIF 1063 1053 END_2D 1064 1054 … … 1082 1072 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 1083 1073 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) * zrhdt11074 zrhdt1 = zrhh(ji,jj,1) - interp3( zdept(ji,jj,1), asp(ji,jj,1), bsp(ji,jj,1), & 1075 & csp(ji,jj,1), dsp(ji,jj,1) ) * 0.25_wp * e3w(ji,jj,1,Kmm) 1076 1077 ! assuming linear profile across the top half surface layer 1078 zhpi(ji,jj,1) = 0.5_wp * e3w(ji,jj,1,Kmm) * zrhdt1 1089 1079 END_2D 1090 1080 1091 1081 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 1092 1082 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) )1083 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 1084 & integ_spline( zdept(ji,jj,jk-1), zdept(ji,jj,jk), & 1085 & asp (ji,jj,jk-1), bsp (ji,jj,jk-1), & 1086 & csp (ji,jj,jk-1), dsp (ji,jj,jk-1) ) 1097 1087 END_3D 1098 1088 … … 1107 1097 ! & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1108 1098 !!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 ) 1099 zsshu_n(ji,jj) = (e1e2u(ji,jj) * ssh(ji,jj,Kmm) + e1e2u(ji+1, jj) * ssh(ji+1,jj,Kmm)) * & 1100 & r1_e1e2u(ji,jj) * umask(ji,jj,1) * 0.5_wp 1101 zsshv_n(ji,jj) = (e1e2v(ji,jj) * ssh(ji,jj,Kmm) + e1e2v(ji+1, jj) * ssh(ji,jj+1,Kmm)) * & 1102 & r1_e1e2v(ji,jj) * vmask(ji,jj,1) * 0.5_wp 1103 END_2D 1116 1104 1117 1105 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) )1106 zu(ji,jj,1) = - ( e3u(ji,jj,1,Kmm) - zsshu_n(ji,jj) ) 1107 zv(ji,jj,1) = - ( e3v(ji,jj,1,Kmm) - zsshv_n(ji,jj) ) 1120 1108 END_2D 1121 1109 1122 1110 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)1111 zu(ji,jj,jk) = zu(ji,jj,jk-1) - e3u(ji,jj,jk,Kmm) 1112 zv(ji,jj,jk) = zv(ji,jj,jk-1) - e3v(ji,jj,jk,Kmm) 1125 1113 END_3D 1126 1114 1127 1115 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)1116 zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * e3u(ji,jj,jk,Kmm) 1117 zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * e3v(ji,jj,jk,Kmm) 1130 1118 END_3D 1131 1119 1132 1120 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) ) )1121 zu(ji,jj,jk) = MIN( zu(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1122 zu(ji,jj,jk) = MAX( zu(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji+1,jj,jk) ) ) 1123 zv(ji,jj,jk) = MIN( zv(ji,jj,jk) , MAX( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1124 zv(ji,jj,jk) = MAX( zv(ji,jj,jk) , MIN( -zdept(ji,jj,jk) , -zdept(ji,jj+1,jk) ) ) 1137 1125 END_3D 1138 1126 1139 1127 1140 1128 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 1129 zpwes = 0._wp; zpwed = 0._wp 1130 zpnss = 0._wp; zpnsd = 0._wp 1131 zuijk = zu(ji,jj,jk) 1132 zvijk = zv(ji,jj,jk) 1133 1134 !!!!! for u equation 1135 IF( jk <= mbku(ji,jj) ) THEN 1136 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 1137 jis = ji + 1; jid = ji 1138 ELSE 1139 jis = ji; jid = ji +1 1140 ENDIF 1141 1142 ! integrate the pressure on the shallow side 1143 jk1 = jk 1144 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 1145 IF( jk1 == mbku(ji,jj) ) THEN 1146 zuijk = -zdept(jis,jj,jk1) 1147 EXIT 1148 ENDIF 1149 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 1150 zpwes = zpwes + & 1151 integ_spline(zdept(jis,jj,jk1), zdeps, & 1152 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 1153 csp(jis,jj,jk1), dsp(jis,jj,jk1)) 1154 jk1 = jk1 + 1 1155 END DO 1156 1157 ! integrate the pressure on the deep side 1158 jk1 = jk 1159 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 1160 IF( jk1 == 1 ) THEN 1161 zdeps = zdept(jid,jj,1) + MIN(zuijk, ssh(jid,jj,Kmm)*znad) 1162 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 1163 bsp(jid,jj,1) , csp(jid,jj,1), & 1164 dsp(jid,jj,1)) * zdeps 1165 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 1166 EXIT 1167 ENDIF 1168 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 1169 zpwed = zpwed + & 1170 integ_spline(zdeps, zdept(jid,jj,jk1), & 1171 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 1172 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) 1173 jk1 = jk1 - 1 1174 END DO 1175 1176 ! update the momentum trends in u direction 1177 zdpdx1 = zcoef0 * r1_e1u(ji,jj) * ( zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk) ) 1178 IF( .NOT.ln_linssh ) THEN 1179 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * & 1180 & ( REAL(jis-jid, wp) * (zpwes + zpwed) + (ssh(ji+1,jj,Kmm)-ssh(ji,jj,Kmm)) ) 1181 ELSE 1182 zdpdx2 = zcoef0 * r1_e1u(ji,jj) * REAL(jis-jid, wp) * (zpwes + zpwed) 1183 ENDIF 1184 IF( ln_wd_il ) THEN 1185 zdpdx1 = zdpdx1 * zcpx(ji,jj) * wdrampu(ji,jj) 1186 zdpdx2 = zdpdx2 * zcpx(ji,jj) * wdrampu(ji,jj) 1187 ENDIF 1188 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + (zdpdx1 + zdpdx2 - zpgu(ji,jj)) * umask(ji,jj,jk) 1152 1189 ENDIF 1153 1190 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) 1191 !!!!! for v equation 1192 IF( jk <= mbkv(ji,jj) ) THEN 1193 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 1194 jjs = jj + 1; jjd = jj 1195 ELSE 1196 jjs = jj ; jjd = jj + 1 1197 ENDIF 1198 1199 ! integrate the pressure on the shallow side 1200 jk1 = jk 1201 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 1202 IF( jk1 == mbkv(ji,jj) ) THEN 1203 zvijk = -zdept(ji,jjs,jk1) 1204 EXIT 1205 ENDIF 1206 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 1207 zpnss = zpnss + & 1208 integ_spline(zdept(ji,jjs,jk1), zdeps, & 1209 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 1210 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) 1211 jk1 = jk1 + 1 1212 END DO 1213 1214 ! integrate the pressure on the deep side 1215 jk1 = jk 1216 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 1217 IF( jk1 == 1 ) THEN 1218 zdeps = zdept(ji,jjd,1) + MIN(zvijk, ssh(ji,jjd,Kmm)*znad) 1219 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 1220 bsp(ji,jjd,1) , csp(ji,jjd,1), & 1221 dsp(ji,jjd,1) ) * zdeps 1222 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 1223 EXIT 1224 ENDIF 1225 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 1226 zpnsd = zpnsd + & 1227 integ_spline(zdeps, zdept(ji,jjd,jk1), & 1228 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 1229 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 1230 jk1 = jk1 - 1 1231 END DO 1232 1233 ! update the momentum trends in v direction 1234 zdpdy1 = zcoef0 * r1_e2v(ji,jj) * ( zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk) ) 1235 IF( .NOT.ln_linssh ) THEN 1236 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * & 1237 ( REAL(jjs-jjd, wp) * (zpnss + zpnsd) + (ssh(ji,jj+1,Kmm)-ssh(ji,jj,Kmm)) ) 1238 ELSE 1239 zdpdy2 = zcoef0 * r1_e2v(ji,jj) * REAL(jjs-jjd, wp) * (zpnss + zpnsd ) 1240 ENDIF 1241 IF( ln_wd_il ) THEN 1242 zdpdy1 = zdpdy1 * zcpy(ji,jj) * wdrampv(ji,jj) 1243 zdpdy2 = zdpdy2 * zcpy(ji,jj) * wdrampv(ji,jj) 1244 ENDIF 1245 1246 pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) + (zdpdy1 + zdpdy2 - zpgv(ji,jj)) * vmask(ji,jj,jk) 1196 1247 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 1248 ! 1264 1249 END_3D … … 1279 1264 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 1280 1265 !!---------------------------------------------------------------------- 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: Linear1266 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(in ) :: fsp, xsp ! value and coordinate 1267 REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT( out) :: asp, bsp, csp, dsp ! coefficients of the interpoated function 1268 INTEGER , INTENT(in ) :: polynomial_type ! 1: cubic spline ; 2: Linear 1284 1269 ! 1285 1270 INTEGER :: ji, jj, jk ! dummy loop indices 1286 INTEGER :: jpi, jpj, jpkm11287 1271 REAL(wp) :: zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 1288 1272 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 ) 1273 REAL(wp) :: zdf(jpk) 1274 !!---------------------------------------------------------------------- 1296 1275 ! 1297 1276 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 1277 DO_2D( 1, 1, 1, 1 ) 1278 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 1279 ! DO jk = 2, jpkm1-1 1280 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1281 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1282 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1283 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1284 ! 1285 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1286 ! 1287 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1288 ! zdf(jk) = 0._wp 1289 ! ELSE 1290 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1291 ! ENDIF 1292 ! END DO 1293 1294 !!Simply geometric average 1295 DO jk = 2, jpk-2 1296 zdf1 = (fsp(ji,jj,jk ) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk ) - xsp(ji,jj,jk-1)) 1297 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk )) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk )) 1298 1299 IF(zdf1 * zdf2 <= 0._wp) THEN 1300 zdf(jk) = 0._wp 1301 ELSE 1302 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1303 ENDIF 1351 1304 END DO 1352 END DO 1305 1306 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1307 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1308 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1309 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - 0.5_wp * zdf(jpk - 2) 1310 1311 DO jk = 1, jpk-2 1312 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1313 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1314 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1315 zddf1 = -2._wp * ztmp1 + ztmp2 1316 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1317 zddf2 = 2._wp * ztmp1 - ztmp2 1318 1319 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1320 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1321 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1322 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1323 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1324 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1325 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1326 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1327 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1328 END DO 1329 END_2D 1353 1330 1354 1331 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 1332 DO_3D( 1, 1, 1, 1, 1, jpk-2 ) 1333 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1334 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1335 1336 dsp(ji,jj,jk) = 0._wp 1337 csp(ji,jj,jk) = 0._wp 1338 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1339 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1340 END_3D 1368 1341 ! 1369 1342 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
