![(please configure the [header_logo] section in trac.ini)](site/little_nemo.jpg)
Changeset 4 for vendor/nemo/current/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
- Timestamp:
- 02/08/12 17:17:04 (12 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
vendor/nemo/current/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90
r1 r4 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 … … 61 61 !!---------------------------------------------------------------------- 62 62 !! NEMO/OPA 3.3 , NEMO Consortium (2010) 63 !! $Id: dynhpg.F90 3 186 2011-11-27 08:16:19Z smasson$63 !! $Id: dynhpg.F90 3294 2012-01-28 16:44:18Z rblod $ 64 64 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 65 65 !!---------------------------------------------------------------------- … … 76 76 !! - Save the trend (l_trddyn=T) 77 77 !!---------------------------------------------------------------------- 78 USE oce, ONLY: tsa ! (tsa) used as 2 3D workspace79 !!80 78 INTEGER, INTENT(in) :: kt ! ocean time-step index 81 79 REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdu, ztrdv … … 85 83 ! 86 84 IF( l_trddyn ) THEN ! Temporary saving of ua and va trends (l_trddyn) 87 ztrdu => tsa(:,:,:,1) 88 ztrdv => tsa(:,:,:,2) 89 ! 85 CALL wrk_alloc( jpi,jpj,jpk, ztrdu, ztrdv ) 90 86 ztrdu(:,:,:) = ua(:,:,:) 91 87 ztrdv(:,:,:) = va(:,:,:) … … 104 100 ztrdv(:,:,:) = va(:,:,:) - ztrdv(:,:,:) 105 101 CALL trd_mod( ztrdu, ztrdv, jpdyn_trd_hpg, 'DYN', kt ) 102 CALL wrk_dealloc( jpi,jpj,jpk, ztrdu, ztrdv ) 106 103 ENDIF 107 104 ! … … 191 188 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 192 189 !!---------------------------------------------------------------------- 193 USE oce, ONLY: tsa ! (tsa) used as 2 3D workspace194 !!195 190 INTEGER, INTENT(in) :: kt ! ocean time-step index 196 191 !! … … 199 194 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 200 195 !!---------------------------------------------------------------------- 201 202 zhpi => tsa(:,:,:,1) 203 zhpj => tsa(:,:,:,2) 196 ! 197 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 204 198 ! 205 199 IF( kt == nit000 ) THEN … … 245 239 END DO 246 240 ! 241 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 242 ! 247 243 END SUBROUTINE hpg_zco 248 244 … … 256 252 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 257 253 !!---------------------------------------------------------------------- 258 USE oce, ONLY: tsa ! (tsa) used as 2 3D workspace259 !!260 254 INTEGER, INTENT(in) :: kt ! ocean time-step index 261 255 !! … … 265 259 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 266 260 !!---------------------------------------------------------------------- 267 268 zhpi => tsa(:,:,:,1) 269 zhpj => tsa(:,:,:,2) 261 ! 262 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 270 263 ! 271 264 IF( kt == nit000 ) THEN … … 343 336 END DO 344 337 ! 345 338 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 339 ! 346 340 END SUBROUTINE hpg_zps 347 341 … … 365 359 !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 366 360 !!---------------------------------------------------------------------- 367 USE oce, ONLY: tsa ! (tsa) used as 2 3D workspace368 !!369 361 INTEGER, INTENT(in) :: kt ! ocean time-step index 370 362 !! … … 373 365 REAL(wp), POINTER, DIMENSION(:,:,:) :: zhpi, zhpj 374 366 !!---------------------------------------------------------------------- 375 376 zhpi => tsa(:,:,:,1) 377 zhpj => tsa(:,:,:,2) 367 ! 368 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zhpj ) 378 369 ! 379 370 IF( kt == nit000 ) THEN … … 432 423 END DO 433 424 ! 425 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zhpj ) 426 ! 434 427 END SUBROUTINE hpg_sco 435 428 … … 442 435 !! Reference: Shchepetkin and McWilliams, J. Geophys. Res., 108(C3), 3090, 2003 443 436 !!---------------------------------------------------------------------- 444 USE oce , ONLY: tsa ! (tsa) used as 2 3D workspace445 !!446 437 INTEGER, INTENT(in) :: kt ! ocean time-step index 447 438 !! … … 458 449 CALL wrk_alloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 459 450 CALL wrk_alloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 460 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k ) 461 zhpi => tsa(:,:,:,1) 462 zhpj => tsa(:,:,:,2) 451 CALL wrk_alloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 463 452 ! 464 453 … … 660 649 CALL wrk_dealloc( jpi, jpj, jpk, dzx , dzy , dzz , dzu , dzv , dzw ) 661 650 CALL wrk_dealloc( jpi, jpj, jpk, drhox, drhoy, drhoz, drhou, drhov, drhow ) 662 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k 651 CALL wrk_dealloc( jpi, jpj, jpk, rho_i, rho_j, rho_k, zhpi, zhpj ) 663 652 ! 664 653 END SUBROUTINE hpg_djc … … 678 667 !! 679 668 !!---------------------------------------------------------------------- 680 USE oce , ONLY: tsa ! (tsa) used as 2 3D workspace681 !!----------------------------------------------------------------------682 !!683 669 INTEGER, PARAMETER :: polynomial_type = 1 ! 1: cubic spline, 2: linear 684 670 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 697 683 !!---------------------------------------------------------------------- 698 684 ! 699 CALL wrk_alloc( jpi, jpj, jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 700 zdeptht => tsa(:,:,:,1) 701 zrhh => tsa(:,:,:,2) 685 CALL wrk_alloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 686 CALL wrk_alloc( jpi,jpj,jpk, zdeptht, zrhh ) 702 687 ! 703 688 IF( kt == nit000 ) THEN … … 717 702 zrhh(:,:,:) = rhd(:,:,:) 718 703 719 ! Preparing vertical density profile for hybrid-sco coordinate704 ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 720 705 DO jj = 1, jpj 721 706 DO ji = 1, jpi … … 732 717 END DO 733 718 719 ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 734 720 DO jj = 1, jpj 735 721 DO ji = 1, jpi … … 753 739 ! Construct the vertical density profile with the 754 740 ! constrained cubic spline interpolation 741 ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 755 742 CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type) 756 743 757 ! Calculate the hydrostatic pressure at T(ji,jj,1)744 ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 758 745 DO jj = 2, jpj 759 746 DO ji = 2, jpi … … 768 755 END DO 769 756 770 ! Calculate the pressure at T(ji,jj,2:jpkm1)757 ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 771 758 DO jk = 2, jpkm1 772 759 DO jj = 2, jpj … … 953 940 END DO 954 941 ! 955 CALL wrk_dealloc( jpi, jpj, jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 942 CALL wrk_dealloc( jpi,jpj,jpk, zhpi, zu, zv, fsp, xsp, asp, bsp, csp, dsp ) 943 CALL wrk_dealloc( jpi,jpj,jpk, zdeptht, zrhh ) 956 944 ! 957 945 END SUBROUTINE hpg_prj … … 964 952 !! 965 953 !! ** 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 954 !! Reference: CJC Kruger, Constrained Cubic Spline Interpoltation 968 955 !! 969 956 !!---------------------------------------------------------------------- … … 987 974 jpkm1 = size(fsp,3) - 1 988 975 989 ! Clean output arrays990 asp = 0.0_wp991 bsp = 0.0_wp992 csp = 0.0_wp993 dsp = 0.0_wp994 976 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 977 IF (polynomial_type == 1) THEN ! Constrained Cubic Spline 978 DO ji = 1, jpi 979 DO jj = 1, jpj 980 !!Fritsch&Butland's method, 1984 (preferred, but more computation) 981 ! DO jk = 2, jpkm1-1 982 ! zdxtmp1 = xsp(ji,jj,jk) - xsp(ji,jj,jk-1) 983 ! zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 984 ! zdf1 = ( fsp(ji,jj,jk) - fsp(ji,jj,jk-1) ) / zdxtmp1 985 ! zdf2 = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp2 986 ! 987 ! zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 988 ! 989 ! IF(zdf1 * zdf2 <= 0._wp) THEN 990 ! zdf(jk) = 0._wp 991 ! ELSE 992 ! zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 993 ! ENDIF 994 ! END DO 995 996 !!Simply geometric average 997 DO jk = 2, jpkm1-1 998 zdf1 = (fsp(ji,jj,jk) - fsp(ji,jj,jk-1)) / (xsp(ji,jj,jk) - xsp(ji,jj,jk-1)) 999 zdf2 = (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / (xsp(ji,jj,jk+1) - xsp(ji,jj,jk)) 1000 1001 IF(zdf1 * zdf2 <= 0._wp) THEN 1002 zdf(jk) = 0._wp 1003 ELSE 1004 zdf(jk) = 2._wp * zdf1 * zdf2 / (zdf1 + zdf2) 1005 ENDIF 1006 END DO 1007 1008 zdf(1) = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / & 1009 & ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - 0.5_wp * zdf(2) 1010 zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 1011 & ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 1012 & 0.5_wp * zdf(jpkm1 - 1) 1003 1013 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) 1014 DO jk = 1, jpkm1 - 1 1015 zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1016 ztmp1 = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 1017 ztmp2 = 6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 1018 zddf1 = -2._wp * ztmp1 + ztmp2 1019 ztmp1 = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 1020 zddf2 = 2._wp * ztmp1 - ztmp2 1021 1022 dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 1023 csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 1024 bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - & 1025 & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 1026 & dsp(ji,jj,jk) * ((xsp(ji,jj,jk+1) + xsp(ji,jj,jk))**2 - & 1027 & xsp(ji,jj,jk+1) * xsp(ji,jj,jk)) 1028 asp(ji,jj,jk) = fsp(ji,jj,jk) - xsp(ji,jj,jk) * (bsp(ji,jj,jk) + & 1029 & (xsp(ji,jj,jk) * (csp(ji,jj,jk) + & 1030 & dsp(ji,jj,jk) * xsp(ji,jj,jk)))) 1031 END DO 1032 END DO 1033 END DO 1034 1035 ELSE IF (polynomial_type == 2) THEN ! Linear 1036 DO ji = 1, jpi 1037 DO jj = 1, jpj 1038 DO jk = 1, jpkm1-1 1039 zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk) 1040 ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 1018 1041 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 1042 dsp(ji,jj,jk) = 0._wp 1043 csp(ji,jj,jk) = 0._wp 1044 bsp(ji,jj,jk) = ztmp1 / zdxtmp 1045 asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 1046 END DO 1047 END DO 1048 END DO 1049 1050 ELSE 1051 CALL ctl_stop( 'invalid polynomial type in cspline' ) 1052 ENDIF 1053 1057 1054 1058 1055 END SUBROUTINE cspline … … 1069 1066 !! extrapolation is also permitted (no value limit) 1070 1067 !! 1071 !! H.Liu, Jan 2009, POL1072 1068 !!---------------------------------------------------------------------- 1073 1069 IMPLICIT NONE
Note: See TracChangeset
for help on using the changeset viewer.
