Changeset 911 for trunk/NEMO/OPA_SRC/DYN

Timestamp:

2008-04-28T11:31:32+02:00 (16 years ago)

Author:

ctlod

Message:

Implementation of the BDY package, see ticket: #126

Location:

trunk/NEMO/OPA_SRC/DYN

Files:

• divcur.F90 (modified) (5 diffs)
• dynnxt.F90 (modified) (4 diffs)
• dynspg_flt.F90 (modified) (4 diffs)
• dynspg_oce.F90 (modified) (1 diff)
• dynspg_ts.F90 (modified) (12 diffs)
• wzvmod.F90 (modified) (5 diffs)

trunk/NEMO/OPA_SRC/DYN/divcur.F90

@@ -14 +14 @@
    USE in_out_manager  ! I/O manager
    USE obc_oce         ! ocean lateral open boundary condition
+   USE bdy_oce        ! Unstructured open boundaries variables
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
 
@@ -77 +78 @@
       !!   8.2  !  00-03  (G. Madec)  no slip accurate
       !!   9.0  !  03-08  (G. Madec)  merged of cur and div, free form, F90
+      !!        !  05-01  (J. Chanut, A. Sellar) unstructured open boundaries
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -134 +136 @@
          ENDIF
 #endif
+#endif         
+#if defined key_bdy || key_bdy_tides
+         ! unstructured open boundaries (div must be zero behind the open boundary)
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
+           END DO
+         END DO
 #endif         
 #if defined key_agrif
@@ -291 +301 @@
       !!   8.1  !  97-08  (J.M. Molines)  Open boundaries
       !!   9.0  !  02-09  (G. Madec, E. Durand)  Free form, F90
+      !!        !  05-01  (J. Chanut) Unstructured open boundaries
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -344 +355 @@
 #endif
 #endif         
+#if defined key_bdy || key_bdy_tides
+         ! unstructured open boundaries (div must be zero behind the open boundary)
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             hdivn(ji,jj,jk)=hdivn(ji,jj,jk)*bdytmask(ji,jj)
+           END DO
+         END DO
+#endif        
 #if defined key_agrif
          if ( .NOT. AGRIF_Root() ) then

trunk/NEMO/OPA_SRC/DYN/dynnxt.F90

@@ -16 +16 @@
    USE obcdyn_bt       ! 2D open boundary condition for momentum (obc_dyn_bt routine)
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE bdy_oce         ! unstructured open boundary conditions
+   USE bdydta          ! unstructured open boundary conditions
+   USE bdydyn          ! unstructured open boundary conditions
    USE dynspg_oce      ! type of surface pressure gradient
    USE lbclnk          ! lateral boundary condition (or mpp link)
@@ -66 +69 @@
       !!        !  02-10  (C. Talandier, A-M. Treguier) Open boundary cond.
       !!   9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
+      !!    "   !  07-07  (D. Storkey) Calls to BDY routines. 
       !!----------------------------------------------------------------------
       !! * Arguments
@@ -186 +190 @@
             END DO
          END DO
+
+         IF( lk_vvl ) THEN
+            ! Unweight velocities prior to updating open boundaries.
+
+            DO jj = 1, jpj                      ! caution: don't use (:,:) for this loop
+               DO ji = 1, jpi                   ! it causes optimization problems on NEC in auto-tasking
+                  ua(ji,jj,jk) = ua(ji,jj,jk) / fse3u(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) / fse3v(ji,jj,jk)
+
+                  un(ji,jj,jk) = un(ji,jj,jk) / fse3u(ji,jj,jk)
+                  vn(ji,jj,jk) = vn(ji,jj,jk) / fse3v(ji,jj,jk)
+
+                  ub(ji,jj,jk) = ub(ji,jj,jk) / zfse3ub(ji,jj,jk)
+                  vb(ji,jj,jk) = vb(ji,jj,jk) / zfse3vb(ji,jj,jk)
+               END DO
+            END DO
+
+         ENDIF
+
 # if defined key_obc
          !                                             ! ===============
@@ -209 +232 @@
 
          IF ( ln_vol_cst ) CALL obc_vol( kt )
