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 3074 – NEMO

Changeset 3074


Ignore:
Timestamp:
2011-11-10T13:43:36+01:00 (12 years ago)
Author:
acc
Message:

Branch dev_NOC_2011_MERGE. #874. Step 11. Merge in changes from the dev_r2802_NOCL_prjhpg. Partially imposed coding standards (more could be done); added missing namelist changes and added first attempt at documentation

Location:
branches/2011/dev_NOC_2011_MERGE
Files:
8 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Chapters/Chap_DYN.tex

    r2986 r3074  
    633633$\bullet$ Rotated axes scheme (rot) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_rot}=true) 
    634634 
    635 Note that expression \eqref{Eq_dynhpg_sco} is used when the variable volume  
     635$\bullet$ Pressure Jacobian scheme (prj) \citep{Thiem_Berntsen_OM06} (\np{ln\_dynhpg\_prj}=true) 
     636 
     637Note that expression \eqref{Eq_dynhpg_sco} is commonly used when the variable volume  
    636638formulation is activated (\key{vvl}) because in that case, even with a flat bottom,  
    637639the coordinate surfaces are not horizontal but follow the free surface  
    638 \citep{Levier2007}. The other pressure gradient options are not yet available. 
     640\citep{Levier2007}. Only the pressure jacobian scheme (\np{ln\_dynhpg\_prj}=true) is available as an  
     641alternative to the default \np{ln\_dynhpg\_sco}=true when \key{vvl} is active.  The pressure Jacobian scheme uses  
     642a constrained cubic spline to reconstruct the density profile across the water column. This method 
     643maintains the monotonicity between the density nodes and is of a higher order than the linear 
     644interpolation method. The pressure can be calculated by analytical integration of the density profile and 
     645a pressure Jacobian method is used to solve the horizontal pressure gradient. This method should 
     646provide a more accurate calculation of the horizontal pressure gradient than the standard scheme. 
    639647 
    640648%-------------------------------------------------------------------------------------------------------------- 
  • branches/2011/dev_NOC_2011_MERGE/DOC/TexFiles/Namelist/namdyn_hpg

    r2540 r3074  
    99   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    1010   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     11   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    1112   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    1213   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/GYRE/EXP00/namelist

    r3072 r3074  
    535535   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    536536   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     537   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    537538   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    538539   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist

    r3072 r3074  
    533533   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    534534   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     535   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    535536   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    536537   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist

    r3072 r3074  
    535535   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    536536   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     537   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    537538   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    538539   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist

    r3072 r3074  
    536536   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    537537   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     538   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    538539   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    539540   ln_dynhpg_imp = .false. !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/CONFIG/POMME/EXP00/namelist

    r3072 r3074  
    535535   ln_hpg_djc  = .false.   !  s-coordinate (Density Jacobian with Cubic polynomial) 
    536536   ln_hpg_rot  = .false.   !  s-coordinate (ROTated axes scheme) 
     537   ln_hpg_prj  = .false.   !  s-coordinate (Pressure Jacobian scheme) 
    537538   rn_gamma    = 0.e0      !  weighting coefficient (wdj scheme) 
    538539   ln_dynhpg_imp = .true.  !  time stepping: semi-implicit time scheme  (T) 
  • branches/2011/dev_NOC_2011_MERGE/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r2715 r3074  
    2727   !!       hpg_djc  : s-coordinate (Density Jacobian with Cubic polynomial) 
    2828   !!       hpg_rot  : s-coordinate (ROTated axes scheme) 
     29   !!       hpg_prj  : s-coordinate (Pressure Jacobian with Cubic polynomial) 
    2930   !!---------------------------------------------------------------------- 
    3031   USE oce             ! ocean dynamics and tracers 
     
    5253   LOGICAL , PUBLIC ::   ln_hpg_djc    = .FALSE.   !: s-coordinate (Density Jacobian with Cubic polynomial) 
    5354   LOGICAL , PUBLIC ::   ln_hpg_rot    = .FALSE.   !: s-coordinate (ROTated axes scheme) 
     55   LOGICAL , PUBLIC ::   ln_hpg_prj    = .FALSE.   !: s-coordinate (Pressure Jacobian scheme) 
    5456   REAL(wp), PUBLIC ::   rn_gamma      = 0._wp     !: weighting coefficient 
    5557   LOGICAL , PUBLIC ::   ln_dynhpg_imp = .FALSE.   !: semi-implicite hpg flag 
    5658 
    57    INTEGER  ::   nhpg  =  0   ! = 0 to 6, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
     59   INTEGER  ::   nhpg  =  0   ! = 0 to 7, type of pressure gradient scheme used ! (deduced from ln_hpg_... flags) 
    5860 
    5961   !! * Substitutions 
     
    100102      CASE (  5 )   ;   CALL hpg_djc    ( kt )      ! s-coordinate (Density Jacobian with Cubic polynomial) 
    101103      CASE (  6 )   ;   CALL hpg_rot    ( kt )      ! s-coordinate (ROTated axes scheme) 
     104      CASE (  7 )   ;   CALL hpg_prj    ( kt )      ! s-coordinate (Pressure Jacobian scheme) 
    102105      END SELECT 
    103106      ! 
     
    129132      !! 
    130133      NAMELIST/namdyn_hpg/ ln_hpg_zco, ln_hpg_zps, ln_hpg_sco, ln_hpg_hel,    & 
    131          &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, rn_gamma  , ln_dynhpg_imp 
     134         &                 ln_hpg_wdj, ln_hpg_djc, ln_hpg_rot, ln_hpg_prj,    & 
     135         &                 rn_gamma  , ln_dynhpg_imp 
    132136      !!---------------------------------------------------------------------- 
    133137      ! 
     
    147151         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    148152         WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot 
     153         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    149154         WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma 
    150155         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp 
    151156      ENDIF 
    152157      ! 
    153       IF( lk_vvl .AND. .NOT. ln_hpg_sco )   & 
    154          &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require the standard jacobian formulation hpg_sco') 
     158      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     159         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl requires:& 
     160                           & the standard jacobian formulation hpg_sco or & 
     161                           & the pressure jacobian formulation hpg_prj') 
    155162      ! 
    156163      !                               ! Set nhpg from ln_hpg_... flags 
     
    162169      IF( ln_hpg_djc )   nhpg = 5 
    163170      IF( ln_hpg_rot )   nhpg = 6 
     171      IF( ln_hpg_prj )   nhpg = 7 
    164172      ! 
    165173      !                               ! Consitency check 
     
    172180      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    173181      IF( ln_hpg_rot )   ioptio = ioptio + 1 
     182      IF( ln_hpg_prj )   ioptio = ioptio + 1 
    174183      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    175184      ! 
     
    10051014   END SUBROUTINE hpg_rot 
    10061015 
     1016    
     1017   SUBROUTINE hpg_prj( kt ) 
     1018      !!--------------------------------------------------------------------- 
     1019      !!                  ***  ROUTINE hpg_prj  *** 
     1020      !! 
     1021      !! ** Method  :   s-coordinate case. 
     1022      !!      Reformulate the horizontal hydrostatical pressure gradient 
     1023      !!      term using Pressure Jacobian. A new correction term  
     1024      !!      is developed to eliminate the sigma-coordinate error. 
     1025      !! 
     1026      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     1027      !!             - Save the trend (l_trddyn=T) 
     1028      !! 
     1029      !!---------------------------------------------------------------------- 
     1030 
     1031      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
     1032      USE oce     , ONLY:   fsde3w_tmp => ta       ! (ta,sa) used as 3D workspace 
     1033      USE oce     , ONLY:      rhd_tmp => sa  
     1034      USE wrk_nemo, ONLY:         zhpi => wrk_3d_3  
     1035      USE wrk_nemo, ONLY:           zu => wrk_3d_4  
     1036      USE wrk_nemo, ONLY:           zv => wrk_3d_5 
     1037      USE wrk_nemo, ONLY:          fsp => wrk_3d_6  
     1038      USE wrk_nemo, ONLY:          xsp => wrk_3d_7 
     1039      USE wrk_nemo, ONLY:          asp => wrk_3d_8 
     1040      USE wrk_nemo, ONLY:          bsp => wrk_3d_9 
     1041      USE wrk_nemo, ONLY:          csp => wrk_3d_10 
     1042      USE wrk_nemo, ONLY:          dsp => wrk_3d_11 
     1043      !! 
     1044      !!---------------------------------------------------------------------- 
     1045      !! 
     1046      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
     1047      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     1048      !! 
     1049      INTEGER  ::   ji, jj, jk, jkk                ! dummy loop indices 
     1050      REAL(wp) ::   zcoef0, znad               ! temporary scalars 
     1051      !! 
     1052      !! The local varialbes for the correction term 
     1053      INTEGER  :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd 
     1054      REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs 
     1055      REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps 
     1056      REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2 
     1057      REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 
     1058      INTEGER  :: bhitwe, bhitns 
     1059      !!---------------------------------------------------------------------- 
     1060 
     1061      IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 
     1062         CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable')   ;   RETURN 
     1063      ENDIF 
     1064 
     1065      IF( kt == nit000 ) THEN 
     1066         IF(lwp) WRITE(numout,*) 
     1067         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 
     1068         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian' 
     1069      ENDIF 
     1070 
     1071      !!---------------------------------------------------------------------- 
     1072      ! Local constant initialization 
     1073      zcoef0 = - grav  
     1074      znad = 0.0_wp 
     1075      IF(lk_vvl) znad = 1._wp 
     1076 
     1077      ! Save fsde3w and rhd 
     1078      fsde3w_tmp(:,:,:) = fsde3w(:,:,:)  
     1079      rhd_tmp(:,:,:)    = rhd(:,:,:)    
     1080 
     1081      ! Clean 3-D work arraies 
     1082      zhpi(:,:,:) = 0. 
     1083       
     1084      !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here 
     1085 
     1086      ! Preparing vertical density profile for hybrid-sco coordinate 
     1087      DO jj = 1, jpj 
     1088        DO ji = 1, jpi    
     1089          jk = mbathy(ji,jj) 
     1090          IF(jk <= 0) THEN; rhd(ji,jj,:) = 0._wp 
     1091            ELSE IF(jk == 1) THEN; rhd(ji,jj, jk+1:jpk) = rhd(ji,jj,jk) 
     1092            ELSE IF(jk < jpkm1) THEN 
     1093            DO jkk = jk+1, jpk 
     1094              rhd(ji,jj,jkk) = interp1(fsde3w(ji,jj,jkk), fsde3w(ji,jj,jkk-1),& 
     1095                                       fsde3w(ji,jj,jkk-2), rhd(ji,jj,jkk-1), rhd(ji,jj,jkk-2)) 
     1096            END DO  
     1097          END IF 
     1098        END DO 
     1099      END DO 
     1100 
     1101      DO jj = 1, jpj 
     1102        DO ji = 1, jpi 
     1103          fsde3w(ji,jj,1) = 0.5_wp * fse3w(ji,jj,1) 
     1104          fsde3w(ji,jj,1) = fsde3w(ji,jj,1) - sshn(ji,jj) * znad 
     1105          DO jk = 2, jpk 
     1106            fsde3w(ji,jj,jk) = fsde3w(ji,jj,jk-1) + fse3w(ji,jj,jk) 
     1107          END DO 
     1108        END DO 
     1109      END DO 
     1110 
     1111      DO jk = 1, jpkm1 
     1112        DO jj = 1, jpj 
     1113          DO ji = 1, jpi 
     1114            fsp(ji,jj,jk) = rhd(ji,jj,jk) 
     1115            xsp(ji,jj,jk) = fsde3w(ji,jj,jk) 
     1116          END DO 
     1117        END DO 
     1118      END DO 
     1119 
     1120      !                 ! Construct the vertical density profile with the  
     1121      !                 !Constrained cubic spline interpolation 
     1122      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
     1123 
     1124      ! Calculate the hydrostatic pressure at T(ji,jj,1) 
     1125      DO jj = 2, jpj 
     1126        DO ji = 2, jpi  
     1127          rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1128                                         bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1129                               * 0.5_wp * fsde3w(ji,jj,1) 
     1130          rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1131 
     1132     !                  ! assuming linear profile across the top half surface layer 
     1133          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1   
     1134        END DO 
     1135      END DO 
     1136 
     1137      ! Calculate the pressure at T(ji,jj,2:jpkm1) 
     1138      DO jk = 2, jpkm1                                   
     1139        DO jj = 2, jpj      
     1140          DO ji = 2, jpi 
     1141            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 
     1142                             integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),& 
     1143                                    asp(ji,jj,jk-1),bsp(ji,jj,jk-1),& 
     1144                                    csp(ji,jj,jk-1),dsp(ji,jj,jk-1)) 
     1145          END DO 
     1146        END DO 
     1147      END DO 
     1148 
     1149      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 
     1150 
     1151      DO jj = 2, jpjm1      
     1152        DO ji = 2, jpim1   
     1153          zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad) 
     1154          zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad) 
     1155        END DO 
     1156      END DO 
     1157 
     1158      DO jk = 2, jpkm1                                   
     1159        DO jj = 2, jpjm1      
     1160          DO ji = 2, jpim1   
     1161            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk) 
     1162            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk) 
     1163          END DO 
     1164        END DO 
     1165      END DO 
     1166                
     1167      !  Start pressure integration 
     1168 
     1169      DO jk = 1, jpkm1                                   
     1170        DO jj = 2, jpjm1      
     1171          DO ji = 2, jpim1   
     1172            pwes = 0._wp; pwed = 0._wp 
     1173            pnss = 0._wp; pnsd = 0._wp 
     1174            zuijk = zu(ji,jj,jk) 
     1175            zvijk = zv(ji,jj,jk) 
     1176 
     1177            !!!!!     for u equation 
     1178            IF( -fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN 
     1179              jis = ji + 1; jid = ji 
     1180            ELSE 
     1181              jis = ji;     jid = ji +1 
     1182            ENDIF 
     1183 
     1184            !     !integrate the pressure on the shallow side 
     1185            jk1 = jk  
     1186            bhitwe = 0 
     1187            IF( jk1 == mbathy(jis,jj) ) THEN 
     1188              bhitwe = 1 
     1189            ELSE 
     1190              DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk ) 
     1191                pwes = pwes + &  
     1192                       integ2(fsde3w(jis,jj,jk1), fsde3w(jis,jj,jk1+1),& 
     1193                              asp(jis,jj,jk1)   , bsp(jis,jj,jk1),     & 
     1194                              csp(jis,jj,jk1)   , dsp(jis,jj,jk1)) 
     1195                jk1 = jk1 + 1 
     1196                IF( jk1 == mbathy(jis,jj) ) THEN 
     1197                  bhitwe = 1 
     1198                  EXIT 
     1199                END IF 
     1200              END DO 
     1201            ENDIF 
     1202             
     1203            IF( bhitwe /= 1 ) THEN 
     1204              pwes = pwes + &  
     1205                     integ2(fsde3w(jis,jj,jk1), -zuijk,      & 
     1206                            asp(jis,jj,jk1), bsp(jis,jj,jk1),& 
     1207                            csp(jis,jj,jk1), dsp(jis,jj,jk1)) 
     1208            ELSE 
     1209               zuijk = -fsde3w(jis,jj,jk1) 
     1210            ENDIF 
     1211 
     1212            !     !integrate the pressure on the deep side 
     1213            jk1 = jk  
     1214            bhitwe = 0 
     1215            IF( jk1 == 1 ) THEN 
     1216              bhitwe = 1 
     1217            ELSE 
     1218              DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk ) 
     1219                pwed = pwed + &  
     1220                       integ2(fsde3w(jid,jj,jk1-1), fsde3w(jid,jj,jk1),& 
     1221                              asp(jid,jj,jk1-1),    bsp(jid,jj,jk1-1), & 
     1222                              csp(jid,jj,jk1-1),    dsp(jid,jj,jk1-1)) 
     1223                jk1 = jk1 - 1 
     1224                IF( jk1 == 1 ) THEN 
     1225                  bhitwe = 1 
     1226                  EXIT 
     1227                END IF 
     1228              END DO 
     1229            ENDIF 
     1230             
     1231            IF( bhitwe /= 1 ) THEN 
     1232              pwed = pwed + &  
     1233                     integ2(-zuijk,             fsde3w(jid,jj,jk1),& 
     1234                             asp(jid,jj,jk1-1), bsp(jid,jj,jk1-1), & 
     1235                             csp(jid,jj,jk1-1), dsp(jid,jj,jk1-1)) 
     1236            ELSE 
     1237              deps  = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad) 
     1238              rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1), asp(jid,jj,1), & 
     1239                                              bsp(jid,jj,1),    csp(jid,jj,1), & 
     1240                                              dsp(jid,jj,1) ) * deps 
     1241              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1242              pwed  = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps 
     1243            ENDIF 
     1244 
     1245            IF( jid > jis ) THEN 
     1246              pe = pwed; pw = pwes 
     1247            ELSE 
     1248              pe = pwes; pw = pwed 
     1249            ENDIF 
     1250 
     1251            dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
     1252            IF( lk_vvl ) THEN 
     1253              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj)))  
     1254             ELSE 
     1255              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe)  
     1256            ENDIF 
     1257 
     1258            ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2) * & 
     1259               &           umask(ji,jj,jk) * tmask(ji,jj,jk) * tmask(ji+1,jj,jk) 
     1260   
     1261            !!!!!     for v equation 
     1262 
     1263            IF( -fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj)) ) THEN 
     1264              jjs = jj + 1; jjd = jj 
     1265            ELSE 
     1266              jjs = jj    ; jjd = jj + 1 
     1267            ENDIF 
     1268 
     1269            !     !integrate the pressure on the shallow side 
     1270            jk1 = jk  
     1271            bhitns = 0 
     1272            IF( jk1 == mbathy(ji,jjs) ) THEN 
     1273              bhitns = 1 
     1274            ELSE 
     1275              DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk ) 
     1276                pnss = pnss + &  
     1277                       integ2(fsde3w(ji,jjs,jk1), fsde3w(ji,jjs,jk1+1),& 
     1278                              asp(ji,jjs,jk1),    bsp(ji,jjs,jk1),     & 
     1279                              csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) ) 
     1280                jk1 = jk1 + 1 
     1281                IF( jk1 == mbathy(ji,jjs) ) THEN 
     1282                  bhitns = 1 
     1283                  EXIT 
     1284                END IF 
     1285              END DO 
     1286            ENDIF 
     1287             
     1288            IF( bhitns /= 1 ) THEN 
     1289              pnss = pnss + &  
     1290                     integ2(fsde3w(ji,jjs,jk1), -zvijk,          & 
     1291                            asp(ji,jjs,jk1),    bsp(ji,jjs,jk1), & 
     1292                            csp(ji,jjs,jk1),    dsp(ji,jjs,jk1) ) 
     1293            ELSE 
     1294               zvijk = -fsde3w(ji,jjs,jk1) 
     1295            ENDIF 
     1296 
     1297            !     !integrate the pressure on the deep side 
     1298            jk1 = jk  
     1299            bhitns = 0 
     1300            IF( jk1 == 1 ) THEN 
     1301              bhitns = 1 
     1302            ELSE 
     1303              DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk ) 
     1304                pnsd = pnsd + &  
     1305                       integ2(fsde3w(ji,jjd,jk1-1), fsde3w(ji,jjd,jk1), & 
     1306                              asp(ji,jjd,jk1-1),    bsp(ji,jjd,jk1-1),  & 
     1307                              csp(ji,jjd,jk1-1),    dsp(ji,jjd,jk1-1) ) 
     1308                jk1 = jk1 - 1 
     1309                IF( jk1 == 1 ) THEN 
     1310                  bhitns = 1 
     1311                  EXIT 
     1312                END IF 
     1313              END DO 
     1314            ENDIF 
     1315             
     1316            IF( bhitns /= 1 ) THEN 
     1317              pnsd = pnsd + &  
     1318                     integ2(-zvijk,            fsde3w(ji,jjd,jk1), & 
     1319                            asp(ji,jjd,jk1-1), bsp(ji,jjd,jk1-1),  & 
     1320                            csp(ji,jjd,jk1-1), dsp(ji,jjd,jk1-1) ) 
     1321            ELSE 
     1322              deps  = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad) 
     1323              rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1), asp(ji,jjd,1), & 
     1324                                              bsp(ji,jjd,1),    csp(ji,jjd,1), & 
     1325                                              dsp(ji,jjd,1) ) * deps 
     1326              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1327              pnsd  = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps 
     1328            ENDIF 
     1329 
     1330            IF( jjd > jjs ) THEN 
     1331              pn = pnsd; ps = pnss 
     1332            ELSE 
     1333              pn = pnss; ps = pnsd 
     1334            ENDIF 
     1335 
     1336            dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     1337            if( lk_vvl ) then 
     1338                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj)))  
     1339              else 
     1340                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn )  
     1341            end if 
     1342 
     1343            va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)* & 
     1344            &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1345 
     1346                     
     1347           END DO 
     1348        END DO 
     1349      END DO 
     1350 
     1351      ! Restore fsde3w and rhd 
     1352      fsde3w(:,:,:) = fsde3w_tmp(:,:,:) 
     1353      rhd(:,:,:)    = rhd_tmp(:,:,:) 
     1354 
     1355      ! 
     1356      IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) )   & 
     1357         CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 
     1358      ! 
     1359   END SUBROUTINE hpg_prj 
     1360 
     1361   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
     1362      !!---------------------------------------------------------------------- 
     1363      !!                 ***  ROUTINE cspline  *** 
     1364      !!        
     1365      !! ** Purpose :   constrained cubic spline interpolation 
     1366      !!           
     1367      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3  
     1368      !! Reference: K.W. Brodlie, A review of mehtods for curve and function 
     1369      !!                          drawing, 1980 
     1370      !! 
     1371      !!---------------------------------------------------------------------- 
     1372      IMPLICIT NONE 
     1373      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate 
     1374      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of  
     1375                                                                    ! the interpoated function 
     1376      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline  
     1377                                                                    ! 2: Linear 
     1378 
     1379      ! Local Variables       
     1380      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     1381      INTEGER  ::   jpi, jpj, jpkm1 
     1382      REAL(wp) ::   zdf1, zdf2, zddf1, zddf2, ztmp1, ztmp2, zdxtmp 
     1383      REAL(wp) ::   zdxtmp1, zdxtmp2, zalpha 
     1384      REAL(wp) ::   zdf(size(fsp,3)) 
     1385      !!---------------------------------------------------------------------- 
     1386 
     1387      jpi   = size(fsp,1) 
     1388      jpj   = size(fsp,2) 
     1389      jpkm1 = size(fsp,3) - 1 
     1390 
     1391      ! Clean output arrays 
     1392      asp = 0.0_wp 
     1393      bsp = 0.0_wp 
     1394      csp = 0.0_wp 
     1395      dsp = 0.0_wp 
     1396       
     1397       
     1398      Do ji = 1, jpi 
     1399      Do jj = 1, jpj 
     1400 
     1401      If (polynomial_type == 1) Then     !Constrained Cubic Spline 
     1402          Do jk = 2, jpkm1-1 
     1403             zdxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)   
     1404             zdxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)   
     1405             zdf1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / zdxtmp1 
     1406             zdf2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / zdxtmp2 
     1407   
     1408             zalpha = ( zdxtmp1 + 2._wp * zdxtmp2 ) / ( zdxtmp1 + zdxtmp2 ) / 3._wp 
     1409              
     1410             If(zdf1 * zdf2 <= 0._wp) Then 
     1411               zdf(jk) = 0._wp 
     1412             Else 
     1413               zdf(jk) = zdf1 * zdf2 / ( ( 1._wp - zalpha ) * zdf1 + zalpha * zdf2 ) 
     1414             End If 
     1415          End Do 
     1416   
     1417          zdf(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - & 
     1418                    &  0.5_wp * zdf(2) 
     1419          zdf(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 
     1420                    &           ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 
     1421                    & 0.5_wp * zdf(jpkm1 - 1) 
     1422   
     1423          Do jk = 1, jpkm1 - 1 
     1424             zdxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1425             ztmp1  = (zdf(jk+1) + 2._wp * zdf(jk)) / zdxtmp 
     1426             ztmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / zdxtmp / zdxtmp 
     1427             zddf1  = -2._wp * ztmp1 + ztmp2  
     1428             ztmp1  = (2._wp * zdf(jk+1) + zdf(jk)) / zdxtmp 
     1429             zddf2  =  2._wp * ztmp1 - ztmp2  
     1430   
     1431             dsp(ji,jj,jk) = (zddf2 - zddf1) / 6._wp / zdxtmp 
     1432             csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * zddf1 - xsp(ji,jj,jk)*zddf2 ) / 2._wp / zdxtmp 
     1433             bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / zdxtmp - &  
     1434                           & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 
     1435                           & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 
     1436                           &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 
     1437                           &                   xsp(ji,jj,jk)**2 ) 
     1438             asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 
     1439                           &                 csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 
     1440                           &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 
     1441          End Do 
     1442 
     1443       Else If (polynomial_type == 2) Then     !Linear 
     1444 
     1445        Do jk = 1, jpkm1-1 
     1446           zdxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1447           ztmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 
     1448   
     1449           dsp(ji,jj,jk) = 0._wp 
     1450           csp(ji,jj,jk) = 0._wp 
     1451           bsp(ji,jj,jk) = ztmp1 / zdxtmp 
     1452           asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 
     1453        End Do 
     1454 
     1455       Else 
     1456        CALL ctl_stop( 'invalid polynomial type in cspline' ) 
     1457      End If 
     1458 
     1459      End Do 
     1460      End Do 
     1461       
     1462   End Subroutine cspline 
     1463 
     1464 
     1465   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)  
     1466      !!---------------------------------------------------------------------- 
     1467      !!                 ***  ROUTINE interp1  *** 
     1468      !!        
     1469      !! ** Purpose :   1-d linear interpolation 
     1470      !!           
     1471      !! ** Method  :   
     1472      !!                interpolation is straight forward 
     1473      !!                extrapolation is also permitted (no value limit)  
     1474      !! 
     1475      !! H.Liu, Jan 2009,  POL  
     1476      !!---------------------------------------------------------------------- 
     1477      IMPLICIT NONE 
     1478      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr    
     1479      REAL(wp)             ::  f ! result of the interpolation (extrapolation) 
     1480      REAL(wp)             ::  zdeltx 
     1481      !!---------------------------------------------------------------------- 
     1482 
     1483      zdeltx = xr - xl 
     1484      IF(abs(zdeltx) <= 10._wp * EPSILON(x)) THEN 
     1485        f = 0.5_wp * (fl + fr) 
     1486       ELSE 
     1487        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / zdeltx 
     1488      END IF 
     1489       
     1490   END FUNCTION interp1 
     1491 
     1492   FUNCTION interp2(x, a, b, c, d)  RESULT(f)  
     1493      !!---------------------------------------------------------------------- 
     1494      !!                 ***  ROUTINE interp1  *** 
     1495      !!        
     1496      !! ** Purpose :   1-d constrained cubic spline interpolation 
     1497      !!           
     1498      !! ** Method  :  cubic spline interpolation 
     1499      !! 
     1500      !!---------------------------------------------------------------------- 
     1501      IMPLICIT NONE 
     1502      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1503      REAL(wp)             ::  f ! value from the interpolation 
     1504      !!---------------------------------------------------------------------- 
     1505 
     1506      f = a + x* ( b + x * ( c + d * x ) )  
     1507 
     1508   END FUNCTION interp2 
     1509 
     1510 
     1511   FUNCTION interp3(x, a, b, c, d)  RESULT(f)  
     1512      !!---------------------------------------------------------------------- 
     1513      !!                 ***  ROUTINE interp1  *** 
     1514      !!        
     1515      !! ** Purpose :   deriavtive of a cubic spline function 
     1516      !!           
     1517      !! ** Method  :   
     1518      !! 
     1519      !!---------------------------------------------------------------------- 
     1520      IMPLICIT NONE 
     1521      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1522      REAL(wp)             ::  f ! value from the interpolation 
     1523      !!---------------------------------------------------------------------- 
     1524 
     1525      f = b + x * ( 2._wp * c + 3._wp * d * x) 
     1526 
     1527   END FUNCTION interp3 
     1528 
     1529    
     1530   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f)  
     1531      !!---------------------------------------------------------------------- 
     1532      !!                 ***  ROUTINE interp1  *** 
     1533      !!        
     1534      !! ** Purpose :   1-d constrained cubic spline integration 
     1535      !!           
     1536      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr  
     1537      !! 
     1538      !!---------------------------------------------------------------------- 
     1539      IMPLICIT NONE 
     1540      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d    
     1541      REAL(wp)             ::  za1, za2, za3       
     1542      REAL(wp)             ::  f                   ! integration result 
     1543      !!---------------------------------------------------------------------- 
     1544 
     1545      za1 = 0.5_wp * b  
     1546      za2 = c / 3.0_wp  
     1547      za3 = 0.25_wp * d  
     1548 
     1549      f  = xr * ( a + xr * ( za1 + xr * ( za2 + za3 * xr ) ) ) - & 
     1550         & xl * ( a + xl * ( za1 + xl * ( za2 + za3 * xl ) ) ) 
     1551 
     1552   END FUNCTION integ2 
     1553 
     1554 
    10071555   !!====================================================================== 
    10081556END MODULE dynhpg 
Note: See TracChangeset for help on using the changeset viewer.