- Timestamp:
- 2012-07-20T11:18:36+02:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90 (modified) (16 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3294 r3434 678 678 REAL(wp) :: zrhdt1 679 679 REAL(wp) :: zdpdx1, zdpdx2, zdpdy1, zdpdy2 680 INTEGER :: zbhitwe, zbhitns 681 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdeptht, zrhh 680 REAL(wp), POINTER, DIMENSION(:,:,:) :: zdept, zrhh 682 681 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp 683 682 !!---------------------------------------------------------------------- 684 683 ! 685 684 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 686 CALL wrk_alloc( jpi,jpj,jpk, zdept ht, zrhh )685 CALL wrk_alloc( jpi,jpj,jpk, zdept, zrhh ) 687 686 ! 688 687 IF( kt == nit000 ) THEN … … 717 716 END DO 718 717 719 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 720 DO jj = 1, jpj 721 DO ji = 1, jpi 722 zdeptht(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 723 zdeptht(ji,jj,1) = zdeptht(ji,jj,1) - sshn(ji,jj) * znad 724 DO jk = 2, jpk 725 zdeptht(ji,jj,jk) = zdeptht(ji,jj,jk-1) + fse3w(ji,jj,jk) 726 END DO 727 END DO 728 END DO 729 730 DO jk = 1, jpkm1 731 DO jj = 1, jpj 732 DO ji = 1, jpi 733 fsp(ji,jj,jk) = zrhh(ji,jj,jk) 734 xsp(ji,jj,jk) = zdeptht(ji,jj,jk) 735 END DO 736 END DO 737 END DO 718 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdept(:,:,:)" 719 DO jj = 1, jpj; DO ji = 1, jpi 720 zdept(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) - sshn(ji,jj) * znad 721 END DO ; END DO 722 723 DO jk = 2, jpk; DO jj = 1, jpj; DO ji = 1, jpi 724 zdept(ji,jj,jk) = zdept(ji,jj,jk-1) + fse3w(ji,jj,jk) 725 END DO ; END DO ; END DO 726 727 fsp(:,:,:) = zrhh(:,:,:) 728 xsp(:,:,:) = zdept(:,:,:) 738 729 739 730 ! Construct the vertical density profile with the … … 745 736 DO jj = 2, jpj 746 737 DO ji = 2, jpi 747 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept ht(ji,jj,1),asp(ji,jj,1), &738 zrhdt1 = zrhh(ji,jj,1) - interp3(zdept(ji,jj,1),asp(ji,jj,1), & 748 739 bsp(ji,jj,1), csp(ji,jj,1), & 749 dsp(ji,jj,1) ) * 0.5_wp * zdeptht(ji,jj,1) 750 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water 740 dsp(ji,jj,1) ) * 0.25_wp * fse3w(ji,jj,1) 751 741 752 742 ! assuming linear profile across the top half surface layer … … 760 750 DO ji = 2, jpi 761 751 zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 762 integ 2(zdeptht(ji,jj,jk-1), zdeptht(ji,jj,jk),&752 integ_spline(zdept(ji,jj,jk-1), zdept(ji,jj,jk),& 763 753 asp(ji,jj,jk-1), bsp(ji,jj,jk-1), & 764 754 csp(ji,jj,jk-1), dsp(ji,jj,jk-1)) … … 793 783 END DO 794 784 785 DO jk = 1, jpkm1 786 DO jj = 2, jpjm1 787 DO ji = 2, jpim1 788 zu(ji,jj,jk) = min(zu(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk))) 789 zu(ji,jj,jk) = max(zu(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji+1,jj,jk))) 790 zv(ji,jj,jk) = min(zv(ji,jj,jk), max(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk))) 791 zv(ji,jj,jk) = max(zv(ji,jj,jk), min(-zdept(ji,jj,jk), -zdept(ji,jj+1,jk))) 792 END DO 793 END DO 794 END DO 795 796 795 797 DO jk = 1, jpkm1 796 798 DO jj = 2, jpjm1 … … 803 805 !!!!! for u equation 804 806 IF( jk <= mbku(ji,jj) ) THEN 805 IF( -zdept ht(ji+1,jj,mbku(ji,jj)) >= -zdeptht(ji,jj,mbku(ji,jj)) ) THEN807 IF( -zdept(ji+1,jj,jk) >= -zdept(ji,jj,jk) ) THEN 806 808 jis = ji + 1; jid = ji 807 809 ELSE … … 811 813 ! integrate the pressure on the shallow side 812 814 jk1 = jk 813 zbhitwe = 0 814 DO WHILE ( -zdeptht(jis,jj,jk1) > zuijk ) 815 DO WHILE ( -zdept(jis,jj,jk1) > zuijk ) 815 816 IF( jk1 == mbku(ji,jj) ) THEN 816 z bhitwe = 1817 zuijk = -zdept(jis,jj,jk1) 817 818 EXIT 818 819 ENDIF 819 zdeps = MIN(zdept ht(jis,jj,jk1+1), -zuijk)820 zdeps = MIN(zdept(jis,jj,jk1+1), -zuijk) 820 821 zpwes = zpwes + & 821 integ 2(zdeptht(jis,jj,jk1), zdeps, &822 integ_spline(zdept(jis,jj,jk1), zdeps, & 822 823 asp(jis,jj,jk1), bsp(jis,jj,jk1), & 823 824 csp(jis,jj,jk1), dsp(jis,jj,jk1)) … … 825 826 END DO 826 827 827 IF(zbhitwe == 1) THEN828 zuijk = -zdeptht(jis,jj,jk1)829 ENDIF830 831 828 ! integrate the pressure on the deep side 832 829 jk1 = jk 833 zbhitwe = 0 834 DO WHILE ( -zdeptht(jid,jj,jk1) < zuijk ) 830 DO WHILE ( -zdept(jid,jj,jk1) < zuijk ) 835 831 IF( jk1 == 1 ) THEN 836 zbhitwe = 1 832 zdeps = zdept(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad) 833 zrhdt1 = zrhh(jid,jj,1) - interp3(zdept(jid,jj,1), asp(jid,jj,1), & 834 bsp(jid,jj,1), csp(jid,jj,1), & 835 dsp(jid,jj,1)) * zdeps 836 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps 837 837 EXIT 838 838 ENDIF 839 zdeps = MAX(zdept ht(jid,jj,jk1-1), -zuijk)839 zdeps = MAX(zdept(jid,jj,jk1-1), -zuijk) 840 840 zpwed = zpwed + & 841 integ 2(zdeps, zdeptht(jid,jj,jk1), &841 integ_spline(zdeps, zdept(jid,jj,jk1), & 842 842 asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 843 843 csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1) ) … … 845 845 END DO 846 846 847 IF( zbhitwe == 1 ) THEN848 zdeps = zdeptht(jid,jj,1) + MIN(zuijk, sshn(jid,jj)*znad)849 zrhdt1 = zrhh(jid,jj,1) - interp3(zdeptht(jid,jj,1), asp(jid,jj,1), &850 bsp(jid,jj,1), csp(jid,jj,1), &851 dsp(jid,jj,1)) * zdeps852 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water853 zpwed = zpwed + 0.5_wp * (zrhh(jid,jj,1) + zrhdt1) * zdeps854 ENDIF855 856 847 ! update the momentum trends in u direction 857 848 … … 870 861 !!!!! for v equation 871 862 IF( jk <= mbkv(ji,jj) ) THEN 872 IF( -zdept ht(ji,jj+1,mbkv(ji,jj)) >= -zdeptht(ji,jj,mbkv(ji,jj)) ) THEN863 IF( -zdept(ji,jj+1,jk) >= -zdept(ji,jj,jk) ) THEN 873 864 jjs = jj + 1; jjd = jj 874 865 ELSE … … 878 869 ! integrate the pressure on the shallow side 879 870 jk1 = jk 880 zbhitns = 0 881 DO WHILE ( -zdeptht(ji,jjs,jk1) > zvijk ) 871 DO WHILE ( -zdept(ji,jjs,jk1) > zvijk ) 882 872 IF( jk1 == mbkv(ji,jj) ) THEN 883 z bhitns = 1873 zvijk = -zdept(ji,jjs,jk1) 884 874 EXIT 885 875 ENDIF 886 zdeps = MIN(zdept ht(ji,jjs,jk1+1), -zvijk)876 zdeps = MIN(zdept(ji,jjs,jk1+1), -zvijk) 887 877 zpnss = zpnss + & 888 integ 2(zdeptht(ji,jjs,jk1), zdeps, &878 integ_spline(zdept(ji,jjs,jk1), zdeps, & 889 879 asp(ji,jjs,jk1), bsp(ji,jjs,jk1), & 890 880 csp(ji,jjs,jk1), dsp(ji,jjs,jk1) ) … … 892 882 END DO 893 883 894 IF(zbhitns == 1) THEN895 zvijk = -zdeptht(ji,jjs,jk1)896 ENDIF897 898 884 ! integrate the pressure on the deep side 899 885 jk1 = jk 900 zbhitns = 0 901 DO WHILE ( -zdeptht(ji,jjd,jk1) < zvijk ) 886 DO WHILE ( -zdept(ji,jjd,jk1) < zvijk ) 902 887 IF( jk1 == 1 ) THEN 903 zbhitns = 1 888 zdeps = zdept(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad) 889 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdept(ji,jjd,1), asp(ji,jjd,1), & 890 bsp(ji,jjd,1), csp(ji,jjd,1), & 891 dsp(ji,jjd,1) ) * zdeps 892 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps 904 893 EXIT 905 894 ENDIF 906 zdeps = MAX(zdept ht(ji,jjd,jk1-1), -zvijk)895 zdeps = MAX(zdept(ji,jjd,jk1-1), -zvijk) 907 896 zpnsd = zpnsd + & 908 integ 2(zdeps, zdeptht(ji,jjd,jk1), &897 integ_spline(zdeps, zdept(ji,jjd,jk1), & 909 898 asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1), & 910 899 csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) … … 912 901 END DO 913 902 914 IF( zbhitns == 1 ) THEN915 zdeps = zdeptht(ji,jjd,1) + MIN(zvijk, sshn(ji,jjd)*znad)916 zrhdt1 = zrhh(ji,jjd,1) - interp3(zdeptht(ji,jjd,1), asp(ji,jjd,1), &917 bsp(ji,jjd,1), csp(ji,jjd,1), &918 dsp(ji,jjd,1) ) * zdeps919 zrhdt1 = MAX(zrhdt1, 1000._wp - rau0) ! no lighter than fresh water920 zpnsd = zpnsd + 0.5_wp * (zrhh(ji,jjd,1) + zrhdt1) * zdeps921 ENDIF922 903 923 904 ! update the momentum trends in v direction … … 941 922 ! 942 923 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 943 CALL wrk_dealloc( jpi,jpj,jpk, zdept ht, zrhh )924 CALL wrk_dealloc( jpi,jpj,jpk, zdept, zrhh ) 944 925 ! 945 926 END SUBROUTINE hpg_prj … … 1121 1102 1122 1103 1123 FUNCTION integ 2(xl, xr, a, b, c, d) RESULT(f)1104 FUNCTION integ_spline(xl, xr, a, b, c, d) RESULT(f) 1124 1105 !!---------------------------------------------------------------------- 1125 1106 !! *** ROUTINE interp1 *** … … 1143 1124 & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 1144 1125 1145 END FUNCTION integ 21126 END FUNCTION integ_spline 1146 1127 1147 1128 1148 1129 !!====================================================================== 1149 1130 END MODULE dynhpg 1131
Note: See TracChangeset
for help on using the changeset viewer.
