- Timestamp:
- 2015-12-16T10:25:22+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_merge_2015/NEMOGCM/NEMO/OPA_SRC/DOM/domzgr.F90
r5836 r6060 2 2 !!============================================================================== 3 3 !! *** MODULE domzgr *** 4 !! Ocean initialization : domain initialization4 !! Ocean domain : definition of the vertical coordinate system 5 5 !!============================================================================== 6 6 !! History : OPA ! 1995-12 (G. Madec) Original code : s vertical coordinate … … 38 38 USE closea ! closed seas 39 39 USE c1d ! 1D vertical configuration 40 ! 40 41 USE in_out_manager ! I/O manager 41 42 USE iom ! I/O library … … 73 74 74 75 !! * Substitutions 75 # include "domzgr_substitute.h90"76 76 # include "vectopt_loop_substitute.h90" 77 77 !!---------------------------------------------------------------------- … … 102 102 INTEGER :: ios 103 103 ! 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav 104 NAMELIST/namzgr/ ln_zco, ln_zps, ln_sco, ln_isfcav, ln_linssh 105 105 !!---------------------------------------------------------------------- 106 106 ! … … 120 120 WRITE(numout,*) 'dom_zgr : vertical coordinate' 121 121 WRITE(numout,*) '~~~~~~~' 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 122 WRITE(numout,*) ' Namelist namzgr : set vertical coordinate' 123 WRITE(numout,*) ' z-coordinate - full steps ln_zco = ', ln_zco 124 WRITE(numout,*) ' z-coordinate - partial steps ln_zps = ', ln_zps 125 WRITE(numout,*) ' s- or hybrid z-s-coordinate ln_sco = ', ln_sco 126 WRITE(numout,*) ' ice shelf cavities ln_isfcav = ', ln_isfcav 127 WRITE(numout,*) ' linear free surface ln_linssh = ', ln_linssh 127 128 ENDIF 129 130 IF( ln_linssh .AND. lwp) WRITE(numout,*) ' linear free surface: the vertical mesh does not change in time' 128 131 129 132 ioptio = 0 ! Check Vertical coordinate options … … 155 158 ! 156 159 IF( nprint == 1 .AND. lwp ) THEN 157 WRITE(numout,*) ' MIN val mbathy ', MINVAL(mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) )160 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 158 161 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 159 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gdep3w_0(:,:,:) )160 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL(e3f_0(:,:,:) ), &161 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL(e3v_0(:,:,:) ), &162 & ' uw', MINVAL( e3uw_0(:,:,:)), ' vw', MINVAL(e3vw_0(:,:,:)), &163 & ' w ', MINVAL(e3w_0(:,:,:) )162 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) ) 163 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), & 164 & ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), & 165 & ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), & 166 & ' w ', MINVAL( e3w_0(:,:,:) ) 164 167 165 168 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 166 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gdep3w_0(:,:,:) )167 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL(e3f_0(:,:,:) ), &168 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL(e3v_0(:,:,:) ), &169 & ' uw', MAXVAL( e3uw_0(:,:,:)), ' vw', MAXVAL( e3vw_0(:,:,:)),&170 & ' w ', MAXVAL(e3w_0(:,:,:) )169 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) ) 170 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), & 171 & ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), & 172 & ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), & 173 & ' w ', MAXVAL( e3w_0(:,:,:) ) 171 174 ENDIF 172 175 ! … … 674 677 !! - update bathy : meter bathymetry (in meters) 675 678 !!---------------------------------------------------------------------- 676 !!677 679 INTEGER :: ji, jj, jl ! dummy loop indices 678 680 INTEGER :: icompt, ibtest, ikmax ! temporary integers 679 681 REAL(wp), POINTER, DIMENSION(:,:) :: zbathy 680 681 682 !!---------------------------------------------------------------------- 682 683 ! … … 775 776 IF(lwp) WRITE(numout,*) ' you can decrease jpk to ', ikmax+1 776 777 ENDIF 777 778 IF( lwp .AND. nprint == 1 ) THEN ! control print779 WRITE(numout,*)780 WRITE(numout,*) ' bathymetric field : number of non-zero T-levels '781 WRITE(numout,*) ' ------------------'782 CALL prihin( mbathy, jpi, jpj, 1, jpi, 1, 1, jpj, 1, 3, numout )783 WRITE(numout,*)784 ENDIF785 778 ! 786 779 CALL wrk_dealloc( jpi, jpj, zbathy ) … … 803 796 !! (min value = 1 over land) 804 797 !!---------------------------------------------------------------------- 805 !!806 798 INTEGER :: ji, jj ! dummy loop indices 807 799 REAL(wp), POINTER, DIMENSION(:,:) :: zmbk … … 835 827 END SUBROUTINE zgr_bot_level 836 828 837 SUBROUTINE zgr_top_level 829 830 SUBROUTINE zgr_top_level 838 831 !!---------------------------------------------------------------------- 839 832 !! *** ROUTINE zgr_bot_level *** … … 847 840 !! (min value = 1) 848 841 !!---------------------------------------------------------------------- 849 !!850 842 INTEGER :: ji, jj ! dummy loop indices 851 843 REAL(wp), POINTER, DIMENSION(:,:) :: zmik … … 881 873 END SUBROUTINE zgr_top_level 882 874 875 883 876 SUBROUTINE zgr_zco 884 877 !!---------------------------------------------------------------------- 885 878 !! *** ROUTINE zgr_zco *** 886 879 !! 887 !! ** Purpose : define the z-coordinate system880 !! ** Purpose : define the reference z-coordinate system 888 881 !! 889 882 !! ** Method : set 3D coord. arrays to reference 1D array … … 895 888 ! 896 889 DO jk = 1, jpk 897 gdept_0 898 gdepw_0 899 gde p3w_0(:,:,jk) = gdepw_1d(jk)900 e3t_0 901 e3u_0 902 e3v_0 903 e3f_0 904 e3w_0 905 e3uw_0 906 e3vw_0 890 gdept_0(:,:,jk) = gdept_1d(jk) 891 gdepw_0(:,:,jk) = gdepw_1d(jk) 892 gde3w_0(:,:,jk) = gdepw_1d(jk) 893 e3t_0 (:,:,jk) = e3t_1d (jk) 894 e3u_0 (:,:,jk) = e3t_1d (jk) 895 e3v_0 (:,:,jk) = e3t_1d (jk) 896 e3f_0 (:,:,jk) = e3t_1d (jk) 897 e3w_0 (:,:,jk) = e3w_1d (jk) 898 e3uw_0 (:,:,jk) = e3w_1d (jk) 899 e3vw_0 (:,:,jk) = e3w_1d (jk) 907 900 END DO 908 901 ! … … 917 910 !! 918 911 !! ** Purpose : the depth and vertical scale factor in partial step 919 !! z-coordinate case912 !! reference z-coordinate case 920 913 !! 921 914 !! ** Method : Partial steps : computes the 3D vertical scale factors … … 957 950 !! Reference : Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 958 951 !!---------------------------------------------------------------------- 959 !!960 952 INTEGER :: ji, jj, jk ! dummy loop indices 961 953 INTEGER :: ik, it, ikb, ikt ! temporary integers 962 LOGICAL :: ll_print ! Allow control print for debugging963 954 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 964 955 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t … … 971 962 IF( nn_timing == 1 ) CALL timing_start('zgr_zps') 972 963 ! 973 CALL wrk_alloc( jpi, jpj, jpk,zprt )964 CALL wrk_alloc( jpi,jpj,jpk, zprt ) 974 965 ! 975 966 IF(lwp) WRITE(numout,*) … … 977 968 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 978 969 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 979 980 ll_print = .FALSE. ! Local variable for debugging981 982 IF(lwp .AND. ll_print) THEN ! control print of the ocean depth983 WRITE(numout,*)984 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)'985 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout )986 ENDIF987 988 970 989 971 ! bathymetry in level (from bathy_meter) … … 1196 1178 IF( MINVAL( gdepw_0(:,:,:) ) < 0._wp ) CALL ctl_stop( ' zgr_zps : e r r o r gdepw_0 < 0' ) 1197 1179 1198 ! Compute gde p3w_0 (vertical sum of e3w)1180 ! Compute gde3w_0 (vertical sum of e3w) 1199 1181 IF ( ln_isfcav ) THEN ! if cavity 1200 WHERE (misfdep == 0)misfdep = 11182 WHERE( misfdep == 0 ) misfdep = 1 1201 1183 DO jj = 1,jpj 1202 1184 DO ji = 1,jpi 1203 gde p3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1)1185 gde3w_0(ji,jj,1) = 0.5_wp * e3w_0(ji,jj,1) 1204 1186 DO jk = 2, misfdep(ji,jj) 1205 gde p3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1187 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1206 1188 END DO 1207 IF (misfdep(ji,jj) .GE. 2) gdep3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj))1189 IF( misfdep(ji,jj) >= 2 ) gde3w_0(ji,jj,misfdep(ji,jj)) = risfdep(ji,jj) + 0.5_wp * e3w_0(ji,jj,misfdep(ji,jj)) 1208 1190 DO jk = misfdep(ji,jj) + 1, jpk 1209 gde p3w_0(ji,jj,jk) = gdep3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk)1191 gde3w_0(ji,jj,jk) = gde3w_0(ji,jj,jk-1) + e3w_0(ji,jj,jk) 1210 1192 END DO 1211 1193 END DO 1212 1194 END DO 1213 1195 ELSE ! no cavity 1214 gde p3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)1196 gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1) 1215 1197 DO jk = 2, jpk 1216 gde p3w_0(:,:,jk) = gdep3w_0(:,:,jk-1) + e3w_0(:,:,jk)1198 gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk) 1217 1199 END DO 1218 1200 END IF 1219 ! ! ================= ! 1220 IF(lwp .AND. ll_print) THEN ! Control print ! 1221 ! ! ================= ! 1222 DO jj = 1,jpj 1223 DO ji = 1, jpi 1224 ik = MAX( mbathy(ji,jj), 1 ) 1225 zprt(ji,jj,1) = e3t_0 (ji,jj,ik) 1226 zprt(ji,jj,2) = e3w_0 (ji,jj,ik) 1227 zprt(ji,jj,3) = e3u_0 (ji,jj,ik) 1228 zprt(ji,jj,4) = e3v_0 (ji,jj,ik) 1229 zprt(ji,jj,5) = e3f_0 (ji,jj,ik) 1230 zprt(ji,jj,6) = gdep3w_0(ji,jj,ik) 1231 END DO 1232 END DO 1233 WRITE(numout,*) 1234 WRITE(numout,*) 'domzgr e3t(mbathy)' ; CALL prihre(zprt(:,:,1),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1235 WRITE(numout,*) 1236 WRITE(numout,*) 'domzgr e3w(mbathy)' ; CALL prihre(zprt(:,:,2),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1237 WRITE(numout,*) 1238 WRITE(numout,*) 'domzgr e3u(mbathy)' ; CALL prihre(zprt(:,:,3),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1239 WRITE(numout,*) 1240 WRITE(numout,*) 'domzgr e3v(mbathy)' ; CALL prihre(zprt(:,:,4),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1241 WRITE(numout,*) 1242 WRITE(numout,*) 'domzgr e3f(mbathy)' ; CALL prihre(zprt(:,:,5),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1243 WRITE(numout,*) 1244 WRITE(numout,*) 'domzgr gdep3w(mbathy)' ; CALL prihre(zprt(:,:,6),jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1245 ENDIF 1246 ! 1247 CALL wrk_dealloc( jpi, jpj, jpk, zprt ) 1201 ! 1202 CALL wrk_dealloc( jpi,jpj,jpk, zprt ) 1248 1203 ! 1249 1204 IF( nn_timing == 1 ) CALL timing_stop('zgr_zps') 1250 1205 ! 1251 1206 END SUBROUTINE zgr_zps 1207 1252 1208 1253 1209 SUBROUTINE zgr_isf … … 1265 1221 !! - bathy and isfdraft are modified 1266 1222 !!---------------------------------------------------------------------- 1267 !!1268 1223 INTEGER :: ji, jj, jk, jl ! dummy loop indices 1269 1224 INTEGER :: ik, it ! temporary integers 1270 1225 INTEGER :: id, jd, nprocd 1271 1226 INTEGER :: icompt, ibtest, ibtestim1, ibtestip1, ibtestjm1, ibtestjp1 ! (ISF) 1272 LOGICAL :: ll_print ! Allow control print for debugging1273 1227 REAL(wp) :: ze3tp , ze3wp ! Last ocean level thickness at T- and W-points 1274 1228 REAL(wp) :: zdepwp, zdepth ! Ajusted ocean depth to avoid too small e3t … … 1281 1235 !!--------------------------------------------------------------------- 1282 1236 ! 1283 IF( nn_timing == 1 ) CALL timing_start('zgr_isf')1284 ! 1285 CALL wrk_alloc( jpi, jpj,zbathy, zmask, zrisfdep)1286 CALL wrk_alloc( jpi, jpj,zmisfdep, zmbathy )1237 IF( nn_timing == 1 ) CALL timing_start('zgr_isf') 1238 ! 1239 CALL wrk_alloc( jpi,jpj, zbathy, zmask, zrisfdep) 1240 CALL wrk_alloc( jpi,jpj, zmisfdep, zmbathy ) 1287 1241 1288 1242 … … 1300 1254 WHERE( 0._wp < risfdep(:,:) .AND. risfdep(:,:) >= zdepth ) misfdep(:,:) = jk+1 1301 1255 END DO 1302 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) .GT.0._wp)1303 risfdep(:,:) = 0. ;misfdep(:,:) = 11256 WHERE (risfdep(:,:) <= e3t_1d(1) .AND. risfdep(:,:) > 0._wp) 1257 risfdep(:,:) = 0. ; misfdep(:,:) = 1 1304 1258 END WHERE 1305 1259 … … 1308 1262 ! run the bathy check 10 times to be sure all the modif in the bathy or iceshelf draft are compatible together 1309 1263 DO jl = 1, 10 1310 WHERE (bathy(:,:) .EQ.risfdep(:,:) )1264 WHERE (bathy(:,:) == risfdep(:,:) ) 1311 1265 misfdep(:,:) = 0 ; risfdep(:,:) = 0._wp 1312 1266 mbathy (:,:) = 0 ; bathy (:,:) = 0._wp … … 1315 1269 misfdep(:,:) = 0; risfdep(:,:) = 0._wp 1316 1270 mbathy (:,:) = 0; bathy (:,:) = 0._wp 1317 END WHERE1271 END WHERE 1318 1272 IF( lk_mpp ) THEN 1319 1273 zbathy(:,:) = FLOAT( misfdep(:,:) ) … … 1360 1314 ! find the minimum change option: 1361 1315 ! test bathy 1362 IF (risfdep(ji,jj) .GT.1) THEN1316 IF (risfdep(ji,jj) > 1) THEN 1363 1317 zbathydiff =ABS(bathy(ji,jj) - (gdepw_1d(mbathy (ji,jj)+1) & 1364 1318 & + MIN( e3zps_min, e3t_1d(mbathy (ji,jj)+1)*e3zps_rat ))) … … 1366 1320 & - MIN( e3zps_min, e3t_1d(misfdep(ji,jj)-1)*e3zps_rat ))) 1367 1321 1368 IF (bathy(ji,jj) .GT. risfdep(ji,jj) .AND. mbathy(ji,jj) .LT.misfdep(ji,jj)) THEN1322 IF (bathy(ji,jj) > risfdep(ji,jj) .AND. mbathy(ji,jj) < misfdep(ji,jj)) THEN 1369 1323 IF (zbathydiff .LE. zrisfdepdiff) THEN 1370 1324 bathy(ji,jj) = gdepw_1d(mbathy(ji,jj)) + MIN( e3zps_min, e3t_1d(mbathy(ji,jj)+1)*e3zps_rat ) … … 1752 1706 CALL wrk_dealloc( jpi, jpj, zmask, zbathy, zrisfdep ) 1753 1707 CALL wrk_dealloc( jpi, jpj, zmisfdep, zmbathy ) 1754 1755 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf')1756 1708 ! 1709 IF( nn_timing == 1 ) CALL timing_stop('zgr_isf') 1710 ! 1757 1711 END SUBROUTINE 1712 1758 1713 1759 1714 SUBROUTINE zgr_sco … … 1801 1756 !! 1802 1757 !!---------------------------------------------------------------------- 1803 !1804 1758 INTEGER :: ji, jj, jk, jl ! dummy loop argument 1805 1759 INTEGER :: iip1, ijp1, iim1, ijm1 ! temporary integers … … 1810 1764 REAL(wp), POINTER, DIMENSION(:,: ) :: ztmpi1, ztmpi2, ztmpj1, ztmpj2 1811 1765 REAL(wp), POINTER, DIMENSION(:,: ) :: zenv, ztmp, zmsk, zri, zrj, zhbat 1812 1766 !! 1813 1767 NAMELIST/namzgr_sco/ln_s_sh94, ln_s_sf12, ln_sigcrit, rn_sbot_min, rn_sbot_max, rn_hc, rn_rmax,rn_theta, & 1814 1768 & rn_thetb, rn_bb, rn_alpha, rn_efold, rn_zs, rn_zb_a, rn_zb_b 1815 1769 !!---------------------------------------------------------------------- 1816 1770 ! 1817 1771 IF( nn_timing == 1 ) CALL timing_start('zgr_sco') 1818 1772 ! 1819 CALL wrk_alloc( jpi, jpj,zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 )1773 CALL wrk_alloc( jpi,jpj, zenv, ztmp, zmsk, zri, zrj, zhbat , ztmpi1, ztmpi2, ztmpj1, ztmpj2 ) 1820 1774 ! 1821 1775 REWIND( numnam_ref ) ! Namelist namzgr_sco in reference namelist : Sigma-stretching parameters … … 1876 1830 DO jj = 1, jpj 1877 1831 DO ji = 1, jpi 1878 IF( bathy(ji,jj) == 0._wp ) THEN 1879 iip1 = MIN( ji+1, jpi ) 1880 ijp1 = MIN( jj+1, jpj ) 1881 iim1 = MAX( ji-1, 1 ) 1882 ijm1 = MAX( jj-1, 1 ) 1883 IF( (bathy(iip1,jj) + bathy(iim1,jj) + bathy(ji,ijp1) + bathy(ji,ijm1) + & 1884 & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1885 zenv(ji,jj) = rn_sbot_min 1886 ENDIF 1832 IF( bathy(ji,jj) == 0._wp ) THEN 1833 iip1 = MIN( ji+1, jpi ) 1834 ijp1 = MIN( jj+1, jpj ) 1835 iim1 = MAX( ji-1, 1 ) 1836 ijm1 = MAX( jj-1, 1 ) 1837 !!gm BUG fix see ticket #1617 1838 IF( ( + bathy(iim1,ijm1) + bathy(ji,ijp1) + bathy(iip1,ijp1) & 1839 & + bathy(iim1,jj ) + bathy(iip1,jj ) & 1840 & + bathy(iim1,ijm1) + bathy(ji,ijm1) + bathy(iip1,ijp1) ) > 0._wp ) zenv(ji,jj) = rn_sbot_min 1841 !!gm 1842 !!gm IF( ( bathy(iip1,jj ) + bathy(iim1,jj ) + bathy(ji,ijp1 ) + bathy(ji,ijm1) + & 1843 !!gm & bathy(iip1,ijp1) + bathy(iim1,ijm1) + bathy(iip1,ijp1) + bathy(iim1,ijm1)) > 0._wp ) THEN 1844 !!gm zenv(ji,jj) = rn_sbot_min 1845 !!gm ENDIF 1846 !!gm end 1887 1847 ENDIF 1888 1848 END DO … … 1975 1935 ENDIF 1976 1936 ! 1977 IF(lwp) THEN ! Control print1978 WRITE(numout,*)1979 WRITE(numout,*) ' domzgr: hbatt field; ocean depth in meters'1980 WRITE(numout,*)1981 CALL prihre( hbatt(1,1), jpi, jpj, 1, jpi, 1, 1, jpj, 1, 0._wp, numout )1982 IF( nprint == 1 ) THEN1983 WRITE(numout,*) ' bathy MAX ', MAXVAL( bathy(:,:) ), ' MIN ', MINVAL( bathy(:,:) )1984 WRITE(numout,*) ' hbatt MAX ', MAXVAL( hbatt(:,:) ), ' MIN ', MINVAL( hbatt(:,:) )1985 ENDIF1986 ENDIF1987 1988 1937 ! ! ============================== 1989 1938 ! ! hbatu, hbatv, hbatf fields … … 2080 2029 CALL lbc_lnk( e3uw_0, 'U', 1._wp ) 2081 2030 CALL lbc_lnk( e3vw_0, 'V', 1._wp ) 2082 2083 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2084 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2085 ! 2086 where (e3t_0 (:,:,:).eq.0.0) e3t_0(:,:,:) = 1.0 2087 where (e3u_0 (:,:,:).eq.0.0) e3u_0(:,:,:) = 1.0 2088 where (e3v_0 (:,:,:).eq.0.0) e3v_0(:,:,:) = 1.0 2089 where (e3f_0 (:,:,:).eq.0.0) e3f_0(:,:,:) = 1.0 2090 where (e3w_0 (:,:,:).eq.0.0) e3w_0(:,:,:) = 1.0 2091 where (e3uw_0 (:,:,:).eq.0.0) e3uw_0(:,:,:) = 1.0 2092 where (e3vw_0 (:,:,:).eq.0.0) e3vw_0(:,:,:) = 1.0 2031 ! 2032 WHERE( e3t_0 (:,:,:) == 0._wp ) e3t_0 (:,:,:) = 1._wp 2033 WHERE( e3u_0 (:,:,:) == 0._wp ) e3u_0 (:,:,:) = 1._wp 2034 WHERE( e3v_0 (:,:,:) == 0._wp ) e3v_0 (:,:,:) = 1._wp 2035 WHERE( e3f_0 (:,:,:) == 0._wp ) e3f_0 (:,:,:) = 1._wp 2036 WHERE( e3w_0 (:,:,:) == 0._wp ) e3w_0 (:,:,:) = 1._wp 2037 WHERE( e3uw_0(:,:,:) == 0._wp ) e3uw_0(:,:,:) = 1._wp 2038 WHERE( e3vw_0(:,:,:) == 0._wp ) e3vw_0(:,:,:) = 1._wp 2093 2039 2094 2040 #if defined key_agrif 2095 ! Ensure meaningful vertical scale factors in ghost lines/columns 2096 IF( .NOT. Agrif_Root() ) THEN 2097 ! 2098 IF((nbondi == -1).OR.(nbondi == 2)) THEN 2099 e3u_0(1,:,:) = e3u_0(2,:,:) 2100 ENDIF 2101 ! 2102 IF((nbondi == 1).OR.(nbondi == 2)) THEN 2103 e3u_0(nlci-1,:,:) = e3u_0(nlci-2,:,:) 2104 ENDIF 2105 ! 2106 IF((nbondj == -1).OR.(nbondj == 2)) THEN 2107 e3v_0(:,1,:) = e3v_0(:,2,:) 2108 ENDIF 2109 ! 2110 IF((nbondj == 1).OR.(nbondj == 2)) THEN 2111 e3v_0(:,nlcj-1,:) = e3v_0(:,nlcj-2,:) 2112 ENDIF 2113 ! 2114 ENDIF 2041 IF( .NOT. Agrif_Root() ) THEN ! Ensure meaningful vertical scale factors in ghost lines/columns 2042 IF( nbondi == -1 .OR. nbondi == 2 ) e3u_0( 1 , : ,:) = e3u_0( 2 , : ,:) 2043 IF( nbondi == 1 .OR. nbondi == 2 ) e3u_0(nlci-1, : ,:) = e3u_0(nlci-2, : ,:) 2044 IF( nbondj == -1 .OR. nbondj == 2 ) e3v_0( : , 1 ,:) = e3v_0( : , 2 ,:) 2045 IF( nbondj == 1 .OR. nbondj == 2 ) e3v_0( : ,nlcj-1,:) = e3v_0( : ,nlcj-2,:) 2046 ENDIF 2115 2047 #endif 2116 2048 2117 fsdept(:,:,:) = gdept_0 (:,:,:) 2118 fsdepw(:,:,:) = gdepw_0 (:,:,:) 2119 fsde3w(:,:,:) = gdep3w_0(:,:,:) 2120 fse3t (:,:,:) = e3t_0 (:,:,:) 2121 fse3u (:,:,:) = e3u_0 (:,:,:) 2122 fse3v (:,:,:) = e3v_0 (:,:,:) 2123 fse3f (:,:,:) = e3f_0 (:,:,:) 2124 fse3w (:,:,:) = e3w_0 (:,:,:) 2125 fse3uw(:,:,:) = e3uw_0 (:,:,:) 2126 fse3vw(:,:,:) = e3vw_0 (:,:,:) 2049 !!gm I don't like that HERE we are supposed to set the reference coordinate (i.e. _0 arrays) 2050 !!gm and only that !!!!! 2051 !!gm THIS should be removed from here ! 2052 gdept_n(:,:,:) = gdept_0(:,:,:) 2053 gdepw_n(:,:,:) = gdepw_0(:,:,:) 2054 gde3w_n(:,:,:) = gde3w_0(:,:,:) 2055 e3t_n (:,:,:) = e3t_0 (:,:,:) 2056 e3u_n (:,:,:) = e3u_0 (:,:,:) 2057 e3v_n (:,:,:) = e3v_0 (:,:,:) 2058 e3f_n (:,:,:) = e3f_0 (:,:,:) 2059 e3w_n (:,:,:) = e3w_0 (:,:,:) 2060 e3uw_n (:,:,:) = e3uw_0 (:,:,:) 2061 e3vw_n (:,:,:) = e3vw_0 (:,:,:) 2062 !!gm and obviously in the following, use the _0 arrays until the end of this subroutine 2063 !! gm end 2127 2064 !! 2128 2065 ! HYBRID : … … 2130 2067 DO ji = 1, jpi 2131 2068 DO jk = 1, jpkm1 2132 IF( scobot(ji,jj) >= fsdept(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk )2133 END DO 2134 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 02069 IF( scobot(ji,jj) >= gdept_n(ji,jj,jk) ) mbathy(ji,jj) = MAX( 2, jk ) 2070 END DO 2071 IF( scobot(ji,jj) == 0._wp ) mbathy(ji,jj) = 0 2135 2072 END DO 2136 2073 END DO … … 2141 2078 WRITE(numout,*) ' MIN val mbathy ', MINVAL( mbathy(:,:) ), ' MAX ', MAXVAL( mbathy(:,:) ) 2142 2079 WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), & 2143 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde p3w_0(:,:,:) )2144 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 2145 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 2146 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 2080 & ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ' , MINVAL( gde3w_0(:,:,:) ) 2081 WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0 (:,:,:) ), ' f ' , MINVAL( e3f_0 (:,:,:) ), & 2082 & ' u ', MINVAL( e3u_0 (:,:,:) ), ' u ' , MINVAL( e3v_0 (:,:,:) ), & 2083 & ' uw', MINVAL( e3uw_0 (:,:,:) ), ' vw' , MINVAL( e3vw_0 (:,:,:) ), & 2147 2084 & ' w ', MINVAL( e3w_0 (:,:,:) ) 2148 2085 2149 2086 WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), & 2150 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde p3w_0(:,:,:) )2151 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 2152 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 2153 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 2087 & ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ' , MAXVAL( gde3w_0(:,:,:) ) 2088 WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0 (:,:,:) ), ' f ' , MAXVAL( e3f_0 (:,:,:) ), & 2089 & ' u ', MAXVAL( e3u_0 (:,:,:) ), ' u ' , MAXVAL( e3v_0 (:,:,:) ), & 2090 & ' uw', MAXVAL( e3uw_0 (:,:,:) ), ' vw' , MAXVAL( e3vw_0 (:,:,:) ), & 2154 2091 & ' w ', MAXVAL( e3w_0 (:,:,:) ) 2155 2092 ENDIF … … 2183 2120 END DO 2184 2121 ENDIF 2185 2186 !================================================================================2187 ! check the coordinate makes sense2188 !================================================================================2122 ! 2123 !================================================================================ 2124 ! check the coordinate makes sense 2125 !================================================================================ 2189 2126 DO ji = 1, jpi 2190 2127 DO jj = 1, jpj 2191 2128 ! 2192 2129 IF( hbatt(ji,jj) > 0._wp) THEN 2193 2130 DO jk = 1, mbathy(ji,jj) 2194 2131 ! check coordinate is monotonically increasing 2195 IF ( fse3w(ji,jj,jk) <= 0._wp .OR. fse3t(ji,jj,jk) <= 0._wp ) THEN2132 IF (e3w_n(ji,jj,jk) <= 0._wp .OR. e3t_n(ji,jj,jk) <= 0._wp ) THEN 2196 2133 WRITE(ctmp1,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2197 2134 WRITE(numout,*) 'ERROR zgr_sco : e3w or e3t =< 0 at point (i,j,k)= ', ji, jj, jk 2198 WRITE(numout,*) 'e3w', fse3w(ji,jj,:)2199 WRITE(numout,*) 'e3t', fse3t(ji,jj,:)2135 WRITE(numout,*) 'e3w',e3w_n(ji,jj,:) 2136 WRITE(numout,*) 'e3t',e3t_n(ji,jj,:) 2200 2137 CALL ctl_stop( ctmp1 ) 2201 2138 ENDIF 2202 2139 ! and check it has never gone negative 2203 IF( fsdepw(ji,jj,jk) < 0._wp .OR. fsdept(ji,jj,jk) < 0._wp ) THEN2140 IF( gdepw_n(ji,jj,jk) < 0._wp .OR. gdept_n(ji,jj,jk) < 0._wp ) THEN 2204 2141 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2205 2142 WRITE(numout,*) 'ERROR zgr_sco : gdepw or gdept =< 0 at point (i,j,k)= ', ji, jj, jk 2206 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2207 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2143 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2144 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2208 2145 CALL ctl_stop( ctmp1 ) 2209 2146 ENDIF 2210 2147 ! and check it never exceeds the total depth 2211 IF( fsdepw(ji,jj,jk) > hbatt(ji,jj) ) THEN2148 IF( gdepw_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2212 2149 WRITE(ctmp1,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2213 2150 WRITE(numout,*) 'ERROR zgr_sco : gdepw > hbatt at point (i,j,k)= ', ji, jj, jk 2214 WRITE(numout,*) 'gdepw', fsdepw(ji,jj,:)2151 WRITE(numout,*) 'gdepw',gdepw_n(ji,jj,:) 2215 2152 CALL ctl_stop( ctmp1 ) 2216 2153 ENDIF 2217 2154 END DO 2218 2155 ! 2219 2156 DO jk = 1, mbathy(ji,jj)-1 2220 2157 ! and check it never exceeds the total depth 2221 IF( fsdept(ji,jj,jk) > hbatt(ji,jj) ) THEN2158 IF( gdept_n(ji,jj,jk) > hbatt(ji,jj) ) THEN 2222 2159 WRITE(ctmp1,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2223 2160 WRITE(numout,*) 'ERROR zgr_sco : gdept > hbatt at point (i,j,k)= ', ji, jj, jk 2224 WRITE(numout,*) 'gdept', fsdept(ji,jj,:)2161 WRITE(numout,*) 'gdept',gdept_n(ji,jj,:) 2225 2162 CALL ctl_stop( ctmp1 ) 2226 2163 ENDIF 2227 2164 END DO 2228 2229 2165 ENDIF 2230 2231 2166 END DO 2232 2167 END DO … … 2238 2173 END SUBROUTINE zgr_sco 2239 2174 2240 !!====================================================================== 2175 2241 2176 SUBROUTINE s_sh94() 2242 2243 2177 !!---------------------------------------------------------------------- 2244 2178 !! *** ROUTINE s_sh94 *** … … 2251 2185 !! Reference : Song and Haidvogel 1994. 2252 2186 !!---------------------------------------------------------------------- 2253 !2254 2187 INTEGER :: ji, jj, jk ! dummy loop argument 2255 2188 REAL(wp) :: zcoeft, zcoefw ! temporary scalars … … 2257 2190 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2258 2191 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2259 2260 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2261 CALL wrk_alloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2192 !!---------------------------------------------------------------------- 2193 2194 CALL wrk_alloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2195 CALL wrk_alloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2262 2196 2263 2197 z_gsigw3 = 0._wp ; z_gsigt3 = 0._wp ; z_gsi3w3 = 0._wp … … 2265 2199 z_esigtu3 = 0._wp ; z_esigtv3 = 0._wp ; z_esigtf3 = 0._wp 2266 2200 z_esigwu3 = 0._wp ; z_esigwv3 = 0._wp 2267 2201 ! 2268 2202 DO ji = 1, jpi 2269 2203 DO jj = 1, jpj 2270 2204 ! 2271 2205 IF( hbatt(ji,jj) > rn_hc ) THEN !deep water, stretched sigma 2272 2206 DO jk = 1, jpk … … 2297 2231 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2298 2232 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2299 gdept_0 2300 gdepw_0 2301 gde p3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft )2233 gdept_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigt3(ji,jj,jk)+rn_hc*zcoeft ) 2234 gdepw_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsigw3(ji,jj,jk)+rn_hc*zcoefw ) 2235 gde3w_0(ji,jj,jk) = ( scosrf(ji,jj) + (hbatt(ji,jj)-rn_hc)*z_gsi3w3(ji,jj,jk)+rn_hc*zcoeft ) 2302 2236 END DO 2303 2237 ! … … 2331 2265 END DO 2332 2266 END DO 2333 2334 CALL wrk_dealloc( jpi, jpj, jpk,z_gsigw3, z_gsigt3, z_gsi3w3 )2335 CALL wrk_dealloc( jpi, jpj, jpk,z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 )2336 2267 ! 2268 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2269 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2270 ! 2337 2271 END SUBROUTINE s_sh94 2338 2272 2273 2339 2274 SUBROUTINE s_sf12 2340 2341 2275 !!---------------------------------------------------------------------- 2342 2276 !! *** ROUTINE s_sf12 *** … … 2354 2288 !! Reference : Siddorn and Furner 2012 (submitted Ocean modelling). 2355 2289 !!---------------------------------------------------------------------- 2356 !2357 2290 INTEGER :: ji, jj, jk ! dummy loop argument 2358 2291 REAL(wp) :: zsmth ! smoothing around critical depth … … 2361 2294 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_gsigw3, z_gsigt3, z_gsi3w3 2362 2295 REAL(wp), POINTER, DIMENSION(:,:,:) :: z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 2363 2296 !!---------------------------------------------------------------------- 2364 2297 ! 2365 2298 CALL wrk_alloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) … … 2425 2358 2426 2359 DO jk = 1, jpk 2427 gdept_0 2428 gdepw_0 2429 gde p3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk)2360 gdept_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigt3(ji,jj,jk) 2361 gdepw_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsigw3(ji,jj,jk) 2362 gde3w_0(ji,jj,jk) = (scosrf(ji,jj)+hbatt(ji,jj))*z_gsi3w3(ji,jj,jk) 2430 2363 END DO 2431 2364 … … 2454 2387 e3f_0(ji,jj,jk)=(scosrf(ji,jj)+hbatf(ji,jj))*z_esigtf3(ji,jj,jk) 2455 2388 ! 2456 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk)2389 e3w_0 (ji,jj,jk)=hbatt(ji,jj)*z_esigw3(ji,jj,jk) 2457 2390 e3uw_0(ji,jj,jk)=hbatu(ji,jj)*z_esigwu3(ji,jj,jk) 2458 2391 e3vw_0(ji,jj,jk)=hbatv(ji,jj)*z_esigwv3(ji,jj,jk) … … 2467 2400 CALL lbc_lnk(e3uw_0,'T',1.) ; CALL lbc_lnk(e3vw_0,'T',1.) 2468 2401 ! 2469 ! ! ============= 2470 2471 CALL wrk_dealloc( jpi, jpj, jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2472 CALL wrk_dealloc( jpi, jpj, jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2473 2402 CALL wrk_dealloc( jpi,jpj,jpk, z_gsigw3, z_gsigt3, z_gsi3w3 ) 2403 CALL wrk_dealloc( jpi,jpj,jpk, z_esigt3, z_esigw3, z_esigtu3, z_esigtv3, z_esigtf3, z_esigwu3, z_esigwv3 ) 2404 ! 2474 2405 END SUBROUTINE s_sf12 2475 2406 2407 2476 2408 SUBROUTINE s_tanh() 2477 2478 2409 !!---------------------------------------------------------------------- 2479 2410 !! *** ROUTINE s_tanh*** … … 2485 2416 !! Reference : Madec, Lott, Delecluse and Crepon, 1996. JPO, 26, 1393-1408. 2486 2417 !!---------------------------------------------------------------------- 2487 2488 INTEGER :: ji, jj, jk ! dummy loop argument 2418 INTEGER :: ji, jj, jk ! dummy loop argument 2489 2419 REAL(wp) :: zcoeft, zcoefw ! temporary scalars 2490 2491 2420 REAL(wp), POINTER, DIMENSION(:) :: z_gsigw, z_gsigt, z_gsi3w 2492 2421 REAL(wp), POINTER, DIMENSION(:) :: z_esigt, z_esigw 2493 2494 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2495 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2422 !!---------------------------------------------------------------------- 2423 2424 CALL wrk_alloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2425 CALL wrk_alloc( jpk, z_esigt, z_esigw ) 2496 2426 2497 2427 z_gsigw = 0._wp ; z_gsigt = 0._wp ; z_gsi3w = 0._wp … … 2523 2453 zcoeft = ( REAL(jk,wp) - 0.5_wp ) / REAL(jpkm1,wp) 2524 2454 zcoefw = ( REAL(jk,wp) - 1.0_wp ) / REAL(jpkm1,wp) 2525 gdept_0 2526 gdepw_0 2527 gde p3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft )2455 gdept_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigt(jk) + hift(:,:)*zcoeft ) 2456 gdepw_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsigw(jk) + hift(:,:)*zcoefw ) 2457 gde3w_0(:,:,jk) = ( scosrf(:,:) + (hbatt(:,:)-hift(:,:))*z_gsi3w(jk) + hift(:,:)*zcoeft ) 2528 2458 END DO 2529 2459 !!gm: e3uw, e3vw can be suppressed (modif in dynzdf, dynzdf_iso, zdfbfr) (save 2 3D arrays) … … 2542 2472 END DO 2543 2473 END DO 2544 2545 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w)2546 CALL wrk_dealloc( jpk, z_esigt, z_esigw)2547 2474 ! 2475 CALL wrk_dealloc( jpk, z_gsigw, z_gsigt, z_gsi3w ) 2476 CALL wrk_dealloc( jpk, z_esigt, z_esigw ) 2477 ! 2548 2478 END SUBROUTINE s_tanh 2479 2549 2480 2550 2481 FUNCTION fssig( pk ) RESULT( pf ) … … 2618 2549 REAL(wp), INTENT(in ) :: pk1(jpk) ! continuous "k" coordinate 2619 2550 REAL(wp) :: p_gamma(jpk) ! stretched coordinate 2620 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth2621 REAL(wp), INTENT(in ) :: pzs ! surface box depth2622 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter2623 REAL(wp) :: za1,za2,za3 ! local variables2624 REAL(wp) :: zn1,zn2 ! local variables2625 REAL(wp) :: za,zb,zx ! local variables2626 integer :: jk2627 !!----------------------------------------------------------------------2628 ! 2629 2630 zn1 = 1. /(jpk-1.)2631 zn2 = 1. - zn12632 2551 REAL(wp), INTENT(in ) :: pzb ! Bottom box depth 2552 REAL(wp), INTENT(in ) :: pzs ! surface box depth 2553 REAL(wp), INTENT(in ) :: psmth ! Smoothing parameter 2554 ! 2555 INTEGER :: jk ! dummy loop index 2556 REAL(wp) :: za1,za2,za3 ! local scalar 2557 REAL(wp) :: zn1,zn2 ! - - 2558 REAL(wp) :: za,zb,zx ! - - 2559 !!---------------------------------------------------------------------- 2560 ! 2561 zn1 = 1._wp / REAL( jpkm1, wp ) 2562 zn2 = 1._wp - zn1 2563 ! 2633 2564 za1 = (rn_alpha+2.0_wp)*zn1**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn1**(rn_alpha+2.0_wp) 2634 2565 za2 = (rn_alpha+2.0_wp)*zn2**(rn_alpha+1.0_wp)-(rn_alpha+1.0_wp)*zn2**(rn_alpha+2.0_wp) 2635 2566 za3 = (zn2**3.0_wp - za2)/( zn1**3.0_wp - za1) 2636 2567 ! 2637 2568 za = pzb - za3*(pzs-za1)-za2 2638 2569 za = za/( zn2-0.5_wp*(za2+zn2**2.0_wp) - za3*(zn1-0.5_wp*(za1+zn1**2.0_wp) ) ) 2639 2570 zb = (pzs - za1 - za*( zn1-0.5_wp*(za1+zn1**2.0_wp ) ) ) / (zn1**3.0_wp - za1) 2640 2571 zx = 1.0_wp-za/2.0_wp-zb 2641 2572 ! 2642 2573 DO jk = 1, jpk 2643 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + &2644 &zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- &2645 &(rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) )2574 p_gamma(jk) = za*(pk1(jk)*(1.0_wp-pk1(jk)/2.0_wp))+zb*pk1(jk)**3.0_wp + & 2575 & zx*( (rn_alpha+2.0_wp)*pk1(jk)**(rn_alpha+1.0_wp)- & 2576 & (rn_alpha+1.0_wp)*pk1(jk)**(rn_alpha+2.0_wp) ) 2646 2577 p_gamma(jk) = p_gamma(jk)*psmth+pk1(jk)*(1.0_wp-psmth) 2647 ENDDO 2648 2578 END DO 2649 2579 ! 2650 2580 END FUNCTION fgamma
Note: See TracChangeset
for help on using the changeset viewer.
