- Timestamp:
- 2013-11-20T17:28:04+01:00 (10 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2013/dev_MERGE_2013/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r3680 r4292 9 9 !! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations 10 10 !! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b 11 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 11 13 !!--------------------------------------------------------------------- 12 14 #if defined key_dynspg_ts || defined key_esopa … … 16 18 !! dyn_spg_ts : compute surface pressure gradient trend using a time- 17 19 !! splitting scheme and add to the general trend 18 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file19 20 !!---------------------------------------------------------------------- 20 21 USE oce ! ocean dynamics and tracers … … 24 25 USE phycst ! physical constants 25 26 USE domvvl ! variable volume 26 USE zdfbfr ! bottom friction27 27 USE dynvor ! vorticity term 28 USE obc_oce ! Lateral open boundary condition29 USE obc_par ! open boundary condition parameters30 USE obcdta ! open boundary condition data31 USE obcfla ! Flather open boundary condition32 28 USE bdy_par ! for lk_bdy 33 29 USE bdy_oce ! Lateral open boundary condition 34 USE bdy dta! open boundary condition data30 USE bdytides ! open boundary condition data 35 31 USE bdydyn2d ! open boundary conditions on barotropic variables 36 USE sbctide 37 USE updtide 32 USE sbctide ! tides 33 USE updtide ! tide potential 38 34 USE lib_mpp ! distributed memory computing library 39 35 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 41 37 USE in_out_manager ! I/O manager 42 38 USE iom ! IOM library 39 USE restart ! only for lrst_oce 43 40 USE zdf_oce ! Vertical diffusion 44 41 USE wrk_nemo ! Memory Allocation 45 USE timing ! Timing 42 USE timing ! Timing 43 USE sbcapr ! surface boundary condition: atmospheric pressure 44 USE dynadv, ONLY: ln_dynadv_vec 45 #if defined key_agrif 46 USE agrif_opa_interp ! agrif 47 #endif 46 48 47 49 … … 49 51 PRIVATE 50 52 51 PUBLIC dyn_spg_ts ! routine called by step.F90 52 PUBLIC ts_rst ! routine called by istate.F90 53 PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F90 54 55 53 PUBLIC dyn_spg_ts ! routine called in dynspg.F90 54 PUBLIC dyn_spg_ts_alloc ! " " " " 55 PUBLIC dyn_spg_ts_init ! " " " " 56 57 ! Potential namelist parameters below to be read in dyn_spg_ts_init 58 LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.TRUE. !: Forward integration of barotropic sub-stepping 59 LOGICAL, PRIVATE, PARAMETER :: ln_bt_av=.TRUE. !: Time averaging of barotropic variables 60 LOGICAL, PRIVATE, PARAMETER :: ln_bt_nn_auto=.FALSE. !: Set number of iterations automatically 61 INTEGER, PRIVATE, PARAMETER :: nn_bt_flt=1 !: Filter choice 62 REAL(wp), PRIVATE, PARAMETER :: rn_bt_cmax=0.8_wp !: Max. courant number (used if ln_bt_nn_auto=T) 63 ! End namelist parameters 64 65 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_baro <= 2.5 nn_baro 66 REAL(wp),SAVE :: rdtbt ! Barotropic time step 67 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 69 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 70 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 71 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 56 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 57 74 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 58 75 59 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocity 60 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocity 76 ! Arrays below are saved to allow testing of the "no time averaging" option 77 ! If this option is not retained, these could be replaced by temporary arrays 78 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 79 ubb_e, ub_e, & 80 vbb_e, vb_e 61 81 62 82 !! * Substitutions … … 64 84 # include "vectopt_loop_substitute.h90" 65 85 !!---------------------------------------------------------------------- 66 !! NEMO/OPA 4.0 , NEMO Consortium (2011)67 !! $Id $86 !! NEMO/OPA 3.5 , NEMO Consortium (2013) 87 !! $Id: dynspg_ts.F90 68 88 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 69 89 !!---------------------------------------------------------------------- … … 74 94 !! *** routine dyn_spg_ts_alloc *** 75 95 !!---------------------------------------------------------------------- 76 ALLOCATE( ftnw (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) , & 77 & ftsw (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc ) 78 ! 96 INTEGER :: ierr(3) 97 !!---------------------------------------------------------------------- 98 ierr(:) = 0 99 100 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 101 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 102 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 103 104 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 105 106 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 107 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 108 109 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 110 79 111 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 80 112 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays') … … 82 114 END FUNCTION dyn_spg_ts_alloc 83 115 84 85 116 SUBROUTINE dyn_spg_ts( kt ) 86 117 !!---------------------------------------------------------------------- 87 !! *** routine dyn_spg_ts ***88 118 !! 89 !! ** Purpose : Compute the now trend due to the surface pressure 90 !! gradient in case of free surface formulation with time-splitting. 91 !! Add it to the general trend of momentum equation. 119 !! ** Purpose : 120 !! -Compute the now trend due to the explicit time stepping 121 !! of the quasi-linear barotropic system. Barotropic variables are 122 !! advanced from internal time steps "n" to "n+1" (if ln_bt_cen=F) 123 !! or from "n-1" to "n+1" time steps (if ln_bt_cen=T) with a 124 !! generalized forward-backward (see ref. below) time stepping. 125 !! -Update the free surface at step "n+1" (ssha, zsshu_a, zsshv_a). 126 !! -Compute barotropic advective velocities at step "n" to be used 127 !! to advect tracers latter on. These are compliant with discrete 128 !! continuity equation taken at the baroclinic time steps, thus 129 !! ensuring tracers conservation. 92 130 !! 93 !! ** Method : Free surface formulation with time-splitting 94 !! -1- Save the vertically integrated trend. This general trend is 95 !! held constant over the barotropic integration. 96 !! The Coriolis force is removed from the general trend as the 97 !! surface gradient and the Coriolis force are updated within 98 !! the barotropic integration. 99 !! -2- Barotropic loop : updates of sea surface height (ssha_e) and 100 !! barotropic velocity (ua_e and va_e) through barotropic 101 !! momentum and continuity integration. Barotropic former 102 !! variables are time averaging over the full barotropic cycle 103 !! (= 2 * baroclinic time step) and saved in uX_b 104 !! and vX_b (X specifying after, now or before). 105 !! -3- The new general trend becomes : 106 !! ua = ua - sum_k(ua)/H + ( un_b - ub_b ) 131 !! ** Method : 107 132 !! 108 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend 133 !! ** Action : - Update barotropic velocities: ua_b, va_b 134 !! - Update trend (ua,va) with barotropic component 135 !! - Update ssha, zsshu_a, zsshv_a 136 !! - Update barotropic advective velocity at kt=now 109 137 !! 110 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 138 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 139 !! The regional oceanic modeling system (ROMS): 140 !! a split-explicit, free-surface, 141 !! topography-following-coordinate oceanic model. 142 !! Ocean Modelling, 9, 347-404. 111 143 !!--------------------------------------------------------------------- 112 144 ! 113 145 INTEGER, INTENT(in) :: kt ! ocean time-step index 114 146 ! 115 INTEGER :: ji, jj, jk, jn ! dummy loop indices 116 INTEGER :: icycle ! local scalar 117 INTEGER :: ikbu, ikbv ! local scalar 118 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 119 REAL(wp) :: z1_8, zx1, zy1 ! - - 120 REAL(wp) :: z1_4, zx2, zy2 ! - - 121 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 122 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 123 REAL(wp) :: ua_btm, va_btm ! - - 124 ! 125 REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 126 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 127 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 147 LOGICAL :: ll_fw_start ! if true, forward integration 148 LOGICAL :: ll_init ! if true, special startup of 2d equations 149 INTEGER :: ji, jj, jk, jn ! dummy loop indices 150 INTEGER :: ikbu, ikbv, noffset ! local integers 151 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 152 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 153 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 154 REAL(wp) :: zu_spg, zv_spg ! - - 155 REAL(wp) :: zhura, zhvra ! - - 156 REAL(wp) :: za0, za1, za2, za3 ! - - 157 ! 158 REAL(wp), POINTER, DIMENSION(:,:) :: zun_e, zvn_e, zsshp2_e 159 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 160 REAL(wp), POINTER, DIMENSION(:,:) :: zu_sum, zv_sum, zwx, zwy, zhdiv 161 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 162 REAL(wp), POINTER, DIMENSION(:,:) :: zhur_b, zhvr_b 163 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 164 REAL(wp), POINTER, DIMENSION(:,:) :: zht, zhf 128 165 !!---------------------------------------------------------------------- 129 166 ! 130 167 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 131 168 ! 132 CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 133 CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 134 CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 135 ! 136 IF( kt == nit000 ) THEN !* initialisation 169 ! !* Allocate temporay arrays 170 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 171 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 172 CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 173 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 174 CALL wrk_alloc( jpi, jpj, zhur_b, zhvr_b ) 175 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 176 CALL wrk_alloc( jpi, jpj, zht, zhf ) 177 ! 178 ! !* Local constant initialization 179 z1_12 = 1._wp / 12._wp 180 z1_8 = 0.125_wp 181 z1_4 = 0.25_wp 182 z1_2 = 0.5_wp 183 zraur = 1._wp / rau0 184 ! 185 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 186 z2dt_bf = rdt 187 ELSE 188 z2dt_bf = 2.0_wp * rdt 189 ENDIF 190 z1_2dt_b = 1.0_wp / z2dt_bf 191 ! 192 ll_init = ln_bt_av ! if no time averaging, then no specific restart 193 ll_fw_start = .FALSE. 194 ! 195 ! time offset in steps for bdy data update 196 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 197 ! 198 IF( kt == nit000 ) THEN !* initialisation 137 199 ! 138 200 IF(lwp) WRITE(numout,*) … … 141 203 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2*nn_baro 142 204 ! 143 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b 144 ! 145 ua_e (:,:) = un_b (:,:) 146 va_e (:,:) = vn_b (:,:) 147 hu_e (:,:) = hu (:,:) 148 hv_e (:,:) = hv (:,:) 149 hur_e (:,:) = hur (:,:) 150 hvr_e (:,:) = hvr (:,:) 151 IF( ln_dynvor_een ) THEN 152 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 205 IF (neuler==0) ll_init=.TRUE. 206 ! 207 IF (ln_bt_fw.OR.(neuler==0)) THEN 208 ll_fw_start=.TRUE. 209 noffset = 0 210 ELSE 211 ll_fw_start=.FALSE. 212 ENDIF 213 ! 214 ! Set averaging weights and cycle length: 215 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 216 ! 217 IF ((neuler/=0).AND.(ln_bt_fw)) CALL ts_rst( nit000, 'READ' ) 218 ! 219 ENDIF 220 ! 221 ! Set arrays to remove/compute coriolis trend. 222 ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 223 ! Note that these arrays are also used during barotropic loop. These are however frozen 224 ! although they should be updated in variable volume case. Not a big approximation. 225 ! To remove this approximation, copy lines below inside barotropic loop 226 ! and update depths at T-F points (ht and hf resp.) at each barotropic time step 227 ! 228 IF ( kt == nit000 .OR. lk_vvl ) THEN 229 IF ( ln_dynvor_een ) THEN 230 ! JC: Simplification needed below: define ht_0 even when volume is fixed 231 IF (lk_vvl) THEN 232 zht(:,:) = (ht_0(:,:) + sshn(:,:)) * tmask(:,:,1) 233 ELSE 234 zht(:,:) = 0. 235 DO jk = 1, jpkm1 236 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 237 END DO 238 ENDIF 239 240 DO jj = 1, jpjm1 241 DO ji = 1, jpim1 242 zwz(ji,jj) = ( zht(ji ,jj+1) + zht(ji+1,jj+1) + & 243 & zht(ji ,jj ) + zht(ji+1,jj ) ) & 244 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 245 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 246 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 247 END DO 248 END DO 249 CALL lbc_lnk( zwz, 'F', 1._wp ) 250 zwz(:,:) = ff(:,:) * zwz(:,:) 251 252 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 153 253 DO jj = 2, jpj 154 254 DO ji = fs_2, jpi ! vector opt. 155 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 156 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 157 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 158 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 159 END DO 160 END DO 161 ENDIF 162 ! 163 ENDIF 164 165 ! !* Local constant initialization 166 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step 167 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic 168 ! time step (euler timestep) 169 z1_8 = 0.125_wp ! coefficient for vorticity estimates 170 z1_4 = 0.25_wp 171 zraur = 1._wp / rau0 ! 1 / volumic mass 172 ! 173 zhdiv(:,:) = 0._wp ! barotropic divergence 174 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 175 zv_sld = 0._wp ; zv_asp = 0._wp 176 177 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 178 z2dt_bf = rdt 179 ELSE 180 z2dt_bf = 2.0_wp * rdt 181 ENDIF 182 255 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 256 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 257 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 258 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 259 END DO 260 END DO 261 ELSE 262 zwz(:,:) = 0._wp 263 zht(:,:) = 0. 264 IF ( .not. ln_sco ) THEN 265 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 266 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 267 ! ENDIF 268 ! zht(:,:) = gdepw_0(:,:,jk+1) 269 ELSE 270 zht(:,:) = hbatf(:,:) 271 END IF 272 273 DO jj = 1, jpjm1 274 zht(:,jj) = zht(:,jj)*(1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 275 END DO 276 277 DO jk = 1, jpkm1 278 DO jj = 1, jpjm1 279 zht(:,jj) = zht(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 280 END DO 281 END DO 282 CALL lbc_lnk( zht, 'F', 1._wp ) 283 ! JC: TBC. hf should be greater than 0 284 DO jj = 1, jpj 285 DO ji = 1, jpi 286 IF( zht(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zht(ji,jj) ! zht is actually hf here but it saves an array 287 END DO 288 END DO 289 zwz(:,:) = ff(:,:) * zwz(:,:) 290 ENDIF 291 ENDIF 292 ! 293 ! If forward start at previous time step, and centered integration, 294 ! then update averaging weights: 295 IF ((.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN 296 ll_fw_start=.FALSE. 297 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 298 ENDIF 299 300 ! before inverse water column height at u- and v- points 301 IF( lk_vvl ) THEN 302 zhur_b(:,:) = 0. 303 zhvr_b(:,:) = 0. 304 DO jk = 1, jpk 305 zhur_b(:,:) = zhur_b(:,:) + fse3u_b(:,:,jk) * umask(:,:,jk) 306 zhvr_b(:,:) = zhvr_b(:,:) + fse3v_b(:,:,jk) * vmask(:,:,jk) 307 END DO 308 zhur_b(:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1. - umask(:,:,1) ) 309 zhvr_b(:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1. - vmask(:,:,1) ) 310 ELSE 311 zhur_b(:,:) = hur(:,:) 312 zhvr_b(:,:) = hvr(:,:) 313 ENDIF 314 183 315 ! ----------------------------------------------------------------------------- 184 316 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 185 317 ! ----------------------------------------------------------------------------- 186 318 ! 319 ! Some vertical sums (at now and before time steps) below could be suppressed 320 ! if one swap barotropic arrays somewhere 321 ! 187 322 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 188 ! ! -------------------------- 189 zu a(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp190 zv a(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp323 ! ! -------------------------------------------------- 324 zu_frc(:,:) = 0._wp ; ub_b(:,:) = 0._wp ; un_b(:,:) = 0._wp 325 zv_frc(:,:) = 0._wp ; vb_b(:,:) = 0._wp ; vn_b(:,:) = 0._wp 191 326 ! 192 327 DO jk = 1, jpkm1 … … 198 333 DO ji = 1, jpi 199 334 #endif 200 ! ! now trend 201 zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 202 zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 203 ! ! now velocity 204 zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk) 205 zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk) 206 ! 207 #if defined key_vvl 208 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk) 209 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk) 210 #else 211 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 212 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 213 #endif 335 ! ! now trend: 336 zu_frc(ji,jj) = zu_frc(ji,jj) + fse3u(ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 337 zv_frc(ji,jj) = zv_frc(ji,jj) + fse3v(ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 338 ! ! now bt transp: 339 un_b(ji,jj) = un_b(ji,jj) + fse3u(ji,jj,jk) * un(ji,jj,jk) * umask(ji,jj,jk) 340 vn_b(ji,jj) = vn_b(ji,jj) + fse3v(ji,jj,jk) * vn(ji,jj,jk) * vmask(ji,jj,jk) 341 ! ! before bt transp: 342 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) 343 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk) 214 344 END DO 215 345 END DO 216 346 END DO 217 347 ! 348 zu_frc(:,:) = zu_frc(:,:) * hur(:,:) 349 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:) 350 ! 351 IF( lk_vvl ) THEN 352 ub_b(:,:) = ub_b(:,:) * zhur_b(:,:) 353 vb_b(:,:) = vb_b(:,:) * zhvr_b(:,:) 354 ELSE 355 ub_b(:,:) = ub_b(:,:) * hur(:,:) 356 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 357 ENDIF 358 ! 218 359 ! !* baroclinic momentum trend (remove the vertical mean trend) 219 DO jk = 1, jpkm1 ! -------------------------- 360 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 220 361 DO jj = 2, jpjm1 221 362 DO ji = fs_2, fs_jpim1 ! vector opt. 222 ua(ji,jj,jk) = ua(ji,jj,jk) - zu a(ji,jj) * hur(ji,jj)223 va(ji,jj,jk) = va(ji,jj,jk) - zv a(ji,jj) * hvr(ji,jj)363 ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 364 va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 224 365 END DO 225 366 END DO 226 367 END DO 227 228 ! !* barotropic Coriolis trends * H (vorticity scheme dependent) 229 ! ! ---------------------------==== 230 zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport 231 zwy(:,:) = zvn(:,:) * e1v(:,:) 368 ! !* barotropic Coriolis trends (vorticity scheme dependent) 369 ! ! -------------------------------------------------------- 370 zwx(:,:) = un_b(:,:) * e2u(:,:) ! now transport 371 zwy(:,:) = vn_b(:,:) * e1v(:,:) 232 372 ! 233 373 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 239 379 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 240 380 ! energy conserving formulation for planetary vorticity term 241 z cu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )242 z cv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )243 END DO 244 END DO 245 ! 246 ELSEIF ( ln_dynvor_ens ) THEN 381 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 382 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 383 END DO 384 END DO 385 ! 386 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 247 387 DO jj = 2, jpjm1 248 388 DO ji = fs_2, fs_jpim1 ! vector opt. 249 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 250 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 251 zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 252 zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 253 END DO 254 END DO 255 ! 256 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 389 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 390 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 391 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 392 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 393 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 394 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 395 END DO 396 END DO 397 ! 398 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 257 399 DO jj = 2, jpjm1 258 400 DO ji = fs_2, fs_jpim1 ! vector opt. 259 zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 260 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 261 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) & 262 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 263 END DO 264 END DO 265 ! 266 ENDIF 267 401 zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 402 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 403 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 404 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 405 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 406 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 407 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 408 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 409 END DO 410 END DO 411 ! 412 ENDIF 413 ! 414 un_b (:,:) = un_b(:,:) * hur(:,:) ! Revert now transport to barotropic velocities 415 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 268 416 ! !* Right-Hand-Side of the barotropic momentum equation 269 417 ! ! ---------------------------------------------------- 270 IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient418 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 271 419 DO jj = 2, jpjm1 272 420 DO ji = fs_2, fs_jpim1 ! vector opt. 273 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 274 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 275 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 276 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 277 END DO 278 END DO 279 ENDIF 280 281 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 421 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 422 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 423 END DO 424 END DO 425 ENDIF 426 427 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 282 428 DO ji = fs_2, fs_jpim1 283 zu a(ji,jj) = zua(ji,jj) - zcu(ji,jj)284 zv a(ji,jj) = zva(ji,jj) - zcv(ji,jj)429 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1) 430 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1) 285 431 END DO 286 END DO 287 288 289 ! ! Remove barotropic contribution of bottom friction 290 ! ! from the barotropic transport trend 291 zcoef = -1._wp * z1_2dt_b 292 293 IF(ln_bfrimp) THEN 294 ! ! Remove the bottom stress trend from 3-D sea surface level gradient 295 ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction 296 DO jj = 2, jpjm1 297 DO ji = fs_2, fs_jpim1 298 ikbu = mbku(ji,jj) 299 ikbv = mbkv(ji,jj) 300 ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 301 va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 302 303 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 304 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 305 END DO 306 END DO 307 308 ELSE 309 310 # if defined key_vectopt_loop 311 DO jj = 1, 1 312 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 313 # else 314 DO jj = 2, jpjm1 315 DO ji = 2, jpim1 316 # endif 317 ! Apply stability criteria for bottom friction 318 !RBbug for vvl and external mode we may need to use varying fse3 319 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 320 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 321 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 322 END DO 323 END DO 324 325 IF( lk_vvl ) THEN 326 DO jj = 2, jpjm1 327 DO ji = fs_2, fs_jpim1 ! vector opt. 328 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 329 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 330 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 331 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 332 END DO 333 END DO 334 ELSE 335 DO jj = 2, jpjm1 336 DO ji = fs_2, fs_jpim1 ! vector opt. 337 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 338 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 339 END DO 340 END DO 341 ENDIF 342 END IF ! end (ln_bfrimp) 343 344 345 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 346 ! ! -------------------------- 347 zua(:,:) = zua(:,:) * hur(:,:) 348 zva(:,:) = zva(:,:) * hvr(:,:) 349 ! 350 IF( lk_vvl ) THEN 351 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 352 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 353 ELSE 354 ub_b(:,:) = ub_b(:,:) * hur(:,:) 355 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 356 ENDIF 357 432 END DO 433 ! 434 ! ! Add bottom stress contribution from baroclinic velocities: 435 IF (ln_bt_fw) THEN 436 DO jj = 2, jpjm1 437 DO ji = fs_2, fs_jpim1 ! vector opt. 438 ikbu = mbku(ji,jj) 439 ikbv = mbkv(ji,jj) 440 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 441 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 442 END DO 443 END DO 444 ELSE 445 DO jj = 2, jpjm1 446 DO ji = fs_2, fs_jpim1 ! vector opt. 447 ikbu = mbku(ji,jj) 448 ikbv = mbkv(ji,jj) 449 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 450 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 451 END DO 452 END DO 453 ENDIF 454 ! 455 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 456 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 457 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 458 ! 459 IF (ln_bt_fw) THEN ! Add wind forcing 460 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:) 461 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:) 462 ELSE 463 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:) 464 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:) 465 ENDIF 466 ! 467 IF ( ln_apr_dyn ) THEN ! Add atm pressure forcing 468 IF (ln_bt_fw) THEN 469 DO jj = 2, jpjm1 470 DO ji = fs_2, fs_jpim1 ! vector opt. 471 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj) 472 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj) 473 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 474 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 475 END DO 476 END DO 477 ELSE 478 DO jj = 2, jpjm1 479 DO ji = fs_2, fs_jpim1 ! vector opt. 480 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 481 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj) 482 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 483 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj) 484 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 485 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg 486 END DO 487 END DO 488 ENDIF 489 ENDIF 490 ! !* Right-Hand-Side of the barotropic ssh equation 491 ! ! ----------------------------------------------- 492 ! ! Surface net water flux and rivers 493 IF (ln_bt_fw) THEN 494 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) ) 495 ELSE 496 zssh_frc(:,:) = zraur * z1_2 * (emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:)) 497 ENDIF 498 #if defined key_asminc 499 ! ! Include the IAU weighted SSH increment 500 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 501 zssh_frc(:,:) = zssh_frc(:,:) + ssh_iau(:,:) 502 ENDIF 503 #endif 504 ! 358 505 ! ----------------------------------------------------------------------- 359 ! Phase 2 : Integration of the barotropic equations with time splitting506 ! Phase 2 : Integration of the barotropic equations 360 507 ! ----------------------------------------------------------------------- 361 508 ! 362 509 ! ! ==================== ! 363 510 ! ! Initialisations ! 511 ! ! ==================== ! 512 ! Initialize barotropic variables: 513 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 514 sshn_e (:,:) = sshn (:,:) 515 zun_e (:,:) = un_b (:,:) 516 zvn_e (:,:) = vn_b (:,:) 517 ELSE ! CENTERED integration: start from BEFORE fields 518 sshn_e (:,:) = sshb (:,:) 519 zun_e (:,:) = ub_b (:,:) 520 zvn_e (:,:) = vb_b (:,:) 521 ENDIF 522 ! 523 ! Initialize depths: 524 IF ( lk_vvl.AND.(.NOT.ln_bt_fw) ) THEN 525 hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 526 hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 527 hur_e (:,:) = zhur_b(:,:) 528 hvr_e (:,:) = zhvr_b(:,:) 529 ELSE 530 hu_e (:,:) = hu (:,:) 531 hv_e (:,:) = hv (:,:) 532 hur_e (:,:) = hur (:,:) 533 hvr_e (:,:) = hvr (:,:) 534 ENDIF 535 ! 536 IF (.NOT.lk_vvl) THEN ! Depths at jn+0.5: 537 zhup2_e (:,:) = hu(:,:) 538 zhvp2_e (:,:) = hv(:,:) 539 ENDIF 540 ! 541 ! Initialize sums: 542 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form) 543 va_b (:,:) = 0._wp 544 ssha (:,:) = 0._wp ! Sum for after averaged sea level 545 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop 546 zv_sum(:,:) = 0._wp 364 547 ! ! ==================== ! 365 icycle = 2 * nn_baro ! Number of barotropic sub time-step 366 367 ! ! Start from NOW field 368 hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points 369 hv_e (:,:) = hv (:,:) 370 hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points 371 hvr_e (:,:) = hvr (:,:) 372 !RBbug zsshb_e(:,:) = sshn (:,:) 373 zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now) 374 sshn_e (:,:) = sshn (:,:) 375 376 zun_e (:,:) = un_b (:,:) ! barotropic velocity (external) 377 zvn_e (:,:) = vn_b (:,:) 378 zub_e (:,:) = un_b (:,:) 379 zvb_e (:,:) = vn_b (:,:) 380 381 zu_sum (:,:) = un_b (:,:) ! summation 382 zv_sum (:,:) = vn_b (:,:) 383 zssh_sum(:,:) = sshn (:,:) 384 385 #if defined key_obc 386 ! set ssh corrections to 0 387 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 388 IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp 389 IF( lp_obc_west ) sshfow_b(:,:) = 0._wp 390 IF( lp_obc_south ) sshfos_b(:,:) = 0._wp 391 IF( lp_obc_north ) sshfon_b(:,:) = 0._wp 392 #endif 393 394 ! ! ==================== ! 395 DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1) 548 DO jn = 1, icycle ! sub-time-step loop ! 396 549 ! ! ==================== ! 397 z2dt_e = 2. * ( rdt / nn_baro )398 IF( jn == 1 ) z2dt_e = rdt / nn_baro399 400 550 ! !* Update the forcing (BDY and tides) 401 551 ! ! ------------------ 402 IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) 403 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 404 IF ( ln_tide_pot .AND. lk_tide) CALL upd_tide( kt, jn ) 405 406 ! !* after ssh_e 552 ! Update only tidal forcing at open boundaries 553 #if defined key_tide 554 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1) ) 555 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset ) 556 #endif 557 ! 558 ! Set extrapolation coefficients for predictor step: 559 IF ((jn<3).AND.ll_init) THEN ! Forward 560 za1 = 1._wp 561 za2 = 0._wp 562 za3 = 0._wp 563 ELSE ! AB3-AM4 Coefficients: bet=0.281105 564 za1 = 1.781105_wp ! za1 = 3/2 + bet 565 za2 = -1.06221_wp ! za2 = -(1/2 + 2*bet) 566 za3 = 0.281105_wp ! za3 = bet 567 ENDIF 568 569 ! Extrapolate barotropic velocities at step jit+0.5: 570 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 571 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 572 573 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 574 ! ! ------------------ 575 ! Extrapolate Sea Level at step jit+0.5: 576 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 577 ! 578 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 579 DO ji = 2, fs_jpim1 ! Vector opt. 580 zwx(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 581 & * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) & 582 & + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 583 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 584 & * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) & 585 & + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 586 END DO 587 END DO 588 CALL lbc_lnk( zwx, 'U', 1._wp ) ; CALL lbc_lnk( zwy, 'V', 1._wp ) 589 ! 590 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 591 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 592 ENDIF 593 ! !* after ssh 407 594 ! ! ----------- 408 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 595 ! One should enforce volume conservation at open boundaries here 596 ! considering fluxes below: 597 ! 598 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 599 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 600 DO jj = 2, jpjm1 409 601 DO ji = fs_2, fs_jpim1 ! vector opt. 410 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & 411 & - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) & 412 & + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) & 413 & - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 414 END DO 415 END DO 416 ! 417 #if defined key_obc 418 ! ! OBC : zhdiv must be zero behind the open boundary 419 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 420 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east 421 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west 422 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north 423 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south 602 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 603 & + zwy(ji,jj) - zwy(ji,jj-1) & 604 & ) / ( e1t(ji,jj) * e2t(ji,jj) ) 605 END DO 606 END DO 607 ! 608 ! Sum over sub-time-steps to compute advective velocities 609 za2 = wgtbtp2(jn) 610 zu_sum (:,:) = zu_sum (:,:) + za2 * ua_e (:,:) * zhup2_e (:,:) 611 zv_sum (:,:) = zv_sum (:,:) + za2 * va_e (:,:) * zhvp2_e (:,:) 612 ! 613 ! Set next sea level: 614 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 615 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 616 617 #if defined key_bdy 618 ! Duplicate sea level across open boundaries (this is only cosmetic if lk_vvl=.false.) 619 IF (lk_bdy) CALL bdy_ssh( ssha_e ) 424 620 #endif 425 #if defined key_ bdy426 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask621 #if defined key_agrif 622 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 427 623 #endif 428 ! 429 DO jj = 2, jpjm1 ! leap-frog on ssh_e 430 DO ji = fs_2, fs_jpim1 ! vector opt. 431 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 432 END DO 433 END DO 434 435 ! !* after barotropic velocities (vorticity scheme dependent) 436 ! ! --------------------------- 437 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 438 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 624 ! 625 ! Sea Surface Height at u-,v-points (vvl case only) 626 IF ( lk_vvl ) THEN 627 DO jj = 2, jpjm1 628 DO ji = 2, jpim1 ! NO Vector Opt. 629 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 630 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha_e(ji ,jj) & 631 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha_e(ji+1,jj) ) 632 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 633 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha_e(ji,jj ) & 634 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha_e(ji,jj+1) ) 635 END DO 636 END DO 637 CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) 638 ENDIF 639 ! 640 ! Half-step back interpolation of SSH for surface pressure computation: 641 !---------------------------------------------------------------------- 642 IF ((jn==1).AND.ll_init) THEN 643 za0=1._wp ! Forward-backward 644 za1=0._wp 645 za2=0._wp 646 za3=0._wp 647 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 648 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 649 za1=-0.1666666666666_wp ! za1 = gam 650 za2= 0.0833333333333_wp ! za2 = eps 651 za3= 0._wp 652 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 653 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 654 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 655 za2=0.088_wp ! za2 = gam 656 za3=0.013_wp ! za3 = eps 657 ENDIF 658 659 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 660 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 661 662 ! 663 ! Compute associated depths at U and V points: 664 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN 665 ! 666 DO jj = 2, jpjm1 667 DO ji = 2, jpim1 668 zx1 = z1_2 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 669 & * ( e1t(ji ,jj) * e2t(ji ,jj) * zsshp2_e(ji ,jj) & 670 & + e1t(ji+1,jj) * e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 671 zy1 = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 672 & * ( e1t(ji,jj ) * e2t(ji,jj ) * zsshp2_e(ji,jj ) & 673 & + e1t(ji,jj+1) * e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 674 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 675 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 676 END DO 677 END DO 678 ENDIF 679 ! 680 ! Add Coriolis trend: 681 ! zwz array below or triads normally depend on sea level with key_vvl and should be updated 682 ! at each time step. We however keep them constant here for optimization. 683 ! Recall that zwx and zwy arrays hold fluxes at this stage: 684 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 685 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 439 686 ! 440 687 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 441 688 DO jj = 2, jpjm1 442 689 DO ji = fs_2, fs_jpim1 ! vector opt. 443 ! surface pressure gradient444 IF( lk_vvl) THEN445 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &446 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)447 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &448 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)449 ELSE450 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)451 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)452 ENDIF453 ! add tidal astronomical forcing454 IF ( ln_tide_pot .AND. lk_tide ) THEN455 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj)456 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj)457 ENDIF458 ! energy conserving formulation for planetary vorticity term459 690 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) 460 691 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 461 692 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 462 693 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 463 zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 464 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 465 ! after velocities with implicit bottom friction 466 467 IF( ln_bfrimp ) THEN ! implicit bottom friction 468 ! A new method to implement the implicit bottom friction. 469 ! H. Liu 470 ! Sept 2011 471 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 472 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 473 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 474 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 475 ! 476 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 477 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 478 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 479 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 480 ! 481 ELSE 482 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) & 483 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 484 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) & 485 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 486 ENDIF 694 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 695 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 487 696 END DO 488 697 END DO … … 491 700 DO jj = 2, jpjm1 492 701 DO ji = fs_2, fs_jpim1 ! vector opt. 493 ! surface pressure gradient 494 IF( lk_vvl) THEN 495 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 496 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 497 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 498 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 499 ELSE 500 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 501 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 502 ENDIF 503 ! add tidal astronomical forcing 504 IF ( ln_tide_pot .AND. lk_tide ) THEN 505 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 506 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 507 ENDIF 508 ! enstrophy conserving formulation for planetary vorticity term 509 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 510 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 511 zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 512 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 513 ! after velocities with implicit bottom friction 514 IF( ln_bfrimp ) THEN 515 ! A new method to implement the implicit bottom friction. 516 ! H. Liu 517 ! Sept 2011 518 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 519 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 520 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 521 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 522 ! 523 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 524 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 525 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 526 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 527 ! 528 ELSE 529 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) & 530 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 531 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) & 532 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 533 ENDIF 702 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 704 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 706 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 534 708 END DO 535 709 END DO … … 538 712 DO jj = 2, jpjm1 539 713 DO ji = fs_2, fs_jpim1 ! vector opt. 540 ! surface pressure gradient 541 IF( lk_vvl) THEN 542 zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 543 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 544 zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 545 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 546 ELSE 547 zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 548 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 549 ENDIF 550 ! add tidal astronomical forcing 551 IF ( ln_tide_pot .AND. lk_tide ) THEN 552 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 553 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 554 ENDIF 555 ! energy/enstrophy conserving formulation for planetary vorticity term 556 zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 557 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 558 zv_cor = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 559 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 560 ! after velocities with implicit bottom friction 561 IF( ln_bfrimp ) THEN 562 ! A new method to implement the implicit bottom friction. 563 ! H. Liu 564 ! Sept 2011 565 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 566 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 567 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 568 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 569 ! 570 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 571 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 572 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 573 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 574 ! 575 ELSE 576 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) & 577 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 578 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) & 579 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 580 ENDIF 714 zu_trd(ji,jj) = + z1_12 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 715 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 716 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 717 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 718 zv_trd(ji,jj) = - z1_12 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 719 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 720 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 721 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 581 722 END DO 582 723 END DO 583 724 ! 584 725 ENDIF 585 ! !* domain lateral boundary 586 ! ! ----------------------- 587 588 ! OBC open boundaries 589 IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 590 591 ! BDY open boundaries 592 #if defined key_bdy 593 pssh => sshn_e 594 phur => hur_e 595 phvr => hvr_e 596 pu2d => ua_e 597 pv2d => va_e 598 599 IF( lk_bdy ) CALL bdy_dyn2d( kt ) 600 #endif 601 602 ! 603 CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries 604 CALL lbc_lnk( va_e , 'V', -1. ) 605 CALL lbc_lnk( ssha_e, 'T', 1. ) 606 607 zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps 608 zv_sum (:,:) = zv_sum (:,:) + va_e (:,:) 609 zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 610 611 ! !* Time filter and swap 612 ! ! -------------------- 613 IF( jn == 1 ) THEN ! Swap only (1st Euler time step) 614 zsshb_e(:,:) = sshn_e(:,:) 615 zub_e (:,:) = zun_e (:,:) 616 zvb_e (:,:) = zvn_e (:,:) 617 sshn_e (:,:) = ssha_e(:,:) 618 zun_e (:,:) = ua_e (:,:) 619 zvn_e (:,:) = va_e (:,:) 620 ELSE ! Swap + Filter 621 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 622 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) 623 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) 624 sshn_e (:,:) = ssha_e(:,:) 625 zun_e (:,:) = ua_e (:,:) 626 zvn_e (:,:) = va_e (:,:) 627 ENDIF 628 629 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 630 ! ! ------------------ 631 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 632 DO ji = 1, fs_jpim1 ! Vector opt. 633 zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 634 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 635 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 636 zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 637 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 638 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 639 END DO 640 END DO 641 CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions 642 CALL lbc_lnk( zsshvn_e, 'V', 1. ) 643 ! 644 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 645 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 726 ! 727 ! Add tidal astronomical forcing if defined 728 IF ( lk_tide.AND.ln_tide_pot ) THEN 729 DO jj = 2, jpjm1 730 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 732 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 733 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 734 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 735 END DO 736 END DO 737 ENDIF 738 ! 739 ! Add bottom stresses: 740 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 741 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 742 ! 743 ! Surface pressure trend: 744 DO jj = 2, jpjm1 745 DO ji = fs_2, fs_jpim1 ! vector opt. 746 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 749 zwx(ji,jj) = zu_spg 750 zwy(ji,jj) = zv_spg 751 END DO 752 END DO 753 ! 754 ! Set next velocities: 755 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN ! Vector form 756 DO jj = 2, jpjm1 757 DO ji = fs_2, fs_jpim1 ! vector opt. 758 ua_e(ji,jj) = ( zun_e(ji,jj) & 759 & + rdtbt * ( zwx(ji,jj) & 760 & + zu_trd(ji,jj) & 761 & + zu_frc(ji,jj) ) & 762 & ) * umask(ji,jj,1) 763 764 va_e(ji,jj) = ( zvn_e(ji,jj) & 765 & + rdtbt * ( zwy(ji,jj) & 766 & + zv_trd(ji,jj) & 767 & + zv_frc(ji,jj) ) & 768 & ) * vmask(ji,jj,1) 769 END DO 770 END DO 771 772 ELSE ! Flux form 773 DO jj = 2, jpjm1 774 DO ji = fs_2, fs_jpim1 ! vector opt. 775 776 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 777 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 779 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 780 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 782 & + hu(ji,jj) * zu_frc(ji,jj) ) & 783 & ) * zhura 784 785 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) & 786 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 787 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 788 & + hv(ji,jj) * zv_frc(ji,jj) ) & 789 & ) * zhvra 790 END DO 791 END DO 792 ENDIF 793 ! 794 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 795 ! ! ---------------------------------------------- 796 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 797 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 646 798 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 647 799 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 648 800 ! 649 801 ENDIF 802 ! !* domain lateral boundary 803 ! ! ----------------------- 804 ! 805 CALL lbc_lnk( ua_e , 'U', -1._wp ) ! local domain boundaries 806 CALL lbc_lnk( va_e , 'V', -1._wp ) 807 808 #if defined key_bdy 809 810 pssh => ssha_e 811 phur => hur_e 812 phvr => hvr_e 813 pua2d => ua_e 814 pva2d => va_e 815 pub2d => zun_e 816 pvb2d => zvn_e 817 818 IF( lk_bdy ) CALL bdy_dyn2d( kt ) ! open boundaries 819 #endif 820 #if defined key_agrif 821 IF( .NOT.Agrif_Root() ) CALL agrif_dyn_ts( kt, jn ) ! Agrif 822 #endif 823 ! !* Swap 824 ! ! ---- 825 ubb_e (:,:) = ub_e (:,:) 826 ub_e (:,:) = zun_e (:,:) 827 zun_e (:,:) = ua_e (:,:) 828 ! 829 vbb_e (:,:) = vb_e (:,:) 830 vb_e (:,:) = zvn_e (:,:) 831 zvn_e (:,:) = va_e (:,:) 832 ! 833 sshbb_e(:,:) = sshb_e(:,:) 834 sshb_e (:,:) = sshn_e(:,:) 835 sshn_e (:,:) = ssha_e(:,:) 836 837 ! !* Sum over whole bt loop 838 ! ! ---------------------- 839 za1 = wgtbtp1(jn) 840 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities 841 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 842 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 843 ELSE ! Sum transports 844 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 845 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) 846 ENDIF 847 ! ! Sum sea level 848 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:) 650 849 ! ! ==================== ! 651 850 END DO ! end loop ! 652 851 ! ! ==================== ! 653 654 #if defined key_obc655 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????656 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)657 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)658 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)659 #endif660 661 852 ! ----------------------------------------------------------------------------- 662 853 ! Phase 3. update the general trend with the barotropic trend 663 854 ! ----------------------------------------------------------------------------- 664 855 ! 665 ! !* Time average ==> after barotropic u, v, ssh 666 zcoef = 1._wp / ( 2 * nn_baro + 1 ) 667 zu_sum(:,:) = zcoef * zu_sum (:,:) 668 zv_sum(:,:) = zcoef * zv_sum (:,:) 669 ! 670 ! !* update the general momentum trend 671 DO jk=1,jpkm1 672 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 673 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 856 ! At this stage ssha holds a time averaged value 857 ! ! Sea Surface Height at u-,v- and f-points 858 IF( lk_vvl ) THEN ! (required only in key_vvl case) 859 DO jj = 1, jpjm1 860 DO ji = 1, jpim1 ! NO Vector Opt. 861 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) / ( e1u(ji ,jj) * e2u(ji ,jj) ) & 862 & * ( e1t(ji ,jj) * e2t(ji ,jj) * ssha(ji ,jj) & 863 & + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) ) 864 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) / ( e1v(ji,jj ) * e2v(ji,jj ) ) & 865 & * ( e1t(ji,jj ) * e2t(ji,jj ) * ssha(ji,jj ) & 866 & + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) ) 867 END DO 868 END DO 869 CALL lbc_lnk( zsshu_a, 'U', 1._wp ) ; CALL lbc_lnk( zsshv_a, 'V', 1._wp ) ! Boundary conditions 870 ENDIF 871 ! 872 ! Set advection velocity correction: 873 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 874 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 875 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 876 ELSE 877 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zu_sum(:,:)) * hur(:,:) 878 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zv_sum(:,:)) * hvr(:,:) 879 END IF 880 881 IF (ln_bt_fw) THEN ! Save integrated transport for next computation 882 ub2_b(:,:) = zu_sum(:,:) 883 vb2_b(:,:) = zv_sum(:,:) 884 ENDIF 885 ! 886 ! Update barotropic trend: 887 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN 888 DO jk=1,jpkm1 889 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b 890 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * z1_2dt_b 891 END DO 892 ELSE 893 hu_e (:,:) = umask(:,:,1) / ( zhur_b(:,:) + 1._wp - umask(:,:,1) ) 894 hv_e (:,:) = vmask(:,:,1) / ( zhvr_b(:,:) + 1._wp - vmask(:,:,1) ) 895 DO jk=1,jpkm1 896 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_e(:,:) ) * z1_2dt_b 897 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_e(:,:) ) * z1_2dt_b 898 END DO 899 ! Save barotropic velocities not transport: 900 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) ) 901 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) ) 902 ENDIF 903 ! 904 DO jk = 1, jpkm1 905 ! Correct velocities: 906 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) )*umask(:,:,jk) 907 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) )*vmask(:,:,jk) 908 ! 674 909 END DO 675 un_b (:,:) = zu_sum(:,:)676 vn_b (:,:) = zv_sum(:,:)677 sshn_b(:,:) = zcoef * zssh_sum(:,:)678 910 ! 679 911 ! !* write time-spliting arrays in the restart 680 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 681 ! 682 CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 683 CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 684 CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 912 IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' ) 913 ! 914 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 915 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 916 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 917 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 918 CALL wrk_dealloc( jpi, jpj, zhur_b, zhvr_b ) 919 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 920 CALL wrk_dealloc( jpi, jpj, zht, zhf ) 685 921 ! 686 922 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') … … 688 924 END SUBROUTINE dyn_spg_ts 689 925 926 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) 927 !!--------------------------------------------------------------------- 928 !! *** ROUTINE ts_wgt *** 929 !! 930 !! ** Purpose : Set time-splitting weights for temporal averaging (or not) 931 !!---------------------------------------------------------------------- 932 LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true. 933 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 934 INTEGER, INTENT(inout) :: jpit ! cycle length 935 REAL(wp), DIMENSION(3*nn_baro), INTENT(inout) :: zwgt1, & ! Primary weights 936 zwgt2 ! Secondary weights 937 938 INTEGER :: jic, jn, ji ! temporary integers 939 REAL(wp) :: za1, za2 940 !!---------------------------------------------------------------------- 941 942 zwgt1(:) = 0._wp 943 zwgt2(:) = 0._wp 944 945 ! Set time index when averaged value is requested 946 IF (ll_fw) THEN 947 jic = nn_baro 948 ELSE 949 jic = 2 * nn_baro 950 ENDIF 951 952 ! Set primary weights: 953 IF (ll_av) THEN 954 ! Define simple boxcar window for primary weights 955 ! (width = nn_baro, centered around jic) 956 SELECT CASE ( nn_bt_flt ) 957 CASE( 0 ) ! No averaging 958 zwgt1(jic) = 1._wp 959 jpit = jic 960 961 CASE( 1 ) ! Boxcar, width = nn_baro 962 DO jn = 1, 3*nn_baro 963 za1 = ABS(float(jn-jic))/float(nn_baro) 964 IF (za1 < 0.5_wp) THEN 965 zwgt1(jn) = 1._wp 966 jpit = jn 967 ENDIF 968 ENDDO 969 970 CASE( 2 ) ! Boxcar, width = 2 * nn_baro 971 DO jn = 1, 3*nn_baro 972 za1 = ABS(float(jn-jic))/float(nn_baro) 973 IF (za1 < 1._wp) THEN 974 zwgt1(jn) = 1._wp 975 jpit = jn 976 ENDIF 977 ENDDO 978 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' ) 979 END SELECT 980 981 ELSE ! No time averaging 982 zwgt1(jic) = 1._wp 983 jpit = jic 984 ENDIF 985 986 ! Set secondary weights 987 DO jn = 1, jpit 988 DO ji = jn, jpit 989 zwgt2(jn) = zwgt2(jn) + zwgt1(ji) 990 END DO 991 END DO 992 993 ! Normalize weigths: 994 za1 = 1._wp / SUM(zwgt1(1:jpit)) 995 za2 = 1._wp / SUM(zwgt2(1:jpit)) 996 DO jn = 1, jpit 997 zwgt1(jn) = zwgt1(jn) * za1 998 zwgt2(jn) = zwgt2(jn) * za2 999 END DO 1000 ! 1001 END SUBROUTINE ts_wgt 690 1002 691 1003 SUBROUTINE ts_rst( kt, cdrw ) … … 698 1010 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 699 1011 ! 700 INTEGER :: ji, jk ! dummy loop indices701 1012 !!---------------------------------------------------------------------- 702 1013 ! 703 1014 IF( TRIM(cdrw) == 'READ' ) THEN 704 IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 705 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued 706 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 1015 CALL iom_get( numror, jpdom_autoglo, 'ub2_b' , ub2_b (:,:) ) 1016 CALL iom_get( numror, jpdom_autoglo, 'vb2_b' , vb2_b (:,:) ) 1017 IF( .NOT.ln_bt_av .AND. iom_varid( numror, 'sshbb_e', ldstop = .FALSE. ) > 0) THEN 1018 CALL iom_get( numror, jpdom_autoglo, 'sshbb_e' , sshbb_e(:,:) ) 1019 CALL iom_get( numror, jpdom_autoglo, 'ubb_e' , ubb_e(:,:) ) 1020 CALL iom_get( numror, jpdom_autoglo, 'vbb_e' , vbb_e(:,:) ) 1021 CALL iom_get( numror, jpdom_autoglo, 'sshb_e' , sshb_e(:,:) ) 1022 CALL iom_get( numror, jpdom_autoglo, 'ub_e' , ub_e(:,:) ) 1023 CALL iom_get( numror, jpdom_autoglo, 'vb_e' , vb_e(:,:) ) 707 1024 ELSE 708 un_b (:,:) = 0._wp 709 vn_b (:,:) = 0._wp 710 ! vertical sum 711 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 712 DO jk = 1, jpkm1 713 DO ji = 1, jpij 714 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 715 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 716 END DO 717 END DO 718 ELSE ! No vector opt. 719 DO jk = 1, jpkm1 720 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 721 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 722 END DO 723 ENDIF 724 un_b (:,:) = un_b(:,:) * hur(:,:) 725 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 726 ENDIF 727 728 ! Vertically integrated velocity (before) 729 IF (neuler/=0) THEN 730 ub_b (:,:) = 0._wp 731 vb_b (:,:) = 0._wp 732 733 ! vertical sum 734 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 735 DO jk = 1, jpkm1 736 DO ji = 1, jpij 737 ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 738 vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 739 END DO 740 END DO 741 ELSE ! No vector opt. 742 DO jk = 1, jpkm1 743 ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 744 vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 745 END DO 746 ENDIF 747 748 IF( lk_vvl ) THEN 749 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 750 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 751 ELSE 752 ub_b(:,:) = ub_b(:,:) * hur(:,:) 753 vb_b(:,:) = vb_b(:,:) * hvr(:,:) 754 ENDIF 755 ELSE ! neuler==0 756 ub_b (:,:) = un_b (:,:) 757 vb_b (:,:) = vn_b (:,:) 758 ENDIF 759 760 IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 761 CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered ssh 762 ELSE 763 sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value 764 ENDIF 1025 sshbb_e = sshn_b ! ACC GUESS WORK 1026 ubb_e = ub_b 1027 vbb_e = vb_b 1028 sshb_e = sshn_b 1029 ub_e = ub_b 1030 vb_e = vb_b 1031 ENDIF 1032 ! 765 1033 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 766 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh 767 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop 768 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) ! 1034 CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) ) 1035 CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) ) 1036 ! 1037 IF (.NOT.ln_bt_av) THEN 1038 CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) ) 1039 CALL iom_rstput( kt, nitrst, numrow, 'ubb_e' , ubb_e(:,:) ) 1040 CALL iom_rstput( kt, nitrst, numrow, 'vbb_e' , vbb_e(:,:) ) 1041 CALL iom_rstput( kt, nitrst, numrow, 'sshb_e' , sshb_e(:,:) ) 1042 CALL iom_rstput( kt, nitrst, numrow, 'ub_e' , ub_e(:,:) ) 1043 CALL iom_rstput( kt, nitrst, numrow, 'vb_e' , vb_e(:,:) ) 1044 ENDIF 769 1045 ENDIF 770 1046 ! 771 1047 END SUBROUTINE ts_rst 772 1048 1049 SUBROUTINE dyn_spg_ts_init( kt ) 1050 !!--------------------------------------------------------------------- 1051 !! *** ROUTINE dyn_spg_ts_init *** 1052 !! 1053 !! ** Purpose : Set time splitting options 1054 !!---------------------------------------------------------------------- 1055 INTEGER , INTENT(in) :: kt ! ocean time-step 1056 ! 1057 INTEGER :: ji ,jj, jk 1058 REAL(wp) :: zxr2, zyr2, zcmax 1059 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zht 1060 !! 1061 ! NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 1062 ! & nn_baro, rn_bt_cmax, nn_bt_flt 1063 !!---------------------------------------------------------------------- 1064 ! REWIND( numnam ) !* Namelist namsplit: split-explicit free surface 1065 ! READ ( numnam, namsplit ) 1066 ! ! Max courant number for ext. grav. waves 1067 ! 1068 CALL wrk_alloc( jpi, jpj, zcu, zht ) 1069 ! 1070 ! JC: Simplification needed below: define ht_0 even when volume is fixed 1071 IF (lk_vvl) THEN 1072 zht(:,:) = ht_0(:,:) * tmask(:,:,1) 1073 ELSE 1074 zht(:,:) = 0. 1075 DO jk = 1, jpkm1 1076 zht(:,:) = zht(:,:) + fse3t_n(:,:,jk) * tmask(:,:,jk) 1077 END DO 1078 ENDIF 1079 1080 DO jj = 1, jpj 1081 DO ji =1, jpi 1082 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1083 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1084 zcu(ji,jj) = sqrt(grav*zht(ji,jj)*(zxr2 + zyr2) ) 1085 END DO 1086 END DO 1087 1088 zcmax = MAXVAL(zcu(:,:)) 1089 IF( lk_mpp ) CALL mpp_max( zcmax ) 1090 1091 ! Estimate number of iterations to satisfy a max courant number=0.8 1092 IF (ln_bt_nn_auto) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1093 1094 rdtbt = rdt / FLOAT(nn_baro) 1095 zcmax = zcmax * rdtbt 1096 ! Print results 1097 IF(lwp) WRITE(numout,*) 1098 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1099 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1100 IF( ln_bt_nn_auto ) THEN 1101 IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.true. Automatically set nn_baro ' 1102 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1103 ELSE 1104 IF(lwp) WRITE(numout,*) ' ln_ts_nn_auto=.false.: Use nn_baro in namelist ' 1105 ENDIF 1106 IF(lwp) WRITE(numout,*) ' nn_baro = ', nn_baro 1107 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rdtbt 1108 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1109 1110 IF(ln_bt_av) THEN 1111 IF(lwp) WRITE(numout,*) ' ln_bt_av=.true. => Time averaging over nn_baro time steps is on ' 1112 ELSE 1113 IF(lwp) WRITE(numout,*) ' ln_bt_av=.false. => No time averaging of barotropic variables ' 1114 ENDIF 1115 ! 1116 ! 1117 IF(ln_bt_fw) THEN 1118 IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables ' 1119 ELSE 1120 IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centered integration of barotropic variables ' 1121 ENDIF 1122 ! 1123 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1124 SELECT CASE ( nn_bt_flt ) 1125 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1126 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1127 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1128 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1129 END SELECT 1130 ! 1131 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN 1132 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1133 ENDIF 1134 IF ( zcmax>0.9_wp ) THEN 1135 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1136 ENDIF 1137 ! 1138 CALL wrk_dealloc( jpi, jpj, zcu, zht ) 1139 ! 1140 END SUBROUTINE dyn_spg_ts_init 1141 773 1142 #else 774 !!---------------------------------------------------------------------- 775 !! Default case : Empty module No standart free surface cst volume 776 !!---------------------------------------------------------------------- 1143 !!--------------------------------------------------------------------------- 1144 !! Default case : Empty module No standard free surface constant volume 1145 !!--------------------------------------------------------------------------- 1146 1147 USE par_kind 1148 LOGICAL, PUBLIC, PARAMETER :: ln_bt_fw=.FALSE. ! Forward integration of barotropic sub-stepping 777 1149 CONTAINS 778 1150 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function … … 787 1159 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 788 1160 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw 789 END SUBROUTINE ts_rst 1161 END SUBROUTINE ts_rst 1162 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine 1163 INTEGER , INTENT(in) :: kt ! ocean time-step 1164 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt 1165 END SUBROUTINE dyn_spg_ts_init 790 1166 #endif 791 1167
Note: See TracChangeset
for help on using the changeset viewer.
