- Timestamp:
- 2016-01-08T10:35:19+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/DYN/sshwzv.F90
r4486 r6225 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 iom ! I/O library 24 USE restart ! only for lrst_oce 25 USE in_out_manager ! I/O manager 26 USE prtctl ! Print control 27 USE phycst 28 USE lbclnk ! ocean lateral boundary condition (or mpp link) 29 USE lib_mpp ! MPP library 30 USE bdy_oce 31 USE bdy_par 32 USE bdydyn2d ! bdy_ssh routine 33 USE diaar5, ONLY: lk_diaar5 34 USE iom 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 35 27 #if defined key_agrif 36 USE agrif_opa_update37 28 USE agrif_opa_interp 38 29 #endif 39 30 #if defined key_asminc 40 USE asminc ! Assimilation increment 41 #endif 42 USE wrk_nemo ! Memory Allocation 43 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 44 42 45 43 IMPLICIT NONE … … 51 49 52 50 !! * Substitutions 53 # include "domzgr_substitute.h90"54 51 # include "vectopt_loop_substitute.h90" 55 52 !!---------------------------------------------------------------------- … … 70 67 !! by the time step. 71 68 !! 72 !! ** action : ssha :after sea surface height69 !! ** action : ssha, after sea surface height 73 70 !! 74 71 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 75 72 !!---------------------------------------------------------------------- 76 ! 77 REAL(wp), POINTER, DIMENSION(:,: ) :: zhdiv 78 INTEGER, INTENT(in) :: kt ! time step 73 INTEGER, INTENT(in) :: kt ! time step 79 74 ! 80 INTEGER :: jk ! dummy loop indice 81 REAL(wp) :: z2dt, z1_rau0 ! local scalars 82 !!---------------------------------------------------------------------- 83 ! 84 IF( nn_timing == 1 ) CALL timing_start('ssh_nxt') 85 ! 86 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 ) 87 83 ! 88 84 IF( kt == nit000 ) THEN 89 !90 85 IF(lwp) WRITE(numout,*) 91 86 IF(lwp) WRITE(numout,*) 'ssh_nxt : after sea surface height' 92 87 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 93 ! 94 ENDIF 95 ! 96 CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity 88 ENDIF 89 ! 90 CALL div_hor( kt ) ! Horizontal divergence 97 91 ! 98 92 z2dt = 2._wp * rdt ! set time step size (Euler/Leapfrog) … … 104 98 zhdiv(:,:) = 0._wp 105 99 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports 106 zhdiv(:,:) = zhdiv(:,:) + fse3t_n(:,:,jk) * hdivn(:,:,jk)100 zhdiv(:,:) = zhdiv(:,:) + e3t_n(:,:,jk) * hdivn(:,:,jk) 107 101 END DO 108 102 ! ! Sea surface elevation time stepping … … 110 104 ! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations. 111 105 ! 112 z1_rau0 = 0.5_wp * r1_rau0 113 ssha(:,:) = ( sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) 114 115 #if ! defined key_dynspg_ts 116 ! These lines are not necessary with time splitting since 117 ! boundary condition on sea level is set during ts loop 118 #if defined key_agrif 119 CALL agrif_ssh( kt ) 120 #endif 121 #if defined key_bdy 122 IF (lk_bdy) THEN 123 CALL lbc_lnk( ssha, 'T', 1. ) ! Not sure that's necessary 124 CALL bdy_ssh( ssha ) ! Duplicate sea level across open boundaries 125 ENDIF 126 #endif 127 #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 128 125 129 126 #if defined key_asminc 130 ! ! Include the IAU weighted SSH increment 131 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 132 128 CALL ssh_asm_inc( kt ) 133 129 ssha(:,:) = ssha(:,:) + z2dt * ssh_iau(:,:) 134 130 ENDIF 135 131 #endif 136 137 132 ! !------------------------------! 138 133 ! ! outputs ! 139 134 ! !------------------------------! 140 CALL iom_put( "ssh" , sshn ) ! sea surface height141 CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height142 135 ! 143 136 IF(ln_ctl) CALL prt_ctl( tab2d_1=ssha, clinfo1=' ssha - : ', mask1=tmask, ovlap=1 ) … … 165 158 !! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling. 166 159 !!---------------------------------------------------------------------- 167 ! 168 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 169 164 REAL(wp), POINTER, DIMENSION(:,: ) :: z2d 170 165 REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d, zhdiv 171 ! 172 INTEGER :: ji, jj, jk ! dummy loop indices 173 REAL(wp) :: z1_2dt ! local scalars 174 !!---------------------------------------------------------------------- 175 176 IF( nn_timing == 1 ) CALL timing_start('wzv') 166 !!---------------------------------------------------------------------- 167 ! 168 IF( nn_timing == 1 ) CALL timing_start('wzv') 177 169 ! 178 170 IF( kt == nit000 ) THEN 179 !180 171 IF(lwp) WRITE(numout,*) 181 172 IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity ' … … 183 174 ! 184 175 wn(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all) 185 !186 176 ENDIF 187 177 ! !------------------------------! … … 199 189 DO jj = 2, jpjm1 200 190 DO ji = fs_2, fs_jpim1 ! vector opt. 201 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) ) 202 192 END DO 203 193 END DO … … 208 198 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 209 199 ! computation of w 210 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk) + zhdiv(:,:,jk)&211 & + 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) 212 202 END DO 213 203 ! IF( ln_vvl_layer ) wn(:,:,:) = 0.e0 … … 216 206 DO jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence 217 207 ! computation of w 218 wn(:,:,jk) = wn(:,:,jk+1) - ( fse3t_n(:,:,jk) * hdivn(:,:,jk)&219 & + 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) 220 210 END DO 221 211 ENDIF 222 212 223 213 #if defined key_bdy 224 IF (lk_bdy) THEN214 IF( lk_bdy ) THEN 225 215 DO jk = 1, jpkm1 226 216 wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:) … … 229 219 #endif 230 220 ! 231 ! !------------------------------!232 ! ! outputs !233 ! !------------------------------!234 CALL iom_put( "woce", wn ) ! vertical velocity235 IF( lk_diaar5 ) THEN ! vertical mass transport & its square value236 CALL wrk_alloc( jpi, jpj, z2d )237 CALL wrk_alloc( jpi, jpj, jpk, z3d )238 ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.239 z2d(:,:) = rau0 * e12t(:,:)240 DO jk = 1, jpk241 z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)242 END DO243 CALL iom_put( "w_masstr" , z3d )244 CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )245 CALL wrk_dealloc( jpi, jpj, z2d )246 CALL wrk_dealloc( jpi, jpj, jpk, z3d )247 ENDIF248 !249 221 IF( nn_timing == 1 ) CALL timing_stop('wzv') 250 251 222 ! 252 223 END SUBROUTINE wzv 224 253 225 254 226 SUBROUTINE ssh_swp( kt ) … … 272 244 !!---------------------------------------------------------------------- 273 245 INTEGER, INTENT(in) :: kt ! ocean time-step index 246 ! 247 REAL(wp) :: zcoef ! local scalar 274 248 !!---------------------------------------------------------------------- 275 249 ! … … 281 255 IF(lwp) WRITE(numout,*) '~~~~~~~ ' 282 256 ENDIF 283 284 # if defined key_dynspg_ts 285 IF( ( neuler == 0 .AND. kt == nit000 ) .OR. ln_bt_fw ) THEN !** Euler time-stepping: no filter 286 # else 287 IF ( neuler == 0 .AND. kt == nit000 ) THEN !** Euler time-stepping at first time-step : no filter 288 #endif 289 sshb(:,:) = sshn(:,:) ! before <-- now 290 sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) 291 ELSE !** Leap-Frog time-stepping: Asselin filter + swap 292 sshb(:,:) = sshn(:,:) + atfp * ( sshb(:,:) - 2 * sshn(:,:) + ssha(:,:) ) ! before <-- now filtered 293 IF( lk_vvl ) sshb(:,:) = sshb(:,:) - atfp * rdt / rau0 * ( emp_b(:,:) - emp(:,:) ) * tmask(:,:,1) 294 sshn(:,:) = ssha(:,:) ! now <-- after 295 ENDIF 296 ! 297 ! Update velocity at AGRIF zoom boundaries 298 #if defined key_agrif 299 IF ( .NOT.Agrif_Root() ) CALL Agrif_Update_Dyn( kt ) 300 #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 301 274 ! 302 275 IF(ln_ctl) CALL prt_ctl( tab2d_1=sshb, clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) 303 276 ! 304 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp')277 IF( nn_timing == 1 ) CALL timing_stop('ssh_swp') 305 278 ! 306 279 END SUBROUTINE ssh_swp
Note: See TracChangeset
for help on using the changeset viewer.