Changeset 3223 for branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC
- Timestamp:
- 2011-12-15T14:13:05+01:00 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r3186 r3223 14 14 !! - ! 2005-11 (G. Madec) style & small optimisation 15 15 !! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase 16 !! 3.4 ! 2011-11 ( A. Coward, H. Liu) introduction of prj scheme;17 !! ! suppression of hel, wdj and rot options16 !! 3.4 ! 2011-11 (H. Liu) hpg_prj: Original code for s-coordinates 17 !! ! (A. Coward) suppression of hel, wdj and rot options 18 18 !!---------------------------------------------------------------------- 19 19 … … 717 717 zrhh(:,:,:) = rhd(:,:,:) 718 718 719 ! Preparing vertical density profile for hybrid-sco coordinate719 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 720 720 DO jj = 1, jpj 721 721 DO ji = 1, jpi … … 732 732 END DO 733 733 734 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 734 735 DO jj = 1, jpj 735 736 DO ji = 1, jpi … … 753 754 ! Construct the vertical density profile with the 754 755 ! constrained cubic spline interpolation 756 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 755 757 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 756 758 757 ! Calculate the hydrostatic pressure at T(ji,jj,1)759 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 758 760 DO jj = 2, jpj 759 761 DO ji = 2, jpi … … 768 770 END DO 769 771 770 ! Calculate the pressure at T(ji,jj,2:jpkm1)772 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 771 773 DO jk = 2, jpkm1 772 774 DO jj = 2, jpj … … 964 966 !! 965 967 !! ** Method : f(x) = asp + bsp*x + csp*x^2 + dsp*x^3 966 !! Reference: K.W. Brodlie, A review of mehtods for curve and function 967 !! drawing, 1980 968 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 968 969 !! 969 970 !!---------------------------------------------------------------------- … … 987 988 jpkm1 = size(fsp,3) - 1 988 989 989 ! Clean output arrays990 asp = 0.0_wp991 bsp = 0.0_wp992 csp = 0.0_wp993 dsp = 0.0_wp994 990 995 DO ji = 1, jpi 996 DO jj = 1, jpj 997 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 998 DO jk = 2, jpkm1-1 999 zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 1000 zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1001 zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 1002 zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 991 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 992 DO ji = 1, jpi 993 DO jj = 1, jpj 994 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 995 ! DO jk = 2, jpkm1-1 996 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 997 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 998 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 999 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 1000 ! 1001 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1002 ! 1003 ! IF(zdf1 * zdf2 <= 0._wp) THEN 1004 ! zdf(jk) = 0._wp 1005 ! ELSE 1006 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1007 ! ENDIF 1008 ! END DO 1009 1010 !!Simply geometric average 1011 DO jk = 2, jpkm1-1 1012 zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1)) 1013 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk)) 1014 1015 IF(zdf1 * zdf2 <= 0._wp) THEN 1016 zdf(jk) = 0._wp 1017 ELSE 1018 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1019 ENDIF 1020 END DO 1021 1022 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1023 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1024 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1025 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1026 & 0.5_wp * zdf(jpkm1 - 1) 1003 1027 1004 zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 1005 1006 IF(zdf1 * zdf2 <= 0._wp) THEN 1007 zdf(jk) = 0._wp 1008 ELSE 1009 zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 1010 ENDIF 1011 END DO 1012 1013 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1014 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1015 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1016 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1017 & 0.5_wp * zdf(jpkm1 - 1) 1028 DO jk = 1, jpkm1 - 1 1029 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1030 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1031 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1032 zddf1 = -2._wp * ztmp1 + ztmp2 1033 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1034 zddf2 = 2._wp * ztmp1 - ztmp2 1035 1036 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1037 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1038 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1039 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1040 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1041 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1042 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1043 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1044 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1045 END DO 1046 END DO 1047 END DO 1048 1049 ELSE IF (polynomial_type == 2) THEN ! Linear 1050 DO ji = 1, jpi 1051 DO jj = 1, jpj 1052 DO jk = 1, jpkm1-1 1053 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1054 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1018 1055 1019 DO jk = 1, jpkm1 - 1 1020 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1021 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1022 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1023 zddf1 = -2._wp * ztmp1 + ztmp2 1024 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1025 zddf2 = 2._wp * ztmp1 - ztmp2 1026 1027 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1028 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1029 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1030 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1031 & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 1032 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 1033 & xsp(ji,jj,jk)**2 ) 1034 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 1035 & csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 1036 & dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 1037 END DO 1038 1039 ELSE IF (polynomial_type == 2) THEN ! Linear 1040 1041 DO jk = 1, jpkm1-1 1042 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1043 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1044 1045 dsp(ji,jj,jk) = 0._wp 1046 csp(ji,jj,jk) = 0._wp 1047 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1048 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1049 END DO 1050 1051 ELSE 1052 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1053 ENDIF 1054 1055 END DO 1056 END DO 1056 dsp(ji,jj,jk) = 0._wp 1057 csp(ji,jj,jk) = 0._wp 1058 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1059 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1060 END DO 1061 END DO 1062 END DO 1063 1064 ELSE 1065 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1066 ENDIF 1067 1057 1068 1058 1069 END SUBROUTINE cspline … … 1069 1080 !! extrapolation is also permitted (no value limit) 1070 1081 !! 1071 !! H.Liu, Jan 2009, POL1072 1082 !!---------------------------------------------------------------------- 1073 1083 IMPLICIT NONE
Note: See TracChangeset
for help on using the changeset viewer.