Changeset 14060
- Timestamp:
- 2020-12-03T17:03:42+01:00 (3 years ago)
- Location:
- NEMO/trunk
- Files:
-
- 2 edited
-
cfgs/SHARED/namelist_ref (modified) (1 diff)
-
src/OCE/DYN/dynhpg.F90 (modified) (9 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/cfgs/SHARED/namelist_ref
r14056 r14060 1013 1013 ln_hpg_isf = .false. ! s-coordinate (sco ) adapted to isf 1014 1014 ln_hpg_djc = .false. ! s-coordinate (Density Jacobian with Cubic polynomial) 1015 ln_hpg_djc_vnh = .true. ! hor. bc type for djc scheme (T=von Neumann, F=linear extrapolation) 1016 ln_hpg_djc_vnv = .true. ! vert. bc type for djc scheme (T=von Neumann, F=linear extrapolation) 1015 1017 ln_hpg_prj = .false. ! s-coordinate (Pressure Jacobian scheme) 1016 1018 / -
NEMO/trunk/src/OCE/DYN/dynhpg.F90
r13982 r14060 17 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 18 !! 3.6 ! 2014-11 (P. Mathiot) hpg_isf: original code for ice shelf cavity 19 !! 4.2 ! 2020-12 (M. Bell, A. Young) hpg_djc: revised djc scheme 19 20 !!---------------------------------------------------------------------- 20 21 … … 72 73 INTEGER, PARAMETER :: np_isf = 5 ! s-coordinate similar to sco modify for isf 73 74 ! 74 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 75 INTEGER, PUBLIC :: nhpg !: type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) (PUBLIC for TAM) 76 ! 77 LOGICAL :: ln_hpg_djc_vnh, ln_hpg_djc_vnv ! flag to specify hpg_djc boundary condition type 78 REAL(wp), PUBLIC :: aco_bc_hor, bco_bc_hor, aco_bc_vrt, bco_bc_vrt !: coefficients for hpg_djc hor and vert boundary conditions 75 79 76 80 !! * Substitutions … … 156 160 !! 157 161 NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, & 158 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf 162 & ln_hpg_djc, ln_hpg_prj, ln_hpg_isf, & 163 & ln_hpg_djc_vnh, ln_hpg_djc_vnv 159 164 !!---------------------------------------------------------------------- 160 165 ! … … 179 184 ENDIF 180 185 ! 181 IF( ln_hpg_djc ) & 182 & CALL ctl_stop('dyn_hpg_init : Density Jacobian: Cubic polynominal method', & 183 & ' currently disabled (bugs under investigation).' , & 184 & ' Please select either ln_hpg_sco or ln_hpg_prj instead' ) 185 ! 186 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf) ) & 186 IF( .NOT.ln_linssh .AND. .NOT.(ln_hpg_sco.OR.ln_hpg_prj.OR.ln_hpg_isf.OR.ln_hpg_djc) ) & 187 187 & CALL ctl_stop('dyn_hpg_init : non-linear free surface requires either ', & 188 188 & ' the standard jacobian formulation hpg_sco or ' , & … … 220 220 ENDIF 221 221 ! 222 IF ( ln_hpg_djc ) THEN 223 IF (ln_hpg_djc_vnh) THEN ! Von Neumann boundary condition 224 IF(lwp) WRITE(numout,*) ' horizontal bc: von Neumann ' 225 aco_bc_hor = 6.0_wp/5.0_wp 226 bco_bc_hor = 7.0_wp/15.0_wp 227 ELSE ! Linear extrapolation 228 IF(lwp) WRITE(numout,*) ' horizontal bc: linear extrapolation' 229 aco_bc_hor = 3.0_wp/2.0_wp 230 bco_bc_hor = 1.0_wp/2.0_wp 231 END IF 232 IF (ln_hpg_djc_vnv) THEN ! Von Neumann boundary condition 233 IF(lwp) WRITE(numout,*) ' vertical bc: von Neumann ' 234 aco_bc_vrt = 6.0_wp/5.0_wp 235 bco_bc_vrt = 7.0_wp/15.0_wp 236 ELSE ! Linear extrapolation 237 IF(lwp) WRITE(numout,*) ' vertical bc: linear extrapolation' 238 aco_bc_vrt = 3.0_wp/2.0_wp 239 bco_bc_vrt = 1.0_wp/2.0_wp 240 END IF 241 END IF 222 242 END SUBROUTINE dyn_hpg_init 223 243 … … 631 651 !! 632 652 INTEGER :: ji, jj, jk ! dummy loop indices 653 INTEGER :: iktb, iktt ! jk indices at tracer points for top and bottom points 633 654 REAL(wp) :: zcoef0, zep, cffw ! temporary scalars 634 REAL(wp) :: z1_10, cffu, cffx ! " " 635 REAL(wp) :: z1_12, cffv, cffy ! " " 655 REAL(wp) :: z_grav_10, z1_12 656 REAL(wp) :: cffu, cffx ! " " 657 REAL(wp) :: cffv, cffy ! " " 636 658 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables 637 659 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zhpi, zhpj 638 REAL(wp), DIMENSION(jpi,jpj,jpk) :: dzx, dzy, dzz, dzu, dzv, dzw 639 REAL(wp), DIMENSION(jpi,jpj,jpk) :: drhox, drhoy, drhoz, drhou, drhov, drhow 640 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rho_i, rho_j, rho_k 660 REAL(wp), DIMENSION(jpi,jpj,jpk) :: rhd_opt 661 662 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdzx, zdzy, zdzz ! Primitive grid differences ('delta_xyz') 663 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdz_i, zdz_j, zdz_k ! Harmonic average of primitive grid differences ('d_xyz') 664 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrhox, zdrhoy, zdrhoz ! Primitive rho differences ('delta_rho') 665 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdrho_i, zdrho_j, zdrho_k ! Harmonic average of primitive rho differences ('d_rho') 666 REAL(wp), DIMENSION(jpi,jpj,jpk) :: z_rho_i, z_rho_j, z_rho_k ! Face intergrals 667 REAL(wp), DIMENSION(jpi,jpj) :: zz_dz_i, zz_dz_j, zz_drho_i, zz_drho_j ! temporary arrays 641 668 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zcpx, zcpy !W/D pressure filter 642 669 !!---------------------------------------------------------------------- … … 691 718 ! Local constant initialization 692 719 zcoef0 = - grav * 0.5_wp 693 z1_10 = 1._wp / 10._wp 694 z1_12 = 1._wp / 12._wp 720 z_grav_10 = grav / 10._wp 721 z1_12 = 1.0_wp / 12._wp 722 723 724 !! To be removed once combined with KERNEL08. KERNEL08 action removes znad from the hpg module, spg is calculated in spg moduel instead. The final code will just have rhd(:,:,:), no rhd_opt(:,:,:) 725 IF( .NOT. ln_linssh ) THEN 726 rhd_opt(:,:,:) = rhd(:,:,:) + 1.0_wp ! for vvl option 727 ELSE 728 rhd_opt(:,:,:) = rhd(:,:,:) 729 END IF 695 730 696 731 !---------------------------------------------------------------------------------------- 697 ! compute and store in provisional arrays elementary vertical and horizontal differences732 ! 1. compute and store elementary vertical differences in provisional arrays 698 733 !---------------------------------------------------------------------------------------- 699 734 700 !!bug gm Not a true bug, but... dzz=e3w for dzx, dzy verify what it is really 701 702 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 703 drhoz(ji,jj,jk) = rhd (ji ,jj ,jk) - rhd (ji,jj,jk-1) 704 dzz (ji,jj,jk) = gde3w(ji ,jj ,jk) - gde3w(ji,jj,jk-1) 705 drhox(ji,jj,jk) = rhd (ji+1,jj ,jk) - rhd (ji,jj,jk ) 706 dzx (ji,jj,jk) = gde3w(ji+1,jj ,jk) - gde3w(ji,jj,jk ) 707 drhoy(ji,jj,jk) = rhd (ji ,jj+1,jk) - rhd (ji,jj,jk ) 708 dzy (ji,jj,jk) = gde3w(ji ,jj+1,jk) - gde3w(ji,jj,jk ) 735 !!bug gm Not a true bug, but... zdzz=e3w for zdzx, zdzy verify what it is really 736 737 DO_3D( 1, 1, 1, 1, 2, jpkm1 ) 738 zdrhoz(ji,jj,jk) = rhd_opt (ji ,jj ,jk) - rhd_opt (ji,jj,jk-1) 739 zdzz (ji,jj,jk) = - gde3w(ji ,jj ,jk) + gde3w(ji,jj,jk-1) 709 740 END_3D 710 741 711 742 !------------------------------------------------------------------------- 712 ! compute harmonic averages using eq. 5.18743 ! 2. compute harmonic averages for vertical differences using eq. 5.18 713 744 !------------------------------------------------------------------------- 714 745 zep = 1.e-15 715 746 716 !!bug gm drhoz not defined at level 1 and used (jk-1 with jk=2) 717 !!bug gm idem for drhox, drhoy et ji=jpi and jj=jpj 718 719 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 720 cffw = 2._wp * drhoz(ji ,jj ,jk) * drhoz(ji,jj,jk-1) 721 722 cffu = 2._wp * drhox(ji+1,jj ,jk) * drhox(ji,jj,jk ) 723 cffx = 2._wp * dzx (ji+1,jj ,jk) * dzx (ji,jj,jk ) 724 725 cffv = 2._wp * drhoy(ji ,jj+1,jk) * drhoy(ji,jj,jk ) 726 cffy = 2._wp * dzy (ji ,jj+1,jk) * dzy (ji,jj,jk ) 727 747 !! mb zdrho_k, zdz_k, zdrho_i, zdz_i, zdrho_j, zdz_j re-centred about the point (ji,jj,jk) 748 zdrho_k(:,:,:) = 0._wp 749 zdz_k (:,:,:) = 0._wp 750 751 DO_3D( 1, 1, 1, 1, 2, jpk-2 ) 752 cffw = 2._wp * zdrhoz(ji ,jj ,jk) * zdrhoz(ji,jj,jk+1) 728 753 IF( cffw > zep) THEN 729 drhow(ji,jj,jk) = 2._wp * drhoz(ji,jj,jk) * drhoz(ji,jj,jk-1) & 730 & / ( drhoz(ji,jj,jk) + drhoz(ji,jj,jk-1) ) 754 zdrho_k(ji,jj,jk) = cffw / ( zdrhoz(ji,jj,jk) + zdrhoz(ji,jj,jk+1) ) 755 ENDIF 756 zdz_k(ji,jj,jk) = 2._wp * zdzz(ji,jj,jk) * zdzz(ji,jj,jk+1) & 757 & / ( zdzz(ji,jj,jk) + zdzz(ji,jj,jk+1) ) 758 END_3D 759 760 !---------------------------------------------------------------------------------- 761 ! 3. apply boundary conditions at top and bottom using 5.36-5.37 762 !---------------------------------------------------------------------------------- 763 764 ! mb for sea-ice shelves we will need to re-write this upper boundary condition in the same form as the lower boundary condition 765 zdrho_k(:,:,1) = aco_bc_vrt * ( rhd_opt (:,:,2) - rhd_opt (:,:,1) ) - bco_bc_vrt * zdrho_k(:,:,2) 766 zdz_k (:,:,1) = aco_bc_vrt * (-gde3w(:,:,2) + gde3w(:,:,1) ) - bco_bc_vrt * zdz_k (:,:,2) 767 768 DO_2D( 1, 1, 1, 1 ) 769 IF ( mbkt(ji,jj)>1 ) THEN 770 iktb = mbkt(ji,jj) 771 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) 772 zdz_k (ji,jj,iktb) = aco_bc_vrt * (-gde3w(ji,jj,iktb) + gde3w(ji,jj,iktb-1) ) - bco_bc_vrt * zdz_k (ji,jj,iktb-1) 773 END IF 774 END_2D 775 776 !-------------------------------------------------------------- 777 ! 4. Compute side face integrals 778 !------------------------------------------------------------- 779 780 !! ssh replaces e3w_n ; gde3w is a depth; the formulae involve heights 781 !! rho_k stores grav * FX / rho_0 782 783 !-------------------------------------------------------------- 784 ! 4. a) Upper half of top-most grid box, compute and store 785 !------------------------------------------------------------- 786 ! *** AY note: ssh(ji,jj,Kmm) + gde3w(ji,jj,1) = e3w(ji,jj,1) 787 DO_2D( 0, 1, 0, 1) 788 z_rho_k(ji,jj,1) = grav * ( ssh(ji,jj,Kmm) + gde3w(ji,jj,1) ) & 789 & * ( rhd_opt(ji,jj,1) & 790 & + 0.5_wp * ( rhd_opt (ji,jj,2) - rhd_opt (ji,jj,1) ) & 791 & * ( ssh (ji,jj,Kmm) + gde3w(ji,jj,1) ) & 792 & / ( - gde3w(ji,jj,2) + gde3w(ji,jj,1) ) ) 793 END_2D 794 795 !-------------------------------------------------------------- 796 ! 4. b) Interior faces, compute and store 797 !------------------------------------------------------------- 798 799 DO_3D( 0, 1, 0, 1, 2, jpkm1 ) 800 z_rho_k(ji,jj,jk) = zcoef0 * ( rhd_opt (ji,jj,jk) + rhd_opt (ji,jj,jk-1) ) & 801 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) ) & 802 & + z_grav_10 * ( & 803 & ( zdrho_k (ji,jj,jk) - zdrho_k (ji,jj,jk-1) ) & 804 & * ( - gde3w(ji,jj,jk) + gde3w(ji,jj,jk-1) - z1_12 * ( zdz_k (ji,jj,jk) + zdz_k (ji,jj,jk-1) ) ) & 805 & - ( zdz_k (ji,jj,jk) - zdz_k (ji,jj,jk-1) ) & 806 & * ( rhd_opt (ji,jj,jk) - rhd_opt (ji,jj,jk-1) - z1_12 * ( zdrho_k(ji,jj,jk) + zdrho_k(ji,jj,jk-1) ) ) & 807 & ) 808 END_3D 809 810 !---------------------------------------------------------------------------------------- 811 ! 5. compute and store elementary horizontal differences in provisional arrays 812 !---------------------------------------------------------------------------------------- 813 814 DO_3D( 1, 0, 1, 0, 1, jpkm1 ) 815 zdrhox(ji,jj,jk) = rhd_opt (ji+1,jj ,jk) - rhd_opt (ji,jj,jk ) 816 zdzx (ji,jj,jk) = - gde3w(ji+1,jj ,jk) + gde3w(ji,jj,jk ) 817 zdrhoy(ji,jj,jk) = rhd_opt (ji ,jj+1,jk) - rhd_opt (ji,jj,jk ) 818 zdzy (ji,jj,jk) = - gde3w(ji ,jj+1,jk) + gde3w(ji,jj,jk ) 819 END_3D 820 821 CALL lbc_lnk_multi( 'dynhpg', zdrhox, 'U', 1., zdzx, 'U', 1., zdrhoy, 'V', 1., zdzy, 'V', 1. ) 822 823 !------------------------------------------------------------------------- 824 ! 6. compute harmonic averages using eq. 5.18 825 !------------------------------------------------------------------------- 826 827 DO_3D( 0, 1, 0, 1, 1, jpkm1 ) 828 cffu = 2._wp * zdrhox(ji-1,jj ,jk) * zdrhox(ji,jj,jk ) 829 IF( cffu > zep ) THEN 830 zdrho_i(ji,jj,jk) = cffu / ( zdrhox(ji-1,jj,jk) + zdrhox(ji,jj,jk) ) 731 831 ELSE 732 drhow(ji,jj,jk) = 0._wp 733 ENDIF 734 735 dzw(ji,jj,jk) = 2._wp * dzz(ji,jj,jk) * dzz(ji,jj,jk-1) & 736 & / ( dzz(ji,jj,jk) + dzz(ji,jj,jk-1) ) 737 738 IF( cffu > zep ) THEN 739 drhou(ji,jj,jk) = 2._wp * drhox(ji+1,jj,jk) * drhox(ji,jj,jk) & 740 & / ( drhox(ji+1,jj,jk) + drhox(ji,jj,jk) ) 832 zdrho_i(ji,jj,jk ) = 0._wp 833 ENDIF 834 835 cffx = 2._wp * zdzx (ji-1,jj ,jk) * zdzx (ji,jj,jk ) 836 IF( cffx > zep ) THEN 837 zdz_i(ji,jj,jk) = cffx / ( zdzx(ji-1,jj,jk) + zdzx(ji,jj,jk) ) 741 838 ELSE 742 drhou(ji,jj,jk) = 0._wp743 ENDIF 744 745 IF( cffx > zep ) THEN746 dzu(ji,jj,jk) = 2._wp * dzx(ji+1,jj,jk) * dzx(ji,jj,jk) &747 & / ( dzx(ji+1,jj,jk) + dzx(ji,jj,jk) )839 zdz_i(ji,jj,jk) = 0._wp 840 ENDIF 841 842 cffv = 2._wp * zdrhoy(ji ,jj-1,jk) * zdrhoy(ji,jj,jk ) 843 IF( cffv > zep ) THEN 844 zdrho_j(ji,jj,jk) = cffv / ( zdrhoy(ji,jj-1,jk) + zdrhoy(ji,jj,jk) ) 748 845 ELSE 749 dzu(ji,jj,jk) = 0._wp750 ENDIF 751 752 IF( cffv > zep ) THEN753 drhov(ji,jj,jk) = 2._wp * drhoy(ji,jj+1,jk) * drhoy(ji,jj,jk) &754 & / ( drhoy(ji,jj+1,jk) + drhoy(ji,jj,jk) )846 zdrho_j(ji,jj,jk) = 0._wp 847 ENDIF 848 849 cffy = 2._wp * zdzy (ji ,jj-1,jk) * zdzy (ji,jj,jk ) 850 IF( cffy > zep ) THEN 851 zdz_j(ji,jj,jk) = cffy / ( zdzy(ji,jj-1,jk) + zdzy(ji,jj,jk) ) 755 852 ELSE 756 drhov(ji,jj,jk) = 0._wp 757 ENDIF 758 759 IF( cffy > zep ) THEN 760 dzv(ji,jj,jk) = 2._wp * dzy(ji,jj+1,jk) * dzy(ji,jj,jk) & 761 & / ( dzy(ji,jj+1,jk) + dzy(ji,jj,jk) ) 762 ELSE 763 dzv(ji,jj,jk) = 0._wp 764 ENDIF 765 766 END_3D 853 zdz_j(ji,jj,jk) = 0._wp 854 ENDIF 855 END_3D 856 857 !!! Note that zdzx, zdzy, zdzz, zdrhox, zdrhoy and zdrhoz should NOT be used beyond this point 767 858 768 859 !---------------------------------------------------------------------------------- 769 ! apply boundary conditions at top and bottomusing 5.36-5.37860 ! 6B. apply boundary conditions at side boundaries using 5.36-5.37 770 861 !---------------------------------------------------------------------------------- 771 drhow(:,:, 1 ) = 1.5_wp * ( drhoz(:,:, 2 ) - drhoz(:,:, 1 ) ) - 0.5_wp * drhow(:,:, 2 ) 772 drhou(:,:, 1 ) = 1.5_wp * ( drhox(:,:, 2 ) - drhox(:,:, 1 ) ) - 0.5_wp * drhou(:,:, 2 ) 773 drhov(:,:, 1 ) = 1.5_wp * ( drhoy(:,:, 2 ) - drhoy(:,:, 1 ) ) - 0.5_wp * drhov(:,:, 2 ) 774 775 drhow(:,:,jpk) = 1.5_wp * ( drhoz(:,:,jpk) - drhoz(:,:,jpkm1) ) - 0.5_wp * drhow(:,:,jpkm1) 776 drhou(:,:,jpk) = 1.5_wp * ( drhox(:,:,jpk) - drhox(:,:,jpkm1) ) - 0.5_wp * drhou(:,:,jpkm1) 777 drhov(:,:,jpk) = 1.5_wp * ( drhoy(:,:,jpk) - drhoy(:,:,jpkm1) ) - 0.5_wp * drhov(:,:,jpkm1) 778 862 863 DO jk = 1, jpkm1 864 zz_drho_i(:,:) = zdrho_i(:,:,jk) 865 zz_dz_i (:,:) = zdz_i (:,:,jk) 866 zz_drho_j(:,:) = zdrho_j(:,:,jk) 867 zz_dz_j (:,:) = zdz_j (:,:,jk) 868 DO_2D( 0, 1, 0, 1) 869 ! Walls coming from left: should check from 2 to jpi-1 (and jpj=2-jpj) 870 IF (ji < jpi) THEN 871 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 872 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) 873 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) 874 END IF 875 END IF 876 ! Walls coming from right: should check from 3 to jpi (and jpj=2-jpj) 877 IF (ji > 2) THEN 878 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 879 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) 880 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) 881 END IF 882 END IF 883 ! Walls coming from left: should check from 2 to jpj-1 (and jpi=2-jpi) 884 IF (jj < jpj) THEN 885 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 886 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) 887 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) 888 END IF 889 END IF 890 ! Walls coming from right: should check from 3 to jpj (and jpi=2-jpi) 891 IF (jj > 2) THEN 892 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 893 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) 894 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) 895 END IF 896 END IF 897 END_2D 898 zdrho_i(:,:,jk) = zz_drho_i(:,:) 899 zdz_i (:,:,jk) = zz_dz_i (:,:) 900 zdrho_j(:,:,jk) = zz_drho_j(:,:) 901 zdz_j (:,:,jk) = zz_dz_j (:,:) 902 END DO 779 903 780 904 !-------------------------------------------------------------- 781 ! Upper half of top-most grid box, compute and store905 ! 7. Calculate integrals on side faces 782 906 !------------------------------------------------------------- 783 907 784 !!bug gm : e3w-gde3w(:,:,:) = 0.5*e3w .... and gde3w(:,:,2)-gde3w(:,:,1)=e3w(:,:,2,Kmm) .... to be verified 785 ! true if gde3w(:,:,:) is really defined as the sum of the e3w scale factors as, it seems to me, it should be 786 787 DO_2D( 0, 0, 0, 0 ) 788 rho_k(ji,jj,1) = -grav * ( e3w(ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 789 & * ( rhd(ji,jj,1) & 790 & + 0.5_wp * ( rhd (ji,jj,2) - rhd (ji,jj,1) ) & 791 & * ( e3w (ji,jj,1,Kmm) - gde3w(ji,jj,1) ) & 792 & / ( gde3w(ji,jj,2) - gde3w(ji,jj,1) ) ) 793 END_2D 794 795 !!bug gm : here also, simplification is possible 796 !!bug gm : optimisation: 1/10 and 1/12 the division should be done before the loop 797 798 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 799 800 rho_k(ji,jj,jk) = zcoef0 * ( rhd (ji,jj,jk) + rhd (ji,jj,jk-1) ) & 801 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) ) & 802 & - grav * z1_10 * ( & 803 & ( drhow (ji,jj,jk) - drhow (ji,jj,jk-1) ) & 804 & * ( gde3w(ji,jj,jk) - gde3w(ji,jj,jk-1) - z1_12 * ( dzw (ji,jj,jk) + dzw (ji,jj,jk-1) ) ) & 805 & - ( dzw (ji,jj,jk) - dzw (ji,jj,jk-1) ) & 806 & * ( rhd (ji,jj,jk) - rhd (ji,jj,jk-1) - z1_12 * ( drhow(ji,jj,jk) + drhow(ji,jj,jk-1) ) ) & 807 & ) 808 809 rho_i(ji,jj,jk) = zcoef0 * ( rhd (ji+1,jj,jk) + rhd (ji,jj,jk) ) & 810 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) & 811 & - grav* z1_10 * ( & 812 & ( drhou (ji+1,jj,jk) - drhou (ji,jj,jk) ) & 813 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzu (ji+1,jj,jk) + dzu (ji,jj,jk) ) ) & 814 & - ( dzu (ji+1,jj,jk) - dzu (ji,jj,jk) ) & 815 & * ( rhd (ji+1,jj,jk) - rhd (ji,jj,jk) - z1_12 * ( drhou(ji+1,jj,jk) + drhou(ji,jj,jk) ) ) & 816 & ) 817 818 rho_j(ji,jj,jk) = zcoef0 * ( rhd (ji,jj+1,jk) + rhd (ji,jj,jk) ) & 819 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) & 820 & - grav* z1_10 * ( & 821 & ( drhov (ji,jj+1,jk) - drhov (ji,jj,jk) ) & 822 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) - z1_12 * ( dzv (ji,jj+1,jk) + dzv (ji,jj,jk) ) ) & 823 & - ( dzv (ji,jj+1,jk) - dzv (ji,jj,jk) ) & 824 & * ( rhd (ji,jj+1,jk) - rhd (ji,jj,jk) - z1_12 * ( drhov(ji,jj+1,jk) + drhov(ji,jj,jk) ) ) & 825 & ) 826 827 END_3D 828 CALL lbc_lnk_multi( 'dynhpg', rho_k, 'W', 1.0_wp, rho_i, 'U', 1.0_wp, rho_j, 'V', 1.0_wp ) 908 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 909 ! two -ve signs cancel in next two lines (within zcoef0 and because gde3w is a depth not a height) 910 z_rho_i(ji,jj,jk) = zcoef0 * ( rhd_opt (ji+1,jj,jk) + rhd_opt (ji,jj,jk) ) & 911 & * ( gde3w(ji+1,jj,jk) - gde3w(ji,jj,jk) ) 912 IF ( umask(ji-1, jj, jk) > 0.5 .OR. umask(ji+1, jj, jk) > 0.5 ) THEN 913 z_rho_i(ji,jj,jk) = z_rho_i(ji,jj,jk) - z_grav_10 * ( & 914 & ( zdrho_i (ji+1,jj,jk) - zdrho_i (ji,jj,jk) ) & 915 & * ( - gde3w(ji+1,jj,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_i (ji+1,jj,jk) + zdz_i (ji,jj,jk) ) ) & 916 & - ( zdz_i (ji+1,jj,jk) - zdz_i (ji,jj,jk) ) & 917 & * ( rhd_opt (ji+1,jj,jk) - rhd_opt (ji,jj,jk) - z1_12 * ( zdrho_i(ji+1,jj,jk) + zdrho_i(ji,jj,jk) ) ) & 918 & ) 919 END IF 920 921 z_rho_j(ji,jj,jk) = zcoef0 * ( rhd_opt (ji,jj+1,jk) + rhd_opt (ji,jj,jk) ) & 922 & * ( gde3w(ji,jj+1,jk) - gde3w(ji,jj,jk) ) 923 IF ( vmask(ji, jj-1, jk) > 0.5 .OR. vmask(ji, jj+1, jk) > 0.5 ) THEN 924 z_rho_j(ji,jj,jk) = z_rho_j(ji,jj,jk) - z_grav_10 * ( & 925 & ( zdrho_j (ji,jj+1,jk) - zdrho_j (ji,jj,jk) ) & 926 & * ( - gde3w(ji,jj+1,jk) + gde3w(ji,jj,jk) - z1_12 * ( zdz_j (ji,jj+1,jk) + zdz_j (ji,jj,jk) ) ) & 927 & - ( zdz_j (ji,jj+1,jk) - zdz_j (ji,jj,jk) ) & 928 & * ( rhd_opt (ji,jj+1,jk) - rhd_opt (ji,jj,jk) - z1_12 * ( zdrho_j(ji,jj+1,jk) + zdrho_j(ji,jj,jk) ) ) & 929 & ) 930 END IF 931 END_3D 932 933 !-------------------------------------------------------------- 934 ! 8. Integrate in the vertical 935 !------------------------------------------------------------- 829 936 830 937 ! --------------- … … 832 939 ! --------------- 833 940 DO_2D( 0, 0, 0, 0 ) 834 zhpi(ji,jj,1) = ( rho_k(ji+1,jj ,1) - rho_k(ji,jj,1) -rho_i(ji,jj,1) ) * r1_e1u(ji,jj)835 zhpj(ji,jj,1) = ( rho_k(ji ,jj+1,1) - rho_k(ji,jj,1) -rho_j(ji,jj,1) ) * r1_e2v(ji,jj)941 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) 942 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) 836 943 IF( ln_wd_il ) THEN 837 944 zhpi(ji,jj,1) = zhpi(ji,jj,1) * zcpx(ji,jj) … … 848 955 DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 849 956 ! hydrostatic pressure gradient along s-surfaces 850 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) &851 & + ( ( rho_k(ji+1,jj,jk) - rho_k(ji,jj,jk ) )&852 & - ( rho_i(ji ,jj,jk) - rho_i(ji,jj,jk-1) ) ) * r1_e1u(ji,jj)853 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) &854 & + ( ( rho_k(ji,jj+1,jk) - rho_k(ji,jj,jk ) )&855 & -( rho_j(ji,jj ,jk) - rho_j(ji,jj,jk-1) ) ) * r1_e2v(ji,jj)957 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) & 958 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji+1,jj,jk ) ) & 959 & - ( z_rho_i(ji,jj,jk) - z_rho_i(ji ,jj,jk-1) ) ) * r1_e1u(ji,jj) 960 zhpj(ji,jj,jk) = zhpj(ji,jj,jk-1) & 961 & + ( ( z_rho_k(ji,jj,jk) - z_rho_k(ji,jj+1,jk ) ) & 962 & -( z_rho_j(ji,jj,jk) - z_rho_j(ji,jj ,jk-1) ) ) * r1_e2v(ji,jj) 856 963 IF( ln_wd_il ) THEN 857 964 zhpi(ji,jj,jk) = zhpi(ji,jj,jk) * zcpx(ji,jj)
Note: See TracChangeset
for help on using the changeset viewer.
