MODULE dynspg_ts !!====================================================================== !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code !! - ! 2005-11 (V. Garnier, G. Madec) optimization !! - ! 2006-08 (S. Masson) distributed restart using iom !! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines !! - ! 2008-01 (R. Benshila) change averaging method !! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation !! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations !! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b !!--------------------------------------------------------------------- #if defined key_dynspg_ts || defined key_esopa !!---------------------------------------------------------------------- !! 'key_dynspg_ts' free surface cst volume with time splitting !!---------------------------------------------------------------------- !! dyn_spg_ts : compute surface pressure gradient trend using a time- !! splitting scheme and add to the general trend !! ts_rst : read/write the time-splitting restart fields in the ocean restart file !!---------------------------------------------------------------------- USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE sbc_oce ! surface boundary condition: ocean USE dynspg_oce ! surface pressure gradient variables USE phycst ! physical constants USE domvvl ! variable volume USE zdfbfr ! bottom friction USE dynvor ! vorticity term USE obc_oce ! Lateral open boundary condition USE obc_par ! open boundary condition parameters USE obcdta ! open boundary condition data USE obcfla ! Flather open boundary condition USE bdy_par ! for lk_bdy USE bdy_oce ! Lateral open boundary condition USE bdydta ! open boundary condition data USE bdydyn2d ! open boundary conditions on barotropic variables USE sbctide USE updtide USE lib_mpp ! distributed memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE in_out_manager ! I/O manager USE iom ! IOM library USE restart ! only for lrst_oce USE zdf_oce ! Vertical diffusion USE wrk_nemo ! Memory Allocation USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC dyn_spg_ts ! routine called by step.F90 PUBLIC ts_rst ! routine called by istate.F90 PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F90 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocity REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocity !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/OPA 4.0 , NEMO Consortium (2011) !! $Id$ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dyn_spg_ts_alloc() !!---------------------------------------------------------------------- !! *** routine dyn_spg_ts_alloc *** !!---------------------------------------------------------------------- ALLOCATE( ftnw (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) , & & ftsw (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc ) ! IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') ! END FUNCTION dyn_spg_ts_alloc SUBROUTINE dyn_spg_ts( kt ) !!---------------------------------------------------------------------- !! *** routine dyn_spg_ts *** !! !! ** Purpose : Compute the now trend due to the surface pressure !! gradient in case of free surface formulation with time-splitting. !! Add it to the general trend of momentum equation. !! !! ** Method : Free surface formulation with time-splitting !! -1- Save the vertically integrated trend. This general trend is !! held constant over the barotropic integration. !! The Coriolis force is removed from the general trend as the !! surface gradient and the Coriolis force are updated within !! the barotropic integration. !! -2- Barotropic loop : updates of sea surface height (ssha_e) and !! barotropic velocity (ua_e and va_e) through barotropic !! momentum and continuity integration. Barotropic former !! variables are time averaging over the full barotropic cycle !! (= 2 * baroclinic time step) and saved in uX_b !! and vX_b (X specifying after, now or before). !! -3- The new general trend becomes : !! ua = ua - sum_k(ua)/H + ( un_b - ub_b ) !! !! ** Action : - Update (ua,va) with the surf. pressure gradient trend !! !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL !!--------------------------------------------------------------------- ! INTEGER, INTENT(in) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: icycle ! local scalar INTEGER :: ikbu, ikbv ! local scalar REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars REAL(wp) :: z1_8, zx1, zy1 ! - - REAL(wp) :: z1_4, zx2, zy2 ! - - REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - REAL(wp) :: ua_btm, va_btm ! - - ! REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum !!---------------------------------------------------------------------- ! IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') ! CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) ! IF( kt == nit000 ) THEN !* initialisation ! IF(lwp) WRITE(numout,*) IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend' IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting' IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2*nn_baro ! CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b ! ua_e (:,:) = un_b (:,:) va_e (:,:) = vn_b (:,:) hu_e (:,:) = hu (:,:) hv_e (:,:) = hv (:,:) hur_e (:,:) = hur (:,:) hvr_e (:,:) = hvr (:,:) IF( ln_dynvor_een ) THEN ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp DO jj = 2, jpj DO ji = fs_2, jpi ! vector opt. ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp END DO END DO ENDIF ! ENDIF ! !* Local constant initialization z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic ! time step (euler timestep) z1_8 = 0.125_wp ! coefficient for vorticity estimates z1_4 = 0.25_wp zraur = 1._wp / rau0 ! 1 / volumic mass ! zhdiv(:,:) = 0._wp ! barotropic divergence zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) zv_sld = 0._wp ; zv_asp = 0._wp IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction z2dt_bf = rdt ELSE z2dt_bf = 2.0_wp * rdt ENDIF ! ----------------------------------------------------------------------------- ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) ! ----------------------------------------------------------------------------- ! ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) ! ! -------------------------- zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp ! DO jk = 1, jpkm1 #if defined key_vectopt_loop DO jj = 1, 1 !Vector opt. => forced unrolling DO ji = 1, jpij #else DO jj = 1, jpj DO ji = 1, jpi #endif ! ! now trend zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) ! ! now velocity zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk) zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk) ! #if defined key_vvl ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk) vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk) #else ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) #endif END DO END DO END DO ! !* baroclinic momentum trend (remove the vertical mean trend) DO jk = 1, jpkm1 ! -------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) END DO END DO END DO ! !* barotropic Coriolis trends * H (vorticity scheme dependent) ! ! ---------------------------==== zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport zwy(:,:) = zvn(:,:) * e1v(:,:) ! IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj) zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) ! energy conserving formulation for planetary vorticity term zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) END DO END DO ! ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) END DO END DO ! ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) zcv(ji,jj) = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) END DO END DO ! ENDIF ! !* Right-Hand-Side of the barotropic momentum equation ! ! ---------------------------------------------------- IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) END DO END DO ENDIF DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend DO ji = fs_2, fs_jpim1 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) END DO END DO ! ! Remove barotropic contribution of bottom friction ! ! from the barotropic transport trend zcoef = -1._wp * z1_2dt_b IF(ln_bfrimp) THEN ! ! Remove the bottom stress trend from 3-D sea surface level gradient ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ikbu = mbku(ji,jj) ikbv = mbkv(ji,jj) ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm END DO END DO ELSE # if defined key_vectopt_loop DO jj = 1, 1 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) # else DO jj = 2, jpjm1 DO ji = 2, jpim1 # endif ! Apply stability criteria for bottom friction !RBbug for vvl and external mode we may need to use varying fse3 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) END DO END DO IF( lk_vvl ) THEN DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) END DO END DO ELSE DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) END DO END DO ENDIF END IF ! end (ln_bfrimp) ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) ! ! -------------------------- zua(:,:) = zua(:,:) * hur(:,:) zva(:,:) = zva(:,:) * hvr(:,:) ! IF( lk_vvl ) THEN ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) ELSE ub_b(:,:) = ub_b(:,:) * hur(:,:) vb_b(:,:) = vb_b(:,:) * hvr(:,:) ENDIF ! ----------------------------------------------------------------------- ! Phase 2 : Integration of the barotropic equations with time splitting ! ----------------------------------------------------------------------- ! ! ! ==================== ! ! ! Initialisations ! ! ! ==================== ! icycle = 2 * nn_baro ! Number of barotropic sub time-step ! ! Start from NOW field hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points hv_e (:,:) = hv (:,:) hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points hvr_e (:,:) = hvr (:,:) !RBbug zsshb_e(:,:) = sshn (:,:) zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now) sshn_e (:,:) = sshn (:,:) zun_e (:,:) = un_b (:,:) ! barotropic velocity (external) zvn_e (:,:) = vn_b (:,:) zub_e (:,:) = un_b (:,:) zvb_e (:,:) = vn_b (:,:) zu_sum (:,:) = un_b (:,:) ! summation zv_sum (:,:) = vn_b (:,:) zssh_sum(:,:) = sshn (:,:) #if defined key_obc ! set ssh corrections to 0 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp IF( lp_obc_west ) sshfow_b(:,:) = 0._wp IF( lp_obc_south ) sshfos_b(:,:) = 0._wp IF( lp_obc_north ) sshfon_b(:,:) = 0._wp #endif ! ! ==================== ! DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1) ! ! ==================== ! z2dt_e = 2. * ( rdt / nn_baro ) IF( jn == 1 ) z2dt_e = rdt / nn_baro ! !* Update the forcing (BDY and tides) ! ! ------------------ IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) IF ( ln_tide_pot ) CALL upd_tide( kt, jn ) ! !* after ssh_e ! ! ----------- DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports DO ji = fs_2, fs_jpim1 ! vector opt. zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & & - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) & & + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) & & - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) END DO END DO ! #if defined key_obc ! ! OBC : zhdiv must be zero behind the open boundary !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south #endif #if defined key_bdy zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask #endif ! DO jj = 2, jpjm1 ! leap-frog on ssh_e DO ji = fs_2, fs_jpim1 ! vector opt. ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) END DO END DO ! !* after barotropic velocities (vorticity scheme dependent) ! ! --------------------------- zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) ! IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! surface pressure gradient IF( lk_vvl) THEN zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) ELSE zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) ENDIF ! add tidal astronomical forcing IF ( ln_tide_pot ) THEN zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) ENDIF ! energy conserving formulation for planetary vorticity term zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) ! after velocities with implicit bottom friction IF( ln_bfrimp ) THEN ! implicit bottom friction ! A new method to implement the implicit bottom friction. ! H. Liu ! Sept 2011 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) ! va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) ! ELSE ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ENDIF END DO END DO ! ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! surface pressure gradient IF( lk_vvl) THEN zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) ELSE zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) ENDIF ! add tidal astronomical forcing IF ( ln_tide_pot ) THEN zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) ENDIF ! enstrophy conserving formulation for planetary vorticity term zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) ! after velocities with implicit bottom friction IF( ln_bfrimp ) THEN ! A new method to implement the implicit bottom friction. ! H. Liu ! Sept 2011 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) ! va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) ! ELSE ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ENDIF END DO END DO ! ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! surface pressure gradient IF( lk_vvl) THEN zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) ELSE zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) ENDIF ! add tidal astronomical forcing IF ( ln_tide_pot ) THEN zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) ENDIF ! energy/enstrophy conserving formulation for planetary vorticity term zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) zv_cor = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) ! after velocities with implicit bottom friction IF( ln_bfrimp ) THEN ! A new method to implement the implicit bottom friction. ! H. Liu ! Sept 2011 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) ! va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) ! ELSE ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ENDIF END DO END DO ! ENDIF ! !* domain lateral boundary ! ! ----------------------- ! OBC open boundaries IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) ! BDY open boundaries #if defined key_bdy pssh => sshn_e phur => hur_e phvr => hvr_e pu2d => ua_e pv2d => va_e IF( lk_bdy ) CALL bdy_dyn2d( kt ) #endif ! CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries CALL lbc_lnk( va_e , 'V', -1. ) CALL lbc_lnk( ssha_e, 'T', 1. ) zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps zv_sum (:,:) = zv_sum (:,:) + va_e (:,:) zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) ! !* Time filter and swap ! ! -------------------- IF( jn == 1 ) THEN ! Swap only (1st Euler time step) zsshb_e(:,:) = sshn_e(:,:) zub_e (:,:) = zun_e (:,:) zvb_e (:,:) = zvn_e (:,:) sshn_e (:,:) = ssha_e(:,:) zun_e (:,:) = ua_e (:,:) zvn_e (:,:) = va_e (:,:) ELSE ! Swap + Filter zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) sshn_e (:,:) = ssha_e(:,:) zun_e (:,:) = ua_e (:,:) zvn_e (:,:) = va_e (:,:) ENDIF IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) ! ! ------------------ DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points DO ji = 1, fs_jpim1 ! Vector opt. zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) END DO END DO CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions CALL lbc_lnk( zsshvn_e, 'V', 1. ) ! hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) ! ENDIF ! ! ==================== ! END DO ! end loop ! ! ! ==================== ! #if defined key_obc IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ????? IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:) IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:) IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:) #endif ! ----------------------------------------------------------------------------- ! Phase 3. update the general trend with the barotropic trend ! ----------------------------------------------------------------------------- ! ! !* Time average ==> after barotropic u, v, ssh zcoef = 1._wp / ( 2 * nn_baro + 1 ) zu_sum(:,:) = zcoef * zu_sum (:,:) zv_sum(:,:) = zcoef * zv_sum (:,:) ! ! !* update the general momentum trend DO jk=1,jpkm1 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b END DO un_b (:,:) = zu_sum(:,:) vn_b (:,:) = zv_sum(:,:) sshn_b(:,:) = zcoef * zssh_sum(:,:) ! ! !* write time-spliting arrays in the restart IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) ! CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) ! IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') ! END SUBROUTINE dyn_spg_ts SUBROUTINE ts_rst( kt, cdrw ) !!--------------------------------------------------------------------- !! *** ROUTINE ts_rst *** !! !! ** Purpose : Read or write time-splitting arrays in restart file !!---------------------------------------------------------------------- INTEGER , INTENT(in) :: kt ! ocean time-step CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag ! INTEGER :: ji, jk ! dummy loop indices !!---------------------------------------------------------------------- ! IF( TRIM(cdrw) == 'READ' ) THEN IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop ELSE un_b (:,:) = 0._wp vn_b (:,:) = 0._wp ! vertical sum IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) END DO ENDIF un_b (:,:) = un_b(:,:) * hur(:,:) vn_b (:,:) = vn_b(:,:) * hvr(:,:) ENDIF ! Vertically integrated velocity (before) IF (neuler/=0) THEN ub_b (:,:) = 0._wp vb_b (:,:) = 0._wp ! vertical sum IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) END DO ENDIF IF( lk_vvl ) THEN ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) ELSE ub_b(:,:) = ub_b(:,:) * hur(:,:) vb_b(:,:) = vb_b(:,:) * hvr(:,:) ENDIF ELSE ! neuler==0 ub_b (:,:) = un_b (:,:) vb_b (:,:) = vn_b (:,:) ENDIF IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered ssh ELSE sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value ENDIF ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) ! ENDIF ! END SUBROUTINE ts_rst #else !!---------------------------------------------------------------------- !! Default case : Empty module No standart free surface cst volume !!---------------------------------------------------------------------- CONTAINS INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function dyn_spg_ts_alloc = 0 END FUNCTION dyn_spg_ts_alloc SUBROUTINE dyn_spg_ts( kt ) ! Empty routine INTEGER, INTENT(in) :: kt WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt END SUBROUTINE dyn_spg_ts SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine INTEGER , INTENT(in) :: kt ! ocean time-step CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw END SUBROUTINE ts_rst #endif !!====================================================================== END MODULE dynspg_ts