Changeset 1438 for trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2009-05-11T16:34:47+02:00 (15 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (27 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r1241 r1438 1 1 MODULE dynspg_ts 2 2 !!====================================================================== 3 !! *** MODULE dynspg_ts *** 4 !! Ocean dynamics: surface pressure gradient trend 5 !!====================================================================== 6 !! History : 9.0 ! 04-12 (L. Bessieres, G. Madec) Original code 7 !! " " ! 05-11 (V. Garnier, G. Madec) optimization 8 !! 9.0 ! 06-08 (S. Masson) distributed restart using iom 9 !! " ! 08-01 (R. Benshila) change averaging method 10 !! " ! 07-07 (D. Storkey) calls to BDY routines 11 !!--------------------------------------------------------------------- 3 !! History : 1.0 ! 04-12 (L. Bessieres, G. Madec) Original code 4 !! - ! 05-11 (V. Garnier, G. Madec) optimization 5 !! - ! 06-08 (S. Masson) distributed restart using iom 6 !! 2.0 ! 07-07 (D. Storkey) calls to BDY routines 7 !! - ! 08-01 (R. Benshila) change averaging method 8 !!--------------------------------------------------------------------- 12 9 #if defined key_dynspg_ts || defined key_esopa 13 10 !!---------------------------------------------------------------------- … … 19 16 !! ts_rst : read/write the time-splitting restart fields in the ocean restart file 20 17 !!---------------------------------------------------------------------- 21 !! * Modules used22 18 USE oce ! ocean dynamics and tracers 23 19 USE dom_oce ! ocean space and time domain … … 49 45 PUBLIC ts_rst ! routine called by istate.F90 50 46 51 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne, & ! triad of coriolis parameter 52 & ftsw, ftse ! (only used with een vorticity scheme) 53 47 REAL(wp), DIMENSION(jpi,jpj) :: ftnw, ftne ! triad of coriolis parameter 48 REAL(wp), DIMENSION(jpi,jpj) :: ftsw, ftse ! (only used with een vorticity scheme) 54 49 55 50 !! * Substitutions 56 51 # include "domzgr_substitute.h90" 57 52 # include "vectopt_loop_substitute.h90" 58 !!---------------------------------------------------------------------- 59 !! OPA 9.0 , LOCEAN-IPSL (2005)53 !!------------------------------------------------------------------------- 54 !! NEMO/OPA 3.2 , LOCEAN-IPSL (2009) 60 55 !! $Id$ 61 !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt62 !!---------------------------------------------------------------------- 56 !! Software is governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 57 !!------------------------------------------------------------------------- 63 58 64 59 CONTAINS … … 85 80 !! (= 2 * baroclinic time step) and saved in zsshX_b, zuX_b 86 81 !! and zvX_b (X specifying after, now or before). 87 !! -3- Update of sea surface height from time averaged barotropic 88 !! variables. 89 !! - apply lateral boundary conditions on sshn. 90 !! -4- The new general trend becomes : 91 !! ua = ua - sum_k(ua)/H + ( zua_b - sum_k(ub) )/H 82 !! -3- The new general trend becomes : 83 !! ua = ua - sum_k(ua)/H + ( zua_e - sum_k(ub) )/H 92 84 !! 93 85 !! ** Action : - Update (ua,va) with the surf. pressure gradient trend … … 96 88 !!--------------------------------------------------------------------- 97 89 INTEGER, INTENT( in ) :: kt ! ocean time-step index 98 99 !! * Local declarations 90 !! 100 91 INTEGER :: ji, jj, jk, jit ! dummy loop indices 101 92 INTEGER :: icycle ! temporary scalar … … 107 98 zcu, zcv, zwx, zwy, zhdiv, & ! temporary arrays 108 99 zua, zva, zub, zvb, & ! " " 109 z ssha_b, zua_b, zva_b, & ! " "110 zub_e, zvb_e, 111 zu n_e, zvn_e ! " "100 zua_e, zva_e, zssha_e, & ! " " 101 zub_e, zvb_e, zsshb_e, & ! " " 102 zu_sum, zv_sum, zssh_sum 112 103 !! Variable volume 113 104 REAL(wp), DIMENSION(jpi,jpj) :: & 114 105 zspgu_1, zspgv_1, zsshun_e, zsshvn_e ! 2D workspace 115 REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfse3un_e, zfse3vn_e ! 3D workspace116 106 !!---------------------------------------------------------------------- 117 107 118 108 ! Arrays initialization 119 109 ! --------------------- 120 zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0121 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0122 110 zhdiv(:,:) = 0.e0 123 111 … … 133 121 ! ! sshb, sshn, sshb_b, sshn_b, un_b, vn_b 134 122 135 ssha_e(:,:) = sshn(:,:) 136 ua_e(:,:) = un_b(:,:) 137 va_e(:,:) = vn_b(:,:) 138 hu_e(:,:) = hu(:,:) 139 hv_e(:,:) = hv(:,:) 140 123 zssha_e(:,:) = sshn(:,:) 124 zua_e (:,:) = un_e(:,:) 125 zva_e (:,:) = vn_e(:,:) 126 hu_e (:,:) = hu (:,:) 127 hv_e (:,:) = hv (:,:) 141 128 IF( ln_dynvor_een ) THEN 142 129 ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0 … … 170 157 zspgv_1(ji,jj) = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) & 171 158 & - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) / e2v(ji,jj) 172 END DO 159 END DO 173 160 END DO 174 161 … … 193 180 zfact2 = 0.5 * 0.5 194 181 zraur = 1. / rauw ! 1 / volumic mass of pure water 195 182 196 183 ! ----------------------------------------------------------------------------- 197 184 ! Phase 1 : Coupling between general trend and barotropic estimates (1st step) … … 215 202 zva(ji,1) = zva(ji,1) + fse3v(ji,1,jk) * vmask(ji,1,jk) * va(ji,1,jk) 216 203 ! ! Vertically integrated transports (before) 217 zub(ji,1) = zub(ji,1) + fse3u (ji,1,jk) * ub(ji,1,jk)218 zvb(ji,1) = zvb(ji,1) + fse3v (ji,1,jk) * vb(ji,1,jk)204 zub(ji,1) = zub(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk) 205 zvb(ji,1) = zvb(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk) 219 206 ! ! Planetary vorticity transport fluxes (now) 220 207 zwx(ji,1) = zwx(ji,1) + e2u(ji,1) * fse3u(ji,1,jk) * un(ji,1,jk) … … 228 215 zva(:,:) = zva(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk) 229 216 ! ! Vertically integrated transports (before) 230 zub(:,:) = zub(:,:) + fse3u (:,:,jk) * ub(:,:,jk)231 zvb(:,:) = zvb(:,:) + fse3v (:,:,jk) * vb(:,:,jk)217 zub(:,:) = zub(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk) 218 zvb(:,:) = zvb(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk) 232 219 ! ! Planetary vorticity (now) 233 220 zwx(:,:) = zwx(:,:) + e2u(:,:) * fse3u(:,:,jk) * un(:,:,jk) … … 304 291 !---------------- 305 292 ! Number of iteration of the barotropic loop 306 icycle = 3 * nn_baro / 2293 icycle = 2 * nn_baro + 1 307 294 308 295 ! variables for the barotropic equations 309 sshb_e (:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now) 310 sshn_e (:,:) = sshn_b(:,:) 311 zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now) 312 zvb_e (:,:) = vn_b (:,:) 313 zun_e (:,:) = un_b (:,:) 314 zvn_e (:,:) = vn_b (:,:) 315 zssha_b(:,:) = 0.e0 316 zua_b (:,:) = 0.e0 317 zva_b (:,:) = 0.e0 318 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 319 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 320 IF( lk_vvl ) THEN 321 zsshun_e(:,:) = sshu (:,:) ! (barotropic) sea surface height (now) at u-point 322 zsshvn_e(:,:) = sshv (:,:) ! (barotropic) sea surface height (now) at v-point 323 zfse3un_e(:,:,:) = fse3u(:,:,:) ! (barotropic) scale factors at u-point 324 zfse3un_e(:,:,:) = fse3v(:,:,:) ! (barotropic) scale factors at v-point 296 zu_sum (:,:) = 0.e0 297 zv_sum (:,:) = 0.e0 298 zssh_sum(:,:) = 0.e0 299 hu_e (:,:) = hu (:,:) ! (barotropic) ocean depth at u-point 300 hv_e (:,:) = hv (:,:) ! (barotropic) ocean depth at v-point 301 zsshb_e (:,:) = sshn_e(:,:) ! (barotropic) sea surface height (before and now) 302 ! vertical sum 303 un_e (:,:) = 0.e0 304 vn_e (:,:) = 0.e0 305 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 306 DO jk = 1, jpkm1 307 DO ji = 1, jpij 308 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 309 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 310 END DO 311 END DO 312 ELSE ! No vector opt. 313 DO jk = 1, jpkm1 314 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 315 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 316 END DO 325 317 ENDIF 318 zub_e (:,:) = un_e(:,:) 319 zvb_e (:,:) = vn_e(:,:) 320 326 321 327 322 ! set ssh corrections to 0 … … 352 347 DO jj = 2, jpjm1 353 348 DO ji = fs_2, fs_jpim1 ! vector opt. 354 zhdiv(ji,jj) = ( e2u(ji ,jj ) * zun_e(ji ,jj) &355 & -e2u(ji-1,jj ) * zun_e(ji-1,jj) &356 & +e1v(ji ,jj ) * zvn_e(ji ,jj) &357 & -e1v(ji ,jj-1) * zvn_e(ji ,jj-1) ) &349 zhdiv(ji,jj) = ( e2u(ji ,jj ) * un_e(ji ,jj) & 350 & -e2u(ji-1,jj ) * un_e(ji-1,jj) & 351 & +e1v(ji ,jj ) * vn_e(ji ,jj) & 352 & -e1v(ji ,jj-1) * vn_e(ji ,jj-1) ) & 358 353 & / (e1t(ji,jj)*e2t(ji,jj)) 359 354 END DO … … 370 365 371 366 #if defined key_bdy 372 373 374 375 376 367 DO jj = 1, jpj 368 DO ji = 1, jpi 369 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj) 370 END DO 371 END DO 377 372 #endif 378 373 … … 381 376 DO jj = 2, jpjm1 382 377 DO ji = fs_2, fs_jpim1 ! vector opt. 383 ssha_e(ji,jj) = (sshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) &384 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)378 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 379 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 385 380 END DO 386 381 END DO … … 388 383 ! evolution of the barotropic transport ( following the vorticity scheme used) 389 384 ! ---------------------------------------------------------------------------- 390 zwx(:,:) = e2u(:,:) * zun_e(:,:)391 zwy(:,:) = e1v(:,:) * zvn_e(:,:)385 zwx(:,:) = e2u(:,:) * un_e(:,:) 386 zwy(:,:) = e1v(:,:) * vn_e(:,:) 392 387 393 388 IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme … … 396 391 ! surface pressure gradient 397 392 IF( lk_vvl) THEN 398 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &399 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)400 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &401 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)393 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 394 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 395 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 396 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 402 397 ELSE 403 398 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 412 407 zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 413 408 ! after transports 414 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)415 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)409 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 410 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 416 411 END DO 417 412 END DO … … 422 417 ! surface pressure gradient 423 418 IF( lk_vvl) THEN 424 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &425 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)426 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &427 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)419 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 420 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 421 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 422 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 428 423 ELSE 429 424 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 432 427 ! enstrophy conserving formulation for planetary vorticity term 433 428 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & 434 429 + zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj) 435 430 zx1 =-zfact1 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) & 436 431 + zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj) 437 432 zcubt = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) 438 433 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 439 434 ! after transports 440 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)441 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)435 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 436 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 442 437 END DO 443 438 END DO … … 449 444 ! surface pressure gradient 450 445 IF( lk_vvl) THEN 451 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj) &452 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hu_e(ji,jj) / e1u(ji,jj)453 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji,jj+1) &454 & - ( rhd(ji,jj,1) + 1 ) * sshn_e(ji,jj) ) * hv_e(ji,jj) / e2v(ji,jj)446 zspgu = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) & 447 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj) 448 zspgv = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) & 449 & - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj) 455 450 ELSE 456 451 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) … … 463 458 & + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) 464 459 ! after transports 465 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)466 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)460 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 461 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 467 462 END DO 468 463 END DO … … 477 472 478 473 ! ... Boundary conditions on ua_e, va_e, ssha_e 479 CALL lbc_lnk( ua_e , 'U', -1. )480 CALL lbc_lnk( va_e , 'V', -1. )481 CALL lbc_lnk( ssha_e, 'T', 1. )474 CALL lbc_lnk( zua_e , 'U', -1. ) 475 CALL lbc_lnk( zva_e , 'V', -1. ) 476 CALL lbc_lnk( zssha_e, 'T', 1. ) 482 477 483 478 ! temporal sum 484 479 !------------- 485 IF( jit >= nn_baro / 2 ) THEN 486 zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) 487 zua_b (:,:) = zua_b (:,:) + ua_e (:,:) 488 zva_b (:,:) = zva_b (:,:) + va_e (:,:) 489 ENDIF 480 zu_sum (:,:) = zu_sum (:,:) + zua_e (:,:) 481 zv_sum (:,:) = zv_sum (:,:) + zva_e (:,:) 482 zssh_sum(:,:) = zssh_sum(:,:) + zssha_e(:,:) 490 483 491 484 ! Time filter and swap of dynamics arrays 492 485 ! --------------------------------------- 493 486 IF( jit == 1 ) THEN ! Euler (forward) time stepping 494 sshb_e (:,:) = sshn_e(:,:)495 zub_e (:,:) = zun_e(:,:)496 zvb_e (:,:) = zvn_e(:,:)497 sshn_e (:,:) = ssha_e(:,:)498 zun_e (:,:) =ua_e (:,:)499 zvn_e (:,:) =va_e (:,:)487 zsshb_e(:,:) = sshn_e (:,:) 488 zub_e (:,:) = un_e (:,:) 489 zvb_e (:,:) = vn_e (:,:) 490 sshn_e (:,:) = zssha_e(:,:) 491 un_e (:,:) = zua_e (:,:) 492 vn_e (:,:) = zva_e (:,:) 500 493 ELSE ! Asselin filtering 501 sshb_e (:,:) = atfp * ( sshb_e(:,:) +ssha_e(:,:) ) + atfp1 * sshn_e(:,:)502 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)503 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)504 sshn_e (:,:) = ssha_e(:,:)505 zun_e (:,:) =ua_e (:,:)506 zvn_e (:,:) =va_e (:,:)494 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * sshn_e(:,:) 495 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * un_e (:,:) 496 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * vn_e (:,:) 497 sshn_e (:,:) = zssha_e(:,:) 498 un_e (:,:) = zua_e (:,:) 499 vn_e (:,:) = zva_e (:,:) 507 500 ENDIF 508 501 509 IF( lk_vvl ) THEN ! Variable volume 502 IF( lk_vvl ) THEN ! Variable volume ! needed for BDY ??? 510 503 511 504 ! Sea surface elevation … … 528 521 529 522 ! Boundaries conditions 530 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; 531 532 ! Scale factors at before and after time step523 CALL lbc_lnk( zsshun_e, 'U', 1. ) ; CALL lbc_lnk( zsshvn_e, 'V', 1. ) 524 525 ! Ocean depth at U- and V-points 533 526 ! ------------------------------------------- 534 CALL dom_vvl_sf( zsshun_e, 'U', zfse3un_e ) ; CALL dom_vvl_sf( zsshvn_e, 'V', zfse3vn_e ) 535 536 ! Ocean depth at U- and V-points 537 hu_e(:,:) = 0.e0 538 hv_e(:,:) = 0.e0 539 540 DO jk = 1, jpk 541 hu_e(:,:) = hu_e(:,:) + zfse3un_e(:,:,jk) * umask(:,:,jk) 542 hv_e(:,:) = hv_e(:,:) + zfse3vn_e(:,:,jk) * vmask(:,:,jk) 543 END DO 544 545 ENDIF ! End variable volume case 546 527 hu_e(:,:) = hu_0(:,:) + zsshun_e(:,:) 528 hv_e(:,:) = hv_0(:,:) + zsshvn_e(:,:) 529 530 ! 531 ENDIF 547 532 ! ! ==================== ! 548 533 END DO ! end loop ! 549 534 ! ! ==================== ! 550 535 551 552 536 ! Time average of after barotropic variables 553 zcoef = 1.e0 / ( nn_baro + 1 ) 554 zssha_b(:,:) = zcoef * zssha_b(:,:) 555 zua_b (:,:) = zcoef * zua_b (:,:) 556 zva_b (:,:) = zcoef * zva_b (:,:) 537 zcoef = 1.e0 / ( 2 * nn_baro + 1 ) 538 un_e (:,:) = zcoef * zu_sum (:,:) 539 vn_e (:,:) = zcoef * zv_sum (:,:) 540 sshn_e(:,:) = zcoef * zssh_sum(:,:) 541 557 542 #if defined key_obc 558 543 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) … … 561 546 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:) 562 547 #endif 563 564 565 ! ---------------------------------------------------------------------------566 ! Phase 3 : Update sea surface height from time averaged barotropic variables567 ! ---------------------------------------------------------------------------568 569 570 ! Horizontal divergence of time averaged barotropic transports571 !-------------------------------------------------------------572 IF( .NOT. lk_vvl ) THEN573 DO jj = 2, jpjm1574 DO ji = fs_2, fs_jpim1 ! vector opt.575 zhdiv(ji,jj) = ( e2u(ji,jj) * un_b(ji,jj) - e2u(ji-1,jj ) * un_b(ji-1,jj ) &576 & +e1v(ji,jj) * vn_b(ji,jj) - e1v(ji ,jj-1) * vn_b(ji ,jj-1) ) &577 & / ( e1t(ji,jj) * e2t(ji,jj) )578 END DO579 END DO580 ENDIF581 582 #if defined key_obc && ! defined key_vvl583 ! open boundaries (div must be zero behind the open boundary)584 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column585 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east586 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west587 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north588 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south589 #endif590 591 #if defined key_bdy592 DO jj = 1, jpj593 DO ji = 1, jpi594 zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)595 END DO596 END DO597 #endif598 599 ! sea surface height600 !-------------------601 IF( lk_vvl ) THEN602 sshbb(:,:) = sshb(:,:)603 sshb (:,:) = sshn(:,:)604 sshn (:,:) = ssha(:,:)605 ELSE606 sshb (:,:) = sshn(:,:)607 sshn(:,:) = ( sshb_b(:,:) - z2dt_b * ( zraur * emp(:,:) + zhdiv(:,:) ) ) * tmask(:,:,1)608 ENDIF609 610 ! ... Boundary conditions on sshn611 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. )612 613 548 614 549 ! ----------------------------------------------------------------------------- 615 ! Phase 4. Coupling between general trend and barotropic estimates - (2nd step)550 ! Phase 3. Coupling between general trend and barotropic estimates - (2nd step) 616 551 ! ----------------------------------------------------------------------------- 617 552 618 ! Swap on time averaged barotropic variables 619 ! ------------------------------------------ 620 sshb_b(:,:) = sshn_b (:,:) 621 IF ( neuler == 0 .AND. kt == nit000 ) zssha_b(:,:) = sshn(:,:) 622 sshn_b(:,:) = zssha_b(:,:) 623 un_b (:,:) = zua_b (:,:) 624 vn_b (:,:) = zva_b (:,:) 625 553 554 626 555 ! add time averaged barotropic coriolis and surface pressure gradient 627 556 ! terms to the general momentum trend 628 557 ! -------------------------------------------------------------------- 629 558 DO jk=1,jpkm1 630 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( zua_b(:,:) - zub(:,:) ) / z2dt_b631 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( zva_b(:,:) - zvb(:,:) ) / z2dt_b559 ua(:,:,jk) = ua(:,:,jk) + hur(:,:) * ( un_e(:,:) - zub(:,:) ) / z2dt_b 560 va(:,:,jk) = va(:,:,jk) + hvr(:,:) * ( vn_e(:,:) - zvb(:,:) ) / z2dt_b 632 561 END DO 633 562 … … 636 565 IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' ) 637 566 638 ! print sum trends (used for debugging)639 IF( ln_ctl ) CALL prt_ctl( tab2d_1=sshn, clinfo1=' ssh : ', mask1=tmask )640 567 ! 641 568 END SUBROUTINE dyn_spg_ts … … 655 582 ! 656 583 IF( TRIM(cdrw) == 'READ' ) THEN 657 IF( iom_varid( numror, 'sshn ', ldstop = .FALSE. ) > 0 ) THEN658 CALL iom_get( numror, jpdom_autoglo, 'ssh b' , sshb(:,:) )659 CALL iom_get( numror, jpdom_autoglo, ' sshn' , sshn(:,:) )660 IF( neuler == 0 ) sshb(:,:) = sshn(:,:)584 IF( iom_varid( numror, 'sshn_e', ldstop = .FALSE. ) > 0 ) THEN 585 CALL iom_get( numror, jpdom_autoglo, 'sshn_e', sshn_e(:,:) ) ! free surface and 586 CALL iom_get( numror, jpdom_autoglo, 'un_e' , un_e (:,:) ) ! horizontal transports issued 587 CALL iom_get( numror, jpdom_autoglo, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 661 588 ELSE 662 IF( nn_rstssh == 1 ) THEN 663 sshb(:,:) = 0.e0 664 sshn(:,:) = 0.e0 665 ENDIF 666 ENDIF 667 IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN 668 CALL iom_get( numror, jpdom_autoglo, 'sshb_b', sshb_b(:,:) ) ! free surface issued 669 CALL iom_get( numror, jpdom_autoglo, 'sshn_b', sshn_b(:,:) ) ! from time-splitting loop 670 CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! horizontal transports issued 671 CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 672 IF( neuler == 0 ) sshb_b(:,:) = sshn_b(:,:) 673 ELSE 674 sshb_b(:,:) = sshb(:,:) 675 sshn_b(:,:) = sshn(:,:) 676 un_b (:,:) = 0.e0 677 vn_b (:,:) = 0.e0 589 sshn_e(:,:) = sshn(:,:) 590 un_e (:,:) = 0.e0 591 vn_e (:,:) = 0.e0 678 592 ! vertical sum 679 593 IF( lk_vopt_loop ) THEN ! vector opt., forced unroll 680 594 DO jk = 1, jpkm1 681 595 DO ji = 1, jpij 682 un_ b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)683 vn_ b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)596 un_e(ji,1) = un_e(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk) 597 vn_e(ji,1) = vn_e(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk) 684 598 END DO 685 599 END DO 686 600 ELSE ! No vector opt. 687 601 DO jk = 1, jpkm1 688 un_ b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)689 vn_ b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)602 un_e(:,:) = un_e(:,:) + fse3u(:,:,jk) * un(:,:,jk) 603 vn_e(:,:) = vn_e(:,:) + fse3v(:,:,jk) * vn(:,:,jk) 690 604 END DO 691 605 ENDIF 692 606 ENDIF 693 607 ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN 694 CALL iom_rstput( kt, nitrst, numrow, 'sshb' , sshb (:,:) ) 695 CALL iom_rstput( kt, nitrst, numrow, 'sshn' , sshn (:,:) ) 696 CALL iom_rstput( kt, nitrst, numrow, 'sshb_b', sshb_b(:,:) ) ! free surface issued 697 CALL iom_rstput( kt, nitrst, numrow, 'sshn_b', sshn_b(:,:) ) ! from barotropic loop 698 CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! horizontal transports issued 699 CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! from barotropic loop 608 CALL iom_rstput( kt, nitrst, numrow, 'sshn_e', sshn_e(:,:) ) ! free surface and 609 CALL iom_rstput( kt, nitrst, numrow, 'un_e' , un_e (:,:) ) ! horizontal transports issued 610 CALL iom_rstput( kt, nitrst, numrow, 'vn_e' , vn_e (:,:) ) ! from barotropic loop 700 611 ENDIF 701 612 !
Note: See TracChangeset
for help on using the changeset viewer.
