- Timestamp:
- 2019-11-22T15:29:17+01:00 (5 years ago)
- Location:
- NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
OCE/DYN/sshwzv.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src
- Property svn:mergeinfo deleted
-
NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/DYN/sshwzv.F90
r11414 r11949 10 10 !! 3.3 ! 2011-10 (M. Leclair) split former ssh_wzv routine and remove all vvl related work 11 11 !! 4.0 ! 2018-12 (A. Coward) add mixed implicit/explicit advection 12 !! 4.1 ! 2019-08 (A. Coward, D. Storkey) Rename ssh_nxt -> ssh_atf. Now only does time filtering. 12 13 !!---------------------------------------------------------------------- 13 14 14 15 !!---------------------------------------------------------------------- 15 16 !! ssh_nxt : after ssh 16 !! ssh_ swp : filter ans swapthe ssh arrays17 !! ssh_atf : time filter the ssh arrays 17 18 !! wzv : compute now vertical velocity 18 19 !!---------------------------------------------------------------------- … … 44 45 PUBLIC wzv ! called by step.F90 45 46 PUBLIC wAimp ! called by step.F90 46 PUBLIC ssh_ swp! called by step.F9047 PUBLIC ssh_atf ! called by step.F90 47 48 48 49 !! * Substitutions … … 55 56 CONTAINS 56 57 57 SUBROUTINE ssh_nxt( kt )58 SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa ) 58 59 !!---------------------------------------------------------------------- 59 60 !! *** ROUTINE ssh_nxt *** 60 61 !! 61 !! ** Purpose : compute the after ssh (ssh a)62 !! ** Purpose : compute the after ssh (ssh(Kaa)) 62 63 !! 63 64 !! ** Method : - Using the incompressibility hypothesis, the ssh increment … … 65 66 !! by the time step. 66 67 !! 67 !! ** action : ssh a, after sea surface height68 !! ** action : ssh(:,:,Kaa), after sea surface height 68 69 !! 69 70 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 70 71 !!---------------------------------------------------------------------- 71 INTEGER, INTENT(in) :: kt ! time step 72 INTEGER , INTENT(in ) :: kt ! time step 73 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index 74 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 72 75 ! 73 76 INTEGER :: jk ! dummy loop indice … … 92 95 ! !------------------------------! 93 96 IF(ln_wd_il) THEN 94 CALL wad_lmt( sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)95 ENDIF 96 97 CALL div_hor( kt )! Horizontal divergence97 CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt, Kmm, uu, vv ) 98 ENDIF 99 100 CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence 98 101 ! 99 102 zhdiv(:,:) = 0._wp 100 103 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 101 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)104 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 102 105 END DO 103 106 ! ! Sea surface elevation time stepping … … 105 108 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 106 109 ! 107 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)110 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 108 111 ! 109 112 #if defined key_agrif 110 CALL agrif_ssh( kt )113 Kbb_a = Kbb; Kmm_a = Kmm; Krhs_a = Kaa; CALL agrif_ssh( kt ) 111 114 #endif 112 115 ! 113 116 IF ( .NOT.ln_dynspg_ts ) THEN 114 117 IF( ln_bdy ) THEN 115 CALL lbc_lnk( 'sshwzv', ssha, 'T', 1. ) ! Not sure that's necessary116 CALL bdy_ssh( ssha) ! Duplicate sea level across open boundaries118 CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. ) ! Not sure that's necessary 119 CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries 117 120 ENDIF 118 121 ENDIF … … 121 124 ! !------------------------------! 122 125 ! 123 IF(ln_ctl) CALL prt_ctl( tab2d_1= ssha, clinfo1=' ssha- : ', mask1=tmask )126 IF(ln_ctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask ) 124 127 ! 125 128 IF( ln_timing ) CALL timing_stop('ssh_nxt') … … 128 131 129 132 130 SUBROUTINE wzv( kt )133 SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa ) 131 134 !!---------------------------------------------------------------------- 132 135 !! *** ROUTINE wzv *** … … 139 142 !! The boundary conditions are w=0 at the bottom (no flux) and. 140 143 !! 141 !! ** action : wn: now vertical velocity144 !! ** action : pww : now vertical velocity 142 145 !! 143 146 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 144 147 !!---------------------------------------------------------------------- 145 INTEGER, INTENT(in) :: kt ! time step 148 INTEGER , INTENT(in) :: kt ! time step 149 INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices 150 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! now vertical velocity 146 151 ! 147 152 INTEGER :: ji, jj, jk ! dummy loop indices … … 157 162 IF(lwp) WRITE(numout,*) '~~~~~ ' 158 163 ! 159 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)164 pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 160 165 ENDIF 161 166 ! !------------------------------! … … 179 184 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 180 185 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 181 ! ! Same question holds for hdiv n. Perhaps just for security186 ! ! Same question holds for hdiv. Perhaps just for security 182 187 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 183 188 ! computation of w 184 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) &185 & + z1_2dt * ( e3t _a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)186 END DO 187 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0189 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) & 190 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 191 END DO 192 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 188 193 DEALLOCATE( zhdiv ) 189 194 ELSE ! z_star and linear free surface cases 190 195 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 196 ! computation of w 192 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) &193 & + z1_2dt * ( e3t _a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk)197 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 198 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 194 199 END DO 195 200 ENDIF … … 197 202 IF( ln_bdy ) THEN 198 203 DO jk = 1, jpkm1 199 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)204 pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 200 205 END DO 201 206 ENDIF … … 203 208 #if defined key_agrif 204 209 IF( .NOT. AGRIF_Root() ) THEN 205 IF ((nbondi == 1).OR.(nbondi == 2)) wn(nlci-1 , : ,:) = 0.e0 ! east206 IF ((nbondi == -1).OR.(nbondi == 2)) wn(2 , : ,:) = 0.e0 ! west207 IF ((nbondj == 1).OR.(nbondj == 2)) wn(: ,nlcj-1 ,:) = 0.e0 ! north208 IF ((nbondj == -1).OR.(nbondj == 2)) wn(: ,2 ,:) = 0.e0 ! south210 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 211 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 212 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 213 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 209 214 ENDIF 210 215 #endif … … 215 220 216 221 217 SUBROUTINE ssh_swp( kt ) 218 !!---------------------------------------------------------------------- 219 !! *** ROUTINE ssh_nxt *** 220 !! 221 !! ** Purpose : achieve the sea surface height time stepping by 222 !! applying Asselin time filter and swapping the arrays 223 !! ssha already computed in ssh_nxt 222 SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh ) 223 !!---------------------------------------------------------------------- 224 !! *** ROUTINE ssh_atf *** 225 !! 226 !! ** Purpose : Apply Asselin time filter to now SSH. 224 227 !! 225 228 !! ** Method : - apply Asselin time fiter to now ssh (excluding the forcing 226 229 !! from the filter, see Leclair and Madec 2010) and swap : 227 !! sshn = ssha + atfp * ( sshb -2 sshn + ssha)230 !! pssh(:,:,Kmm) = pssh(:,:,Kaa) + atfp * ( pssh(:,:,Kbb) -2 pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 228 231 !! - atfp * rdt * ( emp_b - emp ) / rau0 229 !! sshn = ssha 230 !! 231 !! ** action : - sshb, sshn : before & now sea surface height 232 !! ready for the next time step 232 !! 233 !! ** action : - pssh(:,:,Kmm) time filtered 233 234 !! 234 235 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 235 236 !!---------------------------------------------------------------------- 236 INTEGER, INTENT(in) :: kt ! ocean time-step index 237 INTEGER , INTENT(in ) :: kt ! ocean time-step index 238 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! ocean time level indices 239 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! SSH field 237 240 ! 238 241 REAL(wp) :: zcoef ! local scalar 239 242 !!---------------------------------------------------------------------- 240 243 ! 241 IF( ln_timing ) CALL timing_start('ssh_ swp')244 IF( ln_timing ) CALL timing_start('ssh_atf') 242 245 ! 243 246 IF( kt == nit000 ) THEN 244 247 IF(lwp) WRITE(numout,*) 245 IF(lwp) WRITE(numout,*) 'ssh_ swp : Asselin time filter and swapof sea surface height'248 IF(lwp) WRITE(numout,*) 'ssh_atf : Asselin time filter of sea surface height' 246 249 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 247 250 ENDIF 248 251 ! !== Euler time-stepping: no filter, just swap ==! 249 IF ( neuler == 0 .AND. kt == nit000 ) THEN 250 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 251 ! 252 ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==! 253 ! ! before <-- now filtered 254 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 255 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 252 IF ( .NOT.( neuler == 0 .AND. kt == nit000 ) ) THEN ! Only do time filtering for leapfrog timesteps 253 ! ! filtered "now" field 254 pssh(:,:,Kmm) = pssh(:,:,Kmm) + atfp * ( pssh(:,:,Kbb) - 2 * pssh(:,:,Kmm) + pssh(:,:,Kaa) ) 255 IF( .NOT.ln_linssh ) THEN ! "now" <-- with forcing removed 256 256 zcoef = atfp * rdt * r1_rau0 257 sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) &257 pssh(:,:,Kmm) = pssh(:,:,Kmm) - zcoef * ( emp_b(:,:) - emp (:,:) & 258 258 & - rnf_b(:,:) + rnf (:,:) & 259 259 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 260 260 ENDIF 261 sshn(:,:) = ssha(:,:) ! now <-- after 262 ENDIF 263 ! 264 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask ) 265 ! 266 IF( ln_timing ) CALL timing_stop('ssh_swp') 267 ! 268 END SUBROUTINE ssh_swp 269 270 SUBROUTINE wAimp( kt ) 261 ENDIF 262 ! 263 IF(ln_ctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kmm), clinfo1=' pssh(:,:,Kmm) - : ', mask1=tmask ) 264 ! 265 IF( ln_timing ) CALL timing_stop('ssh_atf') 266 ! 267 END SUBROUTINE ssh_atf 268 269 SUBROUTINE wAimp( kt, Kmm ) 271 270 !!---------------------------------------------------------------------- 272 271 !! *** ROUTINE wAimp *** … … 277 276 !! ** Method : - 278 277 !! 279 !! ** action : w n: now vertical velocity (to be handled explicitly)278 !! ** action : ww : now vertical velocity (to be handled explicitly) 280 279 !! : wi : now vertical velocity (for implicit treatment) 281 280 !! … … 285 284 !!---------------------------------------------------------------------- 286 285 INTEGER, INTENT(in) :: kt ! time step 286 INTEGER, INTENT(in) :: Kmm ! time level index 287 287 ! 288 288 INTEGER :: ji, jj, jk ! dummy loop indices … … 308 308 DO jj = 2, jpjm1 309 309 DO ji = 2, fs_jpim1 ! vector opt. 310 z1_e3t = 1._wp / e3t _n(ji,jj,jk)310 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 311 311 ! 2*rdt and not r2dt (for restartability) 312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( w n(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &313 & + ( MAX( e2u(ji ,jj)*e3u _n(ji ,jj,jk)*un(ji ,jj,jk) + un_td(ji ,jj,jk), 0._wp ) - &314 & MIN( e2u(ji-1,jj)*e3u _n(ji-1,jj,jk)*un(ji-1,jj,jk) + un_td(ji-1,jj,jk), 0._wp ) ) &312 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 313 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm) + un_td(ji ,jj,jk), 0._wp ) - & 314 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm) + un_td(ji-1,jj,jk), 0._wp ) ) & 315 315 & * r1_e1e2t(ji,jj) & 316 & + ( MAX( e1v(ji,jj )*e3v _n(ji,jj ,jk)*vn(ji,jj ,jk) + vn_td(ji,jj ,jk), 0._wp ) - &317 & MIN( e1v(ji,jj-1)*e3v _n(ji,jj-1,jk)*vn(ji,jj-1,jk) + vn_td(ji,jj-1,jk), 0._wp ) ) &316 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm) + vn_td(ji,jj ,jk), 0._wp ) - & 317 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm) + vn_td(ji,jj-1,jk), 0._wp ) ) & 318 318 & * r1_e1e2t(ji,jj) & 319 319 & ) * z1_e3t … … 325 325 DO jj = 2, jpjm1 326 326 DO ji = 2, fs_jpim1 ! vector opt. 327 z1_e3t = 1._wp / e3t _n(ji,jj,jk)327 z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm) 328 328 ! 2*rdt and not r2dt (for restartability) 329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( w n(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &330 & + ( MAX( e2u(ji ,jj)*e3u _n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - &331 & MIN( e2u(ji-1,jj)*e3u _n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) &329 Cu_adv(ji,jj,jk) = 2._wp * rdt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 330 & + ( MAX( e2u(ji ,jj)*e3u(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 331 & MIN( e2u(ji-1,jj)*e3u(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 332 332 & * r1_e1e2t(ji,jj) & 333 & + ( MAX( e1v(ji,jj )*e3v _n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - &334 & MIN( e1v(ji,jj-1)*e3v _n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) &333 & + ( MAX( e1v(ji,jj )*e3v(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 334 & MIN( e1v(ji,jj-1)*e3v(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 335 335 & * r1_e1e2t(ji,jj) & 336 336 & ) * z1_e3t … … 366 366 zcff = MIN(1._wp, zcff) 367 367 ! 368 wi(ji,jj,jk) = zcff * w n(ji,jj,jk)369 w n(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)368 wi(ji,jj,jk) = zcff * ww(ji,jj,jk) 369 ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 370 370 ! 371 371 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient below and in stp_ctl … … 381 381 CALL iom_put("wimp",wi) 382 382 CALL iom_put("wi_cff",Cu_adv) 383 CALL iom_put("wexp",w n)383 CALL iom_put("wexp",ww) 384 384 ! 385 385 IF( ln_timing ) CALL timing_stop('wAimp')
Note: See TracChangeset
for help on using the changeset viewer.
