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 3068 for branches/2011/dev_r2802_NOCL_prjhpg – NEMO

Ignore:
Timestamp:
2011-11-09T14:29:55+01:00 (12 years ago)
Author:
hliu
Message:

updated pressure-jacobian horizontal pressure gradient code. only one file (dynhpg.F90) modified. 1) The code has been cleaned and restructured. 2)small modifications for hybrid s-sigma coordinate

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r2873 r3068  
    10501050      !! 
    10511051      !! The local varialbes for the correction term 
    1052       INTEGER  :: jke, jkw, jkn, jks, jk1 
     1052      INTEGER  :: jke, jkw, jkn, jks, jk1, jis, jid, jjs, jjd 
    10531053      REAL(wp) :: zze, zzw, zzn, zzs, rre, rrw, rrn, rrs 
    1054       REAL(wp) :: zuijk, zvijk, pe, pw, pn, ps, dze, dzw, dzn, dzs 
     1054      REAL(wp) :: zuijk, zvijk, pe, pw, pwes, pwed, pn, ps, pnss, pnsd, deps 
    10551055      REAL(wp) :: rhde, rhdw, rhdn, rhds,rhdt1, rhdt2 
    10561056      REAL(wp) :: dpdx1, dpdx2, dpdy1, dpdy2 
    … …  
    10811081      zhpi(:,:,:) = 0. 
    10821082       
    1083  
    10841083      !how to add vector opt.? N.B., jpi&jpi rather than jpim1&jpjm1 are needed here 
    10851084 
    … …  
    11181117      END DO 
    11191118 
    1120       ! Constrained cubic spline interpolation 
    1121  
     1119      !                 ! Construct the vertical density profile with the  
     1120      !                 !Constrained cubic spline interpolation 
    11221121      CALL cspline(fsp,xsp,asp,bsp,csp,dsp,polynomial_type)       
    11231122 
    1124  
    1125  
    1126       ! Calculate the pressure at T(ji,jj,1) 
     1123      ! Calculate the hydrostatic pressure at T(ji,jj,1) 
    11271124      DO jj = 2, jpj 
    11281125        DO ji = 2, jpi  
    … …  
    11301127                                         bsp(ji,jj,1),csp(ji,jj,1),dsp(ji,jj,1)) & 
    11311128                               * 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 
     1129          rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1130 
     1131     !                  ! assuming linear profile across the top half surface layer 
     1132          zhpi(ji,jj,1) =  0.5_wp * fse3w(ji,jj,1) * rhdt1   
    11341133        END DO 
    11351134      END DO 
    11361135 
    1137  
    1138  
    1139 !      print*, 'max&min fse3w=',maxval(fse3w), minval(fse3w) 
    1140 !      print*, 'max&min rhd===',maxval(rhd), minval(rhd) 
    1141              
    11421136      ! Calculate the pressure at T(ji,jj,2:jpkm1) 
    11431137      DO jk = 2, jpkm1                                   
    … …  
    11561150      DO jj = 2, jpjm1      
    11571151        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) 
     1152          zu(ji,jj,1) = - ( 0.5_wp * fse3uw(ji,jj,1) - sshu_n(ji,jj) * znad) 
     1153          zv(ji,jj,1) = - ( 0.5_wp * fse3vw(ji,jj,1) - sshv_n(ji,jj) * znad) 
    11601154        END DO 
    11611155      END DO 
    … …  
    11641158        DO jj = 2, jpjm1      
    11651159          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) 
     1160            zu(ji,jj,jk) = zu(ji,jj,jk-1)- fse3uw(ji,jj,jk) 
     1161            zv(ji,jj,jk) = zv(ji,jj,jk-1)- fse3vw(ji,jj,jk) 
    11681162          END DO 
    11691163        END DO 
    11701164      END DO 
    11711165                
     1166      !  Start pressure integration 
     1167 
    11721168      DO jk = 1, jpkm1                                   
    11731169        DO jj = 2, jpjm1      
    11741170          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     
     1171            pwes = 0._wp; pwed = 0._wp 
     1172            pnss = 0._wp; pnsd = 0._wp 
    11881173            zuijk = zu(ji,jj,jk) 
    11891174            zvijk = zv(ji,jj,jk) 
    11901175 
    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)) 
     1176            !!!!!     for u equation 
     1177            IF(-fsde3w(ji+1,jj,mbathy(ji+1,jj)) >= -fsde3w(ji,jj,mbathy(ji,jj))) THEN 
     1178              jis = ji + 1; jid = ji 
     1179            ELSE 
     1180              jis = ji;     jid = ji +1 
     1181            ENDIF 
     1182 
     1183            !     !integrate the pressure on the shallow side 
     1184            jk1 = jk  
     1185            bhitwe = 0 
     1186            IF(jk1 == mbathy(jis,jj)) THEN 
     1187              bhitwe = 1 
     1188            ELSE 
     1189              DO WHILE ( -fsde3w(jis,jj,jk1+1) > zuijk ) 
     1190                pwes = pwes + &  
     1191                     integ2(fsde3w(jis,jj,jk1),fsde3w(jis,jj,jk1+1),& 
     1192                            asp(jis,jj,jk1),bsp(jis,jj,jk1),& 
     1193                            csp(jis,jj,jk1),dsp(jis,jj,jk1)) 
     1194                jk1 = jk1 + 1 
     1195                IF(jk1 == mbathy(jis,jj)) THEN 
     1196                  bhitwe = 1 
    12061197                  EXIT 
    1207                 ENDIF 
     1198                END IF 
    12081199              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)) 
     1200            ENDIF 
     1201             
     1202            IF(bhitwe /= 1) THEN 
     1203              pwes = pwes + &  
     1204                   integ2(fsde3w(jis,jj,jk1),-zuijk,& 
     1205                          asp(jis,jj,jk1),bsp(jis,jj,jk1),& 
     1206                          csp(jis,jj,jk1),dsp(jis,jj,jk1)) 
     1207            ELSE 
     1208               zuijk = -fsde3w(jis,jj,jk1) 
     1209            ENDIF 
     1210 
     1211            !     !integrate the pressure on the deep side 
     1212            jk1 = jk  
     1213            bhitwe = 0 
     1214            IF(jk1 == 1) THEN 
     1215              bhitwe = 1 
     1216            ELSE 
     1217              DO WHILE ( -fsde3w(jid,jj,jk1-1) < zuijk ) 
     1218                pwed = pwed + &  
     1219                     integ2(fsde3w(jid,jj,jk1-1),fsde3w(jid,jj,jk1),& 
     1220                            asp(jid,jj,jk1-1),bsp(jid,jj,jk1-1),& 
     1221                            csp(jid,jj,jk1-1),dsp(jid,jj,jk1-1)) 
     1222                jk1 = jk1 - 1 
     1223                IF(jk1 == 1) THEN 
     1224                  bhitwe = 1 
    12351225                  EXIT 
    1236                 ENDIF 
     1226                END IF 
    12371227              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  
    13071228            ENDIF 
    1308  
     1229             
     1230            IF(bhitwe /= 1) THEN 
     1231              pwed = pwed + &  
     1232                   integ2(-zuijk,fsde3w(jid,jj,jk1),& 
     1233                          asp(jid,jj,jk1-1),bsp(jid,jj,jk1-1),& 
     1234                          csp(jid,jj,jk1-1),dsp(jid,jj,jk1-1)) 
     1235            ELSE 
     1236              deps = fsde3w(jid,jj,1) + min(zuijk, sshn(jid,jj)*znad) 
     1237              rhdt1 = rhd(jid,jj,1) - interp3(fsde3w(jid,jj,1),asp(jid,jj,1), & 
     1238                                      bsp(jid,jj,1),csp(jid,jj,1),dsp(jid,jj,1)) * deps 
     1239              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1240              pwed = pwed + 0.5_wp * (rhd(jid,jj,1) + rhdt1) * deps 
     1241            ENDIF 
     1242 
     1243            IF(jid > jis) THEN 
     1244              pe = pwed; pw = pwes 
     1245            ELSE 
     1246              pe = pwes; pw = pwed 
     1247            ENDIF 
    13091248 
    13101249            dpdx1 = zcoef0 / e1u(ji,jj) * (zhpi(ji+1,jj,jk) - zhpi(ji,jj,jk)) 
    … …  
    13171256            ua(ji,jj,jk) = ua(ji,jj,jk) + (dpdx1 + dpdx2)*& 
    13181257               &           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)) 
     1258   
     1259            !!!!!     for v equation 
     1260 
     1261            IF(-fsde3w(ji,jj+1,mbathy(ji,jj+1)) >= -fsde3w(ji,jj,mbathy(ji,jj))) THEN 
     1262              jjs = jj + 1; jjd = jj 
     1263            ELSE 
     1264              jjs = jj    ; jjd = jj + 1 
     1265            ENDIF 
     1266 
     1267            !     !integrate the pressure on the shallow side 
     1268            jk1 = jk  
     1269            bhitns = 0 
     1270            IF(jk1 == mbathy(ji,jjs)) THEN 
     1271              bhitns = 1 
     1272            ELSE 
     1273              DO WHILE ( -fsde3w(ji,jjs,jk1+1) > zvijk ) 
     1274                pnss = pnss + &  
     1275                     integ2(fsde3w(ji,jjs,jk1),fsde3w(ji,jjs,jk1+1),& 
     1276                            asp(ji,jjs,jk1),bsp(ji,jjs,jk1),& 
     1277                            csp(ji,jjs,jk1),dsp(ji,jjs,jk1)) 
     1278                jk1 = jk1 + 1 
     1279                IF(jk1 == mbathy(ji,jjs)) THEN 
     1280                  bhitns = 1 
    13361281                  EXIT 
    1337                 ENDIF 
     1282                END IF 
    13381283              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)) 
     1284            ENDIF 
     1285             
     1286            IF(bhitns /= 1) THEN 
     1287              pnss = pnss + &  
     1288                   integ2(fsde3w(ji,jjs,jk1),-zvijk,& 
     1289                          asp(ji,jjs,jk1),bsp(ji,jjs,jk1),& 
     1290                          csp(ji,jjs,jk1),dsp(ji,jjs,jk1)) 
     1291            ELSE 
     1292               zvijk = -fsde3w(ji,jjs,jk1) 
     1293            ENDIF 
     1294 
     1295            !     !integrate the pressure on the deep side 
     1296            jk1 = jk  
     1297            bhitns = 0 
     1298            IF(jk1 == 1) THEN 
     1299              bhitns = 1 
     1300            ELSE 
     1301              DO WHILE ( -fsde3w(ji,jjd,jk1-1) < zvijk ) 
     1302                pnsd = pnsd + &  
     1303                     integ2(fsde3w(ji,jjd,jk1-1),fsde3w(ji,jjd,jk1),& 
     1304                            asp(ji,jjd,jk1-1),bsp(ji,jjd,jk1-1),& 
     1305                            csp(ji,jjd,jk1-1),dsp(ji,jjd,jk1-1)) 
     1306                jk1 = jk1 - 1 
     1307                IF(jk1 == 1) THEN 
     1308                  bhitns = 1 
    13671309                  EXIT 
    1368                 ENDIF 
     1310                END IF 
    13691311              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  
     1312            ENDIF 
     1313             
     1314            IF(bhitns /= 1) THEN 
     1315              pnsd = pnsd + &  
     1316                   integ2(-zvijk,fsde3w(ji,jjd,jk1),& 
     1317                          asp(ji,jjd,jk1-1),bsp(ji,jjd,jk1-1),& 
     1318                          csp(ji,jjd,jk1-1),dsp(ji,jjd,jk1-1)) 
    13811319            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) 
     1320              deps = fsde3w(ji,jjd,1) + min(zvijk, sshn(ji,jjd)*znad) 
     1321              rhdt1 = rhd(ji,jjd,1) - interp3(fsde3w(ji,jjd,1),asp(ji,jjd,1), & 
     1322                                      bsp(ji,jjd,1),csp(ji,jjd,1),dsp(ji,jjd,1)) * deps 
     1323              rhdt1 = max(rhdt1, 1000._wp - rau0)        ! no lighter than fresh water 
     1324              pnsd = pnsd + 0.5_wp * (rhd(ji,jjd,1) + rhdt1) * deps 
     1325            ENDIF 
     1326 
     1327            IF(jjd > jjs) THEN 
     1328              pn = pnsd; ps = pnss 
     1329            ELSE 
     1330              pn = pnss; ps = pnsd 
     1331            ENDIF 
     1332 
     1333            dpdy1 = zcoef0 / e2v(ji,jj) * (zhpi(ji,jj+1,jk) - zhpi(ji,jj,jk)) 
     1334            if(lk_vvl) then 
     1335                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn + (sshn(ji,jj+1)-sshn(ji,jj)))  
     1336              else 
     1337                dpdy2 = zcoef0 / e2v(ji,jj) * (ps + pn )  
     1338            end if 
     1339 
     1340            va(ji,jj,jk) = va(ji,jj,jk) + (dpdy1 + dpdy2)*& 
     1341            &              vmask(ji,jj,jk)*tmask(ji,jj,jk)*tmask(ji,jj+1,jk) 
    14531342 
    14541343                     
    … …  
    14581347 
    14591348      ! Restore fsde3w and rhd 
    1460  
    14611349      fsde3w(:,:,:) = fsde3w_tmp(:,:,:) 
    14621350      rhd(:,:,:)    = rhd_tmp(:,:,:) 
    … …  
    14661354         CALL ctl_stop('dyn:hpg_prj: failed to release workspace arrays') 
    14671355      ! 
    1468  
    14691356   END SUBROUTINE hpg_prj 
    14701357 
Note: See TracChangeset for help on using the changeset viewer.