Changeset 1502 for trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2009-07-20T17:20:23+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (17 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r1438 r1502 1 1 MODULE dynspg_ts 2 2 !!====================================================================== 3 !! History : 1.0 ! 04-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 05-11 (V. Garnier, G. Madec) optimization 5 !! - ! 06-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 07-07 (D. Storkey) calls to BDY routines 7 !! - ! 08-01 (R. Benshila) change averaging method 3 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 2005-11 (V. Garnier, G. Madec) optimization 5 !! - ! 2006-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines 7 !! - ! 2008-01 (R. Benshila) change averaging method 8 !! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation 8 9 !!--------------------------------------------------------------------- 9 10 #if defined key_dynspg_ts || defined key_esopa 10 11 !!---------------------------------------------------------------------- 11 12 !! 'key_dynspg_ts' free surface cst volume with time splitting 12 !!----------------------------------------------------------------------13 13 !!---------------------------------------------------------------------- 14 14 !! dyn_spg_ts : compute surface pressure gradient trend using a time- … … 45 45 PUBLIC ts_rst ! routine called by istate.F90 46 46 47 47 48 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter 48 49 REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme) 50 51 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: un_b, vn_b ! averaged velocity 49 52 50 53 !! * Substitutions … … 66 69 !! gradient in case of free surface formulation with time-splitting. 67 70 !! Add it to the general trend of momentum equation. 68 !! Compute the free surface.69 71 !! 70 72 !! ** Method : Free surface formulation with time-splitting … … 75 77 !! the barotropic integration. 76 78 !! -2- Barotropic loop : updates of sea surface height (ssha_e) and 77 !! barotropic transports(ua_e and va_e) through barotropic79 !! barotropic velocity (ua_e and va_e) through barotropic 78 80 !! momentum and continuity integration. Barotropic former 79 81 !! variables are time averaging over the full barotropic cycle 80 !! (= 2 * baroclinic time step) and saved in z sshX_b, zuX_b82 !! (= 2 * baroclinic time step) and saved in zuX_b 81 83 !! and zvX_b (X specifying after, now or before). 82 84 !! -3- The new general trend becomes : 83 !! ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H85 !! ua = ua - sum_k(ua)/H + ( ua_e - sum_k(ub) ) 84 86 !! 85 87 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 87 89 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 88 90 !!--------------------------------------------------------------------- 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 90 !! 91 INTEGER :: ji, jj, jk, jit ! dummy loop indices 92 INTEGER :: icycle ! temporary scalar 93 REAL(wp) :: & 94 zraur, zcoef, z2dt_e, z2dt_b, zfac25, & ! temporary scalars 95 zfact1, zspgu, zcubt, zx1, zy1, & ! " " 96 zfact2, zspgv, zcvbt, zx2, zy2 ! " " 97 REAL(wp), DIMENSION(jpi,jpj) :: & 98 zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays 99 zua, zva, zub, zvb, & ! " " 100 zua_e, zva_e, zssha_e, & ! " " 101 zub_e, zvb_e, zsshb_e, & ! " " 102 zu_sum, zv_sum, zssh_sum 103 !! Variable volume 104 REAL(wp), DIMENSION(jpi,jpj) :: & 105 zspgu_1, zspgv_1, zsshun_e, zsshvn_e ! 2D workspace 91 INTEGER, INTENT(in) :: kt ! ocean time-step index 92 !! 93 INTEGER :: ji, jj, jk, jn ! dummy loop indices 94 INTEGER :: icycle ! temporary scalar 95 REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b, & ! temporary scalars 96 z1_8, zspgu, zcubt, zx1, zy1, & ! " " 97 z1_4, zspgv, zcvbt, zx2, zy2 ! " " 98 REAL(wp), DIMENSION(jpi,jpj) :: zhdiv, zsshb_e 99 !! 100 REAL(wp), DIMENSION(jpi,jpj) :: zsshun_e, zsshvn_e ! 2D workspace 101 !! 102 REAL(wp), DIMENSION(jpi,jpj) :: zcu, zwx, zua, zun, zub ! 2D workspace 103 REAL(wp), DIMENSION(jpi,jpj) :: zcv, zwy, zva, zvn, zvb ! - - 104 REAL(wp), DIMENSION(jpi,jpj) :: zun_e, zub_e, zu_sum ! 2D workspace 105 REAL(wp), DIMENSION(jpi,jpj) :: zvn_e, zvb_e, zv_sum ! - - 106 REAL(wp), DIMENSION(jpi,jpj) :: zssh_sum ! - - 106 107 !!---------------------------------------------------------------------- 107 108 108 ! Arrays initialization 109 ! --------------------- 110 zhdiv(:,:) = 0.e0 111 112 113 IF( kt == nit000 ) THEN 109 IF( kt == nit000 ) THEN !* initialisation 114 110 ! 115 111 IF(lwp) WRITE(numout,*) … … 117 113 IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting' 118 114 IF(lwp) WRITE(numout,*) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2*nn_baro 119 115 ! 120 116 CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 122 123 zssha_e(:,:) = sshn(:,:) 124 zua_e (:,:) = un_e(:,:) 125 zva_e (:,:) = vn_e(:,:) 126 hu_e (:,:) = hu (:,:) 127 hv_e (:,:) = hv (:,:) 117 ! ! un_b, vn_b 118 ! 119 ua_e (:,:) = un_b (:,:) 120 va_e (:,:) = vn_b (:,:) 121 hu_e (:,:) = hu (:,:) 122 hv_e (:,:) = hv (:,:) 123 hur_e (:,:) = hur (:,:) 124 hvr_e (:,:) = hvr (:,:) 128 125 IF( ln_dynvor_een ) THEN 129 126 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 … … 140 137 ENDIF 141 138 142 ! Substract the surface pressure gradient from the momentum 143 ! --------------------------------------------------------- 144 IF( lk_vvl ) THEN ! Variable volume 145 146 ! 0. Local initialization 147 ! ----------------------- 148 zspgu_1(:,:) = 0.e0 149 zspgv_1(:,:) = 0.e0 150 151 ! 1. Surface pressure gradient (now) 152 ! ---------------------------- 153 DO jj = 2, jpjm1 154 DO ji = fs_2, fs_jpim1 ! vector opt. 155 zspgu_1(ji,jj) = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 156 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e1u(ji,jj) 157 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 158 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 159 END DO 160 END DO 161 162 ! 2. Substract the surface pressure trend from the general trend 163 ! ------------------------------------------------------ 164 DO jk = 1, jpkm1 165 DO jj = 2, jpjm1 166 DO ji = fs_2, fs_jpim1 ! vector opt. 167 ua(ji,jj,jk) = ua(ji,jj,jk) - zspgu_1(ji,jj) 168 va(ji,jj,jk) = va(ji,jj,jk) - zspgv_1(ji,jj) 169 END DO 170 END DO 171 END DO 172 173 ENDIF 174 175 ! Local constant initialization 176 ! -------------------------------- 139 ! !* Local constant initialization 177 140 z2dt_b = 2.0 * rdt ! baroclinic time step 178 IF ( neuler == 0 .AND. kt == nit000 ) z2dt_b = rdt 179 zfact1 = 0.5 * 0.25 ! coefficient for vorticity estimates 180 zfact2 = 0.5 * 0.5 141 z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates 142 z1_4 = 0.5 * 0.5 181 143 zraur = 1. / rauw ! 1 / volumic mass of pure water 144 ! 145 zhdiv(:,:) = 0.e0 ! barotropic divergence 182 146 183 147 ! ----------------------------------------------------------------------------- 184 148 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) 185 149 ! ----------------------------------------------------------------------------- 186 187 ! Vertically integrated quantities 188 ! -------------------------------- 189 zua(:,:) = 0.e0 190 zva(:,:) = 0.e0 191 zub(:,:) = 0.e0 192 zvb(:,:) = 0.e0 193 zwx(:,:) = 0.e0 194 zwy(:,:) = 0.e0 195 196 ! vertical sum 197 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 198 DO jk = 1, jpkm1 150 ! 151 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 152 ! ! -------------------------- 153 zua(:,:) = 0.e0 ; zun(:,:) = 0.e0 ; zub(:,:) = 0.e0 154 zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ; zvb(:,:) = 0.e0 155 ! 156 DO jk = 1, jpkm1 157 #if defined key_vectopt_loop 158 DO jj = 1, 1 !Vector opt. => forced unrolling 199 159 DO ji = 1, jpij 200 ! ! Vertically integrated momentum trends 201 zua(ji,1) = zua(ji,1) + fse3u(ji,1,jk) * umask(ji,1,jk) * ua(ji,1,jk) 202 zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 203 ! ! Vertically integrated transports (before) 204 zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 205 zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 206 ! ! Planetary vorticity transport fluxes (now) 207 zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) 208 zwy(ji,1) = zwy(ji,1) + e1v(ji,1) * fse3v(ji,1,jk) * vn(ji,1,jk) 209 END DO 210 END DO 211 ELSE ! No vector opt. 212 DO jk = 1, jpkm1 213 ! ! Vertically integrated momentum trends 214 zua(:,:) = zua(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk) 215 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 216 ! ! Vertically integrated transports (before) 217 zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 218 zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 219 ! ! Planetary vorticity (now) 220 zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) 221 zwy(:,:) = zwy(:,:) + e1v(:,:) * fse3v(:,:,jk) * vn(:,:,jk) 222 END DO 223 ENDIF 224 160 #else 161 DO jj = 1, jpj 162 DO ji = 1, jpi 163 #endif 164 ! ! now trend 165 zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk) 166 zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk) 167 ! ! now velocity 168 zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk) 169 zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk) 170 ! ! before velocity 171 zub(ji,jj) = zub(ji,jj) + fse3u_b(ji,jj,jk) * ub(ji,jj,jk) 172 zvb(ji,jj) = zvb(ji,jj) + fse3v_b(ji,jj,jk) * vb(ji,jj,jk) 173 END DO 174 END DO 175 END DO 176 177 ! !* baroclinic momentum trend (remove the vertical mean trend) 178 DO jk = 1, jpkm1 ! -------------------------- 179 DO jj = 2, jpjm1 180 DO ji = fs_2, fs_jpim1 ! vector opt. 181 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 182 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 183 END DO 184 END DO 185 END DO 186 187 ! !* barotropic Coriolis trends * H (vorticity scheme dependent) 188 ! ! ---------------------------==== 189 zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport 190 zwy(:,:) = zvn(:,:) * e1v(:,:) 191 ! 225 192 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 226 193 DO jj = 2, jpjm1 … … 231 198 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 232 199 ! energy conserving formulation for planetary vorticity term 233 zcu(ji,jj) = z fact2* ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )234 zcv(ji,jj) =-z fact2* ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )200 zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) 201 zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 235 202 END DO 236 203 END DO … … 239 206 DO jj = 2, jpjm1 240 207 DO ji = fs_2, fs_jpim1 ! vector opt. 241 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 242 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 243 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 244 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 208 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 209 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 245 210 zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 246 211 zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) … … 249 214 ! 250 215 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 251 zfac25 = 0.25252 216 DO jj = 2, jpjm1 253 217 DO ji = fs_2, fs_jpim1 ! vector opt. 254 zcu(ji,jj) = + zfac25 / e1u(ji,jj) & 255 & * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 256 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 257 zcv(ji,jj) = - zfac25 / e2v(ji,jj) & 258 & * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 259 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 218 zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 219 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 220 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) & 221 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 260 222 END DO 261 223 END DO … … 263 225 ENDIF 264 226 265 266 ! Remove barotropic trend from general momentum trend 267 ! --------------------------------------------------- 268 DO jk = 1 , jpkm1 269 DO jj = 2, jpjm1 270 DO ji = fs_2, fs_jpim1 ! vector opt. 271 ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj) 272 va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj) 273 END DO 274 END DO 275 END DO 276 277 ! Remove coriolis term from barotropic trend 278 ! ------------------------------------------ 279 DO jj = 2, jpjm1 227 ! !* Right-Hand-Side of the barotropic momentum equation 228 ! ! ---------------------------------------------------- 229 IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient 230 DO jj = 2, jpjm1 231 DO ji = fs_2, fs_jpim1 ! vector opt. 232 zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) & 233 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 234 zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 235 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 236 END DO 237 END DO 238 ENDIF 239 240 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 280 241 DO ji = fs_2, fs_jpim1 281 242 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) … … 284 245 END DO 285 246 247 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 248 ! ! -------------------------- 249 zua(:,:) = zua(:,:) * hur(:,:) 250 zva(:,:) = zva(:,:) * hvr(:,:) 251 ! 252 IF( lk_vvl ) THEN 253 zub(:,:) = zub(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) ) 254 zvb(:,:) = zvb(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) ) 255 ELSE 256 zub(:,:) = zub(:,:) * hur(:,:) 257 zvb(:,:) = zvb(:,:) * hvr(:,:) 258 ENDIF 259 286 260 ! ----------------------------------------------------------------------- 287 261 ! Phase 2 : Integration of the barotropic equations with time splitting 288 262 ! ----------------------------------------------------------------------- 289 290 ! Initialisations 291 !---------------- 292 ! Number of iteration of the barotropic loop 293 icycle = 2 * nn_baro + 1 294 295 ! variables for the barotropic equations 296 zu_sum (:,:) = 0.e0 297 zv_sum (:,:) = 0.e0 298 zssh_sum(:,:) = 0.e0 299 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 300 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 301 zsshb_e (:,:) = sshn_e(:,:) ! (barotropic) sea surface height (before and now) 302 ! vertical sum 303 un_e (:,:) = 0.e0 304 vn_e (:,:) = 0.e0 305 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 306 DO jk = 1, jpkm1 307 DO ji = 1, jpij 308 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 309 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 310 END DO 311 END DO 312 ELSE ! No vector opt. 313 DO jk = 1, jpkm1 314 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 315 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 316 END DO 317 ENDIF 318 zub_e (:,:) = un_e(:,:) 319 zvb_e (:,:) = vn_e(:,:) 320 321 263 ! 264 ! ! ==================== ! 265 ! ! Initialisations ! 266 ! ! ==================== ! 267 icycle = 2 * nn_baro ! Number of barotropic sub time-step 268 269 ! ! Start from NOW field 270 hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points 271 hv_e (:,:) = hv (:,:) 272 hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points 273 hvr_e (:,:) = hvr (:,:) 274 zsshb_e(:,:) = sshn (:,:) ! sea surface height (before and now) 275 sshn_e (:,:) = sshn (:,:) 276 277 zun_e (:,:) = un_b (:,:) ! barotropic velocity (external) 278 zvn_e (:,:) = vn_b (:,:) 279 zub_e (:,:) = un_b (:,:) 280 zvb_e (:,:) = vn_b (:,:) 281 282 zu_sum (:,:) = un_b (:,:) ! summation 283 zv_sum (:,:) = vn_b (:,:) 284 zssh_sum(:,:) = sshn (:,:) 285 286 #if defined key_obc 322 287 ! set ssh corrections to 0 323 288 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 324 #if defined key_obc325 289 IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0 326 290 IF( lp_obc_west ) sshfow_b(:,:) = 0.e0 … … 329 293 #endif 330 294 331 ! Barotropic integration over 2 baroclinic time steps 332 ! --------------------------------------------------- 333 334 ! ! ==================== ! 335 DO jit = 1, icycle ! sub-time-step loop ! 336 ! ! ==================== ! 295 ! ! ==================== ! 296 DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1) 297 ! ! ==================== ! 337 298 z2dt_e = 2. * ( rdt / nn_baro ) 338 IF ( jit == 1 ) z2dt_e = rdt / nn_baro 339 340 ! Time interpolation of open boundary condition data 341 IF( lk_obc ) CALL obc_dta_bt( kt, jit ) 342 IF( lk_bdy .OR. ln_bdy_tides) CALL bdy_dta_bt( kt, jit+1 ) 343 344 ! Horizontal divergence of barotropic transports 345 !-------------------------------------------------- 346 zhdiv(:,:) = 0.e0 347 DO jj = 2, jpjm1 348 DO ji = fs_2, fs_jpim1 ! vector opt. 349 zhdiv(ji,jj) = ( e2u(ji ,jj ) * un_e(ji ,jj) & 350 & -e2u(ji-1,jj ) * un_e(ji-1,jj) & 351 & +e1v(ji ,jj ) * vn_e(ji ,jj) & 352 & -e1v(ji ,jj-1) * vn_e(ji ,jj-1) ) & 353 & / (e1t(ji,jj)*e2t(ji,jj)) 354 END DO 355 END DO 356 299 IF( jn == 1 ) z2dt_e = rdt / nn_baro 300 301 ! !* Update the forcing (OBC, BDY and tides) 302 ! ! ------------------ 303 IF( lk_obc ) CALL obc_dta_bt( kt, jn ) 304 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dta_bt( kt, jn+1 ) 305 306 ! !* after ssh_e 307 ! ! ----------- 308 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 309 DO ji = fs_2, fs_jpim1 ! vector opt. 310 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & 311 & - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) & 312 & + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) & 313 & - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) ) 314 END DO 315 END DO 316 ! 357 317 #if defined key_obc 358 ! open boundaries (div must be zero behind the open boundary)359 360 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 )= 0.e0 ! east361 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 )= 0.e0 ! west318 ! ! OBC : zhdiv must be zero behind the open boundary 319 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 320 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east 321 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west 362 322 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 363 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 )= 0.e0 ! south323 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south 364 324 #endif 365 366 325 #if defined key_bdy 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 370 END DO 371 END DO 326 zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask 372 327 #endif 373 374 ! Sea surface height from the barotropic system 375 !---------------------------------------------- 376 DO jj = 2, jpjm1 377 DO ji = fs_2, fs_jpim1 ! vector opt. 378 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 379 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 380 END DO 381 END DO 382 383 ! evolution of the barotropic transport ( following the vorticity scheme used) 384 ! ---------------------------------------------------------------------------- 385 zwx(:,:) = e2u(:,:) * un_e(:,:) 386 zwy(:,:) = e1v(:,:) * vn_e(:,:) 387 388 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 328 ! 329 DO jj = 2, jpjm1 ! leap-frog on ssh_e 330 DO ji = fs_2, fs_jpim1 ! vector opt. 331 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 332 END DO 333 END DO 334 335 ! !* after barotropic velocities (vorticity scheme dependent) 336 ! ! --------------------------- 337 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 338 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 339 ! 340 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 389 341 DO jj = 2, jpjm1 390 342 DO ji = fs_2, fs_jpim1 ! vector opt. 391 343 ! surface pressure gradient 392 344 IF( lk_vvl) THEN 393 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj )&394 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)395 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1)&396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)345 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 346 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 347 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 348 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 397 349 ELSE 398 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj)/ e1u(ji,jj)399 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj)/ e2v(ji,jj)350 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 351 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 400 352 ENDIF 401 353 ! energy conserving formulation for planetary vorticity term … … 404 356 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj) 405 357 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 406 zcubt = z fact2 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2)407 zcvbt =-z fact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2)408 ! after transports409 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)410 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)358 zcubt = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj) 359 zcvbt =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 360 ! after barotropic velocity 361 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 362 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 411 363 END DO 412 364 END DO 413 365 ! 414 ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme 366 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 367 DO jj = 2, jpjm1 368 DO ji = fs_2, fs_jpim1 ! vector opt. 369 ! surface pressure gradient 370 IF( lk_vvl) THEN 371 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 372 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 373 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 374 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 375 ELSE 376 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 377 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 378 ENDIF 379 ! enstrophy conserving formulation for planetary vorticity term 380 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) 381 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj) 382 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj) 383 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 384 ! after barotropic velocity 385 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 386 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 387 END DO 388 END DO 389 ! 390 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 415 391 DO jj = 2, jpjm1 416 392 DO ji = fs_2, fs_jpim1 ! vector opt. 417 393 ! surface pressure gradient 418 394 IF( lk_vvl) THEN 419 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj )&420 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)421 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1)&422 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)395 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj) 397 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 398 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj) 423 399 ELSE 424 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 425 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 426 ENDIF 427 ! enstrophy conserving formulation for planetary vorticity term 428 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 429 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 430 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 431 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 432 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 433 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 434 ! after transports 435 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 436 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 437 END DO 438 END DO 439 ! 440 ELSEIF ( ln_dynvor_een ) THEN ! energy and enstrophy conserving scheme 441 zfac25 = 0.25 442 DO jj = 2, jpjm1 443 DO ji = fs_2, fs_jpim1 ! vector opt. 444 ! surface pressure gradient 445 IF( lk_vvl) THEN 446 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 447 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 448 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 449 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 450 ELSE 451 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 452 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 400 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj) 401 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 453 402 ENDIF 454 403 ! energy/enstrophy conserving formulation for planetary vorticity term 455 zcubt = + z fac25/ e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &456 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1))457 zcvbt = - z fac25/ e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &458 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ))459 ! after transports460 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)461 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)404 zcubt = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 405 & + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj) 406 zcvbt = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) & 407 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 408 ! after barotropic velocity 409 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 410 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 462 411 END DO 463 412 END DO 464 413 ! 465 414 ENDIF 466 467 ! Flather's boundary condition for the barotropic loop : 468 ! - Update sea surface height on each open boundary 469 ! - Correct the barotropic transports 470 IF( lk_obc ) CALL obc_fla_ts 471 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dyn_fla 472 473 ! ... Boundary conditions on ua_e, va_e, ssha_e 474 CALL lbc_lnk( zua_e , 'U', -1. ) 475 CALL lbc_lnk( zva_e , 'V', -1. ) 476 CALL lbc_lnk( zssha_e, 'T', 1. ) 477 478 ! temporal sum 479 !------------- 480 zu_sum (:,:) = zu_sum (:,:) + zua_e (:,:) 481 zv_sum (:,:) = zv_sum (:,:) + zva_e (:,:) 482 zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 483 484 ! Time filter and swap of dynamics arrays 485 ! --------------------------------------- 486 IF( jit == 1 ) THEN ! Euler (forward) time stepping 487 zsshb_e(:,:) = sshn_e (:,:) 488 zub_e (:,:) = un_e (:,:) 489 zvb_e (:,:) = vn_e (:,:) 490 sshn_e (:,:) = zssha_e(:,:) 491 un_e (:,:) = zua_e (:,:) 492 vn_e (:,:) = zva_e (:,:) 493 ELSE ! Asselin filtering 494 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 495 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * un_e (:,:) 496 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * vn_e (:,:) 497 sshn_e (:,:) = zssha_e(:,:) 498 un_e (:,:) = zua_e (:,:) 499 vn_e (:,:) = zva_e (:,:) 415 ! !* domain lateral boundary 416 ! ! ----------------------- 417 ! ! Flather's boundary condition for the barotropic loop : 418 ! ! - Update sea surface height on each open boundary 419 ! ! - Correct the velocity 420 421 IF( lk_obc ) CALL obc_fla_ts 422 IF( lk_bdy .OR. ln_bdy_tides ) CALL bdy_dyn_fla( sshn_e ) 423 ! 424 CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries 425 CALL lbc_lnk( va_e , 'V', -1. ) 426 CALL lbc_lnk( ssha_e, 'T', 1. ) 427 428 zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps 429 zv_sum (:,:) = zv_sum (:,:) + va_e (:,:) 430 zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:) 431 432 ! !* Time filter and swap 433 ! ! -------------------- 434 IF( jn == 1 ) THEN ! Swap only (1st Euler time step) 435 zsshb_e(:,:) = sshn_e(:,:) 436 zub_e (:,:) = zun_e (:,:) 437 zvb_e (:,:) = zvn_e (:,:) 438 sshn_e (:,:) = ssha_e(:,:) 439 zun_e (:,:) = ua_e (:,:) 440 zvn_e (:,:) = va_e (:,:) 441 ELSE ! Swap + Filter 442 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 443 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) 444 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) 445 sshn_e (:,:) = ssha_e(:,:) 446 zun_e (:,:) = ua_e (:,:) 447 zvn_e (:,:) = va_e (:,:) 500 448 ENDIF 501 449 502 IF( lk_vvl ) THEN ! Variable volume ! needed for BDY ??? 503 504 ! Sea surface elevation 505 ! --------------------- 506 DO jj = 1, jpjm1 507 DO ji = 1,jpim1 508 509 ! Sea Surface Height at u-point before 510 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 511 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 512 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 513 514 ! Sea Surface Height at v-point before 515 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 516 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 517 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 518 519 END DO 520 END DO 521 522 ! Boundaries conditions 523 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 524 525 ! Ocean depth at U- and V-points 526 ! ------------------------------------------- 527 hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 528 hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 529 450 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 451 ! ! ------------------ 452 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 453 DO ji = 1, fs_jpim1 ! Vector opt. 454 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 455 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 456 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 457 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 458 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 459 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 460 END DO 461 END DO 462 CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions 463 CALL lbc_lnk( zsshvn_e, 'V', 1. ) 464 ! 465 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 466 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 467 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) ) 468 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) ) 530 469 ! 531 470 ENDIF … … 534 473 ! ! ==================== ! 535 474 536 ! Time average of after barotropic variables537 zcoef = 1.e0 / ( 2 * nn_baro + 1 )538 un_e (:,:) = zcoef * zu_sum (:,:)539 vn_e (:,:) = zcoef * zv_sum (:,:)540 sshn_e(:,:) = zcoef * zssh_sum(:,:)541 542 475 #if defined key_obc 543 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 476 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ????? 544 477 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:) 545 478 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:) … … 548 481 549 482 ! ----------------------------------------------------------------------------- 550 ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step)483 ! Phase 3. update the general trend with the barotropic trend 551 484 ! ----------------------------------------------------------------------------- 552 553 554 555 ! add time averaged barotropic coriolis and surface pressure gradient 556 ! terms to the general momentum trend 557 ! -------------------------------------------------------------------- 485 ! 486 ! !* Time average ==> after barotropic u, v, ssh 487 zcoef = 1.e0 / ( 2 * nn_baro + 1 ) 488 un_b (:,:) = zcoef * zu_sum (:,:) 489 vn_b (:,:) = zcoef * zv_sum (:,:) 490 sshn_b(:,:) = zcoef * zssh_sum(:,:) 491 ! 492 ! !* update the general momentum trend 558 493 DO jk=1,jpkm1 559 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b560 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b494 ua(:,:,jk) = ua(:,:,jk) + ( un_b(:,:) - zub(:,:) ) / z2dt_b 495 va(:,:,jk) = va(:,:,jk) + ( vn_b(:,:) - zvb(:,:) ) / z2dt_b 561 496 END DO 562 563 ! write filtered free surface arrays in restart file 564 ! -------------------------------------------------- 497 ! 498 ! !* write time-spliting arrays in the restart 565 499 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 566 567 500 ! 568 501 END SUBROUTINE dyn_spg_ts … … 582 515 ! 583 516 IF( TRIM(cdrw) == 'READ' ) THEN 584 IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 585 CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) ) ! free surface and 586 CALL iom_get( numror, jpdom_autoglo, 'un_e' , un_e (:,:) ) ! horizontal transports issued 587 CALL iom_get( numror, jpdom_autoglo, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 517 IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN 518 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued 519 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 588 520 ELSE 589 sshn_e(:,:) = sshn(:,:) 590 un_e (:,:) = 0.e0 591 vn_e (:,:) = 0.e0 521 un_b (:,:) = 0.e0 522 vn_b (:,:) = 0.e0 592 523 ! vertical sum 593 524 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 594 525 DO jk = 1, jpkm1 595 526 DO ji = 1, jpij 596 un_ e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)597 vn_ e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)527 un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 528 vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 598 529 END DO 599 530 END DO 600 531 ELSE ! No vector opt. 601 532 DO jk = 1, jpkm1 602 un_ e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk)603 vn_ e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk)533 un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk) 534 vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 604 535 END DO 605 536 ENDIF 537 un_b (:,:) = un_b(:,:) * hur(:,:) 538 vn_b (:,:) = vn_b(:,:) * hvr(:,:) 606 539 ENDIF 607 540 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 608 CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) ) ! free surface and 609 CALL iom_rstput( kt, nitrst, numrow, 'un_e' , un_e (:,:) ) ! horizontal transports issued 610 CALL iom_rstput( kt, nitrst, numrow, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 541 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity issued 542 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 611 543 ENDIF 612 544 !
Note: See TracChangeset
for help on using the changeset viewer.