+
+      ENDIF
+
+      !                                                ! ===============
+      DO jk = 1, jpkm1                                 ! Horizontal slab
+         !                                             ! ===============
+# elif defined key_bdy || key_bdy_tides
+         !                                             ! ===============
+      END DO                                           !   End of slab
+      !                                                ! ===============
+      ! Update (ua,va) along open boundaries (for exp or ts options).
+      IF ( lk_dynspg_exp .or. lk_dynspg_ts ) THEN
+  
+         CALL bdy_dyn_frs( kt )
+
+         IF ( ln_bdy_fla ) THEN
+
+            ua_e(:,:)=0.0
+            va_e(:,:)=0.0
+
+            ! Set these variables for use in bdy_dyn_fla
+            hu_e(:,:) = hu(:,:)
+            hv_e(:,:) = hv(:,:)
+
+            DO jk = 1, jpkm1
+               !! Vertically integrated momentum trends
+               ua_e(:,:) = ua_e(:,:) + fse3u(:,:,jk) * umask(:,:,jk) * ua(:,:,jk)
+               va_e(:,:) = va_e(:,:) + fse3v(:,:,jk) * vmask(:,:,jk) * va(:,:,jk)
+            END DO
+
+            DO jk = 1 , jpkm1
+               ua(:,:,jk) = ua(:,:,jk) - ua_e(:,:) * hur(:,:)
+               va(:,:,jk) = va(:,:,jk) - va_e(:,:) * hvr(:,:)
+            END DO
+
+            CALL bdy_dta_bt( kt+1, 0)
+            CALL bdy_dyn_fla
+
+         ENDIF
+
+         DO jk = 1 , jpkm1
+            ua(:,:,jk) = ua(:,:,jk) + ua_e(:,:) * hur(:,:)
+            va(:,:,jk) = va(:,:,jk) + va_e(:,:) * hvr(:,:)
+         END DO
 
       ENDIF

trunk/NEMO/OPA_SRC/DYN/dynspg_flt.F90

@@ -11 +11 @@
    !!            " "  !  05-11  (V. Garnier) Surface pressure gradient organization
    !!            " "  !  06-07  (S. Masson)  distributed restart using iom
+   !!            " "  !  05-01  (J.Chanut, A.Sellar) Calls to BDY routines. 
    !!----------------------------------------------------------------------
 #if defined key_dynspg_flt   ||   defined key_esopa  
@@ -36 +37 @@
    USE obcdyn          ! ocean open boundary condition (obc_dyn routines)
    USE obcvol          ! ocean open boundary condition (obc_vol routines)
+   USE bdy_oce         ! Unstructured open boundaries condition
+   USE bdydyn          ! Unstructured open boundaries condition (bdy_dyn routine) 
+   USE bdyvol          ! Unstructured open boundaries condition (bdy_vol routine)
    USE cla_dynspg      ! cross land advection
    USE in_out_manager  ! I/O manager
@@ -249 +253 @@
       CALL obc_vol( kt )      ! Correction of the barotropic componant velocity to control the volume of the system
 #endif
+#if defined key_bdy
+      ! Update velocities on unstructured boundary using the Flow Relaxation Scheme
+      CALL bdy_dyn_frs( kt )
+
+      IF (ln_bdy_vol) THEN
+      ! Correction of the barotropic component velocity to control the volume of the system
+        CALL bdy_vol( kt )
+      ENDIF
+#endif
 #if defined key_agrif
       CALL Agrif_dyn( kt )    ! Update velocities on each coarse/fine interfaces 
@@ -377 +390 @@
             spgu(ji,jj) = z2dt * ztdgu * obcumask(ji,jj)
             spgv(ji,jj) = z2dt * ztdgv * obcvmask(ji,jj)
+#elif defined key_bdy
+            ! caution : grad D = 0 along open boundaries
+            ! Remark: The filtering force could be reduced here in the FRS zone
+            !         by multiplying spgu/spgv by (1-alpha) ??  
+            spgu(ji,jj) = z2dt * ztdgu * bdyumask(ji,jj)
+            spgv(ji,jj) = z2dt * ztdgv * bdyvmask(ji,jj)           
 #else
             spgu(ji,jj) = z2dt * ztdgu

trunk/NEMO/OPA_SRC/DYN/dynspg_oce.F90

