New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
sshwzv_trc.F90 in branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC – NEMO

source: branches/2011/dev_r2802_TOP_substepping/NEMOGCM/NEMO/TOP_SRC/sshwzv_trc.F90 @ 2830

Last change on this file since 2830 was 2830, checked in by kpedwards, 13 years ago

Updates to average physics variables for TOP substepping.
File size: 9.3 KB
Line 
1MODULE sshwzv_trc   
2   !!==============================================================================
3   !!                       ***  MODULE  sshwzv  ***
4   !! Ocean dynamics : sea surface height and vertical velocity
5   !!==============================================================================
6   !! History :  3.1  !  2009-02  (G. Madec, M. Leclair)  Original code
7   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  modified LF-RA
8   !!             -   !  2010-05  (K. Mogensen, A. Weaver, M. Martin, D. Lea) Assimilation interface
9   !!             -   !  2010-09  (D.Storkey and E.O'Dea) bug fixes for BDY module
10   !!----------------------------------------------------------------------
11
12   !!----------------------------------------------------------------------
13   !!   ssh_wzv        : after ssh & now vertical velocity
14   !!   ssh_nxt        : filter ans swap the ssh arrays
15   !!----------------------------------------------------------------------
16   USE oce             ! ocean dynamics and tracers variables
17   USE dom_oce         ! ocean space and time domain variables
18   USE sbc_oce         ! surface boundary condition: ocean
19   USE domvvl          ! Variable volume
20   USE divcur          ! hor. divergence and curl      (div & cur routines)
21   USE iom             ! I/O library
22   USE restart         ! only for lrst_oce
23   USE in_out_manager  ! I/O manager
24   USE prtctl          ! Print control
25   USE phycst
26   USE lbclnk          ! ocean lateral boundary condition (or mpp link)
27   USE lib_mpp         ! MPP library
28   USE obc_par         ! open boundary cond. parameter
29   USE obc_oce
30   USE bdy_oce
31   USE diaar5, ONLY:   lk_diaar5
32   USE iom
33   USE sbcrnf, ONLY: h_rnf, nk_rnf   ! River runoff
34   USE trc, ONLY:  nn_dttrc, nittrc000 !TOP_TRC substepping
35#if defined key_agrif
36   USE agrif_opa_update
37   USE agrif_opa_interp
38#endif
39
40   IMPLICIT NONE
41   PRIVATE
42
43   PUBLIC   ssh_wzv_trc    ! called by step.F90
44
45   !! * Substitutions
46#  include "domzgr_substitute.h90"
47#  include "vectopt_loop_substitute.h90"
48   !!----------------------------------------------------------------------
49   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
50   !! $Id: sshwzv.F90 2715 2011-03-30 15:58:35Z rblod $
51   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE ssh_wzv_trc( kt ) 
56      !!----------------------------------------------------------------------
57      !!                ***  ROUTINE ssh_wzv  ***
58      !!                   
59      !! ** Purpose :   compute the after ssh (ssha), the now vertical velocity
60      !!              and update the now vertical coordinate (lk_vvl=T).
61      !!
62      !! ** Method  : - Using the incompressibility hypothesis, the vertical
63      !!      velocity is computed by integrating the horizontal divergence 
64      !!      from the bottom to the surface minus the scale factor evolution.
65      !!        The boundary conditions are w=0 at the bottom (no flux) and.
66      !!
67      !! ** action  :   ssha    : after sea surface height
68      !!                wn      : now vertical velocity
69      !!                sshu_a, sshv_a, sshf_a  : after sea surface height (lk_vvl=T)
70      !!                hu, hv, hur, hvr        : ocean depth and its inverse at u-,v-points
71      !!
72      !! Reference  : Leclair, M., and G. Madec, 2009, Ocean Modelling.
73      !!----------------------------------------------------------------------
74      USE wrk_nemo, ONLY:   wrk_in_use, wrk_not_released
75      USE oce     , ONLY:   z3d   => ta                           ! ta used as 3D workspace
76      USE wrk_nemo, ONLY:   zhdiv => wrk_2d_1 , z2d => wrk_2d_2   ! 2D workspace
77      !
78      INTEGER, INTENT(in) ::   kt   ! time step
79      !
80      INTEGER  ::   ji, jj, jk   ! dummy loop indices
81      REAL(wp) ::   zcoefu, zcoefv, zcoeff, z2dt, z1_2dt, z1_rau0   ! local scalars
82      !!----------------------------------------------------------------------
83
84      IF( wrk_in_use(2, 1,2) ) THEN
85         CALL ctl_stop('ssh_wzv: requested workspace arrays unavailable')   ;   RETURN
86      ENDIF
87
88      IF( kt == nittrc000 ) THEN
89         !
90         IF(lwp) WRITE(numout,*)
91         IF(lwp) WRITE(numout,*) 'ssh_wzv_trc : after sea surface height and now vertical velocity '
92         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
93         !
94         wn(:,:,jpk) = 0._wp                  ! bottom boundary condition: w=0 (set once for all)
95         !
96         !
97      ENDIF
98
99      !                                           !------------------------------------------!
100      IF( lk_vvl ) THEN                           !  Regridding: Update Now Vertical coord.  !   (only in vvl case)
101         !                                        !------------------------------------------!
102         DO jk = 1, jpkm1
103            fsdept(:,:,jk) = fsdept_n(:,:,jk)         ! now local depths stored in fsdep. arrays
104            fsdepw(:,:,jk) = fsdepw_n(:,:,jk)
105            fsde3w(:,:,jk) = fsde3w_n(:,:,jk)
106            !
107            fse3t (:,:,jk) = fse3t_n (:,:,jk)         ! vertical scale factors stored in fse3. arrays
108            fse3u (:,:,jk) = fse3u_n (:,:,jk)
109            fse3v (:,:,jk) = fse3v_n (:,:,jk)
110            fse3f (:,:,jk) = fse3f_n (:,:,jk)
111            fse3w (:,:,jk) = fse3w_n (:,:,jk)
112            fse3uw(:,:,jk) = fse3uw_n(:,:,jk)
113            fse3vw(:,:,jk) = fse3vw_n(:,:,jk)
114         END DO
115         !
116         hu(:,:) = hu_0(:,:) + sshu_n(:,:)            ! now ocean depth (at u- and v-points)
117         hv(:,:) = hv_0(:,:) + sshv_n(:,:)
118         !                                            ! now masked inverse of the ocean depth (at u- and v-points)
119         hur(:,:) = umask(:,:,1) / ( hu(:,:) + 1._wp - umask(:,:,1) )
120         hvr(:,:) = vmask(:,:,1) / ( hv(:,:) + 1._wp - vmask(:,:,1) )
121         !
122      ENDIF
123      !
124      CALL div_cur( kt )                              ! Horizontal divergence & Relative vorticity
125      !
126      z2dt = 2._wp * rdt                              ! set time step size (Euler/Leapfrog)
127      IF( neuler == 0 .AND. kt == nittrc000 )   z2dt = rdt
128
129      !                                           !------------------------------!
130      !                                           !   After Sea Surface Height   !
131      !                                           !------------------------------!
132      zhdiv(:,:) = 0._wp
133      DO jk = 1, jpkm1                                 ! Horizontal divergence of barotropic transports
134        zhdiv(:,:) = zhdiv(:,:) + fse3t(:,:,jk) * hdivn(:,:,jk)
135      END DO
136      !                                                ! Sea surface elevation time stepping
137      ! In forward Euler time stepping case, the same formulation as in the leap-frog case can be used
138      ! because emp_b field is initialized with the vlaues of emp field. Hence, 0.5 * ( emp + emp_b ) = emp
139      z1_rau0 = 0.5 / rau0
140      ssha(:,:) = (  sshb(:,:) - z2dt * ( z1_rau0 * ( emp_b(:,:) + emp(:,:) ) + zhdiv(:,:) )  ) * tmask(:,:,1)
141
142#if defined key_agrif
143      CALL agrif_ssh( kt )
144#endif
145#if defined key_obc
146      IF( Agrif_Root() ) THEN
147         ssha(:,:) = ssha(:,:) * obctmsk(:,:)
148         CALL lbc_lnk( ssha, 'T', 1. )                 ! absolutly compulsory !! (jmm)
149      ENDIF
150#endif
151#if defined key_bdy
152      ssha(:,:) = ssha(:,:) * bdytmask(:,:)
153      CALL lbc_lnk( ssha, 'T', 1. ) 
154#endif
155
156      !                                                ! Sea Surface Height at u-,v- and f-points (vvl case only)
157      IF( lk_vvl ) THEN                                ! (required only in key_vvl case)
158         DO jj = 1, jpjm1
159            DO ji = 1, jpim1      ! NO Vector Opt.
160               sshu_a(ji,jj) = 0.5  * umask(ji,jj,1) / ( e1u(ji  ,jj) * e2u(ji  ,jj) )                   &
161                  &                                  * ( e1t(ji  ,jj) * e2t(ji  ,jj) * ssha(ji  ,jj)     &
162                  &                                    + e1t(ji+1,jj) * e2t(ji+1,jj) * ssha(ji+1,jj) )
163               sshv_a(ji,jj) = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj  ) * e2v(ji,jj  ) )                   &
164                  &                                  * ( e1t(ji,jj  ) * e2t(ji,jj  ) * ssha(ji,jj  )     &
165                  &                                    + e1t(ji,jj+1) * e2t(ji,jj+1) * ssha(ji,jj+1) )
166            END DO
167         END DO
168         CALL lbc_lnk( sshu_a, 'U', 1. )   ;   CALL lbc_lnk( sshv_a, 'V', 1. )      ! Boundaries conditions
169      ENDIF
170     
171
172      !                                           !------------------------------!
173      !                                           !     Now Vertical Velocity    !
174      !                                           !------------------------------!
175      z1_2dt = 1.e0 / z2dt
176      DO jk = jpkm1, 1, -1                             ! integrate from the bottom the hor. divergence
177         ! - ML - need 3 lines here because replacement of fse3t by its expression yields too long lines otherwise
178         wn(:,:,jk) = wn(:,:,jk+1) -   fse3t_n(:,:,jk) * hdivn(:,:,jk)        &
179            &                      - ( fse3t_a(:,:,jk) - fse3t_b(:,:,jk) )    &
180            &                         * tmask(:,:,jk) * z1_2dt
181#if defined key_bdy
182         wn(:,:,jk) = wn(:,:,jk) * bdytmask(:,:)
183#endif
184      END DO
185
186      !
187      IF( wrk_not_released(2, 1,2) )   CALL ctl_stop('ssh_wzv: failed to release workspace arrays')
188      !
189   END SUBROUTINE ssh_wzv_trc
190
191
192   !!======================================================================
193END MODULE sshwzv_trc
Note: See TracBrowser for help on using the repository browser.