- Timestamp:
- 2020-09-14T17:40:34+02:00 (4 years ago)
- Location:
- NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/OCE/DYN/dynspg_ts.F90 (modified) (38 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev @HEADext/AGRIF5 ^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 9 # SETTE 10 ^/utils/CI/sette@13382 sette
-
- Property svn:externals
-
NEMO/branches/2019/dev_r11351_fldread_with_XIOS/src/OCE/DYN/dynspg_ts.F90
r11405 r13463 1 1 MODULE dynspg_ts 2 2 3 !! Includes ROMS wd scheme with diagnostic outputs ; un and uaupdates are commented out !3 !! Includes ROMS wd scheme with diagnostic outputs ; puu(:,:,:,Kmm) and puu(:,:,:,Krhs) updates are commented out ! 4 4 5 5 !!====================================================================== … … 31 31 USE dom_oce ! ocean space and time domain 32 32 USE sbc_oce ! surface boundary condition: ocean 33 USE isf_oce ! ice shelf variable (fwfisf) 33 34 USE zdf_oce ! vertical physics: variables 34 35 USE zdfdrg ! vertical physics: top/bottom drag coef. 35 USE sbcisf ! ice shelf variable (fwfisf)36 36 USE sbcapr ! surface boundary condition: atmospheric pressure 37 37 USE dynadv , ONLY: ln_dynadv_vec … … 44 44 USE bdytides ! open boundary condition data 45 45 USE bdydyn2d ! open boundary conditions on barotropic variables 46 USE sbctide ! tides 47 USE updtide ! tide potential 46 USE tide_mod ! 48 47 USE sbcwave ! surface wave 49 USE diatmb ! Top,middle,bottom output50 48 #if defined key_agrif 51 49 USE agrif_oce_interp ! agrif … … 62 60 USE iom ! IOM library 63 61 USE restart ! only for lrst_oce 64 USE diatmb ! Top,middle,bottom output 62 63 USE iom ! to remove 65 64 66 65 IMPLICIT NONE … … 73 72 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 74 73 ! 75 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_ baro <= 2.5 nn_baro76 REAL(wp),SAVE :: r dtbt! Barotropic time step74 INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e 75 REAL(wp),SAVE :: rDt_e ! Barotropic time step 77 76 ! 78 77 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 ! 1st & 2nd weights used in time filtering of barotropic fields … … 87 86 88 87 !! * Substitutions 89 # include "vectopt_loop_substitute.h90" 88 # include "do_loop_substitute.h90" 89 # include "domzgr_substitute.h90" 90 90 !!---------------------------------------------------------------------- 91 91 !! NEMO/OCE 4.0 , NEMO Consortium (2018) … … 103 103 ierr(:) = 0 104 104 ! 105 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 106 ! 105 ALLOCATE( wgtbtp1(3*nn_e), wgtbtp2(3*nn_e), zwz(jpi,jpj), STAT=ierr(1) ) 107 106 IF( ln_dynvor_een .OR. ln_dynvor_eeT ) & 108 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 109 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 107 & ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) ) 110 108 ! 111 109 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) … … 119 117 120 118 121 SUBROUTINE dyn_spg_ts( kt )119 SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) 122 120 !!---------------------------------------------------------------------- 123 121 !! … … 134 132 !! 135 133 !! ** Action : 136 !! -Update the filtered free surface at step "n+1" : ssha137 !! -Update filtered barotropic velocities at step "n+1" : ua_b, va_b134 !! -Update the filtered free surface at step "n+1" : pssh(:,:,Kaa) 135 !! -Update filtered barotropic velocities at step "n+1" : puu_b(:,:,:,Kaa), vv_b(:,:,:,Kaa) 138 136 !! -Compute barotropic advective fluxes at step "n" : un_adv, vn_adv 139 137 !! These are used to advect tracers and are compliant with discrete 140 138 !! continuity equation taken at the baroclinic time steps. This 141 139 !! ensures tracers conservation. 142 !! - ( ua, va) momentum trend updated with barotropic component.140 !! - (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) momentum trend updated with barotropic component. 143 141 !! 144 142 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 145 143 !!--------------------------------------------------------------------- 146 INTEGER, INTENT(in) :: kt ! ocean time-step index 144 INTEGER , INTENT( in ) :: kt ! ocean time-step index 145 INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices 146 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation 147 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels 147 148 ! 148 149 INTEGER :: ji, jj, jk, jn ! dummy loop indices 149 150 LOGICAL :: ll_fw_start ! =T : forward integration 150 151 LOGICAL :: ll_init ! =T : special startup of 2d equations 151 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 152 INTEGER :: ikbu, iktu, noffset ! local integers 153 INTEGER :: ikbv, iktv ! - - 154 REAL(wp) :: r1_2dt_b, z2dt_bf ! local scalars 155 REAL(wp) :: zx1, zx2, zu_spg, zhura, z1_hu ! - - 156 REAL(wp) :: zy1, zy2, zv_spg, zhvra, z1_hv ! - - 152 INTEGER :: noffset ! local integers : time offset for bdy update 153 REAL(wp) :: r1_Dt_b, z1_hu, z1_hv ! local scalars 157 154 REAL(wp) :: za0, za1, za2, za3 ! - - 158 REAL(wp) :: zmdi, zztmp , z1_ht ! - - 159 REAL(wp), DIMENSION(jpi,jpj) :: zsshp2_e, zhf 160 REAL(wp), DIMENSION(jpi,jpj) :: zwx, zu_trd, zu_frc, zssh_frc 161 REAL(wp), DIMENSION(jpi,jpj) :: zwy, zv_trd, zv_frc, zhdiv 162 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhust_e, zhtp2_e 163 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zhvst_e 155 REAL(wp) :: zztmp, zldg ! - - 156 REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - - 157 REAL(wp) :: zun_save, zvn_save ! - - 158 REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg, zssh_frc 159 REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg 160 REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e 161 REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e 164 162 REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points 163 REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes 164 REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v 165 165 ! 166 166 REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True. … … 172 172 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points 173 173 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask 174 REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep 174 175 !!---------------------------------------------------------------------- 175 176 ! … … 178 179 IF( ln_wd_dl ) ALLOCATE( ztwdmask(jpi,jpj), zuwdmask(jpi,jpj), zvwdmask(jpi,jpj), zuwdav2(jpi,jpj), zvwdav2(jpi,jpj)) 179 180 ! 180 zmdi=1.e+20 ! missing data indicator for masking181 !182 181 zwdramp = r_rn_wdmin1 ! simplest ramp 183 182 ! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp 184 ! ! reciprocal of baroclinic time step 185 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 186 ELSE ; z2dt_bf = 2.0_wp * rdt 187 ENDIF 188 r1_2dt_b = 1.0_wp / z2dt_bf 183 ! ! inverse of baroclinic time step 184 r1_Dt_b = 1._wp / rDt 189 185 ! 190 186 ll_init = ln_bt_av ! if no time averaging, then no specific restart 191 187 ll_fw_start = .FALSE. 192 188 ! ! time offset in steps for bdy data update 193 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_ baro189 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_e 194 190 ELSE ; noffset = 0 195 191 ENDIF … … 202 198 IF(lwp) WRITE(numout,*) 203 199 ! 204 IF( neuler == 0) ll_init=.TRUE.205 ! 206 IF( ln_bt_fw .OR. neuler == 0) THEN200 IF( l_1st_euler ) ll_init=.TRUE. 201 ! 202 IF( ln_bt_fw .OR. l_1st_euler ) THEN 207 203 ll_fw_start =.TRUE. 208 204 noffset = 0 … … 210 206 ll_fw_start =.FALSE. 211 207 ENDIF 212 ! 213 ! Set averaging weights and cycle length: 208 ! ! Set averaging weights and cycle length: 214 209 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 215 210 ! 216 ENDIF 217 ! 218 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 219 DO jj = 2, jpjm1 220 DO ji = fs_2, fs_jpim1 ! vector opt. 221 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 222 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 223 END DO 224 END DO 225 ELSE ! bottom friction only 226 DO jj = 2, jpjm1 227 DO ji = fs_2, fs_jpim1 ! vector opt. 228 zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 229 zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 230 END DO 231 END DO 232 ENDIF 233 ! 234 ! Set arrays to remove/compute coriolis trend. 235 ! Do it once at kt=nit000 if volume is fixed, else at each long time step. 236 ! Note that these arrays are also used during barotropic loop. These are however frozen 237 ! although they should be updated in the variable volume case. Not a big approximation. 238 ! To remove this approximation, copy lines below inside barotropic loop 239 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 240 ! 241 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 242 ! 243 SELECT CASE( nvor_scheme ) 244 CASE( np_EEN ) != EEN scheme using e3f (energy & enstrophy scheme) 245 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 246 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 247 DO jj = 1, jpjm1 248 DO ji = 1, jpim1 249 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 250 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 251 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 252 END DO 253 END DO 254 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 255 DO jj = 1, jpjm1 256 DO ji = 1, jpim1 257 zwz(ji,jj) = ( ht_n (ji ,jj+1) + ht_n (ji+1,jj+1) & 258 & + ht_n (ji ,jj ) + ht_n (ji+1,jj ) ) & 259 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 260 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 261 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 262 END DO 263 END DO 264 END SELECT 265 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 266 ! 267 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 268 DO jj = 2, jpj 269 DO ji = 2, jpi 270 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 271 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 272 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 273 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 274 END DO 275 END DO 276 ! 277 CASE( np_EET ) != EEN scheme using e3t (energy conserving scheme) 278 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 279 DO jj = 2, jpj 280 DO ji = 2, jpi 281 z1_ht = ssmask(ji,jj) / ( ht_n(ji,jj) + 1._wp - ssmask(ji,jj) ) 282 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 283 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 284 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 285 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 286 END DO 287 END DO 288 ! 289 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 290 ! 291 zwz(:,:) = 0._wp 292 zhf(:,:) = 0._wp 293 294 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 295 !!gm A priori a better value should be something like : 296 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 297 !!gm divided by the sum of the corresponding mask 298 !!gm 299 !! 300 IF( .NOT.ln_sco ) THEN 301 302 !!gm agree the JC comment : this should be done in a much clear way 303 304 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 305 ! Set it to zero for the time being 306 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 307 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 308 ! ENDIF 309 ! zhf(:,:) = gdepw_0(:,:,jk+1) 310 ! 311 ELSE 312 ! 313 !zhf(:,:) = hbatf(:,:) 314 DO jj = 1, jpjm1 315 DO ji = 1, jpim1 316 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 317 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 318 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 319 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 320 END DO 321 END DO 322 ENDIF 323 ! 324 DO jj = 1, jpjm1 325 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 326 END DO 327 ! 328 DO jk = 1, jpkm1 329 DO jj = 1, jpjm1 330 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 331 END DO 332 END DO 333 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 334 ! JC: TBC. hf should be greater than 0 335 DO jj = 1, jpj 336 DO ji = 1, jpi 337 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) ! zhf is actually hf here but it saves an array 338 END DO 339 END DO 340 zwz(:,:) = ff_f(:,:) * zwz(:,:) 341 END SELECT 342 ENDIF 343 ! 344 ! If forward start at previous time step, and centered integration, 345 ! then update averaging weights: 346 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 347 ll_fw_start=.FALSE. 348 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 349 ENDIF 350 211 ELSEIF( kt == nit000 + 1 ) THEN !* initialisation 2nd time-step 212 ! 213 IF( .NOT.ln_bt_fw ) THEN 214 ! If we did an Euler timestep on the first timestep we need to reset ll_fw_start 215 ! and the averaging weights. We don't have an easy way of telling whether we did 216 ! an Euler timestep on the first timestep (because l_1st_euler is reset to .false. 217 ! at the end of the first timestep) so just do this in all cases. 218 ll_fw_start = .FALSE. 219 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 220 ENDIF 221 ! 222 ENDIF 223 ! 351 224 ! ----------------------------------------------------------------------------- 352 225 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 354 227 ! 355 228 ! 356 ! !* e3*d/dt(Ua) (Vertically integrated) 357 ! ! -------------------------------------------------- 358 zu_frc(:,:) = 0._wp 359 zv_frc(:,:) = 0._wp 360 ! 361 DO jk = 1, jpkm1 362 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 363 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 229 ! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends) 230 ! ! --------------------------- ! 231 DO jk = 1 , jpk 232 ze3u(:,:,jk) = e3u(:,:,jk,Kmm) 233 ze3v(:,:,jk) = e3v(:,:,jk,Kmm) 364 234 END DO 365 235 ! 366 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 367 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 368 ! 369 ! 370 ! !* baroclinic momentum trend (remove the vertical mean trend) 371 DO jk = 1, jpkm1 ! ----------------------------------------------------------- 372 DO jj = 2, jpjm1 373 DO ji = fs_2, fs_jpim1 ! vector opt. 374 ua(ji,jj,jk) = ua(ji,jj,jk) - zu_frc(ji,jj) * umask(ji,jj,jk) 375 va(ji,jj,jk) = va(ji,jj,jk) - zv_frc(ji,jj) * vmask(ji,jj,jk) 376 END DO 377 END DO 236 zu_frc(:,:) = SUM( ze3u(:,:,:) * uu(:,:,:,Krhs) * umask(:,:,:) , DIM=3 ) * r1_hu(:,:,Kmm) 237 zv_frc(:,:) = SUM( ze3v(:,:,:) * vv(:,:,:,Krhs) * vmask(:,:,:) , DIM=3 ) * r1_hv(:,:,Kmm) 238 ! 239 ! 240 ! != U(Krhs) => baroclinic trend =! (remove its vertical mean) 241 DO jk = 1, jpkm1 ! ----------------------------- ! 242 uu(:,:,jk,Krhs) = ( uu(:,:,jk,Krhs) - zu_frc(:,:) ) * umask(:,:,jk) 243 vv(:,:,jk,Krhs) = ( vv(:,:,jk,Krhs) - zv_frc(:,:) ) * vmask(:,:,jk) 378 244 END DO 379 245 … … 381 247 !!gm Is it correct to do so ? I think so... 382 248 383 384 ! !* barotropic Coriolis trends (vorticity scheme dependent) 385 ! ! -------------------------------------------------------- 386 ! 387 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 388 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 389 ! 390 SELECT CASE( nvor_scheme ) 391 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 392 DO jj = 2, jpjm1 393 DO ji = 2, jpim1 ! vector opt. 394 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * r1_hu_n(ji,jj) & 395 & * ( e1e2t(ji+1,jj)*ht_n(ji+1,jj)*ff_t(ji+1,jj) * ( vn_b(ji+1,jj) + vn_b(ji+1,jj-1) ) & 396 & + e1e2t(ji ,jj)*ht_n(ji ,jj)*ff_t(ji ,jj) * ( vn_b(ji ,jj) + vn_b(ji ,jj-1) ) ) 397 ! 398 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * r1_hv_n(ji,jj) & 399 & * ( e1e2t(ji,jj+1)*ht_n(ji,jj+1)*ff_t(ji,jj+1) * ( un_b(ji,jj+1) + un_b(ji-1,jj+1) ) & 400 & + e1e2t(ji,jj )*ht_n(ji,jj )*ff_t(ji,jj ) * ( un_b(ji,jj ) + un_b(ji-1,jj ) ) ) 401 END DO 402 END DO 403 ! 404 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 405 DO jj = 2, jpjm1 406 DO ji = fs_2, fs_jpim1 ! vector opt. 407 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 408 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 409 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 410 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 411 ! energy conserving formulation for planetary vorticity term 412 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 413 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 414 END DO 415 END DO 416 ! 417 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 418 DO jj = 2, jpjm1 419 DO ji = fs_2, fs_jpim1 ! vector opt. 420 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 421 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 422 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 423 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 424 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 425 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 426 END DO 427 END DO 428 ! 429 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 430 DO jj = 2, jpjm1 431 DO ji = fs_2, fs_jpim1 ! vector opt. 432 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 433 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 434 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 435 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 436 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 437 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 438 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 439 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 440 END DO 441 END DO 442 ! 443 END SELECT 444 ! 445 ! !* Right-Hand-Side of the barotropic momentum equation 446 ! ! ---------------------------------------------------- 447 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 448 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 449 DO jj = 2, jpjm1 450 DO ji = 2, jpim1 451 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji+1,jj) ) > & 452 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 453 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji+1,jj) + ht_0(ji+1,jj) ) & 454 & > rn_wdmin1 + rn_wdmin2 455 ll_tmp2 = ( ABS( sshn(ji+1,jj) - sshn(ji ,jj)) > 1.E-12 ).AND.( & 456 & MAX( sshn(ji,jj) , sshn(ji+1,jj) ) > & 457 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 458 IF(ll_tmp1) THEN 459 zcpx(ji,jj) = 1.0_wp 460 ELSEIF(ll_tmp2) THEN 461 ! no worries about sshn(ji+1,jj) - sshn(ji ,jj) = 0, it won't happen ! here 462 zcpx(ji,jj) = ABS( (sshn(ji+1,jj) + ht_0(ji+1,jj) - sshn(ji,jj) - ht_0(ji,jj)) & 463 & / (sshn(ji+1,jj) - sshn(ji ,jj)) ) 464 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 465 ELSE 466 zcpx(ji,jj) = 0._wp 467 ENDIF 468 ! 469 ll_tmp1 = MIN( sshn(ji,jj) , sshn(ji,jj+1) ) > & 470 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 471 & MAX( sshn(ji,jj) + ht_0(ji,jj) , sshn(ji,jj+1) + ht_0(ji,jj+1) ) & 472 & > rn_wdmin1 + rn_wdmin2 473 ll_tmp2 = ( ABS( sshn(ji,jj) - sshn(ji,jj+1)) > 1.E-12 ).AND.( & 474 & MAX( sshn(ji,jj) , sshn(ji,jj+1) ) > & 475 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 476 477 IF(ll_tmp1) THEN 478 zcpy(ji,jj) = 1.0_wp 479 ELSE IF(ll_tmp2) THEN 480 ! no worries about sshn(ji,jj+1) - sshn(ji,jj ) = 0, it won't happen ! here 481 zcpy(ji,jj) = ABS( (sshn(ji,jj+1) + ht_0(ji,jj+1) - sshn(ji,jj) - ht_0(ji,jj)) & 482 & / (sshn(ji,jj+1) - sshn(ji,jj )) ) 483 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 484 ELSE 485 zcpy(ji,jj) = 0._wp 486 ENDIF 487 END DO 488 END DO 489 ! 490 DO jj = 2, jpjm1 491 DO ji = 2, jpim1 492 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 493 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 494 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 495 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 496 END DO 497 END DO 498 ! 499 ELSE 500 ! 501 DO jj = 2, jpjm1 502 DO ji = fs_2, fs_jpim1 ! vector opt. 503 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 504 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 505 END DO 506 END DO 507 ENDIF 508 ! 509 ENDIF 510 ! 511 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 512 DO ji = fs_2, fs_jpim1 513 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 514 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 515 END DO 516 END DO 517 ! 518 ! ! Add bottom stress contribution from baroclinic velocities: 519 IF (ln_bt_fw) THEN 520 DO jj = 2, jpjm1 521 DO ji = fs_2, fs_jpim1 ! vector opt. 522 ikbu = mbku(ji,jj) 523 ikbv = mbkv(ji,jj) 524 zwx(ji,jj) = un(ji,jj,ikbu) - un_b(ji,jj) ! NOW bottom baroclinic velocities 525 zwy(ji,jj) = vn(ji,jj,ikbv) - vn_b(ji,jj) 526 END DO 527 END DO 249 ! != remove 2D Coriolis and pressure gradient trends =! 250 ! ! ------------------------------------------------- ! 251 ! 252 IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient 253 ! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes 254 ! 255 ! !* 2D Coriolis trends 256 zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes 257 zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column 258 ! 259 CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in 260 & zu_trd, zv_trd ) ! ==>> out 261 ! 262 IF( .NOT.ln_linssh ) THEN !* surface pressure gradient (variable volume only) 263 ! 264 IF( ln_wd_il ) THEN ! W/D : limiter applied to spgspg 265 CALL wad_spg( pssh(:,:,Kmm), zcpx, zcpy ) ! Calculating W/D gravity filters, zcpx and zcpy 266 DO_2D( 0, 0, 0, 0 ) 267 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) & 268 & * r1_e1u(ji,jj) * zcpx(ji,jj) * wdrampu(ji,jj) !jth 269 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) & 270 & * r1_e2v(ji,jj) * zcpy(ji,jj) * wdrampv(ji,jj) !jth 271 END_2D 272 ELSE ! now suface pressure gradient 273 DO_2D( 0, 0, 0, 0 ) 274 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( pssh(ji+1,jj ,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e1u(ji,jj) 275 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( pssh(ji ,jj+1,Kmm) - pssh(ji ,jj ,Kmm) ) * r1_e2v(ji,jj) 276 END_2D 277 ENDIF 278 ! 279 ENDIF 280 ! 281 DO_2D( 0, 0, 0, 0 ) 282 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 283 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 284 END_2D 285 ! 286 ! != Add bottom stress contribution from baroclinic velocities =! 287 ! ! ----------------------------------------------------------- ! 288 CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v ) ! also provide the barotropic drag coefficients 289 ! != Add atmospheric pressure forcing =! 290 ! ! ---------------------------------- ! 291 IF( ln_apr_dyn ) THEN 292 IF( ln_bt_fw ) THEN ! FORWARD integration: use kt+1/2 pressure (NOW+1/2) 293 DO_2D( 0, 0, 0, 0 ) 294 zu_frc(ji,jj) = zu_frc(ji,jj) + grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 295 zv_frc(ji,jj) = zv_frc(ji,jj) + grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 296 END_2D 297 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 pressure (NOW) 298 zztmp = grav * r1_2 299 DO_2D( 0, 0, 0, 0 ) 300 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 301 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 302 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 303 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 304 END_2D 305 ENDIF 306 ENDIF 307 ! 308 ! != Add atmospheric pressure forcing =! 309 ! ! ---------------------------------- ! 310 IF( ln_bt_fw ) THEN ! Add wind forcing 311 DO_2D( 0, 0, 0, 0 ) 312 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rho0 * utau(ji,jj) * r1_hu(ji,jj,Kmm) 313 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rho0 * vtau(ji,jj) * r1_hv(ji,jj,Kmm) 314 END_2D 528 315 ELSE 529 DO jj = 2, jpjm1 530 DO ji = fs_2, fs_jpim1 ! vector opt. 531 ikbu = mbku(ji,jj) 532 ikbv = mbkv(ji,jj) 533 zwx(ji,jj) = ub(ji,jj,ikbu) - ub_b(ji,jj) ! BEFORE bottom baroclinic velocities 534 zwy(ji,jj) = vb(ji,jj,ikbv) - vb_b(ji,jj) 535 END DO 536 END DO 537 ENDIF 538 ! 539 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 540 IF( ln_wd_il ) THEN 541 zztmp = -1._wp / rdtbt 542 DO jj = 2, jpjm1 543 DO ji = fs_2, fs_jpim1 ! vector opt. 544 zu_frc(ji,jj) = zu_frc(ji,jj) + & 545 & MAX(r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ), zztmp ) * zwx(ji,jj) * wdrampu(ji,jj) 546 zv_frc(ji,jj) = zv_frc(ji,jj) + & 547 & MAX(r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ), zztmp ) * zwy(ji,jj) * wdrampv(ji,jj) 548 END DO 549 END DO 550 ELSE 551 DO jj = 2, jpjm1 552 DO ji = fs_2, fs_jpim1 ! vector opt. 553 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zwx(ji,jj) 554 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zwy(ji,jj) 555 END DO 556 END DO 316 zztmp = r1_rho0 * r1_2 317 DO_2D( 0, 0, 0, 0 ) 318 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu(ji,jj,Kmm) 319 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv(ji,jj,Kmm) 320 END_2D 321 ENDIF 322 ! 323 ! !----------------! 324 ! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain) 325 ! !----------------! 326 ! != Net water flux forcing applied to a water column =! 327 ! ! --------------------------------------------------- ! 328 IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2) 329 zssh_frc(:,:) = r1_rho0 * ( emp(:,:) - rnf(:,:) + fwfisf_cav(:,:) + fwfisf_par(:,:) ) 330 ELSE ! CENTRED integration: use kt-1/2 + kt+1/2 fluxes (NOW) 331 zztmp = r1_rho0 * r1_2 332 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) & 333 & - rnf(:,:) - rnf_b(:,:) & 334 & + fwfisf_cav(:,:) + fwfisf_cav_b(:,:) & 335 & + fwfisf_par(:,:) + fwfisf_par_b(:,:) ) 336 ENDIF 337 ! != Add Stokes drift divergence =! (if exist) 338 IF( ln_sdw ) THEN ! ----------------------------- ! 339 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:) 340 ENDIF 341 ! 342 ! ! ice sheet coupling 343 IF ( ln_isf .AND. ln_isfcpl ) THEN 344 ! 345 ! ice sheet coupling 346 IF( ln_rstart .AND. kt == nit000 ) THEN 347 zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_ssh(:,:) 348 END IF 349 ! 350 ! conservation option 351 IF( ln_isfcpl_cons ) THEN 352 zssh_frc(:,:) = zssh_frc(:,:) + risfcpl_cons_ssh(:,:) 353 END IF 354 ! 557 355 END IF 558 356 ! 559 IF( ln_isfcav ) THEN ! Add TOP stress contribution from baroclinic velocities:560 IF( ln_bt_fw ) THEN561 DO jj = 2, jpjm1562 DO ji = fs_2, fs_jpim1 ! vector opt.563 iktu = miku(ji,jj)564 iktv = mikv(ji,jj)565 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities566 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj)567 END DO568 END DO569 ELSE570 DO jj = 2, jpjm1571 DO ji = fs_2, fs_jpim1 ! vector opt.572 iktu = miku(ji,jj)573 iktv = mikv(ji,jj)574 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities575 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj)576 END DO577 END DO578 ENDIF579 !580 ! Note that the "unclipped" top friction parameter is used even with explicit drag581 DO jj = 2, jpjm1582 DO ji = fs_2, fs_jpim1 ! vector opt.583 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_hu_n(ji,jj) * r1_2 * ( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zwx(ji,jj)584 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_hv_n(ji,jj) * r1_2 * ( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zwy(ji,jj)585 END DO586 END DO587 ENDIF588 !589 IF( ln_bt_fw ) THEN ! Add wind forcing590 DO jj = 2, jpjm1591 DO ji = fs_2, fs_jpim1 ! vector opt.592 zu_frc(ji,jj) = zu_frc(ji,jj) + r1_rau0 * utau(ji,jj) * r1_hu_n(ji,jj)593 zv_frc(ji,jj) = zv_frc(ji,jj) + r1_rau0 * vtau(ji,jj) * r1_hv_n(ji,jj)594 END DO595 END DO596 ELSE597 zztmp = r1_rau0 * r1_2598 DO jj = 2, jpjm1599 DO ji = fs_2, fs_jpim1 ! vector opt.600 zu_frc(ji,jj) = zu_frc(ji,jj) + zztmp * ( utau_b(ji,jj) + utau(ji,jj) ) * r1_hu_n(ji,jj)601 zv_frc(ji,jj) = zv_frc(ji,jj) + zztmp * ( vtau_b(ji,jj) + vtau(ji,jj) ) * r1_hv_n(ji,jj)602 END DO603 END DO604 ENDIF605 !606 IF( ln_apr_dyn ) THEN ! Add atm pressure forcing607 IF( ln_bt_fw ) THEN608 DO jj = 2, jpjm1609 DO ji = fs_2, fs_jpim1 ! vector opt.610 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj)611 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj)612 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg613 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg614 END DO615 END DO616 ELSE617 zztmp = grav * r1_2618 DO jj = 2, jpjm1619 DO ji = fs_2, fs_jpim1 ! vector opt.620 zu_spg = zztmp * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) &621 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj)622 zv_spg = zztmp * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) &623 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj)624 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg625 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg626 END DO627 END DO628 ENDIF629 ENDIF630 ! !* Right-Hand-Side of the barotropic ssh equation631 ! ! -----------------------------------------------632 ! ! Surface net water flux and rivers633 IF (ln_bt_fw) THEN634 zssh_frc(:,:) = r1_rau0 * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) )635 ELSE636 zztmp = r1_rau0 * r1_2637 zssh_frc(:,:) = zztmp * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) &638 & + fwfisf(:,:) + fwfisf_b(:,:) )639 ENDIF640 !641 IF( ln_sdw ) THEN ! Stokes drift divergence added if necessary642 zssh_frc(:,:) = zssh_frc(:,:) + div_sd(:,:)643 ENDIF644 !645 357 #if defined key_asminc 646 ! ! Include the IAU weighted SSH increment 358 ! != Add the IAU weighted SSH increment =! 359 ! ! ------------------------------------ ! 647 360 IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN 648 361 zssh_frc(:,:) = zssh_frc(:,:) - ssh_iau(:,:) 649 362 ENDIF 650 363 #endif 651 ! ! *Fill boundary data arrays for AGRIF364 ! != Fill boundary data arrays for AGRIF 652 365 ! ! ------------------------------------ 653 366 #if defined key_agrif … … 671 384 vb_e (:,:) = 0._wp 672 385 ENDIF 673 386 ! 387 IF( ln_linssh ) THEN ! mid-step ocean depth is fixed (hup2_e=hu_n=hu_0) 388 zhup2_e(:,:) = hu(:,:,Kmm) 389 zhvp2_e(:,:) = hv(:,:,Kmm) 390 zhtp2_e(:,:) = ht(:,:) 391 ENDIF 674 392 ! 675 393 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 676 sshn_e(:,:) = sshn(:,:)677 un_e (:,:) = un_b(:,:)678 vn_e (:,:) = vn_b(:,:)679 ! 680 hu_e (:,:) = hu _n(:,:)681 hv_e (:,:) = hv _n(:,:)682 hur_e (:,:) = r1_hu _n(:,:)683 hvr_e (:,:) = r1_hv _n(:,:)394 sshn_e(:,:) = pssh(:,:,Kmm) 395 un_e (:,:) = puu_b(:,:,Kmm) 396 vn_e (:,:) = pvv_b(:,:,Kmm) 397 ! 398 hu_e (:,:) = hu(:,:,Kmm) 399 hv_e (:,:) = hv(:,:,Kmm) 400 hur_e (:,:) = r1_hu(:,:,Kmm) 401 hvr_e (:,:) = r1_hv(:,:,Kmm) 684 402 ELSE ! CENTRED integration: start from BEFORE fields 685 sshn_e(:,:) = sshb(:,:) 686 un_e (:,:) = ub_b(:,:) 687 vn_e (:,:) = vb_b(:,:) 688 ! 689 hu_e (:,:) = hu_b(:,:) 690 hv_e (:,:) = hv_b(:,:) 691 hur_e (:,:) = r1_hu_b(:,:) 692 hvr_e (:,:) = r1_hv_b(:,:) 693 ENDIF 694 ! 695 ! 403 sshn_e(:,:) = pssh(:,:,Kbb) 404 un_e (:,:) = puu_b(:,:,Kbb) 405 vn_e (:,:) = pvv_b(:,:,Kbb) 406 ! 407 hu_e (:,:) = hu(:,:,Kbb) 408 hv_e (:,:) = hv(:,:,Kbb) 409 hur_e (:,:) = r1_hu(:,:,Kbb) 410 hvr_e (:,:) = r1_hv(:,:,Kbb) 411 ENDIF 696 412 ! 697 413 ! Initialize sums: 698 ua_b (:,:) = 0._wp ! After barotropic velocities (or transport if flux form)699 va_b (:,:) = 0._wp700 ssha (:,:) = 0._wp ! Sum for after averaged sea level414 puu_b (:,:,Kaa) = 0._wp ! After barotropic velocities (or transport if flux form) 415 pvv_b (:,:,Kaa) = 0._wp 416 pssh (:,:,Kaa) = 0._wp ! Sum for after averaged sea level 701 417 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 702 418 vn_adv(:,:) = 0._wp … … 714 430 ! 715 431 l_full_nf_update = jn == icycle ! false: disable full North fold update (performances) for jn = 1 to icycle-1 716 ! ! ------------------ 717 ! !* Update the forcing (BDY and tides) 718 ! ! ------------------ 719 ! Update only tidal forcing at open boundaries 720 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 721 IF( ln_tide_pot .AND. ln_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 722 ! 723 ! Set extrapolation coefficients for predictor step: 432 ! 433 ! !== Update the forcing ==! (BDY and tides) 434 ! 435 IF( ln_bdy .AND. ln_tide ) CALL bdy_dta_tides( kt, kit=jn, pt_offset= REAL(noffset+1,wp) ) 436 ! Update tide potential at the beginning of current time substep 437 IF( ln_tide_pot .AND. ln_tide ) THEN 438 zt0substep = REAL(nsec_day, wp) - 0.5_wp*rn_Dt + (jn + noffset - 1) * rn_Dt / REAL(nn_e, wp) 439 CALL upd_tide(zt0substep, Kmm) 440 END IF 441 ! 442 ! !== extrapolation at mid-step ==! (jn+1/2) 443 ! 444 ! !* Set extrapolation coefficients for predictor step: 724 445 IF ((jn<3).AND.ll_init) THEN ! Forward 725 446 za1 = 1._wp … … 731 452 za3 = 0.281105_wp ! za3 = bet 732 453 ENDIF 733 734 ! Extrapolate barotropic velocities at step jit+0.5: 454 ! 455 ! !* Extrapolate barotropic velocities at mid-step (jn+1/2) 456 !-- m+1/2 m m-1 m-2 --! 457 !-- u = (3/2+beta) u -(1/2+2beta) u + beta u --! 458 !-------------------------------------------------------------------------! 735 459 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 736 460 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) … … 739 463 ! ! ------------------ 740 464 ! Extrapolate Sea Level at step jit+0.5: 465 !-- m+1/2 m m-1 m-2 --! 466 !-- ssh = (3/2+beta) ssh -(1/2+2beta) ssh + beta ssh --! 467 !--------------------------------------------------------------------------------! 741 468 zsshp2_e(:,:) = za1 * sshn_e(:,:) + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 742 469 743 ! set wetting & drying mask at tracer points for this barotropic sub-step 744 IF ( ln_wd_dl ) THEN 745 ! 746 IF ( ln_wd_dl_rmp ) THEN 747 DO jj = 1, jpj 748 DO ji = 1, jpi ! vector opt. 749 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 750 ! IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 751 ztwdmask(ji,jj) = 1._wp 752 ELSE IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 753 ztwdmask(ji,jj) = (tanh(50._wp*( ( zsshp2_e(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1)) ) 754 ELSE 755 ztwdmask(ji,jj) = 0._wp 756 END IF 757 END DO 758 END DO 759 ELSE 760 DO jj = 1, jpj 761 DO ji = 1, jpi ! vector opt. 762 IF ( zsshp2_e(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 763 ztwdmask(ji,jj) = 1._wp 764 ELSE 765 ztwdmask(ji,jj) = 0._wp 766 ENDIF 767 END DO 768 END DO 769 ENDIF 770 ! 771 ENDIF 470 ! set wetting & drying mask at tracer points for this barotropic mid-step 471 IF( ln_wd_dl ) CALL wad_tmsk( zsshp2_e, ztwdmask ) 772 472 ! 773 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 774 DO ji = 2, fs_jpim1 ! Vector opt. 775 zwx(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 776 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 777 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 778 zwy(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 779 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 780 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 781 END DO 782 END DO 783 CALL lbc_lnk_multi( 'dynspg_ts', zwx, 'U', 1._wp, zwy, 'V', 1._wp ) 473 ! ! ocean t-depth at mid-step 474 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 784 475 ! 785 zhup2_e(:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 786 zhvp2_e(:,:) = hv_0(:,:) + zwy(:,:) 787 zhtp2_e(:,:) = ht_0(:,:) + zsshp2_e(:,:) 788 ELSE 789 zhup2_e(:,:) = hu_n(:,:) 790 zhvp2_e(:,:) = hv_n(:,:) 791 zhtp2_e(:,:) = ht_n(:,:) 792 ENDIF 793 ! !* after ssh 794 ! ! ----------- 795 ! 796 ! Enforce volume conservation at open boundaries: 476 ! ! ocean u- and v-depth at mid-step (separate DO-loops remove the need of a lbc_lnk) 477 DO_2D( 1, 1, 1, 0 ) 478 zhup2_e(ji,jj) = hu_0(ji,jj) + r1_2 * r1_e1e2u(ji,jj) & 479 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 480 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 481 END_2D 482 DO_2D( 1, 0, 1, 1 ) 483 zhvp2_e(ji,jj) = hv_0(ji,jj) + r1_2 * r1_e1e2v(ji,jj) & 484 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 485 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 486 END_2D 487 ! 488 ENDIF 489 ! 490 ! !== after SSH ==! (jn+1) 491 ! 492 ! ! update (ua_e,va_e) to enforce volume conservation at open boundaries 493 ! ! values of zhup2_e and zhvp2_e on the halo are not needed in bdy_vol2d 797 494 IF( ln_bdy .AND. ln_vol ) CALL bdy_vol2d( kt, jn, ua_e, va_e, zhup2_e, zhvp2_e ) 798 ! 799 zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 800 zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 495 ! 496 ! ! resulting flux at mid-step (not over the full domain) 497 zhU(1:jpim1,1:jpj ) = e2u(1:jpim1,1:jpj ) * ua_e(1:jpim1,1:jpj ) * zhup2_e(1:jpim1,1:jpj ) ! not jpi-column 498 zhV(1:jpi ,1:jpjm1) = e1v(1:jpi ,1:jpjm1) * va_e(1:jpi ,1:jpjm1) * zhvp2_e(1:jpi ,1:jpjm1) ! not jpj-row 801 499 ! 802 500 #if defined key_agrif 803 501 ! Set fluxes during predictor step to ensure volume conservation 804 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 805 IF((nbondi == -1).OR.(nbondi == 2)) THEN 806 DO jj = 1, jpj 807 zwx(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * e2u(2:nbghostcells+1,jj) 808 zwy(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * e1v(2:nbghostcells+1,jj) 809 END DO 810 ENDIF 811 IF((nbondi == 1).OR.(nbondi == 2)) THEN 812 DO jj=1,jpj 813 zwx(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * e2u(nlci-nbghostcells-1:nlci-2,jj) 814 zwy(nlci-nbghostcells :nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * e1v(nlci-nbghostcells :nlci-1,jj) 815 END DO 816 ENDIF 817 IF((nbondj == -1).OR.(nbondj == 2)) THEN 818 DO ji=1,jpi 819 zwy(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * e1v(ji,2:nbghostcells+1) 820 zwx(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * e2u(ji,2:nbghostcells+1) 821 END DO 822 ENDIF 823 IF((nbondj == 1).OR.(nbondj == 2)) THEN 824 DO ji=1,jpi 825 zwy(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * e1v(ji,nlcj-nbghostcells-1:nlcj-2) 826 zwx(ji,nlcj-nbghostcells :nlcj-1) = ubdy_n(ji,1:nbghostcells) * e2u(ji,nlcj-nbghostcells :nlcj-1) 827 END DO 828 ENDIF 829 ENDIF 502 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) CALL agrif_dyn_ts_flux( jn, zhU, zhV ) 830 503 #endif 831 IF( ln_wd_il ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 832 833 IF ( ln_wd_dl ) THEN 834 ! 835 ! un_e and vn_e are set to zero at faces where the direction of the flow is from dry cells 836 ! 837 DO jj = 1, jpjm1 838 DO ji = 1, jpim1 839 IF ( zwx(ji,jj) > 0.0 ) THEN 840 zuwdmask(ji, jj) = ztwdmask(ji ,jj) 841 ELSE 842 zuwdmask(ji, jj) = ztwdmask(ji+1,jj) 843 END IF 844 zwx(ji, jj) = zuwdmask(ji,jj)*zwx(ji, jj) 845 un_e(ji,jj) = zuwdmask(ji,jj)*un_e(ji,jj) 846 847 IF ( zwy(ji,jj) > 0.0 ) THEN 848 zvwdmask(ji, jj) = ztwdmask(ji, jj ) 849 ELSE 850 zvwdmask(ji, jj) = ztwdmask(ji, jj+1) 851 END IF 852 zwy(ji, jj) = zvwdmask(ji,jj)*zwy(ji,jj) 853 vn_e(ji,jj) = zvwdmask(ji,jj)*vn_e(ji,jj) 854 END DO 855 END DO 504 IF( ln_wd_il ) CALL wad_lmt_bt(zhU, zhV, sshn_e, zssh_frc, rDt_e) !!gm wad_lmt_bt use of lbc_lnk on zhU, zhV 505 506 IF( ln_wd_dl ) THEN ! un_e and vn_e are set to zero at faces where 507 ! ! the direction of the flow is from dry cells 508 CALL wad_Umsk( ztwdmask, zhU, zhV, un_e, vn_e, zuwdmask, zvwdmask ) ! not jpi colomn for U, not jpj row for V 856 509 ! 857 510 ENDIF 858 859 ! Sum over sub-time-steps to compute advective velocities 860 za2 = wgtbtp2(jn) 861 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 862 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 863 864 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc = True) 865 IF ( ln_wd_dl_bc ) THEN 866 zuwdav2(:,:) = zuwdav2(:,:) + za2 * zuwdmask(:,:) 867 zvwdav2(:,:) = zvwdav2(:,:) + za2 * zvwdmask(:,:) 868 END IF 869 870 ! Set next sea level: 871 DO jj = 2, jpjm1 872 DO ji = fs_2, fs_jpim1 ! vector opt. 873 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 874 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 875 END DO 876 END DO 877 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 878 879 CALL lbc_lnk( 'dynspg_ts', ssha_e, 'T', 1._wp ) 880 511 ! 512 ! 513 ! Compute Sea Level at step jit+1 514 !-- m+1 m m+1/2 --! 515 !-- ssh = ssh - delta_t' * [ frc + div( flux ) ] --! 516 !-------------------------------------------------------------------------! 517 DO_2D( 0, 0, 0, 0 ) 518 zhdiv = ( zhU(ji,jj) - zhU(ji-1,jj) + zhV(ji,jj) - zhV(ji,jj-1) ) * r1_e1e2t(ji,jj) 519 ssha_e(ji,jj) = ( sshn_e(ji,jj) - rDt_e * ( zssh_frc(ji,jj) + zhdiv ) ) * ssmask(ji,jj) 520 END_2D 521 ! 522 CALL lbc_lnk_multi( 'dynspg_ts', ssha_e, 'T', 1._wp, zhU, 'U', -1._wp, zhV, 'V', -1._wp ) 523 ! 881 524 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 882 525 IF( ln_bdy ) CALL bdy_ssh( ssha_e ) … … 884 527 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 885 528 #endif 529 ! 530 ! ! Sum over sub-time-steps to compute advective velocities 531 za2 = wgtbtp2(jn) ! zhU, zhV hold fluxes extrapolated at jn+0.5 532 un_adv(:,:) = un_adv(:,:) + za2 * zhU(:,:) * r1_e2u(:,:) 533 vn_adv(:,:) = vn_adv(:,:) + za2 * zhV(:,:) * r1_e1v(:,:) 534 ! sum over sub-time-steps to decide which baroclinic velocities to set to zero (zuwdav2 is only used when ln_wd_dl_bc=True) 535 IF ( ln_wd_dl_bc ) THEN 536 zuwdav2(1:jpim1,1:jpj ) = zuwdav2(1:jpim1,1:jpj ) + za2 * zuwdmask(1:jpim1,1:jpj ) ! not jpi-column 537 zvwdav2(1:jpi ,1:jpjm1) = zvwdav2(1:jpi ,1:jpjm1) + za2 * zvwdmask(1:jpi ,1:jpjm1) ! not jpj-row 538 END IF 539 ! 886 540 ! 887 541 ! Sea Surface Height at u-,v-points (vvl case only) 888 542 IF( .NOT.ln_linssh ) THEN 889 DO jj = 2, jpjm1 890 DO ji = 2, jpim1 ! NO Vector Opt. 891 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 892 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 893 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 894 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 895 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 896 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 897 END DO 898 END DO 899 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) 543 DO_2D( 0, 0, 0, 0 ) 544 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 545 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 546 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 547 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 548 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 549 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 550 END_2D 900 551 ENDIF 901 ! 902 ! Half-step back interpolation of SSH for surface pressure computation: 903 !---------------------------------------------------------------------- 904 IF ((jn==1).AND.ll_init) THEN 905 za0=1._wp ! Forward-backward 906 za1=0._wp 907 za2=0._wp 908 za3=0._wp 909 ELSEIF ((jn==2).AND.ll_init) THEN ! AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 910 za0= 1.0833333333333_wp ! za0 = 1-gam-eps 911 za1=-0.1666666666666_wp ! za1 = gam 912 za2= 0.0833333333333_wp ! za2 = eps 913 za3= 0._wp 914 ELSE ! AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 915 IF (rn_bt_alpha==0._wp) THEN 916 za0=0.614_wp ! za0 = 1/2 + gam + 2*eps 917 za1=0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 918 za2=0.088_wp ! za2 = gam 919 za3=0.013_wp ! za3 = eps 920 ELSE 921 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 922 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 923 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 924 za1 = 1._wp - za0 - zgamma - zepsilon 925 za2 = zgamma 926 za3 = zepsilon 927 ENDIF 928 ENDIF 929 ! 552 ! 553 ! Half-step back interpolation of SSH for surface pressure computation at step jit+1/2 554 !-- m+1/2 m+1 m m-1 m-2 --! 555 !-- ssh' = za0 * ssh + za1 * ssh + za2 * ssh + za3 * ssh --! 556 !------------------------------------------------------------------------------------------! 557 CALL ts_bck_interp( jn, ll_init, za0, za1, za2, za3 ) ! coeficients of the interpolation 930 558 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 931 559 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 932 933 IF( ln_wd_il ) THEN ! Calculating and applying W/D gravity filters 934 DO jj = 2, jpjm1 935 DO ji = 2, jpim1 936 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 937 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 938 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) ) & 939 & > rn_wdmin1 + rn_wdmin2 940 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji+1,jj)) > 1.E-12 ).AND.( & 941 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji+1,jj) ) > & 942 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 943 944 IF(ll_tmp1) THEN 945 zcpx(ji,jj) = 1.0_wp 946 ELSE IF(ll_tmp2) THEN 947 ! no worries about zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj) = 0, it won't happen ! here 948 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + ht_0(ji+1,jj) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 949 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji ,jj)) ) 950 ELSE 951 zcpx(ji,jj) = 0._wp 952 ENDIF 953 ! 954 ll_tmp1 = MIN( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 955 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 956 & MAX( zsshp2_e(ji,jj) + ht_0(ji,jj) , zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) ) & 957 & > rn_wdmin1 + rn_wdmin2 958 ll_tmp2 = (ABS(zsshp2_e(ji,jj) - zsshp2_e(ji,jj+1)) > 1.E-12 ).AND.( & 959 & MAX( zsshp2_e(ji,jj) , zsshp2_e(ji,jj+1) ) > & 960 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 961 962 IF(ll_tmp1) THEN 963 zcpy(ji,jj) = 1.0_wp 964 ELSEIF(ll_tmp2) THEN 965 ! no worries about zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj ) = 0, it won't happen ! here 966 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + ht_0(ji,jj+1) - zsshp2_e(ji,jj) - ht_0(ji,jj)) & 967 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj )) ) 968 ELSE 969 zcpy(ji,jj) = 0._wp 970 ENDIF 971 END DO 972 END DO 973 ENDIF 974 ! 975 ! Compute associated depths at U and V points: 976 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 977 ! 978 DO jj = 2, jpjm1 979 DO ji = 2, jpim1 980 zx1 = r1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 981 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 982 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 983 zy1 = r1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 984 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 985 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 986 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 987 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 988 END DO 989 END DO 990 ! 560 ! 561 ! ! Surface pressure gradient 562 zldg = ( 1._wp - rn_scal_load ) * grav ! local factor 563 DO_2D( 0, 0, 0, 0 ) 564 zu_spg(ji,jj) = - zldg * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 565 zv_spg(ji,jj) = - zldg * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 566 END_2D 567 IF( ln_wd_il ) THEN ! W/D : gravity filters applied on pressure gradient 568 CALL wad_spg( zsshp2_e, zcpx, zcpy ) ! Calculating W/D gravity filters 569 zu_spg(2:jpim1,2:jpjm1) = zu_spg(2:jpim1,2:jpjm1) * zcpx(2:jpim1,2:jpjm1) 570 zv_spg(2:jpim1,2:jpjm1) = zv_spg(2:jpim1,2:jpjm1) * zcpy(2:jpim1,2:jpjm1) 991 571 ENDIF 992 572 ! … … 994 574 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 995 575 ! at each time step. We however keep them constant here for optimization. 996 ! Recall that zwx and zwy arrays hold fluxes at this stage: 997 ! zwx(:,:) = e2u(:,:) * ua_e(:,:) * zhup2_e(:,:) ! fluxes at jn+0.5 998 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 999 ! 1000 SELECT CASE( nvor_scheme ) 1001 CASE( np_ENT ) ! energy conserving scheme (t-point) 1002 DO jj = 2, jpjm1 1003 DO ji = 2, jpim1 ! vector opt. 1004 1005 z1_hu = ssumask(ji,jj) / ( zhup2_e(ji,jj) + 1._wp - ssumask(ji,jj) ) 1006 z1_hv = ssvmask(ji,jj) / ( zhvp2_e(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1007 1008 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1009 & * ( e1e2t(ji+1,jj)*zhtp2_e(ji+1,jj)*ff_t(ji+1,jj) * ( va_e(ji+1,jj) + va_e(ji+1,jj-1) ) & 1010 & + e1e2t(ji ,jj)*zhtp2_e(ji ,jj)*ff_t(ji ,jj) * ( va_e(ji ,jj) + va_e(ji ,jj-1) ) ) 1011 ! 1012 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1013 & * ( e1e2t(ji,jj+1)*zhtp2_e(ji,jj+1)*ff_t(ji,jj+1) * ( ua_e(ji,jj+1) + ua_e(ji-1,jj+1) ) & 1014 & + e1e2t(ji,jj )*zhtp2_e(ji,jj )*ff_t(ji,jj ) * ( ua_e(ji,jj ) + ua_e(ji-1,jj ) ) ) 1015 END DO 1016 END DO 1017 ! 1018 CASE( np_ENE, np_MIX ) ! energy conserving scheme (f-point) 1019 DO jj = 2, jpjm1 1020 DO ji = fs_2, fs_jpim1 ! vector opt. 1021 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 1022 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1023 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 1024 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1025 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1026 zv_trd(ji,jj) =-r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1027 END DO 1028 END DO 1029 ! 1030 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1031 DO jj = 2, jpjm1 1032 DO ji = fs_2, fs_jpim1 ! vector opt. 1033 zy1 = r1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 1034 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 1035 zx1 = - r1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 1036 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 1037 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1038 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1039 END DO 1040 END DO 1041 ! 1042 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1043 DO jj = 2, jpjm1 1044 DO ji = fs_2, fs_jpim1 ! vector opt. 1045 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 1046 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 1047 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 1048 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 1049 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 1050 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 1051 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 1052 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 1053 END DO 1054 END DO 1055 ! 1056 END SELECT 576 ! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated) 577 CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd ) 1057 578 ! 1058 579 ! Add tidal astronomical forcing if defined 1059 580 IF ( ln_tide .AND. ln_tide_pot ) THEN 1060 DO jj = 2, jpjm1 1061 DO ji = fs_2, fs_jpim1 ! vector opt. 1062 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 1063 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 1064 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 1065 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg 1066 END DO 1067 END DO 581 DO_2D( 0, 0, 0, 0 ) 582 zu_trd(ji,jj) = zu_trd(ji,jj) + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 583 zv_trd(ji,jj) = zv_trd(ji,jj) + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 584 END_2D 1068 585 ENDIF 1069 586 ! … … 1071 588 !jth do implicitly instead 1072 589 IF ( .NOT. ll_wd ) THEN ! Revert to explicit for bit comparison tests in non wad runs 1073 DO jj = 2, jpjm1 1074 DO ji = fs_2, fs_jpim1 ! vector opt. 1075 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 1076 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 1077 END DO 1078 END DO 1079 ENDIF 1080 ! 1081 ! Surface pressure trend: 1082 IF( ln_wd_il ) THEN 1083 DO jj = 2, jpjm1 1084 DO ji = 2, jpim1 1085 ! Add surface pressure gradient 1086 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1087 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1088 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg * zcpx(ji,jj) 1089 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg * zcpy(ji,jj) 1090 END DO 1091 END DO 1092 ELSE 1093 DO jj = 2, jpjm1 1094 DO ji = fs_2, fs_jpim1 ! vector opt. 1095 ! Add surface pressure gradient 1096 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 1097 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 1098 zwx(ji,jj) = (1._wp - rn_scal_load) * zu_spg 1099 zwy(ji,jj) = (1._wp - rn_scal_load) * zv_spg 1100 END DO 1101 END DO 1102 END IF 1103 590 DO_2D( 0, 0, 0, 0 ) 591 zu_trd(ji,jj) = zu_trd(ji,jj) + zCdU_u(ji,jj) * un_e(ji,jj) * hur_e(ji,jj) 592 zv_trd(ji,jj) = zv_trd(ji,jj) + zCdU_v(ji,jj) * vn_e(ji,jj) * hvr_e(ji,jj) 593 END_2D 594 ENDIF 1104 595 ! 1105 596 ! Set next velocities: 597 ! Compute barotropic speeds at step jit+1 (h : total height of the water colomn) 598 !-- VECTOR FORM 599 !-- m+1 m / m+1/2 \ --! 600 !-- u = u + delta_t' * \ (1-r)*g * grad_x( ssh') - f * k vect u + frc / --! 601 !-- --! 602 !-- FLUX FORM --! 603 !-- m+1 __1__ / m m / m+1/2 m+1/2 m+1/2 n \ \ --! 604 !-- u = m+1 | h * u + delta_t' * \ h * (1-r)*g * grad_x( ssh') - h * f * k vect u + h * frc / | --! 605 !-- h \ / --! 606 !------------------------------------------------------------------------------------------------------------------------! 1106 607 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 1107 DO jj = 2, jpjm1 1108 DO ji = fs_2, fs_jpim1 ! vector opt. 1109 ua_e(ji,jj) = ( un_e(ji,jj) & 1110 & + rdtbt * ( zwx(ji,jj) & 1111 & + zu_trd(ji,jj) & 1112 & + zu_frc(ji,jj) ) & 1113 & ) * ssumask(ji,jj) 1114 1115 va_e(ji,jj) = ( vn_e(ji,jj) & 1116 & + rdtbt * ( zwy(ji,jj) & 1117 & + zv_trd(ji,jj) & 1118 & + zv_frc(ji,jj) ) & 1119 & ) * ssvmask(ji,jj) 1120 1121 END DO 1122 END DO 608 DO_2D( 0, 0, 0, 0 ) 609 ua_e(ji,jj) = ( un_e(ji,jj) & 610 & + rDt_e * ( zu_spg(ji,jj) & 611 & + zu_trd(ji,jj) & 612 & + zu_frc(ji,jj) ) & 613 & ) * ssumask(ji,jj) 614 615 va_e(ji,jj) = ( vn_e(ji,jj) & 616 & + rDt_e * ( zv_spg(ji,jj) & 617 & + zv_trd(ji,jj) & 618 & + zv_frc(ji,jj) ) & 619 & ) * ssvmask(ji,jj) 620 END_2D 1123 621 ! 1124 622 ELSE !* Flux form 1125 DO jj = 2, jpjm1 1126 DO ji = fs_2, fs_jpim1 ! vector opt. 1127 1128 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 1129 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 1130 1131 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 1132 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 1133 1134 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 1135 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 1136 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 1137 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 1138 & ) * zhura 1139 1140 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 1141 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 1142 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 1143 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 1144 & ) * zhvra 1145 END DO 1146 END DO 623 DO_2D( 0, 0, 0, 0 ) 624 ! ! hu_e, hv_e hold depth at jn, zhup2_e, zhvp2_e hold extrapolated depth at jn+1/2 625 ! ! backward interpolated depth used in spg terms at jn+1/2 626 zhu_bck = hu_0(ji,jj) + r1_2*r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 627 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) * ssumask(ji,jj) 628 zhv_bck = hv_0(ji,jj) + r1_2*r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 629 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) * ssvmask(ji,jj) 630 ! ! inverse depth at jn+1 631 z1_hu = ssumask(ji,jj) / ( hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - ssumask(ji,jj) ) 632 z1_hv = ssvmask(ji,jj) / ( hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - ssvmask(ji,jj) ) 633 ! 634 ua_e(ji,jj) = ( hu_e (ji,jj) * un_e (ji,jj) & 635 & + rDt_e * ( zhu_bck * zu_spg (ji,jj) & ! 636 & + zhup2_e(ji,jj) * zu_trd (ji,jj) & ! 637 & + hu(ji,jj,Kmm) * zu_frc (ji,jj) ) ) * z1_hu 638 ! 639 va_e(ji,jj) = ( hv_e (ji,jj) * vn_e (ji,jj) & 640 & + rDt_e * ( zhv_bck * zv_spg (ji,jj) & ! 641 & + zhvp2_e(ji,jj) * zv_trd (ji,jj) & ! 642 & + hv(ji,jj,Kmm) * zv_frc (ji,jj) ) ) * z1_hv 643 END_2D 1147 644 ENDIF 1148 645 !jth implicit bottom friction: 1149 646 IF ( ll_wd ) THEN ! revert to explicit for bit comparison tests in non wad runs 1150 DO jj = 2, jpjm1 1151 DO ji = fs_2, fs_jpim1 ! vector opt. 1152 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rdtbt * zCdU_u(ji,jj) * hur_e(ji,jj)) 1153 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rdtbt * zCdU_v(ji,jj) * hvr_e(ji,jj)) 1154 END DO 1155 END DO 1156 ENDIF 1157 1158 1159 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 1160 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 1161 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 1162 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 1163 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 1164 ! 1165 ENDIF 1166 ! !* domain lateral boundary 1167 CALL lbc_lnk_multi( 'dynspg_ts', ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 1168 ! 647 DO_2D( 0, 0, 0, 0 ) 648 ua_e(ji,jj) = ua_e(ji,jj) /(1.0 - rDt_e * zCdU_u(ji,jj) * hur_e(ji,jj)) 649 va_e(ji,jj) = va_e(ji,jj) /(1.0 - rDt_e * zCdU_v(ji,jj) * hvr_e(ji,jj)) 650 END_2D 651 ENDIF 652 653 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 654 hu_e (2:jpim1,2:jpjm1) = hu_0(2:jpim1,2:jpjm1) + zsshu_a(2:jpim1,2:jpjm1) 655 hur_e(2:jpim1,2:jpjm1) = ssumask(2:jpim1,2:jpjm1) / ( hu_e(2:jpim1,2:jpjm1) + 1._wp - ssumask(2:jpim1,2:jpjm1) ) 656 hv_e (2:jpim1,2:jpjm1) = hv_0(2:jpim1,2:jpjm1) + zsshv_a(2:jpim1,2:jpjm1) 657 hvr_e(2:jpim1,2:jpjm1) = ssvmask(2:jpim1,2:jpjm1) / ( hv_e(2:jpim1,2:jpjm1) + 1._wp - ssvmask(2:jpim1,2:jpjm1) ) 658 ENDIF 659 ! 660 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 661 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp & 662 & , hu_e , 'U', 1._wp, hv_e , 'V', 1._wp & 663 & , hur_e, 'U', 1._wp, hvr_e, 'V', 1._wp ) 664 ELSE 665 CALL lbc_lnk_multi( 'dynspg_ts', ua_e , 'U', -1._wp, va_e , 'V', -1._wp ) 666 ENDIF 1169 667 ! ! open boundaries 1170 668 IF( ln_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) … … 1190 688 za1 = wgtbtp1(jn) 1191 689 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 1192 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:)1193 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:)690 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) 691 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) 1194 692 ELSE ! Sum transports 1195 693 IF ( .NOT.ln_wd_dl ) THEN 1196 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:)1197 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:)694 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) 695 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) 1198 696 ELSE 1199 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)1200 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)697 puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:) 698 pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:) 1201 699 END IF 1202 700 ENDIF 1203 701 ! ! Sum sea level 1204 ssha(:,:) = ssha(:,:) + za1 * ssha_e(:,:)702 pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:) 1205 703 1206 704 ! ! ==================== ! … … 1213 711 ! Set advection velocity correction: 1214 712 IF (ln_bt_fw) THEN 1215 zwx(:,:) = un_adv(:,:) 1216 zwy(:,:) = vn_adv(:,:) 1217 IF( .NOT.( kt == nit000 .AND. neuler==0 ) ) THEN 1218 un_adv(:,:) = r1_2 * ( ub2_b(:,:) + zwx(:,:) - atfp * un_bf(:,:) ) 1219 vn_adv(:,:) = r1_2 * ( vb2_b(:,:) + zwy(:,:) - atfp * vn_bf(:,:) ) 1220 ! 1221 ! Update corrective fluxes for next time step: 1222 un_bf(:,:) = atfp * un_bf(:,:) + (zwx(:,:) - ub2_b(:,:)) 1223 vn_bf(:,:) = atfp * vn_bf(:,:) + (zwy(:,:) - vb2_b(:,:)) 713 IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN 714 DO_2D( 1, 1, 1, 1 ) 715 zun_save = un_adv(ji,jj) 716 zvn_save = vn_adv(ji,jj) 717 ! ! apply the previously computed correction 718 un_adv(ji,jj) = r1_2 * ( ub2_b(ji,jj) + zun_save - rn_atfp * un_bf(ji,jj) ) 719 vn_adv(ji,jj) = r1_2 * ( vb2_b(ji,jj) + zvn_save - rn_atfp * vn_bf(ji,jj) ) 720 ! ! Update corrective fluxes for next time step 721 un_bf(ji,jj) = rn_atfp * un_bf(ji,jj) + ( zun_save - ub2_b(ji,jj) ) 722 vn_bf(ji,jj) = rn_atfp * vn_bf(ji,jj) + ( zvn_save - vb2_b(ji,jj) ) 723 ! ! Save integrated transport for next computation 724 ub2_b(ji,jj) = zun_save 725 vb2_b(ji,jj) = zvn_save 726 END_2D 1224 727 ELSE 1225 un_bf(:,:) = 0._wp 1226 vn_bf(:,:) = 0._wp 1227 END IF 1228 ! Save integrated transport for next computation 1229 ub2_b(:,:) = zwx(:,:) 1230 vb2_b(:,:) = zwy(:,:) 728 un_bf(:,:) = 0._wp ! corrective fluxes for next time step set to zero 729 vn_bf(:,:) = 0._wp 730 ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation 731 vb2_b(:,:) = vn_adv(:,:) 732 END IF 1231 733 ENDIF 1232 734 … … 1236 738 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 1237 739 DO jk=1,jpkm1 1238 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * r1_2dt_b1239 va(:,:,jk) = va(:,:,jk) + ( va_b(:,:) - vb_b(:,:) ) * r1_2dt_b740 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b 741 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b 1240 742 END DO 1241 743 ELSE 1242 ! At this stage, ssha has been corrected: compute new depths at velocity points 1243 DO jj = 1, jpjm1 1244 DO ji = 1, jpim1 ! NO Vector Opt. 1245 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 1246 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1247 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1248 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 1249 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1250 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1251 END DO 1252 END DO 744 ! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points 745 DO_2D( 1, 0, 1, 0 ) 746 zsshu_a(ji,jj) = r1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 747 & * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) & 748 & + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) 749 zsshv_a(ji,jj) = r1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 750 & * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) & 751 & + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) 752 END_2D 1253 753 CALL lbc_lnk_multi( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1254 754 ! 1255 755 DO jk=1,jpkm1 1256 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * r1_2dt_b1257 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * r1_2dt_b756 puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b 757 pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b 1258 758 END DO 1259 759 ! Save barotropic velocities not transport: 1260 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )1261 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )760 puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 761 pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 1262 762 ENDIF 1263 763 … … 1265 765 ! Correct velocities so that the barotropic velocity equals (un_adv, vn_adv) (in all cases) 1266 766 DO jk = 1, jpkm1 1267 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:)*r1_hu_n(:,:) - un_b(:,:) ) * umask(:,:,jk)1268 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:)*r1_hv_n(:,:) - vn_b(:,:) ) * vmask(:,:,jk)767 puu(:,:,jk,Kmm) = ( puu(:,:,jk,Kmm) + un_adv(:,:)*r1_hu(:,:,Kmm) - puu_b(:,:,Kmm) ) * umask(:,:,jk) 768 pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk) 1269 769 END DO 1270 770 1271 771 IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN 1272 772 DO jk = 1, jpkm1 1273 un(:,:,jk) = ( un_adv(:,:)*r1_hu_n(:,:) &1274 & + zuwdav2(:,:)*( un(:,:,jk) - un_adv(:,:)*r1_hu_n(:,:)) ) * umask(:,:,jk)1275 vn(:,:,jk) = ( vn_adv(:,:)*r1_hv_n(:,:) &1276 & + zvwdav2(:,:)*( vn(:,:,jk) - vn_adv(:,:)*r1_hv_n(:,:)) ) * vmask(:,:,jk)773 puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) & 774 & + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk) 775 pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) & 776 & + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk) 1277 777 END DO 1278 778 END IF 1279 779 1280 780 1281 CALL iom_put( "ubar", un_adv(:,:)*r1_hu _n(:,:) ) ! barotropic i-current1282 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv _n(:,:) ) ! barotropic i-current781 CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current 782 CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current 1283 783 ! 1284 784 #if defined key_agrif … … 1303 803 IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 ) 1304 804 ! 1305 IF( ln_diatmb ) THEN 1306 CALL iom_put( "baro_u" , un_b*ssumask(:,:)+zmdi*(1.-ssumask(:,:) ) ) ! Barotropic U Velocity 1307 CALL iom_put( "baro_v" , vn_b*ssvmask(:,:)+zmdi*(1.-ssvmask(:,:) ) ) ! Barotropic V Velocity 1308 ENDIF 805 CALL iom_put( "baro_u" , puu_b(:,:,Kmm) ) ! Barotropic U Velocity 806 CALL iom_put( "baro_v" , pvv_b(:,:,Kmm) ) ! Barotropic V Velocity 1309 807 ! 1310 808 END SUBROUTINE dyn_spg_ts … … 1320 818 LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true. 1321 819 INTEGER, INTENT(inout) :: jpit ! cycle length 1322 REAL(wp), DIMENSION(3*nn_ baro), INTENT(inout) :: zwgt1, & ! Primary weights820 REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights 1323 821 zwgt2 ! Secondary weights 1324 822 … … 1332 830 ! Set time index when averaged value is requested 1333 831 IF (ll_fw) THEN 1334 jic = nn_ baro832 jic = nn_e 1335 833 ELSE 1336 jic = 2 * nn_ baro834 jic = 2 * nn_e 1337 835 ENDIF 1338 836 … … 1340 838 IF (ll_av) THEN 1341 839 ! Define simple boxcar window for primary weights 1342 ! (width = nn_ baro, centered around jic)840 ! (width = nn_e, centered around jic) 1343 841 SELECT CASE ( nn_bt_flt ) 1344 842 CASE( 0 ) ! No averaging … … 1346 844 jpit = jic 1347 845 1348 CASE( 1 ) ! Boxcar, width = nn_ baro1349 DO jn = 1, 3*nn_ baro1350 za1 = ABS(float(jn-jic))/float(nn_ baro)846 CASE( 1 ) ! Boxcar, width = nn_e 847 DO jn = 1, 3*nn_e 848 za1 = ABS(float(jn-jic))/float(nn_e) 1351 849 IF (za1 < 0.5_wp) THEN 1352 850 zwgt1(jn) = 1._wp … … 1355 853 ENDDO 1356 854 1357 CASE( 2 ) ! Boxcar, width = 2 * nn_ baro1358 DO jn = 1, 3*nn_ baro1359 za1 = ABS(float(jn-jic))/float(nn_ baro)855 CASE( 2 ) ! Boxcar, width = 2 * nn_e 856 DO jn = 1, 3*nn_e 857 za1 = ABS(float(jn-jic))/float(nn_e) 1360 858 IF (za1 < 1._wp) THEN 1361 859 zwgt1(jn) = 1._wp … … 1401 899 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 1402 900 ! ! --------------- 1403 IF( ln_rstart .AND. ln_bt_fw .AND. ( neuler/=0) ) THEN !* Read the restart file1404 IF(lrxios) CALL iom_swap( TRIM(crxios_context) )1405 CALL iom_get( numror, jpdom_auto glo, 'ub2_b' , ub2_b (:,:), ldxios = lrxios )1406 CALL iom_get( numror, jpdom_auto glo, 'vb2_b' , vb2_b (:,:), ldxios = lrxios )1407 CALL iom_get( numror, jpdom_auto glo, 'un_bf' , un_bf (:,:), ldxios = lrxios )1408 CALL iom_get( numror, jpdom_auto glo, 'vn_bf' , vn_bf (:,:), ldxios = lrxios )901 IF( ln_rstart .AND. ln_bt_fw .AND. (.NOT.l_1st_euler) ) THEN !* Read the restart file 902 IF(lrxios) CALL iom_swap( TRIM(crxios_context) 903 CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 904 CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 905 CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 906 CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 1409 907 IF( .NOT.ln_bt_av ) THEN 1410 CALL iom_get( numror, jpdom_auto glo, 'sshbb_e' , sshbb_e(:,:), ldxios = lrxios )1411 CALL iom_get( numror, jpdom_auto glo, 'ubb_e' , ubb_e(:,:), ldxios = lrxios )1412 CALL iom_get( numror, jpdom_auto glo, 'vbb_e' , vbb_e(:,:), ldxios = lrxios )1413 CALL iom_get( numror, jpdom_auto glo, 'sshb_e' , sshb_e(:,:), ldxios = lrxios )1414 CALL iom_get( numror, jpdom_auto glo, 'ub_e' , ub_e(:,:), ldxios = lrxios )1415 CALL iom_get( numror, jpdom_auto glo, 'vb_e' , vb_e(:,:), ldxios = lrxios )908 CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios ) 909 CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 910 CALL iom_get( numror, jpdom_auto, 'vbb_e' , vbb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 911 CALL iom_get( numror, jpdom_auto, 'sshb_e' , sshb_e(:,:), cd_type = 'T', psgn = 1._wp, ldxios = lrxios ) 912 CALL iom_get( numror, jpdom_auto, 'ub_e' , ub_e(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 913 CALL iom_get( numror, jpdom_auto, 'vb_e' , vb_e(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 1416 914 ENDIF 1417 915 #if defined key_agrif 1418 916 ! Read time integrated fluxes 1419 917 IF ( .NOT.Agrif_Root() ) THEN 1420 CALL iom_get( numror, jpdom_auto glo, 'ub2_i_b' , ub2_i_b(:,:), ldxios = lrxios )1421 CALL iom_get( numror, jpdom_auto glo, 'vb2_i_b' , vb2_i_b(:,:), ldxios = lrxios )918 CALL iom_get( numror, jpdom_auto, 'ub2_i_b' , ub2_i_b(:,:), cd_type = 'U', psgn = -1._wp, ldxios = lrxios ) 919 CALL iom_get( numror, jpdom_auto, 'vb2_i_b' , vb2_i_b(:,:), cd_type = 'V', psgn = -1._wp, ldxios = lrxios ) 1422 920 ENDIF 1423 921 #endif … … 1479 977 ! Max courant number for ext. grav. waves 1480 978 ! 1481 DO jj = 1, jpj 1482 DO ji =1, jpi 1483 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1484 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1485 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 1486 END DO 1487 END DO 1488 ! 1489 zcmax = MAXVAL( zcu(:,:) ) 979 DO_2D( 0, 0, 0, 0 ) 980 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 981 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 982 zcu(ji,jj) = SQRT( grav * MAX(ht_0(ji,jj),0._wp) * (zxr2 + zyr2) ) 983 END_2D 984 ! 985 zcmax = MAXVAL( zcu(Nis0:Nie0,Njs0:Nje0) ) 1490 986 CALL mpp_max( 'dynspg_ts', zcmax ) 1491 987 1492 988 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1493 IF( ln_bt_auto ) nn_ baro = CEILING( rdt / rn_bt_cmax * zcmax)989 IF( ln_bt_auto ) nn_e = CEILING( rn_Dt / rn_bt_cmax * zcmax) 1494 990 1495 r dtbt = rdt / REAL( nn_baro, wp )1496 zcmax = zcmax * r dtbt991 rDt_e = rn_Dt / REAL( nn_e , wp ) 992 zcmax = zcmax * rDt_e 1497 993 ! Print results 1498 994 IF(lwp) WRITE(numout,*) … … 1500 996 IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~' 1501 997 IF( ln_bt_auto ) THEN 1502 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_ baro'998 IF(lwp) WRITE(numout,*) ' ln_ts_auto =.true. Automatically set nn_e ' 1503 999 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1504 1000 ELSE 1505 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_ baro in namelist nn_baro = ', nn_baro1001 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e 1506 1002 ENDIF 1507 1003 1508 1004 IF(ln_bt_av) THEN 1509 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_ barotime steps is on '1005 IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on ' 1510 1006 ELSE 1511 1007 IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables ' … … 1527 1023 SELECT CASE ( nn_bt_flt ) 1528 1024 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1529 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_ baro'1530 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_ baro'1025 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_e' 1026 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_e' 1531 1027 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1, or 2' ) 1532 1028 END SELECT 1533 1029 ! 1534 1030 IF(lwp) WRITE(numout,*) ' ' 1535 IF(lwp) WRITE(numout,*) ' nn_ baro = ', nn_baro1536 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', r dtbt1031 IF(lwp) WRITE(numout,*) ' nn_e = ', nn_e 1032 IF(lwp) WRITE(numout,*) ' Barotropic time step [s] is :', rDt_e 1537 1033 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1538 1034 ! … … 1546 1042 ENDIF 1547 1043 IF( zcmax>0.9_wp ) THEN 1548 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_ baro!' )1044 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_e !' ) 1549 1045 ENDIF 1550 1046 ! … … 1581 1077 END SUBROUTINE dyn_spg_ts_init 1582 1078 1079 1080 SUBROUTINE dyn_cor_2D_init( Kmm ) 1081 !!--------------------------------------------------------------------- 1082 !! *** ROUTINE dyn_cor_2D_init *** 1083 !! 1084 !! ** Purpose : Set time splitting options 1085 !! Set arrays to remove/compute coriolis trend. 1086 !! Do it once during initialization if volume is fixed, else at each long time step. 1087 !! Note that these arrays are also used during barotropic loop. These are however frozen 1088 !! although they should be updated in the variable volume case. Not a big approximation. 1089 !! To remove this approximation, copy lines below inside barotropic loop 1090 !! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 1091 !! 1092 !! Compute zwz = f / ( height of the water colomn ) 1093 !!---------------------------------------------------------------------- 1094 INTEGER, INTENT(in) :: Kmm ! Time index 1095 INTEGER :: ji ,jj, jk ! dummy loop indices 1096 REAL(wp) :: z1_ht 1097 REAL(wp), DIMENSION(jpi,jpj) :: zhf 1098 !!---------------------------------------------------------------------- 1099 ! 1100 SELECT CASE( nvor_scheme ) 1101 CASE( np_EEN ) != EEN scheme using e3f energy & enstrophy scheme 1102 SELECT CASE( nn_een_e3f ) !* ff_f/e3 at F-point 1103 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 1104 DO_2D( 1, 0, 1, 0 ) 1105 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 1106 & ht(ji ,jj ) + ht(ji+1,jj ) ) * 0.25_wp 1107 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1108 END_2D 1109 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 1110 DO_2D( 1, 0, 1, 0 ) 1111 zwz(ji,jj) = ( ht (ji ,jj+1) + ht (ji+1,jj+1) & 1112 & + ht (ji ,jj ) + ht (ji+1,jj ) ) & 1113 & / ( MAX( 1._wp, ssmask(ji ,jj+1) + ssmask(ji+1,jj+1) & 1114 & + ssmask(ji ,jj ) + ssmask(ji+1,jj ) ) ) 1115 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff_f(ji,jj) / zwz(ji,jj) 1116 END_2D 1117 END SELECT 1118 CALL lbc_lnk( 'dynspg_ts', zwz, 'F', 1._wp ) 1119 ! 1120 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1121 DO_2D( 0, 1, 0, 1 ) 1122 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 1123 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 1124 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 1125 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 1126 END_2D 1127 ! 1128 CASE( np_EET ) != EEN scheme using e3t energy conserving scheme 1129 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 1130 DO_2D( 0, 1, 0, 1 ) 1131 z1_ht = ssmask(ji,jj) / ( ht(ji,jj) + 1._wp - ssmask(ji,jj) ) 1132 ftne(ji,jj) = ( ff_f(ji-1,jj ) + ff_f(ji ,jj ) + ff_f(ji ,jj-1) ) * z1_ht 1133 ftnw(ji,jj) = ( ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) + ff_f(ji ,jj ) ) * z1_ht 1134 ftse(ji,jj) = ( ff_f(ji ,jj ) + ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) ) * z1_ht 1135 ftsw(ji,jj) = ( ff_f(ji ,jj-1) + ff_f(ji-1,jj-1) + ff_f(ji-1,jj ) ) * z1_ht 1136 END_2D 1137 ! 1138 CASE( np_ENE, np_ENS , np_MIX ) != all other schemes (ENE, ENS, MIX) except ENT ! 1139 ! 1140 zwz(:,:) = 0._wp 1141 zhf(:,:) = 0._wp 1142 1143 !!gm assume 0 in both cases (which is almost surely WRONG ! ) as hvatf has been removed 1144 !!gm A priori a better value should be something like : 1145 !!gm zhf(i,j) = masked sum of ht(i,j) , ht(i+1,j) , ht(i,j+1) , (i+1,j+1) 1146 !!gm divided by the sum of the corresponding mask 1147 !!gm 1148 !! 1149 IF( .NOT.ln_sco ) THEN 1150 1151 !!gm agree the JC comment : this should be done in a much clear way 1152 1153 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 1154 ! Set it to zero for the time being 1155 ! IF( rn_hmin < 0._wp ) THEN ; jk = - INT( rn_hmin ) ! from a nb of level 1156 ! ELSE ; jk = MINLOC( gdepw_0, mask = gdepw_0 > rn_hmin, dim = 1 ) ! from a depth 1157 ! ENDIF 1158 ! zhf(:,:) = gdepw_0(:,:,jk+1) 1159 ! 1160 ELSE 1161 ! 1162 !zhf(:,:) = hbatf(:,:) 1163 DO_2D( 1, 0, 1, 0 ) 1164 zhf(ji,jj) = ( ht_0 (ji,jj ) + ht_0 (ji+1,jj ) & 1165 & + ht_0 (ji,jj+1) + ht_0 (ji+1,jj+1) ) & 1166 & / MAX( ssmask(ji,jj ) + ssmask(ji+1,jj ) & 1167 & + ssmask(ji,jj+1) + ssmask(ji+1,jj+1) , 1._wp ) 1168 END_2D 1169 ENDIF 1170 ! 1171 DO jj = 1, jpjm1 1172 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 1173 END DO 1174 ! 1175 DO jk = 1, jpkm1 1176 DO jj = 1, jpjm1 1177 zhf(:,jj) = zhf(:,jj) + e3f(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 1178 END DO 1179 END DO 1180 CALL lbc_lnk( 'dynspg_ts', zhf, 'F', 1._wp ) 1181 ! JC: TBC. hf should be greater than 0 1182 DO_2D( 1, 1, 1, 1 ) 1183 IF( zhf(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zhf(ji,jj) 1184 END_2D 1185 zwz(:,:) = ff_f(:,:) * zwz(:,:) 1186 END SELECT 1187 1188 END SUBROUTINE dyn_cor_2d_init 1189 1190 1191 1192 SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd ) 1193 !!--------------------------------------------------------------------- 1194 !! *** ROUTINE dyn_cor_2d *** 1195 !! 1196 !! ** Purpose : Compute u and v coriolis trends 1197 !!---------------------------------------------------------------------- 1198 INTEGER :: ji ,jj ! dummy loop indices 1199 REAL(wp) :: zx1, zx2, zy1, zy2, z1_hu, z1_hv ! - - 1200 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pht, phu, phv, punb, pvnb, zhU, zhV 1201 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: zu_trd, zv_trd 1202 !!---------------------------------------------------------------------- 1203 SELECT CASE( nvor_scheme ) 1204 CASE( np_ENT ) ! enstrophy conserving scheme (f-point) 1205 DO_2D( 0, 0, 0, 0 ) 1206 z1_hu = ssumask(ji,jj) / ( phu(ji,jj) + 1._wp - ssumask(ji,jj) ) 1207 z1_hv = ssvmask(ji,jj) / ( phv(ji,jj) + 1._wp - ssvmask(ji,jj) ) 1208 zu_trd(ji,jj) = + r1_4 * r1_e1e2u(ji,jj) * z1_hu & 1209 & * ( e1e2t(ji+1,jj)*pht(ji+1,jj)*ff_t(ji+1,jj) * ( pvnb(ji+1,jj) + pvnb(ji+1,jj-1) ) & 1210 & + e1e2t(ji ,jj)*pht(ji ,jj)*ff_t(ji ,jj) * ( pvnb(ji ,jj) + pvnb(ji ,jj-1) ) ) 1211 ! 1212 zv_trd(ji,jj) = - r1_4 * r1_e1e2v(ji,jj) * z1_hv & 1213 & * ( e1e2t(ji,jj+1)*pht(ji,jj+1)*ff_t(ji,jj+1) * ( punb(ji,jj+1) + punb(ji-1,jj+1) ) & 1214 & + e1e2t(ji,jj )*pht(ji,jj )*ff_t(ji,jj ) * ( punb(ji,jj ) + punb(ji-1,jj ) ) ) 1215 END_2D 1216 ! 1217 CASE( np_ENE , np_MIX ) ! energy conserving scheme (t-point) ENE or MIX 1218 DO_2D( 0, 0, 0, 0 ) 1219 zy1 = ( zhV(ji,jj-1) + zhV(ji+1,jj-1) ) * r1_e1u(ji,jj) 1220 zy2 = ( zhV(ji,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1221 zx1 = ( zhU(ji-1,jj) + zhU(ji-1,jj+1) ) * r1_e2v(ji,jj) 1222 zx2 = ( zhU(ji ,jj) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1223 ! energy conserving formulation for planetary vorticity term 1224 zu_trd(ji,jj) = r1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 1225 zv_trd(ji,jj) = - r1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) 1226 END_2D 1227 ! 1228 CASE( np_ENS ) ! enstrophy conserving scheme (f-point) 1229 DO_2D( 0, 0, 0, 0 ) 1230 zy1 = r1_8 * ( zhV(ji ,jj-1) + zhV(ji+1,jj-1) & 1231 & + zhV(ji ,jj ) + zhV(ji+1,jj ) ) * r1_e1u(ji,jj) 1232 zx1 = - r1_8 * ( zhU(ji-1,jj ) + zhU(ji-1,jj+1) & 1233 & + zhU(ji ,jj ) + zhU(ji ,jj+1) ) * r1_e2v(ji,jj) 1234 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 1235 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) 1236 END_2D 1237 ! 1238 CASE( np_EET , np_EEN ) ! energy & enstrophy scheme (using e3t or e3f) 1239 DO_2D( 0, 0, 0, 0 ) 1240 zu_trd(ji,jj) = + r1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zhV(ji ,jj ) & 1241 & + ftnw(ji+1,jj) * zhV(ji+1,jj ) & 1242 & + ftse(ji,jj ) * zhV(ji ,jj-1) & 1243 & + ftsw(ji+1,jj) * zhV(ji+1,jj-1) ) 1244 zv_trd(ji,jj) = - r1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zhU(ji-1,jj+1) & 1245 & + ftse(ji,jj+1) * zhU(ji ,jj+1) & 1246 & + ftnw(ji,jj ) * zhU(ji-1,jj ) & 1247 & + ftne(ji,jj ) * zhU(ji ,jj ) ) 1248 END_2D 1249 ! 1250 END SELECT 1251 ! 1252 END SUBROUTINE dyn_cor_2D 1253 1254 1255 SUBROUTINE wad_tmsk( pssh, ptmsk ) 1256 !!---------------------------------------------------------------------- 1257 !! *** ROUTINE wad_lmt *** 1258 !! 1259 !! ** Purpose : set wetting & drying mask at tracer points 1260 !! for the current barotropic sub-step 1261 !! 1262 !! ** Method : ??? 1263 !! 1264 !! ** Action : ptmsk : wetting & drying t-mask 1265 !!---------------------------------------------------------------------- 1266 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pssh ! 1267 REAL(wp), DIMENSION(jpi,jpj), INTENT( out) :: ptmsk ! 1268 ! 1269 INTEGER :: ji, jj ! dummy loop indices 1270 !!---------------------------------------------------------------------- 1271 ! 1272 IF( ln_wd_dl_rmp ) THEN 1273 DO_2D( 1, 1, 1, 1 ) 1274 IF ( pssh(ji,jj) + ht_0(ji,jj) > 2._wp * rn_wdmin1 ) THEN 1275 ! IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin2 ) THEN 1276 ptmsk(ji,jj) = 1._wp 1277 ELSEIF( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN 1278 ptmsk(ji,jj) = TANH( 50._wp*( ( pssh(ji,jj) + ht_0(ji,jj) - rn_wdmin1 )*r_rn_wdmin1) ) 1279 ELSE 1280 ptmsk(ji,jj) = 0._wp 1281 ENDIF 1282 END_2D 1283 ELSE 1284 DO_2D( 1, 1, 1, 1 ) 1285 IF ( pssh(ji,jj) + ht_0(ji,jj) > rn_wdmin1 ) THEN ; ptmsk(ji,jj) = 1._wp 1286 ELSE ; ptmsk(ji,jj) = 0._wp 1287 ENDIF 1288 END_2D 1289 ENDIF 1290 ! 1291 END SUBROUTINE wad_tmsk 1292 1293 1294 SUBROUTINE wad_Umsk( pTmsk, phU, phV, pu, pv, pUmsk, pVmsk ) 1295 !!---------------------------------------------------------------------- 1296 !! *** ROUTINE wad_lmt *** 1297 !! 1298 !! ** Purpose : set wetting & drying mask at tracer points 1299 !! for the current barotropic sub-step 1300 !! 1301 !! ** Method : ??? 1302 !! 1303 !! ** Action : ptmsk : wetting & drying t-mask 1304 !!---------------------------------------------------------------------- 1305 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pTmsk ! W & D t-mask 1306 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: phU, phV, pu, pv ! ocean velocities and transports 1307 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: pUmsk, pVmsk ! W & D u- and v-mask 1308 ! 1309 INTEGER :: ji, jj ! dummy loop indices 1310 !!---------------------------------------------------------------------- 1311 ! 1312 DO_2D( 1, 1, 1, 0 ) 1313 IF ( phU(ji,jj) > 0._wp ) THEN ; pUmsk(ji,jj) = pTmsk(ji ,jj) 1314 ELSE ; pUmsk(ji,jj) = pTmsk(ji+1,jj) 1315 ENDIF 1316 phU(ji,jj) = pUmsk(ji,jj)*phU(ji,jj) 1317 pu (ji,jj) = pUmsk(ji,jj)*pu (ji,jj) 1318 END_2D 1319 ! 1320 DO_2D( 1, 0, 1, 1 ) 1321 IF ( phV(ji,jj) > 0._wp ) THEN ; pVmsk(ji,jj) = pTmsk(ji,jj ) 1322 ELSE ; pVmsk(ji,jj) = pTmsk(ji,jj+1) 1323 ENDIF 1324 phV(ji,jj) = pVmsk(ji,jj)*phV(ji,jj) 1325 pv (ji,jj) = pVmsk(ji,jj)*pv (ji,jj) 1326 END_2D 1327 ! 1328 END SUBROUTINE wad_Umsk 1329 1330 1331 SUBROUTINE wad_spg( pshn, zcpx, zcpy ) 1332 !!--------------------------------------------------------------------- 1333 !! *** ROUTINE wad_sp *** 1334 !! 1335 !! ** Purpose : 1336 !!---------------------------------------------------------------------- 1337 INTEGER :: ji ,jj ! dummy loop indices 1338 LOGICAL :: ll_tmp1, ll_tmp2 1339 REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn 1340 REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy 1341 !!---------------------------------------------------------------------- 1342 DO_2D( 0, 0, 0, 0 ) 1343 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1344 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. & 1345 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji+1,jj) + ht_0(ji+1,jj) ) & 1346 & > rn_wdmin1 + rn_wdmin2 1347 ll_tmp2 = ( ABS( pshn(ji+1,jj) - pshn(ji ,jj)) > 1.E-12 ).AND.( & 1348 & MAX( pshn(ji,jj) , pshn(ji+1,jj) ) > & 1349 & MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) + rn_wdmin1 + rn_wdmin2 ) 1350 IF(ll_tmp1) THEN 1351 zcpx(ji,jj) = 1.0_wp 1352 ELSEIF(ll_tmp2) THEN 1353 ! no worries about pshn(ji+1,jj) - pshn(ji ,jj) = 0, it won't happen ! here 1354 zcpx(ji,jj) = ABS( (pshn(ji+1,jj) + ht_0(ji+1,jj) - pshn(ji,jj) - ht_0(ji,jj)) & 1355 & / (pshn(ji+1,jj) - pshn(ji ,jj)) ) 1356 zcpx(ji,jj) = max(min( zcpx(ji,jj) , 1.0_wp),0.0_wp) 1357 ELSE 1358 zcpx(ji,jj) = 0._wp 1359 ENDIF 1360 ! 1361 ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1362 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) .AND. & 1363 & MAX( pshn(ji,jj) + ht_0(ji,jj) , pshn(ji,jj+1) + ht_0(ji,jj+1) ) & 1364 & > rn_wdmin1 + rn_wdmin2 1365 ll_tmp2 = ( ABS( pshn(ji,jj) - pshn(ji,jj+1)) > 1.E-12 ).AND.( & 1366 & MAX( pshn(ji,jj) , pshn(ji,jj+1) ) > & 1367 & MAX( -ht_0(ji,jj) , -ht_0(ji,jj+1) ) + rn_wdmin1 + rn_wdmin2 ) 1368 1369 IF(ll_tmp1) THEN 1370 zcpy(ji,jj) = 1.0_wp 1371 ELSE IF(ll_tmp2) THEN 1372 ! no worries about pshn(ji,jj+1) - pshn(ji,jj ) = 0, it won't happen ! here 1373 zcpy(ji,jj) = ABS( (pshn(ji,jj+1) + ht_0(ji,jj+1) - pshn(ji,jj) - ht_0(ji,jj)) & 1374 & / (pshn(ji,jj+1) - pshn(ji,jj )) ) 1375 zcpy(ji,jj) = MAX( 0._wp , MIN( zcpy(ji,jj) , 1.0_wp ) ) 1376 ELSE 1377 zcpy(ji,jj) = 0._wp 1378 ENDIF 1379 END_2D 1380 1381 END SUBROUTINE wad_spg 1382 1383 1384 1385 SUBROUTINE dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b ,pvv_b, pu_RHSi, pv_RHSi, pCdU_u, pCdU_v ) 1386 !!---------------------------------------------------------------------- 1387 !! *** ROUTINE dyn_drg_init *** 1388 !! 1389 !! ** Purpose : - add the baroclinic top/bottom drag contribution to 1390 !! the baroclinic part of the barotropic RHS 1391 !! - compute the barotropic drag coefficients 1392 !! 1393 !! ** Method : computation done over the INNER domain only 1394 !!---------------------------------------------------------------------- 1395 INTEGER , INTENT(in ) :: Kbb, Kmm ! ocean time level indices 1396 REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(in ) :: puu, pvv ! ocean velocities and RHS of momentum equation 1397 REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(in ) :: puu_b, pvv_b ! barotropic velocities at main time levels 1398 REAL(wp), DIMENSION(jpi,jpj) , INTENT(inout) :: pu_RHSi, pv_RHSi ! baroclinic part of the barotropic RHS 1399 REAL(wp), DIMENSION(jpi,jpj) , INTENT( out) :: pCdU_u , pCdU_v ! barotropic drag coefficients 1400 ! 1401 INTEGER :: ji, jj ! dummy loop indices 1402 INTEGER :: ikbu, ikbv, iktu, iktv 1403 REAL(wp) :: zztmp 1404 REAL(wp), DIMENSION(jpi,jpj) :: zu_i, zv_i 1405 !!---------------------------------------------------------------------- 1406 ! 1407 ! !== Set the barotropic drag coef. ==! 1408 ! 1409 IF( ln_isfcav ) THEN ! top+bottom friction (ocean cavities) 1410 1411 DO_2D( 0, 0, 0, 0 ) 1412 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) + rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) 1413 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) + rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) 1414 END_2D 1415 ELSE ! bottom friction only 1416 DO_2D( 0, 0, 0, 0 ) 1417 pCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) 1418 pCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) 1419 END_2D 1420 ENDIF 1421 ! 1422 ! !== BOTTOM stress contribution from baroclinic velocities ==! 1423 ! 1424 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities 1425 1426 DO_2D( 0, 0, 0, 0 ) 1427 ikbu = mbku(ji,jj) 1428 ikbv = mbkv(ji,jj) 1429 zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm) 1430 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm) 1431 END_2D 1432 ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities 1433 1434 DO_2D( 0, 0, 0, 0 ) 1435 ikbu = mbku(ji,jj) 1436 ikbv = mbkv(ji,jj) 1437 zu_i(ji,jj) = puu(ji,jj,ikbu,Kbb) - puu_b(ji,jj,Kbb) 1438 zv_i(ji,jj) = pvv(ji,jj,ikbv,Kbb) - pvv_b(ji,jj,Kbb) 1439 END_2D 1440 ENDIF 1441 ! 1442 IF( ln_wd_il ) THEN ! W/D : use the "clipped" bottom friction !!gm explain WHY, please ! 1443 zztmp = -1._wp / rDt_e 1444 DO_2D( 0, 0, 0, 0 ) 1445 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + zu_i(ji,jj) * wdrampu(ji,jj) * MAX( & 1446 & r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) , zztmp ) 1447 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + zv_i(ji,jj) * wdrampv(ji,jj) * MAX( & 1448 & r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp ) 1449 END_2D 1450 ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation) 1451 1452 DO_2D( 0, 0, 0, 0 ) 1453 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj) 1454 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj) 1455 END_2D 1456 END IF 1457 ! 1458 ! !== TOP stress contribution from baroclinic velocities ==! (no W/D case) 1459 ! 1460 IF( ln_isfcav ) THEN 1461 ! 1462 IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity 1463 1464 DO_2D( 0, 0, 0, 0 ) 1465 iktu = miku(ji,jj) 1466 iktv = mikv(ji,jj) 1467 zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm) 1468 zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm) 1469 END_2D 1470 ELSE ! CENTRED integration: use BEFORE top baroclinic velocity 1471 1472 DO_2D( 0, 0, 0, 0 ) 1473 iktu = miku(ji,jj) 1474 iktv = mikv(ji,jj) 1475 zu_i(ji,jj) = puu(ji,jj,iktu,Kbb) - puu_b(ji,jj,Kbb) 1476 zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb) 1477 END_2D 1478 ENDIF 1479 ! 1480 ! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation) 1481 1482 DO_2D( 0, 0, 0, 0 ) 1483 pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj) 1484 pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj) 1485 END_2D 1486 ! 1487 ENDIF 1488 ! 1489 END SUBROUTINE dyn_drg_init 1490 1491 SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in 1492 & za0, za1, za2, za3 ) ! ==> out 1493 !!---------------------------------------------------------------------- 1494 INTEGER ,INTENT(in ) :: jn ! index of sub time step 1495 LOGICAL ,INTENT(in ) :: ll_init ! 1496 REAL(wp),INTENT( out) :: za0, za1, za2, za3 ! Half-step back interpolation coefficient 1497 ! 1498 REAL(wp) :: zepsilon, zgamma ! - - 1499 !!---------------------------------------------------------------------- 1500 ! ! set Half-step back interpolation coefficient 1501 IF ( jn==1 .AND. ll_init ) THEN !* Forward-backward 1502 za0 = 1._wp 1503 za1 = 0._wp 1504 za2 = 0._wp 1505 za3 = 0._wp 1506 ELSEIF( jn==2 .AND. ll_init ) THEN !* AB2-AM3 Coefficients; bet=0 ; gam=-1/6 ; eps=1/12 1507 za0 = 1.0833333333333_wp ! za0 = 1-gam-eps 1508 za1 =-0.1666666666666_wp ! za1 = gam 1509 za2 = 0.0833333333333_wp ! za2 = eps 1510 za3 = 0._wp 1511 ELSE !* AB3-AM4 Coefficients; bet=0.281105 ; eps=0.013 ; gam=0.0880 1512 IF( rn_bt_alpha == 0._wp ) THEN ! Time diffusion 1513 za0 = 0.614_wp ! za0 = 1/2 + gam + 2*eps 1514 za1 = 0.285_wp ! za1 = 1/2 - 2*gam - 3*eps 1515 za2 = 0.088_wp ! za2 = gam 1516 za3 = 0.013_wp ! za3 = eps 1517 ELSE ! no time diffusion 1518 zepsilon = 0.00976186_wp - 0.13451357_wp * rn_bt_alpha 1519 zgamma = 0.08344500_wp - 0.51358400_wp * rn_bt_alpha 1520 za0 = 0.5_wp + zgamma + 2._wp * rn_bt_alpha + 2._wp * zepsilon 1521 za1 = 1._wp - za0 - zgamma - zepsilon 1522 za2 = zgamma 1523 za3 = zepsilon 1524 ENDIF 1525 ENDIF 1526 END SUBROUTINE ts_bck_interp 1527 1528 1583 1529 !!====================================================================== 1584 1530 END MODULE dynspg_ts
Note: See TracChangeset
for help on using the changeset viewer.
