- Timestamp:
- 2016-07-19T10:38:35+02:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/NERC/dev_r5549_BDY_ZEROGRAD/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r5436 r6808 1 1 MODULE dynspg_ts 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend, split-explicit scheme 2 5 !!====================================================================== 3 6 !! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code … … 11 14 !! 3.5 ! 2013-07 (J. Chanut) Switch to Forward-backward time stepping 12 15 !! 3.6 ! 2013-11 (A. Coward) Update for z-tilde compatibility 16 !! 3.7 ! 2015-11 (J. Chanut) free surface simplification 13 17 !!--------------------------------------------------------------------- 14 #if defined key_dynspg_ts || defined key_esopa 18 15 19 !!---------------------------------------------------------------------- 16 !! 'key_dynspg_ts' split explicit free surface17 !! ----------------------------------------------------------------------18 !! dyn_spg_ts : compute surface pressure gradient trend using a time-19 !! splitting scheme and add to the general trend20 !! dyn_spg_ts : compute surface pressure gradient trend using a time-splitting scheme 21 !! dyn_spg_ts_init: initialisation of the time-splitting scheme 22 !! ts_wgt : set time-splitting weights for temporal averaging (or not) 23 !! ts_rst : read/write time-splitting fields in restart file 20 24 !!---------------------------------------------------------------------- 21 25 USE oce ! ocean dynamics and tracers 22 26 USE dom_oce ! ocean space and time domain 23 27 USE sbc_oce ! surface boundary condition: ocean 28 USE zdf_oce ! Bottom friction coefts 24 29 USE sbcisf ! ice shelf variable (fwfisf) 25 USE dynspg_oce ! surface pressure gradient variables 30 USE sbcapr ! surface boundary condition: atmospheric pressure 31 USE dynadv , ONLY: ln_dynadv_vec 26 32 USE phycst ! physical constants 27 33 USE dynvor ! vorticity term 34 USE wet_dry ! wetting/drying flux limter 28 35 USE bdy_par ! for lk_bdy 29 USE bdytides ! open boundary condition data 36 USE bdytides ! open boundary condition data 30 37 USE bdydyn2d ! open boundary conditions on barotropic variables 31 38 USE sbctide ! tides 32 39 USE updtide ! tide potential 40 ! 41 USE in_out_manager ! I/O manager 33 42 USE lib_mpp ! distributed memory computing library 34 43 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 35 44 USE prtctl ! Print control 36 USE in_out_manager ! I/O manager37 45 USE iom ! IOM library 38 46 USE restart ! only for lrst_oce 39 USE zdf_oce ! Bottom friction coefts40 47 USE wrk_nemo ! Memory Allocation 41 48 USE timing ! Timing 42 USE sbcapr ! surface boundary condition: atmospheric pressure 43 USE dynadv, ONLY: ln_dynadv_vec 49 USE diatmb ! Top,middle,bottom output 44 50 #if defined key_agrif 45 51 USE agrif_opa_interp ! agrif … … 49 55 #endif 50 56 57 51 58 IMPLICIT NONE 52 59 PRIVATE … … 60 67 REAL(wp),SAVE :: rdtbt ! Barotropic time step 61 68 62 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: & 63 wgtbtp1, & ! Primary weights used for time filtering of barotropic variables 64 wgtbtp2 ! Secondary weights used for time filtering of barotropic variables 65 66 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz ! ff/h at F points 67 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter 68 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme) 69 70 ! Arrays below are saved to allow testing of the "no time averaging" option 71 ! If this option is not retained, these could be replaced by temporary arrays 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshbb_e, sshb_e, & ! Instantaneous barotropic arrays 73 ubb_e, ub_e, & 74 vbb_e, vb_e 69 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:) :: wgtbtp1, wgtbtp2 !: 1st & 2nd weights used in time filtering of barotropic fields 70 71 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zwz !: ff/h at F points 72 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne !: triad of coriolis parameter 73 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse !: (only used with een vorticity scheme) 74 75 !! Time filtered arrays at baroclinic time step: 76 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step 75 77 76 78 !! * Substitutions 77 # include "domzgr_substitute.h90"78 79 # include "vectopt_loop_substitute.h90" 79 80 !!---------------------------------------------------------------------- … … 91 92 !!---------------------------------------------------------------------- 92 93 ierr(:) = 0 93 94 ALLOCATE( sshb_e(jpi,jpj), sshbb_e(jpi,jpj), & 95 & ub_e(jpi,jpj) , vb_e(jpi,jpj) , & 96 & ubb_e(jpi,jpj) , vbb_e(jpi,jpj) , STAT= ierr(1) ) 97 98 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT= ierr(2) ) 99 100 IF( ln_dynvor_een .or. ln_dynvor_een_old ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 101 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(3) ) 102 103 dyn_spg_ts_alloc = MAXVAL(ierr(:)) 104 94 ! 95 ALLOCATE( wgtbtp1(3*nn_baro), wgtbtp2(3*nn_baro), zwz(jpi,jpj), STAT=ierr(1) ) 96 ! 97 IF( ln_dynvor_een ) ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , & 98 & ftsw(jpi,jpj) , ftse(jpi,jpj) , STAT=ierr(2) ) 99 ! 100 ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) ) 101 ! 102 dyn_spg_ts_alloc = MAXVAL( ierr(:) ) 103 ! 105 104 IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc ) 106 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn spg_oce_alloc: failed to allocate arrays')105 IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dyn_spg_ts_alloc: failed to allocate arrays') 107 106 ! 108 107 END FUNCTION dyn_spg_ts_alloc 108 109 109 110 110 SUBROUTINE dyn_spg_ts( kt ) 111 111 !!---------------------------------------------------------------------- 112 112 !! 113 !! ** Purpose : 114 !! -Compute the now trend due to the explicit time stepping115 !! of the quasi-linear barotropic system.113 !! ** Purpose : - Compute the now trend due to the explicit time stepping 114 !! of the quasi-linear barotropic system, and add it to the 115 !! general momentum trend. 116 116 !! 117 !! ** Method : 117 !! ** Method : - split-explicit schem (time splitting) : 118 118 !! Barotropic variables are advanced from internal time steps 119 119 !! "n" to "n+1" if ln_bt_fw=T … … 129 129 !! continuity equation taken at the baroclinic time steps. This 130 130 !! ensures tracers conservation. 131 !! - Update 3d trend (ua, va)with barotropic component.131 !! - (ua, va) momentum trend updated with barotropic component. 132 132 !! 133 !! References : Shchepetkin, A.F. and J.C. McWilliams, 2005: 134 !! The regional oceanic modeling system (ROMS): 135 !! a split-explicit, free-surface, 136 !! topography-following-coordinate oceanic model. 137 !! Ocean Modelling, 9, 347-404. 133 !! References : Shchepetkin and McWilliams, Ocean Modelling, 2005. 138 134 !!--------------------------------------------------------------------- 139 !140 135 INTEGER, INTENT(in) :: kt ! ocean time-step index 141 136 ! 142 137 LOGICAL :: ll_fw_start ! if true, forward integration 143 138 LOGICAL :: ll_init ! if true, special startup of 2d equations 139 LOGICAL :: ll_tmp1, ll_tmp2 ! local logical variables used in W/D 144 140 INTEGER :: ji, jj, jk, jn ! dummy loop indices 145 141 INTEGER :: ikbu, ikbv, noffset ! local integers 142 INTEGER :: iktu, iktv ! local integers 143 REAL(wp) :: zmdi 146 144 REAL(wp) :: zraur, z1_2dt_b, z2dt_bf ! local scalars 147 REAL(wp) :: zx1, zy1, zx2, zy2 ! - -148 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 149 REAL(wp) :: zu_spg, zv_spg ! - -150 REAL(wp) :: zhura, zhvra 151 REAL(wp) :: za0, za1, za2, za3 152 ! 153 REAL(wp), POINTER, DIMENSION(:,:) :: z un_e, zvn_e, zsshp2_e145 REAL(wp) :: zx1, zy1, zx2, zy2 ! - - 146 REAL(wp) :: z1_12, z1_8, z1_4, z1_2 ! - - 147 REAL(wp) :: zu_spg, zv_spg ! - - 148 REAL(wp) :: zhura, zhvra ! - - 149 REAL(wp) :: za0, za1, za2, za3 ! - - 150 ! 151 REAL(wp), POINTER, DIMENSION(:,:) :: zsshp2_e 154 152 REAL(wp), POINTER, DIMENSION(:,:) :: zu_trd, zv_trd, zu_frc, zv_frc, zssh_frc 155 REAL(wp), POINTER, DIMENSION(:,:) :: z u_sum, zv_sum, zwx, zwy, zhdiv153 REAL(wp), POINTER, DIMENSION(:,:) :: zwx, zwy, zhdiv 156 154 REAL(wp), POINTER, DIMENSION(:,:) :: zhup2_e, zhvp2_e, zhust_e, zhvst_e 157 155 REAL(wp), POINTER, DIMENSION(:,:) :: zsshu_a, zsshv_a 158 156 REAL(wp), POINTER, DIMENSION(:,:) :: zhf 157 REAL(wp), POINTER, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef. 158 REAL(wp), POINTER, DIMENSION(:,:) :: wduflt1, wdvflt1 ! Wetting/Dying velocity filter coef. 159 159 !!---------------------------------------------------------------------- 160 160 ! … … 162 162 ! 163 163 ! !* Allocate temporary arrays 164 CALL wrk_alloc( jpi, jpj, zsshp2_e, zhdiv ) 165 CALL wrk_alloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 166 CALL wrk_alloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc) 167 CALL wrk_alloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 168 CALL wrk_alloc( jpi, jpj, zsshu_a, zsshv_a ) 169 CALL wrk_alloc( jpi, jpj, zhf ) 170 ! 164 CALL wrk_alloc( jpi,jpj, zsshp2_e, zhdiv ) 165 CALL wrk_alloc( jpi,jpj, zu_trd, zv_trd) 166 CALL wrk_alloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc) 167 CALL wrk_alloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e) 168 CALL wrk_alloc( jpi,jpj, zsshu_a, zsshv_a ) 169 CALL wrk_alloc( jpi,jpj, zhf ) 170 IF( ln_wd ) CALL wrk_alloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 171 ! 172 zmdi=1.e+20 ! missing data indicator for masking 171 173 ! !* Local constant initialization 172 174 z1_12 = 1._wp / 12._wp … … 175 177 z1_2 = 0.5_wp 176 178 zraur = 1._wp / rau0 177 ! 178 IF( kt == nit000 .AND. neuler == 0 ) THEN ! reciprocal of baroclinic time step 179 z2dt_bf = rdt 180 ELSE 181 z2dt_bf = 2.0_wp * rdt 179 ! ! reciprocal of baroclinic time step 180 IF( kt == nit000 .AND. neuler == 0 ) THEN ; z2dt_bf = rdt 181 ELSE ; z2dt_bf = 2.0_wp * rdt 182 182 ENDIF 183 183 z1_2dt_b = 1.0_wp / z2dt_bf 184 184 ! 185 ll_init = ln_bt_av! if no time averaging, then no specific restart185 ll_init = ln_bt_av ! if no time averaging, then no specific restart 186 186 ll_fw_start = .FALSE. 187 ! 188 ! time offset in steps for bdy data update 189 IF (.NOT.ln_bt_fw) THEN ; noffset=-2*nn_baro ; ELSE ; noffset = 0 ; ENDIF 187 ! ! time offset in steps for bdy data update 188 IF( .NOT.ln_bt_fw ) THEN ; noffset = - nn_baro 189 ELSE ; noffset = 0 190 ENDIF 190 191 ! 191 192 IF( kt == nit000 ) THEN !* initialisation … … 196 197 IF(lwp) WRITE(numout,*) 197 198 ! 198 IF (neuler==0)ll_init=.TRUE.199 ! 200 IF (ln_bt_fw.OR.(neuler==0)) THEN201 ll_fw_start=.TRUE.202 noffset= 0199 IF( neuler == 0 ) ll_init=.TRUE. 200 ! 201 IF( ln_bt_fw .OR. neuler == 0 ) THEN 202 ll_fw_start =.TRUE. 203 noffset = 0 203 204 ELSE 204 ll_fw_start=.FALSE.205 ll_fw_start =.FALSE. 205 206 ENDIF 206 207 ! 207 208 ! Set averaging weights and cycle length: 208 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) 209 ! 209 CALL ts_wgt( ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2 ) 210 210 ! 211 211 ENDIF … … 218 218 ! and update depths at T-F points (ht and zhf resp.) at each barotropic time step 219 219 ! 220 IF ( kt == nit000 .OR. lk_vvl ) THEN 221 IF ( ln_dynvor_een_old ) THEN 222 DO jj = 1, jpjm1 223 DO ji = 1, jpim1 224 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 225 & ht(ji ,jj ) + ht(ji+1,jj ) ) / 4._wp 226 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 227 END DO 228 END DO 220 IF( kt == nit000 .OR. .NOT.ln_linssh ) THEN 221 IF( ln_dynvor_een ) THEN !== EEN scheme ==! 222 SELECT CASE( nn_een_e3f ) !* ff/e3 at F-point 223 CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4) 224 DO jj = 1, jpjm1 225 DO ji = 1, jpim1 226 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 227 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) * 0.25_wp 228 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 229 END DO 230 END DO 231 CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask) 232 DO jj = 1, jpjm1 233 DO ji = 1, jpim1 234 zwz(ji,jj) = ( ht_n(ji ,jj+1) + ht_n(ji+1,jj+1) + & 235 & ht_n(ji ,jj ) + ht_n(ji+1,jj ) ) & 236 & / ( MAX( 1._wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 237 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 238 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = ff(ji,jj) / zwz(ji,jj) 239 END DO 240 END DO 241 END SELECT 229 242 CALL lbc_lnk( zwz, 'F', 1._wp ) 230 zwz(:,:) = ff(:,:) * zwz(:,:) 231 243 ! 232 244 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 233 245 DO jj = 2, jpj 234 DO ji = fs_2, jpi ! vector opt.246 DO ji = 2, jpi 235 247 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 236 248 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) … … 239 251 END DO 240 252 END DO 241 ELSE IF ( ln_dynvor_een ) THEN 242 DO jj = 1, jpjm1 243 DO ji = 1, jpim1 244 zwz(ji,jj) = ( ht(ji ,jj+1) + ht(ji+1,jj+1) + & 245 & ht(ji ,jj ) + ht(ji+1,jj ) ) & 246 & / ( MAX( 1.0_wp, tmask(ji ,jj+1, 1) + tmask(ji+1,jj+1, 1) + & 247 & tmask(ji ,jj , 1) + tmask(ji+1,jj , 1) ) ) 248 IF( zwz(ji,jj) /= 0._wp ) zwz(ji,jj) = 1._wp / zwz(ji,jj) 249 END DO 250 END DO 251 CALL lbc_lnk( zwz, 'F', 1._wp ) 252 zwz(:,:) = ff(:,:) * zwz(:,:) 253 254 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 255 DO jj = 2, jpj 256 DO ji = fs_2, jpi ! vector opt. 257 ftne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) 258 ftnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) 259 ftse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) 260 ftsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) 261 END DO 262 END DO 263 ELSE 253 ! 254 ELSE !== all other schemes (ENE, ENS, MIX) 264 255 zwz(:,:) = 0._wp 265 zhf(:,:) = 0. 256 zhf(:,:) = 0._wp 266 257 IF ( .not. ln_sco ) THEN 258 259 !!gm agree the JC comment : this should be done in a much clear way 260 267 261 ! JC: It not clear yet what should be the depth at f-points over land in z-coordinate case 268 262 ! Set it to zero for the time being … … 276 270 277 271 DO jj = 1, jpjm1 278 zhf(:,jj) = zhf(:,jj) *(1._wp- umask(:,jj,1) * umask(:,jj+1,1))272 zhf(:,jj) = zhf(:,jj) * (1._wp- umask(:,jj,1) * umask(:,jj+1,1)) 279 273 END DO 280 274 281 275 DO jk = 1, jpkm1 282 276 DO jj = 1, jpjm1 283 zhf(:,jj) = zhf(:,jj) + fse3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk)277 zhf(:,jj) = zhf(:,jj) + e3f_n(:,jj,jk) * umask(:,jj,jk) * umask(:,jj+1,jk) 284 278 END DO 285 279 END DO … … 297 291 ! If forward start at previous time step, and centered integration, 298 292 ! then update averaging weights: 299 IF ( (.NOT.ln_bt_fw).AND.((neuler==0).AND.(kt==nit000+1))) THEN293 IF (.NOT.ln_bt_fw .AND.( neuler==0 .AND. kt==nit000+1 ) ) THEN 300 294 ll_fw_start=.FALSE. 301 295 CALL ts_wgt(ln_bt_av, ll_fw_start, icycle, wgtbtp1, wgtbtp2) … … 313 307 ! 314 308 DO jk = 1, jpkm1 315 zu_frc(:,:) = zu_frc(:,:) + fse3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk)316 zv_frc(:,:) = zv_frc(:,:) + fse3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk)309 zu_frc(:,:) = zu_frc(:,:) + e3u_n(:,:,jk) * ua(:,:,jk) * umask(:,:,jk) 310 zv_frc(:,:) = zv_frc(:,:) + e3v_n(:,:,jk) * va(:,:,jk) * vmask(:,:,jk) 317 311 END DO 318 312 ! 319 zu_frc(:,:) = zu_frc(:,:) * hur(:,:)320 zv_frc(:,:) = zv_frc(:,:) * hvr(:,:)313 zu_frc(:,:) = zu_frc(:,:) * r1_hu_n(:,:) 314 zv_frc(:,:) = zv_frc(:,:) * r1_hv_n(:,:) 321 315 ! 322 316 ! … … 332 326 ! !* barotropic Coriolis trends (vorticity scheme dependent) 333 327 ! ! -------------------------------------------------------- 334 zwx(:,:) = un_b(:,:) * hu (:,:) * e2u(:,:) ! now fluxes335 zwy(:,:) = vn_b(:,:) * hv (:,:) * e1v(:,:)328 zwx(:,:) = un_b(:,:) * hu_n(:,:) * e2u(:,:) ! now fluxes 329 zwy(:,:) = vn_b(:,:) * hv_n(:,:) * e1v(:,:) 336 330 ! 337 331 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme 338 332 DO jj = 2, jpjm1 339 333 DO ji = fs_2, fs_jpim1 ! vector opt. 340 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)341 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)342 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) /e2v(ji,jj)343 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) /e2v(ji,jj)334 zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 335 zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 336 zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 337 zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 344 338 ! energy conserving formulation for planetary vorticity term 345 339 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) … … 352 346 DO ji = fs_2, fs_jpim1 ! vector opt. 353 347 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 354 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)348 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 355 349 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 356 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)350 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 357 351 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 358 352 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 360 354 END DO 361 355 ! 362 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old) THEN ! enstrophy and energy conserving scheme356 ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme 363 357 DO jj = 2, jpjm1 364 358 DO ji = fs_2, fs_jpim1 ! vector opt. 365 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &366 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &367 & + ftse(ji,jj ) * zwy(ji ,jj-1) &368 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )369 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &370 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &371 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &372 & + ftne(ji,jj ) * zwx(ji ,jj ) )359 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 360 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 361 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 362 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 363 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 364 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 365 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 366 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 373 367 END DO 374 368 END DO … … 378 372 ! !* Right-Hand-Side of the barotropic momentum equation 379 373 ! ! ---------------------------------------------------- 380 IF( lk_vvl ) THEN ! Variable volume : remove surface pressure gradient 381 DO jj = 2, jpjm1 382 DO ji = fs_2, fs_jpim1 ! vector opt. 383 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) / e1u(ji,jj) 384 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) / e2v(ji,jj) 385 END DO 386 END DO 374 IF( .NOT.ln_linssh ) THEN ! Variable volume : remove surface pressure gradient 375 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 376 wduflt1(:,:) = 1.0_wp 377 wdvflt1(:,:) = 1.0_wp 378 DO jj = 2, jpjm1 379 DO ji = 2, jpim1 380 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 381 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji+1,jj) + bathy(ji+1,jj)) & 382 & > rn_wdmin1 + rn_wdmin2 383 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji+1,jj)) > MAX(-bathy(ji,jj), -bathy(ji+1,jj)) & 384 & + rn_wdmin1 + rn_wdmin2 385 IF(ll_tmp1) THEN 386 zcpx(ji,jj) = 1.0_wp 387 ELSEIF(ll_tmp2) THEN 388 ! no worries about sshn(ji+1,jj)-sshn(ji,jj) = 0, it won't happen here 389 zcpx(ji,jj) = ABS((sshn(ji+1,jj) + bathy(ji+1,jj) - sshn(ji,jj) - bathy(ji,jj)) & 390 & /(sshn(ji+1,jj) - sshn(ji,jj))) 391 ELSE 392 zcpx(ji,jj) = 0._wp 393 wduflt1(ji,jj) = 0.0_wp 394 END IF 395 396 ll_tmp1 = MIN(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 397 & .and. MAX(sshn(ji,jj) + bathy(ji,jj), sshn(ji,jj+1) + bathy(ji,jj+1)) & 398 & > rn_wdmin1 + rn_wdmin2 399 ll_tmp2 = MAX(sshn(ji,jj), sshn(ji,jj+1)) > MAX(-bathy(ji,jj), -bathy(ji,jj+1)) & 400 & + rn_wdmin1 + rn_wdmin2 401 IF(ll_tmp1) THEN 402 zcpy(ji,jj) = 1.0_wp 403 ELSEIF(ll_tmp2) THEN 404 ! no worries about sshn(ji,jj+1)-sshn(ji,jj) = 0, it won't happen here 405 zcpy(ji,jj) = ABS((sshn(ji,jj+1) + bathy(ji,jj+1) - sshn(ji,jj) - bathy(ji,jj)) & 406 & /(sshn(ji,jj+1) - sshn(ji,jj))) 407 ELSE 408 zcpy(ji,jj) = 0._wp 409 wdvflt1(ji,jj) = 0.0_wp 410 ENDIF 411 412 END DO 413 END DO 414 415 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 416 417 DO jj = 2, jpjm1 418 DO ji = 2, jpim1 419 zu_trd(ji,jj) = ( zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) & 420 & * r1_e1u(ji,jj) ) * zcpx(ji,jj) * wduflt1(ji,jj) 421 zv_trd(ji,jj) = ( zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) & 422 & * r1_e2v(ji,jj) ) * zcpy(ji,jj) * wdvflt1(ji,jj) 423 END DO 424 END DO 425 426 ELSE 427 428 DO jj = 2, jpjm1 429 DO ji = fs_2, fs_jpim1 ! vector opt. 430 zu_trd(ji,jj) = zu_trd(ji,jj) - grav * ( sshn(ji+1,jj ) - sshn(ji ,jj ) ) * r1_e1u(ji,jj) 431 zv_trd(ji,jj) = zv_trd(ji,jj) - grav * ( sshn(ji ,jj+1) - sshn(ji ,jj ) ) * r1_e2v(ji,jj) 432 END DO 433 END DO 434 ENDIF 435 387 436 ENDIF 388 437 389 438 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 390 439 DO ji = fs_2, fs_jpim1 391 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * umask(ji,jj,1)392 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * vmask(ji,jj,1)440 zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj) 441 zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj) 393 442 END DO 394 443 END DO … … 416 465 ! 417 466 ! Note that the "unclipped" bottom friction parameter is used even with explicit drag 418 zu_frc(:,:) = zu_frc(:,:) + hur(:,:) * bfrua(:,:) * zwx(:,:) 419 zv_frc(:,:) = zv_frc(:,:) + hvr(:,:) * bfrva(:,:) * zwy(:,:) 467 IF( ln_wd ) THEN 468 zu_frc(:,:) = zu_frc(:,:) + MAX(r1_hu_n(:,:) * bfrua(:,:),-1._wp / rdtbt) * zwx(:,:) 469 zv_frc(:,:) = zv_frc(:,:) + MAX(r1_hv_n(:,:) * bfrva(:,:),-1._wp / rdtbt) * zwy(:,:) 470 ELSE 471 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * bfrua(:,:) * zwx(:,:) 472 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * bfrva(:,:) * zwy(:,:) 473 END IF 474 ! 475 ! ! Add top stress contribution from baroclinic velocities: 476 IF (ln_bt_fw) THEN 477 DO jj = 2, jpjm1 478 DO ji = fs_2, fs_jpim1 ! vector opt. 479 iktu = miku(ji,jj) 480 iktv = mikv(ji,jj) 481 zwx(ji,jj) = un(ji,jj,iktu) - un_b(ji,jj) ! NOW top baroclinic velocities 482 zwy(ji,jj) = vn(ji,jj,iktv) - vn_b(ji,jj) 483 END DO 484 END DO 485 ELSE 486 DO jj = 2, jpjm1 487 DO ji = fs_2, fs_jpim1 ! vector opt. 488 iktu = miku(ji,jj) 489 iktv = mikv(ji,jj) 490 zwx(ji,jj) = ub(ji,jj,iktu) - ub_b(ji,jj) ! BEFORE top baroclinic velocities 491 zwy(ji,jj) = vb(ji,jj,iktv) - vb_b(ji,jj) 492 END DO 493 END DO 494 ENDIF 495 ! 496 ! Note that the "unclipped" top friction parameter is used even with explicit drag 497 zu_frc(:,:) = zu_frc(:,:) + r1_hu_n(:,:) * tfrua(:,:) * zwx(:,:) 498 zv_frc(:,:) = zv_frc(:,:) + r1_hv_n(:,:) * tfrva(:,:) * zwy(:,:) 420 499 ! 421 500 IF (ln_bt_fw) THEN ! Add wind forcing 422 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * hur(:,:)423 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * hvr(:,:)501 zu_frc(:,:) = zu_frc(:,:) + zraur * utau(:,:) * r1_hu_n(:,:) 502 zv_frc(:,:) = zv_frc(:,:) + zraur * vtau(:,:) * r1_hv_n(:,:) 424 503 ELSE 425 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * hur(:,:)426 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * hvr(:,:)504 zu_frc(:,:) = zu_frc(:,:) + zraur * z1_2 * ( utau_b(:,:) + utau(:,:) ) * r1_hu_n(:,:) 505 zv_frc(:,:) = zv_frc(:,:) + zraur * z1_2 * ( vtau_b(:,:) + vtau(:,:) ) * r1_hv_n(:,:) 427 506 ENDIF 428 507 ! … … 431 510 DO jj = 2, jpjm1 432 511 DO ji = fs_2, fs_jpim1 ! vector opt. 433 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) /e1u(ji,jj)434 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) /e2v(ji,jj)512 zu_spg = grav * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) ) * r1_e1u(ji,jj) 513 zv_spg = grav * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) ) * r1_e2v(ji,jj) 435 514 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 436 515 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 441 520 DO ji = fs_2, fs_jpim1 ! vector opt. 442 521 zu_spg = grav * z1_2 * ( ssh_ib (ji+1,jj ) - ssh_ib (ji,jj) & 443 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) /e1u(ji,jj)522 & + ssh_ibb(ji+1,jj ) - ssh_ibb(ji,jj) ) * r1_e1u(ji,jj) 444 523 zv_spg = grav * z1_2 * ( ssh_ib (ji ,jj+1) - ssh_ib (ji,jj) & 445 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) /e2v(ji,jj)524 & + ssh_ibb(ji ,jj+1) - ssh_ibb(ji,jj) ) * r1_e2v(ji,jj) 446 525 zu_frc(ji,jj) = zu_frc(ji,jj) + zu_spg 447 526 zv_frc(ji,jj) = zv_frc(ji,jj) + zv_spg … … 454 533 ! ! Surface net water flux and rivers 455 534 IF (ln_bt_fw) THEN 456 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + rdivisf *fwfisf(:,:) )535 zssh_frc(:,:) = zraur * ( emp(:,:) - rnf(:,:) + fwfisf(:,:) ) 457 536 ELSE 458 537 zssh_frc(:,:) = zraur * z1_2 * ( emp(:,:) + emp_b(:,:) - rnf(:,:) - rnf_b(:,:) & 459 & + rdivisf * ( fwfisf(:,:) + fwfisf_b(:,:) ))538 & + fwfisf(:,:) + fwfisf_b(:,:) ) 460 539 ENDIF 461 540 #if defined key_asminc … … 465 544 ENDIF 466 545 #endif 467 ! !* Fill boundary data arrays withAGRIF468 ! ! ------------------------------------ -546 ! !* Fill boundary data arrays for AGRIF 547 ! ! ------------------------------------ 469 548 #if defined key_agrif 470 549 IF( .NOT.Agrif_Root() ) CALL agrif_dta_ts( kt ) … … 487 566 vb_e (:,:) = 0._wp 488 567 ENDIF 568 569 IF( ln_wd ) THEN !preserve the positivity of water depth 570 !ssh[b,n,a] should have already been processed for this 571 sshbb_e(:,:) = MAX(sshbb_e(:,:), rn_wdmin1 - bathy(:,:)) 572 sshb_e(:,:) = MAX(sshb_e(:,:) , rn_wdmin1 - bathy(:,:)) 573 ENDIF 489 574 ! 490 575 IF (ln_bt_fw) THEN ! FORWARD integration: start from NOW fields 491 sshn_e(:,:) = sshn(:,:)492 zun_e (:,:) = un_b(:,:)493 zvn_e (:,:) = vn_b(:,:)494 ! 495 hu_e (:,:) = hu(:,:)496 hv_e (:,:) = hv(:,:)497 hur_e (:,:) = hur(:,:)498 hvr_e (:,:) = hvr(:,:)576 sshn_e(:,:) = sshn(:,:) 577 un_e (:,:) = un_b(:,:) 578 vn_e (:,:) = vn_b(:,:) 579 ! 580 hu_e (:,:) = hu_n(:,:) 581 hv_e (:,:) = hv_n(:,:) 582 hur_e (:,:) = r1_hu_n(:,:) 583 hvr_e (:,:) = r1_hv_n(:,:) 499 584 ELSE ! CENTRED integration: start from BEFORE fields 500 sshn_e(:,:) = sshb(:,:)501 zun_e (:,:) = ub_b(:,:)502 zvn_e (:,:) = vb_b(:,:)503 ! 504 hu_e (:,:) = hu_b(:,:)505 hv_e (:,:) = hv_b(:,:)506 hur_e (:,:) = hur_b(:,:)507 hvr_e (:,:) = hvr_b(:,:)585 sshn_e(:,:) = sshb(:,:) 586 un_e (:,:) = ub_b(:,:) 587 vn_e (:,:) = vb_b(:,:) 588 ! 589 hu_e (:,:) = hu_b(:,:) 590 hv_e (:,:) = hv_b(:,:) 591 hur_e (:,:) = r1_hu_b(:,:) 592 hvr_e (:,:) = r1_hv_b(:,:) 508 593 ENDIF 509 594 ! … … 514 599 va_b (:,:) = 0._wp 515 600 ssha (:,:) = 0._wp ! Sum for after averaged sea level 516 zu_sum(:,:) = 0._wp ! Sum for now transport issued from ts loop517 zv_sum(:,:) = 0._wp601 un_adv(:,:) = 0._wp ! Sum for now transport issued from ts loop 602 vn_adv(:,:) = 0._wp 518 603 ! ! ==================== ! 519 604 DO jn = 1, icycle ! sub-time-step loop ! … … 523 608 ! Update only tidal forcing at open boundaries 524 609 #if defined key_tide 525 IF ( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset=(noffset+1))526 IF ( ln_tide_pot .AND. lk_tide ) CALL upd_tide( kt, kit=jn, koffset=noffset)610 IF( lk_bdy .AND. lk_tide ) CALL bdy_dta_tides( kt, kit=jn, time_offset= noffset+1 ) 611 IF( ln_tide_pot .AND. lk_tide ) CALL upd_tide ( kt, kit=jn, time_offset= noffset ) 527 612 #endif 528 613 ! … … 539 624 540 625 ! Extrapolate barotropic velocities at step jit+0.5: 541 ua_e(:,:) = za1 * zun_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:)542 va_e(:,:) = za1 * zvn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:)543 544 IF( lk_vvl ) THEN!* Update ocean depth (variable volume case only)626 ua_e(:,:) = za1 * un_e(:,:) + za2 * ub_e(:,:) + za3 * ubb_e(:,:) 627 va_e(:,:) = za1 * vn_e(:,:) + za2 * vb_e(:,:) + za3 * vbb_e(:,:) 628 629 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 545 630 ! ! ------------------ 546 631 ! Extrapolate Sea Level at step jit+0.5: … … 549 634 DO jj = 2, jpjm1 ! Sea Surface Height at u- & v-points 550 635 DO ji = 2, fs_jpim1 ! Vector opt. 551 zwx(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &552 & * ( e1 2t(ji ,jj) * zsshp2_e(ji ,jj) &553 & + e1 2t(ji+1,jj) * zsshp2_e(ji+1,jj) )554 zwy(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &555 & * ( e1 2t(ji,jj ) * zsshp2_e(ji,jj ) &556 & + e1 2t(ji,jj+1) * zsshp2_e(ji,jj+1) )636 zwx(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 637 & * ( e1e2t(ji ,jj) * zsshp2_e(ji ,jj) & 638 & + e1e2t(ji+1,jj) * zsshp2_e(ji+1,jj) ) 639 zwy(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 640 & * ( e1e2t(ji,jj ) * zsshp2_e(ji,jj ) & 641 & + e1e2t(ji,jj+1) * zsshp2_e(ji,jj+1) ) 557 642 END DO 558 643 END DO … … 561 646 zhup2_e (:,:) = hu_0(:,:) + zwx(:,:) ! Ocean depth at U- and V-points 562 647 zhvp2_e (:,:) = hv_0(:,:) + zwy(:,:) 648 IF( ln_wd ) THEN 649 zhup2_e(:,:) = MAX(zhup2_e (:,:), rn_wdmin1) 650 zhvp2_e(:,:) = MAX(zhvp2_e (:,:), rn_wdmin1) 651 END IF 563 652 ELSE 564 zhup2_e (:,:) = hu (:,:)565 zhvp2_e (:,:) = hv (:,:)653 zhup2_e (:,:) = hu_n(:,:) 654 zhvp2_e (:,:) = hv_n(:,:) 566 655 ENDIF 567 656 ! !* after ssh … … 574 663 ! 575 664 #if defined key_agrif 576 ! Set fluxes during predictor step to ensure 577 ! volume conservation 578 IF( (.NOT.Agrif_Root()).AND.ln_bt_fw ) THEN 665 ! Set fluxes during predictor step to ensure volume conservation 666 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 579 667 IF((nbondi == -1).OR.(nbondi == 2)) THEN 580 668 DO jj=1,jpj … … 599 687 ENDIF 600 688 #endif 689 IF( ln_wd ) CALL wad_lmt_bt(zwx, zwy, sshn_e, zssh_frc, rdtbt) 601 690 ! 602 691 ! Sum over sub-time-steps to compute advective velocities 603 692 za2 = wgtbtp2(jn) 604 zu_sum (:,:) = zu_sum (:,:) + za2 * zwx (:,:) / e2u(:,:)605 zv_sum (:,:) = zv_sum (:,:) + za2 * zwy (:,:) / e1v(:,:)693 un_adv(:,:) = un_adv(:,:) + za2 * zwx(:,:) * r1_e2u(:,:) 694 vn_adv(:,:) = vn_adv(:,:) + za2 * zwy(:,:) * r1_e1v(:,:) 606 695 ! 607 696 ! Set next sea level: … … 609 698 DO ji = fs_2, fs_jpim1 ! vector opt. 610 699 zhdiv(ji,jj) = ( zwx(ji,jj) - zwx(ji-1,jj) & 611 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1 2t(ji,jj)700 & + zwy(ji,jj) - zwy(ji,jj-1) ) * r1_e1e2t(ji,jj) 612 701 END DO 613 702 END DO 614 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1) 703 ssha_e(:,:) = ( sshn_e(:,:) - rdtbt * ( zssh_frc(:,:) + zhdiv(:,:) ) ) * ssmask(:,:) 704 IF( ln_wd ) ssha_e(:,:) = MAX(ssha_e(:,:), rn_wdmin1 - bathy(:,:)) 615 705 CALL lbc_lnk( ssha_e, 'T', 1._wp ) 616 706 617 707 #if defined key_bdy 618 ! Duplicate sea level across open boundaries (this is only cosmetic if l k_vvl=.false.)619 IF (lk_bdy)CALL bdy_ssh( ssha_e )708 ! Duplicate sea level across open boundaries (this is only cosmetic if linssh=T) 709 IF( lk_bdy ) CALL bdy_ssh( ssha_e ) 620 710 #endif 621 711 #if defined key_agrif 622 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn )712 IF( .NOT.Agrif_Root() ) CALL agrif_ssh_ts( jn ) 623 713 #endif 624 714 ! 625 715 ! Sea Surface Height at u-,v-points (vvl case only) 626 IF ( lk_vvl) THEN716 IF( .NOT.ln_linssh ) THEN 627 717 DO jj = 2, jpjm1 628 718 DO ji = 2, jpim1 ! NO Vector Opt. 629 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj)&630 & * ( e1 2t(ji ,jj )* ssha_e(ji ,jj ) &631 & + e1 2t(ji+1,jj )* ssha_e(ji+1,jj ) )632 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj)&633 & * ( e1 2t(ji ,jj )* ssha_e(ji ,jj ) &634 & + e1 2t(ji ,jj+1)* ssha_e(ji ,jj+1) )719 zsshu_a(ji,jj) = z1_2 * ssumask(ji,jj) * r1_e1e2u(ji,jj) & 720 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 721 & + e1e2t(ji+1,jj ) * ssha_e(ji+1,jj ) ) 722 zsshv_a(ji,jj) = z1_2 * ssvmask(ji,jj) * r1_e1e2v(ji,jj) & 723 & * ( e1e2t(ji ,jj ) * ssha_e(ji ,jj ) & 724 & + e1e2t(ji ,jj+1) * ssha_e(ji ,jj+1) ) 635 725 END DO 636 726 END DO … … 656 746 za3=0.013_wp ! za3 = eps 657 747 ENDIF 658 748 ! 659 749 zsshp2_e(:,:) = za0 * ssha_e(:,:) + za1 * sshn_e (:,:) & 660 750 & + za2 * sshb_e(:,:) + za3 * sshbb_e(:,:) 661 751 IF( ln_wd ) THEN ! Calculating and applying W/D gravity filters 752 wduflt1(:,:) = 1._wp 753 wdvflt1(:,:) = 1._wp 754 DO jj = 2, jpjm1 755 DO ji = 2, jpim1 756 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 757 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji+1,jj) + bathy(ji+1,jj) ) & 758 & > rn_wdmin1 + rn_wdmin2 759 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji+1,jj) ) > MAX( -bathy(ji,jj), -bathy(ji+1,jj) ) & 760 & + rn_wdmin1 + rn_wdmin2 761 IF(ll_tmp1) THEN 762 zcpx(ji,jj) = 1._wp 763 ELSE IF(ll_tmp2) THEN 764 ! no worries about zsshp2_e(ji+1,jj)-zsshp2_e(ji,jj) = 0, it won't happen here 765 zcpx(ji,jj) = ABS( (zsshp2_e(ji+1,jj) + bathy(ji+1,jj) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 766 & / (zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj)) ) 767 ELSE 768 zcpx(ji,jj) = 0._wp 769 wduflt1(ji,jj) = 0._wp 770 END IF 771 772 ll_tmp1 = MIN( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 773 & .AND. MAX( zsshp2_e(ji,jj) + bathy(ji,jj), zsshp2_e(ji,jj+1) + bathy(ji,jj+1) ) & 774 & > rn_wdmin1 + rn_wdmin2 775 ll_tmp2 = MAX( zsshp2_e(ji,jj), zsshp2_e(ji,jj+1) ) > MAX( -bathy(ji,jj), -bathy(ji,jj+1) ) & 776 & + rn_wdmin1 + rn_wdmin2 777 IF(ll_tmp1) THEN 778 zcpy(ji,jj) = 1._wp 779 ELSE IF(ll_tmp2) THEN 780 ! no worries about zsshp2_e(ji,jj+1)-zsshp2_e(ji,jj) = 0, it won't happen here 781 zcpy(ji,jj) = ABS( (zsshp2_e(ji,jj+1) + bathy(ji,jj+1) - zsshp2_e(ji,jj) - bathy(ji,jj)) & 782 & / (zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj)) ) 783 ELSE 784 zcpy(ji,jj) = 0._wp 785 wdvflt1(ji,jj) = 0._wp 786 END IF 787 END DO 788 END DO 789 CALL lbc_lnk( zcpx, 'U', 1._wp ) ; CALL lbc_lnk( zcpy, 'V', 1._wp ) 790 ENDIF 662 791 ! 663 792 ! Compute associated depths at U and V points: 664 IF ( lk_vvl.AND.(.NOT.ln_dynadv_vec) ) THEN793 IF( .NOT.ln_linssh .AND. .NOT.ln_dynadv_vec ) THEN !* Vector form 665 794 ! 666 795 DO jj = 2, jpjm1 667 796 DO ji = 2, jpim1 668 zx1 = z1_2 * umask(ji ,jj,1) * r1_e12u(ji ,jj) &669 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj) &670 & + e1 2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) )671 zy1 = z1_2 * vmask(ji ,jj,1) * r1_e12v(ji ,jj ) &672 & * ( e1 2t(ji ,jj ) * zsshp2_e(ji ,jj ) &673 & + e1 2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) )797 zx1 = z1_2 * ssumask(ji ,jj) * r1_e1e2u(ji ,jj) & 798 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj) & 799 & + e1e2t(ji+1,jj ) * zsshp2_e(ji+1,jj ) ) 800 zy1 = z1_2 * ssvmask(ji ,jj) * r1_e1e2v(ji ,jj ) & 801 & * ( e1e2t(ji ,jj ) * zsshp2_e(ji ,jj ) & 802 & + e1e2t(ji ,jj+1) * zsshp2_e(ji ,jj+1) ) 674 803 zhust_e(ji,jj) = hu_0(ji,jj) + zx1 675 804 zhvst_e(ji,jj) = hv_0(ji,jj) + zy1 676 805 END DO 677 806 END DO 807 808 IF( ln_wd ) THEN 809 zhust_e(:,:) = MAX(zhust_e (:,:), rn_wdmin1 ) 810 zhvst_e(:,:) = MAX(zhvst_e (:,:), rn_wdmin1 ) 811 END IF 812 678 813 ENDIF 679 814 ! 680 815 ! Add Coriolis trend: 681 ! zwz array below or triads normally depend on sea level with key_vvland should be updated816 ! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated 682 817 ! at each time step. We however keep them constant here for optimization. 683 818 ! Recall that zwx and zwy arrays hold fluxes at this stage: … … 685 820 ! zwy(:,:) = e1v(:,:) * va_e(:,:) * zhvp2_e(:,:) 686 821 ! 687 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN 822 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==! 688 823 DO jj = 2, jpjm1 689 824 DO ji = fs_2, fs_jpim1 ! vector opt. 690 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) /e1u(ji,jj)691 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)692 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) /e2v(ji,jj)693 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)825 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) * r1_e1u(ji,jj) 826 zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 827 zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) * r1_e2v(ji,jj) 828 zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 694 829 zu_trd(ji,jj) = z1_4 * ( zwz(ji ,jj-1) * zy1 + zwz(ji,jj) * zy2 ) 695 830 zv_trd(ji,jj) =-z1_4 * ( zwz(ji-1,jj ) * zx1 + zwz(ji,jj) * zx2 ) … … 697 832 END DO 698 833 ! 699 ELSEIF ( ln_dynvor_ens ) THEN 834 ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==! 700 835 DO jj = 2, jpjm1 701 836 DO ji = fs_2, fs_jpim1 ! vector opt. 702 837 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 703 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) /e1u(ji,jj)838 & + zwy(ji ,jj ) + zwy(ji+1,jj ) ) * r1_e1u(ji,jj) 704 839 zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 705 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) /e2v(ji,jj)840 & + zwx(ji ,jj ) + zwx(ji ,jj+1) ) * r1_e2v(ji,jj) 706 841 zu_trd(ji,jj) = zy1 * ( zwz(ji ,jj-1) + zwz(ji,jj) ) 707 842 zv_trd(ji,jj) = zx1 * ( zwz(ji-1,jj ) + zwz(ji,jj) ) … … 709 844 END DO 710 845 ! 711 ELSEIF ( ln_dynvor_een .or. ln_dynvor_een_old ) THEN!== energy and enstrophy conserving scheme ==!846 ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==! 712 847 DO jj = 2, jpjm1 713 848 DO ji = fs_2, fs_jpim1 ! vector opt. 714 zu_trd(ji,jj) = + z1_12 /e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) &715 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) &716 & + ftse(ji,jj ) * zwy(ji ,jj-1) &717 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )718 zv_trd(ji,jj) = - z1_12 /e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) &719 & + ftse(ji,jj+1) * zwx(ji ,jj+1) &720 & + ftnw(ji,jj ) * zwx(ji-1,jj ) &721 & + ftne(ji,jj ) * zwx(ji ,jj ) )849 zu_trd(ji,jj) = + z1_12 * r1_e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) & 850 & + ftnw(ji+1,jj) * zwy(ji+1,jj ) & 851 & + ftse(ji,jj ) * zwy(ji ,jj-1) & 852 & + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) 853 zv_trd(ji,jj) = - z1_12 * r1_e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) & 854 & + ftse(ji,jj+1) * zwx(ji ,jj+1) & 855 & + ftnw(ji,jj ) * zwx(ji-1,jj ) & 856 & + ftne(ji,jj ) * zwx(ji ,jj ) ) 722 857 END DO 723 858 END DO … … 729 864 DO jj = 2, jpjm1 730 865 DO ji = fs_2, fs_jpim1 ! vector opt. 731 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) /e1u(ji,jj)732 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) /e2v(ji,jj)866 zu_spg = grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) * r1_e1u(ji,jj) 867 zv_spg = grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) * r1_e2v(ji,jj) 733 868 zu_trd(ji,jj) = zu_trd(ji,jj) + zu_spg 734 869 zv_trd(ji,jj) = zv_trd(ji,jj) + zv_spg … … 738 873 ! 739 874 ! Add bottom stresses: 740 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * zun_e(:,:) * hur_e(:,:) 741 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * zvn_e(:,:) * hvr_e(:,:) 875 zu_trd(:,:) = zu_trd(:,:) + bfrua(:,:) * un_e(:,:) * hur_e(:,:) 876 zv_trd(:,:) = zv_trd(:,:) + bfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 877 ! 878 ! Add top stresses: 879 zu_trd(:,:) = zu_trd(:,:) + tfrua(:,:) * un_e(:,:) * hur_e(:,:) 880 zv_trd(:,:) = zv_trd(:,:) + tfrva(:,:) * vn_e(:,:) * hvr_e(:,:) 742 881 ! 743 882 ! Surface pressure trend: 744 DO jj = 2, jpjm1 745 DO ji = fs_2, fs_jpim1 ! vector opt. 746 ! Add surface pressure gradient 747 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) / e1u(ji,jj) 748 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) / e2v(ji,jj) 749 zwx(ji,jj) = zu_spg 750 zwy(ji,jj) = zv_spg 751 END DO 752 END DO 883 884 IF( ln_wd ) THEN 885 DO jj = 2, jpjm1 886 DO ji = 2, jpim1 887 ! Add surface pressure gradient 888 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 889 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 890 zwx(ji,jj) = zu_spg * zcpx(ji,jj) 891 zwy(ji,jj) = zv_spg * zcpy(ji,jj) 892 END DO 893 END DO 894 ELSE 895 DO jj = 2, jpjm1 896 DO ji = fs_2, fs_jpim1 ! vector opt. 897 ! Add surface pressure gradient 898 zu_spg = - grav * ( zsshp2_e(ji+1,jj) - zsshp2_e(ji,jj) ) * r1_e1u(ji,jj) 899 zv_spg = - grav * ( zsshp2_e(ji,jj+1) - zsshp2_e(ji,jj) ) * r1_e2v(ji,jj) 900 zwx(ji,jj) = zu_spg 901 zwy(ji,jj) = zv_spg 902 END DO 903 END DO 904 END IF 905 753 906 ! 754 907 ! Set next velocities: 755 IF( ln_dynadv_vec .OR. (.NOT. lk_vvl) ) THEN !Vector form908 IF( ln_dynadv_vec .OR. ln_linssh ) THEN !* Vector form 756 909 DO jj = 2, jpjm1 757 910 DO ji = fs_2, fs_jpim1 ! vector opt. 758 ua_e(ji,jj) = ( zun_e(ji,jj) &911 ua_e(ji,jj) = ( un_e(ji,jj) & 759 912 & + rdtbt * ( zwx(ji,jj) & 760 913 & + zu_trd(ji,jj) & 761 914 & + zu_frc(ji,jj) ) & 762 & ) * umask(ji,jj,1)763 764 va_e(ji,jj) = ( zvn_e(ji,jj) &915 & ) * ssumask(ji,jj) 916 917 va_e(ji,jj) = ( vn_e(ji,jj) & 765 918 & + rdtbt * ( zwy(ji,jj) & 766 919 & + zv_trd(ji,jj) & 767 920 & + zv_frc(ji,jj) ) & 768 & ) * vmask(ji,jj,1)769 END DO 770 END DO 771 772 ELSE !Flux form921 & ) * ssvmask(ji,jj) 922 END DO 923 END DO 924 ! 925 ELSE !* Flux form 773 926 DO jj = 2, jpjm1 774 927 DO ji = fs_2, fs_jpim1 ! vector opt. 775 928 776 zhura = umask(ji,jj,1)/(hu_0(ji,jj) + zsshu_a(ji,jj) + 1._wp - umask(ji,jj,1)) 777 zhvra = vmask(ji,jj,1)/(hv_0(ji,jj) + zsshv_a(ji,jj) + 1._wp - vmask(ji,jj,1)) 778 779 ua_e(ji,jj) = ( hu_e(ji,jj) * zun_e(ji,jj) & 929 IF( ln_wd ) THEN 930 zhura = MAX(hu_0(ji,jj) + zsshu_a(ji,jj), rn_wdmin1) 931 zhvra = MAX(hv_0(ji,jj) + zsshv_a(ji,jj), rn_wdmin1) 932 ELSE 933 zhura = hu_0(ji,jj) + zsshu_a(ji,jj) 934 zhvra = hv_0(ji,jj) + zsshv_a(ji,jj) 935 END IF 936 zhura = ssumask(ji,jj)/(zhura + 1._wp - ssumask(ji,jj)) 937 zhvra = ssvmask(ji,jj)/(zhvra + 1._wp - ssvmask(ji,jj)) 938 939 ua_e(ji,jj) = ( hu_e(ji,jj) * un_e(ji,jj) & 780 940 & + rdtbt * ( zhust_e(ji,jj) * zwx(ji,jj) & 781 941 & + zhup2_e(ji,jj) * zu_trd(ji,jj) & 782 & + hu(ji,jj) * zu_frc(ji,jj) ) &942 & + hu_n(ji,jj) * zu_frc(ji,jj) ) & 783 943 & ) * zhura 784 944 785 va_e(ji,jj) = ( hv_e(ji,jj) * zvn_e(ji,jj) &945 va_e(ji,jj) = ( hv_e(ji,jj) * vn_e(ji,jj) & 786 946 & + rdtbt * ( zhvst_e(ji,jj) * zwy(ji,jj) & 787 947 & + zhvp2_e(ji,jj) * zv_trd(ji,jj) & 788 & + hv(ji,jj) * zv_frc(ji,jj) ) &948 & + hv_n(ji,jj) * zv_frc(ji,jj) ) & 789 949 & ) * zhvra 790 950 END DO … … 792 952 ENDIF 793 953 ! 794 IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only) 795 ! ! ---------------------------------------------- 796 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 797 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 798 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 799 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 954 IF( .NOT.ln_linssh ) THEN !* Update ocean depth (variable volume case only) 955 IF( ln_wd ) THEN 956 hu_e (:,:) = MAX(hu_0(:,:) + zsshu_a(:,:), rn_wdmin1) 957 hv_e (:,:) = MAX(hv_0(:,:) + zsshv_a(:,:), rn_wdmin1) 958 ELSE 959 hu_e (:,:) = hu_0(:,:) + zsshu_a(:,:) 960 hv_e (:,:) = hv_0(:,:) + zsshv_a(:,:) 961 END IF 962 hur_e(:,:) = ssumask(:,:) / ( hu_e(:,:) + 1._wp - ssumask(:,:) ) 963 hvr_e(:,:) = ssvmask(:,:) / ( hv_e(:,:) + 1._wp - ssvmask(:,:) ) 800 964 ! 801 965 ENDIF 802 ! !* domain lateral boundary 803 ! ! ----------------------- 804 ! 966 ! !* domain lateral boundary 805 967 CALL lbc_lnk_multi( ua_e, 'U', -1._wp, va_e , 'V', -1._wp ) 806 968 ! 807 969 #if defined key_bdy 808 809 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, zun_e, zvn_e, hur_e, hvr_e, ssha_e )970 ! ! open boundaries 971 IF( lk_bdy ) CALL bdy_dyn2d( jn, ua_e, va_e, un_e, vn_e, hur_e, hvr_e, ssha_e ) 810 972 #endif 811 973 #if defined key_agrif … … 815 977 ! ! ---- 816 978 ubb_e (:,:) = ub_e (:,:) 817 ub_e (:,:) = zun_e(:,:)818 zun_e(:,:) = ua_e (:,:)979 ub_e (:,:) = un_e (:,:) 980 un_e (:,:) = ua_e (:,:) 819 981 ! 820 982 vbb_e (:,:) = vb_e (:,:) 821 vb_e (:,:) = zvn_e(:,:)822 zvn_e(:,:) = va_e (:,:)983 vb_e (:,:) = vn_e (:,:) 984 vn_e (:,:) = va_e (:,:) 823 985 ! 824 986 sshbb_e(:,:) = sshb_e(:,:) … … 829 991 ! ! ---------------------- 830 992 za1 = wgtbtp1(jn) 831 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN ! Sum velocities993 IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities 832 994 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) 833 995 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) 834 ELSE 996 ELSE ! Sum transports 835 997 ua_b (:,:) = ua_b (:,:) + za1 * ua_e (:,:) * hu_e (:,:) 836 998 va_b (:,:) = va_b (:,:) + za1 * va_e (:,:) * hv_e (:,:) … … 845 1007 ! ----------------------------------------------------------------------------- 846 1008 ! 847 ! At this stage ssha holds a time averaged value848 ! ! Sea Surface Height at u-,v- and f-points849 IF( lk_vvl ) THEN ! (required only in key_vvl case)850 DO jj = 1, jpjm1851 DO ji = 1, jpim1 ! NO Vector Opt.852 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e12u(ji,jj) &853 & * ( e12t(ji ,jj) * ssha(ji ,jj) &854 & + e12t(ji+1,jj) * ssha(ji+1,jj) )855 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e12v(ji,jj) &856 & * ( e12t(ji,jj ) * ssha(ji,jj ) &857 & + e12t(ji,jj+1) * ssha(ji,jj+1) )858 END DO859 END DO860 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions861 ENDIF862 !863 1009 ! Set advection velocity correction: 864 IF (((kt==nit000).AND.(neuler==0)).OR.(.NOT.ln_bt_fw)) THEN 865 un_adv(:,:) = zu_sum(:,:)*hur(:,:) 866 vn_adv(:,:) = zv_sum(:,:)*hvr(:,:) 1010 zwx(:,:) = un_adv(:,:) 1011 zwy(:,:) = vn_adv(:,:) 1012 IF( ( kt == nit000 .AND. neuler==0 ) .OR. .NOT.ln_bt_fw ) THEN 1013 un_adv(:,:) = zwx(:,:) * r1_hu_n(:,:) 1014 vn_adv(:,:) = zwy(:,:) * r1_hv_n(:,:) 867 1015 ELSE 868 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + z u_sum(:,:)) * hur(:,:)869 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + z v_sum(:,:)) * hvr(:,:)1016 un_adv(:,:) = z1_2 * ( ub2_b(:,:) + zwx(:,:) ) * r1_hu_n(:,:) 1017 vn_adv(:,:) = z1_2 * ( vb2_b(:,:) + zwy(:,:) ) * r1_hv_n(:,:) 870 1018 END IF 871 1019 872 IF (ln_bt_fw) THEN ! Save integrated transport for next computation873 ub2_b(:,:) = z u_sum(:,:)874 vb2_b(:,:) = z v_sum(:,:)1020 IF( ln_bt_fw ) THEN ! Save integrated transport for next computation 1021 ub2_b(:,:) = zwx(:,:) 1022 vb2_b(:,:) = zwy(:,:) 875 1023 ENDIF 876 1024 ! 877 1025 ! Update barotropic trend: 878 IF (( ln_dynadv_vec ).OR. (.NOT. lk_vvl)) THEN1026 IF( ln_dynadv_vec .OR. ln_linssh ) THEN 879 1027 DO jk=1,jpkm1 880 1028 ua(:,:,jk) = ua(:,:,jk) + ( ua_b(:,:) - ub_b(:,:) ) * z1_2dt_b … … 882 1030 END DO 883 1031 ELSE 1032 ! At this stage, ssha has been corrected: compute new depths at velocity points 1033 DO jj = 1, jpjm1 1034 DO ji = 1, jpim1 ! NO Vector Opt. 1035 zsshu_a(ji,jj) = z1_2 * umask(ji,jj,1) * r1_e1e2u(ji,jj) & 1036 & * ( e1e2t(ji ,jj) * ssha(ji ,jj) & 1037 & + e1e2t(ji+1,jj) * ssha(ji+1,jj) ) 1038 zsshv_a(ji,jj) = z1_2 * vmask(ji,jj,1) * r1_e1e2v(ji,jj) & 1039 & * ( e1e2t(ji,jj ) * ssha(ji,jj ) & 1040 & + e1e2t(ji,jj+1) * ssha(ji,jj+1) ) 1041 END DO 1042 END DO 1043 CALL lbc_lnk_multi( zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions 1044 ! 884 1045 DO jk=1,jpkm1 885 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b886 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b1046 ua(:,:,jk) = ua(:,:,jk) + r1_hu_n(:,:) * ( ua_b(:,:) - ub_b(:,:) * hu_b(:,:) ) * z1_2dt_b 1047 va(:,:,jk) = va(:,:,jk) + r1_hv_n(:,:) * ( va_b(:,:) - vb_b(:,:) * hv_b(:,:) ) * z1_2dt_b 887 1048 END DO 888 1049 ! Save barotropic velocities not transport: 889 ua_b (:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - umask(:,:,1) )890 va_b (:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - vmask(:,:,1) )1050 ua_b(:,:) = ua_b(:,:) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) ) 1051 va_b(:,:) = va_b(:,:) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) ) 891 1052 ENDIF 892 1053 ! 893 1054 DO jk = 1, jpkm1 894 1055 ! Correct velocities: 895 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) *umask(:,:,jk)896 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) *vmask(:,:,jk)1056 un(:,:,jk) = ( un(:,:,jk) + un_adv(:,:) - un_b(:,:) ) * umask(:,:,jk) 1057 vn(:,:,jk) = ( vn(:,:,jk) + vn_adv(:,:) - vn_b(:,:) ) * vmask(:,:,jk) 897 1058 ! 898 1059 END DO 1060 ! 1061 CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current 1062 CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic i-current 899 1063 ! 900 1064 #if defined key_agrif 901 1065 ! Save time integrated fluxes during child grid integration 902 ! (used to update coarse grid transports) 903 ! Useless with 2nd order momentum schemes 904 ! 905 IF ( (.NOT.Agrif_Root()).AND.(ln_bt_fw) ) THEN 906 IF ( Agrif_NbStepint().EQ.0 ) THEN 907 ub2_i_b(:,:) = 0.e0 908 vb2_i_b(:,:) = 0.e0 1066 ! (used to update coarse grid transports at next time step) 1067 ! 1068 IF( .NOT.Agrif_Root() .AND. ln_bt_fw ) THEN 1069 IF( Agrif_NbStepint() == 0 ) THEN 1070 ub2_i_b(:,:) = 0._wp 1071 vb2_i_b(:,:) = 0._wp 909 1072 END IF 910 1073 ! … … 913 1076 vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:) 914 1077 ENDIF 915 !916 !917 1078 #endif 918 !919 1079 ! !* write time-spliting arrays in the restart 920 IF(lrst_oce .AND.ln_bt_fw) CALL ts_rst( kt, 'WRITE' ) 921 ! 922 CALL wrk_dealloc( jpi, jpj, zsshp2_e, zhdiv ) 923 CALL wrk_dealloc( jpi, jpj, zu_trd, zv_trd, zun_e, zvn_e ) 924 CALL wrk_dealloc( jpi, jpj, zwx, zwy, zu_sum, zv_sum, zssh_frc, zu_frc, zv_frc ) 925 CALL wrk_dealloc( jpi, jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 926 CALL wrk_dealloc( jpi, jpj, zsshu_a, zsshv_a ) 927 CALL wrk_dealloc( jpi, jpj, zhf ) 928 ! 1080 IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' ) 1081 ! 1082 CALL wrk_dealloc( jpi,jpj, zsshp2_e, zhdiv ) 1083 CALL wrk_dealloc( jpi,jpj, zu_trd, zv_trd ) 1084 CALL wrk_dealloc( jpi,jpj, zwx, zwy, zssh_frc, zu_frc, zv_frc ) 1085 CALL wrk_dealloc( jpi,jpj, zhup2_e, zhvp2_e, zhust_e, zhvst_e ) 1086 CALL wrk_dealloc( jpi,jpj, zsshu_a, zsshv_a ) 1087 CALL wrk_dealloc( jpi,jpj, zhf ) 1088 IF( ln_wd ) CALL wrk_dealloc( jpi, jpj, zcpx, zcpy, wduflt1, wdvflt1 ) 1089 ! 1090 IF ( ln_diatmb ) THEN 1091 CALL iom_put( "baro_u" , un_b*umask(:,:,1)+zmdi*(1-umask(:,:,1 ) ) ) ! Barotropic U Velocity 1092 CALL iom_put( "baro_v" , vn_b*vmask(:,:,1)+zmdi*(1-vmask(:,:,1 ) ) ) ! Barotropic V Velocity 1093 ENDIF 929 1094 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 930 1095 ! 931 1096 END SUBROUTINE dyn_spg_ts 1097 932 1098 933 1099 SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2) … … 1008 1174 END SUBROUTINE ts_wgt 1009 1175 1176 1010 1177 SUBROUTINE ts_rst( kt, cdrw ) 1011 1178 !!--------------------------------------------------------------------- … … 1061 1228 END SUBROUTINE ts_rst 1062 1229 1063 SUBROUTINE dyn_spg_ts_init( kt ) 1230 1231 SUBROUTINE dyn_spg_ts_init 1064 1232 !!--------------------------------------------------------------------- 1065 1233 !! *** ROUTINE dyn_spg_ts_init *** … … 1067 1235 !! ** Purpose : Set time splitting options 1068 1236 !!---------------------------------------------------------------------- 1069 INTEGER , INTENT(in) :: kt ! ocean time-step 1070 ! 1071 INTEGER :: ji ,jj 1072 INTEGER :: ios ! Local integer output status for namelist read 1073 REAL(wp) :: zxr2, zyr2, zcmax 1074 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1075 !! 1076 NAMELIST/namsplit/ ln_bt_fw, ln_bt_av, ln_bt_nn_auto, & 1077 & nn_baro, rn_bt_cmax, nn_bt_flt 1237 INTEGER :: ji ,jj ! dummy loop indices 1238 REAL(wp) :: zxr2, zyr2, zcmax ! local scalar 1239 REAL(wp), POINTER, DIMENSION(:,:) :: zcu 1078 1240 !!---------------------------------------------------------------------- 1079 1241 ! 1080 REWIND( numnam_ref ) ! Namelist namsplit in reference namelist : time splitting parameters 1081 READ ( numnam_ref, namsplit, IOSTAT = ios, ERR = 901) 1082 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in reference namelist', lwp ) 1083 1084 REWIND( numnam_cfg ) ! Namelist namsplit in configuration namelist : time splitting parameters 1085 READ ( numnam_cfg, namsplit, IOSTAT = ios, ERR = 902 ) 1086 902 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsplit in configuration namelist', lwp ) 1087 IF(lwm) WRITE ( numond, namsplit ) 1088 ! 1089 ! ! Max courant number for ext. grav. waves 1090 ! 1091 CALL wrk_alloc( jpi, jpj, zcu ) 1092 ! 1093 IF (lk_vvl) THEN 1094 DO jj = 1, jpj 1095 DO ji =1, jpi 1096 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1097 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1098 zcu(ji,jj) = sqrt(grav*ht_0(ji,jj)*(zxr2 + zyr2) ) 1099 END DO 1242 ! Max courant number for ext. grav. waves 1243 ! 1244 CALL wrk_alloc( jpi,jpj, zcu ) 1245 ! 1246 DO jj = 1, jpj 1247 DO ji =1, jpi 1248 zxr2 = r1_e1t(ji,jj) * r1_e1t(ji,jj) 1249 zyr2 = r1_e2t(ji,jj) * r1_e2t(ji,jj) 1250 zcu(ji,jj) = SQRT( grav * ht_0(ji,jj) * (zxr2 + zyr2) ) 1100 1251 END DO 1101 ELSE 1102 DO jj = 1, jpj 1103 DO ji =1, jpi 1104 zxr2 = 1./(e1t(ji,jj)*e1t(ji,jj)) 1105 zyr2 = 1./(e2t(ji,jj)*e2t(ji,jj)) 1106 zcu(ji,jj) = sqrt(grav*ht(ji,jj)*(zxr2 + zyr2) ) 1107 END DO 1108 END DO 1109 ENDIF 1110 1111 zcmax = MAXVAL(zcu(:,:)) 1252 END DO 1253 ! 1254 zcmax = MAXVAL( zcu(:,:) ) 1112 1255 IF( lk_mpp ) CALL mpp_max( zcmax ) 1113 1256 1114 1257 ! Estimate number of iterations to satisfy a max courant number= rn_bt_cmax 1115 IF (ln_bt_nn_auto)nn_baro = CEILING( rdt / rn_bt_cmax * zcmax)1258 IF( ln_bt_auto ) nn_baro = CEILING( rdt / rn_bt_cmax * zcmax) 1116 1259 1117 rdtbt = rdt / FLOAT(nn_baro)1260 rdtbt = rdt / REAL( nn_baro , wp ) 1118 1261 zcmax = zcmax * rdtbt 1119 1262 ! Print results … … 1121 1264 IF(lwp) WRITE(numout,*) 'dyn_spg_ts : split-explicit free surface' 1122 1265 IF(lwp) WRITE(numout,*) '~~~~~~~~~~' 1123 IF( ln_bt_ nn_auto ) THEN1124 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.true. Automatically set nn_baro '1266 IF( ln_bt_auto ) THEN 1267 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.true. Automatically set nn_baro ' 1125 1268 IF(lwp) WRITE(numout,*) ' Max. courant number allowed: ', rn_bt_cmax 1126 1269 ELSE 1127 IF(lwp) WRITE(numout,*) ' ln_ts_ nn_auto=.false.: Use nn_baro in namelist '1270 IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_baro in namelist ' 1128 1271 ENDIF 1129 1272 … … 1143 1286 #if defined key_agrif 1144 1287 ! Restrict the use of Agrif to the forward case only 1145 IF ((.NOT.ln_bt_fw ).AND.(.NOT.Agrif_Root()))CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' )1288 IF( .NOT.ln_bt_fw .AND. .NOT.Agrif_Root() ) CALL ctl_stop( 'AGRIF not implemented if ln_bt_fw=.FALSE.' ) 1146 1289 #endif 1147 1290 ! 1148 1291 IF(lwp) WRITE(numout,*) ' Time filter choice, nn_bt_flt: ', nn_bt_flt 1149 1292 SELECT CASE ( nn_bt_flt ) 1150 CASE( 0 ); IF(lwp) WRITE(numout,*) ' Dirac'1151 CASE( 1 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro'1152 CASE( 2 ); IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro'1153 1293 CASE( 0 ) ; IF(lwp) WRITE(numout,*) ' Dirac' 1294 CASE( 1 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = nn_baro' 1295 CASE( 2 ) ; IF(lwp) WRITE(numout,*) ' Boxcar: width = 2*nn_baro' 1296 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt: should 0,1,2' ) 1154 1297 END SELECT 1155 1298 ! … … 1159 1302 IF(lwp) WRITE(numout,*) ' Maximum Courant number is :', zcmax 1160 1303 ! 1161 IF ((.NOT.ln_bt_av).AND.(.NOT.ln_bt_fw)) THEN1304 IF( .NOT.ln_bt_av .AND. .NOT.ln_bt_fw ) THEN 1162 1305 CALL ctl_stop( 'dynspg_ts ERROR: No time averaging => only forward integration is possible' ) 1163 1306 ENDIF 1164 IF 1307 IF( zcmax>0.9_wp ) THEN 1165 1308 CALL ctl_stop( 'dynspg_ts ERROR: Maximum Courant number is greater than 0.9: Inc. nn_baro !' ) 1166 1309 ENDIF 1167 1310 ! 1168 CALL wrk_dealloc( jpi, jpj,zcu )1311 CALL wrk_dealloc( jpi,jpj, zcu ) 1169 1312 ! 1170 1313 END SUBROUTINE dyn_spg_ts_init 1171 1314 1172 #else1173 !!---------------------------------------------------------------------------1174 !! Default case : Empty module No split explicit free surface1175 !!---------------------------------------------------------------------------1176 CONTAINS1177 INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function1178 dyn_spg_ts_alloc = 01179 END FUNCTION dyn_spg_ts_alloc1180 SUBROUTINE dyn_spg_ts( kt ) ! Empty routine1181 INTEGER, INTENT(in) :: kt1182 WRITE(*,*) 'dyn_spg_ts: You should not have seen this print! error?', kt1183 END SUBROUTINE dyn_spg_ts1184 SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine1185 INTEGER , INTENT(in) :: kt ! ocean time-step1186 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag1187 WRITE(*,*) 'ts_rst : You should not have seen this print! error?', kt, cdrw1188 END SUBROUTINE ts_rst1189 SUBROUTINE dyn_spg_ts_init( kt ) ! Empty routine1190 INTEGER , INTENT(in) :: kt ! ocean time-step1191 WRITE(*,*) 'dyn_spg_ts_init : You should not have seen this print! error?', kt1192 END SUBROUTINE dyn_spg_ts_init1193 #endif1194 1195 1315 !!====================================================================== 1196 1316 END MODULE dynspg_ts 1197 1198 1199
Note: See TracChangeset
for help on using the changeset viewer.