@@ -46 +46 @@
          un_b  , vn_b                     ! vertically integrated horizontal velocities (now)
       REAL(wp), PUBLIC, DIMENSION(jpi,jpj) :: & ! variables of the explicit barotropic loop
-         sshn_e, ssha_e,               &  ! sea surface heigth (now,after)
-         ua_e  , va_e                     ! vertically integrated horizontal velocities (after)
+         sshb_e, sshn_e, ssha_e,       &  ! sea surface heigth (before,now,after)
+         ua_e  , va_e,                 &  ! vertically integrated horizontal velocities (after)
+         hu_e  , hv_e                     ! depth arrays for the barotropic solution
 #endif
 

trunk/NEMO/OPA_SRC/DYN/dynspg_ts.F90

@@ -7 +7 @@
    !!             " "  !  05-11  (V. Garnier, G. Madec)  optimization
    !!             9.0  !  06-08  (S. Masson)  distributed restart using iom
-   !!             9.0  !  08-01  (R. Benshila)  change averaging method
+   !!              "   !  08-01  (R. Benshila)  change averaging method
+   !!              "   !  07-07  (D. Storkey) calls to BDY routines
    !!---------------------------------------------------------------------
 #if defined key_dynspg_ts   ||   defined key_esopa
@@ -30 +31 @@
    USE obc_oce         ! Lateral open boundary condition
    USE obc_par         ! open boundary condition parameters
+   USE bdy_oce         ! unstructured open boundaries
+   USE bdy_par         ! unstructured open boundaries
+   USE bdydta          ! unstructured open boundaries
+   USE bdydyn          ! unstructured open boundaries
+   USE bdytides        ! tidal forcing at unstructured open boundaries.
    USE lib_mpp         ! distributed memory computing library
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
@@ -102 +108 @@
          zua, zva, zub, zvb,                &  !     "        "
          zssha_b, zua_b, zva_b,             &  !     "        "
-         zsshb_e, zub_e, zvb_e,             &  !     "        "
+         zub_e, zvb_e,                      &  !     "        "
          zun_e, zvn_e                          !     "        "
       !! Variable volume
       REAL(wp), DIMENSION(jpi,jpj)     ::   &
-         zspgu_1, zspgv_1, zsshun_e, zsshvn_e, hu_e, hv_e         ! 2D workspace
+         zspgu_1, zspgv_1, zsshun_e, zsshvn_e                     ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zfse3un_e, zfse3vn_e  ! 3D workspace
       !!----------------------------------------------------------------------
@@ -130 +136 @@
          ua_e(:,:)   = un_b(:,:)
          va_e(:,:)   = vn_b(:,:)
+         hu_e(:,:)   = hu(:,:)
+         hv_e(:,:)   = hv(:,:)
 
          IF( ln_dynvor_een ) THEN
@@ -300 +308 @@
 
       ! variables for the barotropic equations
