New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 3223 – NEMO

Changeset 3223


Ignore:
Timestamp:
2011-12-15T14:13:05+01:00 (12 years ago)
Author:
hliu
Message:

updated hpg_prj in dynhpg.F90: 1) some optimisation in c-spline subroutine to save some computational amount. 2) added some explanations text in Subroutine hpg_prj(). 3) modifications on the head information about hpg_prj

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/2011/dev_NEMO_MERGE_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r3186 r3223  
    1414   !!             -   !  2005-11  (G. Madec) style & small optimisation 
    1515   !!            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 options 
     16   !!            3.4  !  2011-11  (H. Liu) hpg_prj: Original code for s-coordinates 
     17   !!                 !           (A. Coward) suppression of hel, wdj and rot options 
    1818   !!---------------------------------------------------------------------- 
    1919 
    … …  
    717717      zrhh(:,:,:) = rhd(:,:,:) 
    718718       
    719       ! Preparing vertical density profile for hybrid-sco coordinate 
     719      ! Preparing vertical density profile "zrhh(:,:,:)" for hybrid-sco coordinate 
    720720      DO jj = 1, jpj 
    721721        DO ji = 1, jpi    
    … …  
    732732      END DO 
    733733 
     734      ! Transfer the depth of "T(:,:,:)" to vertical coordinate "zdeptht(:,:,:)" 
    734735      DO jj = 1, jpj 
    735736        DO ji = 1, jpi 
    … …  
    753754      ! Construct the vertical density profile with the  
    754755      ! constrained cubic spline interpolation 
     756      ! rho(z) = asp + bsp*z + csp*z^2 + dsp*z^3 
    755757      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
    756758 
    757       ! Calculate the hydrostatic pressure at T(ji,jj,1) 
     759      ! Integrate the hydrostatic pressure "zhpi(:,:,:)" at "T(ji,jj,1)" 
    758760      DO jj = 2, jpj 
    759761        DO ji = 2, jpi  
    … …  
    768770      END DO 
    769771 
    770       ! Calculate the pressure at T(ji,jj,2:jpkm1) 
     772      ! Calculate the pressure "zhpi(:,:,:)" at "T(ji,jj,2:jpkm1)" 
    771773      DO jk = 2, jpkm1                                   
    772774        DO jj = 2, jpj      
    … …  
    964966      !!           
    965967      !! ** 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 
    968969      !! 
    969970      !!---------------------------------------------------------------------- 
    … …  
    987988      jpkm1 = size(fsp,3) - 1 
    988989 
    989       ! Clean output arrays 
    990       asp = 0.0_wp 
    991       bsp = 0.0_wp 
    992       csp = 0.0_wp 
    993       dsp = 0.0_wp 
    994990       
    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) 
    10031027    
    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) 
    10181055    
    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 
    10571068       
    10581069   END SUBROUTINE cspline 
    … …  
    10691080      !!                extrapolation is also permitted (no value limit)  
    10701081      !! 
    1071       !! H.Liu, Jan 2009,  POL  
    10721082      !!---------------------------------------------------------------------- 
    10731083      IMPLICIT NONE 
Note: See TracChangeset for help on using the changeset viewer.