Changeset 10978 for NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
- Timestamp:
- 2019-05-15T09:41:30+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/DYN/sshwzv.F90
r10969 r10978 54 54 CONTAINS 55 55 56 SUBROUTINE ssh_nxt( kt, Kbb, Kmm )56 SUBROUTINE ssh_nxt( kt, Kbb, Kmm, pssh, Kaa ) 57 57 !!---------------------------------------------------------------------- 58 58 !! *** ROUTINE ssh_nxt *** 59 59 !! 60 !! ** Purpose : compute the after ssh (ssh a)60 !! ** Purpose : compute the after ssh (ssh(Kaa)) 61 61 !! 62 62 !! ** Method : - Using the incompressibility hypothesis, the ssh increment … … 64 64 !! by the time step. 65 65 !! 66 !! ** action : ssh a, after sea surface height66 !! ** action : ssh(:,:,Kaa), after sea surface height 67 67 !! 68 68 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 69 69 !!---------------------------------------------------------------------- 70 INTEGER, INTENT(in) :: kt ! time step 71 INTEGER, INTENT(in) :: Kbb, Kmm ! time level index 70 INTEGER , INTENT(in ) :: kt ! time step 71 INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index 72 REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height 72 73 ! 73 74 INTEGER :: jk ! dummy loop indice … … 92 93 ! !------------------------------! 93 94 IF(ln_wd_il) THEN 94 CALL wad_lmt( sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt)95 CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 95 96 ENDIF 96 97 … … 99 100 zhdiv(:,:) = 0._wp 100 101 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 101 zhdiv(:,:) = zhdiv(:,:) + e3t _n(:,:,jk) * hdivn(:,:,jk)102 zhdiv(:,:) = zhdiv(:,:) + e3t(:,:,jk,Kmm) * hdiv(:,:,jk) 102 103 END DO 103 104 ! ! Sea surface elevation time stepping … … 105 106 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 106 107 ! 107 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:)108 pssh(:,:,Kaa) = ( pssh(:,:,Kbb) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 108 109 ! 109 110 #if defined key_agrif … … 113 114 IF ( .NOT.ln_dynspg_ts ) THEN 114 115 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 boundaries116 CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1. ) ! Not sure that's necessary 117 CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries 117 118 ENDIF 118 119 ENDIF … … 121 122 ! !------------------------------! 122 123 ! 123 IF(ln_ctl) CALL prt_ctl( tab2d_1= ssha, clinfo1=' ssha- : ', mask1=tmask )124 IF(ln_ctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask ) 124 125 ! 125 126 IF( ln_timing ) CALL timing_stop('ssh_nxt') … … 128 129 129 130 130 SUBROUTINE wzv( kt )131 SUBROUTINE wzv( kt, Kbb, Kmm, pww, Kaa ) 131 132 !!---------------------------------------------------------------------- 132 133 !! *** ROUTINE wzv *** … … 139 140 !! The boundary conditions are w=0 at the bottom (no flux) and. 140 141 !! 141 !! ** action : wn: now vertical velocity142 !! ** action : pww : now vertical velocity 142 143 !! 143 144 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 144 145 !!---------------------------------------------------------------------- 145 INTEGER, INTENT(in) :: kt ! time step 146 INTEGER , INTENT(in) :: kt ! time step 147 INTEGER , INTENT(in) :: Kbb, Kmm, Kaa ! time level indices 148 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! now vertical velocity 146 149 ! 147 150 INTEGER :: ji, jj, jk ! dummy loop indices … … 157 160 IF(lwp) WRITE(numout,*) '~~~~~ ' 158 161 ! 159 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)162 pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 160 163 ENDIF 161 164 ! !------------------------------! … … 179 182 CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.) ! - ML - Perhaps not necessary: not used for horizontal "connexions" 180 183 ! ! Is it problematic to have a wrong vertical velocity in boundary cells? 181 ! ! Same question holds for hdiv n. Perhaps just for security184 ! ! Same question holds for hdiv. Perhaps just for security 182 185 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 183 186 ! 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.e0187 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) + zhdiv(:,:,jk) & 188 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 189 END DO 190 ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0 188 191 DEALLOCATE( zhdiv ) 189 192 ELSE ! z_star and linear free surface cases 190 193 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 191 194 ! 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)195 pww(:,:,jk) = pww(:,:,jk+1) - ( e3t(:,:,jk,Kmm) * hdiv(:,:,jk) & 196 & + z1_2dt * ( e3t(:,:,jk,Kaa) - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk) 194 197 END DO 195 198 ENDIF … … 197 200 IF( ln_bdy ) THEN 198 201 DO jk = 1, jpkm1 199 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)202 pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:) 200 203 END DO 201 204 ENDIF … … 203 206 #if defined key_agrif 204 207 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 ! south208 IF ((nbondi == 1).OR.(nbondi == 2)) pww(nlci-1 , : ,:) = 0.e0 ! east 209 IF ((nbondi == -1).OR.(nbondi == 2)) pww(2 , : ,:) = 0.e0 ! west 210 IF ((nbondj == 1).OR.(nbondj == 2)) pww(: ,nlcj-1 ,:) = 0.e0 ! north 211 IF ((nbondj == -1).OR.(nbondj == 2)) pww(: ,2 ,:) = 0.e0 ! south 209 212 ENDIF 210 213 #endif … … 215 218 216 219 217 SUBROUTINE ssh_swp( kt )220 SUBROUTINE ssh_swp( kt, Kbb, Kmm, Kaa ) 218 221 !!---------------------------------------------------------------------- 219 222 !! *** ROUTINE ssh_nxt *** … … 221 224 !! ** Purpose : achieve the sea surface height time stepping by 222 225 !! applying Asselin time filter and swapping the arrays 223 !! ssh aalready computed in ssh_nxt226 !! ssh(:,:,Kaa) already computed in ssh_nxt 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 !! ssh n = ssha + atfp * ( sshb -2 sshn + ssha)230 !! ssh(:,:,Kmm) = ssh(:,:,Kaa) + atfp * ( ssh(:,:,Kbb) -2 ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 228 231 !! - atfp * rdt * ( emp_b - emp ) / rau0 229 !! ssh n = ssha230 !! 231 !! ** action : - ssh b, sshn: before & now sea surface height232 !! ssh(:,:,Kmm) = ssh(:,:,Kaa) 233 !! 234 !! ** action : - ssh(:,:,Kbb), ssh(:,:,Kmm) : before & now sea surface height 232 235 !! ready for the next time step 233 236 !! 234 237 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 235 238 !!---------------------------------------------------------------------- 236 INTEGER, INTENT(in) :: kt ! ocean time-step index 239 INTEGER, INTENT(in) :: kt ! ocean time-step index 240 INTEGER, INTENT(in) :: Kbb, Kmm, Kaa ! ocean time-step index 237 241 ! 238 242 REAL(wp) :: zcoef ! local scalar … … 248 252 ! !== Euler time-stepping: no filter, just swap ==! 249 253 IF ( neuler == 0 .AND. kt == nit000 ) THEN 250 ssh n(:,:) = ssha(:,:) ! now <-- after (before already = now)254 ssh(:,:,Kmm) = ssh(:,:,Kaa) ! now <-- after (before already = now) 251 255 ! 252 256 ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==! 253 257 ! ! before <-- now filtered 254 ssh b(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) )258 ssh(:,:,Kbb) = ssh(:,:,Kmm) + atfp * ( ssh(:,:,Kbb) - 2 * ssh(:,:,Kmm) + ssh(:,:,Kaa) ) 255 259 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 256 260 zcoef = atfp * rdt * r1_rau0 257 ssh b(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) &261 ssh(:,:,Kbb) = ssh(:,:,Kbb) - zcoef * ( emp_b(:,:) - emp (:,:) & 258 262 & - rnf_b(:,:) + rnf (:,:) & 259 263 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 260 264 ENDIF 261 ssh n(:,:) = ssha(:,:) ! now <-- after262 ENDIF 263 ! 264 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssh b, clinfo1=' sshb- : ', mask1=tmask )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 ) 265 269 ! 266 270 IF( ln_timing ) CALL timing_stop('ssh_swp') … … 268 272 END SUBROUTINE ssh_swp 269 273 270 SUBROUTINE wAimp( kt )274 SUBROUTINE wAimp( kt, Kmm ) 271 275 !!---------------------------------------------------------------------- 272 276 !! *** ROUTINE wAimp *** … … 277 281 !! ** Method : - 278 282 !! 279 !! ** action : w n: now vertical velocity (to be handled explicitly)283 !! ** action : ww : now vertical velocity (to be handled explicitly) 280 284 !! : wi : now vertical velocity (for implicit treatment) 281 285 !! … … 283 287 !!---------------------------------------------------------------------- 284 288 INTEGER, INTENT(in) :: kt ! time step 289 INTEGER, INTENT(in) :: Kmm ! time level index 285 290 ! 286 291 INTEGER :: ji, jj, jk ! dummy loop indices … … 305 310 DO jj = 2, jpjm1 306 311 DO ji = 2, fs_jpim1 ! vector opt. 307 z1_e3w = 1._wp / e3w _n(ji,jj,jk)308 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( w n(ji,jj,jk) , 0._wp ) - MIN( wn(ji,jj,jk+1) , 0._wp ) ) &309 & + ( MAX( e2u(ji ,jj)*e3uw _n(ji ,jj,jk)*un(ji ,jj,jk), 0._wp ) - &310 & MIN( e2u(ji-1,jj)*e3uw _n(ji-1,jj,jk)*un(ji-1,jj,jk), 0._wp ) ) &312 z1_e3w = 1._wp / e3w(ji,jj,jk,Kmm) 313 Cu_adv(ji,jj,jk) = r2dt * ( ( MAX( ww(ji,jj,jk) , 0._wp ) - MIN( ww(ji,jj,jk+1) , 0._wp ) ) & 314 & + ( MAX( e2u(ji ,jj)*e3uw(ji ,jj,jk,Kmm)*uu(ji ,jj,jk,Kmm), 0._wp ) - & 315 & MIN( e2u(ji-1,jj)*e3uw(ji-1,jj,jk,Kmm)*uu(ji-1,jj,jk,Kmm), 0._wp ) ) & 311 316 & * r1_e1e2t(ji,jj) & 312 & + ( MAX( e1v(ji,jj )*e3vw _n(ji,jj ,jk)*vn(ji,jj ,jk), 0._wp ) - &313 & MIN( e1v(ji,jj-1)*e3vw _n(ji,jj-1,jk)*vn(ji,jj-1,jk), 0._wp ) ) &317 & + ( MAX( e1v(ji,jj )*e3vw(ji,jj ,jk,Kmm)*vv(ji,jj ,jk,Kmm), 0._wp ) - & 318 & MIN( e1v(ji,jj-1)*e3vw(ji,jj-1,jk,Kmm)*vv(ji,jj-1,jk,Kmm), 0._wp ) ) & 314 319 & * r1_e1e2t(ji,jj) & 315 320 & ) * z1_e3w … … 338 343 zcff = MIN(1._wp, zcff) 339 344 ! 340 wi(ji,jj,jk) = zcff * w n(ji,jj,jk)341 w n(ji,jj,jk) = ( 1._wp - zcff ) * wn(ji,jj,jk)345 wi(ji,jj,jk) = zcff * ww(ji,jj,jk) 346 ww(ji,jj,jk) = ( 1._wp - zcff ) * ww(ji,jj,jk) 342 347 ! 343 348 Cu_adv(ji,jj,jk) = zcff ! Reuse array to output coefficient … … 351 356 CALL iom_put("wimp",wi) 352 357 CALL iom_put("wi_cff",Cu_adv) 353 CALL iom_put("wexp",w n)358 CALL iom_put("wexp",ww) 354 359 ! 355 360 IF( ln_timing ) CALL timing_stop('wAimp')
Note: See TracChangeset
for help on using the changeset viewer.
