- Timestamp:
- 2015-04-07T16:22:54+02:00 (9 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r5151_UKMO_ISF/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5189 r5200 365 365 !! - bathy : meter bathymetry (in meters) 366 366 !!---------------------------------------------------------------------- 367 INTEGER :: ji, jj, j l, jk ! dummy loop indices367 INTEGER :: ji, jj, jk ! dummy loop indices 368 368 INTEGER :: inum ! temporary logical unit 369 369 INTEGER :: ierror ! error flag … … 973 973 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 974 974 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 975 REAL(wp) :: zmax ! Maximum depth976 975 REAL(wp) :: zdiff ! temporary scalar 977 REAL(wp) :: z refdep! temporary scalar976 REAL(wp) :: zmax ! temporary scalar 978 977 REAL(wp), POINTER, DIMENSION(:,:,:) :: zprt 979 978 !!--------------------------------------------------------------------- … … 1014 1013 END DO 1015 1014 1016 IF ( ln_isfcav ) CALL zgr_isf1017 1018 1015 ! Scale factors and depth at T- and W-points 1019 1016 DO jk = 1, jpk ! intitialization to the reference z-coordinate … … 1023 1020 e3w_0 (:,:,jk) = e3w_1d (jk) 1024 1021 END DO 1025 ! 1026 DO jj = 1, jpj 1027 DO ji = 1, jpi 1028 ik = mbathy(ji,jj) 1029 IF( ik > 0 ) THEN ! ocean point only 1030 ! max ocean level case 1031 IF( ik == jpkm1 ) THEN 1032 zdepwp = bathy(ji,jj) 1033 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1034 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1035 e3t_0(ji,jj,ik ) = ze3tp 1036 e3t_0(ji,jj,ik+1) = ze3tp 1037 e3w_0(ji,jj,ik ) = ze3wp 1038 e3w_0(ji,jj,ik+1) = ze3tp 1039 gdepw_0(ji,jj,ik+1) = zdepwp 1040 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1041 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1042 ! 1043 ELSE ! standard case 1044 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1045 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1022 1023 ! Bathy, iceshelf draft, scale factor and depth at T- and W- points in case of isf 1024 IF ( ln_isfcav ) CALL zgr_isf 1025 1026 ! Scale factors and depth at T- and W-points 1027 IF ( .NOT. ln_isfcav ) THEN 1028 DO jj = 1, jpj 1029 DO ji = 1, jpi 1030 ik = mbathy(ji,jj) 1031 IF( ik > 0 ) THEN ! ocean point only 1032 ! max ocean level case 1033 IF( ik == jpkm1 ) THEN 1034 zdepwp = bathy(ji,jj) 1035 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1036 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1037 e3t_0(ji,jj,ik ) = ze3tp 1038 e3t_0(ji,jj,ik+1) = ze3tp 1039 e3w_0(ji,jj,ik ) = ze3wp 1040 e3w_0(ji,jj,ik+1) = ze3tp 1041 gdepw_0(ji,jj,ik+1) = zdepwp 1042 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1043 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1044 ! 1045 ELSE ! standard case 1046 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1047 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1048 ENDIF 1049 !gm Bug? check the gdepw_1d 1050 ! ... on ik 1051 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1052 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1053 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1054 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) & 1055 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ) 1056 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) & 1057 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) ) 1058 ! ... on ik+1 1059 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1060 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1061 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik) 1046 1062 ENDIF 1047 !gm Bug? check the gdepw_1d1048 ! ... on ik1049 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) &1050 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) &1051 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) ))1052 e3t_0 (ji,jj,ik) = e3t_1d (ik) * ( gdepw_0 (ji,jj,ik+1) - gdepw_1d(ik) ) &1053 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )1054 e3w_0(ji,jj,ik) = 0.5_wp * ( gdepw_0(ji,jj,ik+1) + gdepw_1d(ik+1) - 2._wp * gdepw_1d(ik) ) &1055 & * ( e3w_1d(ik) / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) )1056 ! ... on ik+11057 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1058 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik)1059 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + e3t_0(ji,jj,ik)1060 1063 ENDIF 1061 END IF1062 END DO 1063 END DO1064 !1065 it = 01066 DO jj = 1, jpj1067 DO ji = 1, jpi1068 ik = mbathy(ji,jj)1069 IF( ik > 0 ) THEN ! ocean point only1070 e3tp (ji,jj) = e3t_0(ji,jj,ik)1071 e3wp (ji,jj) = e3w_0(ji,jj,ik)1072 ! test1073 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik )1074 IF( zdiff <= 0._wp .AND. lwp ) THEN1075 it = it + 11076 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj1077 WRITE(numout,*) ' bathy = ', bathy(ji,jj)1078 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff1079 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik )1064 END DO 1065 END DO 1066 ! 1067 it = 0 1068 DO jj = 1, jpj 1069 DO ji = 1, jpi 1070 ik = mbathy(ji,jj) 1071 IF( ik > 0 ) THEN ! ocean point only 1072 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1073 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1074 ! test 1075 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1076 IF( zdiff <= 0._wp .AND. lwp ) THEN 1077 it = it + 1 1078 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1079 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1080 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1081 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1082 ENDIF 1080 1083 ENDIF 1081 ENDIF 1082 END DO 1083 END DO 1084 ! 1085 IF ( ln_isfcav ) THEN 1086 ! (ISF) Definition of e3t, u, v, w for ISF case 1087 DO jj = 1, jpj 1088 DO ji = 1, jpi 1089 ik = misfdep(ji,jj) 1090 IF( ik > 1 ) THEN ! ice shelf point only 1091 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1092 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1093 !gm Bug? check the gdepw_0 1094 ! ... on ik 1095 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1096 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1097 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1098 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1099 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1100 1101 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1102 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1103 ENDIF 1104 ! ... on ik / ik-1 1105 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1106 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1107 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1108 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1109 ENDIF 1110 END DO 1111 END DO 1112 ! 1113 it = 0 1114 DO jj = 1, jpj 1115 DO ji = 1, jpi 1116 ik = misfdep(ji,jj) 1117 IF( ik > 1 ) THEN ! ice shelf point only 1118 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1119 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1120 ! test 1121 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1122 IF( zdiff <= 0. .AND. lwp ) THEN 1123 it = it + 1 1124 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1125 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1126 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1127 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1128 ENDIF 1129 ENDIF 1130 END DO 1131 END DO 1084 END DO 1085 END DO 1132 1086 END IF 1133 ! END (ISF) 1134 1087 ! 1135 1088 ! Scale factors and depth at U-, V-, UW and VW-points 1136 1089 DO jk = 1, jpk ! initialisation to z-scale factors … … 1275 1228 !!---------------------------------------------------------------------- 1276 1229 !! 1277 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1278 INTEGER :: ik, it ! temporary integers 1279 INTEGER :: id, jd, nprocd 1230 INTEGER :: ji, jj, jl, jk ! dummy loop indices 1231 INTEGER :: ik, it ! temporary integers 1280 1232 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1281 LOGICAL :: ll_print ! Allow control print for debugging 1233 REAL(wp) :: zdepth ! Ajusted ocean depth to avoid too small e3t 1234 REAL(wp) :: zmax ! Maximum and minimum depth 1235 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar 1282 1236 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1283 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t 1284 REAL(wp) :: zmax, zmin ! Maximum and minimum depth 1237 REAL(wp) :: zdepwp ! Ajusted ocean depth to avoid too small e3t 1285 1238 REAL(wp) :: zdiff ! temporary scalar 1286 REAL(wp) :: zrefdep ! temporary scalar1287 REAL(wp) :: zbathydiff, zrisfdepdiff ! isf temporary scalar1288 1239 REAL(wp), POINTER, DIMENSION(:,:) :: zrisfdep, zbathy, zmask ! 2D workspace (ISH) 1289 1240 INTEGER , POINTER, DIMENSION(:,:) :: zmbathy, zmisfdep ! 2D workspace (ISH) … … 1759 1710 ENDIF 1760 1711 1712 ! compute scale factor and depth at T- and W- points 1713 DO jj = 1, jpj 1714 DO ji = 1, jpi 1715 ik = mbathy(ji,jj) 1716 IF( ik > 0 ) THEN ! ocean point only 1717 ! max ocean level case 1718 IF( ik == jpkm1 ) THEN 1719 zdepwp = bathy(ji,jj) 1720 ze3tp = bathy(ji,jj) - gdepw_1d(ik) 1721 ze3wp = 0.5_wp * e3w_1d(ik) * ( 1._wp + ( ze3tp/e3t_1d(ik) ) ) 1722 e3t_0(ji,jj,ik ) = ze3tp 1723 e3t_0(ji,jj,ik+1) = ze3tp 1724 e3w_0(ji,jj,ik ) = ze3wp 1725 e3w_0(ji,jj,ik+1) = ze3tp 1726 gdepw_0(ji,jj,ik+1) = zdepwp 1727 gdept_0(ji,jj,ik ) = gdept_1d(ik-1) + ze3wp 1728 gdept_0(ji,jj,ik+1) = gdept_0(ji,jj,ik) + ze3tp 1729 ! 1730 ELSE ! standard case 1731 IF( bathy(ji,jj) <= gdepw_1d(ik+1) ) THEN ; gdepw_0(ji,jj,ik+1) = bathy(ji,jj) 1732 ELSE ; gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1733 ENDIF 1734 ! gdepw_0(ji,jj,ik+1) = gdepw_1d(ik+1) 1735 !gm Bug? check the gdepw_1d 1736 ! ... on ik 1737 gdept_0(ji,jj,ik) = gdepw_1d(ik) + ( gdepw_0(ji,jj,ik+1) - gdepw_1d(ik) ) & 1738 & * ((gdept_1d( ik ) - gdepw_1d(ik) ) & 1739 & / ( gdepw_1d( ik+1) - gdepw_1d(ik) )) 1740 e3t_0 (ji,jj,ik ) = gdepw_0(ji,jj,ik+1) - gdepw_1d(ik ) 1741 e3w_0 (ji,jj,ik ) = gdept_0(ji,jj,ik ) - gdept_1d(ik-1) 1742 ! ... on ik+1 1743 e3w_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1744 e3t_0 (ji,jj,ik+1) = e3t_0 (ji,jj,ik) 1745 ENDIF 1746 ENDIF 1747 END DO 1748 END DO 1749 ! 1750 it = 0 1751 DO jj = 1, jpj 1752 DO ji = 1, jpi 1753 ik = mbathy(ji,jj) 1754 IF( ik > 0 ) THEN ! ocean point only 1755 e3tp (ji,jj) = e3t_0(ji,jj,ik) 1756 e3wp (ji,jj) = e3w_0(ji,jj,ik) 1757 ! test 1758 zdiff= gdepw_0(ji,jj,ik+1) - gdept_0(ji,jj,ik ) 1759 IF( zdiff <= 0._wp .AND. lwp ) THEN 1760 it = it + 1 1761 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1762 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 1763 WRITE(numout,*) ' gdept_0 = ', gdept_0(ji,jj,ik), ' gdepw_0 = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1764 WRITE(numout,*) ' e3tp = ', e3t_0 (ji,jj,ik), ' e3wp = ', e3w_0 (ji,jj,ik ) 1765 ENDIF 1766 ENDIF 1767 END DO 1768 END DO 1769 ! 1770 ! (ISF) Definition of e3t, u, v, w for ISF case 1771 DO jj = 1, jpj 1772 DO ji = 1, jpi 1773 ik = misfdep(ji,jj) 1774 IF( ik > 1 ) THEN ! ice shelf point only 1775 IF( risfdep(ji,jj) < gdepw_1d(ik) ) risfdep(ji,jj)= gdepw_1d(ik) 1776 gdepw_0(ji,jj,ik) = risfdep(ji,jj) 1777 !gm Bug? check the gdepw_0 1778 ! ... on ik 1779 gdept_0(ji,jj,ik) = gdepw_1d(ik+1) - ( gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) ) & 1780 & * ( gdepw_1d(ik+1) - gdept_1d(ik) ) & 1781 & / ( gdepw_1d(ik+1) - gdepw_1d(ik) ) 1782 e3t_0 (ji,jj,ik ) = gdepw_1d(ik+1) - gdepw_0(ji,jj,ik) 1783 e3w_0 (ji,jj,ik+1) = gdept_1d(ik+1) - gdept_0(ji,jj,ik) 1784 1785 IF( ik + 1 == mbathy(ji,jj) ) THEN ! ice shelf point only (2 cell water column) 1786 e3w_0 (ji,jj,ik+1) = gdept_0(ji,jj,ik+1) - gdept_0(ji,jj,ik) 1787 ENDIF 1788 ! ... on ik / ik-1 1789 e3w_0 (ji,jj,ik ) = 2._wp * (gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik)) 1790 e3t_0 (ji,jj,ik-1) = gdepw_0(ji,jj,ik) - gdepw_1d(ik-1) 1791 ! The next line isn't required and doesn't affect results - included for consistency with bathymetry code 1792 gdept_0(ji,jj,ik-1) = gdept_1d(ik-1) 1793 ENDIF 1794 END DO 1795 END DO 1796 1797 it = 0 1798 DO jj = 1, jpj 1799 DO ji = 1, jpi 1800 ik = misfdep(ji,jj) 1801 IF( ik > 1 ) THEN ! ice shelf point only 1802 e3tp (ji,jj) = e3t_0(ji,jj,ik ) 1803 e3wp (ji,jj) = e3w_0(ji,jj,ik+1 ) 1804 ! test 1805 zdiff= gdept_0(ji,jj,ik) - gdepw_0(ji,jj,ik ) 1806 IF( zdiff <= 0. .AND. lwp ) THEN 1807 it = it + 1 1808 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 1809 WRITE(numout,*) ' risfdep = ', risfdep(ji,jj) 1810 WRITE(numout,*) ' gdept = ', gdept_0(ji,jj,ik), ' gdepw = ', gdepw_0(ji,jj,ik+1), ' zdiff = ', zdiff 1811 WRITE(numout,*) ' e3tp = ', e3tp(ji,jj), ' e3wp = ', e3wp(ji,jj) 1812 ENDIF 1813 ENDIF 1814 END DO 1815 END DO 1816 1761 1817 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1762 1818 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy )
Note: See TracChangeset
for help on using the changeset viewer.
