Changeset 15035
- Timestamp:
- 2021-06-21T13:07:02+02:00 (3 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0-HEAD_DJC_backwards_port/src/OCE/DYN/dynhpg.F90
r13095 r15035 71 71 ! 72 72 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 73 LOGICAL, PUBLIC :: ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 74 REAL(wp), PUBLIC :: aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt 73 75 74 76 !! * Substitutions … … 147 149 !! 148 150 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 149 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 151 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, & 152 & ln_hpg_djc_vN_hor, ln_hpg_djc_vN_vrt 150 153 !!---------------------------------------------------------------------- 151 154 ! … … 172 175 ENDIF 173 176 ! 174 IF( ln_hpg_djc ) & 175 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method', & 176 & ' currently disabled (bugs under investigation).' , & 177 & ' Please select either ln_hpg_sco or ln_hpg_prj instead' ) 178 ! 179 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 177 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) ) & 180 178 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 181 179 & ' the standard jacobian formulation hpg_sco or ' , & … … 213 211 ENDIF 214 212 ! 213 IF (ln_hpg_djc_vN_hor) THEN ! Von Neumann boundary condition 214 aco_bc_hor = 6.0_wp/5.0_wp 215 bco_bc_hor = 7.0_wp/15.0_wp 216 ELSE ! Linear extrapolation 217 aco_bc_hor = 3.0_wp/2.0_wp 218 bco_bc_hor = 1.0_wp/2.0_wp 219 END IF 220 IF (ln_hpg_djc_vN_vrt) THEN ! Von Neumann boundary condition 221 aco_bc_vrt = 6.0_wp/5.0_wp 222 bco_bc_vrt = 7.0_wp/15.0_wp 223 ELSE ! Linear extrapolation 224 aco_bc_vrt = 3.0_wp/2.0_wp 225 bco_bc_vrt = 1.0_wp/2.0_wp 226 END IF 215 227 IF ( .NOT. ln_isfcav ) THEN !--- no ice shelf load 216 228 riceload(:,:) = 0._wp … … 673 685 !! 674 686 INTEGER :: ji, jj, jk ! dummy loop indices 687 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 675 688 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 676 REAL(wp) :: z1_10, cffu, cffx ! " " 677 REAL(wp) :: z1_12, cffv, cffy ! " " 689 REAL(wp) :: z_grav_10, z1_12 690 REAL(wp) :: cffu, cffx ! " " 691 REAL(wp) :: cffv, cffy ! " " 678 692 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 679 693 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 680 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dzx, dzy, dzz, dzu, dzv, dzw 681 REAL(wp), DIMENSION(jpi,jpj,jpk) :: drhox, drhoy, drhoz, drhou, drhov, drhow 682 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rho_i, rho_j, rho_k 694 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rhd_opt 695 696 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 697 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 698 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 699 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 700 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 701 REAL(wp), DIMENSION(jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 683 702 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 684 703 !!---------------------------------------------------------------------- … … 735 754 ! Local constant initialization 736 755 zcoef0 = - grav * 0.5_wp 737 z1_10 = 1._wp / 10._wp 738 z1_12 = 1._wp / 12._wp 739 756 z_grav_10 = grav / 10._wp 757 z1_12 = 1.0_wp / 12._wp 758 759 IF( .NOT. ln_linssh ) THEN 760 rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option 761 ELSE 762 rhd_opt(:,:,:) = rhd(:,:,:) 763 END IF 740 764 !---------------------------------------------------------------------------------------- 741 ! compute and store in provisional arrays elementary vertical and horizontal differences765 ! 1. compute and store elementary vertical differences in provisional arrays 742 766 !---------------------------------------------------------------------------------------- 743 767 744 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 768 !!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really 769 !! zdzz, zdzx and zdzy changed to heights rather than depths; lower bounds of jj and ji changed from 2 to 1 745 770 746 771 DO jk = 2, jpkm1 772 DO jj = 1, jpj 773 DO ji = 1, jpi 774 zdrhoz(ji,jj,jk) = rhd_opt (ji ,jj ,jk) - rhd_opt (ji,jj,jk-1) 775 zdzz (ji,jj,jk) = - gde3w_n(ji ,jj ,jk) + gde3w_n(ji,jj,jk-1) 776 END DO 777 END DO 778 END DO 779 780 !------------------------------------------------------------------------- 781 ! 2. compute harmonic averages for vertical differences using eq. 5.18 782 !------------------------------------------------------------------------- 783 zep = 1.e-15 784 785 !!bug gm zdrhoz not defined at level 1 and used (jk-1 with jk=2) ; this issue has now been addressed 786 !!bug gm idem for zdrhox, zdrhoy et ji=jpi and jj=jpj; idem 787 788 !! zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 789 !! DO loops broken up so that zdrho_k and zdz_k are calculated only for jk = 2 to jpk - 2 790 zdrho_k(:,:,:) = 0._wp 791 zdz_k (:,:,:) = 0._wp 792 793 DO jk = 2, jpk - 2 794 DO jj = 1, jpj 795 DO ji = 1, jpi 796 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 797 IF( cffw > zep) THEN 798 zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 799 ENDIF 800 801 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 802 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 803 804 END DO 805 END DO 806 END DO 807 808 !---------------------------------------------------------------------------------- 809 ! 3. apply boundary conditions at top and bottom using 5.36-5.37 810 !---------------------------------------------------------------------------------- 811 812 ! for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 813 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt (:,:,2) - rhd_opt (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 814 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w_n(:,:,2) + gde3w_n(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 815 816 DO jj = 1, jpj 817 DO ji = 1, jpi 818 IF ( mbkt(ji,jj)>1 ) THEN 819 iktb = mbkt(ji,jj) 820 zdrho_k(ji,jj,iktb) = aco_bc_vrt * ( rhd_opt(ji,jj,iktb) - rhd_opt(ji,jj,iktb-1) ) - bco_bc_vrt * zdrho_k(ji,jj,iktb-1) 821 zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w_n(ji,jj,iktb) + gde3w_n(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) 822 END IF 823 END DO 824 END DO 825 826 !-------------------------------------------------------------- 827 ! 4. Compute side face integrals 828 !------------------------------------------------------------- 829 830 !! sshn replaces e3w_n ; gde3w_n is a depth; the formulae involve heights 831 !! rho_k stores grav * FX / rho_0 832 833 !-------------------------------------------------------------- 834 ! 4. a) Upper half of top-most grid box, compute and store 835 !------------------------------------------------------------- 836 ! Concerns that zdrho_k might be oddly defined (just -1.5rho) for single celled columns are resolved by the fact that z_rho_k is defined explicity for the surface layer 837 DO jj = 2, jpj 838 DO ji = 2, jpi 839 z_rho_k(ji,jj,1) = grav * ( sshn(ji,jj) + gde3w_n(ji,jj,1) ) & ! *** AY sshn included in djc but not in sco 840 & * ( rhd_opt(ji,jj,1) & 841 & + 0.5_wp * ( rhd_opt (ji,jj,2) - rhd_opt (ji,jj,1) ) & 842 & * ( sshn (ji,jj ) + gde3w_n(ji,jj,1) ) & 843 & / ( - gde3w_n(ji,jj,2) + gde3w_n(ji,jj,1) ) ) 844 END DO 845 END DO 846 847 !-------------------------------------------------------------- 848 ! 4. b) Interior faces, compute and store 849 !------------------------------------------------------------- 850 851 DO jk = 2, jpkm1 852 DO jj = 2, jpj 853 DO ji = 2, jpi 854 z_rho_k(ji,jj,jk) = zcoef0 * ( rhd_opt (ji,jj,jk) + rhd_opt (ji,jj,jk-1) ) & 855 & * ( - gde3w_n(ji,jj,jk) + gde3w_n(ji,jj,jk-1) ) & 856 & + z_grav_10 * ( & 857 & ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) & 858 & * ( - gde3w_n(ji,jj,jk) + gde3w_n(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) & 859 & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) & 860 & * ( rhd_opt (ji,jj,jk) - rhd_opt (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) & 861 & ) 862 END DO 863 END DO 864 END DO 865 866 !---------------------------------------------------------------------------------------- 867 ! 5. compute and store elementary horizontal differences in provisional arrays 868 !---------------------------------------------------------------------------------------- 869 870 DO jk = 1, jpkm1 871 DO jj = 1, jpjm1 872 DO ji = 1, fs_jpim1 ! vector opt. 873 zdrhox(ji,jj,jk) = rhd_opt (ji+1,jj ,jk) - rhd_opt (ji,jj,jk ) 874 zdzx (ji,jj,jk) = - gde3w_n(ji+1,jj ,jk) + gde3w_n(ji,jj,jk ) 875 zdrhoy(ji,jj,jk) = rhd_opt (ji ,jj+1,jk) - rhd_opt (ji,jj,jk ) 876 zdzy (ji,jj,jk) = - gde3w_n(ji ,jj+1,jk) + gde3w_n(ji,jj,jk ) 877 END DO 878 END DO 879 END DO 880 881 CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 882 883 !------------------------------------------------------------------------- 884 ! 6. compute harmonic averages using eq. 5.18 885 !------------------------------------------------------------------------- 886 887 DO jk = 1, jpkm1 888 DO jj = 2, jpj 889 DO ji = fs_2, jpi ! vector opt. 890 891 cffu = 2._wp * zdrhox(ji-1,jj ,jk) * zdrhox(ji,jj,jk ) 892 IF( cffu > zep ) THEN 893 zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 894 ELSE 895 zdrho_i(ji,jj,jk ) = 0._wp 896 ENDIF 897 898 cffx = 2._wp * zdzx (ji-1,jj ,jk) * zdzx (ji,jj,jk ) 899 IF( cffx > zep ) THEN 900 zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 901 ELSE 902 zdz_i(ji,jj,jk) = 0._wp 903 ENDIF 904 905 cffv = 2._wp * zdrhoy(ji ,jj-1,jk) * zdrhoy(ji,jj,jk ) 906 IF( cffv > zep ) THEN 907 zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 908 ELSE 909 zdrho_j(ji,jj,jk) = 0._wp 910 ENDIF 911 912 cffy = 2._wp * zdzy (ji ,jj-1,jk) * zdzy (ji,jj,jk ) 913 IF( cffy > zep ) THEN 914 zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 915 ELSE 916 zdz_j(ji,jj,jk) = 0._wp 917 ENDIF 918 919 END DO 920 END DO 921 END DO 922 923 !!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point 924 925 !---------------------------------------------------------------------------------- 926 ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 927 !---------------------------------------------------------------------------------- 928 929 DO jk = 1, jpkm1 930 DO jj = 2, jpj 931 DO ji = 2, jpi ! vector opt. 932 zz_drho_i(ji,jj) = zdrho_i(ji,jj,jk) 933 zz_dz_i (ji,jj) = zdz_i (ji,jj,jk) 934 zz_drho_j(ji,jj) = zdrho_j(ji,jj,jk) 935 zz_dz_j (ji,jj) = zdz_j (ji,jj,jk) 936 ! Walls coming from left should check from 2 to jpi-1 (and jpj=2-jpj) 937 IF (ji < jpi) THEN 938 IF ( umask(ji,jj,jk) > 0.5_wp .AND. tmask(ji-1,jj,jk) < 0.5_wp .AND. umask(ji+1,jj,jk) > 0.5_wp) THEN 939 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt (ji+1,jj,jk) - rhd_opt (ji,jj,jk) ) - bco_bc_hor * zdrho_i(ji+1,jj,jk) 940 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w_n(ji+1,jj,jk) + gde3w_n(ji,jj,jk) ) - bco_bc_hor * zdz_i (ji+1,jj,jk) 941 END IF 942 END IF 943 ! Walls coming from right should check from 3 to jpi (and jpj=2-jpj) 944 IF (ji > 2) THEN 945 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 946 zz_drho_i(ji,jj) = aco_bc_hor * ( rhd_opt (ji,jj,jk) - rhd_opt (ji-1,jj,jk) ) - bco_bc_hor * zdrho_i(ji-1,jj,jk) 947 zz_dz_i (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(ji-1,jj,jk) ) - bco_bc_hor * zdz_i (ji-1,jj,jk) 948 END IF 949 END IF 950 ! Walls coming from left should check from 2 to jpj-1 (and jpi=2-jpi) 951 IF (jj < jpj) THEN 952 IF ( vmask(ji,jj,jk) > 0.5_wp .AND. tmask(ji,jj-1,jk) < 0.5_wp .AND. vmask(ji,jj+1,jk) > 0.5_wp) THEN 953 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt (ji,jj+1,jk) - rhd_opt (ji,jj,jk) ) - bco_bc_hor * zdrho_j(ji,jj+1,jk) 954 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj+1,jk) + gde3w_n(ji,jj,jk) ) - bco_bc_hor * zdz_j (ji,jj+1,jk) 955 END IF 956 END IF 957 ! Walls coming from right should check from 3 to jpj (and jpi=2-jpi) 958 IF (jj > 2) THEN 959 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 960 zz_drho_j(ji,jj) = aco_bc_hor * ( rhd_opt (ji,jj,jk) - rhd_opt (ji,jj-1,jk) ) - bco_bc_hor * zdrho_j(ji,jj-1,jk) 961 zz_dz_j (ji,jj) = aco_bc_hor * (-gde3w_n(ji,jj,jk) + gde3w_n(ji,jj-1,jk) ) - bco_bc_hor * zdz_j (ji,jj-1,jk) 962 END IF 963 END IF 964 965 END DO 966 END DO 967 968 DO jj = 2, jpj 969 DO ji = 2, jpi ! vector opt. 970 zdrho_i(ji,jj,jk) = zz_drho_i(ji,jj) 971 zdz_i (ji,jj,jk) = zz_dz_i (ji,jj) 972 zdrho_j(ji,jj,jk) = zz_drho_j(ji,jj) 973 zdz_j (ji,jj,jk) = zz_dz_j (ji,jj) 974 END DO 975 END DO 976 977 END DO 978 979 !-------------------------------------------------------------- 980 ! 7. Calculate integrals on side faces 981 !------------------------------------------------------------- 982 983 DO jk = 1, jpkm1 747 984 DO jj = 2, jpjm1 748 985 DO ji = fs_2, fs_jpim1 ! vector opt. 749 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 750 dzz (ji,jj,jk) = gde3w_n(ji ,jj ,jk) - gde3w_n(ji,jj,jk-1) 751 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 752 dzx (ji,jj,jk) = gde3w_n(ji+1,jj ,jk) - gde3w_n(ji,jj,jk ) 753 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 754 dzy (ji,jj,jk) = gde3w_n(ji ,jj+1,jk) - gde3w_n(ji,jj,jk ) 755 END DO 756 END DO 757 END DO 758 759 !------------------------------------------------------------------------- 760 ! compute harmonic averages using eq. 5.18 761 !------------------------------------------------------------------------- 762 zep = 1.e-15 763 764 !!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2) 765 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 766 767 DO jk = 2, jpkm1 768 DO jj = 2, jpjm1 769 DO ji = fs_2, fs_jpim1 ! vector opt. 770 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 771 772 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 773 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 774 775 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 776 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 777 778 IF( cffw > zep) THEN 779 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 780 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 781 ELSE 782 drhow(ji,jj,jk) = 0._wp 783 ENDIF 784 785 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 786 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 787 788 IF( cffu > zep ) THEN 789 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 790 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 791 ELSE 792 drhou(ji,jj,jk ) = 0._wp 793 ENDIF 794 795 IF( cffx > zep ) THEN 796 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) & 797 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) ) 798 ELSE 799 dzu(ji,jj,jk) = 0._wp 800 ENDIF 801 802 IF( cffv > zep ) THEN 803 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) & 804 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) ) 805 ELSE 806 drhov(ji,jj,jk) = 0._wp 807 ENDIF 808 809 IF( cffy > zep ) THEN 810 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 811 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 812 ELSE 813 dzv(ji,jj,jk) = 0._wp 814 ENDIF 815 816 END DO 817 END DO 818 END DO 819 820 !---------------------------------------------------------------------------------- 821 ! apply boundary conditions at top and bottom using 5.36-5.37 822 !---------------------------------------------------------------------------------- 823 drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 ) 824 drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 ) 825 drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 ) 826 827 drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 828 drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 829 drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 830 986 987 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 988 z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt (ji+1,jj,jk) + rhd_opt (ji,jj,jk) ) & 989 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) 990 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 991 z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( & 992 & ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) & 993 & * ( - gde3w_n(ji+1,jj,jk) + gde3w_n(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) & 994 & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) & 995 & * ( rhd_opt (ji+1,jj,jk) - rhd_opt (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) & 996 & ) 997 END IF 998 999 z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt (ji,jj+1,jk) + rhd_opt (ji,jj,jk) ) & 1000 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) 1001 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 1002 z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( & 1003 & ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) & 1004 & * ( - gde3w_n(ji,jj+1,jk) + gde3w_n(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) & 1005 & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) & 1006 & * ( rhd_opt (ji,jj+1,jk) - rhd_opt (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) & 1007 & ) 1008 END IF 1009 1010 END DO 1011 END DO 1012 END DO 831 1013 832 1014 !-------------------------------------------------------------- 833 ! Upper half of top-most grid box, compute and store1015 ! 8. Integrate in the vertical 834 1016 !------------------------------------------------------------- 835 836 !!bug gm : e3w-gde3w = 0.5*e3w .... and gde3w(2)-gde3w(1)=e3w(2) .... to be verified837 ! true if gde3w is really defined as the sum of the e3w scale factors as, it seems to me, it should be838 839 DO jj = 2, jpjm1840 DO ji = fs_2, fs_jpim1 ! vector opt.841 rho_k(ji,jj,1) = -grav * ( e3w_n(ji,jj,1) - gde3w_n(ji,jj,1) ) &842 & * ( rhd(ji,jj,1) &843 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) &844 & * ( e3w_n (ji,jj,1) - gde3w_n(ji,jj,1) ) &845 & / ( gde3w_n(ji,jj,2) - gde3w_n(ji,jj,1) ) )846 END DO847 END DO848 849 !!bug gm : here also, simplification is possible850 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop851 852 DO jk = 2, jpkm1853 DO jj = 2, jpjm1854 DO ji = fs_2, fs_jpim1 ! vector opt.855 856 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) &857 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) ) &858 & - grav * z1_10 * ( &859 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) &860 & * ( gde3w_n(ji,jj,jk) - gde3w_n(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) &861 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) &862 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) &863 & )864 865 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) &866 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) ) &867 & - grav* z1_10 * ( &868 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) &869 & * ( gde3w_n(ji+1,jj,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) &870 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) &871 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) &872 & )873 874 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) &875 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) ) &876 & - grav* z1_10 * ( &877 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) &878 & * ( gde3w_n(ji,jj+1,jk) - gde3w_n(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) &879 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) &880 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) &881 & )882 883 END DO884 END DO885 END DO886 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1., rho_i, 'U', 1., rho_j, 'V', 1. )887 1017 888 1018 ! --------------- … … 891 1021 DO jj = 2, jpjm1 892 1022 DO ji = fs_2, fs_jpim1 ! vector opt. 893 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) -rho_i(ji,jj,1) ) * r1_e1u(ji,jj)894 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) -rho_j(ji,jj,1) ) * r1_e2v(ji,jj)1023 zhpi(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji+1,jj ,1) - z_rho_i(ji,jj,1) ) * r1_e1u(ji,jj) 1024 zhpj(ji,jj,1) = ( z_rho_k(ji,jj,1) - z_rho_k(ji ,jj+1,1) - z_rho_j(ji,jj,1) ) * r1_e2v(ji,jj) 895 1025 IF( ln_wd_il ) THEN 896 1026 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 911 1041 ! hydrostatic pressure gradient along s-surfaces 912 1042 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 913 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) ) &914 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj)1043 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk ) ) & 1044 & - ( z_rho_i(ji,jj,jk) - z_rho_i(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 915 1045 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 916 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) ) &917 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj)1046 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk ) ) & 1047 & -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 918 1048 IF( ln_wd_il ) THEN 919 1049 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
Note: See TracChangeset
for help on using the changeset viewer.
