Changeset 11480
- Timestamp:
- 2019-08-29T11:23:25+02:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps
- Files:
-
- 3 deleted
- 38 edited
- 3 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn.F90
r10564 r11480 60 60 CONTAINS 61 61 62 SUBROUTINE ice_dyn( kt )62 SUBROUTINE ice_dyn( kt, Kmm ) 63 63 !!------------------------------------------------------------------- 64 64 !! *** ROUTINE ice_dyn *** … … 73 73 !!-------------------------------------------------------------------- 74 74 INTEGER, INTENT(in) :: kt ! ice time step 75 INTEGER, INTENT(in) :: Kmm ! ocean time level index 75 76 !! 76 77 INTEGER :: ji, jj, jl ! dummy loop indices … … 92 93 tau_icebfr(:,:) = 0._wp 93 94 DO jl = 1, jpl 94 WHERE( h_i_b(:,:,jl) > ht _n(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr95 WHERE( h_i_b(:,:,jl) > ht(:,:) * rn_depfra ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 95 96 END DO 96 97 ENDIF … … 121 122 122 123 CASE ( np_dynALL ) !== all dynamical processes ==! 123 CALL ice_dyn_rhg ( kt )! -- rheology124 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 124 125 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 125 126 CALL ice_dyn_rdgrft( kt ) ! -- ridging/rafting … … 127 128 128 129 CASE ( np_dynRHGADV ) !== no ridge/raft & no corrections ==! 129 CALL ice_dyn_rhg ( kt )! -- rheology130 CALL ice_dyn_rhg ( kt, Kmm ) ! -- rheology 130 131 CALL ice_dyn_adv ( kt ) ; CALL Hbig( zhi_max, zhs_max, zhip_max ) ! -- advection of ice + correction on ice thickness 131 132 CALL Hpiling ! -- simple pile-up (replaces ridging/rafting) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rhg.F90
r10413 r11480 47 47 CONTAINS 48 48 49 SUBROUTINE ice_dyn_rhg( kt )49 SUBROUTINE ice_dyn_rhg( kt, Kmm ) 50 50 !!------------------------------------------------------------------- 51 51 !! *** ROUTINE ice_dyn_rhg *** … … 58 58 !!-------------------------------------------------------------------- 59 59 INTEGER, INTENT(in) :: kt ! ice time step 60 INTEGER, INTENT(in) :: Kmm ! ocean time level index 60 61 !!-------------------------------------------------------------------- 61 62 ! controls … … 76 77 CASE( np_rhgEVP ) ! Elasto-Viscous-Plastic ! 77 78 ! !------------------------! 78 CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )79 CALL ice_dyn_rhg_evp( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 79 80 ! 80 81 END SELECT -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icedyn_rhg_evp.F90
r10555 r11480 56 56 CONTAINS 57 57 58 SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )58 SUBROUTINE ice_dyn_rhg_evp( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) 59 59 !!------------------------------------------------------------------- 60 60 !! *** SUBROUTINE ice_dyn_rhg_evp *** … … 109 109 !!------------------------------------------------------------------- 110 110 INTEGER , INTENT(in ) :: kt ! time step 111 INTEGER , INTENT(in ) :: Kmm ! ocean time level index 111 112 REAL(wp), DIMENSION(:,:), INTENT(inout) :: pstress1_i, pstress2_i, pstress12_i ! 112 113 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! … … 335 336 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 336 337 ! ice-bottom stress at U points 337 zvCr = zaU(ji,jj) * rn_depfra * hu _n(ji,jj)338 zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm) 338 339 zTauU_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 339 340 ! ice-bottom stress at V points 340 zvCr = zaV(ji,jj) * rn_depfra * hv _n(ji,jj)341 zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm) 341 342 zTauV_ib(ji,jj) = rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 342 343 ! ice_bottom stress at T points 343 zvCr = at_i(ji,jj) * rn_depfra * ht _n(ji,jj)344 zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj) 344 345 tau_icebfr(ji,jj) = rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 345 346 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/ICE/icestp.F90
r10998 r11480 95 95 CONTAINS 96 96 97 SUBROUTINE ice_stp( kt, Kbb, ksbc )97 SUBROUTINE ice_stp( kt, Kbb, Kmm, ksbc ) 98 98 !!--------------------------------------------------------------------- 99 99 !! *** ROUTINE ice_stp *** … … 115 115 !! utau, vtau, taum, wndm, qns , qsr, emp , sfx 116 116 !!--------------------------------------------------------------------- 117 INTEGER, INTENT(in) :: kt ! ocean time step118 INTEGER, INTENT(in) :: Kbb ! ocean time level index119 INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled)117 INTEGER, INTENT(in) :: kt ! ocean time step 118 INTEGER, INTENT(in) :: Kbb, Kmm ! ocean time level indices 119 INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk, or Pure Coupled) 120 120 ! 121 121 INTEGER :: jl ! dummy loop index … … 161 161 ! 162 162 IF( ln_icedyn .AND. .NOT.lk_c1d ) & 163 & CALL ice_dyn( kt )! -- Ice dynamics163 & CALL ice_dyn( kt, Kmm ) ! -- Ice dynamics 164 164 ! 165 165 ! !== lateral boundary conditions ==! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90
r11428 r11480 116 116 END DO 117 117 DO jj = 1, jpj 118 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu _a(ibdy1:ibdy2,jj)118 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 119 119 END DO 120 120 ENDIF … … 137 137 END DO 138 138 DO jj=1,jpj 139 zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu _a(ibdy1:ibdy2,jj)139 zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 140 140 END DO 141 141 … … 158 158 END DO 159 159 DO jj = 1, jpj 160 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv _a(ibdy1:ibdy2,jj)160 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a) 161 161 END DO 162 162 DO jk = 1, jpkm1 … … 194 194 END DO 195 195 DO jj = 1, jpj 196 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu _a(ibdy1:ibdy2,jj)196 uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 197 197 END DO 198 198 ENDIF … … 217 217 END DO 218 218 DO jj=1,jpj 219 zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu _a(ibdy1:ibdy2,jj)219 zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu(ibdy1:ibdy2,jj,Krhs_a) 220 220 END DO 221 221 … … 242 242 END DO 243 243 DO jj = 1, jpj 244 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv _a(ibdy1:ibdy2,jj)244 zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv(ibdy1:ibdy2,jj,Krhs_a) 245 245 END DO 246 246 DO jk = 1, jpkm1 … … 278 278 END DO 279 279 DO ji=1,jpi 280 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv _a(ji,jbdy1:jbdy2)280 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 281 281 END DO 282 282 ENDIF … … 302 302 END DO 303 303 DO ji = 1, jpi 304 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv _a(ji,jbdy1:jbdy2)304 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 305 305 END DO 306 306 … … 325 325 END DO 326 326 DO ji = 1, jpi 327 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu _a(ji,jbdy1:jbdy2)327 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a) 328 328 END DO 329 329 … … 362 362 END DO 363 363 DO ji=1,jpi 364 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv _a(ji,jbdy1:jbdy2)364 vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 365 365 END DO 366 366 ENDIF … … 386 386 END DO 387 387 DO ji = 1, jpi 388 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv _a(ji,jbdy1:jbdy2)388 zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv(ji,jbdy1:jbdy2,Krhs_a) 389 389 END DO 390 390 … … 411 411 END DO 412 412 DO ji = 1, jpi 413 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu _a(ji,jbdy1:jbdy2)413 zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu(ji,jbdy1:jbdy2,Krhs_a) 414 414 END DO 415 415 … … 1129 1129 ! 1130 1130 IF( before ) THEN 1131 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu _n(i1:i2,j1:j2) * uu_b(i1:i2,j1:j2,Kmm_a)1131 ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu(i1:i2,j1:j2,Kmm_a) * uu_b(i1:i2,j1:j2,Kmm_a) 1132 1132 ELSE 1133 1133 western_side = (nb == 1).AND.(ndir == 1) … … 1182 1182 ! 1183 1183 IF( before ) THEN 1184 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv _n(i1:i2,j1:j2) * vv_b(i1:i2,j1:j2,Kmm_a)1184 ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv(i1:i2,j1:j2,Kmm_a) * vv_b(i1:i2,j1:j2,Kmm_a) 1185 1185 ELSE 1186 1186 western_side = (nb == 1).AND.(ndir == 1) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_update.F90
r11428 r11480 234 234 ! uu(:,:,:,Krhs_a) = e3u(:,:,:,Kbb_a) 235 235 ! vv(:,:,:,Krhs_a) = e3v(:,:,:,Kbb_a) 236 hu _a(:,:) = hu_n(:,:)237 hv _a(:,:) = hv_n(:,:)236 hu(:,:,Krhs_a) = hu(:,:,Kmm_a) 237 hv(:,:,Krhs_a) = hv(:,:,Kmm_a) 238 238 239 239 ! 1) NOW fields … … 251 251 ! Update total depths: 252 252 ! -------------------- 253 hu _n(:,:) = 0._wp ! Ocean depth at U-points254 hv _n(:,:) = 0._wp ! Ocean depth at V-points253 hu(:,:,Kmm_a) = 0._wp ! Ocean depth at U-points 254 hv(:,:,Kmm_a) = 0._wp ! Ocean depth at V-points 255 255 DO jk = 1, jpkm1 256 hu _n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk)257 hv _n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk)256 hu(:,:,Kmm_a) = hu(:,:,Kmm_a) + e3u(:,:,jk,Kmm_a) * umask(:,:,jk) 257 hv(:,:,Kmm_a) = hv(:,:,Kmm_a) + e3v(:,:,jk,Kmm_a) * vmask(:,:,jk) 258 258 END DO 259 259 ! ! Inverse of the local depth 260 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )261 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )260 r1_hu(:,:,Kmm_a) = ssumask(:,:) / ( hu(:,:,Kmm_a) + 1._wp - ssumask(:,:) ) 261 r1_hv(:,:,Kmm_a) = ssvmask(:,:) / ( hv(:,:,Kmm_a) + 1._wp - ssvmask(:,:) ) 262 262 263 263 … … 276 276 ! Update total depths: 277 277 ! -------------------- 278 hu _b(:,:) = 0._wp ! Ocean depth at U-points279 hv _b(:,:) = 0._wp ! Ocean depth at V-points278 hu(:,:,Kbb_a) = 0._wp ! Ocean depth at U-points 279 hv(:,:,Kbb_a) = 0._wp ! Ocean depth at V-points 280 280 DO jk = 1, jpkm1 281 hu _b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk)282 hv _b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk)281 hu(:,:,Kbb_a) = hu(:,:,Kbb_a) + e3u(:,:,jk,Kbb_a) * umask(:,:,jk) 282 hv(:,:,Kbb_a) = hv(:,:,Kbb_a) + e3v(:,:,jk,Kbb_a) * vmask(:,:,jk) 283 283 END DO 284 284 ! ! Inverse of the local depth 285 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) )286 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )285 r1_hu(:,:,Kbb_a) = ssumask(:,:) / ( hu(:,:,Kbb_a) + 1._wp - ssumask(:,:) ) 286 r1_hv(:,:,Kbb_a) = ssvmask(:,:) / ( hv(:,:,Kbb_a) + 1._wp - ssvmask(:,:) ) 287 287 ENDIF 288 288 ! … … 636 636 IF (western_side) THEN 637 637 DO jj=j1,j2 638 zcor = uu_b(i1-1,jj,Kmm_a) * hu _a(i1-1,jj) * r1_hu_n(i1-1,jj) - uu_b(i1-1,jj,Kmm_a)638 zcor = uu_b(i1-1,jj,Kmm_a) * hu(i1-1,jj,Krhs_a) * r1_hu(i1-1,jj,Kmm_a) - uu_b(i1-1,jj,Kmm_a) 639 639 uu_b(i1-1,jj,Kmm_a) = uu_b(i1-1,jj,Kmm_a) + zcor 640 640 DO jk=1,jpkm1 … … 646 646 IF (eastern_side) THEN 647 647 DO jj=j1,j2 648 zcor = uu_b(i2+1,jj,Kmm_a) * hu _a(i2+1,jj) * r1_hu_n(i2+1,jj) - uu_b(i2+1,jj,Kmm_a)648 zcor = uu_b(i2+1,jj,Kmm_a) * hu(i2+1,jj,Krhs_a) * r1_hu(i2+1,jj,Kmm_a) - uu_b(i2+1,jj,Kmm_a) 649 649 uu_b(i2+1,jj,Kmm_a) = uu_b(i2+1,jj,Kmm_a) + zcor 650 650 DO jk=1,jpkm1 … … 829 829 IF (southern_side) THEN 830 830 DO ji=i1,i2 831 zcor = vv_b(ji,j1-1,Kmm_a) * hv _a(ji,j1-1) * r1_hv_n(ji,j1-1) - vv_b(ji,j1-1,Kmm_a)831 zcor = vv_b(ji,j1-1,Kmm_a) * hv(ji,j1-1,Krhs_a) * r1_hv(ji,j1-1,Kmm_a) - vv_b(ji,j1-1,Kmm_a) 832 832 vv_b(ji,j1-1,Kmm_a) = vv_b(ji,j1-1,Kmm_a) + zcor 833 833 DO jk=1,jpkm1 … … 839 839 IF (northern_side) THEN 840 840 DO ji=i1,i2 841 zcor = vv_b(ji,j2+1,Kmm_a) * hv _a(ji,j2+1) * r1_hv_n(ji,j2+1) - vv_b(ji,j2+1,Kmm_a)841 zcor = vv_b(ji,j2+1,Kmm_a) * hv(ji,j2+1,Krhs_a) * r1_hv(ji,j2+1,Kmm_a) - vv_b(ji,j2+1,Kmm_a) 842 842 vv_b(ji,j2+1,Kmm_a) = vv_b(ji,j2+1,Kmm_a) + zcor 843 843 DO jk=1,jpkm1 … … 869 869 DO jj=j1,j2 870 870 DO ji=i1,i2 871 tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu _n(ji,jj) * e2u(ji,jj)871 tabres(ji,jj) = zrhoy * uu_b(ji,jj,Kmm_a) * hu(ji,jj,Kmm_a) * e2u(ji,jj) 872 872 END DO 873 873 END DO … … 883 883 END DO 884 884 ! 885 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu _n(ji,jj)885 zcorr = (tabres(ji,jj) - spgu(ji,jj)) * r1_hu(ji,jj,Kmm_a) 886 886 DO jk=1,jpkm1 887 887 uu(ji,jj,jk,Kmm_a) = uu(ji,jj,jk,Kmm_a) + zcorr * umask(ji,jj,jk) … … 891 891 IF ( .NOT.ln_dynspg_ts .OR. (ln_dynspg_ts.AND.(.NOT.ln_bt_fw)) ) THEN 892 892 IF ( .NOT.( lk_agrif_fstep .AND. (neuler==0) ) ) THEN ! Add asselin part 893 zcorr = ( tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu_a(ji,jj) ) * r1_hu_b(ji,jj)893 zcorr = (tabres(ji,jj) - uu_b(ji,jj,Kmm_a) * hu(ji,jj,Krhs_a)) * r1_hu(ji,jj,Kbb_a) 894 894 uu_b(ji,jj,Kbb_a) = uu_b(ji,jj,Kbb_a) + atfp * zcorr * umask(ji,jj,1) 895 895 END IF 896 896 ENDIF 897 uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu _n(ji,jj) * umask(ji,jj,1)897 uu_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hu(ji,jj,Kmm_a) * umask(ji,jj,1) 898 898 ! 899 899 ! Correct "before" velocities to hold correct bt component: … … 903 903 END DO 904 904 ! 905 zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu _b(ji,jj)905 zcorr = uu_b(ji,jj,Kbb_a) - spgu(ji,jj) * r1_hu(ji,jj,Kbb_a) 906 906 DO jk=1,jpkm1 907 907 uu(ji,jj,jk,Kbb_a) = uu(ji,jj,jk,Kbb_a) + zcorr * umask(ji,jj,jk) … … 935 935 DO jj=j1,j2 936 936 DO ji=i1,i2 937 tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv _n(ji,jj) * e1v(ji,jj)937 tabres(ji,jj) = zrhox * vv_b(ji,jj,Kmm_a) * hv(ji,jj,Kmm_a) * e1v(ji,jj) 938 938 END DO 939 939 END DO … … 949 949 END DO 950 950 ! 951 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv _n(ji,jj)951 zcorr = (tabres(ji,jj) - spgv(ji,jj)) * r1_hv(ji,jj,Kmm_a) 952 952 DO jk=1,jpkm1 953 953 vv(ji,jj,jk,Kmm_a) = vv(ji,jj,jk,Kmm_a) + zcorr * vmask(ji,jj,jk) … … 957 957 IF ( .NOT.ln_dynspg_ts .OR. ( ln_dynspg_ts .AND. ( .NOT.ln_bt_fw ) ) ) THEN 958 958 IF ( .NOT. ( lk_agrif_fstep .AND. ( neuler==0 ) ) ) THEN ! Add asselin part 959 zcorr = ( tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv_a(ji,jj) ) * r1_hv_b(ji,jj)959 zcorr = (tabres(ji,jj) - vv_b(ji,jj,Kmm_a) * hv(ji,jj,Krhs_a)) * r1_hv(ji,jj,Kbb_a) 960 960 vv_b(ji,jj,Kbb_a) = vv_b(ji,jj,Kbb_a) + atfp * zcorr * vmask(ji,jj,1) 961 961 END IF 962 962 ENDIF 963 vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv _n(ji,jj) * vmask(ji,jj,1)963 vv_b(ji,jj,Kmm_a) = tabres(ji,jj) * r1_hv(ji,jj,Kmm_a) * vmask(ji,jj,1) 964 964 ! 965 965 ! Correct "before" velocities to hold correct bt component: … … 969 969 END DO 970 970 ! 971 zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv _b(ji,jj)971 zcorr = vv_b(ji,jj,Kbb_a) - spgv(ji,jj) * r1_hv(ji,jj,Kbb_a) 972 972 DO jk=1,jpkm1 973 973 vv(ji,jj,jk,Kbb_a) = vv(ji,jj,jk,Kbb_a) + zcorr * vmask(ji,jj,jk) … … 1387 1387 ! 1388 1388 ! Update total depth: 1389 ht _n(i1:i2,j1:j2) = 0._wp1389 ht(i1:i2,j1:j2) = 0._wp 1390 1390 DO jk = 1, jpkm1 1391 ht _n(i1:i2,j1:j2) = ht_n(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk)1391 ht(i1:i2,j1:j2) = ht(i1:i2,j1:j2) + e3t(i1:i2,j1:j2,jk,Kmm_a) * tmask(i1:i2,j1:j2,jk) 1392 1392 END DO 1393 1393 ! … … 1396 1396 gdept(i1:i2,j1:j2,1,Kmm_a) = 0.5_wp * e3w(i1:i2,j1:j2,1,Kmm_a) 1397 1397 gdepw(i1:i2,j1:j2,1,Kmm_a) = 0.0_wp 1398 gde3w(i1:i2,j1:j2,1 ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( ht _n(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2) ) ! Last term in the rhs is ssh1398 gde3w(i1:i2,j1:j2,1 ) = gdept(i1:i2,j1:j2,1,Kmm_a) - ( ht(i1:i2,j1:j2) - ht_0(i1:i2,j1:j2) ) ! Last term in the rhs is ssh 1399 1399 ! 1400 1400 DO jk = 2, jpk … … 1409 1409 gdept(ji,jj,jk,Kmm_a) = zcoef * ( gdepw(ji,jj,jk ,Kmm_a) + 0.5 * e3w(ji,jj,jk ,Kmm_a) ) & 1410 1410 & + ( 1._wp - zcoef ) * ( gdept(ji,jj,jk-1,Kmm_a) + e3w(ji,jj,jk ,Kmm_a) ) 1411 gde3w(ji,jj,jk ) = gdept(ji,jj,jk ,Kmm_a) - ( ht _n(ji,jj)-ht_0(ji,jj) ) ! Last term in the rhs is ssh1411 gde3w(ji,jj,jk ) = gdept(ji,jj,jk ,Kmm_a) - ( ht(ji,jj)-ht_0(ji,jj) ) ! Last term in the rhs is ssh 1412 1412 END DO 1413 1413 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ASM/asminc.F90
r10969 r11480 805 805 ELSE 806 806 ALLOCATE( ztim(jpi,jpj) ) 807 ztim(:,:) = ssh_iau(:,:) / ( ht _n(:,:) + 1.0 - ssmask(:,:) )807 ztim(:,:) = ssh_iau(:,:) / ( ht(:,:) + 1.0 - ssmask(:,:) ) 808 808 DO jk = 1, jpkm1 809 809 phdivn(:,:,jk) = phdivn(:,:,jk) - ztim(:,:) * tmask(:,:,jk) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydta.F90
r10957 r11480 255 255 CALL fld_read( kt=kt, kn_fsbc=1, sd=bf(jstart:jend), & 256 256 & map=nbmap_ptr(jstart:jend), kt_offset=time_offset, jpk_bdy=nb_jpk_bdy, & 257 & fvl=ln_full_vel_array(jbdy) )257 & fvl=ln_full_vel_array(jbdy), Kmm=Kmm ) 258 258 ENDIF 259 259 ! If full velocities in boundary data then split into barotropic and baroclinic data … … 270 270 & + e3u(ii,ij,ik,Kmm) * umask(ii,ij,ik) * dta%u3d(ib,ik) 271 271 END DO 272 dta%u2d(ib) = dta%u2d(ib) * r1_hu _n(ii,ij)272 dta%u2d(ib) = dta%u2d(ib) * r1_hu(ii,ij,Kmm) 273 273 DO ik = 1, jpkm1 274 274 dta%u3d(ib,ik) = dta%u3d(ib,ik) - dta%u2d(ib) … … 284 284 & + e3v(ii,ij,ik,Kmm) * vmask(ii,ij,ik) * dta%v3d(ib,ik) 285 285 END DO 286 dta%v2d(ib) = dta%v2d(ib) * r1_hv _n(ii,ij)286 dta%v2d(ib) = dta%v2d(ib) * r1_hv(ii,ij,Kmm) 287 287 DO ik = 1, jpkm1 288 288 dta%v3d(ib,ik) = dta%v3d(ib,ik) - dta%v2d(ib) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/BDY/bdydyn.F90
r10957 r11480 78 78 zva2d(:,:) = zva2d(:,:) + e3v(:,:,jk,Kaa) * pvv(:,:,jk,Kaa) * vmask(:,:,jk) 79 79 END DO 80 zua2d(:,:) = zua2d(:,:) * r1_hu _a(:,:)81 zva2d(:,:) = zva2d(:,:) * r1_hv _a(:,:)80 zua2d(:,:) = zua2d(:,:) * r1_hu(:,:,Kaa) 81 zva2d(:,:) = zva2d(:,:) * r1_hv(:,:,Kaa) 82 82 83 83 DO jk = 1 , jpkm1 … … 99 99 !------------------------------------------------------- 100 100 101 IF( ll_dyn2d ) CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu _a(:,:), r1_hv_a(:,:), ssh(:,:,Kaa) )101 IF( ll_dyn2d ) CALL bdy_dyn2d( kt, zua2d, zva2d, uu_b(:,:,Kbb), vv_b(:,:,Kbb), r1_hu(:,:,Kaa), r1_hv(:,:,Kaa), ssh(:,:,Kaa) ) 102 102 103 103 IF( ll_dyn3d ) CALL bdy_dyn3d( kt, Kbb, puu, pvv, Kaa ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/C1D/step_c1d.F90
r11001 r11480 107 107 IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs ) ! horizontal & vertical advection 108 108 IF( ln_zdfosm ) CALL tra_osm( kstp, Nnn , ts, Nrhs ) ! OSMOSIS non-local tracer fluxes 109 CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing109 CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing 110 110 CALL eos( ts(:,:,:,:,Nnn), rhd, rhop, gdept_0(:,:,:) ) ! now potential density for zdfmxl 111 IF( ln_zdfnpc ) CALL tra_npc( kstp, Nnn, Nrhs, ts, Naa ) ! applied non penetrative convective adjustment on (t,s)112 CALL tra_ nxt( kstp, Nbb, Nnn, Nrhs, Naa ) ! tracer fields at next time step111 IF( ln_zdfnpc ) CALL tra_npc( kstp, Nnn, Nrhs, ts, Naa ) ! applied non penetrative convective adjustment on (t,s) 112 CALL tra_atf( kstp, Nbb, Nnn, Nrhs, Naa, ts ) ! time filtering of "now" tracer fields 113 113 114 114 … … 124 124 IF( ln_zdfosm ) CALL dyn_osm ( kstp, Nnn , uu, vv, Nrhs ) ! OSMOSIS non-local velocity fluxes 125 125 CALL dyn_zdf ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! vertical diffusion 126 CALL dyn_nxt ( kstp, Nbb, Nnn, Naa ) ! lateral velocity at next time step 127 IF(.NOT.ln_linssh)CALL ssh_swp ( kstp, Nbb, Nnn, Naa ) ! swap of sea surface height 128 129 IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa ) ! swap of vertical scale factors 126 CALL dyn_atf ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v ) ! time filtering of "now" fields 127 IF(.NOT.ln_linssh)CALL ssh_atf ( kstp, Nbb, Nnn, Naa , ssh ) ! time filtering of "now" sea surface height 128 ! 129 ! Swap time levels 130 Nrhs = Nbb 131 Nbb = Nnn 132 Nnn = Naa 133 Naa = Nrhs 134 ! 135 IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa ) ! swap of vertical scale factors 130 136 131 137 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/dom_oce.F90
r11036 r11480 12 12 !! 3.7 ! 2015-11 (G. Madec) introduce surface and scale factor ratio 13 13 !! - ! 2015-11 (G. Madec, A. Coward) time varying zgr by default 14 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename prognostic variables in preparation for new time scheme. 14 15 !!---------------------------------------------------------------------- 15 16 … … 121 122 REAL(wp), PUBLIC, ALLOCATABLE, SAVE , DIMENSION(:,:) :: e1e2f , r1_e1e2f !: associated metrics at f-point 122 123 ! 123 REAL(wp), PUBLIC, ALLOCATABLE, SAVE 124 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ff_f , ff_t !: Coriolis factor at f- & t-points [1/s] 124 125 !!---------------------------------------------------------------------- 125 126 !! vertical coordinate and scale factors … … 149 150 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w 150 151 151 ! ! ref. ! before ! now ! after ! 152 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 , ht_n !: t-depth [m] 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 , hu_b , hu_n , hu_a !: u-depth [m] 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 , hv_b , hv_n , hv_a !: v-depth [m] 155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hu_b , r1_hu_n , r1_hu_a !: inverse of u-depth [1/m] 156 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: r1_hv_b , r1_hv_n , r1_hv_a !: inverse of v-depth [1/m] 152 ! ! reference heights of water column 153 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0 !: t-depth [m] 154 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hu_0 !: u-depth [m] 155 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hv_0 !: v-depth [m] 156 ! time-dependent heights of water column 157 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht !: height of water column at T-points [m] 158 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hu, hv, r1_hu, r1_hv !: height of water column [m] and reciprocal [1/m] 157 159 158 160 INTEGER, PUBLIC :: nla10 !: deepest W level Above ~10m (nlb10 - 1) … … 268 270 & e3t (jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f (jpi,jpj,jpk) , e3w (jpi,jpj,jpk,jpt) , & 269 271 & e3uw_0(jpi,jpj,jpk) , e3vw_0(jpi,jpj,jpk) , & 270 & e3uw (jpi,jpj,jpk,jpt) , e3vw (jpi,jpj,jpk,jpt) , STAT=ierr(5) ) 271 ! 272 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 273 & hu_b(jpi,jpj) , hv_b(jpi,jpj) , r1_hu_b(jpi,jpj) , r1_hv_b(jpi,jpj) , & 274 & ht_n(jpi,jpj) , hu_n(jpi,jpj) , hv_n(jpi,jpj) , r1_hu_n(jpi,jpj) , r1_hv_n(jpi,jpj) , & 275 & hu_a(jpi,jpj) , hv_a(jpi,jpj) , r1_hu_a(jpi,jpj) , r1_hv_a(jpi,jpj) , STAT=ierr(6) ) 272 & e3uw (jpi,jpj,jpk,jpt) , e3vw (jpi,jpj,jpk,jpt) , STAT=ierr(5) ) 273 ! 274 ALLOCATE( ht_0(jpi,jpj) , hu_0(jpi,jpj) , hv_0(jpi,jpj) , & 275 & ht (jpi,jpj) , hu( jpi,jpj,jpt), hv( jpi,jpj,jpt) , r1_hu(jpi,jpj,jpt) , r1_hv(jpi,jpj,jpt) , & 276 & STAT=ierr(6) ) 276 277 ! 277 278 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domain.F90
r10978 r11480 163 163 ! 164 164 ! before ! now ! after ! 165 gdept(:,:,:,Kbb) = gdept_0 ; gdept(:,:,:,Kmm) = gdept_0 ! ---! depth of grid-points166 gdepw(:,:,:,Kbb) = gdepw_0 ; gdepw(:,:,:,Kmm) = gdepw_0 ! ---!165 gdept(:,:,:,Kbb) = gdept_0 ; gdept(:,:,:,Kmm) = gdept_0 ; gdept(:,:,:,Kaa) = gdept_0 ! depth of grid-points 166 gdepw(:,:,:,Kbb) = gdepw_0 ; gdepw(:,:,:,Kmm) = gdepw_0 ; gdepw(:,:,:,Kaa) = gdepw_0 ! 167 167 gde3w = gde3w_0 ! --- ! 168 168 ! … … 171 171 e3v(:,:,:,Kbb) = e3v_0 ; e3v(:,:,:,Kmm) = e3v_0 ; e3v(:,:,:,Kaa) = e3v_0 ! 172 172 e3f = e3f_0 ! --- ! 173 e3w(:,:,:,Kbb) = e3w_0 ; e3w(:,:,:,Kmm) = e3w_0 ! --- !174 e3uw(:,:,:,Kbb) = e3uw_0 ; e3uw(:,:,:,Kmm) = e3uw_0 ! --- !175 e3vw(:,:,:,Kbb) = e3vw_0 ; e3vw(:,:,:,Kmm) = e3vw_0 ! ---!173 e3w(:,:,:,Kbb) = e3w_0 ; e3w(:,:,:,Kmm) = e3w_0 ; e3w(:,:,:,Kaa) = e3w_0 ! 174 e3uw(:,:,:,Kbb) = e3uw_0 ; e3uw(:,:,:,Kmm) = e3uw_0 ; e3uw(:,:,:,Kaa) = e3uw_0 ! 175 e3vw(:,:,:,Kbb) = e3vw_0 ; e3vw(:,:,:,Kmm) = e3vw_0 ; e3vw(:,:,:,Kaa) = e3vw_0 ! 176 176 ! 177 177 z1_hu_0(:,:) = ssumask(:,:) / ( hu_0(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF … … 179 179 ! 180 180 ! before ! now ! after ! 181 ht _n= ht_0 ! ! water column thickness182 hu _b = hu_0 ; hu_n = hu_0 ; hu_a= hu_0 !183 hv _b = hv_0 ; hv_n = hv_0 ; hv_a= hv_0 !184 r1_hu _b = z1_hu_0 ; r1_hu_n = z1_hu_0 ; r1_hu_a= z1_hu_0 ! inverse of water column thickness185 r1_hv _b = z1_hv_0 ; r1_hv_n = z1_hv_0 ; r1_hv_a= z1_hv_0 !181 ht = ht_0 ! ! water column thickness 182 hu(:,:,Kbb) = hu_0 ; hu(:,:,Kmm) = hu_0 ; hu(:,:,Kaa) = hu_0 ! 183 hv(:,:,Kbb) = hv_0 ; hv(:,:,Kmm) = hv_0 ; hv(:,:,Kaa) = hv_0 ! 184 r1_hu(:,:,Kbb) = z1_hu_0 ; r1_hu(:,:,Kmm) = z1_hu_0 ; r1_hu(:,:,Kaa) = z1_hu_0 ! inverse of water column thickness 185 r1_hv(:,:,Kbb) = z1_hv_0 ; r1_hv(:,:,Kmm) = z1_hv_0 ; r1_hv(:,:,Kaa) = z1_hv_0 ! 186 186 ! 187 187 ! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/domvvl.F90
r10978 r11480 13 13 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 14 14 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 !! dom_vvl_sf_ swp: Swap vertical scale factors and update the vertical grid15 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 16 16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 17 17 !! dom_vvl_rst : read/write restart file … … 37 37 PUBLIC dom_vvl_init ! called by domain.F90 38 38 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_ swp! called by step.F9039 PUBLIC dom_vvl_sf_update ! called by step.F90 40 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 41 … … 181 181 ! 182 182 ! !== thickness of the water column !! (ocean portion only) 183 ht _n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....184 hu _b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)185 hu _n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1)186 hv _b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)187 hv _n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1)183 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 184 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 185 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 186 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 187 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 188 188 DO jk = 2, jpkm1 189 ht _n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)190 hu _b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)191 hu _n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)192 hv _b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)193 hv _n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)189 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 190 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 191 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 192 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 193 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 194 194 END DO 195 195 ! 196 196 ! !== inverse of water column thickness ==! (u- and v- points) 197 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF198 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )199 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )200 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )197 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 198 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 199 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 200 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 201 201 202 202 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 550 550 ! *********************************** ! 551 551 552 hu _a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1)553 hv _a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1)552 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 553 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 554 554 DO jk = 2, jpkm1 555 hu _a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk)556 hv _a(:,:) = hv_a(:,:) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)555 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 556 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 557 557 END DO 558 558 ! ! Inverse of the local depth 559 559 !!gm BUG ? don't understand the use of umask_i here ..... 560 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )561 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )560 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 561 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 562 562 ! 563 563 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 566 566 567 567 568 SUBROUTINE dom_vvl_sf_ swp( kt, Kbb, Kmm, Kaa )569 !!---------------------------------------------------------------------- 570 !! *** ROUTINE dom_vvl_sf_ swp***568 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 569 !!---------------------------------------------------------------------- 570 !! *** ROUTINE dom_vvl_sf_update *** 571 571 !! 572 !! ** Purpose : compute time filter and swap of scale factors572 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 573 573 !! compute all depths and related variables for next time step 574 574 !! write outputs and restart file 575 575 !! 576 !! ** Method : - swap of e3t with trick for volume/tracer conservation 576 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 577 577 !! - reconstruct scale factor at other grid points (interpolate) 578 578 !! - recompute depths and water height fields 579 579 !! 580 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step580 !! ** Action : - tilde_e3t_(b/n) ready for next time step 581 581 !! - Recompute: 582 582 !! e3(u/v)_b … … 599 599 IF( ln_linssh ) RETURN ! No calculation in linear free surface 600 600 ! 601 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')601 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 602 602 ! 603 603 IF( kt == nit000 ) THEN 604 604 IF(lwp) WRITE(numout,*) 605 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'606 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'605 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 606 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 607 607 ENDIF 608 608 ! … … 619 619 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 620 620 ENDIF 621 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)622 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)623 624 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)625 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)626 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)627 621 628 622 ! Compute all missing vertical scale factor and depths … … 630 624 ! Horizontal scale factor interpolations 631 625 ! -------------------------------------- 632 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are al lready computed in dynnxt633 ! - JC - hu _b, hv_b, hur_b, hvr_b also626 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are already computed in dynnxt 627 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 634 628 635 629 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) … … 663 657 ! Local depth and Inverse of the local depth of the water 664 658 ! ------------------------------------------------------- 665 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 666 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 667 ! 668 ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 659 ! 660 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 669 661 DO jk = 2, jpkm1 670 ht _n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)662 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 671 663 END DO 672 664 … … 675 667 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 676 668 ! 677 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')678 ! 679 END SUBROUTINE dom_vvl_sf_ swp669 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 670 ! 671 END SUBROUTINE dom_vvl_sf_update 680 672 681 673 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/iscplrst.F90
r10978 r11480 108 108 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 109 109 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 110 hu _b (:,:) = hu_n (:,:)111 hv _b (:,:) = hv_n (:,:)112 r1_hu _b(:,:) = r1_hu_n(:,:)113 r1_hv _b(:,:) = r1_hv_n(:,:)110 hu (:,:,Kbb) = hu (:,:,Kmm) 111 hv (:,:,Kbb) = hv (:,:,Kmm) 112 r1_hu(:,:,Kbb) = r1_hu(:,:,Kmm) 113 r1_hv(:,:,Kbb) = r1_hv(:,:,Kmm) 114 114 ! 115 115 END SUBROUTINE iscpl_stp … … 240 240 ! t-, u- and v- water column thickness 241 241 ! ------------------------------------ 242 ht _n(:,:) = 0._wp ; hu_n(:,:) = 0._wp ; hv_n(:,:) = 0._wp242 ht(:,:) = 0._wp ; hu(:,:,Kmm) = 0._wp ; hv(:,:,Kmm) = 0._wp 243 243 DO jk = 1, jpkm1 244 hu _n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)245 hv _n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)246 ht _n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)244 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 245 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 246 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 247 247 END DO 248 248 ! ! Inverse of the local depth 249 r1_hu _n(:,:) = 1._wp / ( hu_n(:,:) + 1._wp - ssumask(:,:) ) * ssumask(:,:)250 r1_hv _n(:,:) = 1._wp / ( hv_n(:,:) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:)249 r1_hu(:,:,Kmm) = 1._wp / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) * ssumask(:,:) 250 r1_hv(:,:,Kmm) = 1._wp / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) * ssvmask(:,:) 251 251 252 252 END IF -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DOM/istate.F90
r10978 r11480 175 175 END DO 176 176 ! 177 uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu _n(:,:)178 vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv _n(:,:)177 uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm) 178 vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm) 179 179 ! 180 uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu _b(:,:)181 vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv _b(:,:)180 uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb) 181 vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb) 182 182 ! 183 183 END SUBROUTINE istate_init -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/dynspg_ts.F90
r11027 r11480 250 250 DO jj = 1, jpjm1 251 251 DO ji = 1, jpim1 252 zwz(ji,jj) = ( ht _n(ji ,jj+1) + ht_n(ji+1,jj+1) + &253 & ht _n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp252 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 253 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 254 254 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 255 255 END DO … … 258 258 DO jj = 1, jpjm1 259 259 DO ji = 1, jpim1 260 zwz(ji,jj) = ( ht _n (ji ,jj+1) + ht_n(ji+1,jj+1) &261 & + ht _n (ji ,jj ) + ht_n(ji+1,jj ) ) &260 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 261 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & 262 262 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 263 263 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) … … 282 282 DO jj = 2, jpj 283 283 DO ji = 2, jpi 284 z1_ht = ssmask(ji,jj) / ( ht _n(ji,jj) + 1._wp - ssmask(ji,jj) )284 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 285 285 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 286 286 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht … … 367 367 END DO 368 368 ! 369 zu_frc(:,:) = zu_frc(:,:) * r1_hu _n(:,:)370 zv_frc(:,:) = zv_frc(:,:) * r1_hv _n(:,:)369 zu_frc(:,:) = zu_frc(:,:) * r1_hu(:,:,Kmm) 370 zv_frc(:,:) = zv_frc(:,:) * r1_hv(:,:,Kmm) 371 371 ! 372 372 ! … … 388 388 ! ! -------------------------------------------------------- 389 389 ! 390 zwx(:,:) = puu_b(:,:,Kmm) * hu _n(:,:) * e2u(:,:) ! now fluxes391 zwy(:,:) = pvv_b(:,:,Kmm) * hv _n(:,:) * e1v(:,:)390 zwx(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 391 zwy(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) 392 392 ! 393 393 SELECT CASE( nvor_scheme ) … … 395 395 DO jj = 2, jpjm1 396 396 DO ji = 2, jpim1 ! vector opt. 397 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu _n(ji,jj) &398 & * ( e1e2t(ji+1,jj)*ht _n(ji+1,jj)*ff_t(ji+1,jj) * ( pvv_b(ji+1,jj,Kmm) + pvv_b(ji+1,jj-1,Kmm) ) &399 & + e1e2t(ji ,jj)*ht _n(ji ,jj)*ff_t(ji ,jj) * ( pvv_b(ji ,jj,Kmm) + pvv_b(ji ,jj-1,Kmm) ) )397 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu(ji,jj,Kmm) & 398 & * ( e1e2t(ji+1,jj)*ht(ji+1,jj)*ff_t(ji+1,jj) * ( pvv_b(ji+1,jj,Kmm) + pvv_b(ji+1,jj-1,Kmm) ) & 399 & + e1e2t(ji ,jj)*ht(ji ,jj)*ff_t(ji ,jj) * ( pvv_b(ji ,jj,Kmm) + pvv_b(ji ,jj-1,Kmm) ) ) 400 400 ! 401 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv _n(ji,jj) &402 & * ( e1e2t(ji,jj+1)*ht _n(ji,jj+1)*ff_t(ji,jj+1) * ( puu_b(ji,jj+1,Kmm) + puu_b(ji-1,jj+1,Kmm) ) &403 & + e1e2t(ji,jj )*ht _n(ji,jj )*ff_t(ji,jj ) * ( puu_b(ji,jj ,Kmm) + puu_b(ji-1,jj ,Kmm) ) )401 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv(ji,jj,Kmm) & 402 & * ( e1e2t(ji,jj+1)*ht(ji,jj+1)*ff_t(ji,jj+1) * ( puu_b(ji,jj+1,Kmm) + puu_b(ji-1,jj+1,Kmm) ) & 403 & + e1e2t(ji,jj )*ht(ji,jj )*ff_t(ji,jj ) * ( puu_b(ji,jj ,Kmm) + puu_b(ji-1,jj ,Kmm) ) ) 404 404 END DO 405 405 END DO … … 546 546 DO ji = fs_2, fs_jpim1 ! vector opt. 547 547 zu_frc(ji,jj) = zu_frc(ji,jj) + & 548 & MAX(r1_hu _n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj)548 & MAX(r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj) 549 549 zv_frc(ji,jj) = zv_frc(ji,jj) + & 550 & MAX(r1_hv _n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj)550 & MAX(r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj) 551 551 END DO 552 552 END DO … … 554 554 DO jj = 2, jpjm1 555 555 DO ji = fs_2, fs_jpim1 ! vector opt. 556 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu _n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj)557 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv _n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj)556 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 557 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 558 558 END DO 559 559 END DO … … 584 584 DO jj = 2, jpjm1 585 585 DO ji = fs_2, fs_jpim1 ! vector opt. 586 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu _n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)587 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv _n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)586 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj) 587 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj) 588 588 END DO 589 589 END DO … … 593 593 DO jj = 2, jpjm1 594 594 DO ji = fs_2, fs_jpim1 ! vector opt. 595 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu _n(ji,jj)596 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv _n(ji,jj)595 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 596 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) 597 597 END DO 598 598 END DO … … 601 601 DO jj = 2, jpjm1 602 602 DO ji = fs_2, fs_jpim1 ! vector opt. 603 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu _n(ji,jj)604 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv _n(ji,jj)603 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 604 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) 605 605 END DO 606 606 END DO … … 681 681 vn_e (:,:) = pvv_b(:,:,Kmm) 682 682 ! 683 hu_e (:,:) = hu _n(:,:)684 hv_e (:,:) = hv _n(:,:)685 hur_e (:,:) = r1_hu _n(:,:)686 hvr_e (:,:) = r1_hv _n(:,:)683 hu_e (:,:) = hu(:,:,Kmm) 684 hv_e (:,:) = hv(:,:,Kmm) 685 hur_e (:,:) = r1_hu(:,:,Kmm) 686 hvr_e (:,:) = r1_hv(:,:,Kmm) 687 687 ELSE ! CENTRED integration: start from BEFORE fields 688 688 sshn_e(:,:) = pssh(:,:,Kbb) … … 690 690 vn_e (:,:) = pvv_b(:,:,Kbb) 691 691 ! 692 hu_e (:,:) = hu _b(:,:)693 hv_e (:,:) = hv _b(:,:)694 hur_e (:,:) = r1_hu _b(:,:)695 hvr_e (:,:) = r1_hv _b(:,:)692 hu_e (:,:) = hu(:,:,Kbb) 693 hv_e (:,:) = hv(:,:,Kbb) 694 hur_e (:,:) = r1_hu(:,:,Kbb) 695 hvr_e (:,:) = r1_hv(:,:,Kbb) 696 696 ENDIF 697 697 ! … … 790 790 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 791 791 ELSE 792 zhup2_e(:,:) = hu _n(:,:)793 zhvp2_e(:,:) = hv _n(:,:)794 zhtp2_e(:,:) = ht _n(:,:)792 zhup2_e(:,:) = hu(:,:,Kmm) 793 zhvp2_e(:,:) = hv(:,:,Kmm) 794 zhtp2_e(:,:) = ht(:,:) 795 795 ENDIF 796 796 ! !* after ssh … … 1138 1138 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 1139 1139 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1140 & + hu _n(ji,jj) * zu_frc(ji,jj) ) &1140 & + hu(ji,jj,Kmm) * zu_frc(ji,jj) ) & 1141 1141 & ) * zhura 1142 1142 … … 1144 1144 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1145 1145 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1146 & + hv _n(ji,jj) * zv_frc(ji,jj) ) &1146 & + hv(ji,jj,Kmm) * zv_frc(ji,jj) ) & 1147 1147 & ) * zhvra 1148 1148 END DO … … 1257 1257 ! 1258 1258 DO jk=1,jpkm1 1259 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu _n(:,:) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu_b(:,:) ) * r1_2dt_b1260 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv _n(:,:) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv_b(:,:) ) * r1_2dt_b1259 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_2dt_b 1260 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_2dt_b 1261 1261 END DO 1262 1262 ! Save barotropic velocities not transport: … … 1268 1268 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1269 1269 DO jk = 1, jpkm1 1270 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu _n(:,:) - puu_b(:,:,Kmm) ) * umask(:,:,jk)1271 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv _n(:,:) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)1270 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk) 1271 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk) 1272 1272 END DO 1273 1273 1274 1274 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 1275 1275 DO jk = 1, jpkm1 1276 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu _n(:,:) &1277 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu _n(:,:)) ) * umask(:,:,jk)1278 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv _n(:,:) &1279 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv _n(:,:)) ) * vmask(:,:,jk)1276 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 1277 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 1278 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 1279 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 1280 1280 END DO 1281 1281 END IF 1282 1282 1283 1283 1284 CALL iom_put( "ubar", un_adv(:,:)*r1_hu _n(:,:) ) ! barotropic i-current1285 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv _n(:,:) ) ! barotropic i-current1284 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 1285 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current 1286 1286 ! 1287 1287 #if defined key_agrif -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
r11421 r11480 9 9 !! - ! 2010-09 (D.Storkey and E.O'Dea) bug fixes for BDY module 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 11 12 !!---------------------------------------------------------------------- 12 13 13 14 !!---------------------------------------------------------------------- 14 15 !! ssh_nxt : after ssh 15 !! ssh_ swp : filter ans swapthe ssh arrays16 !! ssh_atf : time filter the ssh arrays 16 17 !! wzv : compute now vertical velocity 17 18 !!---------------------------------------------------------------------- … … 43 44 PUBLIC wzv ! called by step.F90 44 45 PUBLIC wAimp ! called by step.F90 45 PUBLIC ssh_ swp! called by step.F9046 PUBLIC ssh_atf ! called by step.F90 46 47 47 48 !! * Substitutions … … 218 219 219 220 220 SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa ) 221 !!---------------------------------------------------------------------- 222 !! *** ROUTINE ssh_nxt *** 223 !! 224 !! ** Purpose : achieve the sea surface height time stepping by 225 !! applying Asselin time filter and swapping the arrays 226 !! ssh(:,:,Kaa) already computed in ssh_nxt 221 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh ) 222 !!---------------------------------------------------------------------- 223 !! *** ROUTINE ssh_atf *** 224 !! 225 !! ** Purpose : Apply Asselin time filter to now SSH. 227 226 !! 228 227 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 229 228 !! from the filter, see Leclair and Madec 2010) and swap : 230 !! ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) +ssh(:,:,Kaa) )229 !! pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 231 230 !! - atfp * rdt * ( emp_b - emp ) / rau0 232 !! ssh(:,:,Kmm) = ssh(:,:,Kaa) 233 !! 234 !! ** action : - ssh(:,:,Kbb), ssh(:,:,Kmm) : before & now sea surface height 235 !! ready for the next time step 231 !! 232 !! ** action : - pssh(:,:,Kmm) time filtered 236 233 !! 237 234 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 238 235 !!---------------------------------------------------------------------- 239 INTEGER, INTENT(in) :: kt ! ocean time-step index 240 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time-step index 236 INTEGER , INTENT(in ) :: kt ! ocean time-step index 237 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 238 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field 241 239 ! 242 240 REAL(wp) :: zcoef ! local scalar 243 241 !!---------------------------------------------------------------------- 244 242 ! 245 IF( ln_timing ) CALL timing_start('ssh_ swp')243 IF( ln_timing ) CALL timing_start('ssh_atf') 246 244 ! 247 245 IF( kt == nit000 ) THEN 248 246 IF(lwp) WRITE(numout,*) 249 IF(lwp) WRITE(numout,*) 'ssh_ swp : Asselin time filter and swapof sea surface height'247 IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height' 250 248 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 251 249 ENDIF 252 250 ! !== Euler time-stepping: no filter, just swap ==! 253 IF ( neuler == 0 .AND. kt == nit000 ) THEN 254 ssh(:,:,Kmm) = ssh(:,:,Kaa) ! now <-- after (before already = now) 255 ! 256 ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==! 257 ! ! before <-- now filtered 258 ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 259 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 251 IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN ! Only do time filtering for leapfrog timesteps 252 ! ! filtered "now" field 253 pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 254 IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed 260 255 zcoef = atfp * rdt * r1_rau0 261 ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * ( emp_b(:,:) - emp (:,:) &256 pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) & 262 257 & - rnf_b(:,:) + rnf (:,:) & 263 258 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 264 259 ENDIF 265 ssh(:,:,Kmm) = ssh(:,:,Kaa) ! now <-- after 266 ENDIF 267 ! 268 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssh(:,:,Kbb), clinfo1=' ssh(:,:,Kbb) - : ', mask1=tmask ) 269 ! 270 IF( ln_timing ) CALL timing_stop('ssh_swp') 271 ! 272 END SUBROUTINE ssh_swp 260 ENDIF 261 ! 262 IF(ln_ctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm) - : ', mask1=tmask ) 263 ! 264 IF( ln_timing ) CALL timing_stop('ssh_atf') 265 ! 266 END SUBROUTINE ssh_atf 273 267 274 268 SUBROUTINE wAimp( kt, Kmm ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/fldread.F90
r10922 r11480 130 130 CONTAINS 131 131 132 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl )132 SUBROUTINE fld_read( kt, kn_fsbc, sd, map, kit, kt_offset, jpk_bdy, fvl, Kmm ) 133 133 !!--------------------------------------------------------------------- 134 134 !! *** ROUTINE fld_read *** … … 153 153 INTEGER , INTENT(in ), OPTIONAL :: jpk_bdy ! number of vertical levels in the BDY data 154 154 LOGICAL , INTENT(in ), OPTIONAL :: fvl ! number of vertical levels in the BDY data 155 INTEGER , INTENT(in ), OPTIONAL :: Kmm ! ocean time level index 155 156 !! 156 157 INTEGER :: itmp ! local variable … … 287 288 ! read after data 288 289 IF( PRESENT(jpk_bdy) ) THEN 289 CALL fld_get( sd(jf), imap, jpk_bdy, fvl )290 CALL fld_get( sd(jf), imap, jpk_bdy, fvl, Kmm ) 290 291 ELSE 291 292 CALL fld_get( sd(jf), imap ) … … 614 615 615 616 616 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl )617 SUBROUTINE fld_get( sdjf, map, jpk_bdy, fvl, Kmm ) 617 618 !!--------------------------------------------------------------------- 618 619 !! *** ROUTINE fld_get *** … … 624 625 INTEGER , INTENT(in), OPTIONAL :: jpk_bdy ! number of vertical levels in the bdy data 625 626 LOGICAL , INTENT(in), OPTIONAL :: fvl ! number of vertical levels in the bdy data 627 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 626 628 ! 627 629 INTEGER :: ipk ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) … … 638 640 IF( PRESENT(jpk_bdy) ) THEN 639 641 IF( sdjf%ln_tint ) THEN ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), & 640 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )642 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 641 643 ELSE ; CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,: ), & 642 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl )644 sdjf%nrec_a(1), map, sdjf%igrd, sdjf%ibdy, jpk_bdy, fvl, Kmm ) 643 645 ENDIF 644 646 ELSE … … 701 703 END SUBROUTINE fld_get 702 704 703 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl )705 SUBROUTINE fld_map( num, clvar, dta, nrec, map, igrd, ibdy, jpk_bdy, fvl, Kmm ) 704 706 !!--------------------------------------------------------------------- 705 707 !! *** ROUTINE fld_map *** … … 718 720 INTEGER , INTENT(in), OPTIONAL :: igrd, ibdy, jpk_bdy ! grid type, set number and number of vertical levels in the bdy data 719 721 LOGICAL , INTENT(in), OPTIONAL :: fvl ! grid type, set number and number of vertical levels in the bdy data 722 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 720 723 INTEGER :: jpkm1_bdy! number of vertical levels in the bdy data minus 1 721 724 !! … … 813 816 814 817 IF ( ln_bdy ) & 815 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta )818 CALL fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 816 819 817 820 ELSE ! boundary data assumed to be on model grid … … 838 841 END SUBROUTINE fld_map 839 842 840 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta )843 SUBROUTINE fld_bdy_interp(dta_read, dta_read_z, dta_read_dz, map, jpk_bdy, igrd, ibdy, fv, dta, fvl, ilendta, Kmm) 841 844 842 845 !!--------------------------------------------------------------------- … … 857 860 INTEGER , INTENT(in) :: igrd, ibdy, jpk_bdy ! number of levels in bdy data 858 861 INTEGER , INTENT(in) :: ilendta ! length of data in file 862 INTEGER , INTENT(in), OPTIONAL :: Kmm ! ocean time level index 859 863 !! 860 864 INTEGER :: ipi ! length of boundary data on local process … … 900 904 SELECT CASE( igrd ) 901 905 CASE(1) 902 IF( ABS( (zh - ht _n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN906 IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 903 907 WRITE(ibstr,"(I10.10)") map%ptr(ib) 904 908 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 905 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:, nfld_Nnn), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj909 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1), ht(zij,zjj), map%ptr(ib), ib, zij, zjj 906 910 ENDIF 907 911 CASE(2) 908 IF( ABS( (zh - hu _n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN912 IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 909 913 WRITE(ibstr,"(I10.10)") map%ptr(ib) 910 914 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 911 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:, nfld_Nnn), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), &912 & hu _n(zij,zjj), map%ptr(ib), ib, zij, zjj, narea-1 , &915 IF(lwp) WRITE(*,*) 'DEPTHU', zh, sum(e3u(zij,zjj,:,Kmm), mask=umask(zij,zjj,:)==1), sum(umask(zij,zjj,:)), & 916 & hu(zij,zjj,Kmm), map%ptr(ib), ib, zij, zjj, narea-1 , & 913 917 & dta_read(map%ptr(ib),1,:) 914 918 ENDIF 915 919 CASE(3) 916 IF( ABS( (zh - hv _n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN920 IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 917 921 WRITE(ibstr,"(I10.10)") map%ptr(ib) 918 922 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') … … 922 926 SELECT CASE( igrd ) 923 927 CASE(1) 924 zl = gdept(zij,zjj,ik, nfld_Nnn) ! if using in step could use fsdept instead of gdept_n?928 zl = gdept(zij,zjj,ik,Kmm) ! if using in step could use fsdept instead of gdept_n? 925 929 CASE(2) 926 930 IF(ln_sco) THEN 927 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?931 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 928 932 ELSE 929 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) )933 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) ) 930 934 ENDIF 931 935 CASE(3) 932 936 IF(ln_sco) THEN 933 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?937 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 934 938 ELSE 935 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )939 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 936 940 ENDIF 937 941 END SELECT … … 941 945 dta(ib,1,ik) = dta_read(map%ptr(ib),1,MAXLOC(dta_read_z(map%ptr(ib),1,:),1)) 942 946 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 943 DO ikk = 1, jpkm1_bdy ! when gdept(ikk, nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)947 DO ikk = 1, jpkm1_bdy ! when gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 944 948 IF( ( (zl-dta_read_z(map%ptr(ib),1,ikk)) * (zl-dta_read_z(map%ptr(ib),1,ikk+1)) <= 0._wp) & 945 949 & .AND. (dta_read_z(map%ptr(ib),1,ikk+1) /= fv_alt)) THEN … … 965 969 ENDDO 966 970 DO ik = 1, ipk ! calculate transport on model grid 967 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik, nfld_Nnn) * umask(zij,zjj,ik)971 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 968 972 ENDDO 969 973 DO ik = 1, ipk ! make transport correction 970 974 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 971 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)975 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 972 976 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 973 IF( ABS(ztrans * r1_hu _n(zij,zjj)) > 0.01_wp ) &977 IF( ABS(ztrans * r1_hu(zij,zjj,Kmm)) > 0.01_wp ) & 974 978 & CALL ctl_warn('fld_bdy_interp: barotropic component of > 0.01 ms-1 found in baroclinic velocities at') 975 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu _n(zij,zjj) * umask(zij,zjj,ik)979 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu(zij,zjj,Kmm) * umask(zij,zjj,ik) 976 980 ENDIF 977 981 ENDDO … … 990 994 ENDDO 991 995 DO ik = 1, ipk ! calculate transport on model grid 992 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik, nfld_Nnn) * vmask(zij,zjj,ik)996 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 993 997 ENDDO 994 998 DO ik = 1, ipk ! make transport correction 995 999 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 996 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1000 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 997 1001 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 998 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv _n(zij,zjj) * vmask(zij,zjj,ik)1002 dta(ib,1,ik) = dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv(zij,zjj,Kmm) * vmask(zij,zjj,ik) 999 1003 ENDIF 1000 1004 ENDDO … … 1025 1029 SELECT CASE( igrd ) 1026 1030 CASE(1) 1027 IF( ABS( (zh - ht _n(zij,zjj)) / ht_n(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN1031 IF( ABS( (zh - ht(zij,zjj)) / ht(zij,zjj)) * tmask(zij,zjj,1) > 0.01_wp ) THEN 1028 1032 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1029 1033 CALL ctl_warn('fld_bdy_interp: T depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1030 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:, nfld_Nnn), mask=tmask(zij,zjj,:)==1), ht_n(zij,zjj), map%ptr(ib), ib, zij, zjj1034 ! IF(lwp) WRITE(*,*) 'DEPTHT', zh, sum(e3t(zij,zjj,:,Kmm), mask=tmask(zij,zjj,:)==1), ht(zij,zjj), map%ptr(ib), ib, zij, zjj 1031 1035 ENDIF 1032 1036 CASE(2) 1033 IF( ABS( (zh - hu _n(zij,zjj)) * r1_hu_n(zij,zjj)) * umask(zij,zjj,1) > 0.01_wp ) THEN1037 IF( ABS( (zh - hu(zij,zjj,Kmm)) * r1_hu(zij,zjj,Kmm)) * umask(zij,zjj,1) > 0.01_wp ) THEN 1034 1038 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1035 1039 CALL ctl_warn('fld_bdy_interp: U depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') 1036 1040 ENDIF 1037 1041 CASE(3) 1038 IF( ABS( (zh - hv _n(zij,zjj)) * r1_hv_n(zij,zjj)) * vmask(zij,zjj,1) > 0.01_wp ) THEN1042 IF( ABS( (zh - hv(zij,zjj,Kmm)) * r1_hv(zij,zjj,Kmm)) * vmask(zij,zjj,1) > 0.01_wp ) THEN 1039 1043 WRITE(ibstr,"(I10.10)") map%ptr(ib) 1040 1044 CALL ctl_warn('fld_bdy_interp: V depths differ between grids at BDY point '//TRIM(ibstr)//' by more than 1%') … … 1044 1048 SELECT CASE( igrd ) ! coded for sco - need zco and zps option using min 1045 1049 CASE(1) 1046 zl = gdept(zij,zjj,ik, nfld_Nnn) ! if using in step could use fsdept instead of gdept_n?1050 zl = gdept(zij,zjj,ik,Kmm) ! if using in step could use fsdept instead of gdept_n? 1047 1051 CASE(2) 1048 1052 IF(ln_sco) THEN 1049 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij+1,zjj,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?1053 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij+1,zjj,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1050 1054 ELSE 1051 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij+1,zjj,ik,nfld_Nnn) )1055 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij+1,zjj,ik,Kmm) ) 1052 1056 ENDIF 1053 1057 CASE(3) 1054 1058 IF(ln_sco) THEN 1055 zl = ( gdept(zij,zjj,ik, nfld_Nnn) + gdept(zij,zjj+1,ik,nfld_Nnn) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n?1059 zl = ( gdept(zij,zjj,ik,Kmm) + gdept(zij,zjj+1,ik,Kmm) ) * 0.5_wp ! if using in step could use fsdept instead of gdept_n? 1056 1060 ELSE 1057 zl = MIN( gdept(zij,zjj,ik, nfld_Nnn), gdept(zij,zjj+1,ik,nfld_Nnn) )1061 zl = MIN( gdept(zij,zjj,ik,Kmm), gdept(zij,zjj+1,ik,Kmm) ) 1058 1062 ENDIF 1059 1063 END SELECT … … 1063 1067 dta(ib,1,ik) = dta_read(ji,jj,MAXLOC(dta_read_z(ji,jj,:),1)) 1064 1068 ELSE ! inbetween : vertical interpolation between ikk & ikk+1 1065 DO ikk = 1, jpkm1_bdy ! when gdept(ikk, nfld_Nnn) < zl < gdept(ikk+1,nfld_Nnn)1069 DO ikk = 1, jpkm1_bdy ! when gdept(ikk,Kmm) < zl < gdept(ikk+1,Kmm) 1066 1070 IF( ( (zl-dta_read_z(ji,jj,ikk)) * (zl-dta_read_z(ji,jj,ikk+1)) <= 0._wp) & 1067 1071 & .AND. (dta_read_z(ji,jj,ikk+1) /= fv_alt)) THEN … … 1089 1093 ENDDO 1090 1094 DO ik = 1, ipk ! calculate transport on model grid 1091 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik, nfld_Nnn) * umask(zij,zjj,ik)1095 ztrans_new = ztrans_new + dta(ib,1,ik) * e3u(zij,zjj,ik,Kmm) * umask(zij,zjj,ik) 1092 1096 ENDDO 1093 1097 DO ik = 1, ipk ! make transport correction 1094 1098 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1095 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)1099 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 1096 1100 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1097 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu _n(zij,zjj) ) * umask(zij,zjj,ik)1101 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hu(zij,zjj,Kmm) ) * umask(zij,zjj,ik) 1098 1102 ENDIF 1099 1103 ENDDO … … 1114 1118 ENDDO 1115 1119 DO ik = 1, ipk ! calculate transport on model grid 1116 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik, nfld_Nnn) * vmask(zij,zjj,ik)1120 ztrans_new = ztrans_new + dta(ib,1,ik) * e3v(zij,zjj,ik,Kmm) * vmask(zij,zjj,ik) 1117 1121 ENDDO 1118 1122 DO ik = 1, ipk ! make transport correction 1119 1123 IF(fvl) THEN ! bdy data are total velocity so adjust bt transport term to match input data 1120 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1124 dta(ib,1,ik) = ( dta(ib,1,ik) + ( ztrans - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 1121 1125 ELSE ! we're just dealing with bc velocity so bt transport term should sum to zero 1122 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv _n(zij,zjj) ) * vmask(zij,zjj,ik)1126 dta(ib,1,ik) = ( dta(ib,1,ik) + ( 0._wp - ztrans_new ) * r1_hv(zij,zjj,Kmm) ) * vmask(zij,zjj,ik) 1123 1127 ENDIF 1124 1128 ENDDO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/SBC/sbcmod.F90
r11027 r11480 442 442 CASE( 1 ) ; CALL sbc_ice_if ( kt, Kbb, Kmm ) ! Ice-cover climatology ("Ice-if" model) 443 443 #if defined key_si3 444 CASE( 2 ) ; CALL ice_stp ( kt, Kbb, nsbc )! SI3 ice model444 CASE( 2 ) ; CALL ice_stp ( kt, Kbb, Kmm, nsbc ) ! SI3 ice model 445 445 #endif 446 446 CASE( 3 ) ; CALL sbc_ice_cice ( kt, nsbc ) ! CICE ice model -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/trasbc.F90
r10985 r11480 233 233 DO jj = 2, jpj 234 234 DO ji = fs_2, fs_jpim1 235 ztim = ssh_iau(ji,jj) / ( ht _n(ji,jj) + 1. - ssmask(ji, jj) )235 ztim = ssh_iau(ji,jj) / ( ht(ji,jj) + 1. - ssmask(ji, jj) ) 236 236 pts(ji,jj,:,jp_tem,Krhs) = pts(ji,jj,:,jp_tem,Krhs) + pts(ji,jj,:,jp_tem,Kmm) * ztim 237 237 pts(ji,jj,:,jp_sal,Krhs) = pts(ji,jj,:,jp_sal,Krhs) + pts(ji,jj,:,jp_sal,Kmm) * ztim -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trddyn.F90
r10946 r11480 123 123 z3dx(:,:,:) = 0._wp ! U.dxU & V.dyV (approximation) 124 124 z3dy(:,:,:) = 0._wp 125 DO jk = 1, jpkm1 ! no mask as u n,vnare masked125 DO jk = 1, jpkm1 ! no mask as uu, vv are masked 126 126 DO jj = 2, jpjm1 127 127 DO ji = 2, jpim1 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRD/trdvor.F90
r10946 r11480 189 189 190 190 ! Average except for Beta.V 191 zudpvor(:,:) = zudpvor(:,:) * r1_hu _n(:,:)192 zvdpvor(:,:) = zvdpvor(:,:) * r1_hv _n(:,:)191 zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm) 192 zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm) 193 193 194 194 ! Curl … … 276 276 END DO 277 277 ! Average of the Curl and Surface mask 278 vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * r1_hu _n(:,:) * fmask(:,:,1)278 vortrd(:,:,jpvor_bev) = vortrd(:,:,jpvor_bev) * r1_hu(:,:,Kmm) * fmask(:,:,1) 279 279 ENDIF 280 280 ! 281 281 ! Average 282 zudpvor(:,:) = zudpvor(:,:) * r1_hu _n(:,:)283 zvdpvor(:,:) = zvdpvor(:,:) * r1_hv _n(:,:)282 zudpvor(:,:) = zudpvor(:,:) * r1_hu(:,:,Kmm) 283 zvdpvor(:,:) = zvdpvor(:,:) * r1_hv(:,:,Kmm) 284 284 ! 285 285 ! Curl … … 342 342 END DO 343 343 344 zuu(:,:) = zuu(:,:) * r1_hu _n(:,:)345 zvv(:,:) = zvv(:,:) * r1_hv _n(:,:)344 zuu(:,:) = zuu(:,:) * r1_hu(:,:,Kmm) 345 zvv(:,:) = zvv(:,:) * r1_hv(:,:,Kmm) 346 346 347 347 ! Curl -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/ZDF/zdfosm.F90
r10955 r11480 489 489 490 490 zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it 491 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht _n(:,:))491 zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:)) 492 492 zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom 493 493 … … 525 525 526 526 zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w(ji,jj,jk,Kmm) ) 527 zhbl_s = MIN(zhbl_s, ht _n(ji,jj))527 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 528 528 529 529 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1 … … 546 546 & * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj) ! ALMG to investigate whether need to include ww here 547 547 548 zhbl_s = MIN(zhbl_s, ht _n(ji,jj))548 zhbl_s = MIN(zhbl_s, ht(ji,jj)) 549 549 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1 550 550 END DO -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/nemogcm.F90
r11421 r11480 473 473 #if defined key_top 474 474 ! ! Passive tracers 475 CALL trc_init ( Nbb, Nnn, Naa )475 CALL trc_init 476 476 #endif 477 477 IF( l_ldfslp ) CALL ldf_slp_init ! slope of lateral mixing -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/oce.F90
r11047 r11480 8 8 !! 3.3 ! 2010-09 (C. Ethe) TRA-TRC merge: add ts, gtsu, gtsv 4D arrays 9 9 !! 3.7 ! 2014-01 (G. Madec) suppression of curl and before hdiv from in-core memory 10 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rename prognostic variables in preparation for new time scheme 10 11 !!---------------------------------------------------------------------- 11 12 USE par_oce ! ocean parameters -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step.F90
r11421 r11480 31 31 !! - ! 2015-11 (J. Chanut) free surface simplification (remove filtered free surface) 32 32 !! 4.0 ! 2017-05 (G. Madec) introduction of the vertical physics manager (zdfphy) 33 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme 33 34 !!---------------------------------------------------------------------- 34 35 … … 265 266 266 267 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 267 ! Set boundary conditions and Swap268 ! Set boundary conditions, time filter and swap time levels 268 269 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 269 270 !!jc1: For agrif, it would be much better to finalize tracers/momentum here (e.g. bdy conditions) and move the swap … … 281 282 !! 282 283 !!jc2: dynnxt must be the latest call. e3t(:,:,:,Nbb) are indeed updated in that routine 283 CALL tra_nxt ( kstp, Nbb, Nnn, Nrhs, Naa ) ! finalize (bcs) tracer fields at next time step and swap 284 CALL dyn_nxt ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa ) ! finalize (bcs) velocities at next time step and swap (always called after tra_nxt) 285 CALL ssh_swp ( kstp, Nbb, Nnn, Naa ) ! swap of sea surface height 286 IF(.NOT.ln_linssh) CALL dom_vvl_sf_swp( kstp, Nbb, Nnn, Naa ) ! swap of vertical scale factors 284 CALL tra_atf ( kstp, Nbb, Nnn, Naa, ts ) ! time filtering of "now" tracer arrays 285 CALL dyn_atf ( kstp, Nbb, Nnn, Naa, uu, vv, e3t, e3u, e3v ) ! time filtering of "now" velocities and scale factors 286 CALL ssh_atf ( kstp, Nbb, Nnn, Naa, ssh ) ! time filtering of "now" sea surface height 287 ! 288 ! Swap time levels 289 Nrhs = Nbb 290 Nbb = Nnn 291 Nnn = Naa 292 Naa = Nrhs 293 ! 294 IF(.NOT.ln_linssh) CALL dom_vvl_sf_update( kstp, Nbb, Nnn, Naa ) ! recompute vertical scale factors 287 295 ! 288 296 IF( ln_diahsb ) CALL dia_hsb ( kstp, Nbb, Nnn ) ! - ML - global conservation diagnostics … … 300 308 ! AGRIF 301 309 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 310 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices 302 311 CALL Agrif_Integrate_ChildGrids( stp ) ! allows to finish all the Child Grids before updating 303 312 304 313 IF( Agrif_NbStepint() == 0 ) THEN 305 Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices306 314 CALL Agrif_update_all( ) ! Update all components 307 315 ENDIF … … 312 320 ! Control 313 321 !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< 314 CALL stp_ctl ( kstp, N nn, indic )322 CALL stp_ctl ( kstp, Nbb, Nnn, indic ) 315 323 316 324 IF( kstp == nit000 ) THEN ! 1st time step only … … 337 345 ! 338 346 END SUBROUTINE stp 339 347 ! 340 348 !!====================================================================== 341 349 END MODULE step -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/step_oce.F90
r10927 r11480 30 30 USE traldf ! lateral mixing (tra_ldf routine) 31 31 USE trazdf ! vertical mixing (tra_zdf routine) 32 USE tra nxt ! time-stepping (tra_nxtroutine)32 USE traatf ! time filtering (tra_atf routine) 33 33 USE tranpc ! non-penetrative convection (tra_npc routine) 34 34 … … 43 43 USE dynspg ! surface pressure gradient (dyn_spg routine) 44 44 45 USE dyn nxt ! time-stepping (dyn_nxtroutine)45 USE dynatf ! time-filtering (dyn_atf routine) 46 46 47 47 USE stopar ! Stochastic parametrization (sto_par routine) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/stpctl.F90
r10965 r11480 42 42 CONTAINS 43 43 44 SUBROUTINE stp_ctl( kt, K mm, kindic )44 SUBROUTINE stp_ctl( kt, Kbb, Kmm, kindic ) 45 45 !!---------------------------------------------------------------------- 46 46 !! *** ROUTINE stp_ctl *** … … 60 60 !!---------------------------------------------------------------------- 61 61 INTEGER, INTENT(in ) :: kt ! ocean time-step index 62 INTEGER, INTENT(in ) :: K mm ! ocean time level index62 INTEGER, INTENT(in ) :: Kbb, Kmm ! ocean time level index 63 63 INTEGER, INTENT(inout) :: kindic ! error indicator 64 64 !! -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/dtadyn.F90
r11027 r11480 46 46 PRIVATE 47 47 48 PUBLIC dta_dyn_init ! called by opa.F90 49 PUBLIC dta_dyn ! called by step.F90 50 PUBLIC dta_dyn_sed_init ! called by opa.F90 51 PUBLIC dta_dyn_sed ! called by step.F90 52 PUBLIC dta_dyn_swp ! called by step.F90 48 PUBLIC dta_dyn_init ! called by nemo_init 49 PUBLIC dta_dyn ! called by nemo_gcm 50 PUBLIC dta_dyn_sed_init ! called by nemo_init 51 PUBLIC dta_dyn_sed ! called by nemo_gcm 52 PUBLIC dta_dyn_atf ! called by nemo_gcm 53 PUBLIC dta_dyn_sf_interp ! called by nemo_gcm 53 54 54 55 CHARACTER(len=100) :: cn_dir !: Root directory for location of ssr files … … 535 536 END SUBROUTINE dta_dyn_sed_init 536 537 537 SUBROUTINE dta_dyn_ swp( kt, Kbb, Kmm, Kaa )538 SUBROUTINE dta_dyn_atf( kt, Kbb, Kmm, Kaa ) 538 539 !!--------------------------------------------------------------------- 539 540 !! *** ROUTINE dta_dyn_swp *** 540 541 !! 541 !! ** Purpose : Swap and the data and compute the vertical scale factor 542 !! at U/V/W pointand the depht 542 !! ** Purpose : Asselin time filter of now SSH 543 543 !!--------------------------------------------------------------------- 544 544 INTEGER, INTENT(in) :: kt ! time step 545 545 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time level indices 546 546 ! 547 !!--------------------------------------------------------------------- 548 549 IF( kt == nit000 ) THEN 550 IF(lwp) WRITE(numout,*) 551 IF(lwp) WRITE(numout,*) 'dta_dyn_atf : Asselin time filter of sea surface height' 552 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ ' 553 ENDIF 554 555 ssh(:,:,Kmm) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) 556 557 !! Do we also need to time filter e3t?? 558 ! 559 END SUBROUTINE dta_dyn_atf 560 561 SUBROUTINE dta_dyn_sf_interp( kt, Kmm ) 562 !!--------------------------------------------------------------------- 563 !! *** ROUTINE dta_dyn_sf_interp *** 564 !! 565 !! ** Purpose : Calculate scale factors at U/V/W points and depths 566 !! given the after e3t field 567 !!--------------------------------------------------------------------- 568 INTEGER, INTENT(in) :: kt ! time step 569 INTEGER, INTENT(in) :: Kmm ! ocean time level indices 570 ! 547 571 INTEGER :: ji, jj, jk 548 572 REAL(wp) :: zcoef 549 573 !!--------------------------------------------------------------------- 550 551 IF( kt == nit000 ) THEN552 IF(lwp) WRITE(numout,*)553 IF(lwp) WRITE(numout,*) 'ssh_swp : Asselin time filter and swap of sea surface height'554 IF(lwp) WRITE(numout,*) '~~~~~~~ '555 ENDIF556 557 ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa)) ! before <-- now filtered558 ssh(:,:,Kmm) = ssh(:,:,Kaa)559 560 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)561 562 ! Reconstruction of all vertical scale factors at now and before time steps563 ! =============================================================================564 574 565 575 ! Horizontal scale factor interpolations … … 571 581 ! ------------------------------------ 572 582 CALL dom_vvl_interpol( e3t(:,:,:,Kmm), e3w (:,:,:,Kmm), 'W' ) 573 574 e3t(:,:,:,Kbb) = e3t(:,:,:,Kmm)575 e3u(:,:,:,Kbb) = e3u(:,:,:,Kmm)576 e3v(:,:,:,Kbb) = e3v(:,:,:,Kmm)577 583 578 584 ! t- and w- points depth … … 592 598 END DO 593 599 ! 594 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm) 595 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm) 596 ! 597 END SUBROUTINE dta_dyn_swp 598 600 END SUBROUTINE dta_dyn_sf_interp 599 601 600 602 SUBROUTINE dta_dyn_ssh( kt, phdivtr, psshb, pemp, pssha, pe3ta ) -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OFF/nemogcm.F90
r11027 r11480 119 119 #else 120 120 CALL dta_dyn ( istp, Nbb, Nnn, Naa ) ! Interpolation of the dynamical fields 121 IF( .NOT.ln_linssh ) CALL dta_dyn_swp( istp, Nbb, Nnn, Naa ) ! swap of sea surface height and vertical scale factors122 121 #endif 123 122 CALL trc_stp ( istp, Nbb, Nnn, Nrhs, Naa ) ! time-stepping 123 #if ! defined key_sed_off 124 IF( .NOT.ln_linssh ) CALL dta_dyn_atf( istp, Nbb, Nnn, Naa ) ! time filter of sea surface height and vertical scale factors 125 #endif 126 ! Swap time levels 127 Nrhs = Nbb 128 Nbb = Nnn 129 Nnn = Naa 130 Naa = Nrhs 131 ! 132 #if ! defined key_sed_off 133 IF( .NOT.ln_linssh ) CALL dta_dyn_sf_interp( istp, Nnn ) ! calculate now grid parameters 134 #endif 124 135 CALL stp_ctl ( istp, indic ) ! Time loop: control and print 125 136 istp = istp + 1 … … 333 344 #endif 334 345 335 CALL trc_init ( Nbb, Nnn, Naa )! Passive tracers initialization346 CALL trc_init ! Passive tracers initialization 336 347 CALL dia_ptr_init ! Poleward TRansports initialization 337 348 -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/PISCES/SED/trcdmp_sed.F90
r10975 r11480 149 149 !!---------------------------------------------------------------------- 150 150 CONTAINS 151 SUBROUTINE trc_dmp_sed( kt )! Empty routine151 SUBROUTINE trc_dmp_sed( kt, Kbb, Kmm, Krhs ) ! Empty routine 152 152 INTEGER, INTENT(in) :: kt 153 INTEGER, INTENT(in) :: Kbb, Kmm, Krhs 153 154 WRITE(*,*) 'trc_dmp_sed: You should not have seen this print! error?', kt 154 155 END SUBROUTINE trc_dmp_sed -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcrad.F90
r10985 r11480 37 37 CONTAINS 38 38 39 SUBROUTINE trc_rad( kt, Kbb, Kmm, Krhs,ptr )39 SUBROUTINE trc_rad( kt, Kbb, Kmm, ptr ) 40 40 !!---------------------------------------------------------------------- 41 41 !! *** ROUTINE trc_rad *** … … 52 52 !! (the total CFC content is not strictly preserved) 53 53 !!---------------------------------------------------------------------- 54 INTEGER, INTENT(in ) :: kt 55 INTEGER, INTENT(in ) :: Kbb, Kmm , Krhs! time level indices56 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr 54 INTEGER, INTENT(in ) :: kt ! ocean time-step index 55 INTEGER, INTENT(in ) :: Kbb, Kmm ! time level indices 56 REAL(wp), DIMENSION(jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! passive tracers and RHS of tracer equation 57 57 ! 58 58 CHARACTER (len=22) :: charout … … 61 61 IF( ln_timing ) CALL timing_start('trc_rad') 62 62 ! 63 IF( ln_age ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_age , jp_age ) ! AGE64 IF( ll_cfc ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_cfc0, jp_cfc1 ) ! CFC model65 IF( ln_c14 ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_c14 , jp_c14 ) ! C1466 IF( ln_pisces ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model67 IF( ln_my_trc ) CALL trc_rad_sms( kt, K mm, Krhs, ptr(:,:,:,:,Kbb), ptr(:,:,:,:,Kmm), jp_myt0, jp_myt1 ) ! MY_TRC model63 IF( ln_age ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_age , jp_age ) ! AGE 64 IF( ll_cfc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_cfc0, jp_cfc1 ) ! CFC model 65 IF( ln_c14 ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_c14 , jp_c14 ) ! C14 66 IF( ln_pisces ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_pcs0, jp_pcs1, cpreserv='Y' ) ! PISCES model 67 IF( ln_my_trc ) CALL trc_rad_sms( kt, Kbb, Kmm, ptr, jp_myt0, jp_myt1 ) ! MY_TRC model 68 68 ! 69 69 IF(ln_ctl) THEN ! print mean trends (used for debugging) 70 70 WRITE(charout, FMT="('rad')") 71 71 CALL prt_ctl_trc_info( charout ) 72 CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,K mm), mask=tmask, clinfo=ctrcnm )72 CALL prt_ctl_trc( tab4d=ptr(:,:,:,:,Kbb), mask=tmask, clinfo=ctrcnm ) 73 73 ENDIF 74 74 ! … … 115 115 116 116 117 SUBROUTINE trc_rad_sms( kt, Kmm, Krhs, ptrb, ptrn, jp_sms0, jp_sms1, cpreserv ) 118 !!----------------------------------------------------------------------------- 119 !! *** ROUTINE trc_rad_sms *** 120 !! 121 !! ** Purpose : "crappy" routine to correct artificial negative 122 !! concentrations due to isopycnal scheme 123 !! 124 !! ** Method : 2 cases : 125 !! - Set negative concentrations to zero while computing 126 !! the corresponding tracer content that is added to the 127 !! tracers. Then, adjust the tracer concentration using 128 !! a multiplicative factor so that the total tracer 129 !! concentration is preserved. 130 !! - simply set to zero the negative CFC concentration 131 !! (the total content of concentration is not strictly preserved) 132 !!-------------------------------------------------------------------------------- 133 INTEGER , INTENT(in ) :: kt ! ocean time-step index 134 INTEGER , INTENT(in ) :: Kmm, Krhs ! time level indices 135 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 136 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra), INTENT(inout) :: ptrb , ptrn ! before and now traceur concentration 137 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 138 ! 139 INTEGER :: ji, ji2, jj, jj2, jk, jn ! dummy loop indices 140 INTEGER :: icnt 141 LOGICAL :: lldebug = .FALSE. ! local logical 142 REAL(wp):: zcoef, zs2rdt, ztotmass 143 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 144 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 145 !!---------------------------------------------------------------------- 146 ! 147 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 148 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 149 ! 150 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 151 ! 152 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 153 154 DO jn = jp_sms0, jp_sms1 155 ztrneg(:,:,jn) = SUM( MIN( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 156 ztrpos(:,:,jn) = SUM( MAX( 0., ptrb(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 157 END DO 158 CALL sum3x3( ztrneg ) 159 CALL sum3x3( ztrpos ) 160 161 DO jn = jp_sms0, jp_sms1 162 ! 163 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input tr(:,:,:,:,Kbb) for trend computation 164 ! 165 DO jk = 1, jpkm1 166 DO jj = 1, jpj 167 DO ji = 1, jpi 168 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 169 ! 170 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 171 IF( ptrb(ji,jj,jk,jn) < 0. ) ptrb(ji,jj,jk,jn) = 0. ! supress negative values 172 IF( ptrb(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 173 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 174 ptrb(ji,jj,jk,jn) = ptrb(ji,jj,jk,jn) * zcoef 175 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 176 gainmass(jn,1) = gainmass(jn,1) - ptrb(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 177 ptrb(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 178 ENDIF 179 ENDIF 180 ! 181 ENDIF 182 END DO 183 END DO 184 END DO 185 ! 186 IF( l_trdtrc ) THEN 187 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 188 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 189 ENDIF 190 ! 191 END DO 192 193 IF( kt == nitend ) THEN 194 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 195 DO jn = jp_sms0, jp_sms1 196 IF( gainmass(jn,1) > 0. ) THEN 197 ztotmass = glob_sum( 'trcrad', ptrb(:,:,:,jn) * cvol(:,:,:) ) 198 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 199 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 200 END IF 201 END DO 202 ENDIF 203 204 DO jn = jp_sms0, jp_sms1 205 ztrneg(:,:,jn) = SUM( MIN( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 206 ztrpos(:,:,jn) = SUM( MAX( 0., ptrn(:,:,:,jn) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 207 END DO 208 CALL sum3x3( ztrneg ) 209 CALL sum3x3( ztrpos ) 210 211 DO jn = jp_sms0, jp_sms1 212 ! 213 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input tr for trend computation 214 ! 215 DO jk = 1, jpkm1 216 DO jj = 1, jpj 217 DO ji = 1, jpi 218 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 219 ! 220 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * tmask(ji,jj,jk) ! really needed? 221 IF( ptrn(ji,jj,jk,jn) < 0. ) ptrn(ji,jj,jk,jn) = 0. ! supress negative values 222 IF( ptrn(ji,jj,jk,jn) > 0. ) THEN ! use positive values to compensate mass gain 223 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptrb > 0 224 ptrn(ji,jj,jk,jn) = ptrn(ji,jj,jk,jn) * zcoef 225 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 226 gainmass(jn,2) = gainmass(jn,2) - ptrn(ji,jj,jk,jn) * cvol(ji,jj,jk) ! we are adding mass... 227 ptrn(ji,jj,jk,jn) = 0. ! limit the compensation to keep positive value 228 ENDIF 229 ENDIF 230 ! 231 ENDIF 232 END DO 233 END DO 234 END DO 235 ! 236 IF( l_trdtrc ) THEN 237 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 238 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 239 ENDIF 240 ! 241 END DO 242 243 IF( kt == nitend ) THEN 244 CALL mpp_sum( 'trcrad', gainmass(:,2) ) 245 DO jn = jp_sms0, jp_sms1 246 IF( gainmass(jn,2) > 0. ) THEN 247 ztotmass = glob_sum( 'trcrad', ptrn(:,:,:,jn) * cvol(:,:,:) ) 248 WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrn, traceur ', jn & 249 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 250 END IF 251 END DO 252 ENDIF 253 254 DEALLOCATE( ztrneg, ztrpos ) 255 ! 256 ELSE !== total CFC content is NOT strictly preserved ==! 257 ! 258 DO jn = jp_sms0, jp_sms1 259 ! 260 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrb(:,:,:,jn) ! save input tr for trend computation 261 ! 262 WHERE( ptrb(:,:,:,jn) < 0. ) ptrb(:,:,:,jn) = 0. 263 ! 264 IF( l_trdtrc ) THEN 265 ztrtrd(:,:,:) = ( ptrb(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 266 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 267 ENDIF 268 ! 269 IF( l_trdtrc ) ztrtrd(:,:,:) = ptrn(:,:,:,jn) ! save input tr for trend computation 270 ! 271 WHERE( ptrn(:,:,:,jn) < 0. ) ptrn(:,:,:,jn) = 0. 272 ! 273 IF( l_trdtrc ) THEN 274 ztrtrd(:,:,:) = ( ptrn(:,:,:,jn) - ztrtrd(:,:,:) ) * zs2rdt 275 CALL trd_tra( kt, Kmm, Krhs, 'TRC', jn, jptra_radn, ztrtrd ) ! standard trend handling 276 ENDIF 277 ! 278 END DO 279 ! 280 ENDIF 117 SUBROUTINE trc_rad_sms( kt, Kbb, Kmm, ptr, jp_sms0, jp_sms1, cpreserv ) 118 !!----------------------------------------------------------------------------- 119 !! *** ROUTINE trc_rad_sms *** 120 !! 121 !! ** Purpose : "crappy" routine to correct artificial negative 122 !! concentrations due to isopycnal scheme 123 !! 124 !! ** Method : 2 cases : 125 !! - Set negative concentrations to zero while computing 126 !! the corresponding tracer content that is added to the 127 !! tracers. Then, adjust the tracer concentration using 128 !! a multiplicative factor so that the total tracer 129 !! concentration is preserved. 130 !! - simply set to zero the negative CFC concentration 131 !! (the total content of concentration is not strictly preserved) 132 !!-------------------------------------------------------------------------------- 133 INTEGER , INTENT(in ) :: kt ! ocean time-step index 134 INTEGER , INTENT(in ) :: Kbb, Kmm ! time level indices 135 INTEGER , INTENT(in ) :: jp_sms0, jp_sms1 ! First & last index of the passive tracer model 136 REAL(wp), DIMENSION (jpi,jpj,jpk,jptra,jpt), INTENT(inout) :: ptr ! before and now traceur concentration 137 CHARACTER( len = 1), OPTIONAL , INTENT(in ) :: cpreserv ! flag to preserve content or not 138 ! 139 INTEGER :: ji, ji2, jj, jj2, jk, jn, jt ! dummy loop indices 140 INTEGER :: icnt, itime 141 LOGICAL :: lldebug = .FALSE. ! local logical 142 REAL(wp):: zcoef, zs2rdt, ztotmass 143 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrneg, ztrpos 144 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrtrd ! workspace arrays 145 !!---------------------------------------------------------------------- 146 ! 147 IF( l_trdtrc ) ALLOCATE( ztrtrd(jpi,jpj,jpk) ) 148 zs2rdt = 1. / ( 2. * rdt * REAL( nn_dttrc, wp ) ) 149 ! 150 DO jt = 1,2 ! Loop over time indices since exactly the same fix is applied to "now" and "after" fields 151 IF( jt == 1 ) itime = Kbb 152 IF( jt == 2 ) itime = Kmm 153 154 IF( PRESENT( cpreserv ) ) THEN !== total tracer concentration is preserved ==! 155 ! 156 ALLOCATE( ztrneg(1:jpi,1:jpj,jp_sms0:jp_sms1), ztrpos(1:jpi,1:jpj,jp_sms0:jp_sms1) ) 157 158 DO jn = jp_sms0, jp_sms1 159 ztrneg(:,:,jn) = SUM( MIN( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the negative values 160 ztrpos(:,:,jn) = SUM( MAX( 0., ptr(:,:,:,jn,itime) ) * cvol(:,:,:), dim = 3 ) ! sum of the positive values 161 END DO 162 CALL sum3x3( ztrneg ) 163 CALL sum3x3( ztrpos ) 164 165 DO jn = jp_sms0, jp_sms1 166 ! 167 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr(:,:,:,:,Kbb) for trend computation 168 ! 169 DO jk = 1, jpkm1 170 DO jj = 1, jpj 171 DO ji = 1, jpi 172 IF( ztrneg(ji,jj,jn) /= 0. ) THEN ! if negative values over the 3x3 box 173 ! 174 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * tmask(ji,jj,jk) ! really needed? 175 IF( ptr(ji,jj,jk,jn,itime) < 0. ) ptr(ji,jj,jk,jn,itime) = 0. ! suppress negative values 176 IF( ptr(ji,jj,jk,jn,itime) > 0. ) THEN ! use positive values to compensate mass gain 177 zcoef = 1. + ztrneg(ji,jj,jn) / ztrpos(ji,jj,jn) ! ztrpos > 0 as ptr > 0 178 ptr(ji,jj,jk,jn,itime) = ptr(ji,jj,jk,jn,itime) * zcoef 179 IF( zcoef < 0. ) THEN ! if the compensation exceed the positive value 180 gainmass(jn,1) = gainmass(jn,1) - ptr(ji,jj,jk,jn,itime) * cvol(ji,jj,jk) ! we are adding mass... 181 ptr(ji,jj,jk,jn,itime) = 0. ! limit the compensation to keep positive value 182 ENDIF 183 ENDIF 184 ! 185 ENDIF 186 END DO 187 END DO 188 END DO 189 ! 190 IF( l_trdtrc ) THEN 191 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 192 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 193 ENDIF 194 ! 195 END DO 196 197 IF( kt == nitend ) THEN 198 CALL mpp_sum( 'trcrad', gainmass(:,1) ) 199 DO jn = jp_sms0, jp_sms1 200 IF( gainmass(jn,1) > 0. ) THEN 201 ztotmass = glob_sum( 'trcrad', ptr(:,:,:,jn,itime) * cvol(:,:,:) ) 202 IF(lwp) WRITE(numout, '(a, i2, a, D23.16, a, D23.16)') 'trcrad ptrb, traceur ', jn & 203 & , ' total mass : ', ztotmass, ', mass gain : ', gainmass(jn,1) 204 END IF 205 END DO 206 ENDIF 207 208 DEALLOCATE( ztrneg, ztrpos ) 209 ! 210 ELSE !== total CFC content is NOT strictly preserved ==! 211 ! 212 DO jn = jp_sms0, jp_sms1 213 ! 214 IF( l_trdtrc ) ztrtrd(:,:,:) = ptr(:,:,:,jn,itime) ! save input tr for trend computation 215 ! 216 WHERE( ptr(:,:,:,jn,itime) < 0. ) ptr(:,:,:,jn,itime) = 0. 217 ! 218 IF( l_trdtrc ) THEN 219 ztrtrd(:,:,:) = ( ptr(:,:,:,jn,itime) - ztrtrd(:,:,:) ) * zs2rdt 220 CALL trd_tra( kt, Kbb, Kmm, 'TRC', jn, jptra_radb, ztrtrd ) ! Asselin-like trend handling 221 ENDIF 222 ! 223 END DO 224 ! 225 ENDIF 226 ! 227 END DO 281 228 ! 282 229 IF( l_trdtrc ) DEALLOCATE( ztrtrd ) … … 289 236 !!---------------------------------------------------------------------- 290 237 CONTAINS 291 SUBROUTINE trc_rad( kt, Kbb, Kmm , Krhs) ! Empty routine238 SUBROUTINE trc_rad( kt, Kbb, Kmm ) ! Empty routine 292 239 INTEGER, INTENT(in) :: kt 293 INTEGER, INTENT(in) :: Kbb, Kmm , Krhs! time level indices240 INTEGER, INTENT(in) :: Kbb, Kmm ! time level indices 294 241 WRITE(*,*) 'trc_rad: You should not have seen this print! error?', kt 295 242 END SUBROUTINE trc_rad -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trcsbc.F90
r10985 r11480 197 197 !! Dummy module : NO passive tracer 198 198 !!---------------------------------------------------------------------- 199 USE par_oce 200 USE par_trc 199 201 CONTAINS 200 202 SUBROUTINE trc_sbc ( kt, Kmm, ptr, Krhs ) ! Empty routine -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/TRP/trctrp.F90
r10985 r11480 20 20 USE trcadv ! advection (trc_adv routine) 21 21 USE trczdf ! vertical diffusion (trc_zdf routine) 22 USE trc nxt ! time-stepping (trc_nxtroutine)22 USE trcatf ! time filtering (trc_atf routine) 23 23 USE trcrad ! positivity (trc_rad routine) 24 24 USE trcsbc ! surface boundary condition (trc_sbc routine) … … 53 53 !! - Update the passive tracers 54 54 !!---------------------------------------------------------------------- 55 INTEGER, INTENT( in ) :: kt ! ocean time-step index56 INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! time level indices55 INTEGER, INTENT( in ) :: kt ! ocean time-step index 56 INTEGER, INTENT( inout ) :: Kbb, Kmm, Krhs, Kaa ! TOP time level indices (swapped in this routine) 57 57 !! --------------------------------------------------------------------- 58 58 ! … … 78 78 #endif 79 79 CALL trc_zdf ( kt, Kbb, Kmm, Krhs, tr, Kaa ) ! vert. mixing & after tracer ==> after 80 CALL trc_nxt ( kt, Kbb, Kmm, Krhs ) ! tracer fields at next time step 81 IF( ln_trcrad ) CALL trc_rad ( kt, Kbb, Kmm, Krhs, tr ) ! Correct artificial negative concentrations 82 IF( ln_trcdmp_clo ) CALL trc_dmp_clo( kt, Kbb, Kmm ) ! internal damping trends on closed seas only 80 CALL trc_atf ( kt, Kbb, Kmm, Kaa , tr ) ! time filtering of "now" tracer fields 81 ! 82 ! Swap TOP time levels (= Nrhs_trc, Nbb_trc etc) 83 Krhs = Kbb 84 Kbb = Kmm 85 Kmm = Kaa 86 Kaa = Krhs 87 ! 88 IF( ln_trcrad ) CALL trc_rad ( kt, Kbb, Kmm, tr ) ! Correct artificial negative concentrations 89 IF( ln_trcdmp_clo ) CALL trc_dmp_clo( kt, Kbb, Kmm ) ! internal damping trends on closed seas only 83 90 84 91 ! … … 87 94 IF( ln_trcdmp ) CALL trc_dmp( kt, Kbb, Kmm, tr, Krhs ) ! internal damping trends 88 95 CALL trc_zdf( kt, Kbb, Kmm, Krhs, tr, Kaa ) ! vert. mixing & after tracer ==> after 89 CALL trc_nxt( kt, Kbb, Kmm, Krhs ) ! tracer fields at next time step 90 IF( ln_trcrad ) CALL trc_rad( kt, Kbb, Kmm, Krhs, tr ) ! Correct artificial negative concentrations 96 CALL trc_atf( kt, Kbb, Kmm, Kaa , tr ) ! time filtering of "now" tracer fields 97 ! 98 ! Swap TOP time levels (= Nrhs_trc, Nbb_trc etc) 99 Krhs = Kbb 100 Kbb = Kmm 101 Kmm = Kaa 102 Kaa = Krhs 103 ! 104 IF( ln_trcrad ) CALL trc_rad( kt, Kbb, Kmm, tr ) ! Correct artificial negative concentrations 91 105 ! 92 106 END IF -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcini.F90
r11027 r11480 40 40 CONTAINS 41 41 42 SUBROUTINE trc_init ( Kbb, Kmm, Kaa )42 SUBROUTINE trc_init 43 43 !!--------------------------------------------------------------------- 44 44 !! *** ROUTINE trc_init *** … … 52 52 !! or read data or analytical formulation 53 53 !!--------------------------------------------------------------------- 54 INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! time level indices55 54 ! 56 55 IF( ln_timing ) CALL timing_start('trc_init') … … 59 58 IF(lwp) WRITE(numout,*) 'trc_init : initial set up of the passive tracers' 60 59 IF(lwp) WRITE(numout,*) '~~~~~~~~' 60 ! 61 ! Initialise time level indices 62 Nbb_trc = 1; Nnn_trc = 2; Naa_trc = 3; Nrhs_trc = Naa_trc 61 63 ! 62 64 CALL trc_ini_ctl ! control … … 71 73 IF(lwp) WRITE(numout,*) 72 74 ! 73 CALL trc_ini_sms( Kmm) ! SMS75 CALL trc_ini_sms( Nnn_trc ) ! SMS 74 76 CALL trc_ini_trp ! passive tracers transport 75 77 CALL trc_ice_ini ! Tracers in sea ice … … 79 81 ENDIF 80 82 ! 81 CALL trc_ini_state( Kbb, Kmm, Kaa) ! passive tracers initialisation : from a restart or from clim83 CALL trc_ini_state( Nbb_trc, Nnn_trc, Naa_trc ) ! passive tracers initialisation : from a restart or from clim 82 84 IF( nn_dttrc /= 1 ) & 83 CALL trc_sub_ini( Kmm) ! Initialize variables for substepping passive tracers84 ! 85 CALL trc_ini_inv( Kmm) ! Inventories85 CALL trc_sub_ini( Nnn_trc ) ! Initialize variables for substepping passive tracers 86 ! 87 CALL trc_ini_inv( Nnn_trc ) ! Inventories 86 88 ! 87 89 IF( ln_timing ) CALL timing_stop('trc_init') -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcstp.F90
r11027 r11480 38 38 39 39 !!---------------------------------------------------------------------- 40 !! time level indices 41 !!---------------------------------------------------------------------- 42 INTEGER, PUBLIC :: Nbb_trc, Nnn_trc, Naa_trc, Nrhs_trc !! used by trc_init 43 44 !!---------------------------------------------------------------------- 40 45 !! NEMO/TOP 4.0 , NEMO Consortium (2018) 41 46 !! $Id$ … … 44 49 CONTAINS 45 50 46 SUBROUTINE trc_stp( kt, Kbb , Kmm, Krhs, Kaa)51 SUBROUTINE trc_stp( kt, Kbb_oce, Kmm_oce, Krhs_oce, Kaa_oce ) 47 52 !!------------------------------------------------------------------- 48 53 !! *** ROUTINE trc_stp *** … … 54 59 !!------------------------------------------------------------------- 55 60 INTEGER, INTENT( in ) :: kt ! ocean time-step index 56 INTEGER, INTENT( in ) :: Kbb , Kmm, Krhs, Kaa! time level indices61 INTEGER, INTENT( in ) :: Kbb_oce, Kmm_oce, Krhs_oce, Kaa_oce ! time level indices 57 62 ! 58 63 INTEGER :: jk, jn ! dummy loop indices … … 76 81 IF( .NOT.ln_linssh ) THEN ! update ocean volume due to ssh temporal evolution 77 82 DO jk = 1, jpk 78 cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm ) * tmask(:,:,jk)83 cvol(:,:,jk) = e1e2t(:,:) * e3t(:,:,jk,Kmm_oce) * tmask(:,:,jk) 79 84 END DO 80 85 IF ( ln_ctl .OR. kt == nitrst .OR. ( ln_check_mass .AND. kt == nitend ) & … … 86 91 IF( l_trcdm2dc ) CALL trc_mean_qsr( kt ) 87 92 ! 88 IF( nn_dttrc /= 1 ) CALL trc_sub_stp( kt, Kbb, Kmm, Krhs ) ! averaging physical variables for sub-stepping 93 IF( nn_dttrc == 1 ) THEN 94 IF(lwp) WRITE(numout,*) "Kbb_oce, Kmm_oce, Kaa_oce, Krhs_oce : ",Kbb_oce, Kmm_oce, Kaa_oce, Krhs_oce 95 IF(lwp) WRITE(numout,*) "Nbb_trc, Nnn_trc, Naa_trc, Nrhs_trc : ",Nbb_trc, Nnn_trc, Naa_trc, Nrhs_trc 96 IF(lwp) CALL FLUSH(numout) 97 CALL mppsync() 98 IF( Kmm_oce /= Nnn_trc .OR. Kaa_oce /= Naa_trc .OR. Krhs_oce /= Nrhs_trc ) THEN 99 ! The nn_dttrc == 1 case depends on the OCE and TRC time indices being the same always. 100 ! If this is not the case then something has gone wrong. 101 CALL ctl_stop( 'trc_stp : nn_dttrc = 1 but OCE and TRC time indices are different! Something has gone wrong.' ) 102 ENDIF 103 ELSE 104 CALL trc_sub_stp( kt, Nbb_trc, Nnn_trc, Nrhs_trc ) ! averaging physical variables for sub-stepping 105 ENDIF 89 106 ! 90 107 IF( MOD( kt , nn_dttrc ) == 0 ) THEN ! only every nn_dttrc time step … … 95 112 ENDIF 96 113 ! 97 tr(:,:,:,:, Krhs) = 0.e0114 tr(:,:,:,:,Nrhs_trc) = 0.e0 98 115 ! 99 116 CALL trc_rst_opn ( kt ) ! Open tracer restart file 100 117 IF( lrst_trc ) CALL trc_rst_cal ( kt, 'WRITE' ) ! calendar 101 CALL trc_wri ( kt, Kmm )! output of passive tracers with iom I/O manager102 CALL trc_sms ( kt, Kbb, Kmm, Krhs )! tracers: sinks and sources103 CALL trc_trp ( kt, Kbb, Kmm, Krhs, Kaa) ! transport of passive tracers118 CALL trc_wri ( kt, Nnn_trc ) ! output of passive tracers with iom I/O manager 119 CALL trc_sms ( kt, Nbb_trc, Nnn_trc, Nrhs_trc ) ! tracers: sinks and sources 120 CALL trc_trp ( kt, Nbb_trc, Nnn_trc, Nrhs_trc, Naa_trc ) ! transport of passive tracers 104 121 IF( kt == nittrc000 ) THEN 105 122 CALL iom_close( numrtr ) ! close input tracer restart file 106 123 IF(lwm) CALL FLUSH( numont ) ! flush namelist output 107 124 ENDIF 108 IF( lrst_trc ) CALL trc_rst_wri ( kt, Kbb, Kmm, Krhs) ! write tracer restart file109 IF( lk_trdmxl_trc ) CALL trd_mxl_trc ( kt, Kmm) ! trends: Mixed-layer110 ! 111 IF( nn_dttrc /= 1 ) CALL trc_sub_reset( kt, Kbb, Kmm, Krhs) ! resetting physical variables when sub-stepping125 IF( lrst_trc ) CALL trc_rst_wri ( kt, Nbb_trc, Nnn_trc, Nrhs_trc ) ! write tracer restart file 126 IF( lk_trdmxl_trc ) CALL trd_mxl_trc ( kt, Nnn_trc ) ! trends: Mixed-layer 127 ! 128 IF( nn_dttrc /= 1 ) CALL trc_sub_reset( kt, Nbb_trc, Nnn_trc, Nrhs_trc ) ! resetting physical variables when sub-stepping 112 129 ! 113 130 ENDIF … … 116 133 ztrai = 0._wp ! content of all tracers 117 134 DO jn = 1, jptra 118 ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn, Kmm) * cvol(:,:,:) )135 ztrai = ztrai + glob_sum( 'trcstp', tr(:,:,:,jn,Nnn_trc) * cvol(:,:,:) ) 119 136 END DO 120 137 IF( lwm ) WRITE(numstr,9300) kt, ztrai / areatot -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/TOP/trcsub.F90
r10969 r11480 92 92 !!------------------------------------------------------------------- 93 93 INTEGER, INTENT( in ) :: kt ! ocean time-step index 94 INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs ! oceantime-level index94 INTEGER, INTENT( in ) :: Kbb, Kmm, Krhs ! time-level index 95 95 ! 96 96 INTEGER :: ji, jj, jk ! dummy loop indices -
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/tests/VORTEX/MY_SRC/domvvl.F90
r11412 r11480 13 13 !! dom_vvl_init : define initial vertical scale factors, depths and column thickness 14 14 !! dom_vvl_sf_nxt : Compute next vertical scale factors 15 !! dom_vvl_sf_ swp: Swap vertical scale factors and update the vertical grid15 !! dom_vvl_sf_update : Swap vertical scale factors and update the vertical grid 16 16 !! dom_vvl_interpol : Interpolate vertical scale factors from one grid point to another 17 17 !! dom_vvl_rst : read/write restart file … … 37 37 PUBLIC dom_vvl_init ! called by domain.F90 38 38 PUBLIC dom_vvl_sf_nxt ! called by step.F90 39 PUBLIC dom_vvl_sf_ swp! called by step.F9039 PUBLIC dom_vvl_sf_update ! called by step.F90 40 40 PUBLIC dom_vvl_interpol ! called by dynnxt.F90 41 41 … … 181 181 ! 182 182 ! !== thickness of the water column !! (ocean portion only) 183 ht _n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) ....184 hu _b(:,:) = e3u(:,:,1,Kbb) * umask(:,:,1)185 hu _n(:,:) = e3u(:,:,1,Kmm) * umask(:,:,1)186 hv _b(:,:) = e3v(:,:,1,Kbb) * vmask(:,:,1)187 hv _n(:,:) = e3v(:,:,1,Kmm) * vmask(:,:,1)183 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) !!gm BUG : this should be 1/2 * e3w(k=1) .... 184 hu(:,:,Kbb) = e3u(:,:,1,Kbb) * umask(:,:,1) 185 hu(:,:,Kmm) = e3u(:,:,1,Kmm) * umask(:,:,1) 186 hv(:,:,Kbb) = e3v(:,:,1,Kbb) * vmask(:,:,1) 187 hv(:,:,Kmm) = e3v(:,:,1,Kmm) * vmask(:,:,1) 188 188 DO jk = 2, jpkm1 189 ht _n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)190 hu _b(:,:) = hu_b(:,:) + e3u(:,:,jk,Kbb) * umask(:,:,jk)191 hu _n(:,:) = hu_n(:,:) + e3u(:,:,jk,Kmm) * umask(:,:,jk)192 hv _b(:,:) = hv_b(:,:) + e3v(:,:,jk,Kbb) * vmask(:,:,jk)193 hv _n(:,:) = hv_n(:,:) + e3v(:,:,jk,Kmm) * vmask(:,:,jk)189 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 190 hu(:,:,Kbb) = hu(:,:,Kbb) + e3u(:,:,jk,Kbb) * umask(:,:,jk) 191 hu(:,:,Kmm) = hu(:,:,Kmm) + e3u(:,:,jk,Kmm) * umask(:,:,jk) 192 hv(:,:,Kbb) = hv(:,:,Kbb) + e3v(:,:,jk,Kbb) * vmask(:,:,jk) 193 hv(:,:,Kmm) = hv(:,:,Kmm) + e3v(:,:,jk,Kmm) * vmask(:,:,jk) 194 194 END DO 195 195 ! 196 196 ! !== inverse of water column thickness ==! (u- and v- points) 197 r1_hu _b(:,:) = ssumask(:,:) / ( hu_b(:,:) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF198 r1_hu _n(:,:) = ssumask(:,:) / ( hu_n(:,:) + 1._wp - ssumask(:,:) )199 r1_hv _b(:,:) = ssvmask(:,:) / ( hv_b(:,:) + 1._wp - ssvmask(:,:) )200 r1_hv _n(:,:) = ssvmask(:,:) / ( hv_n(:,:) + 1._wp - ssvmask(:,:) )197 r1_hu(:,:,Kbb) = ssumask(:,:) / ( hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) ! _i mask due to ISF 198 r1_hu(:,:,Kmm) = ssumask(:,:) / ( hu(:,:,Kmm) + 1._wp - ssumask(:,:) ) 199 r1_hv(:,:,Kbb) = ssvmask(:,:) / ( hv(:,:,Kbb) + 1._wp - ssvmask(:,:) ) 200 r1_hv(:,:,Kmm) = ssvmask(:,:) / ( hv(:,:,Kmm) + 1._wp - ssvmask(:,:) ) 201 201 202 202 ! !== z_tilde coordinate case ==! (Restoring frequencies) … … 550 550 ! *********************************** ! 551 551 552 hu _a(:,:) = e3u(:,:,1,Kaa) * umask(:,:,1)553 hv _a(:,:) = e3v(:,:,1,Kaa) * vmask(:,:,1)552 hu(:,:,Kaa) = e3u(:,:,1,Kaa) * umask(:,:,1) 553 hv(:,:,Kaa) = e3v(:,:,1,Kaa) * vmask(:,:,1) 554 554 DO jk = 2, jpkm1 555 hu _a(:,:) = hu_a(:,:) + e3u(:,:,jk,Kaa) * umask(:,:,jk)556 hv _a(:,:) = hv_a(:,:) + e3v(:,:,jk,Kaa) * vmask(:,:,jk)555 hu(:,:,Kaa) = hu(:,:,Kaa) + e3u(:,:,jk,Kaa) * umask(:,:,jk) 556 hv(:,:,Kaa) = hv(:,:,Kaa) + e3v(:,:,jk,Kaa) * vmask(:,:,jk) 557 557 END DO 558 558 ! ! Inverse of the local depth 559 559 !!gm BUG ? don't understand the use of umask_i here ..... 560 r1_hu _a(:,:) = ssumask(:,:) / ( hu_a(:,:) + 1._wp - ssumask(:,:) )561 r1_hv _a(:,:) = ssvmask(:,:) / ( hv_a(:,:) + 1._wp - ssvmask(:,:) )560 r1_hu(:,:,Kaa) = ssumask(:,:) / ( hu(:,:,Kaa) + 1._wp - ssumask(:,:) ) 561 r1_hv(:,:,Kaa) = ssvmask(:,:) / ( hv(:,:,Kaa) + 1._wp - ssvmask(:,:) ) 562 562 ! 563 563 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_nxt') … … 566 566 567 567 568 SUBROUTINE dom_vvl_sf_ swp( kt, Kbb, Kmm, Kaa )569 !!---------------------------------------------------------------------- 570 !! *** ROUTINE dom_vvl_sf_ swp***568 SUBROUTINE dom_vvl_sf_update( kt, Kbb, Kmm, Kaa ) 569 !!---------------------------------------------------------------------- 570 !! *** ROUTINE dom_vvl_sf_update *** 571 571 !! 572 !! ** Purpose : compute time filter and swap of scale factors572 !! ** Purpose : for z tilde case: compute time filter and swap of scale factors 573 573 !! compute all depths and related variables for next time step 574 574 !! write outputs and restart file 575 575 !! 576 !! ** Method : - swap of e3t with trick for volume/tracer conservation 576 !! ** Method : - swap of e3t with trick for volume/tracer conservation (ONLY FOR Z TILDE CASE) 577 577 !! - reconstruct scale factor at other grid points (interpolate) 578 578 !! - recompute depths and water height fields 579 579 !! 580 !! ** Action : - e3t_(b/n), tilde_e3t_(b/n) and e3(u/v)_nready for next time step580 !! ** Action : - tilde_e3t_(b/n) ready for next time step 581 581 !! - Recompute: 582 582 !! e3(u/v)_b … … 599 599 IF( ln_linssh ) RETURN ! No calculation in linear free surface 600 600 ! 601 IF( ln_timing ) CALL timing_start('dom_vvl_sf_ swp')601 IF( ln_timing ) CALL timing_start('dom_vvl_sf_update') 602 602 ! 603 603 IF( kt == nit000 ) THEN 604 604 IF(lwp) WRITE(numout,*) 605 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_ swp : - time filter and swap of scale factors'606 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~ - interpolate scale factors and compute depths for next time step'605 IF(lwp) WRITE(numout,*) 'dom_vvl_sf_update : - interpolate scale factors and compute depths for next time step' 606 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~' 607 607 ENDIF 608 608 ! … … 619 619 tilde_e3t_n(:,:,:) = tilde_e3t_a(:,:,:) 620 620 ENDIF 621 gdept(:,:,:,Kbb) = gdept(:,:,:,Kmm)622 gdepw(:,:,:,Kbb) = gdepw(:,:,:,Kmm)623 624 e3t(:,:,:,Kmm) = e3t(:,:,:,Kaa)625 e3u(:,:,:,Kmm) = e3u(:,:,:,Kaa)626 e3v(:,:,:,Kmm) = e3v(:,:,:,Kaa)627 621 628 622 ! Compute all missing vertical scale factor and depths … … 631 625 ! -------------------------------------- 632 626 ! - ML - e3u(:,:,:,Kbb) and e3v(:,:,:,Kbb) are allready computed in dynnxt 633 ! - JC - hu _b, hv_b, hur_b, hvr_b also627 ! - JC - hu(:,:,:,Kbb), hv(:,:,:,:,Kbb), hur_b, hvr_b also 634 628 635 629 CALL dom_vvl_interpol( e3u(:,:,:,Kmm), e3f(:,:,:), 'F' ) … … 663 657 ! Local depth and Inverse of the local depth of the water 664 658 ! ------------------------------------------------------- 665 hu_n(:,:) = hu_a(:,:) ; r1_hu_n(:,:) = r1_hu_a(:,:) 666 hv_n(:,:) = hv_a(:,:) ; r1_hv_n(:,:) = r1_hv_a(:,:) 667 ! 668 ht_n(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 659 ! 660 ht(:,:) = e3t(:,:,1,Kmm) * tmask(:,:,1) 669 661 DO jk = 2, jpkm1 670 ht _n(:,:) = ht_n(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk)662 ht(:,:) = ht(:,:) + e3t(:,:,jk,Kmm) * tmask(:,:,jk) 671 663 END DO 672 664 … … 675 667 IF( lrst_oce ) CALL dom_vvl_rst( kt, Kbb, Kmm, 'WRITE' ) 676 668 ! 677 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_ swp')678 ! 679 END SUBROUTINE dom_vvl_sf_ swp669 IF( ln_timing ) CALL timing_stop('dom_vvl_sf_update') 670 ! 671 END SUBROUTINE dom_vvl_sf_update 680 672 681 673 … … 931 923 ELSE 932 924 ! 933 ! ! 934 ! usr_def_istate called here only to get ssh, that is needed to initialize e3t_b and e3t_n 935 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 936 ! usr_def_istate will be called again in istate_init to initialize ts(bn), ssh(bn), u(bn) and v(bn) 925 ! usr_def_istate called here only to get ssh(Kbb) needed to initialize e3t(Kbb) and e3t(Kmm) 926 ! 927 CALL usr_def_istate( gdept_0, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb), ssh(:,:,Kbb) ) 928 ! 929 ! usr_def_istate will be called again in istate_init to initialize ts, ssh, u and v 937 930 ! 938 931 DO jk=1,jpk 939 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) 940 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) &941 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(Kbb) != 0 on land points932 e3t(:,:,jk,Kbb) = e3t_0(:,:,jk) * ( ht_0(:,:) + ssh(:,:,Kbb) ) & 933 & / ( ht_0(:,:) + 1._wp - ssmask(:,:) ) * tmask(:,:,jk) & 934 & + e3t_0(:,:,jk) * ( 1._wp - tmask(:,:,jk) ) ! make sure e3t(:,:,:,Kbb) != 0 on land points 942 935 END DO 943 936 e3t(:,:,:,Kmm) = e3t(:,:,:,Kbb) 944 ssh(:,: ,Kmm) = ssh(:,: ,Kbb)! needed later for gde3w937 ssh(:,:,Kmm) = ssh(:,:,Kbb) ! needed later for gde3w 945 938 ! 946 939 END IF ! end of ll_wd edits
Note: See TracChangeset
for help on using the changeset viewer.