MODULE dynspg_ts !!====================================================================== !! *** MODULE dynspg_ts *** !! Ocean dynamics: surface pressure gradient trend !!====================================================================== !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code !! " " ! 05-11 (V. Garnier, G. Madec) optimization !! 9.0 ! 06-08 (S. Masson) distributed restart using iom !! 9.0 ! 08-01 (R. Benshila) change averaging method !!--------------------------------------------------------------------- #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 !!---------------------------------------------------------------------- !! * Modules used USE oce ! ocean dynamics and tracers USE dom_oce ! ocean space and time domain USE phycst ! physical constants USE ocesbc ! ocean surface boundary condition 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 lib_mpp ! distributed memory computing library USE lbclnk ! ocean lateral boundary conditions (or mpp link) USE prtctl ! Print control USE dynspg_oce ! surface pressure gradient variables USE in_out_manager ! I/O manager USE iom USE restart ! only for lrst_oce USE domvvl ! variable volume IMPLICIT NONE PRIVATE PUBLIC dyn_spg_ts ! routine called by step.F90 PUBLIC ts_rst ! routine called by j-k-i subroutine REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter & ftsw, ftse ! (only used with een vorticity scheme) !! * Substitutions # include "domzgr_substitute.h90" # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! OPA 9.0 , LOCEAN-IPSL (2005) !! $Header: /home/opalod/NEMOCVSROOT/NEMO/OPA_SRC/DYN/dynspg_ts.F90,v 1.16 2007/06/05 10:38:27 opalod Exp $ !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt !!---------------------------------------------------------------------- CONTAINS 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. !! Compute the free surface. !! !! ** 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 transports (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 zsshX_b, zuX_b !! and zvX_b (X specifying after, now or before). !! -3- Update of sea surface height from time averaged barotropic !! variables. !! - apply lateral boundary conditions on sshn. !! -4- The new general trend becomes : !! ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H !! !! ** 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 !! * Local declarations INTEGER :: ji, jj, jk, jit ! dummy loop indices INTEGER :: icycle, ibaro ! temporary scalar REAL(wp) :: & zraur, zcoef, z2dt_e, z2dt_b, zfac25, & ! temporary scalars zfact1, zspgu, zcubt, zx1, zy1, & ! " " zfact2, zspgv, zcvbt, zx2, zy2 ! " " REAL(wp), DIMENSION(jpi,jpj) :: & zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays zua, zva, zub, zvb, & ! " " zssha_b, zua_b, zva_b, & ! " " zsshb_e, zub_e, zvb_e, & ! " " zun_e, zvn_e ! " " !! Variable volume REAL(wp), DIMENSION(jpi,jpj) :: & zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e ! 2D workspace REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace !!---------------------------------------------------------------------- ! Arrays initialization ! --------------------- zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0 zhdiv(:,:) = 0.e0 IF( kt == nit000 ) THEN ! 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 = ', FLOOR( 2*rdt/rdtbt ) CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b ssha_e(:,:) = sshn(:,:) ua_e(:,:) = un_b(:,:) va_e(:,:) = vn_b(:,:) 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 ! Substract the surface pressure gradient from the momentum ! --------------------------------------------------------- IF( lk_vvl ) THEN ! Variable volume ! 0. Local initialization ! ----------------------- zspgu_1(:,:) = 0.e0 zspgv_1(:,:) = 0.e0 ! 1. Surface pressure gradient (now) ! ---------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj) zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) END DO END DO ! 2. Substract the surface pressure trend from the general 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) - zspgu_1(ji,jj) va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) END DO END DO END DO ENDIF ! Local constant initialization ! -------------------------------- z2dt_b = 2.0 * rdt ! baroclinic time step IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt zfact1 = 0.5 * 0.25 ! coefficient for vorticity estimates zfact2 = 0.5 * 0.5 zraur = 1. / rauw ! 1 / volumic mass of pure water ! ----------------------------------------------------------------------------- ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) ! ----------------------------------------------------------------------------- ! Vertically integrated quantities ! -------------------------------- zua(:,:) = 0.e0 zva(:,:) = 0.e0 zub(:,:) = 0.e0 zvb(:,:) = 0.e0 zwx(:,:) = 0.e0 zwy(:,:) = 0.e0 ! vertical sum IF( lk_vopt_loop ) THEN ! vector opt., forced unroll DO jk = 1, jpkm1 DO ji = 1, jpij ! ! Vertically integrated momentum trends zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk) zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) ! ! Vertically integrated transports (before) zub(ji,1) = zub(ji,1) + fse3u(ji,1,jk) * ub(ji,1,jk) zvb(ji,1) = zvb(ji,1) + fse3v(ji,1,jk) * vb(ji,1,jk) ! ! Planetary vorticity transport fluxes (now) zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk) END DO END DO ELSE ! No vector opt. DO jk = 1, jpkm1 ! ! Vertically integrated momentum trends zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) ! ! Vertically integrated transports (before) zub(:,:) = zub(:,:) + fse3u(:,:,jk) * ub(:,:,jk) zvb(:,:) = zvb(:,:) + fse3v(:,:,jk) * vb(:,:,jk) ! ! Planetary vorticity (now) zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) END DO ENDIF 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) = zfact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) zcv(ji,jj) =-zfact2 * ( 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 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 =-zfact1 * ( 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 zfac25 = 0.25 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zcu(ji,jj) = + zfac25 / 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) = - zfac25 / 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 ! Remove barotropic trend from general momentum 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 ! Remove coriolis term from barotropic trend ! ------------------------------------------ DO jj = 2, jpjm1 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 ! ----------------------------------------------------------------------- ! Phase 2 : Integration of the barotropic equations with time splitting ! ----------------------------------------------------------------------- ! Initialisations !---------------- ! Number of iteration of the barotropic loop ibaro = FLOOR( rdt / rdtbt ) icycle = 3 /2 * ibaro ! variables for the barotropic equations zsshb_e(:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now) sshn_e (:,:) = sshn_b(:,:) zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now) zvb_e (:,:) = vn_b (:,:) zun_e (:,:) = un_b (:,:) zvn_e (:,:) = vn_b (:,:) zssha_b(:,:) = 0.e0 zua_b (:,:) = 0.e0 zva_b (:,:) = 0.e0 IF( lk_vvl ) THEN zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point hu_e(:,:) = hu (:,:) ! (barotropic) ocean depth at u-point hv_e(:,:) = hv (:,:) ! (barotropic) ocean depth at v-point zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point ENDIF ! set ssh corrections to 0 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop #if defined key_obc 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 ! Barotropic integration over 2 baroclinic time steps ! --------------------------------------------------- ! ! ==================== ! DO jit = 1, icycle ! sub-time-step loop ! ! ! ==================== ! z2dt_e = 2. * rdtbt IF ( jit == 1 ) z2dt_e = rdtbt ! Time interpolation of open boundary condition data IF( lk_obc ) CALL obc_dta_bt( kt, jit ) ! Horizontal divergence of barotropic transports !-------------------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun_e(ji ,jj) & & -e2u(ji-1,jj ) * zun_e(ji-1,jj) & & +e1v(ji ,jj ) * zvn_e(ji ,jj) & & -e1v(ji ,jj-1) * zvn_e(ji ,jj-1) ) & & / (e1t(ji,jj)*e2t(ji,jj)) END DO END DO #if defined key_obc ! open boundaries (div 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 ! Sea surface height from the barotropic system !---------------------------------------------- DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) END DO END DO ! evolution of the barotropic transport ( following the vorticity scheme used) ! ---------------------------------------------------------------------------- zwx(:,:) = e2u(:,:) * zun_e(:,:) zwy(:,:) = e1v(:,:) * zvn_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 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) ELSE zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(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) zcubt = zfact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) ! after transports ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 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 zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) ELSE zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) ENDIF ! enstrophy conserving formulation for planetary vorticity term zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) ! after transports ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) END DO END DO ! ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme zfac25 = 0.25 DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. ! surface pressure gradient IF( lk_vvl) THEN zspgu = -grav * ( ( rhd(ji+1,jj,1) + 1 ) * sshn_e(ji+1,jj) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj) zspgv = -grav * ( ( rhd(ji,jj+1,1) + 1 ) * sshn_e(ji,jj+1) & & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj) ELSE zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) ENDIF ! energy/enstrophy conserving formulation for planetary vorticity term zcubt = + zfac25 / 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) ) zcvbt = - zfac25 / 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 ) ) ! after transports ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) END DO END DO ! ENDIF ! Flather's boundary condition for the barotropic loop : ! - Update sea surface height on each open boundary ! - Correct the barotropic transports IF( lk_obc ) CALL obc_fla_ts ! ... Boundary conditions on ua_e, va_e, ssha_e CALL lbc_lnk( ua_e , 'U', -1. ) CALL lbc_lnk( va_e , 'V', -1. ) CALL lbc_lnk( ssha_e, 'T', 1. ) ! temporal sum !------------- IF( jit >= ibaro/2 ) THEN zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) zua_b (:,:) = zua_b (:,:) + ua_e (:,:) zva_b (:,:) = zva_b (:,:) + va_e (:,:) ENDIF ! Time filter and swap of dynamics arrays ! --------------------------------------- IF( neuler == 0 .AND. jit == 1 ) THEN ! Euler (forward) time stepping 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 ! Asselin filtering 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 ! Variable volume ! Sea surface elevation ! --------------------- DO jj = 1, jpjm1 DO ji = 1,jpim1 ! Sea Surface Height at u-point before 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) ) ! Sea Surface Height at v-point before 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 ! Boundaries conditions CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) ! Scale factors at before and after time step ! ------------------------------------------- CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ; CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) ! Ocean depth at U- and V-points hu_e(:,:) = 0.e0 hv_e(:,:) = 0.e0 DO jk = 1, jpk hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) END DO ENDIF ! End variable volume case ! ! ==================== ! END DO ! end loop ! ! ! ==================== ! ! Time average of after barotropic variables zcoef = 1.e0 / ( ibaro + 1 ) zssha_b(:,:) = zcoef * zssha_b(:,:) zua_b (:,:) = zcoef * zua_b (:,:) zva_b (:,:) = zcoef * zva_b (:,:) #if defined key_obc IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 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 sea surface height from time averaged barotropic variables ! --------------------------------------------------------------------------- ! Horizontal divergence of time averaged barotropic transports !------------------------------------------------------------- IF( .NOT. lk_vvl ) THEN DO jj = 2, jpjm1 DO ji = fs_2, fs_jpim1 ! vector opt. zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) & & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) & & / ( e1t(ji,jj) * e2t(ji,jj) ) END DO END DO ENDIF #if defined key_obc && ! defined key_vvl ! open boundaries (div 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 ! sea surface height !------------------- IF( lk_vvl ) THEN sshbb(:,:) = sshb(:,:) sshb (:,:) = sshn(:,:) sshn (:,:) = ssha(:,:) ELSE sshb (:,:) = sshn(:,:) sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) ENDIF ! ... Boundary conditions on sshn IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. ) ! ----------------------------------------------------------------------------- ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step) ! ----------------------------------------------------------------------------- ! Swap on time averaged barotropic variables ! ------------------------------------------ sshb_b(:,:) = sshn_b (:,:) IF ( kt == nit000 ) zssha_b(:,:) = sshn(:,:) sshn_b(:,:) = zssha_b(:,:) un_b (:,:) = zua_b (:,:) vn_b (:,:) = zva_b (:,:) ! add time averaged barotropic coriolis and surface pressure gradient ! terms to the general momentum trend ! -------------------------------------------------------------------- DO jk=1,jpkm1 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b END DO ! write filtered free surface arrays in restart file ! -------------------------------------------------- IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) ! print sum trends (used for debugging) IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask ) ! 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, 'sshn', ldstop = .FALSE. ) > 0 ) THEN CALL iom_get( numror, jpdom_autoglo, 'sshb' , sshb(:,:) ) CALL iom_get( numror, jpdom_autoglo, 'sshn' , sshn(:,:) ) IF( neuler == 0 ) sshb(:,:) = sshn(:,:) ELSE IF( nn_rstssh == 1 ) THEN sshb(:,:) = 0.e0 sshn(:,:) = 0.e0 ENDIF ENDIF IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) ) ! free surface issued CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! horizontal transports issued CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) ELSE sshb_b(:,:) = sshb(:,:) sshn_b(:,:) = sshn(:,:) 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 ENDIF ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop ENDIF ! END SUBROUTINE ts_rst #else !!---------------------------------------------------------------------- !! Default case : Empty module No standart free surface cst volume !!---------------------------------------------------------------------- CONTAINS SUBROUTINE dyn_spg_ts( kt ) ! Empty routine 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