MODULE sshwzv !!============================================================================== !! *** MODULE sshwzv *** !! Ocean dynamics : sea surface height and vertical velocity !!============================================================================== !! History : 3.1 ! 2009-02 (G. Madec, M. Leclair) Original code !!---------------------------------------------------------------------- !!---------------------------------------------------------------------- !! ssh_wzv : after ssh & now vertical velocity !! ssh_nxt : filter ans swap the ssh arrays !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers variables USE dom_oce ! ocean space and time domain variables USE sbc_oce ! surface boundary condition: ocean USE domvvl ! Variable volume USE divcur ! hor. divergence and curl (div & cur routines) USE cla_div ! cross land: hor. divergence (div_cla routine) USE iom ! I/O library USE restart ! only for lrst_oce USE in_out_manager ! I/O manager USE prtctl ! Print control USE phycst USE lbclnk ! ocean lateral boundary condition (or mpp link) USE obc_par ! open boundary cond. parameter USE obc_oce USE diaar5, ONLY : lk_diaar5 USE iom IMPLICIT NONE PRIVATE PUBLIC ssh_wzv ! called by step.F90 PUBLIC ssh_nxt ! called by step.F90 !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) !! $Id$ !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ssh_wzv( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE ssh_wzv *** !! !! ** Purpose : compute the after ssh (ssha), the now vertical velocity !! and update the now vertical coordinate (lk_vvl=T). !! !! ** Method : - !! !! - Using the incompressibility hypothesis, the vertical !! velocity is computed by integrating the horizontal divergence !! from the bottom to the surface minus the scale factor evolution. !! The boundary conditions are w=0 at the bottom (no flux) and. !! !! ** action : ssha : after sea surface height !! wn : now vertical velocity !! if lk_vvl=T: sshu_a, sshv_a, sshf_a : after sea surface height !! at u-, v-, f-point s !! hu, hv, hur, hvr : ocean depth and its inverse at u-,v-points !!---------------------------------------------------------------------- USE oce, ONLY : z3d => ta ! use ta as 3D workspace !! INTEGER, INTENT(in) :: kt ! time step !! INTEGER :: ji, jj, jk ! dummy loop indices REAL(wp) :: zcoefu, zcoefv, zcoeff ! temporary scalars REAL(wp) :: z2dt, zraur ! temporary scalars REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ssh_wzv : after sea surface height and now vertical velocity ' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ! wn(:,:,jpk) = 0.e0 ! bottom boundary condition: w=0 (set once for all) ! IF( lk_vvl ) THEN ! before and now Sea SSH at u-, v-, f-points (vvl case only) DO jj = 1, jpjm1 DO ji = 1, jpim1 ! caution: use of Vector Opt. not possible zcoefu = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) zcoefv = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1) sshu_b(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshb(ji ,jj) & & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) ) sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshb(ji,jj ) & & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) ) sshf_b(ji,jj) = zcoeff * ( sshb(ji ,jj) + sshb(ji ,jj+1) & & + sshb(ji+1,jj) + sshb(ji+1,jj+1) ) sshu_n(ji,jj) = zcoefu * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn(ji ,jj) & & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) ) sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn(ji,jj ) & & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) ) sshf_n(ji,jj) = zcoeff * ( sshn(ji ,jj) + sshn(ji ,jj+1) & & + sshn(ji+1,jj) + sshn(ji+1,jj+1) ) END DO END DO CALL lbc_lnk( sshu_b, 'U', 1. ) ; CALL lbc_lnk( sshu_n, 'U', 1. ) CALL lbc_lnk( sshv_b, 'V', 1. ) ; CALL lbc_lnk( sshv_n, 'V', 1. ) CALL lbc_lnk( sshf_b, 'F', 1. ) ; CALL lbc_lnk( sshf_n, 'F', 1. ) ENDIF ! ENDIF ! !------------------------------! IF( lk_vvl ) THEN ! Update Now Vertical coord. ! (only in vvl case) ! !------------------------------! DO jk = 1, jpkm1 fsdept(:,:,jk) = fsdept_n(:,:,jk) ! now local depths stored in fsdep. arrays fsdepw(:,:,jk) = fsdepw_n(:,:,jk) fsde3w(:,:,jk) = fsde3w_n(:,:,jk) ! fse3t (:,:,jk) = fse3t_n (:,:,jk) ! vertical scale factors stored in fse3. arrays fse3u (:,:,jk) = fse3u_n (:,:,jk) fse3v (:,:,jk) = fse3v_n (:,:,jk) fse3f (:,:,jk) = fse3f_n (:,:,jk) fse3w (:,:,jk) = fse3w_n (:,:,jk) fse3uw(:,:,jk) = fse3uw_n(:,:,jk) fse3vw(:,:,jk) = fse3vw_n(:,:,jk) END DO ! ! now ocean depth (at u- and v-points) hu(:,:) = hu_0(:,:) + sshu_n(:,:) hv(:,:) = hv_0(:,:) + sshv_n(:,:) ! ! now masked inverse of the ocean depth (at u- and v-points) hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1.e0 - umask(:,:,1) ) hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1.e0 - vmask(:,:,1) ) ! ENDIF CALL div_cur( kt ) ! Horizontal divergence & Relative vorticity IF( n_cla == 1 ) CALL div_cla( kt ) ! Cross Land Advection (Update Hor. divergence) ! set time step size (Euler/Leapfrog) z2dt = 2. * rdt IF( neuler == 0 .AND. kt == nit000 ) z2dt =rdt zraur = 1. / rau0 ! !------------------------------! ! ! After Sea Surface Height ! ! !------------------------------! zhdiv(:,:) = 0.e0 DO jk = 1, jpkm1 ! Horizontal divergence of barotropic transports zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk) END DO ! ! Sea surface elevation time stepping ssha(:,:) = ( sshb(:,:) - z2dt * ( zraur * ( emp(:,:)-rnf(:,:) ) + zhdiv(:,:) ) ) * tmask(:,:,1) #if defined key_obc IF ( Agrif_Root() ) THEN ssha(:,:) = ssha(:,:) * obctmsk(:,:) CALL lbc_lnk( ssha, 'T', 1. ) ! absolutly compulsory !! (jmm) ENDIF #endif ! ! Sea Surface Height at u-,v- and f-points (vvl case only) IF( lk_vvl ) THEN ! (required only in key_vvl case) DO jj = 1, jpjm1 DO ji = 1, jpim1 ! NO Vector Opt. sshu_a(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) sshv_a(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) sshf_a(ji,jj) = 0.25 * umask(ji,jj,1) * umask (ji,jj+1,1) & & * ( ssha(ji ,jj) + ssha(ji ,jj+1) & & + ssha(ji+1,jj) + ssha(ji+1,jj+1) ) END DO END DO CALL lbc_lnk( sshu_a, 'U', 1. ) ! Boundaries conditions CALL lbc_lnk( sshv_a, 'V', 1. ) CALL lbc_lnk( sshf_a, 'F', 1. ) ENDIF ! !------------------------------! ! ! Now Vertical Velocity ! ! !------------------------------! ! ! integrate from the bottom the hor. divergence DO jk = jpkm1, 1, -1 wn(:,:,jk) = wn(:,:,jk+1) - fse3t_n(:,:,jk) * hdivn(:,:,jk) & & - ( fse3t_a(:,:,jk) & & - fse3t_b(:,:,jk) ) * tmask(:,:,jk) / z2dt END DO ! CALL iom_put( "woce", wn ) ! vertical velocity CALL iom_put( "ssh" , sshn ) ! sea surface height CALL iom_put( "ssh2", sshn(:,:) * sshn(:,:) ) ! square of sea surface height IF( lk_diaar5 ) THEN z2d(:,:) = rau0 * e1t(:,:) * e2t(:,:) DO jk = 1, jpk z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:) END DO CALL iom_put( "w_masstr" , z3d ) ! vertical mass transport CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) ) ! square of vertical mass transport ENDIF ! END SUBROUTINE ssh_wzv SUBROUTINE ssh_nxt( kt ) !!---------------------------------------------------------------------- !! *** ROUTINE ssh_nxt *** !! !! ** Purpose : achieve the sea surface height time stepping by !! applying Asselin time filter and swapping the arrays !! ssha already computed in ssh_wzv !! !! ** Method : - apply Asselin time fiter to now ssh and swap : !! sshn = ssha + atfp * ( sshb -2 sshn + ssha ) !! sshn = ssha !! !! ** action : - sshb, sshn : before & now sea surface height !! ready for the next time step !!---------------------------------------------------------------------- INTEGER, INTENT( in ) :: kt ! ocean time-step index !! INTEGER :: ji, jj ! dummy loop indices !!---------------------------------------------------------------------- IF( kt == nit000 ) THEN IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'ssh_nxt : next sea surface height (Asselin time filter + swap)' IF(lwp) WRITE(numout,*) '~~~~~~~ ' ENDIF ! Time filter and swap of the ssh ! ------------------------------- IF( lk_vvl ) THEN ! Variable volume levels : ssh at t-, u-, v, f-points ! ! ---------------------- ! IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter sshn (:,:) = ssha (:,:) ! now <-- after (before already = now) sshu_n(:,:) = sshu_a(:,:) sshv_n(:,:) = sshv_a(:,:) sshf_n(:,:) = sshf_a(:,:) ELSE ! Leap-Frog time-stepping: Asselin filter + swap DO jj = 1, jpj DO ji = 1, jpi ! before <-- now filtered sshb (ji,jj) = sshn(ji,jj) + atfp * ( sshb (ji,jj) - 2 * sshn (ji,jj) + ssha (ji,jj) ) sshu_b(ji,jj) = sshu_n(ji,jj) + atfp * ( sshu_b(ji,jj) - 2 * sshu_n(ji,jj) + sshu_a(ji,jj) ) sshv_b(ji,jj) = sshv_n(ji,jj) + atfp * ( sshv_b(ji,jj) - 2 * sshv_n(ji,jj) + sshv_a(ji,jj) ) sshf_b(ji,jj) = sshf_n(ji,jj) + atfp * ( sshf_b(ji,jj) - 2 * sshf_n(ji,jj) + sshf_a(ji,jj) ) sshn (ji,jj) = ssha (ji,jj) ! now <-- after sshu_n(ji,jj) = sshu_a(ji,jj) sshv_n(ji,jj) = sshv_a(ji,jj) sshf_n(ji,jj) = sshf_a(ji,jj) END DO END DO ENDIF ! ELSE ! fixed levels : ssh at t-point only ! ! ------------ ! IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler time-stepping at first time-step : no filter sshn(:,:) = ssha(:,:) ! now <-- after (before already = now) ! ELSE ! Leap-Frog time-stepping: Asselin filter + swap DO jj = 1, jpj DO ji = 1, jpi ! before <-- now filtered sshb(ji,jj) = sshn(ji,jj) + atfp * ( sshb(ji,jj) - 2 * sshn(ji,jj) + ssha(ji,jj) ) sshn(ji,jj) = ssha(ji,jj) ! now <-- after END DO END DO ENDIF ! ENDIF ! IF(ln_ctl) CALL prt_ctl(tab2d_1=sshb , clinfo1=' sshb - : ', mask1=tmask, ovlap=1 ) ! END SUBROUTINE ssh_nxt !!====================================================================== END MODULE sshwzv