Changeset 3294 for trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2012-01-28T17:44:18+01:00 (12 years ago)
- File:
-
- 1 edited
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (29 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r2724 r3294 25 25 USE domvvl ! variable volume 26 26 USE zdfbfr ! bottom friction 27 USE obcdta ! open boundary condition data28 USE obcfla ! Flather open boundary condition29 27 USE dynvor ! vorticity term 30 28 USE obc_oce ! Lateral open boundary condition 31 29 USE obc_par ! open boundary condition parameters 32 USE bdy_oce ! unstructured open boundaries 33 USE bdy_par ! unstructured open boundaries 34 USE bdydta ! unstructured open boundaries 35 USE bdydyn ! unstructured open boundaries 36 USE bdytides ! tidal forcing at unstructured open boundaries. 30 USE obcdta ! open boundary condition data 31 USE obcfla ! Flather open boundary condition 32 USE bdy_par ! for lk_bdy 33 USE bdy_oce ! Lateral open boundary condition 34 USE bdydta ! open boundary condition data 35 USE bdydyn2d ! open boundary conditions on barotropic variables 36 USE sbctide 37 USE updtide 37 38 USE lib_mpp ! distributed memory computing library 38 39 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 41 42 USE iom ! IOM library 42 43 USE restart ! only for lrst_oce 43 USE zdf_oce 44 USE zdf_oce ! Vertical diffusion 45 USE wrk_nemo ! Memory Allocation 46 USE timing ! Timing 47 44 48 45 49 IMPLICIT NONE … … 107 111 !! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL 108 112 !!--------------------------------------------------------------------- 109 USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released110 USE wrk_nemo, ONLY: zsshun_e => wrk_2d_1 , zsshb_e => wrk_2d_2 , zhdiv => wrk_2d_3111 USE wrk_nemo, ONLY: zsshvn_e => wrk_2d_4 , zssh_sum => wrk_2d_5112 USE wrk_nemo, ONLY: zcu => wrk_2d_6 , zwx => wrk_2d_7 , zua => wrk_2d_8 , zbfru => wrk_2d_9113 USE wrk_nemo, ONLY: zcv => wrk_2d_10 , zwy => wrk_2d_11 , zva => wrk_2d_12 , zbfrv => wrk_2d_13114 USE wrk_nemo, ONLY: zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17115 USE wrk_nemo, ONLY: zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21116 113 ! 117 114 INTEGER, INTENT(in) :: kt ! ocean time-step index … … 119 116 INTEGER :: ji, jj, jk, jn ! dummy loop indices 120 117 INTEGER :: icycle ! local scalar 121 REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! local scalars 122 REAL(wp) :: z1_8, zx1, zy1 ! - - 123 REAL(wp) :: z1_4, zx2, zy2 ! - - 124 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 125 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 118 INTEGER :: ikbu, ikbv ! local scalar 119 REAL(wp) :: zraur, zcoef, z2dt_e, z1_2dt_b, z2dt_bf ! local scalars 120 REAL(wp) :: z1_8, zx1, zy1 ! - - 121 REAL(wp) :: z1_4, zx2, zy2 ! - - 122 REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - - 123 REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - - 124 REAL(wp) :: ua_btm, va_btm ! - - 125 ! 126 REAL(wp), POINTER, DIMENSION(:,:) :: zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv 127 REAL(wp), POINTER, DIMENSION(:,:) :: zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e 128 REAL(wp), POINTER, DIMENSION(:,:) :: zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum 126 129 !!---------------------------------------------------------------------- 127 128 IF( wrk_in_use(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & 129 & 11,12,13,14,15,16,17,18,19,20,21 ) ) THEN 130 CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' ) ; RETURN 131 ENDIF 132 130 ! 131 IF( nn_timing == 1 ) CALL timing_start('dyn_spg_ts') 132 ! 133 CALL wrk_alloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 134 CALL wrk_alloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 135 CALL wrk_alloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 136 ! 133 137 IF( kt == nit000 ) THEN !* initialisation 134 138 ! … … 147 151 hvr_e (:,:) = hvr (:,:) 148 152 IF( ln_dynvor_een ) THEN 149 ftne(1,:) = 0. e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0153 ftne(1,:) = 0._wp ; ftnw(1,:) = 0._wp ; ftse(1,:) = 0._wp ; ftsw(1,:) = 0._wp 150 154 DO jj = 2, jpj 151 155 DO ji = fs_2, jpi ! vector opt. 152 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3. 153 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3. 154 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3. 155 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3. 156 ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3._wp 157 ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3._wp 158 ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3._wp 159 ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3._wp 156 160 END DO 157 161 END DO … … 160 164 ENDIF 161 165 162 ! !* Local constant initialization 163 z2dt_b = 2.0 * rdt ! baroclinic time step 164 z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates 165 z1_4 = 0.5 * 0.5 166 zraur = 1. / rau0 ! 1 / volumic mass 167 ! 168 zhdiv(:,:) = 0.e0 ! barotropic divergence 169 zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F) 170 zv_sld = 0.e0 ; zv_asp = 0.e0 166 ! !* Local constant initialization 167 z1_2dt_b = 1._wp / ( 2.0_wp * rdt ) ! reciprocal of baroclinic time step 168 IF( neuler == 0 .AND. kt == nit000 ) z1_2dt_b = 1.0_wp / rdt ! reciprocal of baroclinic 169 ! time step (euler timestep) 170 z1_8 = 0.125_wp ! coefficient for vorticity estimates 171 z1_4 = 0.25_wp 172 zraur = 1._wp / rau0 ! 1 / volumic mass 173 ! 174 zhdiv(:,:) = 0._wp ! barotropic divergence 175 zu_sld = 0._wp ; zu_asp = 0._wp ! tides trends (lk_tide=F) 176 zv_sld = 0._wp ; zv_asp = 0._wp 177 178 IF( kt == nit000 .AND. neuler == 0) THEN ! for implicit bottom friction 179 z2dt_bf = rdt 180 ELSE 181 z2dt_bf = 2.0_wp * rdt 182 ENDIF 171 183 172 184 ! ----------------------------------------------------------------------------- … … 176 188 ! !* e3*d/dt(Ua), e3*Ub, e3*Vn (Vertically integrated) 177 189 ! ! -------------------------- 178 zua(:,:) = 0. e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0179 zva(:,:) = 0. e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0190 zua(:,:) = 0._wp ; zun(:,:) = 0._wp ; ub_b(:,:) = 0._wp 191 zva(:,:) = 0._wp ; zvn(:,:) = 0._wp ; vb_b(:,:) = 0._wp 180 192 ! 181 193 DO jk = 1, jpkm1 … … 195 207 ! 196 208 #if defined key_vvl 197 ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)*(1.+sshu_b(ji,jj)*muu(ji,jj,jk)))* ub(ji,jj,jk)198 vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)*(1.+sshv_b(ji,jj)*muv(ji,jj,jk)))* vb(ji,jj,jk)209 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_b(ji,jj,jk)* ub(ji,jj,jk) *umask(ji,jj,jk) 210 vb_b(ji,jj) = vb_b(ji,jj) + fse3v_b(ji,jj,jk)* vb(ji,jj,jk) *vmask(ji,jj,jk) 199 211 #else 200 212 ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk) … … 270 282 DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend 271 283 DO ji = fs_2, fs_jpim1 272 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)273 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)274 END DO284 zua(ji,jj) = zua(ji,jj) - zcu(ji,jj) 285 zva(ji,jj) = zva(ji,jj) - zcv(ji,jj) 286 END DO 275 287 END DO 276 288 … … 278 290 ! ! Remove barotropic contribution of bottom friction 279 291 ! ! from the barotropic transport trend 280 zcoef = -1. / z2dt_b 292 zcoef = -1._wp * z1_2dt_b 293 294 IF(ln_bfrimp) THEN 295 ! ! Remove the bottom stress trend from 3-D sea surface level gradient 296 ! ! and Coriolis forcing in case of 3D semi-implicit bottom friction 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 299 ikbu = mbku(ji,jj) 300 ikbv = mbkv(ji,jj) 301 ua_btm = zcu(ji,jj) * z2dt_bf * hur(ji,jj) * umask (ji,jj,ikbu) 302 va_btm = zcv(ji,jj) * z2dt_bf * hvr(ji,jj) * vmask (ji,jj,ikbv) 303 304 zua(ji,jj) = zua(ji,jj) - bfrua(ji,jj) * ua_btm 305 zva(ji,jj) = zva(ji,jj) - bfrva(ji,jj) * va_btm 306 END DO 307 END DO 308 309 ELSE 310 281 311 # if defined key_vectopt_loop 282 DO jj = 1, 1283 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)312 DO jj = 1, 1 313 DO ji = 1, jpij-jpi ! vector opt. (forced unrolling) 284 314 # else 285 DO jj = 2, jpjm1286 DO ji = 2, jpim1315 DO jj = 2, jpjm1 316 DO ji = 2, jpim1 287 317 # endif 288 318 ! Apply stability criteria for bottom friction 289 319 !RBbug for vvl and external mode we may need to use varying fse3 290 320 !!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx. 291 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 292 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 293 END DO 294 END DO 295 296 IF( lk_vvl ) THEN 297 DO jj = 2, jpjm1 298 DO ji = fs_2, fs_jpim1 ! vector opt. 299 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 300 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) ) 301 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 302 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) ) 303 END DO 304 END DO 305 ELSE 306 DO jj = 2, jpjm1 307 DO ji = fs_2, fs_jpim1 ! vector opt. 308 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 309 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 310 END DO 311 END DO 312 ENDIF 313 321 zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef ) 322 zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef ) 323 END DO 324 END DO 325 326 IF( lk_vvl ) THEN 327 DO jj = 2, jpjm1 328 DO ji = fs_2, fs_jpim1 ! vector opt. 329 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) & 330 & / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1._wp - umask(ji,jj,1) ) 331 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) & 332 & / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1._wp - vmask(ji,jj,1) ) 333 END DO 334 END DO 335 ELSE 336 DO jj = 2, jpjm1 337 DO ji = fs_2, fs_jpim1 ! vector opt. 338 zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj) 339 zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj) 340 END DO 341 END DO 342 ENDIF 343 END IF ! end (ln_bfrimp) 344 345 314 346 ! !* d/dt(Ua), Ub, Vn (Vertical mean velocity) 315 347 ! ! -------------------------- … … 318 350 ! 319 351 IF( lk_vvl ) THEN 320 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )321 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )352 ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 353 vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 322 354 ELSE 323 355 ub_b(:,:) = ub_b(:,:) * hur(:,:) … … 355 387 ! set ssh corrections to 0 356 388 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 357 IF( lp_obc_east ) sshfoe_b(:,:) = 0. e0358 IF( lp_obc_west ) sshfow_b(:,:) = 0. e0359 IF( lp_obc_south ) sshfos_b(:,:) = 0. e0360 IF( lp_obc_north ) sshfon_b(:,:) = 0. e0389 IF( lp_obc_east ) sshfoe_b(:,:) = 0._wp 390 IF( lp_obc_west ) sshfow_b(:,:) = 0._wp 391 IF( lp_obc_south ) sshfos_b(:,:) = 0._wp 392 IF( lp_obc_north ) sshfon_b(:,:) = 0._wp 361 393 #endif 362 394 … … 367 399 IF( jn == 1 ) z2dt_e = rdt / nn_baro 368 400 369 ! !* Update the forcing ( OBC,BDY and tides)401 ! !* Update the forcing (BDY and tides) 370 402 ! ! ------------------ 371 403 IF( lk_obc ) CALL obc_dta_bt ( kt, jn ) 372 IF( lk_bdy ) CALL bdy_dta_fla( kt, jn+1, icycle ) 404 IF( lk_bdy ) CALL bdy_dta ( kt, jit=jn, time_offset=+1 ) 405 IF ( ln_tide_pot ) CALL upd_tide( kt, jn ) 373 406 374 407 ! !* after ssh_e 375 408 ! ! ----------- 376 DO jj = 2, jpjm1 409 DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports 377 410 DO ji = fs_2, fs_jpim1 ! vector opt. 378 411 zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) & … … 386 419 ! ! OBC : zhdiv must be zero behind the open boundary 387 420 !! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 388 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0. e0! east389 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0. e0! west390 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0. e0! north391 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0. e0! south421 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0._wp ! east 422 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0._wp ! west 423 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0._wp ! north 424 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0._wp ! south 392 425 #endif 393 426 #if defined key_bdy … … 403 436 ! !* after barotropic velocities (vorticity scheme dependent) 404 437 ! ! --------------------------- 405 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) 438 zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport 406 439 zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:) 407 440 ! … … 419 452 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 420 453 ENDIF 454 ! add tidal astronomical forcing 455 IF ( ln_tide_pot ) THEN 456 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 457 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 458 ENDIF 421 459 ! energy conserving formulation for planetary vorticity term 422 460 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) … … 427 465 zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj) 428 466 ! after velocities with implicit bottom friction 429 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 430 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 431 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 432 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 467 468 IF( ln_bfrimp ) THEN ! implicit bottom friction 469 ! A new method to implement the implicit bottom friction. 470 ! H. Liu 471 ! Sept 2011 472 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 473 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 474 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 475 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 476 ! 477 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 478 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 479 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 480 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 481 ! 482 ELSE 483 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 484 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 485 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 486 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 487 ENDIF 433 488 END DO 434 489 END DO … … 447 502 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 448 503 ENDIF 504 ! add tidal astronomical forcing 505 IF ( ln_tide_pot ) THEN 506 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 507 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 508 ENDIF 449 509 ! enstrophy conserving formulation for planetary vorticity term 450 510 zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj) … … 453 513 zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj) 454 514 ! after velocities with implicit bottom friction 455 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 456 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 457 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 458 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 515 IF( ln_bfrimp ) THEN 516 ! A new method to implement the implicit bottom friction. 517 ! H. Liu 518 ! Sept 2011 519 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 520 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 521 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 522 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 523 ! 524 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 525 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 526 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 527 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 528 ! 529 ELSE 530 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 531 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 532 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 533 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 534 ENDIF 459 535 END DO 460 536 END DO … … 473 549 zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj) 474 550 ENDIF 551 ! add tidal astronomical forcing 552 IF ( ln_tide_pot ) THEN 553 zu_spg = zu_spg + grav * ( pot_astro(ji+1,jj) - pot_astro(ji,jj) ) / e1u(ji,jj) 554 zv_spg = zv_spg + grav * ( pot_astro(ji,jj+1) - pot_astro(ji,jj) ) / e2v(ji,jj) 555 ENDIF 475 556 ! energy/enstrophy conserving formulation for planetary vorticity term 476 557 zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) & … … 479 560 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj) 480 561 ! after velocities with implicit bottom friction 481 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) & 482 & / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 483 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) & 484 & / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 562 IF( ln_bfrimp ) THEN 563 ! A new method to implement the implicit bottom friction. 564 ! H. Liu 565 ! Sept 2011 566 ua_e(ji,jj) = umask(ji,jj,1) * ( zub_e(ji,jj) + & 567 & z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp ) & 568 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) ) 569 ua_e(ji,jj) = ( ua_e(ji,jj) + z2dt_e * zua(ji,jj) ) * umask(ji,jj,1) 570 ! 571 va_e(ji,jj) = vmask(ji,jj,1) * ( zvb_e(ji,jj) + & 572 & z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp ) & 573 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) ) 574 va_e(ji,jj) = ( va_e(ji,jj) + z2dt_e * zva(ji,jj) ) * vmask(ji,jj,1) 575 ! 576 ELSE 577 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj))) * umask(ji,jj,1) & 578 & / ( 1._wp - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) ) 579 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj))) * vmask(ji,jj,1) & 580 & / ( 1._wp - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) ) 581 ENDIF 485 582 END DO 486 583 END DO 487 584 ! 488 585 ENDIF 489 ! !* domain lateral boundary 490 ! ! ----------------------- 491 ! ! Flather's boundary condition for the barotropic loop : 492 ! ! - Update sea surface height on each open boundary 493 ! ! - Correct the velocity 494 586 ! !* domain lateral boundary 587 ! ! ----------------------- 588 589 ! OBC open boundaries 495 590 IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e ) 496 IF( lk_bdy .OR. ln_tides ) CALL bdy_dyn_fla( sshn_e ) 591 592 ! BDY open boundaries 593 #if defined key_bdy 594 pssh => sshn_e 595 phur => hur_e 596 phvr => hvr_e 597 pu2d => ua_e 598 pv2d => va_e 599 600 IF( lk_bdy ) CALL bdy_dyn2d( kt ) 601 #endif 602 497 603 ! 498 604 CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries … … 526 632 DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points 527 633 DO ji = 1, fs_jpim1 ! Vector opt. 528 zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &529 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &530 & +e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )531 zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &532 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &533 & +e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )634 zsshun_e(ji,jj) = 0.5_wp * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) & 635 & * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) & 636 & + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) ) 637 zsshvn_e(ji,jj) = 0.5_wp * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) & 638 & * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) & 639 & + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) ) 534 640 END DO 535 641 END DO … … 539 645 hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points 540 646 hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:) 541 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1. e0- umask(:,:,1) )542 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1. e0- vmask(:,:,1) )647 hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1._wp - umask(:,:,1) ) 648 hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1._wp - vmask(:,:,1) ) 543 649 ! 544 650 ENDIF … … 559 665 ! 560 666 ! !* Time average ==> after barotropic u, v, ssh 561 zcoef = 1. e0/ ( 2 * nn_baro + 1 )667 zcoef = 1._wp / ( 2 * nn_baro + 1 ) 562 668 zu_sum(:,:) = zcoef * zu_sum (:,:) 563 669 zv_sum(:,:) = zcoef * zv_sum (:,:) … … 565 671 ! !* update the general momentum trend 566 672 DO jk=1,jpkm1 567 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b568 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b673 ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) * z1_2dt_b 674 va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) * z1_2dt_b 569 675 END DO 570 676 un_b (:,:) = zu_sum(:,:) … … 575 681 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 576 682 ! 577 ! 578 IF( wrk_not_released(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, & 579 & 11,12,13,14,15,16,17,18,19,20,21) ) & 580 CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays') 683 CALL wrk_dealloc( jpi, jpj, zsshun_e, zsshvn_e, zsshb_e, zssh_sum, zhdiv ) 684 CALL wrk_dealloc( jpi, jpj, zua, zva, zun, zvn, zun_e, zvn_e, zub_e, zvb_e ) 685 CALL wrk_dealloc( jpi, jpj, zcu, zcv, zwx, zwy, zbfru, zbfrv, zu_sum, zv_sum ) 686 ! 687 IF( nn_timing == 1 ) CALL timing_stop('dyn_spg_ts') 581 688 ! 582 689 END SUBROUTINE dyn_spg_ts … … 600 707 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 601 708 ELSE 602 un_b (:,:) = 0. e0603 vn_b (:,:) = 0. e0709 un_b (:,:) = 0._wp 710 vn_b (:,:) = 0._wp 604 711 ! vertical sum 605 712 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll … … 622 729 ! Vertically integrated velocity (before) 623 730 IF (neuler/=0) THEN 624 ub_b (:,:) = 0. e0625 vb_b (:,:) = 0. e0731 ub_b (:,:) = 0._wp 732 vb_b (:,:) = 0._wp 626 733 627 734 ! vertical sum … … 641 748 642 749 IF( lk_vvl ) THEN 643 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1. e0- umask(:,:,1) )644 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1. e0- vmask(:,:,1) )750 ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1._wp - umask(:,:,1) ) 751 vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1._wp - vmask(:,:,1) ) 645 752 ELSE 646 753 ub_b(:,:) = ub_b(:,:) * hur(:,:)
Note: See TracChangeset
for help on using the changeset viewer.
