Changeset 367 for trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
- Timestamp:
- 2005-12-28T10:25:10+01:00 (18 years ago)
- File:
-
- 1 edited
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90 (modified) (20 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90
r363 r367 17 17 USE phycst ! physical constants 18 18 USE ocesbc ! ocean surface boundary condition 19 USE obcdta ! open boundary condition data 20 USE obcfla ! Flather open boundary condition 19 21 USE dynvor ! vorticity term 20 22 USE obc_oce ! Lateral open boundary condition … … 22 24 USE lbclnk ! ocean lateral boundary conditions (or mpp link) 23 25 USE prtctl ! Print control 26 USE dynspg_oce ! surface pressure gradient variables 24 27 USE in_out_manager ! I/O manager 25 28 … … 29 32 !! * Accessibility 30 33 PUBLIC dyn_spg_ts ! routine called by step.F90 31 32 !! * Module variables33 REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables averaged over the barotropic loop34 sshn_b, sshb_b, & ! sea surface heigth (now, before)35 un_b , vn_b ! vertically integrated horizontal velocities (now)36 34 37 35 !! * Substitutions … … 61 59 !! surface gradient and the Coriolis force are updated within 62 60 !! the barotropic integration. 63 !! -2- Barotropic loop : updates of sea surface height ( zssha_e) and64 !! barotropic transports ( zua_e and zva_e) through barotropic61 !! -2- Barotropic loop : updates of sea surface height (ssha_e) and 62 !! barotropic transports (ua_e and va_e) through barotropic 65 63 !! momentum and continuity integration. Barotropic former 66 64 !! variables are time averaging over the full barotropic cycle … … 97 95 zssha_b, zua_b, zva_b, & ! " " 98 96 zsshb_e, zub_e, zvb_e, & ! " " 99 zsshn_e, zun_e, zvn_e, & ! " " 100 zssha_e, zua_e, zva_e ! " " 97 zun_e, zvn_e ! " " 101 98 REAL(wp), DIMENSION(jpi,jpj),SAVE :: & 102 99 ztnw, ztne, ztsw, ztse … … 105 102 ! Arrays initialization 106 103 ! --------------------- 107 zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0 ; zua_e(:,:) = 0.e0108 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0 ; zva_e(:,:) = 0.e0104 zua_b(:,:) = 0.e0 ; zub_e(:,:) = 0.e0 ; zun_e(:,:) = 0.e0 105 zva_b(:,:) = 0.e0 ; zvb_e(:,:) = 0.e0 ; zvn_e(:,:) = 0.e0 109 106 zhdiv(:,:) = 0.e0 110 107 … … 138 135 ENDIF 139 136 ENDIF 140 zssha_e(:,:) = sshn(:,:)141 zua_e (:,:)= un_b(:,:)142 zva_e (:,:)= vn_b(:,:)137 ssha_e(:,:) = sshn(:,:) 138 ua_e(:,:) = un_b(:,:) 139 va_e(:,:) = vn_b(:,:) 143 140 144 141 IF( ln_dynvor_een ) THEN … … 278 275 ! variables for the barotropic equations 279 276 zsshb_e(:,:) = sshn_b(:,:) ! (barotropic) sea surface height (before and now) 280 zsshn_e(:,:) = sshn_b(:,:) 281 zub_e(:,:) = un_b(:,:) ! barotropic transports issued from the barotropic equations (before and now) 282 zvb_e(:,:) = vn_b(:,:) 283 zun_e(:,:) = un_b(:,:) 284 zvn_e(:,:) = vn_b(:,:) 285 zssha_b(:,:) = sshn(:,:) ! time averaged variables over all sub-timesteps 286 zua_b(:,:) = un_b(:,:) 287 zva_b(:,:) = vn_b(:,:) 277 sshn_e (:,:) = sshn_b(:,:) 278 zub_e (:,:) = un_b (:,:) ! barotropic transports issued from the barotropic equations (before and now) 279 zvb_e (:,:) = vn_b (:,:) 280 zun_e (:,:) = un_b (:,:) 281 zvn_e (:,:) = vn_b (:,:) 282 zssha_b(:,:) = sshn (:,:) ! time averaged variables over all sub-timesteps 283 zua_b (:,:) = un_b (:,:) 284 zva_b (:,:) = vn_b (:,:) 285 286 ! set ssh corrections to 0 287 ! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop 288 #if defined key_obc 289 IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0 290 IF( lp_obc_west ) sshfow_b(:,:) = 0.e0 291 IF( lp_obc_south ) sshfos_b(:,:) = 0.e0 292 IF( lp_obc_north ) sshfon_b(:,:) = 0.e0 293 #endif 288 294 289 295 ! Barotropic integration over 2 baroclinic time steps … … 296 302 z2dt_e = 2. * rdtbt 297 303 IF ( jit == 1 ) z2dt_e = rdtbt 304 305 ! Time interpolation of open boundary condition data 306 IF( lk_obc ) CALL obc_dta_bt( kt, jit ) 298 307 299 308 ! Horizontal divergence of barotropic transports … … 312 321 ! open boundaries (div must be zero behind the open boundary) 313 322 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 314 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east315 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west316 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north317 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south323 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 324 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 325 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 326 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 318 327 #endif 319 328 … … 322 331 DO jj = 2, jpjm1 323 332 DO ji = fs_2, fs_jpim1 ! vector opt. 324 zssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) &333 ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * emp(ji,jj) & 325 334 & + zhdiv(ji,jj) ) ) * tmask(ji,jj,1) 326 335 END DO … … 336 345 DO ji = fs_2, fs_jpim1 ! vector opt. 337 346 ! surface pressure gradient 338 zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)339 zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)347 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 348 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 340 349 ! energy conserving formulation for planetary vorticity term 341 350 zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj) … … 346 355 zcvbt =-zfact2 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) 347 356 ! after transports 348 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)349 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)357 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 358 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 350 359 END DO 351 360 END DO … … 355 364 DO ji = fs_2, fs_jpim1 ! vector opt. 356 365 ! surface pressure gradient 357 zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)358 zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)366 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 367 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 359 368 ! enstrophy conserving formulation for planetary vorticity term 360 369 zy1 = zfact1 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) & … … 365 374 zcvbt = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) 366 375 ! after transports 367 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)368 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)376 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 377 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 369 378 END DO 370 379 END DO … … 375 384 DO ji = fs_2, fs_jpim1 ! vector opt. 376 385 ! surface pressure gradient 377 zspgu = -grav * ( zsshn_e(ji+1,jj) - zsshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj)378 zspgv = -grav * ( zsshn_e(ji,jj+1) - zsshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj)386 zspgu = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) * hu(ji,jj) / e1u(ji,jj) 387 zspgv = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) * hv(ji,jj) / e2v(ji,jj) 379 388 ! energy/enstrophy conserving formulation for planetary vorticity term 380 389 zcubt = + zfac25 / e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) & … … 383 392 & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) ) 384 393 ! after transports 385 zua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1)386 zva_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1)394 ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zcubt + zspgu + zua(ji,jj) ) ) * umask(ji,jj,1) 395 va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zcvbt + zspgv + zva(ji,jj) ) ) * vmask(ji,jj,1) 387 396 END DO 388 397 END DO 389 398 ENDIF 390 399 391 ! ... Boundary conditions on zua_e, zva_e, zssha_e 392 CALL lbc_lnk( zua_e, 'U', -1. ) 393 CALL lbc_lnk( zva_e, 'V', -1. ) 394 CALL lbc_lnk( zssha_e, 'T', 1. ) 400 ! Flather's boundary condition for the barotropic loop : 401 ! - Update sea surface height on each open boundary 402 ! - Correct the barotropic transports 403 IF( lk_obc ) CALL obc_fla_ts( kt ) 404 405 406 ! ... Boundary conditions on ua_e, va_e, ssha_e 407 CALL lbc_lnk( ua_e , 'U', -1. ) 408 CALL lbc_lnk( va_e , 'V', -1. ) 409 CALL lbc_lnk( ssha_e, 'T', 1. ) 395 410 396 411 ! temporal sum 397 412 !------------- 398 zssha_b(:,:) = zssha_b(:,:) + zssha_e(:,:)399 zua_b (:,:) = zua_b (:,:) + zua_e (:,:)400 zva_b (:,:) = zva_b (:,:) + zva_e (:,:)413 zssha_b(:,:) = zssha_b(:,:) + ssha_e(:,:) 414 zua_b (:,:) = zua_b (:,:) + ua_e (:,:) 415 zva_b (:,:) = zva_b (:,:) + va_e (:,:) 401 416 402 417 ! Time filter and swap of dynamics arrays 403 418 ! --------------------------------------- 404 419 IF( neuler == 0 .AND. kt == nit000 ) THEN ! Euler (forward) time stepping 405 zsshb_e(:,:) = zsshn_e(:,:)406 zub_e (:,:) = zun_e 407 zvb_e (:,:) = zvn_e 408 zsshn_e(:,:) = zssha_e(:,:)409 zun_e (:,:) = zua_e (:,:)410 zvn_e (:,:) = zva_e (:,:)420 zsshb_e(:,:) = sshn_e(:,:) 421 zub_e (:,:) = zun_e (:,:) 422 zvb_e (:,:) = zvn_e (:,:) 423 sshn_e (:,:) = ssha_e(:,:) 424 zun_e (:,:) = ua_e (:,:) 425 zvn_e (:,:) = va_e (:,:) 411 426 ELSE ! Asselin filtering 412 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + zssha_e(:,:) ) + atfp1 * zsshn_e(:,:)413 zub_e (:,:) = atfp * ( zub_e (:,:) + zua_e (:,:) ) + atfp1 * zun_e (:,:)414 zvb_e (:,:) = atfp * ( zvb_e (:,:) + zva_e (:,:) ) + atfp1 * zvn_e (:,:)415 zsshn_e(:,:) = zssha_e(:,:)416 zun_e (:,:) = zua_e (:,:)417 zvn_e (:,:) = zva_e (:,:)427 zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:) 428 zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:) 429 zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:) 430 sshn_e (:,:) = ssha_e(:,:) 431 zun_e (:,:) = ua_e (:,:) 432 zvn_e (:,:) = va_e (:,:) 418 433 ENDIF 419 434 … … 426 441 zcoef = 1.e0 / ( FLOAT( icycle +1 ) ) 427 442 zssha_b(:,:) = zcoef * zssha_b(:,:) 428 zua_b (:,:) = zcoef * zua_b (:,:) 429 zva_b (:,:) = zcoef * zva_b (:,:) 443 zua_b (:,:) = zcoef * zua_b (:,:) 444 zva_b (:,:) = zcoef * zva_b (:,:) 445 #if defined key_obc 446 IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) 447 IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:) 448 IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:) 449 IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:) 450 #endif 430 451 431 452 … … 448 469 ! open boundaries (div must be zero behind the open boundary) 449 470 ! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column 450 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0! east451 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0! west471 IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1) = 0.e0 ! east 472 IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1) = 0.e0 ! west 452 473 IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north 453 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0! south474 IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1) = 0.e0 ! south 454 475 #endif 455 476 … … 460 481 461 482 ! ... Boundary conditions on sshn 462 CALL lbc_lnk( sshn, 'T', 1. )483 IF( .NOT. lk_obc ) CALL lbc_lnk( sshn, 'T', 1. ) 463 484 464 485
Note: See TracChangeset
for help on using the changeset viewer.
