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 obcdta ! open boundary condition data USE obcfla ! Flather open boundary condition USE dynvor ! vorticity term USE obc_oce ! Lateral open boundary condition USE obc_par ! open boundary condition parameters USE bdy_oce ! unstructured open boundaries USE bdy_par ! unstructured open boundaries USE bdydta ! unstructured open boundaries USE bdydyn ! unstructured open boundaries USE bdytides ! tidal forcing at unstructured open boundaries. 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 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 !!--------------------------------------------------------------------- USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released USE wrk_nemo, ONLY: zsshun_e => wrk_2d_1 , zsshb_e => wrk_2d_2 , zhdiv => wrk_2d_3 USE wrk_nemo, ONLY: zsshvn_e => wrk_2d_4 , zssh_sum => wrk_2d_5 USE wrk_nemo, ONLY: zcu => wrk_2d_6 , zwx => wrk_2d_7 , zua => wrk_2d_8 , zbfru => wrk_2d_9 USE wrk_nemo, ONLY: zcv => wrk_2d_10 , zwy => wrk_2d_11 , zva => wrk_2d_12 , zbfrv => wrk_2d_13 USE wrk_nemo, ONLY: zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17 USE wrk_nemo, ONLY: zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21 ! INTEGER, INTENT(in) :: kt ! ocean time-step index ! INTEGER :: ji, jj, jk, jn ! dummy loop indices INTEGER :: icycle ! local scalar REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! 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 ! - - !!---------------------------------------------------------------------- IF( wrk_in_use(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & & 11,12,13,14,15,16,17,18,19,20,21 ) ) THEN CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' ) ; RETURN ENDIF 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.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 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. ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. END DO END DO ENDIF ! ENDIF ! !* Local constant initialization z2dt_b = 2.0 * rdt ! baroclinic time step z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates z1_4 = 0.5 * 0.5 zraur = 1. / rau0 ! 1 / volumic mass ! zhdiv(:,:) = 0.e0 ! barotropic divergence zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F) zv_sld = 0.e0 ; zv_asp = 0.e0 ! ----------------------------------------------------------------------------- ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) ! ----------------------------------------------------------------------------- ! ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) ! ! -------------------------- zua(:,:) = 0.e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0 zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0 ! 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_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk) vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(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. / z2dt_b # 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.e0 - 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.e0 - 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 ! !* 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.e0 - umask(:,:,1) ) vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - 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.e0 IF( lp_obc_west ) sshfow_b(:,:) = 0.e0 IF( lp_obc_south ) sshfos_b(:,:) = 0.e0 IF( lp_obc_north ) sshfon_b(:,:) = 0.e0 #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 (OBC, BDY and tides) ! ! ------------------ IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) IF( lk_bdy ) CALL bdy_dta_fla( kt, jn+1, icycle ) ! !* 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.e0 ! east IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! 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 ! 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 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.e0 - 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.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 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 ! 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 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.e0 - 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.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 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 ! 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 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.e0 - 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.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) END DO END DO ! ENDIF ! !* domain lateral boundary ! ! ----------------------- ! ! Flather's boundary condition for the barotropic loop : ! ! - Update sea surface height on each open boundary ! ! - Correct the velocity IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) IF( lk_bdy .OR. ln_tides ) CALL bdy_dyn_fla( sshn_e ) ! 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 * 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 * 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.e0 - umask(:,:,1) ) hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - 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.e0 / ( 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(:,:) ) / z2dt_b va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_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' ) ! ! IF( wrk_not_released(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & & 11,12,13,14,15,16,17,18,19,20,21) ) & CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays') ! 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.e0 vn_b (:,:) = 0.e0 ! 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.e0 vb_b (:,:) = 0.e0 ! 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.e0 - umask(:,:,1) ) vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - 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