- Timestamp:
- 2016-12-01T11:30:29+01:00 (8 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdylib.F90
r6140 r7412 5 5 !!====================================================================== 6 6 !! History : 3.6 ! 2013 (D. Storkey) original code 7 !! 4.0 ! 2014 (T. Lovato) Generalize OBC structure 7 8 !!---------------------------------------------------------------------- 8 #if defined key_bdy9 !!----------------------------------------------------------------------10 !! 'key_bdy' : Unstructured Open Boundary Condition11 9 !!---------------------------------------------------------------------- 12 10 !! bdy_orlanski_2d … … 25 23 PRIVATE 26 24 27 PUBLIC bdy_orlanski_2d ! routine called where? 28 PUBLIC bdy_orlanski_3d ! routine called where? 25 PUBLIC bdy_frs, bdy_spe, bdy_nmn, bdy_orl 26 PUBLIC bdy_orlanski_2d 27 PUBLIC bdy_orlanski_3d 29 28 30 29 !!---------------------------------------------------------------------- 31 !! NEMO/OPA 3.3 , NEMO Consortium (2010)30 !! NEMO/OPA 4.0 , NEMO Consortium (2016) 32 31 !! $Id$ 33 32 !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) 34 33 !!---------------------------------------------------------------------- 35 34 CONTAINS 35 36 SUBROUTINE bdy_frs( idx, pta, dta ) 37 !!---------------------------------------------------------------------- 38 !! *** SUBROUTINE bdy_frs *** 39 !! 40 !! ** Purpose : Apply the Flow Relaxation Scheme for tracers at open boundaries. 41 !! 42 !! Reference : Engedahl H., 1995, Tellus, 365-382. 43 !!---------------------------------------------------------------------- 44 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 45 REAL(wp), DIMENSION(:,:), INTENT(in) :: dta ! OBC external data 46 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! tracer trend 47 !! 48 REAL(wp) :: zwgt ! boundary weight 49 INTEGER :: ib, ik, igrd ! dummy loop indices 50 INTEGER :: ii, ij ! 2D addresses 51 !!---------------------------------------------------------------------- 52 ! 53 IF( nn_timing == 1 ) CALL timing_start('bdy_frs') 54 ! 55 igrd = 1 ! Everything is at T-points here 56 DO ib = 1, idx%nblen(igrd) 57 DO ik = 1, jpkm1 58 ii = idx%nbi(ib,igrd) 59 ij = idx%nbj(ib,igrd) 60 zwgt = idx%nbw(ib,igrd) 61 pta(ii,ij,ik) = ( pta(ii,ij,ik) + zwgt * (dta(ib,ik) - pta(ii,ij,ik) ) ) * tmask(ii,ij,ik) 62 END DO 63 END DO 64 ! 65 IF( nn_timing == 1 ) CALL timing_stop('bdy_frs') 66 ! 67 END SUBROUTINE bdy_frs 68 69 SUBROUTINE bdy_spe( idx, pta, dta ) 70 !!---------------------------------------------------------------------- 71 !! *** SUBROUTINE bdy_spe *** 72 !! 73 !! ** Purpose : Apply a specified value for tracers at open boundaries. 74 !! 75 !!---------------------------------------------------------------------- 76 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 77 REAL(wp), DIMENSION(:,:), INTENT(in) :: dta ! OBC external data 78 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! tracer trend 79 !! 80 REAL(wp) :: zwgt ! boundary weight 81 INTEGER :: ib, ik, igrd ! dummy loop indices 82 INTEGER :: ii, ij ! 2D addresses 83 !!---------------------------------------------------------------------- 84 ! 85 IF( nn_timing == 1 ) CALL timing_start('bdy_spe') 86 ! 87 igrd = 1 ! Everything is at T-points here 88 DO ib = 1, idx%nblenrim(igrd) 89 ii = idx%nbi(ib,igrd) 90 ij = idx%nbj(ib,igrd) 91 DO ik = 1, jpkm1 92 pta(ii,ij,ik) = dta(ib,ik) * tmask(ii,ij,ik) 93 END DO 94 END DO 95 ! 96 IF( nn_timing == 1 ) CALL timing_stop('bdy_spe') 97 ! 98 END SUBROUTINE bdy_spe 99 100 SUBROUTINE bdy_orl( idx, ptb, pta, dta, ll_npo ) 101 !!---------------------------------------------------------------------- 102 !! *** SUBROUTINE bdy_orl *** 103 !! 104 !! ** Purpose : Apply Orlanski radiation for tracers at open boundaries. 105 !! This is a wrapper routine for bdy_orlanski_3d below 106 !! 107 !!---------------------------------------------------------------------- 108 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 109 REAL(wp), DIMENSION(:,:), INTENT(in) :: dta ! OBC external data 110 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: ptb ! before tracer field 111 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pta ! tracer trend 112 LOGICAL, INTENT(in) :: ll_npo ! switch for NPO version 113 !! 114 INTEGER :: igrd ! grid index 115 !!---------------------------------------------------------------------- 116 ! 117 IF( nn_timing == 1 ) CALL timing_start('bdy_orl') 118 ! 119 igrd = 1 ! Everything is at T-points here 120 ! 121 CALL bdy_orlanski_3d( idx, igrd, ptb(:,:,:), pta(:,:,:), dta, ll_npo ) 122 ! 123 IF( nn_timing == 1 ) CALL timing_stop('bdy_orl') 124 ! 125 END SUBROUTINE bdy_orl 36 126 37 127 SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext, ll_npo ) … … 355 445 END SUBROUTINE bdy_orlanski_3d 356 446 357 358 #else 359 !!---------------------------------------------------------------------- 360 !! Dummy module NO Unstruct Open Boundary Conditions 361 !!---------------------------------------------------------------------- 362 CONTAINS 363 SUBROUTINE bdy_orlanski_2d( idx, igrd, phib, phia, phi_ext ) ! Empty routine 364 WRITE(*,*) 'bdy_orlanski_2d: You should not have seen this print! error?', kt 365 END SUBROUTINE bdy_orlanski_2d 366 SUBROUTINE bdy_orlanski_3d( idx, igrd, phib, phia, phi_ext ) ! Empty routine 367 WRITE(*,*) 'bdy_orlanski_3d: You should not have seen this print! error?', kt 368 END SUBROUTINE bdy_orlanski_3d 369 #endif 447 SUBROUTINE bdy_nmn( idx, igrd, phia ) 448 !!---------------------------------------------------------------------- 449 !! *** SUBROUTINE bdy_nmn *** 450 !! 451 !! ** Purpose : Duplicate the value at open boundaries, zero gradient. 452 !! 453 !!---------------------------------------------------------------------- 454 INTEGER, INTENT(in) :: igrd ! grid index 455 REAL(wp), DIMENSION(:,:,:), INTENT(inout) :: phia ! model after 3D field (to be updated) 456 TYPE(OBC_INDEX), INTENT(in) :: idx ! OBC indices 457 !! 458 REAL(wp) :: zcoef, zcoef1, zcoef2 459 REAL(wp), POINTER, DIMENSION(:,:,:) :: pmask ! land/sea mask for field 460 REAL(wp), POINTER, DIMENSION(:,:) :: bdypmask ! land/sea mask for field 461 INTEGER :: ib, ik ! dummy loop indices 462 INTEGER :: ii, ij, ip, jp ! 2D addresses 463 !!---------------------------------------------------------------------- 464 !! 465 IF( nn_timing == 1 ) CALL timing_start('bdy_nmn') 466 ! 467 SELECT CASE(igrd) 468 CASE(1) 469 pmask => tmask(:,:,:) 470 bdypmask => bdytmask(:,:) 471 CASE(2) 472 pmask => umask(:,:,:) 473 bdypmask => bdyumask(:,:) 474 CASE(3) 475 pmask => vmask(:,:,:) 476 bdypmask => bdyvmask(:,:) 477 CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' ) 478 END SELECT 479 DO ib = 1, idx%nblenrim(igrd) 480 ii = idx%nbi(ib,igrd) 481 ij = idx%nbj(ib,igrd) 482 DO ik = 1, jpkm1 483 ! search the sense of the gradient 484 zcoef1 = bdypmask(ii-1,ij )*pmask(ii-1,ij,ik) + bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) 485 zcoef2 = bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik) + bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) 486 IF ( nint(zcoef1+zcoef2) == 0) THEN 487 ! corner **** we probably only want to set the tangentail component for the dynamics here 488 zcoef = pmask(ii-1,ij,ik) + pmask(ii+1,ij,ik) + pmask(ii,ij-1,ik) + pmask(ii,ij+1,ik) 489 IF (zcoef > .5_wp) THEN ! Only set none isolated points. 490 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik) + & 491 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik) + & 492 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik) + & 493 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik) 494 phia(ii,ij,ik) = ( phia(ii,ij,ik) / zcoef ) * pmask(ii,ij,ik) 495 ELSE 496 phia(ii,ij,ik) = phia(ii,ij ,ik) * pmask(ii,ij ,ik) 497 ENDIF 498 ELSEIF ( nint(zcoef1+zcoef2) == 2) THEN 499 ! oblique corner **** we probably only want to set the normal component for the dynamics here 500 zcoef = pmask(ii-1,ij,ik)*bdypmask(ii-1,ij ) + pmask(ii+1,ij,ik)*bdypmask(ii+1,ij ) + & 501 & pmask(ii,ij-1,ik)*bdypmask(ii,ij -1 ) + pmask(ii,ij+1,ik)*bdypmask(ii,ij+1 ) 502 phia(ii,ij,ik) = phia(ii-1,ij ,ik) * pmask(ii-1,ij ,ik)*bdypmask(ii-1,ij ) + & 503 & phia(ii+1,ij ,ik) * pmask(ii+1,ij ,ik)*bdypmask(ii+1,ij ) + & 504 & phia(ii ,ij-1,ik) * pmask(ii ,ij-1,ik)*bdypmask(ii,ij -1 ) + & 505 & phia(ii ,ij+1,ik) * pmask(ii ,ij+1,ik)*bdypmask(ii,ij+1 ) 506 507 phia(ii,ij,ik) = ( phia(ii,ij,ik) / MAX(1._wp, zcoef) ) * pmask(ii,ij,ik) 508 ELSE 509 ip = nint(bdypmask(ii+1,ij )*pmask(ii+1,ij,ik) - bdypmask(ii-1,ij )*pmask(ii-1,ij,ik)) 510 jp = nint(bdypmask(ii ,ij+1)*pmask(ii,ij+1,ik) - bdypmask(ii ,ij-1)*pmask(ii,ij-1,ik)) 511 phia(ii,ij,ik) = phia(ii+ip,ij+jp,ik) * pmask(ii+ip,ij+jp,ik) 512 ENDIF 513 END DO 514 END DO 515 ! 516 IF( nn_timing == 1 ) CALL timing_stop('bdy_nmn') 517 ! 518 END SUBROUTINE bdy_nmn 370 519 371 520 !!======================================================================
Note: See TracChangeset
for help on using the changeset viewer.