Changeset 2873


Ignore:
Timestamp:
2011-09-28T10:25:27+02:00 (9 years ago)
Author:
hliu
Message:

1). added pressure Jacobian horizontal pressure gradient code, 2) added par_AMM7/12.h90, 3) added two arch files for NOCL mobius and ubuntu linux

Location:
branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM
Files:
4 added
2 edited

Legend:

Unmodified
Added
Removed
  • branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/DYN/dynhpg.F90

    r2715 r2873  
    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, rn_gamma  , ln_dynhpg_imp 
    132135      !!---------------------------------------------------------------------- 
    133136      ! 
     
    147150         WRITE(numout,*) '      s-coord. (Density Jacobian: Cubic polynomial)     ln_hpg_djc    = ', ln_hpg_djc 
    148151         WRITE(numout,*) '      s-coord. (ROTated axes scheme)                    ln_hpg_rot    = ', ln_hpg_rot 
     152         WRITE(numout,*) '      s-coord. (Pressure Jacobian: Cubic polynomial)    ln_hpg_prj    = ', ln_hpg_prj 
    149153         WRITE(numout,*) '      weighting coeff. (wdj scheme)                     rn_gamma      = ', rn_gamma 
    150154         WRITE(numout,*) '      time stepping: centered (F) or semi-implicit (T)  ln_dynhpg_imp = ', ln_dynhpg_imp 
    151155      ENDIF 
    152156      ! 
    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') 
     157      IF( lk_vvl .AND. .NOT. (ln_hpg_sco.OR.ln_hpg_prj) )   & 
     158         &   CALL ctl_stop('dyn_hpg_init : variable volume key_vvl require:& 
     159                           & the standard jacobian formulation hpg_sco or & 
     160                           & the pressure jacobian formulation hpg_prj') 
    155161      ! 
    156162      !                               ! Set nhpg from ln_hpg_... flags 
     
    162168      IF( ln_hpg_djc )   nhpg = 5 
    163169      IF( ln_hpg_rot )   nhpg = 6 
     170      IF( ln_hpg_prj )   nhpg = 7 
    164171      ! 
    165172      !                               ! Consitency check 
     
    172179      IF( ln_hpg_djc )   ioptio = ioptio + 1 
    173180      IF( ln_hpg_rot )   ioptio = ioptio + 1 
     181      IF( ln_hpg_prj )   ioptio = ioptio + 1 
    174182      IF( ioptio /= 1 )   CALL ctl_stop( 'NO or several hydrostatic pressure gradient options used' ) 
    175183      ! 
     
    10051013   END SUBROUTINE hpg_rot 
    10061014 
     1015    
     1016   SUBROUTINE hpg_prj( kt ) 
     1017      !!--------------------------------------------------------------------- 
     1018      !!                  ***  ROUTINE hpg_prj  *** 
     1019      !! 
     1020      !! ** Method  :   s-coordinate case. 
     1021      !!      Reformulate the horizontal hydrostatical pressure gradient 
     1022      !!      term using Pressure Jacobian. A new correction term  
     1023      !!      is developed to eliminate the sigma-coordinate error. 
     1024      !! 
     1025      !! ** Action : - Update (ua,va) with the now hydrastatic pressure trend 
     1026      !!             - Save the trend (l_trddyn=T) 
     1027      !! 
     1028      !!---------------------------------------------------------------------- 
     1029 
     1030      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released 
     1031      USE oce     , ONLY:   fsde3w_tmp => ta       ! (ta,sa) used as 3D workspace 
     1032      USE oce     , ONLY:      rhd_tmp => sa  
     1033      USE wrk_nemo, ONLY:         zhpi => wrk_3d_3  
     1034      USE wrk_nemo, ONLY:           zu => wrk_3d_4  
     1035      USE wrk_nemo, ONLY:           zv => wrk_3d_5 
     1036      USE wrk_nemo, ONLY:          fsp => wrk_3d_6  
     1037      USE wrk_nemo, ONLY:          xsp => wrk_3d_7 
     1038      USE wrk_nemo, ONLY:          asp => wrk_3d_8 
     1039      USE wrk_nemo, ONLY:          bsp => wrk_3d_9 
     1040      USE wrk_nemo, ONLY:          csp => wrk_3d_10 
     1041      USE wrk_nemo, ONLY:          dsp => wrk_3d_11 
     1042      !! 
     1043      !!---------------------------------------------------------------------- 
     1044      !! 
     1045      INTEGER, PARAMETER  :: polynomial_type = 1    ! 1: cubic spline, 2: linear 
     1046      INTEGER, INTENT(in) ::   kt    ! ocean time-step index 
     1047      !! 
     1048      INTEGER  ::   ji, jj, jk, jkk                ! dummy loop indices 
     1049      REAL(wp) ::   zcoef0, znad               ! temporary scalars 
     1050      !! 
     1051      !! The local varialbes for the correction term 
     1052      INTEGER  :: jke, jkw, jkn, jks, jk1 
     1053      REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs 
     1054      REAL(wp) :: zuijk, zvijk, pe, pw, pn, ps, dze, dzw, dzn, dzs 
     1055      REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2 
     1056      REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 
     1057      INTEGER  :: bhitwe, bhitns 
     1058      !!---------------------------------------------------------------------- 
     1059 
     1060      IF( wrk_in_use(3, 3,4,5,6,7,8,9,10,11) ) THEN 
     1061         CALL ctl_stop('dyn:hpg_prj: requested workspace arrays unavailable')   ;   RETURN 
     1062      ENDIF 
     1063 
     1064      IF( kt == nit000 ) THEN 
     1065         IF(lwp) WRITE(numout,*) 
     1066         IF(lwp) WRITE(numout,*) 'dyn:hpg_prj : hydrostatic pressure gradient trend' 
     1067         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~   s-coordinate case, cubic spline pressure Jacobian' 
     1068      ENDIF 
     1069 
     1070      !!---------------------------------------------------------------------- 
     1071      ! Local constant initialization 
     1072      zcoef0 = - grav  
     1073      znad = 0.0_wp 
     1074      IF(lk_vvl) znad = 1._wp 
     1075 
     1076      ! Save fsde3w and rhd 
     1077      fsde3w_tmp(:,:,:) = fsde3w(:,:,:)  
     1078      rhd_tmp(:,:,:)    = rhd(:,:,:)    
     1079 
     1080      ! Clean 3-D work arraies 
     1081      zhpi(:,:,:) = 0. 
     1082       
     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      ! Constrained cubic spline interpolation 
     1121 
     1122      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
     1123 
     1124 
     1125 
     1126      ! Calculate the pressure at T(ji,jj,1) 
     1127      DO jj = 2, jpj 
     1128        DO ji = 2, jpi  
     1129          rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1130                                         bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1131                               * 0.5_wp * fsde3w(ji,jj,1) 
     1132          rhdt1 = max(rhdt1, 0.0_wp) 
     1133          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1 
     1134        END DO 
     1135      END DO 
     1136 
     1137 
     1138 
     1139!      print*, 'max&min fse3w=',maxval(fse3w), minval(fse3w) 
     1140!      print*, 'max&min rhd===',maxval(rhd), minval(rhd) 
     1141             
     1142      ! Calculate the pressure at T(ji,jj,2:jpkm1) 
     1143      DO jk = 2, jpkm1                                   
     1144        DO jj = 2, jpj      
     1145          DO ji = 2, jpi 
     1146            zhpi(ji,jj,jk) = zhpi(ji,jj,jk-1) + & 
     1147                             integ2(fsde3w(ji,jj,jk-1),fsde3w(ji,jj,jk),& 
     1148                                    asp(ji,jj,jk-1),bsp(ji,jj,jk-1),& 
     1149                                    csp(ji,jj,jk-1),dsp(ji,jj,jk-1)) 
     1150          END DO 
     1151        END DO 
     1152      END DO 
     1153 
     1154      ! Z coordinate of U(ji,jj,1:jpkm1) and V(ji,jj,1:jpkm1) 
     1155 
     1156      DO jj = 2, jpjm1      
     1157        DO ji = 2, jpim1   
     1158          zu(ji,jj,1) = - ( fse3u(ji,jj,1) - sshu_n(ji,jj) * znad) 
     1159          zv(ji,jj,1) = - ( fse3v(ji,jj,1) - sshv_n(ji,jj) * znad) 
     1160        END DO 
     1161      END DO 
     1162 
     1163      DO jk = 2, jpkm1                                   
     1164        DO jj = 2, jpjm1      
     1165          DO ji = 2, jpim1   
     1166            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3u(ji,jj,jk) 
     1167            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3v(ji,jj,jk) 
     1168          END DO 
     1169        END DO 
     1170      END DO 
     1171                
     1172      DO jk = 1, jpkm1                                   
     1173        DO jj = 2, jpjm1      
     1174          DO ji = 2, jpim1   
     1175            zu(ji,jj,jk) = zu(ji,jj,jk) + 0.5_wp * fse3u(ji,jj,jk) 
     1176            zv(ji,jj,jk) = zv(ji,jj,jk) + 0.5_wp * fse3v(ji,jj,jk) 
     1177          END DO 
     1178        END DO 
     1179      END DO 
     1180 
     1181      DO jk = 1, jpkm1                                   
     1182        DO jj = 2, jpjm1      
     1183          DO ji = 2, jpim1   
     1184                
     1185            pe = 0.; pw = 0.; dze = 0.; dzw = 0. 
     1186            pn = 0.; ps = 0.; dzn = 0.; dzs = 0. 
     1187            bhitwe = 0;      bhitns = 0     
     1188            zuijk = zu(ji,jj,jk) 
     1189            zvijk = zv(ji,jj,jk) 
     1190 
     1191            !! for u equation 
     1192            IF(-fsde3w(ji+1,jj,jk) > zuijk) THEN 
     1193              DO jk1 = jk+1, mbathy(ji+1,jj) 
     1194                IF(-fsde3w(ji+1,jj,jk1) > zuijk) THEN 
     1195                  pe = pe + &  
     1196                       integ2(fsde3w(ji+1,jj,jk1-1),fsde3w(ji+1,jj,jk1),& 
     1197                              asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),& 
     1198                              csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1)) 
     1199                  IF(jk1 == mbathy(ji+1,jj)) bhitwe = 1  
     1200                 ELSE 
     1201                  zze = -fsde3w(ji+1,jj,jk1-1) - zuijk 
     1202                  pe = pe + &  
     1203                  integ2(fsde3w(ji+1,jj,jk1-1),-zuijk,& 
     1204                         asp(ji+1,jj,jk1-1),bsp(ji+1,jj,jk1-1),& 
     1205                         csp(ji+1,jj,jk1-1),dsp(ji+1,jj,jk1-1)) 
     1206                  EXIT 
     1207                ENDIF 
     1208              END DO 
     1209 
     1210              IF(jk == mbathy(ji+1,jj))  bhitwe = 1  
     1211 
     1212              IF(bhitwe == 1) zuijk = -fsde3w(ji+1,jj,mbathy(ji+1,jj))   
     1213 
     1214              DO jk1 = jk-1, 1, -1 
     1215                IF(-fsde3w(ji,jj,jk1) < zuijk) THEN 
     1216                  pw = pw + & 
     1217                           integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),& 
     1218                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),& 
     1219                                   csp(ji,jj,jk1),dsp(ji,jj,jk1)) 
     1220                  IF(jk1 == 1) THEN 
     1221                    rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1222                                           bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1223                                         * (zuijk + fsde3w(ji,jj,1)) 
     1224                    rhdt1 = max(rhdt1,0.0_wp) 
     1225                    rrw = rhdt1 
     1226                    zzw = zuijk + fsde3w(ji,jj,1) 
     1227                    pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1)) 
     1228                  END IF 
     1229                 ELSE 
     1230                  zzw = zuijk + fsde3w(ji,jj,jk1+1) 
     1231                  pw = pw + & 
     1232                           integ2(-zuijk,fsde3w(ji,jj,jk1+1),& 
     1233                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),& 
     1234                                   csp(ji,jj,jk1),dsp(ji,jj,jk1)) 
     1235                  EXIT 
     1236                ENDIF 
     1237              END DO 
     1238 
     1239              IF(jk == 1) THEN 
     1240                rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1241                                       bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1242                                     * (zuijk + fsde3w(ji,jj,1)) 
     1243                rhdt1 = max(rhdt1,0.0_wp) 
     1244                rrw = rhdt1 
     1245                zzw = zuijk + fsde3w(ji,jj,1) 
     1246                pw = pw + min(0.5 * (rhd(ji,jj,1) + rrw) * zzw, zhpi(ji,jj,1)) 
     1247              END IF 
     1248 
     1249              ELSE 
     1250 
     1251              DO jk1 = jk+1, mbathy(ji,jj) 
     1252                IF(-fsde3w(ji,jj,jk1) > zuijk) THEN 
     1253                  pw = pw - & 
     1254                          integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),& 
     1255                                  asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 
     1256                                  csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 
     1257                  IF(jk1 == mbathy(ji,jj)) bhitwe = 1  
     1258                ELSE 
     1259                  zzw = -fsde3w(ji,jj,jk1-1) - zuijk 
     1260                  pw = pw - & 
     1261                           integ2(fsde3w(ji,jj,jk1-1),-zuijk,& 
     1262                                  asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 
     1263                                  csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 
     1264                  EXIT 
     1265                ENDIF 
     1266              END DO 
     1267 
     1268              IF(jk == mbathy(ji,jj)) bhitwe = 1  
     1269 
     1270              IF(bhitwe == 1) zuijk = -fsde3w(ji,jj,mbathy(ji,jj))  
     1271 
     1272              DO jk1 = jk-1, 1, -1 
     1273                IF(-fsde3w(ji+1,jj,jk1) < zuijk) THEN 
     1274                  pe = pe - & 
     1275                           integ2(fsde3w(ji+1,jj,jk1),fsde3w(ji+1,jj,jk1+1),& 
     1276                                   asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),& 
     1277                                   csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1)) 
     1278                  IF(jk1 == 1) THEN 
     1279                    rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), & 
     1280                                           bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) & 
     1281                                         * (zuijk + fsde3w(ji+1,jj,1)) 
     1282                    rhdt1 = max(rhdt1,0.0_wp) 
     1283                    rre = rhdt1 
     1284                    zze = zuijk + fsde3w(ji+1,jj,1) 
     1285                    pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1)) 
     1286                  END IF 
     1287                 ELSE 
     1288                  zze = zuijk + fsde3w(ji+1,jj,jk1+1) 
     1289                  pe = pe - & 
     1290                           integ2(-zuijk,fsde3w(ji+1,jj,jk1+1),& 
     1291                                   asp(ji+1,jj,jk1),bsp(ji+1,jj,jk1),& 
     1292                                   csp(ji+1,jj,jk1),dsp(ji+1,jj,jk1)) 
     1293                  EXIT 
     1294                ENDIF 
     1295              END DO 
     1296 
     1297              IF(jk == 1) THEN 
     1298                 rhdt1 = rhd(ji+1,jj,1) - interp3(fsde3w(ji+1,jj,1),asp(ji+1,jj,1), & 
     1299                                          bsp(ji+1,jj,1),csp(ji+1,jj,1),dsp(ji+1,jj,1)) & 
     1300                                        * (zuijk + fsde3w(ji+1,jj,1)) 
     1301                rhdt1 = max(rhdt1,0.0_wp) 
     1302                rre = rhdt1 
     1303                zze = zuijk + fsde3w(ji+1,jj,1) 
     1304                pe = pe - min(0.5 * (rhd(ji+1,jj,1) + rre) * zze, zhpi(ji+1,jj,1)) 
     1305              END IF 
     1306 
     1307            ENDIF 
     1308 
     1309 
     1310            dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
     1311            IF(lk_vvl) THEN 
     1312              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe + (sshn(ji+1,jj)-sshn(ji,jj)))  
     1313             ELSE 
     1314              dpdx2 = zcoef0 / e1u(ji,jj) * (pw + pe)  
     1315            ENDIF 
     1316 
     1317            ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2)*& 
     1318               &           umask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji+1,jj,jk) 
     1319                 
     1320            !! for v equation 
     1321 
     1322            IF(-fsde3w(ji,jj+1,jk) > zvijk) THEN 
     1323              DO jk1 = jk+1, mbathy(ji,jj+1) 
     1324                IF(-fsde3w(ji,jj+1,jk1) > zvijk) THEN 
     1325                  pn = pn + & 
     1326                           integ2(fsde3w(ji,jj+1,jk1-1),fsde3w(ji,jj+1,jk1),& 
     1327                                   asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),& 
     1328                                   csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1)) 
     1329                  IF(jk1 == mbathy(ji,jj+1)) bhitns = 1  
     1330                 ELSE 
     1331                  zzn = -fsde3w(ji,jj+1,jk1-1) - zvijk 
     1332                  pn = pn + & 
     1333                           integ2(fsde3w(ji,jj+1,jk1-1),-zvijk,& 
     1334                                   asp(ji,jj+1,jk1-1),bsp(ji,jj+1,jk1-1),& 
     1335                                   csp(ji,jj+1,jk1-1),dsp(ji,jj+1,jk1-1)) 
     1336                  EXIT 
     1337                ENDIF 
     1338              END DO 
     1339 
     1340              IF(jk == mbathy(ji,jj+1)) bhitns = 1 
     1341 
     1342 
     1343              IF(bhitns == 1) zvijk = -fsde3w(ji,jj+1,mbathy(ji,jj+1))  
     1344 
     1345 
     1346              DO jk1 = jk-1, 1, -1 
     1347                IF(-fsde3w(ji,jj,jk1) < zvijk) THEN 
     1348                  ps = ps + & 
     1349                           integ2(fsde3w(ji,jj,jk1),fsde3w(ji,jj,jk1+1),& 
     1350                                  asp(ji,jj,jk1),bsp(ji,jj,jk1),& 
     1351                                  csp(ji,jj,jk1),dsp(ji,jj,jk1)) 
     1352                  IF(jk1 == 1) THEN 
     1353                    rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1354                                           bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1355                                         * (zvijk + fsde3w(ji,jj,1)) 
     1356                    rhdt1 = max(rhdt1,0.0_wp) 
     1357                    rrs = rhdt1 
     1358                    zzs = zvijk + fsde3w(ji,jj,1) 
     1359                    ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs ) * zzs, zhpi(ji,jj,1)) 
     1360                  END IF 
     1361                ELSE 
     1362                  zzs = zvijk + fsde3w(ji,jj,jk1+1) 
     1363                  ps = ps + & 
     1364                           integ2(-zvijk,fsde3w(ji,jj,jk1+1),& 
     1365                                   asp(ji,jj,jk1),bsp(ji,jj,jk1),& 
     1366                                   csp(ji,jj,jk1),dsp(ji,jj,jk1)) 
     1367                  EXIT 
     1368                ENDIF 
     1369              END DO 
     1370 
     1371              IF(jk == 1) THEN 
     1372                rhdt1 = rhd(ji,jj,1) - interp3(fsde3w(ji,jj,1),asp(ji,jj,1), & 
     1373                                               bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
     1374                                     * (zvijk + fsde3w(ji,jj,1)) 
     1375                rhdt1 = max(rhdt1,0.0_wp) 
     1376                rrs = rhdt1 
     1377                zzs = zvijk + fsde3w(ji,jj,1) 
     1378                ps = ps + min(0.5 * (rhd(ji,jj,1) + rrs) * zzs, zhpi(ji,jj,1)) 
     1379              END IF 
     1380 
     1381            ELSE 
     1382 
     1383              DO jk1 = jk+1, mbathy(ji,jj) 
     1384                IF(-fsde3w(ji,jj,jk1) > zvijk) THEN 
     1385                  ps = ps - & 
     1386                           integ2(fsde3w(ji,jj,jk1-1),fsde3w(ji,jj,jk1),& 
     1387                                   asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 
     1388                                   csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 
     1389                  IF(jk1 == mbathy(ji,jj)) bhitns = 1 
     1390                 ELSE 
     1391                  zzs = -fsde3w(ji,jj,jk1-1) - zvijk 
     1392                  ps = ps - & 
     1393                           integ2(fsde3w(ji,jj,jk1-1),-zvijk,& 
     1394                                   asp(ji,jj,jk1-1),bsp(ji,jj,jk1-1),& 
     1395                                   csp(ji,jj,jk1-1),dsp(ji,jj,jk1-1)) 
     1396                  EXIT 
     1397                ENDIF 
     1398              END DO 
     1399 
     1400              IF(jk == mbathy(ji,jj)) bhitns = 1  
     1401 
     1402 
     1403              IF(bhitns == 1) zvijk = -fsde3w(ji,jj,mbathy(ji,jj)) 
     1404 
     1405                DO jk1 = jk-1, 1, -1 
     1406                  IF(-fsde3w(ji,jj+1,jk1) < zvijk) THEN 
     1407                    pn = pn - & 
     1408                             integ2(fsde3w(ji,jj+1,jk1),fsde3w(ji,jj+1,jk1+1),& 
     1409                                    asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),& 
     1410                                    csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1)) 
     1411                    IF(jk1 == 1) THEN 
     1412                      rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), & 
     1413                                             bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) & 
     1414                                           * (zvijk + fsde3w(ji,jj+1,1)) 
     1415                      rhdt1 = max(rhdt1,0.0_wp) 
     1416                      rrn = rhdt1 
     1417                      zzn = zvijk + fsde3w(ji,jj+1,1) 
     1418                      pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1)) 
     1419                    END IF 
     1420 
     1421                   ELSE 
     1422                    zzn = zvijk + fsde3w(ji,jj+1,jk1+1) 
     1423                    pn = pn - & 
     1424                             integ2(-zvijk,fsde3w(ji,jj+1,jk1+1),& 
     1425                                     asp(ji,jj+1,jk1),bsp(ji,jj+1,jk1),& 
     1426                                     csp(ji,jj+1,jk1),dsp(ji,jj+1,jk1)) 
     1427                    EXIT 
     1428                  ENDIF 
     1429                END DO 
     1430 
     1431                IF(jk == 1) THEN 
     1432                  rhdt1 = rhd(ji,jj+1,1) - interp3(fsde3w(ji,jj+1,1),asp(ji,jj+1,1), & 
     1433                                          bsp(ji,jj+1,1),csp(ji,jj+1,1),dsp(ji,jj+1,1)) & 
     1434                                        * (zvijk + fsde3w(ji,jj+1,1)) 
     1435                  rhdt1 = max(rhdt1,0.0_wp) 
     1436                  rrn = rhdt1 
     1437                  zzn = zvijk + fsde3w(ji,jj+1,1) 
     1438                  pn = pn - min(0.5 * (rhd(ji,jj+1,1) + rrn) * zzn, zhpi(ji,jj+1,1)) 
     1439                END IF 
     1440 
     1441              ENDIF 
     1442 
     1443 
     1444              dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     1445              if(lk_vvl) then 
     1446                  dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj)))  
     1447                else 
     1448                  dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn )  
     1449              end if 
     1450 
     1451              va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*& 
     1452              &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
     1453 
     1454                     
     1455           END DO 
     1456        END DO 
     1457      END DO 
     1458 
     1459      ! Restore fsde3w and rhd 
     1460 
     1461      fsde3w(:,:,:) = fsde3w_tmp(:,:,:) 
     1462      rhd(:,:,:)    = rhd_tmp(:,:,:) 
     1463 
     1464      ! 
     1465      IF( wrk_not_released(3, 3,4,5,6,7,8,9,10,11) )   & 
     1466         CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 
     1467      ! 
     1468 
     1469   END SUBROUTINE hpg_prj 
     1470 
     1471   SUBROUTINE cspline(fsp, xsp, asp, bsp, csp, dsp, polynomial_type) 
     1472      !!---------------------------------------------------------------------- 
     1473      !!                 ***  ROUTINE cspline  *** 
     1474      !!        
     1475      !! ** Purpose :   constrained cubic spline interpolation 
     1476      !!           
     1477      !! ** Method  :   f(x) = asp + bsp*x + csp*x^2 + dsp*x^3  
     1478      !! Reference: K.W. Brodlie, A review of mehtods for curve and function 
     1479      !!                          drawing, 1980 
     1480      !! 
     1481      !!---------------------------------------------------------------------- 
     1482      IMPLICIT NONE 
     1483      REAL(wp), DIMENSION(:,:,:), INTENT(in)  :: fsp, xsp           ! value and coordinate 
     1484      REAL(wp), DIMENSION(:,:,:), INTENT(out) :: asp, bsp, csp, dsp ! coefficients of  
     1485                                                                    ! the interpoated function 
     1486      INTEGER, INTENT(in) :: polynomial_type                        ! 1: cubic spline  
     1487                                                                    ! 2: Linear 
     1488 
     1489      ! Local Variables       
     1490      INTEGER  ::   ji, jj, jk                 ! dummy loop indices 
     1491      INTEGER  ::   jpi, jpj, jpkm1 
     1492      REAL(wp) :: df1, df2, ddf1, ddf2, tmp1, tmp2, dxtmp 
     1493      REAL(wp) :: dxtmp1, dxtmp2, alpha 
     1494      REAL(wp) :: df(size(fsp,3)) 
     1495      !!---------------------------------------------------------------------- 
     1496 
     1497      jpi   = size(fsp,1) 
     1498      jpj   = size(fsp,2) 
     1499      jpkm1 = size(fsp,3) - 1 
     1500 
     1501      ! Clean output arrays 
     1502      asp = 0.0_wp 
     1503      bsp = 0.0_wp 
     1504      csp = 0.0_wp 
     1505      dsp = 0.0_wp 
     1506       
     1507       
     1508      Do ji = 1, jpi 
     1509      Do jj = 1, jpj 
     1510 
     1511      If (polynomial_type == 1) Then     !Constrained Cubic Spline 
     1512          Do jk = 2, jpkm1-1 
     1513             dxtmp1 = xsp(ji,jj,jk)   - xsp(ji,jj,jk-1)   
     1514             dxtmp2 = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)   
     1515             df1    = ( fsp(ji,jj,jk)   - fsp(ji,jj,jk-1) ) / dxtmp1 
     1516             df2    = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk)   ) / dxtmp2 
     1517   
     1518             alpha = ( dxtmp1 + 2._wp * dxtmp2 ) / ( dxtmp1 + dxtmp2 ) / 3._wp 
     1519              
     1520             If(df1 * df2 <= 0._wp) Then 
     1521               df(jk) = 0._wp 
     1522             Else 
     1523               df(jk) = df1 * df2 / ( ( 1._wp - alpha ) * df1 + alpha * df2 ) 
     1524             End If 
     1525          End Do 
     1526   
     1527          df(1)     = 1.5_wp * ( fsp(ji,jj,2) - fsp(ji,jj,1) ) / ( xsp(ji,jj,2) - xsp(ji,jj,1) ) - & 
     1528                    & 0.5_wp * df(2) 
     1529          df(jpkm1) = 1.5_wp * ( fsp(ji,jj,jpkm1) - fsp(ji,jj,jpkm1-1) ) / & 
     1530                    &          ( xsp(ji,jj,jpkm1) - xsp(ji,jj,jpkm1-1) ) - & 
     1531                    & 0.5_wp * df(jpkm1 - 1) 
     1532   
     1533          Do jk = 1, jpkm1 - 1 
     1534             dxtmp = xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1535             tmp1  = (df(jk+1) + 2._wp * df(jk)) / dxtmp 
     1536             tmp2  =  6._wp * (fsp(ji,jj,jk+1) - fsp(ji,jj,jk)) / dxtmp / dxtmp 
     1537             ddf1  = -2._wp * tmp1 + tmp2  
     1538             tmp1  = (2._wp * df(jk+1) + df(jk)) / dxtmp 
     1539             ddf2  =  2._wp * tmp1 - tmp2  
     1540   
     1541             dsp(ji,jj,jk) = (ddf2 - ddf1) / 6._wp / dxtmp 
     1542             csp(ji,jj,jk) = ( xsp(ji,jj,jk+1) * ddf1 - xsp(ji,jj,jk)*ddf2 ) / 2._wp / dxtmp 
     1543             bsp(ji,jj,jk) = ( fsp(ji,jj,jk+1) - fsp(ji,jj,jk) ) / dxtmp - &  
     1544                           & csp(ji,jj,jk) * ( xsp(ji,jj,jk+1) + xsp(ji,jj,jk) ) - & 
     1545                           & dsp(ji,jj,jk) * ( xsp(ji,jj,jk+1)**2 + & 
     1546                           &                   xsp(ji,jj,jk+1) * xsp(ji,jj,jk) + & 
     1547                           &                   xsp(ji,jj,jk)**2 ) 
     1548             asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) - & 
     1549                           &                 csp(ji,jj,jk) * xsp(ji,jj,jk)**2 - & 
     1550                           &                 dsp(ji,jj,jk) * xsp(ji,jj,jk)**3 
     1551          End Do 
     1552 
     1553       Else If (polynomial_type == 2) Then     !Linear 
     1554 
     1555        Do jk = 1, jpkm1-1 
     1556           dxtmp =xsp(ji,jj,jk+1) - xsp(ji,jj,jk)  
     1557           tmp1 = fsp(ji,jj,jk+1) - fsp(ji,jj,jk) 
     1558   
     1559           dsp(ji,jj,jk) = 0._wp 
     1560           csp(ji,jj,jk) = 0._wp 
     1561           bsp(ji,jj,jk) = tmp1 / dxtmp 
     1562           asp(ji,jj,jk) = fsp(ji,jj,jk) - bsp(ji,jj,jk) * xsp(ji,jj,jk) 
     1563        End Do 
     1564 
     1565       Else 
     1566        CALL ctl_stop( 'invalid polynomial type in cspline' ) 
     1567      End If 
     1568 
     1569      End Do 
     1570      End Do 
     1571       
     1572   End Subroutine cspline 
     1573 
     1574 
     1575   FUNCTION interp1(x, xl, xr, fl, fr)  RESULT(f)  
     1576      !!---------------------------------------------------------------------- 
     1577      !!                 ***  ROUTINE interp1  *** 
     1578      !!        
     1579      !! ** Purpose :   1-d linear interpolation 
     1580      !!           
     1581      !! ** Method  :   
     1582      !!                interpolation is straigt forward 
     1583      !!                extrapolation is also permitted (no value limit)  
     1584      !! 
     1585      !! H.Liu, Jan 2009,  POL  
     1586      !!---------------------------------------------------------------------- 
     1587      IMPLICIT NONE 
     1588      REAL(wp), INTENT(in) ::  x, xl, xr, fl, fr    
     1589      REAL(wp)             ::  f ! result of the interpolation (extrapolation) 
     1590      REAL(wp)             ::  deltx 
     1591      !!---------------------------------------------------------------------- 
     1592 
     1593      deltx = xr - xl 
     1594      IF(abs(deltx) <= 10._wp * EPSILON(x)) THEN 
     1595        f = 0.5_wp * (fl + fr) 
     1596       ELSE 
     1597        f = ( (x - xl ) * fr - ( x - xr ) * fl ) / deltx 
     1598      END IF 
     1599       
     1600   END FUNCTION interp1 
     1601 
     1602   FUNCTION interp2(x, a, b, c, d)  RESULT(f)  
     1603      !!---------------------------------------------------------------------- 
     1604      !!                 ***  ROUTINE interp1  *** 
     1605      !!        
     1606      !! ** Purpose :   1-d constrained cubic spline interpolation 
     1607      !!           
     1608      !! ** Method  :  cubic spline interpolation 
     1609      !! 
     1610      !!---------------------------------------------------------------------- 
     1611      IMPLICIT NONE 
     1612      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1613      REAL(wp)             ::  f ! value from the interpolation 
     1614      !!---------------------------------------------------------------------- 
     1615 
     1616      f = a + x* ( b + x * ( c + d * x ) )  
     1617 
     1618   END FUNCTION interp2 
     1619 
     1620 
     1621   FUNCTION interp3(x, a, b, c, d)  RESULT(f)  
     1622      !!---------------------------------------------------------------------- 
     1623      !!                 ***  ROUTINE interp1  *** 
     1624      !!        
     1625      !! ** Purpose :   deriavtive of a cubic spline function 
     1626      !!           
     1627      !! ** Method  :   
     1628      !! 
     1629      !!---------------------------------------------------------------------- 
     1630      IMPLICIT NONE 
     1631      REAL(wp), INTENT(in) ::  x, a, b, c, d    
     1632      REAL(wp)             ::  f ! value from the interpolation 
     1633      !!---------------------------------------------------------------------- 
     1634 
     1635      f = b + x * ( 2._wp * c + 3._wp * d * x) 
     1636 
     1637   END FUNCTION interp3 
     1638 
     1639    
     1640   FUNCTION integ2(xl, xr, a, b, c, d)  RESULT(f)  
     1641      !!---------------------------------------------------------------------- 
     1642      !!                 ***  ROUTINE interp1  *** 
     1643      !!        
     1644      !! ** Purpose :   1-d constrained cubic spline integration 
     1645      !!           
     1646      !! ** Method  :  integrate polynomial a+bx+cx^2+dx^3 from xl to xr  
     1647      !! 
     1648      !!---------------------------------------------------------------------- 
     1649      IMPLICIT NONE 
     1650      REAL(wp), INTENT(in) ::  xl, xr, a, b, c, d    
     1651      REAL(wp)             ::  a1, a2,a3       
     1652      REAL(wp)             ::  f                   ! integration result 
     1653      !!---------------------------------------------------------------------- 
     1654 
     1655      a1 = 0.5_wp * b  
     1656      a2 = c / 3.0_wp  
     1657      a3 = 0.25_wp * d  
     1658 
     1659      f  = xr * ( a + xr * ( a1 + xr * ( a2 + a3 * xr ) ) ) - & 
     1660         & xl * ( a + xl * ( a1 + xl * ( a2 + a3 * xl ) ) ) 
     1661 
     1662   END FUNCTION integ2 
     1663 
     1664 
    10071665   !!====================================================================== 
    10081666END MODULE dynhpg 
  • branches/2011/dev_r2802_NOCL_prjhpg/NEMOGCM/NEMO/OPA_SRC/par_oce.F90

    r2715 r2873  
    8181   !!--------------------------------------------------------------------- 
    8282#             include "par_POMME_R025.h90" 
     83#elif defined key_amm12  
     84   !!---------------------------------------------------------------------  
     85   !!   'key_amm12'         :   Atlantic Margin Model (~12km)     :      AMM  
     86   !!---------------------------------------------------------------------  
     87#             include "par_AMM12.h90"  
     88#elif defined key_amm7  
     89   !!---------------------------------------------------------------------  
     90   !!   'key_amm7'         :   Atlantic Margin Model (~7km)     :      AMM  
     91   !!---------------------------------------------------------------------  
     92#             include "par_AMM7.h90"  
    8393#else 
    8494   !!--------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.