-      zsshb_e(:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
+      sshb_e (:,:) = sshn_b(:,:)       ! (barotropic) sea surface height (before and now)
       sshn_e (:,:) = sshn_b(:,:)
       zub_e  (:,:) = un_b  (:,:)       ! barotropic transports issued from the barotropic equations (before and now)
@@ -309 +317 @@
       zua_b  (:,:) = 0.e0
       zva_b  (:,:) = 0.e0
+      hu_e   (:,:) = hu   (:,:)     ! (barotropic) ocean depth at u-point
+      hv_e   (:,:) = hv   (:,:)     ! (barotropic) ocean depth at v-point
       IF( lk_vvl ) THEN
          zsshun_e(:,:)    = sshu (:,:)     ! (barotropic) sea surface height (now) at u-point
          zsshvn_e(:,:)    = sshv (:,:)     ! (barotropic) sea surface height (now) at v-point
-         hu_e(:,:)        = hu   (:,:)     ! (barotropic) ocean depth at u-point
-         hv_e(:,:)        = hv   (:,:)     ! (barotropic) ocean depth at v-point
          zfse3un_e(:,:,:) = fse3u(:,:,:)   ! (barotropic) scale factors  at u-point
          zfse3un_e(:,:,:) = fse3v(:,:,:)   ! (barotropic) scale factors  at v-point
@@ -338 +346 @@
          ! Time interpolation of open boundary condition data
          IF( lk_obc )   CALL obc_dta_bt( kt, jit )
+         IF( lk_bdy .or. lk_bdy_tides)   CALL bdy_dta_bt( kt, jit+1 )
 
          ! Horizontal divergence of barotropic transports
          !--------------------------------------------------
+         zhdiv(:,:) = 0.e0
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -360 +370 @@
 #endif
 
+         IF( lk_bdy .or. lk_bdy_tides ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+               END DO
+            END DO
+         ENDIF
+
          ! Sea surface height from the barotropic system
          !----------------------------------------------
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
+               ssha_e(ji,jj) = ( sshb_e(ji,jj) - z2dt_e *  ( zraur * emp(ji,jj)  &
             &            +  zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
             END DO
@@ -457 +475 @@
          !         - Correct the barotropic transports
          IF( lk_obc )   CALL obc_fla_ts
-
+         IF( lk_bdy .or. lk_bdy_tides )   CALL bdy_dyn_fla
 
          ! ... Boundary conditions on ua_e, va_e, ssha_e
@@ -475 +493 @@
          ! ---------------------------------------
          IF( neuler == 0 .AND. jit == 1 ) THEN   ! Euler (forward) time stepping
-            zsshb_e(:,:) = sshn_e(:,:)
+            sshb_e (:,:) = sshn_e(:,:)
             zub_e  (:,:) = zun_e (:,:)
             zvb_e  (:,:) = zvn_e (:,:)
@@ -482 +500 @@
             zvn_e  (:,:) = va_e  (:,:)
          ELSE                                        ! Asselin filtering
-            zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
+            sshb_e (:,:) = atfp * ( sshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
             zub_e  (:,:) = atfp * ( zub_e  (:,:) + ua_e  (:,:) ) + atfp1 * zun_e  (:,:)
             zvb_e  (:,:) = atfp * ( zvb_e  (:,:) + va_e  (:,:) ) + atfp1 * zvn_e  (:,:)
@@ -571 +589 @@
       IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
 #endif
+
+      IF ( lk_bdy .or. lk_bdy_tides ) THEN
+         DO jj = 1, jpj
+           DO ji = 1, jpi
+             zhdiv(ji,jj) = zhdiv(ji,jj)*bdytmask(ji,jj)
+           END DO
+         END DO
+      ENDIF
 
       ! sea surface height

trunk/NEMO/OPA_SRC/DYN/wzvmod.F90

@@ -7 +7 @@
    !!             7.0  !  96-01  (G. Madec)  Statement function for e3
    !!             8.5  !  02-07  (G. Madec)  Free form, F90
+   !!              "   !  07-07  (D. Storkey) Zero zhdiv at open boundary (BDY) 
    !!----------------------------------------------------------------------
    !!   wzv        : Compute the vertical velocity
@@ -18 +19 @@
    USE prtctl          ! Print control
    USE phycst
+   USE bdy_oce         ! unstructured open boundaries
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
 
@@ -54 +56 @@
 
       !! * Local declarations
+      INTEGER  ::     jgrd, jb           ! temporary scalars
       INTEGER  ::           jk           ! dummy loop indices
       !! Variable volume
@@ -88 +91 @@
          ! Horizontal divergence of barotropic transports
          !--------------------------------------------------
+         zhdiv(:,:) = 0.e0
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! vector opt.
@@ -105 +109 @@
          IF( lp_obc_north )   zhdiv(nin0  :nin1  ,njn0p1:njn1p1) = 0.e0    ! north
          IF( lp_obc_south )   zhdiv(nis0  :nis1  ,njs0  :njs1)   = 0.e0    ! south
+#endif
+
+#if defined key_bdy || defined key_bdy_tides
+         jgrd=1 !: tracer grid.
+         DO jb = 1, nblenrim(jgrd)
+           ji = nbi(jb,jgrd)
+           jj = nbj(jb,jgrd)
+           zhdiv(ji, jj) = 0.e0
+           zhdiv(ji, jj) = 0.e0
+           zhdiv(ji, jj) = 0.e0
+         END DO
 #endif