- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r5467 r6808 12 12 13 13 !!---------------------------------------------------------------------- 14 !! ssh_nxt : after ssh 15 !! ssh_swp : filter ans swap the ssh arrays 16 !! wzv : compute now vertical velocity 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and tracers variables 19 USE dom_oce ! ocean space and time domain variables 20 USE sbc_oce ! surface boundary condition: ocean 21 USE domvvl ! Variable volume 22 USE divcur ! hor. divergence and curl (div & cur routines) 23 USE restart ! only for lrst_oce 24 USE in_out_manager ! I/O manager 25 USE prtctl ! Print control 26 USE phycst 27 USE lbclnk ! ocean lateral boundary condition (or mpp link) 28 USE lib_mpp ! MPP library 29 USE bdy_oce 30 USE bdy_par 31 USE bdydyn2d ! bdy_ssh routine 14 !! ssh_nxt : after ssh 15 !! ssh_swp : filter ans swap the ssh arrays 16 !! wzv : compute now vertical velocity 17 !!---------------------------------------------------------------------- 18 USE oce ! ocean dynamics and tracers variables 19 USE dom_oce ! ocean space and time domain variables 20 USE sbc_oce ! surface boundary condition: ocean 21 USE domvvl ! Variable volume 22 USE divhor ! horizontal divergence 23 USE phycst ! physical constants 24 USE bdy_oce ! 25 USE bdy_par ! 26 USE bdydyn2d ! bdy_ssh routine 32 27 #if defined key_agrif 33 USE agrif_opa_update34 28 USE agrif_opa_interp 35 29 #endif 36 30 #if defined key_asminc 37 USE asminc ! Assimilation increment 38 #endif 39 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 31 USE asminc ! Assimilation increment 32 #endif 33 ! 34 USE in_out_manager ! I/O manager 35 USE restart ! only for lrst_oce 36 USE prtctl ! Print control 37 USE lbclnk ! ocean lateral boundary condition (or mpp link) 38 USE lib_mpp ! MPP library 39 USE wrk_nemo ! Memory Allocation 40 USE timing ! Timing 41 USE wet_dry ! Wetting/Drying flux limting 41 42 42 43 IMPLICIT NONE … … 48 49 49 50 !! * Substitutions 50 # include "domzgr_substitute.h90"51 51 # include "vectopt_loop_substitute.h90" 52 52 !!---------------------------------------------------------------------- … … 67 67 !! by the time step. 68 68 !! 69 !! ** action : ssha :after sea surface height69 !! ** action : ssha, after sea surface height 70 70 !! 71 71 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 72 72 !!---------------------------------------------------------------------- 73 ! 74 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 75 INTEGER, INTENT(in) :: kt ! time step 73 INTEGER, INTENT(in) :: kt ! time step 76 74 ! 77 INTEGER :: jk ! dummy loop indice 78 REAL(wp) :: z2dt, z1_rau0 ! local scalars 79 !!---------------------------------------------------------------------- 80 ! 81 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 82 ! 83 CALL wrk_alloc( jpi, jpj, zhdiv ) 75 INTEGER :: jk ! dummy loop indice 76 REAL(wp) :: z2dt, zcoef ! local scalars 77 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv ! 2D workspace 78 !!---------------------------------------------------------------------- 79 ! 80 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 81 ! 82 CALL wrk_alloc( jpi,jpj, zhdiv ) 84 83 ! 85 84 IF( kt == nit000 ) THEN 86 !87 85 IF(lwp) WRITE(numout,*) 88 86 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 89 87 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 90 ! 91 ENDIF 92 ! 93 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 88 ENDIF 89 ! 90 CALL div_hor( kt ) ! Horizontal divergence 94 91 ! 95 92 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) … … 101 98 zhdiv(:,:) = 0._wp 102 99 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 103 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)100 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 104 101 END DO 105 102 ! ! Sea surface elevation time stepping … … 107 104 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 108 105 ! 109 z1_rau0 = 0.5_wp * r1_rau0 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 112 #if ! defined key_dynspg_ts 113 ! These lines are not necessary with time splitting since 114 ! boundary condition on sea level is set during ts loop 115 #if defined key_agrif 116 CALL agrif_ssh( kt ) 117 #endif 118 #if defined key_bdy 119 IF (lk_bdy) THEN 120 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 121 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 122 ENDIF 123 #endif 124 #endif 106 zcoef = 0.5_wp * r1_rau0 107 108 IF(ln_wd) CALL wad_lmt(sshb, zcoef * (emp_b(:,:) + emp(:,:)), z2dt) 109 110 ssha(:,:) = ( sshb(:,:) - z2dt * ( zcoef * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * ssmask(:,:) 111 112 IF ( .NOT.ln_dynspg_ts ) THEN 113 ! These lines are not necessary with time splitting since 114 ! boundary condition on sea level is set during ts loop 115 # if defined key_agrif 116 CALL agrif_ssh( kt ) 117 # endif 118 # if defined key_bdy 119 IF( lk_bdy ) THEN 120 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 121 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 122 ENDIF 123 # endif 124 ENDIF 125 125 126 126 #if defined key_asminc 127 ! ! Include the IAU weighted SSH increment 128 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 127 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN ! Include the IAU weighted SSH increment 129 128 CALL ssh_asm_inc( kt ) 130 129 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 131 130 ENDIF 132 131 #endif 133 134 132 ! !------------------------------! 135 133 ! ! outputs ! … … 160 158 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 161 159 !!---------------------------------------------------------------------- 162 ! 163 INTEGER, INTENT(in) :: kt ! time step 160 INTEGER, INTENT(in) :: kt ! time step 161 ! 162 INTEGER :: ji, jj, jk ! dummy loop indices 163 REAL(wp) :: z1_2dt ! local scalars 164 164 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 165 165 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 166 ! 167 INTEGER :: ji, jj, jk ! dummy loop indices 168 REAL(wp) :: z1_2dt ! local scalars 169 !!---------------------------------------------------------------------- 170 171 IF( nn_timing == 1 ) CALL timing_start('wzv') 166 !!---------------------------------------------------------------------- 167 ! 168 IF( nn_timing == 1 ) CALL timing_start('wzv') 172 169 ! 173 170 IF( kt == nit000 ) THEN 174 !175 171 IF(lwp) WRITE(numout,*) 176 172 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 178 174 ! 179 175 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 180 !181 176 ENDIF 182 177 ! !------------------------------! … … 194 189 DO jj = 2, jpjm1 195 190 DO ji = fs_2, fs_jpim1 ! vector opt. 196 zhdiv(ji,jj,jk) = r1_e1 2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )191 zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) ) 197 192 END DO 198 193 END DO … … 203 198 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 204 199 ! computation of w 205 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)&206 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )) * tmask(:,:,jk)200 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk) & 201 & + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 207 202 END DO 208 203 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 … … 211 206 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 212 207 ! computation of w 213 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk)&214 & + z1_2dt * ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )) * tmask(:,:,jk)208 wn(:,:,jk) = wn(:,:,jk+1) - ( e3t_n(:,:,jk) * hdivn(:,:,jk) & 209 & + z1_2dt * ( e3t_a(:,:,jk) - e3t_b(:,:,jk) ) ) * tmask(:,:,jk) 215 210 END DO 216 211 ENDIF 217 212 218 213 #if defined key_bdy 219 IF (lk_bdy) THEN214 IF( lk_bdy ) THEN 220 215 DO jk = 1, jpkm1 221 216 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 225 220 ! 226 221 IF( nn_timing == 1 ) CALL timing_stop('wzv') 227 228 222 ! 229 223 END SUBROUTINE wzv 224 230 225 231 226 SUBROUTINE ssh_swp( kt ) … … 249 244 !!---------------------------------------------------------------------- 250 245 INTEGER, INTENT(in) :: kt ! ocean time-step index 246 ! 247 REAL(wp) :: zcoef ! local scalar 251 248 !!---------------------------------------------------------------------- 252 249 ! … … 258 255 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 259 256 ENDIF 260 261 # if defined key_dynspg_ts 262 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 263 # else 264 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 265 #endif 266 sshb(:,:) = sshn(:,:) ! before <-- now 267 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 268 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 269 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 270 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) - rnf_b(:,:) + rnf(:,:) ) * ssmask(:,:) 271 sshn(:,:) = ssha(:,:) ! now <-- after 272 ENDIF 273 ! 274 ! Update velocity at AGRIF zoom boundaries 275 #if defined key_agrif 276 IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt ) 277 #endif 257 ! !== Euler time-stepping: no filter, just swap ==! 258 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. & 259 & ( ln_bt_fw .AND. ln_dynspg_ts ) ) THEN 260 sshb(:,:) = sshn(:,:) ! before <-- now 261 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 262 ! 263 ELSE !== Leap-Frog time-stepping: Asselin filter + swap ==! 264 ! ! before <-- now filtered 265 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) 266 IF( .NOT.ln_linssh ) THEN ! before <-- with forcing removed 267 zcoef = atfp * rdt * r1_rau0 268 sshb(:,:) = sshb(:,:) - zcoef * ( emp_b(:,:) - emp (:,:) & 269 & - rnf_b(:,:) + rnf (:,:) & 270 & + fwfisf_b(:,:) - fwfisf(:,:) ) * ssmask(:,:) 271 ENDIF 272 sshn(:,:) = ssha(:,:) ! now <-- after 273 ENDIF 278 274 ! 279 275 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 280 276 ! 281 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')277 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp') 282 278 ! 283 279 END SUBROUTINE ssh_swp
Note: See TracChangeset
for help on using the changeset viewer.