- Timestamp:
- 2018-07-11T10:24:17+02:00 (6 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2018/dev_r9838_ENHANCE04_MLF/src/OCE/DYN/dynspg_ts.F90
r9863 r9923 69 69 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 70 70 ! 71 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_ baro <= 2.5 nn_baro72 REAL(wp),SAVE :: r dtbt ! Barotropic timestep71 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5*nn_e 72 REAL(wp),SAVE :: rDt_e ! external mode time-step 73 73 ! 74 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 81 81 REAL(wp) :: r1_4 = 0.25_wp ! 82 82 REAL(wp) :: r1_2 = 0.5_wp ! 83 84 REAL(wp) :: r1_2dt_b, r2dt_bf ! local scalars85 83 86 84 !! * Substitutions … … 101 99 ierr(:) = 0 102 100 ! 103 ALLOCATE( wgtbtp1(3*nn_ baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) )101 ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) ) 104 102 ! 105 103 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & … … 150 148 INTEGER :: ikbu, iktu, noffset ! local integers 151 149 INTEGER :: ikbv, iktv ! - - 150 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 152 151 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 153 152 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 154 153 REAL(wp) :: za0, za1, za2, za3 ! - - 155 154 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 155 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 156 REAL(wp) :: zload 156 157 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 157 158 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc … … 161 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 162 163 ! 163 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. 164 165 INTEGER :: iwdg, jwdg, kwdg ! short-hand values for the indices of the output point 164 166 165 167 166 REAL(wp) :: zepsilon, zgamma ! - - … … 179 178 zwdramp = r_rn_wdmin1 ! simplest ramp 180 179 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 181 ! ! reciprocal of baroclinic time step182 IF( l_1st_euler ) THEN ; r2dt_bf = rdt183 ELSE ; r2dt_bf = 2.0_wp * rdt184 ENDIF185 r1_2dt_b = 1.0_wp / r2dt_bf186 180 ! 187 181 ll_init = ln_bt_av ! if no time averaging, then no specific restart 188 182 ll_fw_start = .FALSE. 189 183 ! ! time offset in steps for bdy data update 190 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_ baro184 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 191 185 ELSE ; noffset = 0 192 186 ENDIF … … 459 453 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 460 454 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 461 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp)455 zcpx(ji,jj) = MAX( 0._wp , MIN( zcpx(ji,jj) , 1._wp ) ) 462 456 ELSE 463 457 zcpx(ji,jj) = 0._wp … … 536 530 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 537 531 IF( ln_wd_il ) THEN 538 zztmp = -1._wp / r dtbt532 zztmp = -1._wp / rDt_e 539 533 DO jj = 2, jpjm1 540 534 DO ji = fs_2, fs_jpim1 ! vector opt. … … 587 581 DO jj = 2, jpjm1 588 582 DO ji = fs_2, fs_jpim1 ! vector opt. 589 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_r au0 * utau(ji,jj) * r1_hu_n(ji,jj)590 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_r au0 * vtau(ji,jj) * r1_hv_n(ji,jj)583 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu_n(ji,jj) 584 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv_n(ji,jj) 591 585 END DO 592 586 END DO 593 587 ELSE 594 zztmp = r1_r au0 * r1_2588 zztmp = r1_rho0 * r1_2 595 589 DO jj = 2, jpjm1 596 590 DO ji = fs_2, fs_jpim1 ! vector opt. … … 629 623 ! ! Surface net water flux and rivers 630 624 IF (ln_bt_fw) THEN 631 zssh_frc(:,:) = r1_r au0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )625 zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 632 626 ELSE 633 zztmp = r1_r au0 * r1_2627 zztmp = r1_rho0 * r1_2 634 628 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 635 629 & + fwfisf(:,:) + fwfisf_b(:,:) ) … … 818 812 ENDIF 819 813 #endif 820 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, r dtbt)814 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rDt_e) 821 815 822 816 IF ( ln_wd_dl ) THEN … … 864 858 END DO 865 859 END DO 866 ssha_e(:,:) = ( sshn_e(:,:) - r dtbt* ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:)860 ssha_e(:,:) = ( sshn_e(:,:) - rDt_e * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 867 861 868 862 CALL lbc_lnk( ssha_e, 'T', 1._wp ) … … 1068 1062 ENDIF 1069 1063 ! 1070 ! Surface pressure trend: 1064 ! Surface pressure trend 1065 IF( ln_scal_load ) THEN ; zload = 1._wp 1066 ELSE ; zload = 1._wp - rn_load 1067 ENDIF 1071 1068 IF( ln_wd_il ) THEN 1072 1069 DO jj = 2, jpjm1 … … 1075 1072 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1076 1073 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1077 zwx(ji,jj) = (1._wp - rn_scal_load)* zu_spg * zcpx(ji,jj)1078 zwy(ji,jj) = (1._wp - rn_scal_load)* zv_spg * zcpy(ji,jj)1074 zwx(ji,jj) = zload * zu_spg * zcpx(ji,jj) 1075 zwy(ji,jj) = zload * zv_spg * zcpy(ji,jj) 1079 1076 END DO 1080 1077 END DO … … 1085 1082 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1086 1083 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1087 zwx(ji,jj) = (1._wp - rn_scal_load)* zu_spg1088 zwy(ji,jj) = (1._wp - rn_scal_load)* zv_spg1084 zwx(ji,jj) = zload * zu_spg 1085 zwy(ji,jj) = zload * zv_spg 1089 1086 END DO 1090 1087 END DO … … 1097 1094 DO ji = fs_2, fs_jpim1 ! vector opt. 1098 1095 ua_e(ji,jj) = ( un_e(ji,jj) & 1099 & + r dtbt* ( zwx(ji,jj) &1096 & + rDt_e * ( zwx(ji,jj) & 1100 1097 & + zu_trd(ji,jj) & 1101 1098 & + zu_frc(ji,jj) ) & … … 1103 1100 1104 1101 va_e(ji,jj) = ( vn_e(ji,jj) & 1105 & + r dtbt* ( zwy(ji,jj) &1102 & + rDt_e * ( zwy(ji,jj) & 1106 1103 & + zv_trd(ji,jj) & 1107 1104 & + zv_frc(ji,jj) ) & … … 1110 1107 !jth implicit bottom friction: 1111 1108 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1112 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - r dtbt* zCdU_u(ji,jj) * hur_e(ji,jj))1113 va_e(ji,jj) = va_e(ji,jj) /(1.0 - r dtbt* zCdU_v(ji,jj) * hvr_e(ji,jj))1109 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 1110 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 1114 1111 ENDIF 1115 1112 … … 1128 1125 1129 1126 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1130 & + r dtbt* ( zhust_e(ji,jj) * zwx(ji,jj) &1127 & + rDt_e * ( zhust_e(ji,jj) * zwx(ji,jj) & 1131 1128 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1132 1129 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & … … 1134 1131 1135 1132 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1136 & + r dtbt* ( zhvst_e(ji,jj) * zwy(ji,jj) &1133 & + rDt_e * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1137 1134 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1138 1135 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & … … 1202 1199 zwy(:,:) = vn_adv(:,:) 1203 1200 IF( .NOT.l_1st_euler ) THEN 1204 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) )1205 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) )1201 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - rn_atfp * un_bf(:,:) ) 1202 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - rn_atfp * vn_bf(:,:) ) 1206 1203 ! 1207 1204 ! Update corrective fluxes for next time step: 1208 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:))1209 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:))1205 un_bf(:,:) = rn_atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1206 vn_bf(:,:) = rn_atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 1210 1207 ELSE 1211 1208 un_bf(:,:) = 0._wp … … 1222 1219 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1223 1220 DO jk=1,jpkm1 1224 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_ 2dt_b1225 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_ 2dt_b1221 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_Dt 1222 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_Dt 1226 1223 END DO 1227 1224 ELSE … … 1229 1226 DO jj = 1, jpjm1 1230 1227 DO ji = 1, jpim1 ! NO Vector Opt. 1231 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 1232 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1233 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1234 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 1235 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1236 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1228 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1229 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1230 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1231 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1237 1232 END DO 1238 1233 END DO … … 1240 1235 ! 1241 1236 DO jk=1,jpkm1 1242 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_ 2dt_b1243 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_ 2dt_b1237 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_Dt 1238 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_Dt 1244 1239 END DO 1245 1240 ! Save barotropic velocities not transport: … … 1306 1301 LOGICAL , INTENT(in ) :: ll_fw ! forward time splitting =.true. 1307 1302 INTEGER , INTENT(inout) :: jpit ! cycle length 1308 REAL(wp), DIMENSION(3*nn_ baro), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights1303 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights 1309 1304 ! 1310 1305 INTEGER :: jic, jn, ji ! temporary integers … … 1316 1311 ! 1317 1312 ! Set time index when averaged value is requested 1318 IF ( ll_fw ) THEN ; jic = nn_ baro1319 ELSE ; jic = 2 * nn_ baro1313 IF ( ll_fw ) THEN ; jic = nn_e 1314 ELSE ; jic = 2 * nn_e 1320 1315 ENDIF 1321 1316 … … 1323 1318 ! 1324 1319 IF (ll_av) THEN !* Define simple boxcar window for primary weights 1325 ! ! (width = nn_ baro, centered around jic)1320 ! ! (width = nn_e, centered around jic) 1326 1321 SELECT CASE ( nn_bt_flt ) 1327 1322 CASE( 0 ) ! No averaging … … 1329 1324 jpit = jic 1330 1325 ! 1331 CASE( 1 ) ! Boxcar, width = nn_ baro1332 DO jn = 1, 3*nn_ baro1333 za1 = ABS(float(jn-jic))/float(nn_ baro)1326 CASE( 1 ) ! Boxcar, width = nn_e 1327 DO jn = 1, 3*nn_e 1328 za1 = ABS(float(jn-jic))/float(nn_e) 1334 1329 IF ( za1 < 0.5_wp ) THEN 1335 1330 zwgt1(jn) = 1._wp … … 1338 1333 END DO 1339 1334 ! 1340 CASE( 2 ) ! Boxcar, width = 2 * nn_ baro1341 DO jn = 1, 3*nn_ baro1342 za1 = ABS(float(jn-jic))/float(nn_ baro)1335 CASE( 2 ) ! Boxcar, width = 2 * nn_e 1336 DO jn = 1, 3*nn_e 1337 za1 = ABS(float(jn-jic))/float(nn_e) 1343 1338 IF ( za1 < 1._wp ) THEN 1344 1339 zwgt1(jn) = 1._wp … … 1474 1469 1475 1470 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1476 IF( ln_bt_auto ) nn_ baro = CEILING( rdt / rn_bt_cmax * zcmax)1471 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 1477 1472 1478 r dtbt = rdt / REAL( nn_baro, wp )1479 zcmax = zcmax * r dtbt1473 rDt_e = rn_Dt / REAL( nn_e , wp ) 1474 zcmax = zcmax * rDt_e 1480 1475 ! Print results 1481 1476 IF(lwp) WRITE(numout,*) … … 1483 1478 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 1484 1479 IF( ln_bt_auto ) THEN 1485 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_ baro'1480 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e ' 1486 1481 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1487 1482 ELSE 1488 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_ baro in namelist nn_baro = ', nn_baro1483 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 1489 1484 ENDIF 1490 1485 1491 1486 IF(ln_bt_av) THEN 1492 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_ barotime steps is on '1487 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' 1493 1488 ELSE 1494 1489 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' … … 1510 1505 SELECT CASE ( nn_bt_flt ) 1511 1506 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1512 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_ baro'1513 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_ baro'1507 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1508 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1514 1509 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1515 1510 END SELECT 1516 1511 ! 1517 1512 IF(lwp) WRITE(numout,*) ' ' 1518 IF(lwp) WRITE(numout,*) ' nn_ baro = ', nn_baro1519 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt1520 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax1513 IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e 1514 IF(lwp) WRITE(numout,*) ' external mode time step is : rDt_e', rDt_e, ' [s]' 1515 IF(lwp) WRITE(numout,*) ' Maximum Courant number is : ', zcmax 1521 1516 ! 1522 1517 IF(lwp) WRITE(numout,*) ' Time diffusion parameter rn_bt_alpha: ', rn_bt_alpha … … 1529 1524 ENDIF 1530 1525 IF( zcmax>0.9_wp ) THEN 1531 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_ baro!' )1526 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1532 1527 ENDIF 1533 1528 !
Note: See TracChangeset
for help on using the changeset viewer.
