Changeset 1083
- Timestamp:
- 2008-06-05T16:26:59+02:00 (16 years ago)
- Location:
- trunk/NEMO/OPA_SRC/DOM
- Files:
-
- 1 deleted
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DOM/domzgr.F90
r900 r1083 14 14 !! zgr_zps : z-coordinate with partial steps 15 15 !! zgr_sco : s-coordinate 16 !! fssig : sigma coordinate non-dimensional function 16 17 !!--------------------------------------------------------------------- 17 18 !! * Modules used … … 729 730 END SUBROUTINE zgr_zco 730 731 731 732 733 # include "domzgr_zps.h90" 734 735 736 #if ! defined key_zco 732 #if defined key_zco 737 733 !!---------------------------------------------------------------------- 738 !! NOT 'key_zco' : 3D gdep. and e3.734 !! 'key_zco' : "pure" zco (gdep & e3 are 1D arrays) 739 735 !!---------------------------------------------------------------------- 736 SUBROUTINE zgr_zps ! Empty routine 737 END SUBROUTINE zgr_zps 738 SUBROUTINE zgr_sco ! Empty routine 739 END SUBROUTINE zgr_sco 740 741 #else 742 !!---------------------------------------------------------------------- 743 !! Default option : zco, zps and/or sco available (gedp & e3 are 3D arrays) 744 !!---------------------------------------------------------------------- 745 746 SUBROUTINE zgr_zps 747 !!---------------------------------------------------------------------- 748 !! *** ROUTINE zgr_zps *** 749 !! 750 !! ** Purpose : the depth and vertical scale factor in partial step 751 !! z-coordinate case 752 !! 753 !! ** Method : Partial steps : computes the 3D vertical scale factors 754 !! of T-, U-, V-, W-, UW-, VW and F-points that are associated with 755 !! a partial step representation of bottom topography. 756 !! 757 !! The reference depth of model levels is defined from an analytical 758 !! function the derivative of which gives the reference vertical 759 !! scale factors. 760 !! From depth and scale factors reference, we compute there new value 761 !! with partial steps on 3d arrays ( i, j, k ). 762 !! 763 !! w-level: gdepw(i,j,k) = fsdep(k) 764 !! e3w(i,j,k) = dk(fsdep)(k) = fse3(i,j,k) 765 !! t-level: gdept(i,j,k) = fsdep(k+0.5) 766 !! e3t(i,j,k) = dk(fsdep)(k+0.5) = fse3(i,j,k+0.5) 767 !! 768 !! With the help of the bathymetric file ( bathymetry_depth_ORCA_R2.nc), 769 !! we find the mbathy index of the depth at each grid point. 770 !! This leads us to three cases: 771 !! 772 !! - bathy = 0 => mbathy = 0 773 !! - 1 < mbathy < jpkm1 774 !! - bathy > gdepw(jpk) => mbathy = jpkm1 775 !! 776 !! Then, for each case, we find the new depth at t- and w- levels 777 !! and the new vertical scale factors at t-, u-, v-, w-, uw-, vw- 778 !! and f-points. 779 !! 780 !! This routine is given as an example, it must be modified 781 !! following the user s desiderata. nevertheless, the output as 782 !! well as the way to compute the model levels and scale factors 783 !! must be respected in order to insure second order accuracy 784 !! schemes. 785 !! 786 !! c a u t i o n : gdept_0, gdepw_0 and e3._0 are positives 787 !! - - - - - - - gdept, gdepw and e3. are positives 788 !! 789 !! Reference : 790 !! Pacanowsky & Gnanadesikan 1997, Mon. Wea. Rev., 126, 3248-3270. 791 !! 792 !! History : 793 !! 8.5 ! 02-09 (A. Bozec, G. Madec) F90: Free form and module 794 !! 9.0 ! 02-09 (A. de Miranda) rigid-lid + islands 795 !!---------------------------------------------------------------------- 796 !! * Local declarations 797 INTEGER :: ji, jj, jk ! dummy loop indices 798 INTEGER :: & 799 ik, it ! temporary integers 800 801 REAL(wp) :: & 802 ze3tp, ze3wp, & ! Last ocean level thickness at T- and W-points 803 zdepwp, & ! Ajusted ocean depth to avoid too small e3t 804 zdepth, & ! " " 805 zmax, zmin, & ! Maximum and minimum depth 806 zdiff ! temporary scalar 807 808 REAL(wp), DIMENSION(jpi,jpj) :: & 809 zprt ! " " 810 811 LOGICAL :: ll_print ! Allow control print for debugging 812 813 !!--------------------------------------------------------------------- 814 !! OPA8.5, LODYC-IPSL (2002) 815 !!--------------------------------------------------------------------- 816 817 ! Local variable for debugging 818 ll_print=.FALSE. 819 !!! ll_print=.TRUE. 820 821 ! Initialization of constant 822 zmax = gdepw_0(jpk) + e3t_0(jpk) 823 zmin = gdepw_0(4) 824 825 ! Ocean depth 826 IF(lwp .AND. ll_print) THEN 827 WRITE(numout,*) 828 WRITE(numout,*) 'dom_zgr_zps: bathy (in hundred of meters)' 829 CALL prihre( bathy, jpi, jpj, 1,jpi, 1, 1, jpj, 1, 1.e-2, numout ) 830 ENDIF 831 832 IF(lwp) WRITE(numout,*) 833 IF(lwp) WRITE(numout,*) ' zgr_zps : z-coordinate with partial steps' 834 IF(lwp) WRITE(numout,*) ' ~~~~~~~ ' 835 IF(lwp) WRITE(numout,*) ' mbathy is recomputed : bathy_level file is NOT used' 836 837 838 ! bathymetry in level (from bathy_meter) 839 ! =================== 840 841 ! initialize mbathy to the maximum ocean level available 842 mbathy(:,:) = jpkm1 843 844 ! storage of land and island's number (zera and negative values) in mbathy 845 DO jj = 1, jpj 846 DO ji= 1, jpi 847 IF( bathy(ji,jj) <= 0. ) mbathy(ji,jj) = INT( bathy(ji,jj) ) 848 END DO 849 END DO 850 851 ! bounded value of bathy 852 ! minimum depth == 3 levels 853 ! maximum depth == gdepw_0(jpk)+e3t_0(jpk) 854 ! i.e. the last ocean level thickness cannot exceed e3t_0(jpkm1)+e3t_0(jpk) 855 DO jj = 1, jpj 856 DO ji= 1, jpi 857 IF( bathy(ji,jj) <= 0. ) THEN 858 bathy(ji,jj) = 0.e0 859 ELSE 860 bathy(ji,jj) = MAX( bathy(ji,jj), zmin ) 861 bathy(ji,jj) = MIN( bathy(ji,jj), zmax ) 862 ENDIF 863 END DO 864 END DO 865 866 ! Compute mbathy for ocean points (i.e. the number of ocean levels) 867 ! find the number of ocean levels such that the last level thickness 868 ! is larger than the minimum of e3zps_min and e3zps_rat * e3t_0 (where 869 ! e3t_0 is the reference level thickness 870 DO jk = jpkm1, 1, -1 871 zdepth = gdepw_0(jk) + MIN( e3zps_min, e3t_0(jk)*e3zps_rat ) 872 DO jj = 1, jpj 873 DO ji = 1, jpi 874 IF( 0. < bathy(ji,jj) .AND. bathy(ji,jj) <= zdepth ) mbathy(ji,jj) = jk-1 875 END DO 876 END DO 877 END DO 878 879 ! Scale factors and depth at T- and W-points 880 881 ! intitialization to the reference z-coordinate 882 DO jk = 1, jpk 883 gdept(:,:,jk) = gdept_0(jk) 884 gdepw(:,:,jk) = gdepw_0(jk) 885 e3t (:,:,jk) = e3t_0 (jk) 886 e3w (:,:,jk) = e3w_0 (jk) 887 END DO 888 hdept(:,:) = gdept(:,:,2 ) 889 hdepw(:,:) = gdepw(:,:,3 ) 890 891 ! 892 DO jj = 1, jpj 893 DO ji = 1, jpi 894 ik = mbathy(ji,jj) 895 ! ocean point only 896 IF( ik > 0 ) THEN 897 ! max ocean level case 898 IF( ik == jpkm1 ) THEN 899 zdepwp = bathy(ji,jj) 900 ze3tp = bathy(ji,jj) - gdepw_0(ik) 901 ze3wp = 0.5 * e3w_0(ik) * ( 1. + ( ze3tp/e3t_0(ik) ) ) 902 e3t(ji,jj,ik ) = ze3tp 903 e3t(ji,jj,ik+1) = ze3tp 904 e3w(ji,jj,ik ) = ze3wp 905 e3w(ji,jj,ik+1) = ze3tp 906 gdepw(ji,jj,ik+1) = zdepwp 907 gdept(ji,jj,ik ) = gdept_0(ik-1) + ze3wp 908 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + ze3tp 909 ! standard case 910 ELSE 911 !!alex 912 IF( bathy(ji,jj) <= gdepw_0(ik+1) ) THEN 913 gdepw(ji,jj,ik+1) = bathy(ji,jj) 914 ELSE 915 !!alex ctl write(*,*) 'zps',ji,jj,'bathy', bathy(ji,jj), 'depw ',gdepw_0(ik+1) 916 gdepw(ji,jj,ik+1) = gdepw_0(ik+1) 917 ENDIF 918 !!Alex 919 !!Alex gdepw(ji,jj,ik+1) = bathy(ji,jj) 920 !!bug gm verifier les gdepw_0 921 gdept(ji,jj,ik ) = gdepw_0(ik) + ( gdepw(ji,jj,ik+1) - gdepw_0(ik) ) & 922 & * ((gdept_0( ik ) - gdepw_0(ik) ) & 923 & / ( gdepw_0( ik+1) - gdepw_0(ik) )) 924 e3t(ji,jj,ik) = e3t_0(ik) * ( gdepw (ji,jj,ik+1) - gdepw_0(ik) ) & 925 & / ( gdepw_0( ik+1) - gdepw_0(ik) ) 926 e3w(ji,jj,ik) = 0.5 * ( gdepw(ji,jj,ik+1) + gdepw_0(ik+1) - 2.*gdepw_0(ik) ) & 927 & * ( e3w_0(ik) / ( gdepw_0(ik+1) - gdepw_0(ik) ) ) 928 ! ... on ik+1 929 e3w(ji,jj,ik+1) = e3t(ji,jj,ik) 930 e3t(ji,jj,ik+1) = e3t(ji,jj,ik) 931 gdept(ji,jj,ik+1) = gdept(ji,jj,ik) + e3t(ji,jj,ik) 932 ENDIF 933 ENDIF 934 END DO 935 END DO 936 937 it = 0 938 DO jj = 1, jpj 939 DO ji = 1, jpi 940 ik = mbathy(ji,jj) 941 ! ocean point only 942 IF( ik > 0 ) THEN 943 ! bathymetry output 944 hdept(ji,jj) = gdept(ji,jj,ik ) 945 hdepw(ji,jj) = gdepw(ji,jj,ik+1) 946 e3tp (ji,jj) = e3t(ji,jj,ik ) 947 e3wp (ji,jj) = e3w(ji,jj,ik ) 948 ! test 949 zdiff= gdepw(ji,jj,ik+1) - gdept(ji,jj,ik ) 950 IF( zdiff <= 0. .AND. lwp ) THEN 951 it=it+1 952 WRITE(numout,*) ' it = ', it, ' ik = ', ik, ' (i,j) = ', ji, jj 953 WRITE(numout,*) ' bathy = ', bathy(ji,jj) 954 WRITE(numout,*) ' gdept= ', gdept(ji,jj,ik), ' gdepw= ', gdepw(ji,jj,ik+1), & 955 ' zdiff = ', zdiff 956 WRITE(numout,*) ' e3tp = ', e3t(ji,jj,ik ), ' e3wp = ', e3w(ji,jj,ik ) 957 ENDIF 958 ENDIF 959 END DO 960 END DO 961 962 ! Scale factors and depth at U-, V-, UW and VW-points 963 964 ! initialisation to z-scale factors 965 DO jk = 1, jpk 966 e3u (:,:,jk) = e3t_0(jk) 967 e3v (:,:,jk) = e3t_0(jk) 968 e3uw(:,:,jk) = e3w_0(jk) 969 e3vw(:,:,jk) = e3w_0(jk) 970 END DO 971 972 ! Computed as the minimum of neighbooring scale factors 973 DO jk = 1,jpk 974 DO jj = 1, jpjm1 975 DO ji = 1, fs_jpim1 ! vector opt. 976 e3u (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji+1,jj,jk)) 977 e3v (ji,jj,jk) = MIN( e3t(ji,jj,jk), e3t(ji,jj+1,jk)) 978 e3uw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji+1,jj,jk) ) 979 e3vw(ji,jj,jk) = MIN( e3w(ji,jj,jk), e3w(ji,jj+1,jk) ) 980 END DO 981 END DO 982 END DO 983 984 ! Boundary conditions 985 CALL lbc_lnk( e3u , 'U', 1. ) ; CALL lbc_lnk( e3uw, 'U', 1. ) 986 CALL lbc_lnk( e3v , 'V', 1. ) ; CALL lbc_lnk( e3vw, 'V', 1. ) 987 988 ! set to z-scale factor if zero (i.e. along closed boundaries) 989 DO jk = 1, jpk 990 DO jj = 1, jpj 991 DO ji = 1, jpi 992 IF( e3u (ji,jj,jk) == 0.e0 ) e3u (ji,jj,jk) = e3t_0(jk) 993 IF( e3v (ji,jj,jk) == 0.e0 ) e3v (ji,jj,jk) = e3t_0(jk) 994 IF( e3uw(ji,jj,jk) == 0.e0 ) e3uw(ji,jj,jk) = e3w_0(jk) 995 IF( e3vw(ji,jj,jk) == 0.e0 ) e3vw(ji,jj,jk) = e3w_0(jk) 996 END DO 997 END DO 998 END DO 999 1000 ! Scale factor at F-point 1001 1002 ! initialisation to z-scale factors 1003 DO jk = 1, jpk 1004 e3f (:,:,jk) = e3t_0(jk) 1005 END DO 1006 1007 ! Computed as the minimum of neighbooring V-scale factors 1008 DO jk = 1, jpk 1009 DO jj = 1, jpjm1 1010 DO ji = 1, fs_jpim1 ! vector opt. 1011 e3f(ji,jj,jk) = MIN( e3v(ji,jj,jk), e3v(ji+1,jj,jk) ) 1012 END DO 1013 END DO 1014 END DO 1015 ! Boundary conditions 1016 CALL lbc_lnk( e3f, 'F', 1. ) 1017 1018 ! set to z-scale factor if zero (i.e. along closed boundaries) 1019 DO jk = 1, jpk 1020 DO jj = 1, jpj 1021 DO ji = 1, jpi 1022 IF( e3f(ji,jj,jk) == 0.e0 ) e3f(ji,jj,jk) = e3t_0(jk) 1023 END DO 1024 END DO 1025 END DO 1026 !!bug ? gm: must be a do loop with mj0,mj1 1027 1028 ! we duplicate factor scales for jj = 1 and jj = 2 1029 e3t(:,mj0(1),:) = e3t(:,mj0(2),:) 1030 e3w(:,mj0(1),:) = e3w(:,mj0(2),:) 1031 e3u(:,mj0(1),:) = e3u(:,mj0(2),:) 1032 e3v(:,mj0(1),:) = e3v(:,mj0(2),:) 1033 e3f(:,mj0(1),:) = e3f(:,mj0(2),:) 1034 1035 1036 1037 ! Compute gdep3w (vertical sum of e3w) 1038 1039 gdep3w (:,:,1) = 0.5 * e3w (:,:,1) 1040 DO jk = 2, jpk 1041 gdep3w (:,:,jk) = gdep3w (:,:,jk-1) + e3w (:,:,jk) 1042 END DO 1043 1044 ! Control print 1045 9600 FORMAT(9x,' level gdept gdepw e3t e3w ') 1046 9610 FORMAT(10x,i4,4f9.2) 1047 IF(lwp .AND. ll_print) THEN 1048 DO jj = 1,jpj 1049 DO ji = 1, jpi 1050 ik = MAX(mbathy(ji,jj),1) 1051 zprt(ji,jj) = e3t(ji,jj,ik) 1052 END DO 1053 END DO 1054 WRITE(numout,*) 1055 WRITE(numout,*) 'domzgr e3t(mbathy)' 1056 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1057 DO jj = 1,jpj 1058 DO ji = 1, jpi 1059 ik = MAX(mbathy(ji,jj),1) 1060 zprt(ji,jj) = e3w(ji,jj,ik) 1061 END DO 1062 END DO 1063 WRITE(numout,*) 1064 WRITE(numout,*) 'domzgr e3w(mbathy)' 1065 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1066 DO jj = 1,jpj 1067 DO ji = 1, jpi 1068 ik = MAX(mbathy(ji,jj),1) 1069 zprt(ji,jj) = e3u(ji,jj,ik) 1070 END DO 1071 END DO 1072 1073 WRITE(numout,*) 1074 WRITE(numout,*) 'domzgr e3u(mbathy)' 1075 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1076 DO jj = 1,jpj 1077 DO ji = 1, jpi 1078 ik = MAX(mbathy(ji,jj),1) 1079 zprt(ji,jj) = e3v(ji,jj,ik) 1080 END DO 1081 END DO 1082 WRITE(numout,*) 1083 WRITE(numout,*) 'domzgr e3v(mbathy)' 1084 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1085 DO jj = 1,jpj 1086 DO ji = 1, jpi 1087 ik = MAX(mbathy(ji,jj),1) 1088 zprt(ji,jj) = e3f(ji,jj,ik) 1089 END DO 1090 END DO 1091 1092 WRITE(numout,*) 1093 WRITE(numout,*) 'domzgr e3f(mbathy)' 1094 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1095 DO jj = 1,jpj 1096 DO ji = 1, jpi 1097 ik = MAX(mbathy(ji,jj),1) 1098 zprt(ji,jj) = gdep3w(ji,jj,ik) 1099 END DO 1100 END DO 1101 WRITE(numout,*) 1102 WRITE(numout,*) 'domzgr gdep3w(mbathy)' 1103 CALL prihre(zprt,jpi,jpj,1,jpi,1,1,jpj,1,1.e-3,numout) 1104 1105 ENDIF 1106 1107 1108 DO jk = 1,jpk 1109 DO jj = 1, jpj 1110 DO ji = 1, jpi 1111 IF( e3w(ji,jj,jk) <= 0. .or. e3t(ji,jj,jk) <= 0. ) THEN 1112 IF(lwp) THEN 1113 WRITE(numout,*) ' e r r o r : e3w or e3t =< 0 ' 1114 WRITE(numout,*) ' ========= --------------- ' 1115 WRITE(numout,*) 1116 ENDIF 1117 STOP 'domzgr.psteps' 1118 ENDIF 1119 IF( gdepw(ji,jj,jk) < 0. .or. gdept(ji,jj,jk) < 0. ) THEN 1120 IF(lwp) THEN 1121 WRITE(numout,*) ' e r r o r : gdepw or gdept < 0 ' 1122 WRITE(numout,*) ' ========= ------------------ ' 1123 WRITE(numout,*) 1124 ENDIF 1125 STOP 'domzgr.psteps' 1126 ENDIF 1127 END DO 1128 END DO 1129 END DO 1130 1131 IF(lwp) THEN 1132 WRITE(numout,*) ' e3t lev 21 ' 1133 CALL prihre(e3t(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1134 WRITE(numout,*) ' e3w lev 21 ' 1135 CALL prihre(e3w(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1136 WRITE(numout,*) ' e3u lev 21 ' 1137 CALL prihre(e3u(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1138 WRITE(numout,*) ' e3v lev 21 ' 1139 CALL prihre(e3v(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1140 WRITE(numout,*) ' e3f lev 21 ' 1141 CALL prihre(e3f(1,1,21),jpi,jpj,50,59,1,1,5,1,0.,numout) 1142 WRITE(numout,*) ' e3t lev 22 ' 1143 CALL prihre(e3t(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1144 WRITE(numout,*) ' e3w lev 22 ' 1145 CALL prihre(e3w(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1146 WRITE(numout,*) ' e3u lev 22 ' 1147 CALL prihre(e3u(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1148 WRITE(numout,*) ' e3v lev 22 ' 1149 CALL prihre(e3v(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1150 WRITE(numout,*) ' e3f lev 22 ' 1151 CALL prihre(e3f(1,1,22),jpi,jpj,50,59,1,1,5,1,0.,numout) 1152 ENDIF 1153 1154 ! =========== 1155 ! Zoom domain 1156 ! =========== 1157 1158 IF( lzoom ) CALL zgr_bat_zoom 1159 1160 ! ================ 1161 ! Bathymetry check 1162 ! ================ 1163 1164 CALL zgr_bat_ctl 1165 1166 END SUBROUTINE zgr_zps 1167 1168 1169 FUNCTION fssig( pk ) RESULT( pf ) 1170 !!---------------------------------------------------------------------- 1171 !! *** ROUTINE eos_init *** 1172 !! 1173 !! ** Purpose : provide the analytical function in s-coordinate 1174 !! 1175 !! ** Method : the function provide the non-dimensional position of 1176 !! T and W (i.e. between 0 and 1) 1177 !! T-points at integer values (between 1 and jpk) 1178 !! W-points at integer values - 1/2 (between 0.5 and jpk-0.5) 1179 !! 1180 !! Reference : ??? 1181 !!---------------------------------------------------------------------- 1182 REAL(wp), INTENT(in ) :: pk ! continuous "k" coordinate 1183 REAL(wp) :: pf ! sigma value 1184 !!---------------------------------------------------------------------- 1185 ! 1186 pf = ( TANH( theta * ( -(pk-0.5) / REAL(jpkm1) + thetb ) ) & 1187 & - TANH( thetb * theta ) ) & 1188 & * ( COSH( theta ) & 1189 & + COSH( theta * ( 2.e0 * thetb - 1.e0 ) ) ) & 1190 & / ( 2.e0 * SINH( theta ) ) 1191 ! 1192 END FUNCTION fssig 1193 740 1194 741 1195 SUBROUTINE zgr_sco … … 758 1212 !! hbatf = mi( mj( hbatt ) ) 759 1213 !! - Compute gsigt, gsigw, esigt, esigw from an analytical 760 !! function and its derivative given as statementfunction.761 !! gsigt(k) = fssig (k +0.5)762 !! gsigw(k) = fssig (k 763 !! esigt(k) = fsdsig(k +0.5)764 !! esigw(k) = fsdsig(k 1214 !! function and its derivative given as function. 1215 !! gsigt(k) = fssig (k ) 1216 !! gsigw(k) = fssig (k-0.5) 1217 !! esigt(k) = fsdsig(k ) 1218 !! esigw(k) = fsdsig(k-0.5) 765 1219 !! This routine is given as an example, it must be modified 766 1220 !! following the user s desiderata. nevertheless, the output as … … 781 1235 !! * Local declarations 782 1236 INTEGER :: ji, jj, jk, jl 783 REAL(wp) :: fssig, fsdsig, pfkk784 785 1237 INTEGER :: iip1, ijp1, iim1, ijm1 786 REAL(wp) :: & 787 fssigt, fssigw, zcoeft, zcoefw, & 788 zrmax, ztaper 1238 REAL(wp) :: zrmax, ztaper 789 1239 790 1240 REAL(wp), DIMENSION(jpi,jpj) :: & … … 792 1242 793 1243 NAMELIST/nam_zgr_sco/ sbot_max, sbot_min, theta, thetb, r_max 794 !!----------------------------------------------------------------------795 796 ! a. Sigma stretching coefficient797 fssigt(pfkk) = (tanh(theta*(-(pfkk-0.5)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &798 * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))799 fssigw(pfkk) = (tanh(theta*(-(pfkk-1.0)/FLOAT(jpkm1)+thetb))-tanh(thetb*theta)) &800 * (cosh(theta)+cosh(theta*(2.*thetb-1.)))/(2.*sinh(theta))801 802 ! b. Vertical derivative of sigma stretching coefficient803 804 !!bug fsdsig(pfkk)= put here the analytical derivative of fssigt and w ...805 1244 !!---------------------------------------------------------------------- 806 1245 … … 1043 1482 ! 2.1 Coefficients for model level depth at w- and t-levels 1044 1483 1045 !!bug gm : change the def of statment function....1046 1484 DO jk = 1, jpk 1047 gsigw(jk) = -fssig w(FLOAT(jk))1048 gsigt(jk) = -fssig t(FLOAT(jk))1485 gsigw(jk) = -fssig( FLOAT(jk)-0.5 ) 1486 gsigt(jk) = -fssig( FLOAT(jk) ) 1049 1487 END DO 1050 1488 1051 1489 IF( nprint == 1 .AND. lwp ) WRITE(numout,*) 'gsigw 1 jpk ', gsigw(1), gsigw(jpk) 1052 1053 !!org DO jk = 1, jpk1054 !!org gsigw(jk) = -fssig ( FLOAT(jk) )1055 !!org gsigt(jk) = -fssig ( FLOAT(jk)+0.5)1056 !!org END DO1057 1490 1058 1491 ! 2.2 Coefficients for vertical scale factors at w-, t- levels … … 1068 1501 esigw(1)=esigw(2) 1069 1502 esigt(jpk)=esigt(jpkm1) 1070 1071 !!org DO jk = 1, jpk1072 !!org esigw(jk)=fsdsig( FLOAT(jk))1073 !!org esigt(jk)=fsdsig( FLOAT(jk)+0.5)1074 !!org END DO1075 1503 1076 1504 ! 2.3 Coefficients for vertical depth as the sum of e3w scale factors … … 1205 1633 END SUBROUTINE zgr_sco 1206 1634 1207 #else1208 !!----------------------------------------------------------------------1209 !! Default option : Empty routine1210 !!----------------------------------------------------------------------1211 SUBROUTINE zgr_sco1212 END SUBROUTINE zgr_sco1213 1635 #endif 1214 1636
Note: See TracChangeset
for help on using the changeset viewer.